%  t=0:0.01:20;
% plot(t,besselj(0,t))



% lambda wavelength  
% A  aperture      
% z2   distance 

 lambda =1 ;
 %  A  = 10   ;  A=3;
 z2   = 50;
 r=0:lambda/20:lambda*10;

intens=r*0;

u= 2*pi/lambda * z2 * (A/z2 )^2;
v=2*pi/lambda *(A/z2 ) * r;

for m=1:max(size(r)) 
     % evaluate integral
      t=0:0.001:1;

       Func=t.* besselj(0,v(m)*t)   .*    exp (i* 0.5*t.^2 * u);
     %  figure(1); plot(t,Func)  ; drawnow
       
       intensr(m)=sum(real(Func));
       intensi(m)=sum(imag(Func));
       %pause
end    
% figure(2), plot(r,intensr,'b-',r,intensi,'r-')
figure(2), plot(r,intensr.^2 + intensi.^2,'r-')