%slit diffraction from Rayleigh-Sommerfeld based on recursive adaptive Simpson
%quadrature quad(F,lim1,lim2)
%assumes observation position x2,z2, and wavelength=1
%calls RSkernel as a function, where RSkernel is the Rayleigh Sommerfeld
%integration kernel
clear u
slitwidth=50
tol=1e-7;
z2=100;
x2vals=linspace(-slitwidth/2,slitwidth/2,1000);
N=length(x2vals);
tic
for p=1:N
x2=x2vals(p);
u(p)=quad(@(x)RSkernel(x,x2,z2),-slitwidth/2,slitwidth/2,tol);
end;
toc
%plot the intensity, 
plot(x2vals,z2*abs(u).^2,'g')
