% function edgedif
% Edge diffraction integral, Eq. (24) in the notes about diffraction 
% u is normalized position in aperture
u  = -2:.05:8.; % compute integral in 201 positions
%
% first simple sum approximation to integral
e0 =  (u(2)-u(1)).*cumsum (edgeint(u)); 
% the Fresnel integral is (1 - i)/2 for u = 0
e0 =  e0 - e0(find(u==0)) + (1. - i)/2.;
i0 =  .5.*e0.*conj(e0);  % intensity is field amplitude squared
%
% then more accurate numerical integration
e1 =  zeros(1,201);  % preallocate arrays of integral values
for l = 2:length(u)
   e1(l) = e1(l-1) + quad ('edgeint', u(l-1), u(l));
end

% function eint = edgeint (u)
% Integrand needed in edge diffraction integral, Eq. (24) in 
% the note about Fresnel diffraction 
% Computes integrand to be used with function edgedif
% Input:  
%  u  array of values to compute integrand for
%     normalized position in aperture
% Output: 
%  eint  array of corresponding values for the integrand
% eint  =  exp (-0.5.*j.*pi.*u.^2); 
% return

e1 = e1 - e1(find(u==0)) + (1. - i)/2.;
i1 = .5.*e1.*conj(e1);  % intensity is field amplitude squared
plot (u, i0, 'b:', u, i1, 'k');
title ('Edge Diffraction Pattern')
xlabel ('Position')
ylabel ('Intensity')
