################################################
# calculate power absorbed in  region where  
# the refractive indes is that of a particular materials
# (Silicon and Silver) for all frequencies,
# by integrating the absorption over this material.
#
# This file should work for both 2D and 3D simulation
#################################################
clear;
closeall;

# get monitor data
runanalysis;
m="pabs_adv";
Pabs=getresult(m,"Pabs");
n = getresult(m+"::index","index");
n = n.index_x;
f=Pabs.f;
nx=length(Pabs.x);
ny=length(Pabs.y);
nz=length(Pabs.z);
nf=length(Pabs.f);

# define material filters

filter_mat_Si = "Si (Silicon) - Palik";
filter_mat_Ag = "Ag (Silver) - CRC";

filter_real_Si = matrix(nx,ny,nz,nf);
filter_imag_Si = matrix(nx,ny,nz,nf);
filter_real_Ag = matrix(nx,ny,nz,nf);
filter_imag_Ag = matrix(nx,ny,nz,nf);

n_filter_real_Si = real(getfdtdindex(filter_mat_Si,f,min(f),max(f)));
n_filter_imag_Si = imag(getfdtdindex(filter_mat_Si,f,min(f),max(f)));
n_filter_real_Ag = real(getfdtdindex(filter_mat_Ag,f,min(f),max(f)));
n_filter_imag_Ag = imag(getfdtdindex(filter_mat_Ag,f,min(f),max(f)));

#relative difference for almosequal comparator
rel_diff = 1e-15;



for (i=1:nf) {
    filter_real_Si(1:nx,1:ny,1:nz,i) =  almostequal(n_filter_real_Si(2),pinch(real(n),4,2),rel_diff);
    filter_imag_Si(1:nx,1:ny,1:nz,i) =  almostequal(n_filter_imag_Si(2),pinch(imag(n),4,2),rel_diff);
}
filter_Si=filter_real_Si*filter_imag_Si;




for (i=1:nf) {
    filter_real_Ag(1:nx,1:ny,1:nz,i) =  almostequal(n_filter_real_Ag(2),pinch(real(n),4,2),rel_diff);
    filter_imag_Ag(1:nx,1:ny,1:nz,i) =  almostequal(n_filter_imag_Ag(2),pinch(imag(n),4,2),rel_diff);
}
filter_Ag=filter_real_Ag*filter_imag_Ag;

#integrate

Pabs_Si = integrate2(Pabs.Pabs*filter_Si,1:3,Pabs.x,Pabs.y,Pabs.z);
Pabs_Ag = integrate2(Pabs.Pabs*filter_Ag,1:3,Pabs.x,Pabs.y,Pabs.z);

#Plots of filters and total absorption at a wavelength
if (nz>1){ # 3D
image(Pabs.x*1e6,Pabs.y*1e6,pinch(pinch(filter_Si,4,8),3,20),"x (um)","y (um)","filter_Silicon");
image(Pabs.x*1e6,Pabs.y*1e6,pinch(pinch(filter_Ag,4,8),3,20),"x (um)","y (um)","filter_Silver");
image(Pabs.x*1e6,Pabs.y*1e6,pinch(pinch(log10(Pabs.Pabs),4,8),3,20),"x (um)","y (um)","Total Absorption");
}else{ # 2D
image(Pabs.x*1e6,Pabs.y*1e6,pinch(filter_Si,4,8),"x (um)","y (um)","filter_Silicon");
image(Pabs.x*1e6,Pabs.y*1e6,pinch(filter_Ag,4,8),"x (um)","y (um)","filter_Silver");
image(Pabs.x*1e6,Pabs.y*1e6,pinch(log10(Pabs.Pabs),4,8),"x (um)","y (um)","Total Absorption");
}

# plots of absorption spectrum
plot(c/f*1e9,Pabs_Si,"wavelength","Power Absorption","Absorption in Silicon");
plot(c/f*1e9,Pabs_Ag,"wavelength","Power Absorption","Absorption in Silver");
plot(c/f*1e9,Pabs_Si+Pabs_Ag,"wavelength","Power Absorption","Total absorption");