#DESCRIPTION:
#Implements small signal laser model and calculates the RF response and 3dB bandwidth
#based on DC parameters and results from DC simulation.
#
#INPUTS:
#ssTime - time at which to extract DC results (when simulation reached steady state)
#maxFreq - maximum frequency for the RF response
#numPoints - number of frequency points in the RF response

function ss_bandwidth(ssTime,maxFreq,numPoints){
    
    freq = linspace(0,maxFreq,numPoints);
    wfreq = 2*pi*freq;
    
    l = getnamed("TWLM_1","length");
    t = getnamed("TWLM_1","active region thickness");
    w = getnamed("TWLM_1","active region width");
    Gamma = getnamed("TWLM_1","mode confinement factor 1");
    beta = getnamed("TWLM_1","spontaneous emission factor 1");
    vg = c/getnamed("TWLM_1","group index 1");
    g0 = getnamed("TWLM_1","gain coefficient"); #differential gain
    chden_tr = getnamed("TWLM_1","carrier density at transparency");
    Vw = getnamed("TWLM_1","total well thickness");
    eps = getnamed("TWLM_1","gain compression factor");
    taus = 1/getnamed("TWLM_1","well carrier capture rate");
    taue = 1/getnamed("TWLM_1","well carrier escape rate");
    
    pden = getresult("TWLM_1","photon density profile");
    position = pden.getparameter("position");
    time = pden.getparameter("time");
    ssTimeIndex = find(time,ssTime);
    pden = pden.getattribute("photon density (1/m)");
    pden = pden(ssTimeIndex,:); #time index when steady state reached
    pden_avg = (integrate(pden,2,position))/l; #averaged over length
    S0 = pden_avg/t/w;
    
    chden = getresult("TWLM_1","carrier density profile");
    position = chden.getparameter("position");
    chden = chden.getattribute("carrier density (1/m^3)");
    chden = chden(ssTimeIndex,:); #time index when steady state reached
    Rsp = getresult("TWLM_1","spontaneous emission recombination rate profile");
    Rsp = Rsp.getattribute("spontaneous emission rate (1/s)");
    Rsp = Rsp(ssTimeIndex,:);
    Rsp_avg = (integrate(Rsp*chden,2,position))/l; #averaged over length
    chden_avg = (integrate(chden,2,position))/l; #averaged over length
    taup = 1/(Gamma*beta*Rsp_avg/S0 + Gamma*vg*g0*(chden_avg-chden_tr)/(1+eps*S0));
    
    Rsp_avg = (integrate(Rsp,2,position))/l; #averaged over length
    Rnr = getresult("TWLM_1","non-radiative recombination rate profile");
    Rnr = Rnr.getattribute("non-radiative recombination rate (1/s)");
    Rnr = Rnr(ssTimeIndex,:);
    Rnr_avg = (integrate(Rnr,2,position))/l; #averaged over length
    taun = 1/(Rnr_avg+(1-beta)*Rsp_avg);
    
    chi = 1 + taus/taue;
    
    A0 = vg*g0*S0/taup*(1+eps/vg/g0/taun);
    A1 = vg*g0*S0*(1+taus/taup)+eps*S0/taup*(1+taus/taue+taus/taun)+1/taun*(1+eps*S0);
    A2 = (1+eps*S0)*(1+taus/taue+taus/taun)+vg*g0*S0*taus+eps*S0*taus/taup;
    A3 = taus*(1+eps*S0);
    
    RF = Gamma*vg*g0*S0/e/Vw/(A0+1i*A1*wfreq-A2*wfreq^2-1i*A3*wfreq^3);
    RF = RF/RF(1);
    RF = 10*log10(abs(RF));
    
    output = struct;
    output.bandwidth = freq(find(RF,-3));
    output.rf_response = RF;
    output.freq = freq;
    return output;
}