# Michael Kelzenberg
# LUMERICAL script for calculating optical absorption of Si wire array after
# completion of FDTD simulation.  Reports results with and without partial
# spectral averaging.

clear;

#The "datasetname" is a descriptor for each (sub)simulation, e.g., the wavelength.
setname = "Test";

#The "basename" is a unique descriptor for the simulation structure
basename = "Example";

matlab_mat_file=basename + "_" + setname + ".mat";
logfile=basename + "_" + setname + ".txt";

# This is the incident power density (my simulations always use 1 mW/cm2)
incdPwr = 10;  # W / m2

#This is the area of the simulation unit cell
myArea = (7e-6)^2;
#Note: Lumerical, aware of the symmetric BC's for the simulations, seems to 
# "unfold" the simulation data so that monitor readings correspond to the
# full unit cell, not just the 1/4 cell that was actually simulated.

#Get grid coordinates of the simulation:
Pabs_x=getdata("VolumeMonitor","x");
Pabs_y=getdata("VolumeMonitor","y");
Pabs_z=getdata("VolumeMonitor","z");

#Get simulation frequency (i.e., the wavelength)
freq=getdata("VolumeMonitor","f");

#Get the refractive index throughout the simulation volume.
n_x = getdata("RefrMonitor","x");
n_z = getdata("RefrMonitor","z");
n_y = getdata("RefrMonitor","y");
# I use the refractive index data to "mask" the c-Si volume for numerical 
# integration later.  I believe I needed separate x, y, and z vectors because
# FDTD stores refractive index on a different (interlaced) grid than electric
# field data... 

# Number of grid points for electric field data...
nx=length(Pabs_x);
ny=length(Pabs_y);
nz=length(Pabs_z);

#Get the ... |E|^2 data (?) :
E2_pavg=getdata("VolumeMonitor","E2x_pavg") + getdata("VolumeMonitor","E2y_pavg") + getdata("VolumeMonitor","E2z_pavg");

#Get the refractive index data.
n=getdata("RefrMonitor","index_x");

#This code sets up the c-Si "mask":
n_i = (length(n_x) - nx)/2;
n = n( (n_i+1):(n_i+nx), (n_i+1):(n_i+nx), 1:nz,1);
nSi = n( nx/2, ny/2, nz/2);
isSi = n==nSi;

#This applies the mask:
n = n*isSi;

#This gets the non-pavg |E|^2 data:
E2=getelectric("VolumeMonitor");	

#Parameters for calculating optical generation
omega = 2*pi*freq;
epsilon = eps0*n^2;

#This calculates absorbed power.  
Pabs_pavg = 0.5*omega*E2_pavg*imag(epsilon);
Pabs = 0.5*omega*E2*imag(epsilon);

#Convert the power dissipation density to optical generation rate:
Ngen_pavg = Pabs_pavg*1e-6/(6.626e-34*freq); #per cm3 per sec
Ngen = Pabs*1e-6/(6.626e-34*freq);

#Integrate optical generation to calculate absorption (photo)current:
Current_pavg = 1e6*1.61e-19*integrate(Ngen_pavg, 1:3, Pabs_x, Pabs_y, Pabs_z); #A per wire
Current = 1e6*1.61e-19*integrate(Ngen, 1:3, Pabs_x, Pabs_y, Pabs_z);   # (10^6 to go from cm3 to m3) 

# This can be used to visualize the result...
# image(Pabs_x,Pabs_y,Ngen_pavg);

#This might be redundant with above integrations:
IntgPwr_pavg = integrate(Pabs_pavg, 1:3, Pabs_x,Pabs_y,Pabs_z );
IntgPwr = integrate(Pabs, 1:3, Pabs_x,Pabs_y,Pabs_z );

#Now, determine the power absorption based on the (plane) monitors:
toppwra = transmission_pavg("ReflMonitor");
Absfrac_pavg = 1 - toppwra;

toppwr = transmission("ReflMonitor");
Absfrac = 1 - toppwr;

#Usually, the "volume integration" (above) and "plane monitors" give slightly 
# different results as far as the mangitude of total absorption.  I make sure
# to record both here to confirm consistency of the simulation.  

write(logfile,"Processed data.");
write(logfile,"Monitors report:          absorption_pavg=" + num2str(Absfrac_pavg ) +                      "  ,    absorption=" + num2str(Absfrac ) );
write(logfile,"Monitors report:            absPower_pavg=" + num2str(Absfrac_pavg*incdPwr*myArea ) +       " W,  absPower=" + num2str(Absfrac*incdPwr*myArea) + " W" );
write(logfile,"Volume integration reports: absPower_pavg=" + num2str(IntgPwr_pavg) +                       " W,  absPower=" + num2str(IntgPwr) + " W" );
write(logfile,"Volume integration reports:  current_pavg=" + num2str(Current_pavg) +                       " A,  current=" + num2str(Current) + " A" );


write(logfile,"Saving output file: " + matlab_mat_file );

if( fileexists(matlab_mat_file) )
{
	del(matlab_mat_file );
}

matlabsave(matlab_mat_file,Pabs_x,Pabs_y,Pabs_z,Pabs,Pabs_pavg,Ngen_pavg,Ngen,freq,Absfrac_pavg,Absfrac,IntgPwr_pavg,IntgPwr,Current_pavg,Current,n);