%this file compares diffraction for the Fresnel and
%Fraunhoffer approximations

l=1e-6;     %wavelength
Wx=100e-6;  %slit width

%NF1=10;
%z=Wx^2/(l*NF1);
z=.05;       %for Fresnel z=0.0022      
            %for Fraunhoffer z=0.1

xtmp=max([1.5*Wx,l*z/Wx]);
x=linspace(-xtmp,xtmp,101); %position


NF=Wx^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

I2=(sinc(2*Wx*x/(l*z))).^2;     %Fraunhoffer intensity

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


a=max(x)+Wx;
%fraunhoffer approximation
approx_fraun=Wx^2/(2*l*z);

%fresnel approximation
approx_fresnel=(a^4/(4*l*z^3));

[a z]

[approx_fresnel,approx_fraun]