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

//units are in terms of wavelengths
theta=linspace(0,0.001,100); //angles to measure [rad]
R=1000; //aperture radius [wvl]
x=linspace(-R,R,100); //oscillators to add
y=sqrt(R.^2-x.^2); //amount of oscillators in y
sim0=ones(1,length(x)); //initial simultaneity plane (complex phase + ampl at each x)

curv=0e-9;
//dz=curv.*x.^2; //shift of initial plane
dz=x.*sin(-0);
Wi=exp(%i*2*%pi*dz); //phase shifts at each oscillator
sim0=Wi.*sim0; //modified simultaneity plane

nt=1;
Ww=zeros(1,length(theta));

for(th1=theta)
    dz=x.*sin(th1); //shift of path lengths vs x [wvl]
    Wi=exp(%i*2*%pi*dz); //phase shifts at each oscillator
    sim1=Wi.*sim0; //modified simultaneity plane from view of th1
    W=2*inttrap(x,sim1.*y)/(%pi*R^2); //add up all oscillators weighed by their areal density 2D
    //W=inttrap(x,sim1)/(2*R); //add up all oscillators weighed by their areal density 1D
    Ww(nt)=W;
    nt=nt+1;
end

//disp(W)
//plot(theta,abs(Ww));
plot(theta,real(Ww));
//plot(theta,imag(Ww),'g');
