Secs1d package

• Example, a PIN diode

% physical constants and parameters
 secs1d_physical_constants;
 secs1d_silicon_material_properties;
 
 % geometry
 L  = 10e-6;          % [m] 
 xm = L/2;
 
 Nelements = 1000;
 x         = linspace (0, L, Nelements+1)';
 sinodes   = [1:length(x)];
 
 % dielectric constant (silicon)
 er = esir * ones (Nelements, 1);
 
 % doping profile [m^{-3}]
 Na = 1e23 * (x <= xm);
 Nd = 1e23 * (x > xm);
 
 % avoid zero doping
 D  = Nd - Na;  
  
 % initial guess for n, p, V, phin, phip
 V_p = -1;
 V_n =  0;
 
 Fp = V_p * (x <= xm);
 Fn = Fp;
 
 p = abs (D) / 2 .* (1 + sqrt (1 + 4 * (ni./abs(D)) .^2)) .* (x <= xm) + ...
     ni^2 ./ (abs (D) / 2 .* (1 + sqrt (1 + 4 * (ni ./ abs (D)) .^2))) .* (x > xm);
 
 n = abs (D) / 2 .* (1 + sqrt (1 + 4 * (ni ./ abs (D)) .^ 2)) .* (x > xm) + ...
     ni ^ 2 ./ (abs (D) / 2 .* (1 + sqrt (1 + 4 * (ni ./ abs (D)) .^2))) .* (x <= xm);
 
 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 = 1/xbar;                  % [m^{-1}]
 
 % scaling procedure
 l2 = e0 * Vbar / (q * nbar * xbar^2);
 theta = ni / nbar;
 
 xin = x / xbar;
 Din = D / nbar;
 Nain = Na / nbar;
 Ndin = Nd / 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;
 uminnin = uminn / mubar;
 vsatnin = vsatn / (mubar * Ebar);
 
 u0pin = u0p / mubar;
 uminpin = uminp / mubar;
 vsatpin = vsatp / (mubar * Ebar);
 
 Nrefnin = Nrefn / nbar;
 Nrefpin = Nrefp / nbar;
 
 Cnin     = Cn / CAubar;
 Cpin     = Cp / CAubar;
 
 anin     = an / abar;
 apin     = ap / abar;
 Ecritnin = Ecritn / Ebar;
 Ecritpin = Ecritp / Ebar;
 
 % tolerances for convergence checks
 toll  = 1e-3;
 maxit = 1000;
 ptoll = 1e-12;
 pmaxit = 1000;
 
 % solve the problem using the full DD model
 [nout, pout, Vout, Fnout, Fpout, Jnout, Jpout, it, res] = ...
       secs1d_dd_gummel_map (xin, Din, Nain, Ndin, pin, nin, Vin, Fnin, Fpin, ...
                             l2, er, u0nin, uminnin, vsatnin, betan, Nrefnin, ...
 	                       u0pin, uminpin, vsatpin, betap, Nrefpin, theta, ...
 		               tnin, tpin, Cnin, Cpin, anin, apin, ...
 		               Ecritnin, Ecritpin, toll, maxit, 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;
 
 % 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