clear;

#############################
# layer configuration & strain
#############################
n_wells = 6;
x_well = 0.24;
y_well = 0.79;
x_barrier = 0.29;
y_barrier = 0.55;
n_layers = 2*n_wells + 1;
t_layers = matrix(n_layers,1);
t_layers(1) = 17e-9;
t_layers(end) = 17e-9;
si_layers = t_layers; # strain
si_layers(1) = 0.0004;
si_layers(end) = 0.0004;
qwLength = 0;
totalLength = t_layers(1)+t_layers(end);
for (j = 2:n_layers-1) {
    if (mod(j,2) == 0) {  # even: well
        t_layers(j) = 6.4e-9;
        si_layers(j) = -0.01;
        qwLength = qwLength + t_layers(j);
    } else { # odd: barrier
        t_layers(j) = 5.5e-9;
        si_layers(j) = 0.0004;
    }
    totalLength = totalLength + t_layers(j);
}

lambdamin = 1.45e-6;
lambdamax = 1.62e-6;
######################
# Stack properties, except material properties that depend on temperature and charge density
######################
stackProperties = struct;
stackProperties.neff = [
c/lambdamax, 3.1;
c/lambdamin, 3.1];
stackProperties.gamma = 2*h/2/pi/1e-13/e; #intraband relaxation tau = 0.1 ps
stackProperties.length = t_layers;
stackProperties.strain = si_layers;

#########################
# Simulation parameters, except temperature and charge density and material dependent params (like kt)
#########################
simParameters = struct;
simParameters.V = [0, 0; sum(stackProperties.length), 0];
simParameters.phfreq = flip(c/linspace(lambdamin,lambdamax,51),1);

#######################
# Simulation config
#######################
bc = struct;
bc.pmlactive = true;
#bc.hwcutoff = [5.e-4,5e-4];
bc.pmlcutoff = [1e-5,1e-3];
#bc.pmllength = [10e-9,10e-9];
#bc.pmlcoefficient = [0.5+1i*0.5,0.5+1i*0.5,-1+1i*1.4,-1+1i*1.4];
simcfg = struct;
simcfg.bcs = bc;
simcfg.dz = 1e-10;
#simcfg.numeigenvalues = 30;
simcfg.materialdb = struct;

#######################################
#Sweep parameters and loop over them
#######################################
default_mat_db = struct;
temperatureVec = [293,333,373];
cdenVec = 0.8e24; #m^{-3}; this is average density over the entire MQW thickness (density in quantum wells will be larger than this)
emission = cell(length(temperatureVec)*length(cdenVec));
for(temperatureInd = 1:length(temperatureVec)){
    for(cdenInd = 1:length(cdenVec)){
        T_sim = temperatureVec(temperatureInd);
        cden_sim = cdenVec(cdenInd);
        
        ###############################
        # Temperature and charge dependent material properties
        ###############################
        layer_material = cell(n_layers);
        #For default database calculate these at 300K and apply bandgap temperature gradient from Piprek 2000.
        InGaAsP_well = buildmqwmaterial(default_mat_db, 300, "InGaAsP",1-x_well,y_well,"Gamma");
        InGaAsP_barrier = buildmqwmaterial(default_mat_db, 300, "InGaAsP",1-x_barrier,y_barrier,"Gamma");
        #Override bandgap with custom formula from Chuang, Physics of Optoelectronic Devices
        InGaAsP_well.eg = 1.35+0.668*x_well-1.068*y_well+0.758*x_well^2+0.078*y_well^2-0.069*x_well*y_well-0.322*x_well^2*y_well+0.03*x_well*y_well^2;
        InGaAsP_barrier.eg = 1.35+0.668*x_barrier-1.068*y_barrier+0.758*x_barrier^2+0.078*y_barrier^2-0.069*x_barrier*y_barrier-0.322*x_barrier^2*y_barrier+0.03*x_barrier*y_barrier^2;
        #Include bandgap temperature gradient
        InGaAsP_well.eg = InGaAsP_well.eg - 0.28e-3*(T_sim-300);
        InGaAsP_barrier.eg = InGaAsP_barrier.eg - 0.28e-3*(T_sim-300);
        #Override deformation potentials from Chuang 1991
        a_GaAs = -9.77;
        a_InAs = -6;
        a_GaP = -8.83;
        a_InP = -6.31;
        a_GaInAsP_well = a_GaAs*x_well*y_well + a_InP*(1-x_well)*(1-y_well) + a_InAs*(1-x_well)*y_well + a_GaP*x_well*(1-y_well);
        a_GaInAsP_barrier = a_GaAs*x_barrier*y_barrier + a_InP*(1-x_barrier)*(1-y_barrier) + a_InAs*(1-x_barrier)*y_barrier + a_GaP*x_barrier*(1-y_barrier);
        InGaAsP_well.av = -a_GaInAsP_well/3;
        InGaAsP_well.ac = a_GaInAsP_well/3;
        InGaAsP_barrier.av = -a_GaInAsP_barrier*2/3;
        InGaAsP_barrier.ac = a_GaInAsP_barrier*2/3;
        
        #Band gap shrinkage due to carrier-carrier interaction (consider all density is inside QW and none in the barriers).
        dEg = -1e-10*(cden_sim*totalLength/qwLength)^(1/3); #QW density is larger than average density
        InGaAsP_well.eg = InGaAsP_well.eg + dEg;
        
        #Cell for material properties
        for (j = 1:length(layer_material)) {
            layer_material{j} = struct;
            if (mod(j,2) == 0) {  # even: well
                layer_material{j} = InGaAsP_well;
            } else { # odd: barrier
                layer_material{j} = InGaAsP_barrier;
            }
        }
        ?"          well        barrier";
        ?"Band gap  "+num2str(InGaAsP_well.eg)+", "+num2str(InGaAsP_barrier.eg);
        ?"me        "+num2str(InGaAsP_well.me)+", "+num2str(InGaAsP_barrier.me);
        
        stackProperties.material = layer_material;
        # Override band offsets using the ratio of dEc/dEv = 0.4/0.6 for band offsets from Piprek 2000
        dEg = InGaAsP_barrier.eg - InGaAsP_well.eg;
        dEv = 0.6*dEg;
        vb_off = matrix(n_layers,1);
        vb_off(1:2:end) = -dEv;
        stackProperties.vb = struct;
        stackProperties.vb.method = "direct";
        stackProperties.vb.offsets = vb_off;
        
        #########################
        # Temperature and charge dependent simulation parameters
        #########################
        simParameters.T = T_sim;
        simParameters.cden = cden_sim;
        minLCStrained = stackProperties.material{1}.lc*(1.+si_layers(1));
        for(i=1:n_layers){
            minLCStrained = min([minLCStrained,stackProperties.material{i}.lc*(1.+si_layers(i))]);
        }
        simParameters.kt = linspace(0,2*pi/minLCStrained*0.2,75);
        
        #################################################################
        ##################### Call mqw solver ###########################
        #################################################################
        out = mqwgain(stackProperties,simParameters,simcfg);
        
        ####################################################################################
        ##################### Find true band gap from band structure (includes quantum confinement effect)
        ####################################################################################
        Ek = out.bandstructure;
        nCBsubbands = size(Ek.conduction_band,2);
        nVBUsubbands = size(Ek.valence_band_up,2);
        nVBLsubbands = size(Ek.valence_band_lo,2);
        minEc = min(Ek.conduction_band(:,1));
        maxEv = max(Ek.valence_band_up(:,1));
        for (j = 1:nCBsubbands) {
            if (finite(Ek.conduction_band(1,j))) {
                minEc = min([minEc; max(Ek.conduction_band(:,j))]);
            }
        }
        for (j = 1:nVBUsubbands) {
            if (finite(Ek.valence_band_up(1,j))) {
                maxEv = max([maxEv; max(Ek.valence_band_up(:,j))]);
            }
        }
        for (j = 1:nVBLsubbands) {
            if (finite(Ek.valence_band_lo(1,j))) {
                maxEv = max([maxEv; max(Ek.valence_band_lo(:,j))]);
            }
        }
        if (0) { holdoff; }
        Egk0 = minEc - maxEv;
        Egk0_lambda = h*c/Egk0/e;
        ?"Bandgap from E(k=0): "+num2str(Egk0)+" eV ("+num2str(Egk0_lambda*1e6)+" um).";
        
        emission{cdenInd + (temperatureInd-1)*length(cdenVec)} = out.emission;
    }
}

#Plot simulated gain against reference.
piprek2000_fig7;
plot(piprek2000_fig7_20C(:,1),piprek2000_fig7_20C(:,2),"lambda (um)","Gain (1/cm)","","plot type=point,marker style=x");
holdon;
plot(piprek2000_fig7_60C(:,1),piprek2000_fig7_60C(:,2),"","","","plot type=point,marker style=x");
plot(piprek2000_fig7_100C(:,1),piprek2000_fig7_100C(:,2),"","","","plot type=point,marker style=x");
for(temperatureInd = 1:length(temperatureVec)){
    for(cdenInd = 1:length(cdenVec)){
        plot(emission{cdenInd + (temperatureInd-1)*length(cdenVec)}.wavelength*1e6, 
            emission{cdenInd + (temperatureInd-1)*length(cdenVec)}.stimulated_TE/100*totalLength/qwLength,"","","","linewidth=2"); # to 1/cm
    }
}
legend("ref,T=293K","ref,T=333K","ref,T=373K","sim,T=293K","sim,T=333K","sim,T=373K");
holdoff;