Difference between revisions of "Secs1d package"

From Octave
Jump to navigation Jump to search
(The SECS1D Package)
 
m (Rename "Octave-Forge" to "Octave Forge" (https://lists.gnu.org/archive/html/octave-maintainers/2018-08/msg00138.html).)
 
Line 122: Line 122:
 
</code>
 
</code>
  
[[Category:Octave-Forge]]
+
[[Category:Octave Forge]]

Latest revision as of 11:25, 2 September 2018

  • 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