# DESCRIPTION
# {
#     TE mode maximum gain vs temperature for a single QW as defined in the paper
#     1988, Ahn, Valence band mixing effects on the gain and the refractive index change of quantum well lasers
#     
#     NOTES:
#
#       The electric field is zero and pml boundary conditions are used. The assumption is that the band gap does not change with temperature
#       (temperature dependence not clearly stated in the reference).
#
# }

clear;

mqw_material_build_functions;

GaAs = buildBinaryMaterialChuang("GaAs",300,0);
AlGaAs = buildTernaryMaterialChuang("Al(x)Ga(1-x)As",300,0,0.25);
GaAs.me = 0.071; #This gives a better match of the temperature gradient of gain than 0.067 from Chuang (in some older papers of Chuang, like 1991, paper this value was used actually).
AlGaAs.me = GaAs.me;
AlGaAs.gamma1 = GaAs.gamma1;
AlGaAs.gamma2 = GaAs.gamma2;
AlGaAs.gamma3 = GaAs.gamma3;
GaAs.vb = AlGaAs.vb + (AlGaAs.eg-GaAs.eg)*0.43; #dEc=0.57/dEg, equivalently dEv=0.43/dEg (as in the reference)

#Cell for material properties
material = cell(3);
material{1} = AlGaAs;
material{2} = GaAs;
material{3} = AlGaAs;

############### Input arguments for mqwgain command #####################
stackProperties = struct;
stackProperties.neff = [
    340e12, 3.6;
    370e12, 3.6];
stackProperties.gamma = 1.317e-2;
stackProperties.material = material;
qwLength = 100e-10;
barrierLength = 100e-10;
stackProperties.length = [barrierLength,qwLength,barrierLength];
#stackProperties.strain = [0];

lengthConversion = sum(stackProperties.length)/qwLength;

simParameters = struct;
simParameters.cden = [3e24/lengthConversion]; #length over which density is averaged is total length (barriers+QW), while in the reference it is just QW length
#simParameters.V = [0, 0; sum(stackProperties.length), 0];
simParameters.phfreq = linspace(344.32e12,366.81e12,51);
simParameters.kt = linspace(0,2*pi/GaAs.lc*0.12,71);

bc = struct;
bc.pmlactive = true;
#bc.hwcutoff = [5.e-4,5e-4];
bc.pmlcutoff = [1e-3,1e-11];
#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];

simConfig = struct;
simConfig.bcs = bc;
simConfig.materialdb = struct;
#simConfig.dz = 1e-10;
#simConfig.numeigenvalues = 30;

##################### Call mqw solver ###########################
Tvec = [100,200,300,400];
maxgainVsT = zeros(length(Tvec));
for(i=1:length(Tvec)){
    simParameters.T = Tvec(i);
    out = struct;
    out = mqwgain(stackProperties,simParameters,simConfig);
    
    #Convert lengths to compare with reference (we consider total length, not just QW length)
    maxgainVsT(i) = max(out.emission.stimulated_TE*lengthConversion);
}

#Plot main result
plot(Tvec,maxgainVsT);