# DESCRIPTION
# {
#     Gain 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 50 kV/cm applied from 0-300A (total length). Hard wall boundary conditions are used.
#
# }

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.063; #0.067 in Chuang's textbook, but it was not specified in the reference and this gives a better match for conduction bands.
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)
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.T = 300;
simParameters.cden = [1e24,2e24,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;
    300e-10,150e-3];
simParameters.phfreq = linspace(344.32e12,366.81e12,51);
simParameters.kt = linspace(0,2*pi/GaAs.lc*0.09,51); #convergence with size of Brillouin zone is at around 12%, but it is not clear how much was usd in the reference. 9% gives a better match.

bc = struct;
bc.pmlactive = false;
bc.hwcutoff = [1e-2,4e-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];

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

##################### Call mqw solver ###########################
out = struct;
out = mqwgain(stackProperties,simParameters,simConfig);

#Convert lengths (we consider total length, not just QW length)
emission = matrixdataset("emission");
emission.addparameter("wavelength",out.emission.wavelength,"energy",out.emission.energy,"frequency",out.emission.frequency);
emission.addparameter("ndensity",out.emission.ndensity);
emission.addattribute("stimulated_TE",out.emission.stimulated_TE*lengthConversion);
emission.addattribute("stimulated_TM",out.emission.stimulated_TM*lengthConversion);
emission.addattribute("spontaneous_TE",out.emission.spontaneous_TE*lengthConversion);
emission.addattribute("spontaneous_TM",out.emission.spontaneous_TM*lengthConversion);

visualize(emission);