dim=64;
[kx,ky]=meshgrid(-1:2/(dim-1):1);
circ=sqrt(kx.^2+ky.^2)<1;
kx=kx/2;
ky=ky/2;
alp=asin(0.8);
k0=1/sin(alp);
k=256;
kz=sqrt(k0^2-(kx.^2+ky.^2));
Gx=sqrt(k0./kz).*((k0*ky.^2+kz.*kx.^2)./(k0*(kx.^2+ky.^2)));
Gy=sqrt(k0./kz).*((kz-k0).*kx.*ky)./(k0*(kx.^2+ky.^2));
Gz=sqrt(k0./kz).*(kx./k0);
z=-64:1:63;
[r,c]=size(z);
mz=pi*64/(256);
Ex=zeros([k,k,c]);
Ey=zeros([k,k,c]);
Ez=zeros([k,k,c]);
for jj=1:c
Ex(:,:,jj)=fftshift(fft2(exp(1i*kz*z(jj)*mz).*Gx.*circ,k,k));
Ey(:,:,jj)=fftshift(fft2(exp(1i*kz*z(jj)*mz).*Gy.*circ,k,k));
Ez(:,:,jj)=fftshift(fft2(exp(1i*kz*z(jj)*mz).*Gz.*circ,k,k));
end
I =Ex.*conj(Ex)+Ey.*conj(Ey)+Ez.*conj(Ez);
M=I(k/2,:,:);
M1=squeeze(M).';
figure(1)
imagesc((M1));colormap gray