%% Solve the scaled stationary bipolar DD equation system using Newton's method, for a PIN diode geometry 
%% 
%%
%% [n, p, V, Fn, Fp, Jn, Jp, it, res] = secs1d_dd_newton (x, D, Vin, nin, 
%%                                                        pin, l2, er, un, 
%%                                                        up, theta, tn, tp, 
%%                                                        Cn, Cp, toll, maxit)
%%
%%     input: 
%%       x                spatial grid
%%       D                doping profile
%%       Vin              initial guess for electrostatic potential
%%       nin              initial guess for electron concentration
%%       pin              initial guess for hole concentration
%%       l2               scaled Debye length squared
%%       er               relative electric permittivity
%%       un               electron mobility model coefficients
%%       up               electron mobility model coefficients
%%       theta            intrinsic carrier density
%%       tn, tp, Cn, Cp   generation recombination model parameters
%%       toll             tolerance for Gummel iterarion convergence test
%%       maxit            maximum number of Gummel iterarions
%%
%%     output: 
%%       n     electron concentration
%%       p     hole concentration
%%       V     electrostatic potential
%%       Fn    electron Fermi potential
%%       Fp    hole Fermi potential
%%       Jn    electron current density
%%       Jp    hole current density
%%       it    number of Gummel iterations performed
%%       res   total potential increment at each step



 % physical constants and parameters
 pkg load secs1d;
 secs1d_physical_constants;
 secs1d_silicon_material_properties;
 
 % geometry
 L  = 180*10^-6;           % [m] 
 x_pi_jun = 60*10^-6;      %   coord of p-i junction (m)
 x_in_jun = 120*10^-6;     %   coord of i-n junction (m) 
 
 Nelements = 10000;
 x         = linspace (0, L, Nelements+1)';
 sinodes   = [1:length(x)];
 
 % dielectric constant (silicon)
 er = esir * ones (Nelements, 1);
 
 % doping profile [m^{-3}]

 Na = 1e16 * (x <= x_pi_jun);
 Nd = 1e16 * (x > x_in_jun);
 No = 1e10 * (x > x_pi_jun) .* (x <= x_in_jun) ;
 
 % avoid zero doping
 D  = Nd - Na +No; 
  
 % initial guess for n, p, V, phin, phip
 %V_p = -1;
 V_p =  0;  %changed V_p from -1 to 0, to start from initial zero voltage. correct?
 V_n =  0;
 
 Fp = V_p * (x <= x_pi_jun);
 Fn = Fp;
 
 p = abs (D) / 2 .* (1 + sqrt (1 + 4 * (ni./abs(D)) .^2)) .* (x <= x_pi_jun) + ...
     ni .* (x > x_pi_jun) .* (x <= x_in_jun) + ...
     ni^2 ./ (abs (D) / 2 .* (1 + sqrt (1 + 4 * (ni ./ abs (D)) .^2))) .* (x > x_in_jun);
 
 n = abs (D) / 2 .* (1 + sqrt (1 + 4 * (ni ./ abs (D)) .^ 2)) .* (x > x_pi_jun) + ...
     ni .* (x > x_pi_jun) .* (x <= x_in_jun) + ...
     ni ^ 2 ./ (abs (D) / 2 .* (1 + sqrt (1 + 4 * (ni ./ abs (D)) .^2))) .* (x <= x_in_jun);
 
 V = Fn + Vth * log (n / ni);
 
 % scaling factors
 xbar = L;                         % [m]
 nbar = norm(D, 'inf');            % [m^{-3}]
 Vbar = Vth;                       % [V]
 mubar = max(u0n, u0p);            % [m^2 V^{-1} s^{-1}]
 tbar = xbar^2/(mubar*Vbar);       % [s]
 Rbar = nbar/tbar;                 % [m^{-3} s^{-1}]
 Ebar = Vbar/xbar;                 % [V m^{-1}]
 Jbar = q*mubar*nbar*Ebar;         % [A m^{-2}]
 CAubar = Rbar/nbar^3;             % [m^6 s^{-1}]
 abar = xbar^(-1);                 % [m^{-1}]
 
 % scaling procedure
 l2 = e0*Vbar/(q*nbar*xbar^2);    
 theta = ni/nbar;  
 
 xin = x/xbar;
 Din = D/nbar;
 pin = p/nbar;
 nin = n/nbar;
 Vin = V/Vbar;
 Fnin = Vin - log(nin);
 Fpin = Vin + log(pin);
 
 tnin = tn/tbar;
 tpin = tp/tbar;
 
 u0nin = u0n/mubar;
 u0pin = u0p/mubar;
 
 Cnin     = Cn/CAubar;
 Cpin     = Cp/CAubar;
 
 % tolerances for convergence checks
 ptoll = 1e-12;
 pmaxit = 1000;
 
 
 % solve the problem using the full DD model
 
 [n, p, V, Fn, Fp, Jn, Jp, it, res] = secs1d_dd_newton (xin, Din, Vin, nin, 
                                                        pin, l2, er, u0nin,            
                                                        u0pin, theta, tnin, tpin, 
                                                        Cnin, Cpin, ptoll, pmaxit)

 
 % Descaling procedure
 n    = nout*nbar;
 p    = pout*nbar;
 V    = Vout*Vbar;
 Fn   = V - Vth*log(n/ni);
 Fp   = V + Vth*log(p/ni);
 dV   = diff(V);
 dx   = diff(x);
 E    = -dV./dx;
 
 % compute current densities 
 [Bp, Bm] = bimu_bernoulli (dV/Vth);
 Jn       =  q*u0n*Vth .* (n(2:end) .* Bp - n(1:end-1) .* Bm) ./ dx; 
 Jp       = -q*u0p*Vth .* (p(2:end) .* Bm - p(1:end-1) .* Bp) ./ dx;
 Jtot     =  Jn+Jp;
 
 % band structure
 Efn  = -Fn;
 Efp  = -Fp;
 Ec   = Vth*log(Nc./n)+Efn;
 Ev   = -Vth*log(Nv./p)+Efp;

 plot (x, Efn, x, Efp, x, Ec, x, Ev)
 legend ('Efn', 'Efp', 'Ec', 'Ev')
 axis tight
 
 figure(2)
 plot (x, Vout*Vbar)

 figure(3)
 plot(x, pout*nbar)
 
 figure(4)
 plot(x, nout*nbar)
 
 axis tight