#This script reproduces results from table 2 in 
#1997, Chan, Lowery, Semiconductor laser linewidth from the transmission-line laser model
#
#Device parameters are set using GUI (icp project file with the ame name), except grating coupling coefficient
#which is set using script below

clear;

#get/set some device parameters needed for calculation												  
alphai = getnamed('TWLM_1','loss 1')/(10*log10(exp(1))); #waveguide loss
ng = getnamed('TWLM_1','group index 1'); #group index
vg = c/ng; #group velocity
nsp = 2; #population inversion factor
alphaeff = getnamed('TWLM_1','linewidth enhancement factor'); #material linewidth enhancement factor
len = getnamed('TWLM_1','length');
thick = getnamed('TWLM_1','active region thickness');
width = getnamed('TWLM_1','active region width');
k = 2.0/len; #grating coupling coefficient from reference
switchtolayout;
setnamed('TWLM_1','grating coupling real coefficient 1',k);
kl = 1.3; #linewidth correction factor due to longitudinal field enhancement in QWS DFB from reference

Pout = 1/[0.1,0.16,0.23,0.32,0.45,1.06]*1e-3; #power levels from reference [W]
curr = [0.0967187,0.0715625,0.0589844,0.0507813,0.0442188,0.0349219]; #precalculated injection current [A] corresponding to power levels Pout
gth = 10623.7; #precalculated gain at threshold

df1 = zeros(length(curr)); #linewidth from Method 1
df2 = zeros(length(curr)); #linewidth from Method 2
Pfacet = zeros(length(curr)); #facet power
rsp1 = zeros(length(curr)); #total cavity spontaneous recombination rate from Method 1
rsp2 = zeros(length(curr)); #total cavity spontaneous recombination rate from Method 2

#Loop over facet power and calculate linewidths
for(i=1:length(curr)){
    switchtolayout;
    setnamed('DC_1','amplitude',curr(i));
    run; #run laser simulation for target facet power
    
    #Find lasing frequency from spectrum
    spectr = getresult('OSA_1','sum/signal');
    f = c/spectr.getparameter('wavelength');
    spectr = spectr.getattribute('power (dBm)');
    f = f(find(spectr,max(spectr)));
    
    #Calculate total number of photons in cavity
    Snum_dataset = getresult('TWLM_1','photon density profile');
    t = Snum_dataset.getparameter('time');
    starttime = getnamed('OPWM_1','start time');
    x = Snum_dataset.getparameter('position');
    dx = x(2)-x(1);
    Snum_dataset = Snum_dataset.getattribute('photon density (1/m)');
    Snum_dataset = sum(Snum_dataset(find(t>starttime),:),2)*dx;
    Snum = mean(Snum_dataset);
    
    #Find facet power
    P1 = getresult('OPWM_1','sum/power');
    P2 = getresult('OPWM_2','sum/power');
    Pfacet(i) = mean([P1,P2]);
    
    #Linewidth from Method 1
    alpham = (P1+P2)/(vg*h*f*Snum);
    rsp1(i) = vg*(alphai+alpham)*nsp;    
    df1(i) = kl*rsp1(i)/(4*pi*Snum)*(1+alphaeff^2);
    
    #Linewidth from Method 2
    rsp2(i) = vg*gth*nsp;
    df2(i) = kl*rsp2(i)/(4*pi*Snum)*(1+alphaeff^2);
}

#Plot and compare 3 linewidth methods with reference (Table 2)
df_ref_low = [11,18,23,32,45,105];
df_ref_high = [15,22,30,40,55,115];
plot(Pfacet,df1*1e-6,'P facet [W]','linewidth [MHz]','','plot type = point, marker style = x, marker size=10');
holdon;
plot(Pfacet,df2*1e-6,'P facet [W]','linewidth [MHz]','','plot type = point, marker style = +, marker size=10');
plot(Pout,df_ref_low,'P facet [W]','linewidth [MHz]','','plot type = point, marker style = o, marker size=6');
plot(Pout,df_ref_high,'P facet [W]','linewidth [MHz]','','plot type = point, marker style = o, marker size=6');
holdoff;