runanalysis;

# get power by taking differnce of transmission function
T = getresult("transmission","T");
R = getresult("reflection","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");



f=Pabs_standard.f;
lf=length(f);
layer2 = getnamed("layer2","material");
layer3 = getnamed("layer3","material");

n = matrix(3,lf);
n(1,1:lf) = 1+matrix(lf);
n(2,1:lf) = getfdtdindex(layer2,f,min(f),max(f));
n(3,1:lf) = getfdtdindex(layer3,f,min(f),max(f));
d = [1;0.1e-6;1];
theta = 0;
RT_theory = stackrt(n,d,f,theta);


plot(c/f*1e6,1-RT_theory.Rs-RT_theory.Ts,Pabs_advanced.Pabs_total,Pabs_standard.Pabs_total,1-T.T+R.T);
legend("P_abs_theory","Pabs_advanced","Pabs_standard","1-T-R");