#This script reproduces results from table 1 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.6; #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 = [1.05,1.25,1.35,1.6,2.15]/len; #array of grating coupling coefficients from reference
ith = [0.038,0.0362,0.0348,0.0345,0.033]; #precalculated threshold current for each coupling coefficient [A]
gth = [9075.21,7954.15,7053.27,6858.21,5847.71]; #precalculated threshold gain for each coupling coefficient [1/m]
curr = [0.0609375,0.060625,0.06125,0.061875,0.06375]; #precalculated current for 5 mW facet power for each coupling coefficient
kl = [3,2.7,2.5,1.8,1.3]; #linewidth correction factor due to longitudinal field enhancement in QWS DFB from reference for each coupling coefficient

#Initialize variables for caching calculated values
Snum = zeros(length(k)); #Cavity number of photons
alpham = zeros(length(k)); #DFB effective mirror loss
Pfacet = zeros(length(k)); #facet power
rsp1 = zeros(length(k)); #total cavity spontaneous recombination rate for Method 1
rsp2 = zeros(length(k)); #total cavity spontaneous recombination rate for Method 2
df1 = zeros(length(k)); #linewidth using spontaneous emission rate from Method 1
df2 = zeros(length(k)); #linewidth using spontaneous emission rate from Method 2

#Loop over grating coupling coefficients and calculate linewidth for each at 5 mW facet power
for(i=1:length(k)){
    switchtolayout;
    setnamed('TWLM_1','grating coupling real coefficient 1',k(i));
    setnamed('DC_1','amplitude',curr(i)); #set current for 5 mW facet power
    run; #run laser simulation for 5 mW 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(i) = 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(i) = (P1+P2)/(vg*h*f*Snum(i));
    rsp1(i) = vg*(alphai+alpham(i))*nsp;
    df1(i) = kl(i)*rsp1(i)/(4*pi*Snum(i))*(1+alphaeff^2);

    #Linewidth from Method 2
    rsp2(i) = vg*gth(i)*nsp;
    df2(i) = kl(i)*rsp2(i)/(4*pi*Snum(i))*(1+alphaeff^2);
}

#Plot and compare linewidth with reference (Table 1)
#Also compare photon number and mirror loss with reference (Table 1)
Snum_ref = [213603,229972,250981,278798,331159];
alpham_ref = [4556,4140,3892,3500,2898]; #1/m
df_ref = [49.6,39.5,31.9,20,11.1]; #MHz
plot(k*len,Snum,'coupling strength','Photon number','','plot type = point, marker style = x, marker size=6');
holdon;
plot(k*len,Snum_ref,'coupling strength','Photon number','','plot type = point, marker style = o, marker size=6');
holdoff;

plot(k*len,alpham,'coupling strength','mirror loss [1/m]','','plot type = point, marker style = x, marker size=6');
holdon;
plot(k*len,alpham_ref,'coupling strength','mirros loss [1/m]','','plot type = point, marker style = o, marker size=6');
holdoff;

plot(k*len,df1*1e-6,'coupling strength','linewidth [MHz]','','plot type = point, marker style = x, marker size=6');
holdon;
plot(k*len,df2*1e-6,'coupling strength','linewidth [MHz]','','plot type = point, marker style = x, marker size=6');
plot(k*len,df_ref,'coupling strength','linewidth [MHz]','','plot type = point, marker style = o, marker size=6');
holdoff;