#########################################################################
# scriptfile: phase_constrast_analysis;
#
# Description: Script to calculate image field from object field
#              of a phase contrast microscope with a magnification
#              of 4X for an incoherent, unpolarized source. This
#              script assumes that the script
#              phase_contrast_batch_run.lsf has already been run.
#
# Copyright 2007 Lumerical Solutions Inc.
########################################################################

# choose the filename of the template file
filename = "phase_contrast1";

# choose magnification factor
M = 4;

# choose NA of imaging optics
NA_i = 0.8;

# the name of the monitor recording the data
mname = "Monitor2";

# choose the number of azimuthal angles used in phase_contrast_batch_run.lsf
points = 4;

# choose the resolution for far field projections
res = 201; # resolutions is res x res

# choose the resolution for the phase offset plots in -90 to 90 degree range (phase plate)
phase_res = 5;

# loop over the points and polarization states
simulation_counter = 0;
for(i=1:points) {
    for(polarization=0:90:90) {
        simulation_counter = simulation_counter + 1;
        ?"Analyzing simulation " + num2str(simulation_counter) + " of " + num2str(2*points);

        #####################################################
        # calculate the specimen beam
        #####################################################
            load(filename+"_"+num2str(simulation_counter));     

        # image the near field for the first simulation
            if(simulation_counter==1) {
                ## get data from monitor to plot phase of field at object
                f  = getdata(mname,"f");
                Ex = getdata(mname,"Ex");
                E2 = getelectric(mname);
                x = getdata(mname,"x");
                y = getdata(mname,"y");
                image(x*1e6,y*1e6,E2,"x (microns)","y (microns)","Object Field Intensity");
                image(x*1e6,y*1e6,angle(Ex),"x (microns)","y (microns)","Object Field Phase");
            }

		# calculate the field strength of the grating orders
			Es = gratingpolar(mname,1);
			Er = pinch(Es,3,1);
            Et = pinch(Es,3,2);
            Ep = pinch(Es,3,3);
	# normalize such that the sum of the orders is equal to the
        # normalized transmission
            norm_factor = sqrt( transmission(mname));

            Er = Er * norm_factor;
            Et = Et * norm_factor;
            Ep = Ep * norm_factor;       

        #####################################################
        # calculate the reference beam
        #####################################################
            load(filename+"_reference_"+num2str(simulation_counter));  

		# calculate the field strength of the grating orders
            Es_ref = gratingpolar(mname,1);
            Er_ref = pinch(Es_ref,3,1);
            Et_ref = pinch(Es_ref,3,2);
            Ep_ref = pinch(Es_ref,3,3);

        ## calculate the grating order direction cosines
			ux = gratingu1(mname);
            uy = gratingu2(mname);

        # for the reference beam, force only the zeroth order
        # other orders are only present due to numerical error on the interpolation
            gn = gratingn(mname);
            gm = gratingm(mname);
	    filter = 0*Et_ref;
            filter(find(gn,0),find(gm,0)) = 1;
	    Er_ref = Er_ref*filter;
	    Et_ref = Et_ref*filter;
	    Ep_ref = Ep_ref*filter;
	# normalize such that the sum of the orders is equal to the
        # normalized transmission
            norm_factor = sqrt( transmission(mname) / sum(abs(Er_ref)^2 + abs(Et_ref)^2 + abs(Ep_ref)^2)  );
            Er_ref = Er_ref * norm_factor;
            Et_ref = Et_ref * norm_factor; 
            Ep_ref = Ep_ref * norm_factor;

        #####################################################
        # calculate field at image plane
        #####################################################
            lamda=c/f;     # wavelength
            k=2*pi/lamda;

            Ux = meshgridx(ux,uy);
            Uy = meshgridy(ux,uy);

        ## apply magnification factor
            ux = ux/M;
            uy = uy/M;
            Ux=(Ux)/M;
            Uy=(Uy)/M;

        ## filter for propagating waves through aperture
            filter = real(sqrt(Ux^2 + Uy^2)) < NA_i/M;  

            Uxy = sqrt(Ux^2+Uy^2)+1e-20; # add 1e-20 to avoid divide by zero problems
            Uz  = sqrt(1-Uxy^2);
            sin_phi = Uy/Uxy;
            cos_phi = Ux/Uxy;
            cos_theta = Uz;
            sin_theta = Uxy;

        ## calculate magnified field components in image plane
            Ex = (-Ep*sin_phi + Et*cos_phi*cos_theta);
            Ey = ( Ep*cos_phi + Et*sin_phi*cos_theta);
            Ez =  -Et*sin_theta;
            Ex_ref = (-Ep_ref*sin_phi + Et_ref*cos_phi*cos_theta);
            Ey_ref = ( Ep_ref*cos_phi + Et_ref*sin_phi*cos_theta);
            Ez_ref =  -Et_ref*sin_theta;

        ## define image plane 
            npts=200;
            simulation; 
            xmin=getnamed("FDTD","x min")*M;
            xmax=getnamed("FDTD","x max")*M;
            ymin=getnamed("FDTD","y min")*M;
            ymax=getnamed("FDTD","y max")*M;

            x=linspace(xmin,xmax,npts);
            y=linspace(ymin,ymax,npts);

        # calculate the image, uzing chirped z-transform
            kx = ux*k;
            ky = uy*k;
            Ex_image = czt(filter*Ex,kx,ky,x,y);
            Ey_image = czt(filter*Ey,kx,ky,x,y);
            Ez_image = czt(filter*Ez,kx,ky,x,y);
            Ex_ref_image = czt(filter*Ex_ref,kx,ky,x,y);
            Ey_ref_image = czt(filter*Ey_ref,kx,ky,x,y);
            Ez_ref_image = czt(filter*Ez_ref,kx,ky,x,y);

        ## calculate |E|^2 at the image plane
            if(simulation_counter==1) {
                phase_delay = linspace(-90,90,phase_res)*pi/180;
                E2_pc_image = matrix(npts,npts,length(phase_delay));
                E2_image = 0;
                E2_ref_image = 0;
            }
	 # adding up the intensities for all injection angles
            E2_image = E2_image + abs(Ex_image)^2 + abs(Ey_image)^2 + abs(Ez_image)^2;
            E2_ref_image = E2_ref_image + abs(Ex_ref_image)^2 + abs(Ey_ref_image)^2 + abs(Ez_ref_image)^2;
            for(j=1:length(phase_delay)) {
                E2_pc_image(1:npts,1:npts,j) = pinch(E2_pc_image(1:npts,1:npts,j)) 
                                               + abs(Ex_image+exp(1i*phase_delay(j))*Ex_ref_image)^2
                                               + abs(Ey_image+exp(1i*phase_delay(j))*Ey_ref_image)^2
                                               + abs(Ez_image+exp(1i*phase_delay(j))*Ez_ref_image)^2;
            }
    }
}

## plot field at image plane with no re-interference, ie bright field result
   image(x*1e6,y*1e6,E2_image/max(E2_image),"x (microns)","y (microns)","Bright field microscope image");

# image the results for different phase delays
for(j=1:length(phase_delay)) {
    image(x*1e6,y*1e6,pinch(E2_pc_image,3,j)/max(E2_image),"x (microns)","y (microns)","Phase offset: "+num2str(phase_delay(j)*180/pi));
}

# reload the template file
load(filename);