### INPUT

filename="Piprek2000OQE";

### END INPUT

if(!exist('fname')){
    fname = filename;
}
jsonload(fname);

mqwOut = struct;
mqwOut.input_filename = fname;
jsonload(fname+"_opt");
mqwOut.opt_output_filename = fname+"_opt";
mqwOut.mqw_result_filename = common.name+"_mqw";

#############################
# layer configuration & strain
#############################
if(mqw.uncoupled_wells){
    n_wells = 1; #solve only for one well since wells are independent and equal in geometry
}
else{
    n_wells = common.n_wells;
}
x_well = mqw.x_well;
y_well = mqw.y_well;
x_barrier = mqw.x_barrier;
y_barrier = mqw.y_barrier;
n_layers = 2*n_wells+ 1;
t_layers = matrix(n_layers,1);
t_layers(1) = mqw.edge1_barrier_thickness;
t_layers(end) = mqw.edge2_barrier_thickness;
si_layers = t_layers; # strain
si_layers(1) = mqw.barrier_strain;
si_layers(end) = mqw.barrier_strain;
qwLength = 0; #qw length for SQW approximation
totalLength = t_layers(1)+t_layers(end); #total length for SQW approximation
for (j = 2:n_layers-1) {
    if (mod(j,2) == 0) {  # even: well
        t_layers(j) = common.well_thickness;
        si_layers(j) = mqw.well_strain;
        qwLength = qwLength + t_layers(j);
    } else { # odd: barrier
        t_layers(j) = common.barrier_thickness;
        si_layers(j) = mqw.barrier_strain;
    }
    totalLength = totalLength + t_layers(j);
}
mqwLength = 0; #total QW length for full MQW structure (different than qwLength variable if mqw.uncoupled_wells = true).
mqwTotalLength = t_layers(1)+t_layers(end); #total length for full MQW structure
mqw_n_layers = 2*common.n_wells+ 1;
for (j = 2:mqw_n_layers-1) {
    if (mod(j,2) == 0) {  # even: well
        mqwTotalLength = mqwTotalLength + common.well_thickness;
        mqwLength = mqwLength + common.well_thickness;
    } else { # odd: barrier
        mqwTotalLength = mqwTotalLength + common.barrier_thickness;
    }
}

lambdamin = mqw.lambdamin;
lambdamax = mqw.lambdamax;
######################
# Stack properties, except material properties that depend on temperature and charge density
######################
stackProperties = struct;
stackProperties.neff = [
c/lambdamax, real(wgSweep.neff(2));
c/lambdamin, real(wgSweep.neff(2))];
stackProperties.gamma = 2*h/2/pi/mqw.intraband_relaxation_tau/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];
if(common.frequpsampletype < 2){
    simParameters.phfreq = flip(c/linspace(mqw.lambdamin,mqw.lambdamax,mqw.numfreq),1);
}
else{
    simParameters.phfreq = chpts(c/mqw.lambdamax,c/mqw.lambdamin,mqw.numfreq,2);
}

#######################
# Simulation config
#######################
bc = struct;
bc.pmlactive = mqw.pmlactive;
#bc.hwcutoff = [5.e-4,5e-4];
bc.pmlcutoff = mqw.pmlcutoff;
#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 = mqw.dz;
#simcfg.numeigenvalues = 30;
simcfg.materialdb = struct;

#######################################
#Sweep parameters and loop over them 
#######################################
mat_db = struct; #use default database when empty struct
temperatureVec = mqw.temperature;
cdenVec = mqw.charge_density; #m^{-3};
mqwOut.emission = struct;
mqwOut.emission.frequency=simParameters.phfreq;
mqwOut.emission.te = mqwOut.emission.tm=struct;
mqwOut.emission.te.stim=mqwOut.emission.tm.stim=
mqwOut.emission.te.spont=mqwOut.emission.tm.spont=
zeros(length(cdenVec),length(temperatureVec),length(mqwOut.emission.frequency));
mqwOut.Egk0 = zeros(length(cdenVec),length(temperatureVec));
mqwOut.strainedBandgaps = zeros(length(cdenVec),length(temperatureVec),n_layers);
sqw_to_mqw_conversion = 1;
if(mqw.uncoupled_wells){
    sqw_to_mqw_conversion = totalLength/mqwTotalLength*common.n_wells; #transform SQW average density to MQW average density and same for gain
}
#Convert results to represent the physical footprint in case the buffer layers (first and last barrier) are set to be thinner in MQW solver than in the physical footprint for efficiency reasons
#(making them thinner is OK as long as the wavefunction decays sufficiently at the edges of the active region so that the boundary conditions can work properly).
lengthConversion = mqwTotalLength/(mqwTotalLength+common.edge1_barrier_thickness+common.edge2_barrier_thickness-mqw.edge1_barrier_thickness-mqw.edge2_barrier_thickness);
for(temperatureInd = 1:length(temperatureVec)){
    for(cdenInd = 1:length(cdenVec)){
        T_sim = temperatureVec(temperatureInd);
        cden_sim = cdenVec(cdenInd)/sqw_to_mqw_conversion/lengthConversion; #cdenVec is for MQW structure, cden_sim may be in SQW approximation and with thinner buffer layers
        
        ###############################
        # Temperature and charge dependent material properties
        ###############################
        layer_material = cell(n_layers);
        #Calculate these at 300K and apply bandgap temperature gradient from the paper.
        InGaAsP_well = buildmqwmaterial(mat_db, 300, "InGaAsP",1-x_well,y_well,"Gamma");       
        InGaAsP_barrier = buildmqwmaterial(mat_db, 300, "InGaAsP",1-x_barrier,y_barrier,"Gamma");
        #Override bandgap with custom formula
        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;

        #Fitting (if nonzero)
        InGaAsP_well.eg = InGaAsP_well.eg + mqw.well_Eg_shift;
        InGaAsP_barrier.eg = InGaAsP_barrier.eg + mqw.barrier_Eg_shift;
        
        #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 paper ratio of dEc/dEv = 0.4/0.6 for band offsets
        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,81);        
        
        #################################################################
        ##################### Call mqw solver ###########################
        #################################################################
        out = mqwgain(stackProperties,simParameters,simcfg);
                
        ####################################################################################
        ##################### Calculate and print true band gap from bandstructure ###################
        ####################################################################################
        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).";
        
        #export results
        mqwOut.Egk0(cdenInd,temperatureInd) = Egk0;
        mqwOut.emission.te.stim(cdenInd,temperatureInd,:) = out.emission.stimulated_TE*sqw_to_mqw_conversion*lengthConversion;
        mqwOut.emission.te.spont(cdenInd,temperatureInd,:) = out.emission.spontaneous_TE*sqw_to_mqw_conversion*lengthConversion;
        mqwOut.emission.tm.stim(cdenInd,temperatureInd,:) = out.emission.stimulated_TM*sqw_to_mqw_conversion*lengthConversion;
        mqwOut.emission.tm.spont(cdenInd,temperatureInd,:) = out.emission.spontaneous_TM*sqw_to_mqw_conversion*lengthConversion;
        
        ec = out.banddiagram.Ec;
        ev = out.banddiagram.Ev;
        position = out.banddiagram.z;
        layer_distance = 0;
        mqwOut.strainedBandgaps(cdenInd,temperatureInd,1) = ec(1)-ev(1);
        for(i=2; i<=n_layers; i=i+1){
            layer_distance = layer_distance + t_layers(i-1);
            for(j=1; position(j)<layer_distance; j=j+1){
            }
            mqwOut.strainedBandgaps(cdenInd,temperatureInd,i) = ec(j)-ev(j);
        }
    }
}
holdoff;

mqwOut.cdenVec = cdenVec;
mqwOut.temperatureVec=temperatureVec;
mqwOut.qwLength = mqwLength;
mqwOut.totalLength = mqwTotalLength/lengthConversion; #to get physical MQW length in case buffer layers (first and last barrier) were made thinner for simulation speed up
mqwOut.y_well = y_well;
jsonsave(mqwOut.mqw_result_filename,mqwOut);