from pylab import * from rectangular import * from scipy import * def fra (z,wl): z = float64(z) wl = float64(wl) dxi = float64(0.1) xi = arange(-10,10,dxi) w = float64(5) a = (len(xi)/2) m = arange(-a,a) M = len(m) k = (2*pi)/wl U1 = rect(xi/(2*w))/(w*2) U2 =U1*exp((1j*k*(xi**2))/(2*z)) ##plot function subplot(2,2,1) plot(xi,U1) title("Input Wave Field") ft = fftshift(fft(U2)) ##plot fourier transform subplot(2,2,2) plot (xi,ft) title("Fourier Transform") p = m dx = (wl*z)/(dxi*M) x = p*dx ##phase calculation pIs = exp(1j*z*k)/(1j*wl*z) pDs = exp((1j*k*(x**2))/(2*z)) U3 = pIs*pDs*ft ##Brake down of insensity calculation Nf = (w**2)/(wl*z) print Nf X = x/sqrt(wl*z) alpha1 = -1.00*(sqrt(2.))*(sqrt(Nf)+X) alpha2 = (sqrt(2.))*(sqrt(Nf)-X) Falpha1 = special.fresnel(alpha1) Falpha2 = special.fresnel(alpha2) Salpha1 = Falpha1[0] Calpha1 = Falpha1[1] Salpha2 = Falpha2[0] Calpha2 = Falpha2[1] I = (((Calpha2-Calpha1)**2)+((Salpha2-Salpha1)**2))/(2*wl*z) ##plot Inensity subplot(2,2,3) plot(x,I) title("Book Intensity") ## plot calculated Intensity and Intensity subplot(2,2,4) plot (x,U3*U3.conj()) plot (x,I) title("Book Intensity over Calculated Intensity")