##############################################################
# scriptfile: microscopy_imaging.lsf
#
# Description:  This script decomposes light propagating 
#               through a monitor into a ray set.  It then
#               applies a filter to remove all light at 
#               steep angles (set by the NA).  The remaining 
#               light is re-focused onto an image plane.
#
# Copyright 2007 Lumerical Solutions Inc.
##############################################################

##############################################################
# User Settings
NA = 1;         	# choose NA of imaging optics
mname = "Monitor1"; 	# monitor recording the data
res = 201;             # res controls the number of plane waves
                       #    used for decomposition, which is (res x res)
z0 = 0e-6;            # defocus
##############################################################


# get simulation frequency
f  = getdata(mname,"f");  # frequency of monitor
x=getdata(mname,"x");
y=getdata(mname,"y");
lambda=c/f;          		# wavelength
k=2*pi/lambda;       		# k

# decompose nearfield into plane waves using farfield projection
E = farfieldvector3d(mname,1,res,res);
Ex = pinch(E,3,1);
Ey = pinch(E,3,2);
Ez = pinch(E,3,3);
ux = farfieldux(mname,1,res,res);
uy = farfielduy(mname,1,res,res);
Ux = meshgridx(ux,uy);
Uy = meshgridy(ux,uy);
Uz = real(sqrt(1-(Ux^2+Uy^2))) + 1e-5;

# filter for propagating waves through aperture
filter = real(sqrt(Ux^2 + Uy^2)) < NA;  
Ex=filter*Ex;
Ey=filter*Ey;
Ez=filter*Ez;

# define image plane 
x_image=linspace(min(x),max(x),res);
y_image=linspace(min(y),max(y),res);


# calculate the image, uzing chirped z-transform
kx = ux*k;
ky = uy*k;
dkx = kx(2)-kx(1);
dky = ky(2)-ky(1);
Kz = Uz*k;

Ex_image = czt(Ex*exp(1i*Kz*z0)*sqrt(1/Uz/k^2),kx,ky,x_image,y_image)*dkx*dky/(2*pi);
Ey_image = czt(Ey*exp(1i*Kz*z0)*sqrt(1/Uz/k^2),kx,ky,x_image,y_image)*dkx*dky/(2*pi);
Ez_image = czt(Ez*exp(1i*Kz*z0)*sqrt(1/Uz/k^2),kx,ky,x_image,y_image)*dkx*dky/(2*pi);

# calculate |E|^2 at the image plane
E2_image =abs(Ex_image)^2 + abs(Ey_image)^2 + abs(Ez_image)^2;

# plot fields at image plane
image(x_image*1e6,y_image*1e6,E2_image,              
    "x (um)","y (um)","Image. NA="+num2str(NA)+", defocus="+num2str(z0*1e6));

# plot fields in near field
image(x*1e6,y*1e6,getelectric(mname),
    "x (um)","y (um)","Near field");