function [x,uz,Iz]=oneDdiffract(slitwidth,z)
%function [x,uz,Iz]=oneDdiffract(slitwidth,z)
%calculates the one dimensional diffraction pattern of a slit 
%	slitwidth is the width of the slit in wavelengths
%	z is the distance behind the screen in wavelengths


mult=8; %this zero-pads to an input vector that is mult x slitwidth
if z>1000;
   mult=4*2^floor(log(sqrt(z)/4)/log(2)); %for large z, mult may need to grow
end;

X=mult*slitwidth;   %This is full X scale length
N=18;   %I'm using 2^N values in my array

points=([1:2^N]-2^(N-1)-0.5)./(2^N);    %symmetric about zero, with 2^N points
x=X*points;         %scaled x variable
rect=find(abs(x)<=slitwidth/2); %these three lines set a rectangular aperture
aper=zeros(size(x));
aper(rect)=1;
%aper=aper.*(1+cos(2*pi*x*10/slitwidth));  %other interesting aperture functions
%aper=aper.*(1-abs(2*x./slitwidth));

A=fftshift(fft(aper));  
f=([0:2^N-1]-2^(N-1))./X;   %this sets the physical frequency, with zero corresponding to an fftshift'ed array
propfreq=find(abs(f)<1);    %the next three lines mask out non-propagating terms
propmask=zeros(size(f));
propmask(propfreq)=ones(size(propfreq));

Az=A.*exp(j*2*pi*z*sqrt(1-f.^2));   %this does the propagation

uz=ifft(fftshift(Az .*propmask));   %this does the inverse transform - note the fftshift again
Iz=uz .*conj(uz);   %predicted intensity pattern

plot(x,Iz,'r-');
%plot(x(2^(N-1)-2^(N-4):2^(N-1)+2^(N-4)),Iz(2^(N-1)-2^(N-4):2^(N-1)+2^(N-4)),'g');