//Infinite-distance angular spectrum integrator
//For a finite spot circular laser source

theta=0.00001;//observation angle [rad]
R=3e-3; //laser aperture radius [m]
lam=635e-9;//laser wavelength [m]

function W=circ_intg(R,theta,lam)
    x=linspace(-R*cos(theta),R*cos(theta),1000);
    y=sqrt(R^2-((1/cos(theta))*x).^2);
    Wi=exp(%i*2*%pi*(x*sin(theta)/lam));
    dA=(2*y)/cos(theta);
    W=inttrap(x,Wi.*dA)/(%pi*R^2);
endfunction

function W=circ_intg2(R,theta,lam)
    x=linspace(-R,R,1000);
    y=sqrt(R^2-x.^2);
    Wi=cos(2*%pi*(x*sin(theta)/lam));
    W=inttrap(x,2*Wi.*y)/(%pi*R^2);
endfunction

tt=0.00001*(0:500);
Ww=zeros(1,length(tt));
for(nt=1:length(tt))
    W=circ_intg2(R,tt(nt),lam);
    Ww(nt)=W;
end

//disp(W)
plot(tt,real(Ww),'r');
