#########################################################################
#
# Description: This file should be used in conjunction with 
#   propagate periodic simulation file to test the method of projecting
#   fields through homogeneous materials for scattering from
#   periodic structures
#
# Copyright 2010 Lumerical Solutions Inc.
########################################################################

# test it with the FDTD actual simulation or not, 0 to not use test_monitor results
test=1;

# only for not to use the test_monitor results (test=0)
z_project_res = 1000;
z_start = 1e-007;
z_end = 100e-006;

# the name of the monitor recording the data
mname = "transmission";
index = 1;

# choose the desired wavelength
lambda_target = 0.5e-6;

# image the near field at the object plane
E2 = getelectric(mname);
x = getdata(mname,"x");
y = getdata(mname,"y");
z = getdata(mname,"z");
f = getdata(mname,"f");
fpoint = find(c/f,lambda_target);
f = f(fpoint);
E2 = pinch(E2,4,fpoint);
power1 = transmission(mname);
power1 = power1(fpoint)*sourcepower(f); # power in Watts
image(x*1e6,y*1e6,E2,"x (microns)","y (microns)","Object Field Intensity");

# calculate the field strength of the grating orders
E = gratingvector(mname,fpoint);
Ex = pinch(E,3,1);
Ey = pinch(E,3,2);
Ez = pinch(E,3,3);

# calculate wavelength and k-vector
lambda=c/f;     # wavelength
k=2*pi/lambda;  # k-vector

## calculate the grating order direction cosines
# and various derived quantities
ux = gratingu1(mname,fpoint);
uy = gratingu2(mname,fpoint);
kx = k*ux;
ky = k*uy;
Ux = meshgridx(ux,uy);
Uy = meshgridy(ux,uy);
Uxy = sqrt(Ux^2+Uy^2)+1e-5; # add 1e-5 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;

# correct fields for sqrt(cos(theta)) factor
Ex = Ex/(sqrt(cos_theta)+1e-5);
Ey = Ey/(sqrt(cos_theta)+1e-5);
Ez = Ez/(sqrt(cos_theta)+1e-5);

# calculate filter to include only propagating plane waves 
# (ie eliminate evanescent waves)
filter = (Ux^2+Uy^2 < 1);

# make sure tha total power of plane waves is equal to total power of
# the original monitor
area = (max(x)-min(x))*(max(y)-min(y));
power2 = 0.5*index*sqrt(eps0/mu0)*sum(cos_theta*(abs(Ex)^2+abs(Ey)^2+abs(Ez^2)))*area;
norm_factor = sqrt(power1/power2);
Ex = Ex*norm_factor;
Ey = Ey*norm_factor;
Ez = Ez*norm_factor;

###############################################################
# perform test to compare with FDTD simulated result
###############################################################

# get the FDTD result OR only use project
if (test){
xt = getdata("test_monitor","x");
yt = getdata("test_monitor","y");
zt = getdata("test_monitor","z");
E2t = getelectric("test_monitor");
E2t = pinch(E2t,4,fpoint);
Ext = getdata("test_monitor","Ex");
Ext = pinch(Ext,4,fpoint);
tempx = linspace(min(xt),max(xt),length(xt));
tempy = linspace(min(yt),max(yt),length(yt));
tempz = linspace(min(zt),max(zt),length(zt));
E2t = interp(E2t,xt,yt,zt,tempx,tempy,tempz);
Ext = interp(Ext,xt,yt,zt,tempx,tempy,tempz);
xt = tempx;
yt = tempy;
zt = tempz;

E2t = pinch(E2t,2,find(yt,0));
Ext = pinch(Ext,2,find(yt,0));
image(xt*1e6,zt*1e6,E2t,"x (microns)","z (microns)","|E|^2 by FDTD at " + num2str(lambda*1e9) + "nm");
E2p = 0*E2t;
Exp = 0*E2t;

} else {
zt=linspace(z_start, z_end, z_project_res);
xt = getdata("test_monitor","x");
yt = getdata("test_monitor","y");
E2p=matrix(length(xt),length(zt));
Exp=matrix(length(xt),length(zt));
}

# calculate projected result
yt = yt(find(yt,0):(find(yt,0)+1));
for(i=1:length(zt)) {
  ?"calculating projection " + num2str(i) + " of " + num2str(length(zt));
  
  # find phase correction factor
  # the i in front of exp(-ikx) is because of backward propagating the far field (1m) to the near field, where i=exp(i*pi/2).
  z_phase = 1i*exp(1i*k*Uz*filter*(zt(i))-1i*k*1);

  Ex_image = czt(Ex*z_phase*filter,kx,ky,xt,yt);
  Ey_image = czt(Ey*z_phase*filter,kx,ky,xt,yt);
  Ez_image = czt(Ez*z_phase*filter,kx,ky,xt,yt);
  E2p(1:length(xt),i) = pinch(abs(Ex_image)^2+abs(Ey_image)^2+abs(Ez_image)^2 ,2 ,1);
  Exp(1:length(xt),i) = pinch(Ex_image ,2 ,1);
}

image(xt*1e6,zt*1e6,E2p,"x (microns)","z (microns)","|E|^2 by projection at " + num2str(lambda*1e9) + "nm");
if (test) {
plot(zt*1e6,pinch(E2p,1,find(xt,0)),pinch(E2t,1,find(xt,0)),"zt (um)","|E|^2","|E|^2 at " + num2str(lambda*1e9) + "nm");
legend("projection","FDTD");
plot(zt*1e6,pinch(real(Exp),1,find(xt,0)),pinch(real(Ext),1,find(xt,0)),"zt (um)","Ex","Re(Ex) at " + num2str(lambda*1e9) + "nm");
legend("projection","FDTD");
} else {
plot(zt*1e6,pinch(E2p,1,find(xt,0)),"zt (um)","|E|^2","|E|^2 at " + num2str(lambda*1e9) + "nm");
legend("projection");
}