clear all;
cc = 2.99792458E8;
x=linspace(-10e-6,10e-6,401);
y=linspace(-11e-6,11e-6,441);
X=(meshgrid(x,y)).';
Y=meshgrid(y,x);
lambda0=0.5e-6;
f=cc/lambda0;
w=2*pi*f;
k=2*pi/lambda0;
NA=0.2;
kx=linspace(-k,k,201);
ky=linspace(-k,k,221);
Kx=(meshgrid(kx,ky)).';
Ky=meshgrid(ky,kx);
phi=atan2(Ky,Kx);
theta=real(acos(sqrt(1-Kx.^2./k^2-Ky.^2/k^2)));
envelope=exp(-0.5.*(Kx.^2+Ky.^2)./(NA.*k)^2);
Exk=cos(phi).*cos(theta).*envelope;
ny=length(y);
fy1=min(y);
fy2=max(y);
samp_freq_y=1/(ky(2)-ky(1));
wy=exp(1i*2*pi*(fy2-fy1)/(samp_freq_y*ny));
ay=exp(1i*2*pi*fy1/samp_freq_y);
tmp=czt(Exk.',ny,wy,ay);
Exf=(tmp).';
figure;imagesc(abs(Exf));
nx=length(x);
fx1=min(x);
fx2=max(x);
samp_freq_x=1/(kx(2)-kx(1));
wx=exp(1i*2*pi*(fx2-fx1)/(samp_freq_x*nx));
ax=exp(1i*2*pi*fx1/samp_freq_x);
Exf=czt(tmp,nx,wx,ax);
figure;
imagesc(abs(Exf));
