clear;

# This commented block explains how is the length of TWLM set so that it can be discretized into 
# an integer number of sections. The length of each section is determined by the light velocity and
# time step.
#
# ng = getnamed('TWLM_1','group index 1');
# fs = getnamed('::Root Element::','sample rate');
# numSections = 10;
# len = c/ng/fs*numSections;

conf = getnamed('TWLM_1','mode confinement factor 1');
len = getnamed('TWLM_1','length');

############################ real part of complex gain from ONA impulse response transmission
Tr = getsweepresult("sweep","transmission");
f = Tr.getparameter("frequency")/1e12;
cden = Tr.getparameter("cden_init");
Tr = Tr.getattribute("TE transmission");

numcurves = length(cden);

g = log((abs(Tr))^2)/conf/len;
############################


############################ original complex gain directly from TWLM diagnostic output
gain_diag = getresult("TWLM_1","diagnostic/gain/fitted spectrum");
fff = gain_diag.getparameter("frequency")/1e12;

ggg = zeros(length(fff),numcurves);
for (i=1:numcurves) {
    ggg(:,i) = gain_diag.getattribute(i);
}
gggleg = splitstring(gain_diag.getattribute,endl);
############################


############### Plot comparison of real part of TWLM gain and real part of gain from ONA transmission
myleg=cell(2*numcurves);
for (i=1:numcurves) {
    plot(fff,real(ggg(:,i)),"frequency [THz]", "real(g) [1/m]"," ");
    holdon;
    myleg{i} = "TWLM "+ gggleg{i};
}

for (i=1:numcurves){
    plot(f(1:40:end),g(1:40:end,i),"frequency [THz]", "real(g) [1/m]"," ","pen=--");
    holdon;
    myleg{i+numcurves} = "ONA "+num2str(cden(i));
}

legend(myleg);
holdoff;

#expand x-axis to accomodate legend
xmin = fff(1);
xmax = fff(1) + 2*( fff(length(fff)) - fff(1));
setplot("x min",xmin);
setplot("x max",xmax);
setplot("legend position",6);
##########################################