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");