################################################
# Filename: usr_calc_absorption_old.lsf
#
# This script will integrate div P to 
# measure loss in a volume.  This result
# is compared to the loss as calculated with a 
# box of 2D monitors.
#
# Copyright 2007, Lumerical Solutions, Inc.
################################################

mname="mon3D";  	# monitor name
f_point=1;  		# select first frequency point;
f = getdata(mname,"f");
f=f(f_point);		


# Measure total absorbed power by using a
# box of 2D power monitors
Power_absorbed1 = transmission("x2") - transmission("x1") +
               transmission("y2") - transmission("y1") +
               transmission("z2") - transmission("z1") ;

Power_absorbed1=Power_absorbed1(f_point);
?"Total Absorbed Power (2D method): "+num2str(-100*Power_absorbed1)+ " %.";



# Measure total absorbed power by calculating 
# the divergence of P and integrating.
Px=getdata(mname,"Px");
Px=pinch(Px,4,f_point);               # select first frequency point
Py=getdata(mname,"Py");
Py=pinch(Py,4,f_point);
Pz=getdata(mname,"Pz");
Pz=pinch(Pz,4,f_point);
x=getdata(mname,"x");
y=getdata(mname,"y");
z=getdata(mname,"z");
n = size(Px);

# calculate dx.  Note that x2,y2,z2 are not at same locations as x,y,z
dx = x(2:n(1)) - x(1:(n(1)-1));
x2 = 0.5*( x(2:n(1)) + x(1:(n(1)-1)));
dy = y(2:n(2)) - y(1:(n(2)-1));
y2 = 0.5*( y(2:n(2)) + y(1:(n(2)-1)));
dz = z(2:n(3)) - z(1:(n(3)-1));
z2 = 0.5*( z(2:n(3)) + z(1:(n(3)-1)));


# calculate div P
?"dPx_dx"; 
temp = Px(2:n(1),1:n(2),1:n(3))- Px(1:(n(1)-1),1:n(2),1:n(3));  # calculate x component of divergence
temp= temp / meshgrid3dx(dx,y,z);  # divide by dx
dPx_dx=interp(temp,x2,y,z,x,y,z);  # interpolate back to standard grid

?"dPy_dy";
temp = Py(1:n(1),2:n(2),1:n(3))- Py(1:n(1),1:(n(2)-1),1:n(3));
temp= temp/ meshgrid3dy(x,dy,z);
dPy_dy=interp(temp,x,y2,z,x,y,z);

?"dPz_dz";
temp = Pz(1:n(1),1:n(2),2:n(3))- Pz(1:n(1),1:n(2),1:(n(3)-1));
temp= temp/ meshgrid3dz(x,y,dz);
dPz_dz=interp(temp,x,y,z2,x,y,z);



# sum each component, convert to units of loss, and 
# normalize to source power.
loss= -0.5 * real(dPx_dx + dPy_dy + dPz_dz) / sourcepower(f);
loss=pinch(loss);

# calculate total loss in box
Power_absorbed2= 100*integrate(loss,1:3,x,y,z);
?"Total Absorbed Power (3D method): "+num2str(Power_absorbed2)+ " %.";


# # Create integration filter in particle by analytic formula
# rad=0.2e-6;
# x0=0;
# y0=0;
# z0=0;
# X=meshgrid3dx(x,y,z);
# Y=meshgrid3dy(x,y,z);
# Z=meshgrid3dz(x,y,z);
#
# filter= ((X-x0)^2+(Y-y0)^2+(Z-z0)^2)<=rad^2;

# create integration filter in particle from index monitor
nx=getdata("index3d","index_x");
nx=pinch(real(nx));
filter = nx<1.4;   # set filter to 1 for all values of real n < 1.4

# integrate loss over filter region
Power_absorbed_sphere=100*integrate(loss*filter,1:3,x,y,z);
?"Loss in sphere: "+num2str(Power_absorbed_sphere)+ " %.";


image(x*1e9,y*1e9,loss(1:n(1),1:n(2),n(3)/2),"x (nm)","y (nm)","loss per unit volume");
image(x*1e9,y*1e9,filter(1:n(1),1:n(2),n(3)/2),"x (nm)","y (nm)","integration filter");