Qucs-core  0.0.19
diode.cpp
Go to the documentation of this file.
00001 /*
00002  * diode.cpp - diode class implementation
00003  *
00004  * Copyright (C) 2004, 2005, 2006, 2007, 2008 Stefan Jahn <stefan@lkcc.org>
00005  *
00006  * This is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2, or (at your option)
00009  * any later version.
00010  *
00011  * This software is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this package; see the file COPYING.  If not, write to
00018  * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
00019  * Boston, MA 02110-1301, USA.
00020  *
00021  * $Id$
00022  *
00023  */
00024 
00025 #if HAVE_CONFIG_H
00026 # include <config.h>
00027 #endif
00028 
00029 #include "component.h"
00030 #include "device.h"
00031 #include "devstates.h"
00032 #include "diode.h"
00033 
00034 #define NODE_C 0 /* cathode node */
00035 #define NODE_A 1 /* anode node   */
00036 
00037 #define StateVars 1 // state variables
00038 
00039 // state variable indices
00040 #define _UdPrev 0
00041 
00042 // state variable shortcuts
00043 #define UdPrev deviceVar (_UdPrev)
00044 
00045 using namespace qucs;
00046 using namespace qucs::device;
00047 
00048 // Constructor for the diode.
00049 diode::diode () : circuit (2) {
00050   rs = NULL;
00051   type = CIR_DIODE;
00052 }
00053 
00054 // Callback for S-parameter analysis.
00055 void diode::calcSP (nr_double_t frequency) {
00056   nr_double_t gd = getOperatingPoint ("gd");
00057   nr_double_t Cd = getOperatingPoint ("Cd");
00058   nr_complex_t y = 2 * z0 * nr_complex_t (gd, Cd * 2.0 * pi * frequency);
00059   setS (NODE_C, NODE_C, 1.0 / (1.0 + y));
00060   setS (NODE_A, NODE_A, 1.0 / (1.0 + y));
00061   setS (NODE_C, NODE_A, y / (1.0 + y));
00062   setS (NODE_A, NODE_C, y / (1.0 + y));
00063 }
00064 
00065 // Callback for S-parameter noise analysis.
00066 void diode::calcNoiseSP (nr_double_t frequency) {
00067 #if MICHAEL /* shot noise only */
00068   nr_double_t Id = getOperatingPoint ("Id");
00069   nr_double_t Is = getPropertyDouble ("Is") + getPropertyDouble ("Isr");
00070 
00071   // adjust shot noise current if necessary
00072   if (Id < -Is) Id = -Is;
00073 
00074   nr_double_t gd = getOperatingPoint ("gd");
00075   nr_double_t Cd = getOperatingPoint ("Cd");
00076 
00077   nr_complex_t y = rect (gd, Cd * 2.0 * pi * frequency);
00078   nr_complex_t f = 2 * z0 * (Id + 2 * Is) / norm (2 * z0 * y + 1) * QoverkB / T0;
00079   setN (NODE_C, NODE_C, +f); setN (NODE_A, NODE_A, +f);
00080   setN (NODE_C, NODE_A, -f); setN (NODE_A, NODE_C, -f);
00081 #else
00082   setMatrixN (cytocs (calcMatrixCy (frequency) * z0, getMatrixS ()));
00083 #endif
00084 }
00085 
00086 // Computes noise correlation matrix Cy.
00087 matrix diode::calcMatrixCy (nr_double_t frequency) {
00088   // fetch computed operating points
00089   nr_double_t Id = getOperatingPoint ("Id");
00090   nr_double_t Is = getPropertyDouble ("Is") + getPropertyDouble ("Isr");
00091 
00092   // adjust shot noise current if necessary
00093   if (Id < -Is) Id = -Is;
00094 
00095   nr_double_t Kf  = getPropertyDouble ("Kf");
00096   nr_double_t Af  = getPropertyDouble ("Af");
00097   nr_double_t Ffe = getPropertyDouble ("Ffe");
00098 
00099   // build noise current correlation matrix
00100   matrix cy (2);
00101   nr_double_t i = 2 * (Id + 2 * Is) * QoverkB / T0 +    // shot noise
00102     Kf * qucs::pow (fabs (Id), Af) / qucs::pow (frequency, Ffe) / kB / T0; // flicker noise
00103   cy.set (NODE_C, NODE_C, +i); cy.set (NODE_A, NODE_A, +i);
00104   cy.set (NODE_A, NODE_C, -i); cy.set (NODE_C, NODE_A, -i);
00105   return cy;
00106 }
00107 
00108 // Initializes the diode model including temperature and area effects.
00109 void diode::initModel (void) {
00110   // fetch necessary device properties
00111   nr_double_t T  = getPropertyDouble ("Temp");
00112   nr_double_t Tn = getPropertyDouble ("Tnom");
00113   nr_double_t A  = getPropertyDouble ("Area");
00114 
00115   // compute Is temperature and area dependency
00116   nr_double_t Is  = getPropertyDouble ("Is");
00117   nr_double_t N   = getPropertyDouble ("N");
00118   nr_double_t Xti = getPropertyDouble ("Xti");
00119   nr_double_t Eg  = getPropertyDouble ("Eg");
00120   nr_double_t T1, T2;
00121   T2 = celsius2kelvin (T);
00122   T1 = celsius2kelvin (Tn);
00123   Is = pnCurrent_T (T1, T2, Is, Eg, N, Xti);
00124   setScaledProperty ("Is", Is * A);
00125 
00126   // compute Isr temperature and area dependency
00127   nr_double_t Isr = getPropertyDouble ("Isr");
00128   nr_double_t Nr  = getPropertyDouble ("Nr");
00129   Isr = pnCurrent_T (T1, T2, Isr, Eg, Nr, Xti);
00130   setScaledProperty ("Isr", Isr * A);
00131 
00132   // check unphysical parameters
00133   if (Nr < 1.0) {
00134     logprint (LOG_ERROR, "WARNING: Unphysical model parameter Nr = %g in "
00135               "diode `%s'\n", Nr, getName ());
00136   }
00137   if (N < 1.0) {
00138     logprint (LOG_ERROR, "WARNING: Unphysical model parameter N = %g in "
00139               "diode `%s'\n", N, getName ());
00140   }
00141 
00142   // compute Vj temperature dependency
00143   nr_double_t Vj = getPropertyDouble ("Vj");
00144   nr_double_t VjT;
00145   VjT = pnPotential_T (T1,T2, Vj);
00146   setScaledProperty ("Vj", VjT);
00147 
00148   // compute Cj0 temperature and area dependency
00149   nr_double_t Cj0 = getPropertyDouble ("Cj0");
00150   nr_double_t M   = getPropertyDouble ("M");
00151   Cj0 = pnCapacitance_T (T1, T2, M, VjT / Vj, Cj0);
00152   setScaledProperty ("Cj0", Cj0 * A);
00153 
00154   // check unphysical parameters
00155   if (M > 1.0) {
00156     logprint (LOG_ERROR, "WARNING: Unphysical model parameter M = %g in "
00157               "Diode `%s'\n", M, getName ());
00158   }
00159 
00160   // compute Bv temperature dependency
00161   nr_double_t Bv  = getPropertyDouble ("Bv");
00162   nr_double_t Tbv = getPropertyDouble ("Tbv");
00163   nr_double_t DT  = T2 - T1;
00164   Bv = Bv - Tbv * DT;
00165   setScaledProperty ("Bv", Bv);
00166 
00167   // compute Tt temperature dependency
00168   nr_double_t Tt   = getPropertyDouble ("Tt");
00169   nr_double_t Ttt1 = getPropertyDouble ("Ttt1");
00170   nr_double_t Ttt2 = getPropertyDouble ("Ttt2");
00171   Tt = Tt * (1 + Ttt1 * DT + Ttt2 * DT * DT);
00172   setScaledProperty ("Tt", Tt);
00173 
00174   // compute M temperature dependency
00175   nr_double_t Tm1 = getPropertyDouble ("Tm1");
00176   nr_double_t Tm2 = getPropertyDouble ("Tm2");
00177   M = M * (1 + Tm1 * DT + Tm2 * DT * DT);
00178   setScaledProperty ("M", M);
00179 
00180   // compute Rs temperature and area dependency
00181   nr_double_t Rs  = getPropertyDouble ("Rs");
00182   nr_double_t Trs = getPropertyDouble ("Trs");
00183   Rs = Rs * (1 + Trs * DT);
00184   setScaledProperty ("Rs", Rs / A);
00185 }
00186 
00187 // Prepares DC (i.e. HB) analysis.
00188 void diode::prepareDC (void) {
00189   // allocate MNA matrices
00190   allocMatrixMNA ();
00191 
00192   // initialize scalability
00193   initModel ();
00194 
00195   // initialize starting values
00196   Ud = real (getV (NODE_A) - getV (NODE_C));
00197   for (int i = 0; i < deviceStates (); i++) {
00198     deviceState (i);
00199     UdPrev = Ud;
00200   }
00201 
00202   // get device temperature
00203   nr_double_t T = getPropertyDouble ("Temp");
00204 
00205   // possibly insert series resistance
00206   nr_double_t Rs = getScaledProperty ("Rs");
00207   if (Rs != 0.0) {
00208     // create additional circuit if necessary and reassign nodes
00209     rs = splitResistor (this, rs, "Rs", "anode", NODE_A);
00210     rs->setProperty ("Temp", T);
00211     rs->setProperty ("R", Rs);
00212     rs->setProperty ("Controlled", getName ());
00213     rs->initDC ();
00214   }
00215   // no series resistance
00216   else {
00217     disableResistor (this, rs, NODE_A);
00218   }
00219 
00220   // calculate actual breakdown voltage
00221   Bv = getScaledProperty ("Bv");
00222   if (Bv != 0) {
00223     nr_double_t Ibv, Is, tol, Ut, Xbv, Xibv;
00224     Ibv = getPropertyDouble ("Ibv");
00225     Is = getScaledProperty ("Is");
00226     Ut = celsius2kelvin (T) * kBoverQ;
00227     // adjust very small breakdown currents
00228     if (Ibv < Is * Bv / Ut) {
00229       Ibv = Is * Bv / Ut;
00230       Xbv = Bv;
00231       logprint (LOG_ERROR, "WARNING: Increased breakdown current to %g to "
00232                 "match the saturation current %g\n", Ibv, Is);
00233     }
00234     // fit reverse and forward regions
00235     else {
00236       int good = 0;
00237       tol = 1e-3 * Ibv;
00238       Xbv = Bv - Ut * qucs::log (1 + Ibv / Is);
00239       for (int i = 0; i < 25 ; i++) {
00240         Xbv = Bv - Ut * qucs::log (Ibv / Is + 1 - Xbv / Ut);
00241         Xibv = Is * (qucs::exp ((Bv - Xbv) / Ut) - 1 + Xbv / Ut);
00242         if (fabs (Xibv - Ibv) < tol) {
00243           Bv = Xbv;
00244           good = 1;
00245           break;
00246         }
00247       }
00248       if (!good) {
00249         logprint (LOG_ERROR, "WARNING: Unable to fit reverse and forward "
00250                   "diode regions using Bv=%g and Ibv=%g\n", Bv, Ibv);
00251       }
00252     }
00253   }
00254 }
00255 
00256 // Callback for initializing the DC analysis.
00257 void diode::initDC (void) {
00258   deviceStates (StateVars, 1);
00259   doHB = false;
00260   prepareDC ();
00261 }
00262 
00263 // Callback for restarting the DC analysis.
00264 void diode::restartDC (void) {
00265   // apply starting value to previous iteration value
00266   UdPrev = real (getV (NODE_A) - getV (NODE_C));
00267 }
00268 
00269 // Callback for DC analysis.
00270 void diode::calcDC (void) {
00271   // get device properties
00272   nr_double_t Is  = getScaledProperty ("Is");
00273   nr_double_t N   = getPropertyDouble ("N");
00274   nr_double_t Isr = getScaledProperty ("Isr");
00275   nr_double_t Nr  = getPropertyDouble ("Nr");
00276   nr_double_t Ikf = getPropertyDouble ("Ikf");
00277   nr_double_t T   = getPropertyDouble ("Temp");
00278 
00279   nr_double_t Ut, Ieq, Ucrit, gtiny;
00280 
00281   T = celsius2kelvin (T);
00282   Ut = T * kBoverQ;
00283   Ud = real (getV (NODE_A) - getV (NODE_C));
00284 
00285   // critical voltage necessary for bad start values
00286   Ucrit = pnCriticalVoltage (Is, N * Ut);
00287   if (Bv != 0 && Ud < std::min (0.0, -Bv + 10 * N * Ut)) {
00288     nr_double_t V = -(Ud + Bv);
00289     V = pnVoltage (V, -(UdPrev + Bv), Ut * N, Ucrit);
00290     Ud = -(V + Bv);
00291   }
00292   else {
00293     Ud = pnVoltage (Ud, UdPrev, Ut * N, Ucrit);
00294   }
00295   UdPrev = Ud;
00296 
00297   // tiny derivative for little junction voltage
00298   gtiny = (Ud < - 10 * Ut * N && Bv != 0) ? (Is + Isr) : 0;
00299 
00300   if (Ud >= -3 * N * Ut) { // forward region
00301     gd = pnConductance (Ud, Is, Ut * N) + pnConductance (Ud, Isr, Ut * Nr);
00302     Id = pnCurrent (Ud, Is, Ut * N) + pnCurrent (Ud, Isr, Ut * Nr);
00303   }
00304   else if (Bv == 0 || Ud >= -Bv) { // reverse region
00305     nr_double_t a = 3 * N * Ut / (Ud * euler);
00306     a = cubic (a);
00307     Id = -Is * (1 + a);
00308     gd = +Is * 3 * a / Ud;
00309   }
00310   else { // middle region
00311     nr_double_t a = qucs::exp (-(Bv + Ud) / N / Ut);
00312     Id = -Is * a;
00313     gd = +Is * a / Ut / N;
00314   }
00315 
00316   // knee current calculations
00317   if (Ikf != 0.0) {
00318     nr_double_t a = Ikf / (Ikf + Id);
00319     gd *= 0.5 * (2 - Id * a / Ikf) * qucs::sqrt (a);
00320     Id *= qucs::sqrt (a);
00321   }
00322 
00323   Id += gtiny * Ud;
00324   gd += gtiny;
00325 
00326   // HB simulation
00327   if (doHB) {
00328     Ieq = Id;
00329     setGV (NODE_C, -gd * Ud);
00330     setGV (NODE_A, +gd * Ud);
00331   }
00332   // DC and transient simulation
00333   else {
00334     Ieq = Id - Ud * gd;
00335   }
00336 
00337   // fill in I-Vector
00338   setI (NODE_C, +Ieq);
00339   setI (NODE_A, -Ieq);
00340 
00341   // fill in G-Matrix
00342   setY (NODE_C, NODE_C, +gd); setY (NODE_A, NODE_A, +gd);
00343   setY (NODE_C, NODE_A, -gd); setY (NODE_A, NODE_C, -gd);
00344 }
00345 
00346 // Saves operating points (voltages).
00347 void diode::saveOperatingPoints (void) {
00348   nr_double_t Vd = real (getV (NODE_A) - getV (NODE_C));
00349   setOperatingPoint ("Vd", Vd);
00350 }
00351 
00352 // Loads operating points (voltages).
00353 void diode::loadOperatingPoints (void) {
00354   Ud = getOperatingPoint ("Vd");
00355 }
00356 
00357 // Calculates and saves operating points.
00358 void diode::calcOperatingPoints (void) {
00359 
00360   // load operating points
00361   loadOperatingPoints ();
00362 
00363   // get necessary properties
00364   nr_double_t M   = getScaledProperty ("M");
00365   nr_double_t Cj0 = getScaledProperty ("Cj0");
00366   nr_double_t Vj  = getScaledProperty ("Vj");
00367   nr_double_t Fc  = getPropertyDouble ("Fc");
00368   nr_double_t Cp  = getPropertyDouble ("Cp");
00369   nr_double_t Tt  = getScaledProperty ("Tt");
00370 
00371   // calculate capacitances and charges
00372   nr_double_t Cd;
00373   Cd = pnCapacitance (Ud, Cj0, Vj, M, Fc) + Tt * gd + Cp;
00374   Qd = pnCharge (Ud, Cj0, Vj, M, Fc) + Tt * Id + Cp * Ud;
00375 
00376   // save operating points
00377   setOperatingPoint ("gd", gd);
00378   setOperatingPoint ("Id", Id);
00379   setOperatingPoint ("Cd", Cd);
00380 }
00381 
00382 // Callback for initializing the AC analysis.
00383 void diode::initAC (void) {
00384   allocMatrixMNA ();
00385 }
00386 
00387 // Callback for the AC analysis.
00388 void diode::calcAC (nr_double_t frequency) {
00389   nr_double_t gd = getOperatingPoint ("gd");
00390   nr_double_t Cd = getOperatingPoint ("Cd");
00391   nr_complex_t y = nr_complex_t (gd, Cd * 2.0 * pi * frequency);
00392   setY (NODE_C, NODE_C, +y); setY (NODE_A, NODE_A, +y);
00393   setY (NODE_C, NODE_A, -y); setY (NODE_A, NODE_C, -y);
00394 }
00395 
00396 // Callback for the AC noise analysis.
00397 void diode::calcNoiseAC (nr_double_t frequency) {
00398   setMatrixN (calcMatrixCy (frequency));
00399 }
00400 
00401 #define qState 0 // charge state
00402 #define cState 1 // current state
00403 
00404 // Callback for initializing the TR analysis.
00405 void diode::initTR (void) {
00406   setStates (2);
00407   initDC ();
00408 }
00409 
00410 // Callback for the TR analysis.
00411 void diode::calcTR (nr_double_t) {
00412   calcDC ();
00413   saveOperatingPoints ();
00414   calcOperatingPoints ();
00415 
00416   nr_double_t Cd = getOperatingPoint ("Cd");
00417 
00418   transientCapacitance (qState, NODE_A, NODE_C, Cd, Ud, Qd);
00419 }
00420 
00421 // Callback for initializing the HB analysis.
00422 void diode::initHB (int frequencies) {
00423   deviceStates (StateVars, frequencies);
00424   doHB = true;
00425   prepareDC ();
00426   allocMatrixHB ();
00427 }
00428 
00429 // Callback for the HB analysis.
00430 void diode::calcHB (int frequency) {
00431   // set current frequency state
00432   deviceState (frequency);
00433 
00434   // g's (dI/dU) into Y-Matrix and I's into I-Vector
00435   calcDC ();
00436 
00437   // calculate Q and C
00438   saveOperatingPoints ();
00439   calcOperatingPoints ();
00440 
00441   nr_double_t Cd = getOperatingPoint ("Cd");
00442 
00443   // fill in Q's in Q-Vector
00444   setQ (NODE_C, +Qd);
00445   setQ (NODE_A, -Qd);
00446 
00447   setCV (NODE_C, -Cd * Ud);
00448   setCV (NODE_A, +Cd * Ud);
00449 
00450   // fill in C's (dQ/dU) into QV-Matrix
00451   setQV (NODE_C, NODE_C, +Cd); setQV (NODE_A, NODE_A, +Cd);
00452   setQV (NODE_C, NODE_A, -Cd); setQV (NODE_A, NODE_C, -Cd);
00453 }
00454 
00455 // properties
00456 PROP_REQ [] = {
00457   { "Is", PROP_REAL, { 1e-15, PROP_NO_STR }, PROP_POS_RANGE },
00458   { "N", PROP_REAL, { 1, PROP_NO_STR }, PROP_RNGII (1e-6, 100) },
00459   { "M", PROP_REAL, { 0.5, PROP_NO_STR }, PROP_RNGII (0, 2) },
00460   { "Cj0", PROP_REAL, { 10e-15, PROP_NO_STR }, PROP_POS_RANGE },
00461   { "Vj", PROP_REAL, { 0.7, PROP_NO_STR }, PROP_RNGXI (0, 10) },
00462   PROP_NO_PROP };
00463 PROP_OPT [] = {
00464   { "Rs", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00465   { "Isr", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00466   { "Nr", PROP_REAL, { 2, PROP_NO_STR }, PROP_RNGII (0.1, 100) },
00467   { "Bv", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00468   { "Ibv", PROP_REAL, { 1e-3, PROP_NO_STR }, PROP_POS_RANGE },
00469   { "Ikf", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00470   { "Tt", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00471   { "Fc", PROP_REAL, { 0.5, PROP_NO_STR }, PROP_RNGIX (0, 1) },
00472   { "Cp", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00473   { "Kf", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00474   { "Af", PROP_REAL, { 1, PROP_NO_STR }, PROP_POS_RANGE },
00475   { "Ffe", PROP_REAL, { 1, PROP_NO_STR }, PROP_POS_RANGE },
00476   { "Temp", PROP_REAL, { 26.85, PROP_NO_STR }, PROP_MIN_VAL (K) },
00477   { "Xti", PROP_REAL, { 3, PROP_NO_STR }, PROP_POS_RANGE },
00478   { "Eg", PROP_REAL, { EgSi, PROP_NO_STR }, PROP_POS_RANGE },
00479   { "Tbv", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00480   { "Trs", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
00481   { "Ttt1", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
00482   { "Ttt2", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
00483   { "Tm1", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
00484   { "Tm2", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
00485   { "Tnom", PROP_REAL, { 26.85, PROP_NO_STR }, PROP_MIN_VAL (K) },
00486   { "Area", PROP_REAL, { 1, PROP_NO_STR }, PROP_POS_RANGEX },
00487   PROP_NO_PROP };
00488 struct define_t diode::cirdef =
00489   { "Diode", 2, PROP_COMPONENT, PROP_NO_SUBSTRATE, PROP_NONLINEAR, PROP_DEF };