##############################################################
# 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 2013 Lumerical Solutions Inc.
##############################################################

# run simulation;
run;
closeall;  clear;

##############################################################
# User Settings
NAin = 0.8;         	# choose NA of imaging optics, input lens
NAout = 1;         	# choose NA of imaging optics, output. generally leave this at 1
magnification = 1;     # magnification factor
mname = "T"; 	# monitor recording the data
res_ff = 2001;     	# resolution for far field projection
res_image = 1001;      # resolution of near field image
y0 = 0e-6;             # defocus. This specifies the location of the focal plane in global coordinates
index = 1;             # far field refractive index
##############################################################


# get simulation frequency
f  = getdata(mname,"f");  	# frequency of monitor
nf=length(f);
x=getdata(mname,"x");
lamda=c/f;          		# wavelength
k=2*pi*index/lamda;       	# k
E2_image = matrix(res_image,length(f));	# initialize matrix to store E2 data on image plane


for (i=1:nf) {  # loop over frequency
  ?i;

  # decompose nearfield into plane waves using farfield projection
  E = farfieldpolar2d(mname,i,res_ff);
  Er = pinch(E,2,1);
  Et = pinch(E,2,2);
  Ez = pinch(E,2,3);
  theta=farfieldangle(mname,i,res_ff);
  ux=sin(theta*pi/180);
  uy=real(sqrt(1-ux^2))+1e-20;

  # filter for propagating waves through aperture
  filter = abs(ux) < NAin;  
  Er=filter*Er;
  Et=filter*Et;
  Ez=filter*Ez;
  #?farfield2dintegrate(abs(Er)^2+abs(Et)^2+abs(Ez)^2,theta);  # check far field power

  # apply the magnification and output NA
  ux = ux/magnification;
  uy = real(sqrt(1-ux^2))+1e-20;

  filter = abs(ux) < NAout;
  Er=filter*Er*sqrt(magnification);
  Et=filter*Et*sqrt(magnification);
  Ez=filter*Ez*sqrt(magnification);
  #farfield2dintegrate(abs(Er)^2+abs(Et)^2+abs(Ez)^2,asin(ux)*180/pi);  # check far field power

  # return to cartesian coordinates so we can sum the fields
  theta_rad = asin(ux);
  Ex = Er*sin(theta_rad) - Et*cos(theta_rad);
  Ey = Er*cos(theta_rad) + Et*sin(theta_rad);

  # define image plane 
  x_image=linspace(min(x),max(x),res_image)*magnification;

  # calculate the image, uzing chirped z-transform to sum the individual plane waves
  kx = ux*k(i);
  ky = uy*k(i);
  dkx = kx(2)-kx(1);

  Ex_image = czt(Ex*exp(1i*ky*y0)/uy,kx,x_image)*dkx/sqrt(2*pi*k(i));
  Ey_image = czt(Ey*exp(1i*ky*y0)/uy,kx,x_image)*dkx/sqrt(2*pi*k(i));
  Ez_image = czt(Ez*exp(1i*ky*y0)/uy,kx,x_image)*dkx/sqrt(2*pi*k(i));

  # calculate |E|^2 at the image plane
  E2_image(1:res_image,i) =abs(Ex_image)^2 + abs(Ey_image)^2 + abs(Ez_image)^2;
}

## plot fields at image plane 
#plot(x_image*1e6,E2_image,interp(getelectric(mname),getdata(mname,"x"),x_image),"x (um)","E2");
#legend("E2 from microscope image","E2 from monitor");

if(nf==1) {
  plotxy(x_image*1e6,E2_image,x*1e6,getelectric("test1"),"x (um)","E2");
  legend("E2 image plane","E2 from monitor");
} else {
  image(x_image*1e6,c/f*1e6,E2_image,"x (um)","wavelength","E2");
}