# DESCRIPTION:
# Test of mqwindex 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 and pml boundary conditions are used.

clear;

#use custom material properties from the reference instead of our database to make comparison fair
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;
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 mqwindex 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 MATLAB and in the reference it is just QW length
simParameters.V = [
    0,0;
    300e-10,150e-3];
simParameters.phfreq = linspace(344.32e12,392.68e+12,51);
simParameters.kt = linspace(0,2*pi/GaAs.lc*0.3,201);

bc = struct;
bc.pmlactive = true;
#bc.hwcutoff = [1e-2,4e-3];
bc.pmlcutoff = [5e-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 = mqwindex(stackProperties,simParameters,simConfig);

#Extract real part of refractive index into a dataset
index =  matrixdataset("index");
index.addparameter("energy", out.index.energy);
index.addparameter("ndensity", out.index.ndensity);
index.addattribute("index_TE", real(out.index.index_TE)*lengthConversion/stackProperties.neff(1,2));
index.addattribute("index_TM", real(out.index.index_TM)*lengthConversion/stackProperties.neff(1,2));

#Plot real part of refractive index
for(i=1; i<=length(simParameters.cden); i=i+1){
    plot(index.energy,index.index_TE(:,i),"energy [eV]","relative refractive index change [%]","","linewidth=2");
    holdon;
    plot(index.energy,index.index_TM(:,i),"energy [eV]","relative refractive index change [%]","","linewidth=2");
    holdon;
    legend("simulation TE","simulation TM");
    setplot("title","density = " + num2str(simParameters.cden(i)) + " cm^{-3}");
    holdoff;
}

#Extract imaginary part of refractive index into a dataset
imagindex =  matrixdataset("index");
imagindex.addparameter("energy", out.index.energy);
imagindex.addparameter("ndensity", out.index.ndensity);
imagindex.addattribute("index_TE", imag(out.index.index_TE)*lengthConversion/stackProperties.neff(1,2));
imagindex.addattribute("index_TM", imag(out.index.index_TM)*lengthConversion/stackProperties.neff(1,2));
for(i=1; i<=length(simParameters.cden); i=i+1){
    plot(index.energy,imagindex.index_TE(:,i),"energy [eV]","imaginary part of index","","linewidth=2");
    holdon;
    plot(index.energy,imagindex.index_TM(:,i),"energy [eV]","imaginary part of index","","linewidth=2");
    holdon;
    legend("simulation TE","simulation TM");
    setplot("title","density = " + num2str(simParameters.cden(i)) + " cm^{-3}");
    holdoff;
}

#Plot wavefunctions on top of banddiagram.
plot_wavefunction_banddiagram;
plot_wf_on_banddiagram(out,3,4);