from math import cos,sin,pi,sqrt
import matplotlib.pyplot as plt

def simpson(xl, xr, ndiv, fun): #liczy całke
        h=(xr-xl)/ndiv/2.
        integ=0.
        for i in range(0,ndiv):
                integ=integ+2.*h*(fun(xl+i*2.*h)+4.*fun(xl+i*2.*h+h)+fun(xl+i*2.*h+2*h))/6.
        return integ

def real_amp(x):
        r=sqrt(d*d+(y-x)**2)
        return ampl*cos(k*r)/(sqrt(r))

def imag_amp(x):
        r=(d*d+(y-x)**2)**(1./2.)
        return ampl*sin(k*r)/(sqrt(r))


 
ndiv=100 #kroki w simpsonie
ampl=1. #amplituda
lmbd=2 #długosc fali
k=2.*pi/lmbd #liczba falowa

a=1. #polozenie szczelin
d=20. #odleglosc
yl=-30.
yr=30.
ny=1000 #dokładnosc sumowana amplitudy



hy=(yr-yl)/(ny-1)

yarr=[] #os
iarr=[] #ampiltuda z pierwszej szczeliny
jarr=[] #ampiltuda z drugiej szczeliny
karr=[] #suma amplitud

for iy in range(1,ny):
    y=yl+(iy-1.)*hy 
    yarr.append(y)
    i=simpson(-5.*a,-3.*a,ndiv,real_amp)+simpson(-5.*a,-3.*a,ndiv,imag_amp)
    j=simpson(3.*a,5*a,ndiv,real_amp)+simpson(3.*a,5.*a,ndiv,imag_amp)

    iarr.append(i)
    jarr.append(j)
    #yarr.append(y)
    karr.append((i+j)**2) #natezenie
    
#n=2
   
#for p in range(0,ny-1):
#    k=0
 #   k=k+(iarr[p]+jarr[p]) 
 #   karr.append(k)
    
 
plt.xlabel("odleglosci")
plt.title("Dyfrakcja na dwoch szczelinach")    
plt.plot(yarr,karr)