clear; nonorm; # this example only works in the nonorm state

power_integration = 1; # use the "power" returned by the monitors
manual_poynting_integration = 0; # NOTE: needs a LOT more memories in the time monitor

if (power_integration) { # part 1
switchtolayout;
setnamed("time", "output Pz", 0);
setnamed("frequency", "output Pz", 0);
run;

#################################
# power integration - time domain
t = getdata("time","t"); # [Second]
power_time = getdata("time","power"); # [Watt]
energy_time = integrate(power_time,2,t); # [Joule]
plot(t,power_time,"time(s)","power (W)","instantaneous power - time domain");

################################# 
# Energy spectrum integration - frequency domain
f = getdata("frequency","f"); # [Hz]

# time averaged power(w) = 1/2*integral(dS ExH*)
# multiplied by 2 for a pulse to remove the time averaged factor (0.5)
energy_spectrum_f=2*real(getdata("frequency","power")); # [Watt/Hz^2] in nonorm state

# a factor of 2 coming from integraqting the negative frequencies. 
# NOTE that this factor of 2 is not needed if complex fields are used in the simulation
energy_f = integrate(energy_spectrum_f,1,f)*2; # [Watt/Hz], which is equivalently Watt*s=Joule

##################################
# sourcepower command verification - frequency domain
# a factor of 2 from negative frequencies, another factor of 2 to remove the time averaged factor (1/2) that we defined
energy_f_source = 2*2*integrate( sourcepower(f), 1, 2*pi*f) / (2*pi); # [Watt/Hz], which is equivalently Watt*s=Joule


plot(f,energy_spectrum_f,2*sourcepower(f),"frequency (Hz)","Energy per unit freuquency (Joule/Hz)","Energy spectrum - frequency domain");
legend("energy spectrum","2*sourcepower");

?"Ratio of energy_time to energy_f = " + num2str(energy_time/energy_f);
?"Ratio of energy_time to energy_f_sourcepower = " + num2str(energy_time/energy_f_source);
}



# Poynting vector integration. WARNING: it needs more memories to record data in every time step
#######################################################################################
if (manual_poynting_integration) { # part 2
switchtolayout;
setnamed("time", "output Pz", 1);
setnamed("frequency", "output Pz", 1);
run;

# manual integration - time domain
t = getdata("time","t"); # [Second]
x = getdata("time","x"); # [meter]
y = getdata("time","y"); # [meter]
Pz_time = getdata("time","Pz"); # [Watt/meter^2]
power_time_manual = integrate(Pz_time, 1:2, x,y); # [Watt]
energy_time_manual = integrate(power_time_manual, 2, t); # [Joule]

# manual integration - frequency domain
f = getdata("frequency","f"); # [Hz]
Pz_f =getdata("frequency","Pz"); # [Watt/meter^2/Hz^2]
energy_spectrum_f_manual = integrate(real(Pz_f), 1:2, x,y); # [Watt/Hz^2]
# a factor of 2 coming from integraqting the negative frequencies. 
energy_f_manual = integrate(energy_spectrum_f_manual, 2, f)*2; # [Watt/Hz], which is equivalently Watt*s=Joule.

?"Ratio of energy_time_manual to energy_f_manual = " + num2str(energy_time_manual/energy_f_manual);
}