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")