#########################################################
# Calculate the total transmission through monitor with 
# transmission function.  This can be used to check the 
# poynting vector integrals.
mon_name="transmission";
f=getdata(mon_name,"f");
T1=transmission(mon_name);
R1=-transmission("reflection");   #- sign due to direction of power flow



#########################################################
# Reproduce the results of the transmission function by 
# manually integrating the poynting vector.  Normalize this to the
# sourcepower.  
mon_name="transmission";
x=getdata(mon_name,"x");
y=getdata(mon_name,"y");
f=getdata(mon_name,"f");
Pz=getdata(mon_name,"Pz");

Pz=pinch(real(Pz));

T2 = 0.5*integrate(Pz,1:2,x,y);  # power in watts
T2 = T2/sourcepower(f);  # normalize to source power



#########################################################
# Plot T1 and T2 for comparison. They should be exactly the same.
plot(f/1e12,T1*100,T2*100,"frequency (THz)","Power");
legend("T from transmission function",
       "T by integrating Pz");

#########################################################
# Plot T,R, and absorption. 
plot(f/1e12,T1*100,R1*100,(1-T1-R1)*100,"frequency (THz)","Power");
legend("Transmission",
       "Reflection",
       "Absorption");