import numpy as np
import matplotlib.pyplot as plt
plt.interactive(True)


     


'''
	Main program diifraction-1.py
    
    Compute the diffracion pattern from sigle slit
    
    X = 25e-6          - half-width of the slit
    lambda = 0.6e-6 m  - red light wavelength
    z = 1 m            - distance to screen
    y                  - coordinte on the screen
    Range of y:  approx. -5 to +5 cm,
    Compute only for y>0, use symmetry.
    Also, try to use vector operations on arrays,
    and np.sum() summation, 
    this will greatly speed up the integration
    
    Formula for intensity I(y) on the screen:
    I(y) = |f(y)|^2,   where f(y) is a complex-valued sum 
    f(y) = (z^2+y^2)^{-1/2} integral_{-X}^{+X} exp[i*k*delta(x)] dx, 
    delta = sqrt(z^2+(y-x)^2) - z
    k = 2 pi / lambda 
    i = 1j in Python = imaginary unit
    Use 1001 points both on the screen and across the slit
'''

plt.figure(dpi=140)
X, Y, lamb, z = 15e-6, 0.06, 0.6e-6, 1. 
k = 2*np.pi/lamb        # wave number
ik = 1j * k
z2 = z*z
N = 1001
x = np.linspace(-X,X,N); dx = x[1]-x[0]
y = np.linspace(-Y,Y,N)
I = y*0   # initialize intensity array I(y)

print('ik',ik)


# fill I(y) array
for l in range(N//2): 
    yy = y[N//2+l]
    deltas = (np.sqrt(z2+(yy-x)**2) - z)  # array of N floats
    integral  = np.sum(np.exp(ik*deltas)) * dx
    f = integral/(z2+yy**2)**0.5
    intens = abs(f)**2
    I[N//2+l] = intens
    I[N//2-l] = intens
I = I/I.max()
# output
plt.plot(100*y,I)
plt.plot(100*y,y*0,linewidth=0.2,color=(.7,.7,.7))
plt.title('Diffraction pattern of a singe slit '+str(2e6*X)+' um wide')
plt.xlabel('y [cm]')
plt.ylabel('Intensity')
plt.grid()
plt.show()
input('next plot? ')