# run analysis
#runanalysis;

# get power with transmission box
Pabs_box=getresult("trans_box","T");

# get power by spatially integrating the absorption
# use standard technique
m="pabs";
Pabs_standard=getresult(m,"Pabs_total");

# get power by spatially integrating the absorption
# use advanced technique
m="pabs_adv";
Pabs_advanced=getresult(m,"Pabs_total");


# compare with analytic result from the text file 'mie_au_jc_FDTD.txt'
f=Pabs_box.f;
a = readdata("mie_au_jc_FDTD.txt");
r=50e-9;
size_parameter = 2*pi*r*f/c;
sizeparam = pinch(a,2,1);
qabs = pinch(a,2,3);
Qabs_theory= interp(qabs,sizeparam,size_parameter);
Qabs_theory = Qabs_theory*r^2*pi*sourceintensity(f)/sourcepower(f);  # normalize data


# plot absorption spectrum
plot(c/f*1e6,-Pabs_box.T,Pabs_standard.Pabs_total,Pabs_advanced.Pabs_total,Qabs_theory,    
     "wavelength (um)","Power","Absorption");
legend("Net power into box","Standard","Advanced","Analytic");

# plot one cross section of the absorption profile
m="pabs";
Pabs=getresult(m,"Pabs");
image(Pabs.x*1e9,Pabs.z*1e9,
      pinch(pinch(Pabs.Pabs,4,18),2,find(Pabs.y,0)),
      "x (nm)","z (nm)","Absorption profile");