%this file computes the Fresnel diffraction of a slit

l=1e-6;     %wavelength
a=100e-6;   %slit half width
k=2*pi/l;
%NF1=1;
%z=a^2/(l*NF1)
z=1e-3;      

xtmp=max([1.15*a,l*z/a]);
x=linspace(-xtmp,xtmp,1001); %position


NF=a^2/(l*z) %Fresnel number


    a1=-sqrt(2)*(sqrt(NF)+x/sqrt(l*z));
    a2=sqrt(2)*(sqrt(NF)-x/sqrt(l*z));

    C1=mfun('FresnelC',a1);
    C2=mfun('FresnelC',a2);
    S1=mfun('FresnelS',a1); 
    S2=mfun('FresnelS',a2);

    I1=.5*((C2-C1).^2+(S2-S1).^2);  %Fresnel intensity



plot(x,I1/max(I1),'b')
hold
plot([-a,-a],[0,1],'k')
plot([a,a],[0,1],'k')
hold
xlabel('position')
ylabel('I(x)')



%is approximation valid?
((max(x)+a)^4)/(l*4*z^3)