clear; closeall;
run;

####################################
# Standard fourier transform (SFT)

# R vs f
f=getdata("R_SFT","f");
R=-transmission("R_SFT");

plot(f/1e12,R,"f (THz)","reflection","Standard Fourier Transform (SFT)");
legend("SFT");



####################################
# Partial spectral average (PSA)

# SFT
f  = getdata("R_SFT","f");
w  = transpose(2*pi*f);
nf = length(f);
R  = -transmission("R_SFT");    # normalized reflected power
R2 = R*sourcepower(f);          # remove source power norm.  R2 will have units of watts, rather than being unitless
sn = abs(sourcenorm(f))^2;      # source_normalization (the weighting function)
source_power = sourcepower(f);  # source power 


# PSA
f_psa  = getdata("R_PSA","f");
w_psa  = transpose(2*pi*f_psa);
R_psa1 = -transmission_pavg("R_PSA" );
R_psa1 = interp(R_psa1,f_psa,f); #interpolate for easier comparison with SFT
delta  = getdata("R_PSA","delta");


# calculate PSA from SFT, at the same frequencies as the PSA monitor
R_psa2 = matrix(length(f_psa));
for (i=1:length(f_psa)) { 
  h2=delta / ( (w - w_psa(i))^2 + (pi*delta)^2 );

  # manually calculate the partial spectral average.  It is the 
  # integral of the reflection weighted by the source spectrum and lorentz function divided by the 
  # integral of the sourcepower wieghted by the source spectrum and lorentz function.
  R_psa2(i) = integrate( sn*h2*R2            ,1,w) /
              integrate( sn*h2*source_power  ,1,w) ;
}
R_psa2 = interp(R_psa2,f_psa,f); #interpolate for easier comparison with SFT



# Compare results at specific frequency.
fi=530e12;      	# target frequency
fi=f_psa( find(f_psa,fi) );
i=find(f,fi);
fi=f(i);

?"Partial Spectral Average (PSA) at "+num2str(round(fi/1e12))+"THz.";
?"PSA from STF: "+num2str(R_psa2(i));
?"PSA             : "+num2str(R_psa1(i));
?"";

# h2 for plot
h2=delta / ( (w - 2*pi*fi)^2 + (pi*delta)^2 );

# Plot PSA
plot(f/1e12,R,R_psa2,R_psa1,h2/max(h2),sn/max(sn),
     "f (THz)", "reflection","Partial Spectral Average (PSA)");
legend("SFT","PSA from SFT","PSA","Lorentz" ,"Source spectrum");




####################################
# Total spectral average (TSA)

# Calculate R_avg by weighted average of standard fourier transform
f=getdata("R_SFT","f");
w = 2*pi*f;
nf=length(f);
R=-transmission("R_SFT");       # normalized reflected power
R2 = R*sourcepower(f);          # remove source power norm.  R2 will have units of watts, rather than being unitless
sn = abs(sourcenorm(f))^2;      # source_normalization (the weighting function)
source_power = sourcepower(f);  # source power 

# manually calculate the total spectral average.  It is 
# the integral of the reflection weighted by the source spectrum divided by 
# the integral of the sourcepower wieghted by the source spectrum.
R_avg_SFT = integrate( sn*R2          ,1,w) /
            integrate( sn*source_power,1,w) ;


# Calculate R_avg with total spectral average monitor
R_avg_TSA =-transmission_avg("R_TSA");

?"Total Spectral Average (TSA).";
?"avg(SFT): " + num2str(R_avg_SFT);
?"TSA     : " + num2str(R_avg_TSA);

plot(f/1e12,R,matrix(nf)+R_avg_SFT,matrix(nf)+R_avg_TSA,sn/max(sn),
  "f (THz)", "reflection","Total Spectral Average (TSA)");
legend("SFT","avg(SFT)","TSA","Source spectrum");