clear;

switchtolayout;
groupscope("::model");
deleteall;

filename = "Piprek2000OQE";

if(!exist('fname')){
    fname = filename;
}
jsonload(fname);
jsonload(fname+"_opt");
jsonload(fname+"_mqw");
jsonload(fname+"_T="+num2str(twlm.temperature)+"_result_twlm");

#In this simulation use single quantum well with total thickness equal to the original multiple quantum wells
common.well_thickness = common.n_wells*common.well_thickness;# [m] the thickness of each of the quantum wells
common.edge1_barrier_thickness = (common.n_wells-1)*common.barrier_thickness+common.edge1_barrier_thickness;# [m] #the thickness of the edge of the top barrier
common.edge2_barrier_thickness = common.edge1_barrier_thickness; # [m] #the thickness of the edge of the bottom barrier
common.n_wells = 1; # the number of quantum wells

#Additional CHARGE simulation properties
charge = struct;
charge.halfwidth = 8e-6;
charge.ridgehalfwidth = 4e-6;
charge.ridgewidthactual = 57e-6;
charge.length = 269e-6;
charge.topcladdingthickness = 0.2e-6;
charge.bottomcladdingthickness = 0.2e-6;
charge.topsclthickness = 0.1e-6;
charge.bottomsclthickness = 0.1e-6;
charge.bottomcladdingthickness = 0.2e-6;
charge.srhtaun = 20e-9;
charge.srhtaup  = 20e-9;
charge.caun = 0;
kb = 1.38064852e-23; # Boltzmann constant J/K
charge.caup = 1.722491012524473e-039*exp(-60e-3/(kb*twlm.temperature/e)); #m^6/s
charge.forwardbiastable = [0,0.25,0.5,0.6,0.7,0.75,0.78,0.82,0.85,0.86,0.87,0.88,0.9,0.92,0.94,0.96,0.98,1.0];


builddevice;


#############################################################################
####################### Run simulation and postprocess ######################
#############################################################################
run("CHARGE");

#All currents below are in positive (x,z) direction except recombination currents which are positive for positive recombination (sink).
#Electron current is the negative of the actual direction of electron motion.

#Total current
anode = getresult("CHARGE","anode");
V = anode.V_anode;
normlength = getnamed("CHARGE","norm length");
surfacefluxA = getresult("CHARGE::injection flux uniform ridge","surfaceflux");
Itotala = -2*anode.I-surfacefluxA.I*(charge.ridgewidthactual-2*charge.ridgehalfwidth)/charge.ridgehalfwidth/uniform_ridge_length_ratio*normlength; #multiplication factor is to scale up to full width considering that current density is mostly uniform along the ridge
cathode = getresult("CHARGE","cathode");
surfacefluxC = getresult("CHARGE::injection flux uniform substrate","surfaceflux");
Itotalc = 2*cathode.I-surfacefluxC.I*(charge.ridgewidthactual-2*charge.ridgehalfwidth)/charge.ridgehalfwidth/uniform_ridge_length_ratio*normlength; #multiplication factor is to scale up to full width considering that current density is mostly uniform along the ridge
Ivertical = -2*anode.In-surfacefluxA.In*(charge.ridgewidthactual-2*charge.ridgehalfwidth)/charge.ridgehalfwidth/uniform_ridge_length_ratio*normlength+
2*cathode.Ip-surfacefluxC.Ip*(charge.ridgewidthactual-2*charge.ridgehalfwidth)/charge.ridgehalfwidth/uniform_ridge_length_ratio*normlength;

#Recombination currents
Rsrh = getdata("CHARGE","recombination","Rsrh");
Rau = getdata("CHARGE","recombination","Rau");
Ropt = getdata("CHARGE","recombination","Ropt");
Rstim = getdata("CHARGE","recombination","Rstim");
vtx = getdata("CHARGE","charge","vertices");
tri = getdata("CHARGE","charge","elements");
numVtx = size(vtx,1);
Rsrh = pinch(Rsrh(:,1,:,1));
Rau = pinch(Rau(:,1,:,1));
Ropt = pinch(Ropt(:,1,:,1));
Rstim = pinch(Rstim(:,1,:,1));
Rsrh_uniform = Rsrh;
Rau_uniform = Rau;
Ropt_uniform = Ropt;
Rstim_uniform = Rstim;
vtx = vtx(:,[1,3]);
#Set recombination rates to zero at substrate vertices which are not under the ridge to avoid double counting the lateral current leakage
#Rsrh(find(vtx(:,1) > charge.ridgehalfwidth),:) = Rau(find(vtx(:,1) > charge.ridgehalfwidth),:) = Ropt(find(vtx(:,1) > charge.ridgehalfwidth),:) = 0;
Rsrh_uniform(find(vtx(:,1) >= uniform_ridge_length_ratio*charge.ridgehalfwidth),:) = Rau_uniform(find(vtx(:,1) >= uniform_ridge_length_ratio*charge.ridgehalfwidth),:) =
Ropt_uniform(find(vtx(:,1) >= uniform_ridge_length_ratio*charge.ridgehalfwidth),:) = Rstim_uniform(find(vtx(:,1) >= uniform_ridge_length_ratio*charge.ridgehalfwidth),:) = 0;

Isrh = zeros(size(Rsrh,2));
Iau = zeros(size(Rau,2));
Iopt = zeros(size(Ropt,2));
Istim = zeros(size(Rsrh,2));
Isrh_uniform = zeros(size(Rsrh,2));
Iau_uniform = zeros(size(Rau,2));
Iopt_uniform = zeros(size(Ropt,2));
Istim_uniform = zeros(size(Rstim,2));
for(i=1; i<=length(V); i=i+1){
    Isrh(i) = quadtri(tri,vtx,Rsrh(:,i))*e*normlength*1e6*2;
    Iau(i) = quadtri(tri,vtx,Rau(:,i))*e*normlength*1e6*2*
    (charge.caun+charge.caup)/(charge.caun+charge.caup+rsp_coeffs(2));
    Iopt(i) = quadtri(tri,vtx,Ropt(:,i))*e*normlength*1e6*2 +
    quadtri(tri,vtx,Rau(:,i))*e*normlength*1e6*2*
    (rsp_coeffs(2))/(charge.caun+charge.caup+rsp_coeffs(2));
    Istim(i) = quadtri(tri,vtx,Rstim(:,i))*e*normlength*1e6*2;
    
    Isrh_uniform(i) = quadtri(tri,vtx,Rsrh_uniform(:,i))*e*normlength*1e6*(charge.ridgewidthactual-2*charge.ridgehalfwidth)/charge.ridgehalfwidth/uniform_ridge_length_ratio;
    Iau_uniform(i) = quadtri(tri,vtx,Rau_uniform(:,i))*e*normlength*1e6*(charge.ridgewidthactual-2*charge.ridgehalfwidth)/charge.ridgehalfwidth/uniform_ridge_length_ratio*
    (charge.caun+charge.caup)/(charge.caun+charge.caup+rsp_coeffs(2));
    Iopt_uniform(i) = quadtri(tri,vtx,Ropt_uniform(:,i))*e*normlength*1e6*(charge.ridgewidthactual-2*charge.ridgehalfwidth)/charge.ridgehalfwidth/uniform_ridge_length_ratio +
    quadtri(tri,vtx,Rau_uniform(:,i))*e*normlength*1e6*(charge.ridgewidthactual-2*charge.ridgehalfwidth)/charge.ridgehalfwidth/uniform_ridge_length_ratio*
    (rsp_coeffs(2))/(charge.caun+charge.caup+rsp_coeffs(2));
    Istim_uniform(i) = quadtri(tri,vtx,Rstim_uniform(:,i))*e*normlength*1e6*(charge.ridgewidthactual-2*charge.ridgehalfwidth)/charge.ridgehalfwidth/uniform_ridge_length_ratio;
}

plot(V,abs(Itotala),abs(Ivertical),abs(Isrh+Isrh_uniform),abs(Iau+Iau_uniform),abs(Iopt+Iopt_uniform),abs(Istim+Istim_uniform),"V anode","current","","log10y,linewidth=2");
legend("Ianode","Ivertical","Isrh","Iau","Iopt","Istim");

plot(abs(Itotala),V,"current [A]","","","linewidth=2,y axis location=left");
holdon;
plot(abs(Itotala),(Istim+Istim_uniform)/twlmResultOut.IstimOutputPowerRatio,"current [A]","","","linewidth=2,y axis location=right");
legend("voltage [V]","power [W]");
holdoff;

dVdI = (V(2:end)-V(1:end-1))/(abs(Itotala(2:end))-abs(Itotala(1:end-1)));
plot(abs(Itotala(2:end)),(abs(Itotala(1:end-1))+abs(Itotala(2:end)))/2*dVdI,"current","Itotal*dV/dItotal","","linewidth=2");
legend("Itotal*dV/dItotal");