include "../fsem_main_Include.edp";
fsemVERBOSITY = 10;
fsemModUseCpp = true;
verbosity = 0;
UsrNewtonNbrMaxIterCarr = 20;


bool wf = true; // Wait Flag for plots


/****************************************************************
 * Define the mesh. The mesh must be adimentional (I will try to 
 * make that automatice in the future)
 ****************************************************************/
include "meshes/meshunif2d.edp";
include "meshes/meshpn2d.edp";

real Thick = 100e-6 / GloScaX;    // wafer thickness [m]
real Width = 100e-6 / GloScaX;    // wafer width (radius) [m]
real Radius = 1.5e-6 / GloScaX;   // laser radius [m]
real AbsLength = 20e-6 / GloScaX; // Laser absoprtion length [m]
real Xj = 50e-9 / GloScaX;        // junction depth (m) 
real L0 = 50e-6 / GloScaX;        // Lateral doping width

//meshunif2d(GloMesh, Width, Thick, Radius, AbsLength);
meshpn2d(GloMesh, Width, Thick, Xj, L0, Radius, AbsLength);  // define new mesh with Ultra-Shallow Layer (USL)

// Get region number in order to define the P+N junction later on.
int USL1 = GloMesh(1e-9/GloScaX, -1e-9/GloScaX).region;            // Ultra Shallow Layer (USL)
int USL2 = GloMesh(L0-1e-9/GloScaX, 0.0).region;
int USL3 = GloMesh(L0+1e-9/GloScaX, 0.0).region;
int SUB1 = GloMesh(0.0, -Thick).region;
int SUB2 = GloMesh(0.0, -Xj-1e-9/GloScaX).region;

plot(GloMesh, wait=wf);


/*****************************
 * Doping parameters
 *****************************/
real dopingSUB =1e15 * ConvC2M^(-3); // m^-3;
real dopingUSL =1e19 * ConvC2M^(-3); // m^-3;


/***************************************
 * Solve for each doping concentration
 ****************************************/

ParUsrPDOP = 0;
ParUsrNDOP = dopingSUB;
GloVarSemiMOBELE = 500 * ConvC2M^(2) ;
GloVarSemiMOBHOL = 500 * ConvC2M^(2);

/****************
 * Set laser(s)
 ****************/
 
ParUsrLaserN = 1;                          /* number of lasers */
ParUsrLaserE[0] = 1.24 / 0.670;            /* energy [eV] */
ParUsrLaserI[0] = 1.59e5 * ConvC2M^(-2);   /* Irradiance [W m-2] */
ParUsrLaserRefl[0] =  0.303;               /* Reflectance [adim] */
ParUsrLaserR[0] = 1 * ConvU2M;             /* Radius [m] */
  
/*******************
 * Set model flags
 *******************/

ModUsrOptCarrGen=1;   // Switch on optical carrier generation
ModUsrBTBIP=1;        // Band-to-band absorption model (Smith)
ModUsrBGNSchenk=0;
ModUsrBGNSlotboom=1;
ModUsrRecAuger=1;   
ModUsrRecSRH=1;
ModUsrRecSRHConcDepLifeTime=1;
//ModUsrMob = ModCSTMobCst; 
ModUsrMob = ModCSTMobPHUNI;

/****************************
 * Set material parameters
 ***************************/
ParMatxyzMobn = 500  * ConvC2M^(2) ;
ParMatxyzMobp = 500  * ConvC2M^(2) ;

/***********************
 * Set initial value
 ***********************/

GloVarSemiPOT =0.0;
GloVarSemiHOL=1.448141431e+10 * 1.448141431e+10   * ConvC2M^(-6) / ParUsrNDOP;
GloVarSemiELE=ParUsrNDOP;
GloVarTempT = 300;

/***************************
 * Set solver info
 *******************/

UsrNewtonEpsIncrPot = 1e-10;
UsrGummelEpsPot = 1e-6;
UsrGummelEpsCarr = 1e-6;



/**************************************************************
 **************************************************************
 **************************************************************
 * Solve the drift-diffusion equations in slotboom varaibles
 *
 *   UNIFORM DOPING (substrate)
 *
 **************************************************************/

GloVhP1 AdaptFunc;
UsrGummelEpsPot = 1e-5;
UsrGummelEpsCarr = 1e-5;



SolveDDSG(true);  //true because Ntype --> Hole equation should be solved before


  // DO NOT include potential for refinement in the case of a uniform doping (too noisy) 
  plot(GloVarSemiHOL , wait=wf, value=true); 
  AdaptFunc = 0; GloMeshGetAdaptFunc(AdaptFunc[], GloVarSemiCARRGEN[0][], GloVarSemiDDSGRhoP[], GloVarSemiDDSGRhoN[]);
  GloMesh = adaptmesh(GloMesh, AdaptFunc, err=0.001, nbsmooth=5, nbvx=100000, iso=false);

ModUsrEleNonlinear = true;  // next resolution with Newton solver for RhoN ...
ModUsrHolNonlinear = true;  // ... and RhoP
UsrGummelEpsPot = 1e-7;
UsrGummelEpsCarr = 1e-7;

SolveDDSG(true);

  plot(GloVarSemiPOT , wait=wf, value=true);
  plot(GloVarSemiHOL , wait=wf, value=true);
  plot(GloVarSemiELE ,  wait=wf, value=true);


/**************************************************************
 **************************************************************
 **************************************************************
 * Solve the drift-diffusion equations in slotboom variables
 *
 *   P+/N junction
 *
 **************************************************************/


// Define P+N doping 
ParUsrPDOP = dopingUSL * (region==USL1 | region==USL2 | region==USL3 ) + 0.0 * (region==SUB1 | region==SUB2);    
ParUsrNDOP = dopingSUB * (region==SUB1 | region==SUB2) + 0.0 * (region==USL1 | region==USL2 | region==USL3);
GloVarSemiDOPNET = ParUsrNDOP - ParUsrPDOP;


for(int i=0; i<3; i++) {
  
  AdaptFunc = 0; GloMeshGetAdaptFunc(AdaptFunc[], GloVarSemiDOPNET[], GloVarSemiCARRGEN[0][], GloVarSemiDDSGRhoP[], GloVarSemiDDSGRhoN[]);
  GloMesh = adaptmesh(GloMesh, AdaptFunc, err=0.006, nbsmooth=5, nbvx=100000, iso=false);
  
  ParUsrPDOP = dopingUSL * (region==USL1 | region==USL2 | region==USL3 ) + 0.0 * (region==SUB1 | region==SUB2);    
  ParUsrNDOP = dopingSUB * (region==SUB1 | region==SUB2) + 0.0 * (region==USL1 | region==USL2 | region==USL3);
  GloVarSemiDOPNET = ParUsrNDOP - ParUsrPDOP;
  GloVarSemiCARRGEN[0] =  GloVarSemiCARRGEN[0];
  GloVarSemiDDSGRhoP = GloVarSemiDDSGRhoP;
  GloVarSemiDDSGRhoN = GloVarSemiDDSGRhoN;
  plot(GloVarSemiDOPNET, wait=wf, value=true);

 }

//ModUsrEleNonlinear = true;
//ModUsrHolNonlinear = true;
//UsrGummelEpsPot = 1e-7;
//UsrGummelEpsCarr = 1e-7;


SolveDDSG(true);

plot(GloVarSemiPOT , wait=wf, value=true);
plot(GloVarSemiHOL , wait=wf, value=true);
plot(GloVarSemiELE ,  wait=wf, value=true);

// Resolution with Schenk's BGN model

ModUsrBGNSlotboom=0;
ModUsrBGNSchenk=1;
SolveDDSG(true);


plot(GloVarSemiPOT , wait=true, value=true);
plot(GloVarSemiHOL , wait=true, value=true);
plot(GloVarSemiELE ,  wait=true, value=true);

                                                                                                                                       
AdaptFunc = 0; GloMeshGetAdaptFunc(AdaptFunc[], GloVarSemiDOPNET[], GloVarSemiCARRGEN[0][], GloVarSemiDDSGRhoP[], GloVarSemiDDSGRhoN[], GloVarSemiPOT[]);
GloMesh = adaptmesh(GloMesh, AdaptFunc, err=0.006, nbsmooth=5, nbvx=100000, iso=false);

SolveDDSG(true);
                                                                                                                                                        
 
GloSaveGmshPosSPARSEDST(GloMesh, GloVarSemiPOT, 1, "GloVarSemiPOT", "GloVarSemiPOT.pos", 1e6, GloScaX, GloScaY);
GloSaveGmshPosSPARSEDST(GloMesh, GloVarSemiHOL, ConvM2C^(-3), "GloVarSemiHOL", "GloVarSemiHOL.pos", 1e6, GloScaX, GloScaY);
GloSaveGmshPosSPARSEDST(GloMesh, GloVarSemiELE, ConvM2C^(-3), "GloVarSemiELE", "GloVarSemiELE.pos", 1e6, GloScaX, GloScaY);
                                                                                                                                        


/*********************************
 * Test (for debugging purpose)
 *********************************/

x=0;
y=0;
real GloVarSemiPOT00 = GloVarSemiPOT;
real GloVarSemiPOT00saved =  -0.03601889679;
real epspot = abs(GloVarSemiPOT00saved) * 1e-3;
real GloVarSemiELE00 = GloVarSemiELE;
real GloVarSemiELE00saved =   8.882481043e+23 ;
real epsele = abs(GloVarSemiELE00saved) * 1e-3;
real GloVarSemiHOL00 = GloVarSemiHOL;
 real GloVarSemiHOL00saved = 8.872483954e+23;
real epshol = abs(GloVarSemiHOL00saved) * 1e-3;

cout << "GloVarSemiPOT00saved: " << GloVarSemiPOT00saved << endl;
cout << "GloVarSemiPOT00: " << GloVarSemiPOT00 << endl;

cout << "DEBUG: min, max: " << GloVarSemiPOT[].min << ", " << GloVarSemiPOT[].max << endl;
//notice the perfectly flat potential due to equal values of mobp and mobn. if taken different, non zero potential gradient because of hole-electron dissymetry. 

cout << "GloVarSemiELE00saved: " << GloVarSemiELE00saved << endl;
cout << "GloVarSemiELE00: " << GloVarSemiELE00 << endl;
cout << "DEBUG: min, max: " << GloVarSemiELE[].min << ", " << GloVarSemiELE[].max << endl;

cout << "GloVarSemiHOL00saved: " << GloVarSemiHOL00saved << endl;
cout << "GloVarSemiHOL00: " << GloVarSemiHOL00 << endl;
cout << "DEBUG: min, max: " << GloVarSemiHOL[].min << ", " << GloVarSemiHOL[].max << endl;

if(GloVarSemiPOT00 > GloVarSemiPOT00saved-epspot && GloVarSemiPOT00 < GloVarSemiPOT00saved+epspot
&& GloVarSemiELE00 > GloVarSemiELE00saved-epsele && GloVarSemiELE00 < GloVarSemiELE00saved+epsele
&& GloVarSemiHOL00 > GloVarSemiHOL00saved-epshol && GloVarSemiHOL00 < GloVarSemiHOL00saved+epshol


) {
  cout << endl << "*** 007b: Test successfully passed ! ***" << endl;
} 
else {
  cout << endl << "*** ERROR: expected value was not obtained "  << endl;
  exit(0);
}