clear;
closeall;
#########################################################################
# scriptfile: propagate_periodic.fsp
#
# Description: This file can be used to project the fields from a planar
#   x-y FDTD monitor to any desired distance for periodic structures.
#
# Copyright 2010 Lumerical Solutions Inc.
########################################################################

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

# specify the z-position at which to reconstruct the fields
z_project = 99.83e-6;

# specify the index of the material in which the monitor is located
index = 1;


# 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");
power1 = transmission(mname);
power1 = power1*sourcepower(f); # power in Watts

# initializing field component at image
size_correction = 2; # extend space one step outside of simulation region
Ex = matrix(length(x) + size_correction,  length(y) + size_correction, length(f) );
Ey = matrix(length(x) + size_correction,  length(y) + size_correction, length(f) );
Ez = matrix(length(x) + size_correction,  length(y) + size_correction, length(f) );
Hx = matrix(length(x) + size_correction,  length(y) + size_correction, length(f) );
Hy = matrix(length(x) + size_correction,  length(y) + size_correction, length(f) );
Hz = matrix(length(x) + size_correction,  length(y) + size_correction, length(f) );

# calculate the field strength of the grating orders
E = gratingvector(mname,1:length(f));
Exf = pinch(E,3,1);
Eyf = pinch(E,3,2);
Ezf = 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
ux = gratingu1(mname,1:length(f));
uy = gratingu2(mname,1:length(f));

for (i = 1 : length(f)) { # loop over all wavelengths

# field components for each i
Exi = pinch(Exf(:,:,i));
Eyi = pinch(Eyf(:,:,i));
Ezi = pinch(Ezf(:,:,i));

# calculate the grating order direction cosines
# and various derived quantities
kx = k(i)*ux(:,i);
ky = k(i)*uy(:,i);
Ux = meshgridx(ux(:,i),uy(:,i));
Uy = meshgridy(ux(:,i),uy(:,i));
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
Exi = Exi/(sqrt(cos_theta)+1e-5);
Eyi = Eyi/(sqrt(cos_theta)+1e-5);
Ezi = Ezi/(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(Exi)^2+abs(Eyi)^2+abs(Ezi^2)))*area;
norm_factor = sqrt(power1(i)/power2);
Exi = Exi*norm_factor;
Eyi = Eyi*norm_factor;
Ezi = Ezi*norm_factor;

Hxi = index*sqrt(eps0/mu0)*(Uy*Ezi - Uz*Eyi);
Hyi = index*sqrt(eps0/mu0)*(Uz*Exi - Ux*Ezi);
Hzi = index*sqrt(eps0/mu0)*(Ux*Eyi - Uy*Exi);

# Perform the projection find phase correction factor
# extend space one step outside of simulation region
# and ensure a linearly space vector
dx = max(x(2:length(x))-x(1:length(x)-1));
dy = max(y(2:length(y))-y(1:length(y)-1));
if (i==1){ # only need to do this operation once
x = linspace(min(x)-dx,max(x)+dx,length(x)+2);
y = linspace(min(y)-dx,max(y)+dy,length(y)+2);
}

z_phase = 1i*exp(1i*k(i)*Uz*filter*(z_project)-1i*k(i)*1);
Ex_image = czt(Exi*z_phase*filter,kx,ky,x,y);
Ey_image = czt(Eyi*z_phase*filter,kx,ky,x,y);
Ez_image = czt(Ezi*z_phase*filter,kx,ky,x,y);
Hx_image = czt(Hxi*z_phase*filter,kx,ky,x,y);
Hy_image = czt(Hyi*z_phase*filter,kx,ky,x,y);
Hz_image = czt(Hzi*z_phase*filter,kx,ky,x,y);
E2_image = abs(Ex_image)^2+abs(Ey_image)^2+abs(Ez_image)^2;
H2_image = abs(Hx_image)^2+abs(Hy_image)^2+abs(Hz_image)^2;

image(x*1e6,y*1e6,E2_image,"x (microns)","y (microns)","|E|^2 at z="+num2str(z_project*1e6)+"microns, " + num2str(lambda(i)*1e9) + "nm");
image(x*1e6,y*1e6,H2_image,"x (microns)","y (microns)","|H|^2 at z="+num2str(z_project*1e6)+"microns, " + num2str(lambda(i)*1e9) + "nm");

# export results to an fld file to be used as a source
Ex(:,:,i) = Ex_image;
Ey(:,:,i) = Ey_image;
Ez(:,:,i) = Ez_image;
Hx(:,:,i) = Hx_image;
Hy(:,:,i) = Hy_image;
Hz(:,:,i) = Hz_image;
}


# package the data for exporting
EM = rectilineardataset("EM fields",x,y,z);
EM.addparameter("lambda",c/f,"f",f);   
EM.addattribute("E",Ex,Ey,Ez);
EM.addattribute("H",Hx,Hy,Hz);
  
# save dataset to .mat file which can be imported at a later time
matlabsave("EM.mat",EM);