# DESCRIPTION
# {
#     Strained bandstructure for a single QW as defined in the paper
#     19991, Chuang, Efficient band-structure calculations of strained quantum wells
#
#          1. Band gap from lowest CB to highest heavy hole VB matches the reference for 3 different material compositions (3 different strains).
#          1. Band gap from lowest CB to highest light hole VB matches the reference for 3 different material compositions (3 different strains).
#	3. Band gap for "infinite" QW length matches the reference for 3 different material compositions (3 different strains). 
#             This band gap is directly plotted in the banddiagram dataset of mqwgain solver, no need to actually simulate very long QW.
#
#     NOTES:
#
#       The electric field is zero and hard wall boundary conditions are used. The QW length is 6 nm (as in Fig. 5 of the reference).
# }

clear;
mqw_material_build_functions;

x = [0.37,0.47,0.57];
Lwell = [60]*1e-10;
hhpos = [1,1,2]; #Index of HH1 band (highest heavy hole) as extracted from Fig. 4c (tells which is the lowest band: HH or LH) and Fig. 5
lhpos = [3,2,1]; #same for LH1 band (highest light hole)
c1_hh1 = zeros(length(x),length(Lwell));
c1_lh1 = zeros(length(x),length(Lwell));
bandgapInfLength = zeros(length(x));
for(i=1:length(x)){
    for(j=1:length(Lwell)){
        
        GaInAs = buildTernaryMaterialChuang("Ga(x)In(1-x)As",300,0,x(i));
        GaInAsP = buildQuaternaryMaterialInPMatchedChuang("Ga(x)In(1-x)As(y)P(1-y)",300,0,x(i),0.581); #y=0.581 to give Eg = 0.95 eV as in the reference
        #Custom parameters from the reference
        GaInAs.vb = GaInAsP.vb + 0.64*(GaInAsP.eg-GaInAs.eg);
        GaInAs.lc = (6.0584-0.4051*x(i))*1e-10;
        a_GaAs = -9.77;
        a_InAs = -6;
        a_GaInAs = a_GaAs*x(i)+a_InAs*(1-x(i));
        GaInAs.av = -a_GaInAs/3;
        GaInAs.ac = 2/3*a_GaInAs;
        material = cell(3);
        material{1} = GaInAsP;
        material{2} = GaInAs;
        material{3} = GaInAsP;
        
        #
        # Stack properties
        #
        stackProperties = struct;
        stackProperties.neff = [
        340e12, 3.6;
        370e12, 3.6];
        stackProperties.gamma = 1.317e-2;
        stackProperties.material = material;
        stackProperties.length = [100e-10,Lwell(j),100e-10];
        stackProperties.strain = [0,(GaInAsP.lc-GaInAs.lc)/GaInAsP.lc,0];
        
        #
        # Simulation parameters
        #
        simParameters = struct;
        #simParameters.T = 300;
        simParameters.cden = [3e24];
        #simParameters.V = [0, 0; sum(stackProperties.length), 0];
        simParameters.phfreq = linspace(344.32e12,366.81e12,5);
        InP = buildBinaryMaterialChuang("InP",300,0);
        simParameters.kt = linspace(0,2*pi/InP.lc*0.06,31);
        
        #
        # Simulation config
        #
        bc = struct;
        bc.pmlactive = false;
        bc.hwcutoff = [5e-3,5e-3];
        #bc.pmlcutoff = [1e-2,1e-2];
        #bc.pmllength = [100e-10,100e-10];
        #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;
        
        ##################### Call mqw solver ###########################
        out = mqwgain(stackProperties,simParameters, simcfg);
        
        cb = pinch(out.bandstructure.conduction_band);
        vb = pinch(out.bandstructure.valence_band_lo);
        cb_edge = min(out.banddiagram.Ec);
        vb_edge = max(out.banddiagram.Ev);
        kt = out.bandstructure.kt;
        vb_mass = (h/2/pi)^2/((vb(3,:)-2*vb(2,:)+vb(1,:))*e/(kt(2)-kt(1))^2);
        
        #main results
        c1_hh1(i,j) = cb(1,1) - vb(1,hhpos(i));
        c1_lh1(i,j) = cb(1,1) - vb(1,lhpos(i));
        bandgapInfLength(i) = cb_edge - vb_edge; #this does not depend on well length so it will get overwritten for each length iteration
    }
}