################################################
# Filename: usr_absorption.lsf
# Description: This file shows how to calculate
#              loss per unit volume with a set
#              of power and index monitors.  
#              It also shows how to integrate the
#              loss function, and compares the 
#              results with an alternate loss 
#              calculation technique.  
# Copyright 2009 Lumerical Solutions, Inc.
################################################

clear;
closeall;

# get monitor data

runanalysis("Pabs");
E2  =getresult("Pabs::field","E");
n   =getresult("Pabs::index","index");
Pabs=getresult("Pabs","Pabs");


# send datasets to visualizer
#visualize(Pabs);
#visualize(E2);
#visualize(n);


################################################
# image a cross section of the field 
# and index data at y=0, second frequency point
image(E2.x*1e6,E2.z*1e6,pinch(pinch(E2.E2     ,4,2),2,find(E2.y,0)),"x (um)","z (um)","E^2");
image(n.x*1e6,n.z*1e6,pinch(pinch(real(n.index_x),4,2),2,find(n.y,0)),"x (um)","z (um)","real(n)");


# plot cross section of absorption at y=0, second frequency point
image(Pabs.x*1e6,Pabs.z*1e6,pinch(pinch(Pabs.Pabs,4,2),2,find(Pabs.y,0)),"x (um)","z (um)","Absorption");


# get total absorbed power as a function of frequency
Pabs_total = getresult("Pabs","Pabs_total");
?"Wavelength (nm):                       "+num2str(transpose(c/Pabs_total.f)*1e9);  
?"Integrate absorption technique (%):      "+num2str((Pabs_total.Pabs_total)*100);



################################################
# check accuracy of total power absorption by comparing with
# a box of transmission monitors.  Assume no power flows out the sides,
# so only two monitors are required.
T1=-transmission("T1");
T2=-transmission("T2");
?"Net transmission technique (%):  "+num2str(transpose(T1-T2)*100);


################################################
# calculate power absorbed in circular region centered at 
# z=-200nm, with radius 200nm, at lambda=500nm,
# by integrating the absorption over this spherical region.

f = Pabs.f;
x = Pabs.x;
y = Pabs.y;
z = Pabs.z;

# select frequency point from Pabs
fi=find(Pabs.f,c/500e-9);
Pabs_temp=pinch(Pabs.Pabs,4,fi);

# create integration filter
X=meshgrid3dx(x,y,z);
Y=meshgrid3dy(x,y,z);
Z=meshgrid3dz(x,y,z);
filter= (X^2 + Y^2 + (Z+200e-9)^2 ) < (200e-9)^2;

# image cross section of integration filter
image(x*1e6,z*1e6,pinch(filter,2,find(y,0)),"x (um)","z (um)","integration filter");

# integrate power
Pabs_sphere = integrate(Pabs_temp*filter,1:3,x,y,z);

?"";
?"Absorption in sphere at "+num2str(c/f(fi)*1e9)+" nm: "+num2str(transpose(Pabs_sphere)*100);