### INPUT

filename = "Piprek2000OQE";

### END INPUT

### START FUNCTIONS

function spontaneousrecombrate(gsp,Eph,nr){

prefactor = nr^2*Eph^2/pi^2/(h/2/pi/e)^3/c^2;
rate = prefactor*gsp;

result = integrate(rate,2,Eph);
return result;
}

function srh(taun,taup,n,p,ni,Ei,T){
    kbT = 26e-3*T/300;
    n1 = ni*exp(Ei/kbT);
    p1 = ni*exp(-Ei/kbT);
   
    return (n*p-ni^2)/(taup*(n+n1)+taun*(p+p1));
}

function auger(caun,caup,n,p,ni){
    return (caun*n+caup*p)*(n*p-ni^2);
}


function ni_GaInAsP(y,T,Eg){
    #http://www.ioffe.ru/SVA/NSM/Semicond/GaInAsP/bandstr.html#carrier%20concentration 
    #(GaInAsP lattice matched to InP except for Eg which is input)
    
    Nc = 4.82e15*(0.08-0.039*y)^(3.0/2.0)*T^(3.0/2.0)*1e6; #m^{-3}
    Nv = 4.82e15*(0.6-0.18*y)^(3.0/2.0)*T^(3.0/2.0)*1e6; #m^{-3}
    
    kT = 8.617e-5*T; #eV
    ni = sqrt(Nc*Nv)*exp(-Eg/2/kT);
    return ni;
}

function savetwlmgain(filename,cden,freq,gain){

format long;

gainSortMap = sortmap(cden,false);
cden=cden(gainSortMap);
for(j=1:size(gain,2)){
    gain(1:size(gain,1),j) = gain(gainSortMap,j);
}

numcurves = size(gain,1);
numfreq = size(gain,2);

write(filename,'(1,'+num2str(numcurves)+')');

str = '';
for(n=1:numcurves){
    str = str + num2str(cden(n)) + ' ';
}
write(filename,str);

write(filename,'(' + num2str(numfreq) + ',' + num2str(numcurves+1) + ')');

str = '';
for(m=1:numfreq){
    str = str + num2str(freq(m)) + ' ';
    for (n=1:numcurves){
        str = str + num2str(gain(n,m)) + ' ';
    }
    str = str + endl;
}
write(filename,str);

}

function saverecombinationrate(filename,cden,radRecomb,nonRadRecomb){

format long;

numcurves = length(radRecomb);

recombSortMap = sortmap(cden,false);
cden=cden(recombSortMap);
radRecomb = radRecomb(recombSortMap);
nonRadRecomb = nonRadRecomb(recombSortMap);

for(n=1:numcurves){
    write(filename,num2str(cden(n))+' '+num2str(radRecomb(n))+' '+num2str(nonRadRecomb(n)));
}

}

function plotData(projectname, temperature, datafileRoot){

mylegend=cell(2*length(temperature));
for (i=1:length(temperature)){
    
    load(projectname+"_T="+num2str(temperature(i))+".icp");
    x=getsweepdata("sweep","current");
    y=getsweepdata("sweep","power");  
    plot(x,y,"curent [A]","Power [W]","LI Curves");
    holdon;
    refdata = readdata(datafileRoot+"_LI_T="+num2str(temperature(i))+".txt");
    plot(refdata(:,1),refdata(:,2));
    mylegend{2*i-1}= "Lumerical "  +num2str(temperature(i))+" K";
    mylegend{2*i} =  "Reference 1 "+num2str(temperature(i))+" K";
        
}
holdoff;                            
legend(mylegend);

}

#### END FUNCTIONS

new;

if(!exist('fname')){
    fname = filename;
}
jsonload(fname);
jsonload(fname+"_opt");
jsonload(fname+"_mqw");

temperatureInd = find(mqwOut.temperatureVec==twlm.temperature);
Temp = mqwOut.temperatureVec(temperatureInd);  

freq_base = mqwOut.emission.frequency;
cden=mqwOut.cdenVec;
numcden = length(mqwOut.cdenVec);
gsp_TE_base=gsp_TM_base=gst_TE_base=gst_TM_base=
matrix(numcden,length(freq_base));
for(i=1:numcden){
     gsp_TE_base(i,:) = mqwOut.emission.te.spont(i,temperatureInd,:);  
     gsp_TM_base(i,:) = mqwOut.emission.tm.spont(i,temperatureInd,:);
     gst_TE_base(i,:) = mqwOut.emission.te.stim(i,temperatureInd,:);  
     gst_TM_base(i,:) = mqwOut.emission.te.stim(i,temperatureInd,:);
}

#Upsample Frequency
freq = zeros((length(freq_base)-1)*common.frequpsample+1);
for(i=1:length(freq_base)-1){
    df = (freq_base(i+1)-freq_base(i))/common.frequpsample;
    for(j=1:common.frequpsample){
        freq((i-1)*common.frequpsample+j) = freq_base(i)+df*(j-1);
    }
}
freq(end) = freq_base(end);
gsp_TE=gsp_TM=gst_TE=gst_TM=matrix(numcden,length(freq));
if(common.frequpsampletype==0){
    for(i=1:numcden){
        gsp_TE(i,:) = interp(gsp_TE_base(i,:),freq_base,freq);
        gsp_TM(i,:) = interp(gsp_TM_base(i,:),freq_base,freq);
        gst_TE(i,:) = interp(gst_TE_base(i,:),freq_base,freq);
        gst_TM(i,:) = interp(gst_TM_base(i,:),freq_base,freq);
    }
}
else if(common.frequpsampletype==1){
    for(i=1:numcden){
        gsp_TE(i,:) = spline(gsp_TE_base(i,:),freq_base,freq);
        gsp_TM(i,:) = spline(gsp_TM_base(i,:),freq_base,freq);
        gst_TE(i,:) = spline(gst_TE_base(i,:),freq_base,freq);
        gst_TM(i,:) = spline(gst_TM_base(i,:),freq_base,freq);
    }
}
else{
    for(i=1:numcden){
        gsp_TE(i,:) = chebin(transpose(gsp_TE_base(i,:)),freq_base,freq,freq_base(1),freq_base(end));
        gsp_TM(i,:) = chebin(transpose(gsp_TM_base(i,:)),freq_base,freq,freq_base(1),freq_base(end));
        gst_TE(i,:) = chebin(transpose(gst_TE_base(i,:)),freq_base,freq,freq_base(1),freq_base(end));
        gst_TM(i,:) = chebin(transpose(gst_TM_base(i,:)),freq_base,freq,freq_base(1),freq_base(end));
    }
}
   

# Calculate radiative (spontaneous) and non-radiative recombination rates
# Radiative recombination is just from 1 TE mode (call it first TE mode) within the twlm spectrum. 
# Radiative recombination of the other 2 modes (second TE and TM) as well as the recmbination from 
# the first TE mode falling oustide of the twlm spectrum must contribute to nonradiative recombination.
twlmSpectrumStart = twlm.centre_frequency(temperatureInd) - twlm.sample_rate/2;
twlmSpectrumEnd = twlm.centre_frequency(temperatureInd) + twlm.sample_rate/2;
radRecombinationRate = zeros(numcden);

for(i=1:numcden){
    
    # extract the in-INTC-simulation-BW TE spontaneous emission spectrum and photon energy arrays,resavgspec and resengy.
    # the 2/3 factor comes from the fact that there are 2 TE modes: one with the propagation direction  E field in each 
    # orthogonal direction parallel to the mqw and slab waveguide layers and one TM mode with E field and direction of
    # propagation normal to the mqw and slab waveguide layers
    resavgspec = ( 2*gsp_TE(i,find((freq>twlmSpectrumStart) & (freq<twlmSpectrumEnd)))/3 ); #take TE only
                      #+gsp_TM(i,find((freq>twlmSpectrumStart) & (freq<twlmSpectrumEnd))) )/3;
    
    resengy =  h*transpose( freq( find( (freq>twlmSpectrumStart) & (freq<twlmSpectrumEnd) ) ) )/e;                
    
    radRecombinationRate(i) = spontaneousrecombrate( resavgspec,resengy,real(wgSweep.neff(2))  );
}
radRecombinationRate = radRecombinationRate/cden; #1/s

tausrhn = twlm.tausrhn;
tausrhp = twlm.tausrhp;
ni = ni_GaInAsP(mqwOut.y_well,Temp,mqwOut.Egk0(:,temperatureInd));
Ei = twlm.Ei;
kb = 1.38064852e-23; # Boltzmann constant J/K
caun =twlm.C0augn*exp(twlm.Eaaugn/(kb*Temp/e))*(mqwOut.totalLength/mqwOut.qwLength)^2/twlm.qw_density_correction(temperatureInd)^3;
caup =twlm.C0augp*exp(twlm.Eaaugp/(kb*Temp/e))*(mqwOut.totalLength/mqwOut.qwLength)^2/twlm.qw_density_correction(temperatureInd)^3;
nonradRecombinationRate = srh(tausrhn,tausrhp,cden,cden,ni,Ei,Temp)
                        + auger(caun,caup,cden,cden,ni);
#add "nonradiative" contributions of spontaneous recombination, i.e., out of band TE + all of TM                        
for(i=1:numcden){
    nonradRecombinationRate(i) = nonradRecombinationRate(i) + 
        spontaneousrecombrate((2*gsp_TE(i,:)+gsp_TM(i,:))/3,h*transpose(freq)/e,real(wgSweep.neff(2)))-radRecombinationRate(i)*cden(i);
}
nonradRecombinationRate = nonradRecombinationRate/(cden); #1/s

# Calculate diffusion constant
mun = twlm.mun; #m2/(Vs)
mup = twlm.mup;
mu = mun*mup/(mun+mup);
diffConst = mu*kb*Temp/e; #m2/s

# Save data files to import to TWLM model
if(fileexists(common.name+'_gain.dat')){rm(common.name+'_gain.dat');}
savetwlmgain(common.name+'_gain.dat',mqwOut.cdenVec,freq,gst_TE);
if(fileexists(common.name+'_spem.dat')) {rm(common.name+'_spem.dat');}
savetwlmgain(common.name+'_spem.dat',mqwOut.cdenVec,freq,gsp_TE);
if(fileexists(common.name+'_rcmb.dat')){rm(common.name+'_rcmb.dat');}
saverecombinationrate(common.name+'_rcmb.dat',mqwOut.cdenVec,radRecombinationRate,nonradRecombinationRate);

twlmOut=struct;
twlmOut.twlm_template_filename = "twlm.icp";
twlmOut.twlm_parameter_filename=common.name+"_T="+num2str(Temp)+"_param_twlm";
twlmOut.twlm_run_filename=common.name+"_T="+num2str(Temp)+".icp";
twlmOut.twlm_result_filename=common.name+"_T="+num2str(Temp)+"_result_twlm";
twlmOut.frequency = twlm.centre_frequency(temperatureInd);
twlmOut.length=common.active_region_length;
twlmOut.active_region_width=common.active_region_width;
twlmOut.active_region_thickness=mqwOut.totalLength; #input
twlmOut.diffusion_constant=diffConst; #calculated
twlmOut.effective_index=wgSweep.neff(2); #calculated
twlmOut.group_index = wgSweep.ng(2); #calculated
twlmOut.loss = twlm.mode_loss(temperatureInd); 
twlmOut.temperature=mqwOut.temperatureVec;
twlmOut.facet_reflectivity_left=twlm.facet_reflectivity_left;  
twlmOut.facet_reflectivity_right=twlm.facet_reflectivity_right;
twlmOut.mode_confinement_factor = wgSweep.confinement_factor(2);
twlmOut.current_injection_efficiency = twlm.current_injection_efficiency(temperatureInd);

twlm_run_filename=twlmOut.twlm_run_filename;
if(fileexists(twlmOut.twlm_run_filename)) {rm(twlmOut.twlm_run_filename);}
cp(twlmOut.twlm_template_filename,twlmOut.twlm_run_filename);

if(fileexists(twlmOut.twlm_parameter_filename)) {rm(twlmOut.twlm_parameter_filename);}
jsonsave(twlmOut.twlm_parameter_filename,twlmOut);

load(twlmOut.twlm_run_filename);

switchtodesign;

setnamed('TWLM_1','frequency',twlmOut.frequency);
setnamed('TWLM_1','length',twlmOut.length);
setnamed('TWLM_1','active region width',twlmOut.active_region_width);
setnamed('TWLM_1','active region thickness',twlmOut.active_region_thickness);
setnamed('TWLM_1','recombination input parameter','table');
setnamed('TWLM_1','load recombination from file',true);
setnamed('TWLM_1','recombination filename',common.name+'_rcmb.dat');
setnamed('TWLM_1','gain shape','user defined');
setnamed('TWLM_1','load gain from file',true);
setnamed('TWLM_1','gain filename',common.name+'_gain.dat');
setnamed('TWLM_1','diffusion constant',twlmOut.diffusion_constant);
setnamed('TWLM_1','spontaneous emission from gain',false);
setnamed('TWLM_1','load spontaneous emission from file',true);
setnamed('TWLM_1','spontaneous emission filename',common.name+'_spem.dat');
setnamed('TWLM_1','effective index 1',twlmOut.effective_index);
setnamed('TWLM_1','group index 1',twlmOut.group_index);
setnamed('TWLM_1','loss 1',10*log10(exp(1))*twlmOut.loss);
setnamed('TWLM_1','facet reflectivity left 1',twlmOut.facet_reflectivity_left);
setnamed('TWLM_1','facet reflectivity right 1',twlmOut.facet_reflectivity_right);
setnamed('TWLM_1','mode confinement factor 1',twlmOut.mode_confinement_factor);
setnamed('TWLM_1','output carrier density',true);

setnamed('TWLM_1','minimum number of sections',ceil(twlmOut.length/10e-6));

setnamed('TWLM_1','gain fitting number of coefficients',twlm.num_fit_coefficients);
setnamed('TWLM_1','gain fitting tolerance',twlm.fit_tolerance);
setnamed('TWLM_1','gain fitting maximum number of iterations',twlm.num_fit_iterations);
setnamed('TWLM_1','gain fitting rolloff',twlm.fit_rolloff);

setnamed('TWLM_1','spontaneous emission fitting number of coefficients',twlm.num_fit_coefficients);
setnamed('TWLM_1','spontaneous emission fitting tolerance',twlm.fit_tolerance);
setnamed('TWLM_1','spontaneous emission fitting maximum number of iterations',twlm.num_fit_iterations);
setnamed('TWLM_1','spontaneous emission fitting rolloff',twlm.fit_rolloff);

setnamed("TWLM_1","record carrier density profile",true);
setnamed("TWLM_1","record recombination rates profile",true);														 
setnamed('GAIN_1','gain',twlmOut.current_injection_efficiency);

setnamed('DC_1','amplitude',0.2);

setnamed("","simulation input","sample rate");
setnamed("","time window",twlm.time_window);
setnamed("","sample rate",twlm.sample_rate);

sweepName = "sweep";
sweep = getsweep("sweep");
sweepCurrent = getsweep("sweep","current");
myCurrent = struct;
myCurrent.name = sweepCurrent.name;
myCurrent.parameter=sweepCurrent.parameter;
myCurrent.type=sweepCurrent.type;
numcurrents = length(twlm.li_current);
for(i=1:numcurrents){  
    eval("myCurrent.value_"+num2str(i)+"="+"twlm.li_current("+num2str(i)+");");  
}
setsweep(sweepName, "type", "Values");
setsweep(sweepName, "number of points", numcurrents);
setsweep(sweepName, "current", myCurrent);

setnamed("","number of samples",16);
run;
waituntildone;
save(twlmOut.twlm_run_filename);

gain_fitted=getresult("TWLM_1","diagnostic/gain/fitted spectrum");
gain_input=getresult("TWLM_1","diagnostic/gain/measured spectrum");
spem_fitted=getresult("TWLM_1","diagnostic/spontaneous emission/fitted spectrum");
spem_input=getresult("TWLM_1","diagnostic/spontaneous emission/measured spectrum");
savedata(twlmOut.twlm_parameter_filename,gain_fitted,gain_input,spem_fitted,spem_input);

switchtodesign;

setnamed('TWLM_1','load recombination from file',true);
setnamed('TWLM_1','load gain from file',false);
setnamed('TWLM_1','load spontaneous emission from file',false);

setnamed("","time window",twlm.time_window);
save(twlmOut.twlm_run_filename);

runsweep(sweepName);
save(twlmOut.twlm_run_filename);

#Calculate Rstim(I) in units of 1/s/m3:
#Rstim(I<=Ith) = 0,
#Rstim(I>Ith) = c*(J/et-Jth/et).
li = getsweepresult("sweep","power");
current = li.getparameter("current")*twlm.current_injection_efficiency(temperatureInd); #just recombination current
li = li.getattribute("TE power (W)");
current_interp = linspace(current(1),current(end),100);
li_interp = spline(li,current,current_interp);
d2liinterp_di2 = li_interp(3:end)-2*li_interp(2:end-1)+li_interp(1:end-2);
th_index = find(d2liinterp_di2,max(d2liinterp_di2))+1;
Ith = current_interp(th_index);
rst_indices = find(current>=Ith);
sweepfolder = common.name+"_T="+num2str(Temp)+"_"+sweepName;
oldfilename = currentfilename;
cd(sweepfolder);
Rst_avg = zeros(length(rst_indices));
l = getnamed("TWLM_1","length");
t = getnamed("TWLM_1","active region thickness");
w = getnamed("TWLM_1","active region width");
for(i=1;i<=length(rst_indices);i=i+1){
    sweepfilename = sweepName+"_"+num2str(rst_indices(i));
    load(sweepfilename+".icp");
    chden = getresult("TWLM_1","carrier density profile");
    position = chden.getparameter("position");
    chden = chden.getattribute("carrier density (1/m^3)");
    chden = chden(end-end*0.1:end,:); #average over last few time steps
    chden = sum(chden,1)/size(chden,1);
    Rst = getresult("TWLM_1","stimulated emission recombination rate profile");
    Rst = Rst.getattribute("stimulated emission rate (1/s)");
    Rst = Rst(end-end*0.1:end,:);
    Rst = sum(Rst,1)/size(Rst,1);
    Rst_avg(i) = (integrate(Rst*chden,1,position))/l; #averaged over length
}
lincoeff = (Rst_avg(end)-Rst_avg(1))*t*l*w*e/(current(rst_indices(end))-current(rst_indices(1)));
IstimToOutputPower = (Rst_avg(end)-Rst_avg(1))*t*l*w*e/(li(rst_indices(end))-li(rst_indices(1)));

twlmResultOut = struct;
twlmResultOut.Jth = Ith/l/w;
twlmResultOut.lincoeff = lincoeff;
twlmResultOut.IstimOutputPowerRatio = IstimToOutputPower;

load(oldfilename);

if(fileexists(twlmOut.twlm_result_filename)) {rm(twlmOut.twlm_result_filename);}
jsonsave(twlmOut.twlm_result_filename,twlmResultOut);