Qucs-core  0.0.19
bjt.cpp
Go to the documentation of this file.
00001 /*
00002  * bjt.cpp - bipolar junction transistor class implementation
00003  *
00004  * Copyright (C) 2004, 2005, 2006, 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 "bjt.h"
00032 
00033 #define NEWSGP 0
00034 
00035 #define NODE_B 0 /* base node       */
00036 #define NODE_C 1 /* collector node  */
00037 #define NODE_E 2 /* emitter node    */
00038 #define NODE_S 3 /* substrate node  */
00039 
00040 using namespace qucs;
00041 using namespace qucs::device;
00042 
00043 bjt::bjt () : circuit (4) {
00044   cbcx = rb = re = rc = NULL;
00045   type = CIR_BJT;
00046 }
00047 
00048 void bjt::calcSP (nr_double_t frequency) {
00049   // build admittance matrix and convert it to S-parameter matrix
00050   setMatrixS (ytos (calcMatrixY (frequency)));
00051 }
00052 
00053 matrix bjt::calcMatrixY (nr_double_t frequency) {
00054 
00055   // fetch computed operating points
00056   nr_double_t Cbe  = getOperatingPoint ("Cbe");
00057   nr_double_t gbe  = getOperatingPoint ("gpi");
00058   nr_double_t Cbci = getOperatingPoint ("Cbci");
00059   nr_double_t gbc  = getOperatingPoint ("gmu");
00060   nr_double_t Ccs  = getOperatingPoint ("Ccs");
00061 #if NEWSGP
00062   nr_double_t gm   = getOperatingPoint ("gmf");
00063   nr_double_t gmr  = getOperatingPoint ("gmr");
00064 #else
00065   nr_double_t gm   = getOperatingPoint ("gm");
00066   nr_double_t go   = getOperatingPoint ("go");
00067 #endif
00068   nr_double_t Ptf  = getPropertyDouble ("Ptf");
00069   nr_double_t Tf   = getPropertyDouble ("Tf");
00070 
00071   // compute admittance matrix entries
00072   nr_complex_t Ybe = nr_complex_t (gbe, 2.0 * pi * frequency * Cbe);
00073   nr_complex_t Ybc = nr_complex_t (gbc, 2.0 * pi * frequency * Cbci);
00074   nr_complex_t Ycs = nr_complex_t (0.0, 2.0 * pi * frequency * Ccs);
00075 
00076   // admittance matrix entries for "transcapacitance"
00077   nr_complex_t Ybebc = nr_complex_t (0.0, 2.0 * pi * frequency * dQbedUbc);
00078 
00079   // compute influence of excess phase
00080   nr_double_t phase = deg2rad (Ptf) * Tf * 2 * pi * frequency;
00081 #if NEWSGP
00082   nr_complex_t gmf = qucs::polar (gm, -phase);
00083 #else
00084   nr_complex_t gmf = qucs::polar (gm + go, -phase) - go;
00085 #endif
00086 
00087   // build admittance matrix
00088   matrix y (4);
00089 #if NEWSGP
00090   // for some reason this small signal equivalent can't be used
00091   y.set (NODE_B, NODE_B, Ybc + Ybe + Ybebc);
00092   y.set (NODE_B, NODE_C, -Ybc - Ybebc);
00093   y.set (NODE_B, NODE_E, -Ybe);
00094   y.set (NODE_B, NODE_S, 0);
00095   y.set (NODE_C, NODE_B, -Ybc + gmf + gmr);
00096   y.set (NODE_C, NODE_C, Ybc - gmr + Ycs);
00097   y.set (NODE_C, NODE_E, -gmf);
00098   y.set (NODE_C, NODE_S, -Ycs);
00099   y.set (NODE_E, NODE_B, -Ybe - gmf - gmr - Ybebc);
00100   y.set (NODE_E, NODE_C, gmr + Ybebc);
00101   y.set (NODE_E, NODE_E, Ybe + gmf);
00102   y.set (NODE_E, NODE_S, 0);
00103   y.set (NODE_S, NODE_B, 0);
00104   y.set (NODE_S, NODE_C, -Ycs);
00105   y.set (NODE_S, NODE_E, 0);
00106   y.set (NODE_S, NODE_S, Ycs);
00107 #else /* !NEWSGP */
00108   y.set (NODE_B, NODE_B, Ybc + Ybe + Ybebc);
00109   y.set (NODE_B, NODE_C, -Ybc - Ybebc);
00110   y.set (NODE_B, NODE_E, -Ybe);
00111   y.set (NODE_B, NODE_S, 0);
00112   y.set (NODE_C, NODE_B, -Ybc + gmf);
00113   y.set (NODE_C, NODE_C, Ybc + Ycs + go);
00114   y.set (NODE_C, NODE_E, -gmf - go);
00115   y.set (NODE_C, NODE_S, -Ycs);
00116   y.set (NODE_E, NODE_B, -Ybe - gmf - Ybebc);
00117   y.set (NODE_E, NODE_C, -go + Ybebc);
00118   y.set (NODE_E, NODE_E, Ybe + gmf + go);
00119   y.set (NODE_E, NODE_S, 0);
00120   y.set (NODE_S, NODE_B, 0);
00121   y.set (NODE_S, NODE_C, -Ycs);
00122   y.set (NODE_S, NODE_E, 0);
00123   y.set (NODE_S, NODE_S, Ycs);
00124 #endif /* !NEWSGP */
00125   return y;
00126 }
00127 
00128 void bjt::calcNoiseSP (nr_double_t frequency) {
00129   setMatrixN (cytocs (calcMatrixCy (frequency) * z0, getMatrixS ()));
00130 }
00131 
00132 matrix bjt::calcMatrixCy (nr_double_t frequency) {
00133 
00134   // fetch computed operating points
00135   nr_double_t Ibe = fabs (getOperatingPoint ("Ibe"));
00136   nr_double_t Ice = fabs (getOperatingPoint ("Ice"));
00137 
00138   // get model properties
00139   nr_double_t Kf  = getPropertyDouble ("Kf");
00140   nr_double_t Af  = getPropertyDouble ("Af");
00141   nr_double_t Ffe = getPropertyDouble ("Ffe");
00142   nr_double_t Kb  = getPropertyDouble ("Kb");
00143   nr_double_t Ab  = getPropertyDouble ("Ab");
00144   nr_double_t Fb  = getPropertyDouble ("Fb");
00145 
00146   nr_double_t ib = 2 * Ibe * QoverkB / T0 +            // shot noise
00147     (Kf * qucs::pow (Ibe, Af) / qucs::pow (frequency, Ffe) +       // flicker noise
00148      Kb * qucs::pow (Ibe, Ab) / (1 + sqr (frequency / Fb)))  // burst noise
00149     / kB / T0;
00150   nr_double_t ic = 2 * Ice * QoverkB / T0;             // shot noise
00151 
00152   /* build noise current correlation matrix and convert it to
00153      noise-wave correlation matrix */
00154   matrix cy = matrix (4);
00155   cy.set (NODE_B, NODE_B, ib);
00156   cy.set (NODE_B, NODE_E, -ib);
00157   cy.set (NODE_C, NODE_C, ic);
00158   cy.set (NODE_C, NODE_E, -ic);
00159   cy.set (NODE_E, NODE_B, -ib);
00160   cy.set (NODE_E, NODE_C, -ic);
00161   cy.set (NODE_E, NODE_E, ic + ib);
00162   return cy;
00163 }
00164 
00165 void bjt::initModel (void) {
00166   // fetch necessary device properties
00167   nr_double_t T  = getPropertyDouble ("Temp");
00168   nr_double_t Tn = getPropertyDouble ("Tnom");
00169   nr_double_t A  = getPropertyDouble ("Area");
00170 
00171   // compute Is temperature and area dependency
00172   nr_double_t Is  = getPropertyDouble ("Is");
00173   nr_double_t Xti = getPropertyDouble ("Xti");
00174   nr_double_t Eg  = getPropertyDouble ("Eg");
00175   nr_double_t T1, T2, IsT;
00176   T2 = celsius2kelvin (T);
00177   T1 = celsius2kelvin (Tn);
00178   IsT = pnCurrent_T (T1, T2, Is, Eg, 1, Xti);
00179   setScaledProperty ("Is", IsT * A);
00180 
00181   // compute Vje, Vjc and Vjs temperature dependencies
00182   nr_double_t Vje = getPropertyDouble ("Vje");
00183   nr_double_t Vjc = getPropertyDouble ("Vjc");
00184   nr_double_t Vjs = getPropertyDouble ("Vjs");
00185   nr_double_t VjeT, VjcT, VjsT;
00186   VjeT = pnPotential_T (T1,T2, Vje);
00187   VjcT = pnPotential_T (T1,T2, Vjc);
00188   VjsT = pnPotential_T (T1,T2, Vjs);
00189   setScaledProperty ("Vje", VjeT);
00190   setScaledProperty ("Vjc", VjcT);
00191   setScaledProperty ("Vjs", VjsT);
00192 
00193   // compute Bf and Br temperature dependencies
00194   nr_double_t Bf  = getPropertyDouble ("Bf");
00195   nr_double_t Br  = getPropertyDouble ("Br");
00196   nr_double_t Xtb = getPropertyDouble ("Xtb");
00197   nr_double_t F = qucs::exp (Xtb * qucs::log (T2 / T1));
00198   setScaledProperty ("Bf", Bf * F);
00199   setScaledProperty ("Br", Br * F);
00200 
00201   // compute Ise and Isc temperature and area dependencies
00202   nr_double_t Ise = getPropertyDouble ("Ise");
00203   nr_double_t Isc = getPropertyDouble ("Isc");
00204   nr_double_t Ne  = getPropertyDouble ("Ne");
00205   nr_double_t Nc  = getPropertyDouble ("Nc");
00206   nr_double_t G = qucs::log (IsT / Is);
00207   nr_double_t F1 = qucs::exp (G / Ne);
00208   nr_double_t F2 = qucs::exp (G / Nc);
00209   Ise = Ise / F * F1;
00210   Isc = Isc / F * F2;
00211   setScaledProperty ("Ise", Ise * A);
00212   setScaledProperty ("Isc", Isc * A);
00213 
00214   // check unphysical parameters
00215   nr_double_t Nf = getPropertyDouble ("Nf");
00216   nr_double_t Nr = getPropertyDouble ("Nr");
00217   if (Nf < 1.0) {
00218     logprint (LOG_ERROR, "WARNING: Unphysical model parameter Nf = %g in "
00219               "BJT `%s'\n", Nf, getName ());
00220   }
00221   if (Nr < 1.0) {
00222     logprint (LOG_ERROR, "WARNING: Unphysical model parameter Nr = %g in "
00223               "BJT `%s'\n", Nr, getName ());
00224   }
00225   if (Ne < 1.0) {
00226     logprint (LOG_ERROR, "WARNING: Unphysical model parameter Ne = %g in "
00227               "BJT `%s'\n", Ne, getName ());
00228   }
00229   if (Nc < 1.0) {
00230     logprint (LOG_ERROR, "WARNING: Unphysical model parameter Nc = %g in "
00231               "BJT `%s'\n", Nc, getName ());
00232   }
00233 
00234   /* Originally Vtf was expected to be PROP_POS_RANGE, but there are models
00235    * which use negative values. Instead of failing, warn the user.
00236    * \todo Provide a way to silece such warnings
00237    */
00238   nr_double_t Vtf = getPropertyDouble ("Vtf");
00239   if (Vtf < 0.0) {
00240     logprint (LOG_ERROR, "WARNING: Unphysical model parameter Vtf = %g in "
00241               "BJT `%s'\n", Vtf, getName ());
00242   }
00243 
00244   // compute Cje, Cjc and Cjs temperature and area dependencies
00245   nr_double_t Cje = getPropertyDouble ("Cje");
00246   nr_double_t Cjc = getPropertyDouble ("Cjc");
00247   nr_double_t Cjs = getPropertyDouble ("Cjs");
00248   nr_double_t Mje = getPropertyDouble ("Mje");
00249   nr_double_t Mjc = getPropertyDouble ("Mjc");
00250   nr_double_t Mjs = getPropertyDouble ("Mjs");
00251   Cje = pnCapacitance_T (T1, T2, Mje, VjeT / Vje, Cje);
00252   Cjc = pnCapacitance_T (T1, T2, Mjc, VjcT / Vjc, Cjc);
00253   Cjs = pnCapacitance_T (T1, T2, Mjs, VjsT / Vjs, Cjs);
00254   setScaledProperty ("Cje", Cje * A);
00255   setScaledProperty ("Cjc", Cjc * A);
00256   setScaledProperty ("Cjs", Cjs * A);
00257 
00258   // compute Rb, Rc, Re and Rbm area dependencies
00259   nr_double_t Rb  = getPropertyDouble ("Rb");
00260   nr_double_t Re  = getPropertyDouble ("Re");
00261   nr_double_t Rc  = getPropertyDouble ("Rc");
00262   nr_double_t Rbm = getPropertyDouble ("Rbm");
00263   setScaledProperty ("Rb", Rb / A);
00264   setScaledProperty ("Re", Re / A);
00265   setScaledProperty ("Rc", Rc / A);
00266   setScaledProperty ("Rbm", Rbm / A);
00267 
00268   // compute Ikf, Ikr, Irb and Itf area dependencies
00269   nr_double_t Ikf = getPropertyDouble ("Ikf");
00270   nr_double_t Ikr = getPropertyDouble ("Ikr");
00271   nr_double_t Irb = getPropertyDouble ("Irb");
00272   nr_double_t Itf = getPropertyDouble ("Itf");
00273   setScaledProperty ("Ikf", Ikf * A);
00274   setScaledProperty ("Ikr", Ikr * A);
00275   setScaledProperty ("Irb", Irb * A);
00276   setScaledProperty ("Itf", Itf * A);
00277 }
00278 
00279 void bjt::initDC (void) {
00280 
00281   // no transient analysis
00282   doTR = false;
00283 
00284   // allocate MNA matrices
00285   allocMatrixMNA ();
00286 
00287   // initialize scalability
00288   initModel ();
00289 
00290   // apply polarity of BJT
00291   const char * const type = getPropertyString ("Type");
00292   pol = !strcmp (type, "pnp") ? -1 : 1;
00293 
00294   // get simulation temperature
00295   nr_double_t T = getPropertyDouble ("Temp");
00296 
00297   // initialize starting values
00298   restartDC ();
00299 
00300   // disable additional base-collector capacitance
00301   if (deviceEnabled (cbcx)) {
00302     disableCapacitor (this, cbcx);
00303   }
00304 
00305   // possibly insert series resistance at emitter
00306   nr_double_t Re = getScaledProperty ("Re");
00307   if (Re != 0.0) {
00308     // create additional circuit if necessary and reassign nodes
00309     re = splitResistor (this, re, "Re", "emitter", NODE_E);
00310     re->setProperty ("R", Re);
00311     re->setProperty ("Temp", T);
00312     re->setProperty ("Controlled", getName ());
00313     re->initDC ();
00314   }
00315   // no series resistance at emitter
00316   else {
00317     disableResistor (this, re, NODE_E);
00318   }
00319 
00320   // possibly insert series resistance at collector
00321   nr_double_t Rc = getScaledProperty ("Rc");
00322   if (Rc != 0.0) {
00323     // create additional circuit if necessary and reassign nodes
00324     rc = splitResistor (this, rc, "Rc", "collector", NODE_C);
00325     rc->setProperty ("R", Rc);
00326     rc->setProperty ("Temp", T);
00327     rc->setProperty ("Controlled", getName ());
00328     rc->initDC ();
00329   }
00330   // no series resistance at collector
00331   else {
00332     disableResistor (this, rc, NODE_C);
00333   }
00334 
00335   // possibly insert base series resistance
00336   nr_double_t Rb  = getScaledProperty ("Rb");
00337   nr_double_t Rbm = getScaledProperty ("Rbm");
00338   if (Rbm <= 0.0) Rbm = Rb; // Rbm defaults to Rb if zero
00339   if (Rb < Rbm)   Rbm = Rb; // Rbm must be less or equal Rb
00340   setScaledProperty ("Rbm", Rbm);
00341   if (Rbm != 0.0) {
00342     // create additional circuit and reassign nodes
00343     rb = splitResistor (this, rb, "Rbb", "base", NODE_B);
00344     rb->setProperty ("R", Rb);
00345     rb->setProperty ("Temp", T);
00346     rb->setProperty ("Controlled", getName ());
00347     rb->initDC ();
00348   }
00349   // no series resistance at base
00350   else {
00351     disableResistor (this, rb, NODE_B);
00352     Rbb = 0.0;                 // set this operating point
00353     setProperty ("Xcjc", 1.0); // other than 1 is senseless here
00354   }
00355 }
00356 
00357 void bjt::restartDC (void) {
00358   // apply starting values to previous iteration values
00359   UbePrev = real (getV (NODE_B) - getV (NODE_E)) * pol;
00360   UbcPrev = real (getV (NODE_B) - getV (NODE_C)) * pol;
00361 }
00362 
00363 #define cexState 6 // extra excess phase state
00364 
00365 void bjt::calcDC (void) {
00366 
00367   // fetch device model parameters
00368   nr_double_t Is   = getScaledProperty ("Is");
00369   nr_double_t Nf   = getPropertyDouble ("Nf");
00370   nr_double_t Nr   = getPropertyDouble ("Nr");
00371   nr_double_t Vaf  = getPropertyDouble ("Vaf");
00372   nr_double_t Var  = getPropertyDouble ("Var");
00373   nr_double_t Ikf  = getScaledProperty ("Ikf");
00374   nr_double_t Ikr  = getScaledProperty ("Ikr");
00375   nr_double_t Bf   = getScaledProperty ("Bf");
00376   nr_double_t Br   = getScaledProperty ("Br");
00377   nr_double_t Ise  = getScaledProperty ("Ise");
00378   nr_double_t Isc  = getScaledProperty ("Isc");
00379   nr_double_t Ne   = getPropertyDouble ("Ne");
00380   nr_double_t Nc   = getPropertyDouble ("Nc");
00381   nr_double_t Rb   = getScaledProperty ("Rb");
00382   nr_double_t Rbm  = getScaledProperty ("Rbm");
00383   nr_double_t Irb  = getScaledProperty ("Irb");
00384   nr_double_t T    = getPropertyDouble ("Temp");
00385 
00386   nr_double_t Ut, Q1, Q2;
00387   nr_double_t Iben, Ibcn, Ibei, Ibci, Ibc, gbe, gbc, gtiny;
00388   nr_double_t IeqB, IeqC, IeqE, IeqS, UbeCrit, UbcCrit;
00389   nr_double_t gm, go;
00390 
00391   // interpret zero as infinity for these model parameters
00392   Ikf = Ikf > 0 ? 1.0 / Ikf : 0;
00393   Ikr = Ikr > 0 ? 1.0 / Ikr : 0;
00394   Vaf = Vaf > 0 ? 1.0 / Vaf : 0;
00395   Var = Var > 0 ? 1.0 / Var : 0;
00396 
00397   T = celsius2kelvin (T);
00398   Ut = T * kBoverQ;
00399   Ube = real (getV (NODE_B) - getV (NODE_E)) * pol;
00400   Ubc = real (getV (NODE_B) - getV (NODE_C)) * pol;
00401 
00402   // critical voltage necessary for bad start values
00403   UbeCrit = pnCriticalVoltage (Is, Nf * Ut);
00404   UbcCrit = pnCriticalVoltage (Is, Nr * Ut);
00405   UbePrev = Ube = pnVoltage (Ube, UbePrev, Ut * Nf, UbeCrit);
00406   UbcPrev = Ubc = pnVoltage (Ubc, UbcPrev, Ut * Nr, UbcCrit);
00407 
00408   Uce = Ube - Ubc;
00409 
00410   // base-emitter diodes
00411   gtiny = Ube < - 10 * Ut * Nf ? (Is + Ise) : 0;
00412 #if 0
00413   If = pnCurrent (Ube, Is, Ut * Nf);
00414   Ibei = If / Bf;
00415   gif = pnConductance (Ube, Is, Ut * Nf);
00416   gbei = gif / Bf;
00417   Iben = pnCurrent (Ube, Ise, Ut * Ne);
00418   gben = pnConductance (Ube, Ise, Ut * Ne);
00419   Ibe = Ibei + Iben + gtiny * Ube;
00420   gbe = gbei + gben + gtiny;
00421 #else
00422   pnJunctionBIP (Ube, Is, Ut * Nf, If, gif);
00423   Ibei = If / Bf;
00424   gbei = gif / Bf;
00425   pnJunctionBIP (Ube, Ise, Ut * Ne, Iben, gben);
00426   Iben += gtiny * Ube;
00427   gben += gtiny;
00428   Ibe = Ibei + Iben;
00429   gbe = gbei + gben;
00430 #endif
00431 
00432   // base-collector diodes
00433   gtiny = Ubc < - 10 * Ut * Nr ? (Is + Isc) : 0;
00434 #if 0
00435   Ir = pnCurrent (Ubc, Is, Ut * Nr);
00436   Ibci = Ir / Br;
00437   gir = pnConductance (Ubc, Is, Ut * Nr);
00438   gbci = gir / Br;
00439   Ibcn = pnCurrent (Ubc, Isc, Ut * Nc);
00440   gbcn = pnConductance (Ubc, Isc, Ut * Nc);
00441   Ibc = Ibci + Ibcn + gtiny * Ubc;
00442   gbc = gbci + gbcn + gtiny;
00443 #else
00444   pnJunctionBIP (Ubc, Is, Ut * Nr, Ir, gir);
00445   Ibci = Ir / Br;
00446   gbci = gir / Br;
00447   pnJunctionBIP (Ubc, Isc, Ut * Nc, Ibcn, gbcn);
00448   Ibcn += gtiny * Ubc;
00449   gbcn += gtiny;
00450   Ibc = Ibci + Ibcn;
00451   gbc = gbci + gbcn;
00452 #endif
00453 
00454   // compute base charge quantities
00455   Q1 = 1 / (1 - Ubc * Vaf - Ube * Var);
00456   Q2 = If * Ikf + Ir * Ikr;
00457   nr_double_t SArg = 1.0 + 4.0 * Q2;
00458   nr_double_t Sqrt = SArg > 0 ? qucs::sqrt (SArg) : 1;
00459   Qb = Q1 * (1 + Sqrt) / 2;
00460   dQbdUbe = Q1 * (Qb * Var + gif * Ikf / Sqrt);
00461   dQbdUbc = Q1 * (Qb * Vaf + gir * Ikr / Sqrt);
00462 
00463   // during transient analysis only
00464   if (doTR) {
00465     // calculate excess phase influence
00466     If /= Qb;
00467     excessPhase (cexState, If, gif);
00468     If *= Qb;
00469   }
00470 
00471   // compute transfer current
00472   It = (If - Ir) / Qb;
00473 
00474   // compute forward and backward transconductance
00475   gitf = (+gif - It * dQbdUbe) / Qb;
00476   gitr = (-gir - It * dQbdUbc) / Qb;
00477 
00478   // compute old SPICE values
00479   go = -gitr;
00480   gm = +gitf - go;
00481   setOperatingPoint ("gm", gm);
00482   setOperatingPoint ("go", go);
00483 
00484   // calculate current-dependent base resistance
00485   if (Rbm != 0.0) {
00486     if (Irb != 0.0) {
00487       nr_double_t a, b, z;
00488       a = (Ibci + Ibcn + Ibei + Iben) / Irb;
00489       a = std::max (a, NR_TINY); // enforce positive values
00490       z = (qucs::sqrt (1 + 144 / sqr (pi) * a) - 1) / 24 * sqr (pi) / qucs::sqrt (a);
00491       b = qucs::tan (z);
00492       Rbb = Rbm + 3 * (Rb - Rbm) * (b - z) / z / sqr (b);
00493     }
00494     else {
00495       Rbb = Rbm + (Rb - Rbm) / Qb;
00496     }
00497     rb->setScaledProperty ("R", Rbb);
00498     rb->calcDC ();
00499   }
00500 
00501   // compute autonomic current sources
00502   IeqB = Ibe - Ube * gbe;
00503   IeqC = Ibc - Ubc * gbc;
00504 #if NEWSGP
00505   IeqE = It - Ube * gitf - Ubc * gitr;
00506 #else
00507   IeqE = It - Ube * gm - Uce * go;
00508 #endif
00509   IeqS = 0;
00510   setI (NODE_B, (-IeqB - IeqC) * pol);
00511   setI (NODE_C, (+IeqC - IeqE - IeqS) * pol);
00512   setI (NODE_E, (+IeqB + IeqE) * pol);
00513   setI (NODE_S, (+IeqS) * pol);
00514 
00515   // apply admittance matrix elements
00516 #if NEWSGP
00517   setY (NODE_B, NODE_B, gbc + gbe);
00518   setY (NODE_B, NODE_C, -gbc);
00519   setY (NODE_B, NODE_E, -gbe);
00520   setY (NODE_B, NODE_S, 0);
00521   setY (NODE_C, NODE_B, -gbc + gitf + gitr);
00522   setY (NODE_C, NODE_C, gbc - gitr);
00523   setY (NODE_C, NODE_E, -gitf);
00524   setY (NODE_C, NODE_S, 0);
00525   setY (NODE_E, NODE_B, -gbe - gitf - gitr);
00526   setY (NODE_E, NODE_C, gitr);
00527   setY (NODE_E, NODE_E, gbe + gitf);
00528   setY (NODE_E, NODE_S, 0);
00529   setY (NODE_S, NODE_B, 0);
00530   setY (NODE_S, NODE_C, 0);
00531   setY (NODE_S, NODE_E, 0);
00532   setY (NODE_S, NODE_S, 0);
00533 #else
00534   setY (NODE_B, NODE_B, gbc + gbe);
00535   setY (NODE_B, NODE_C, -gbc);
00536   setY (NODE_B, NODE_E, -gbe);
00537   setY (NODE_B, NODE_S, 0);
00538   setY (NODE_C, NODE_B, -gbc + gm);
00539   setY (NODE_C, NODE_C, go + gbc);
00540   setY (NODE_C, NODE_E, -go - gm);
00541   setY (NODE_C, NODE_S, 0);
00542   setY (NODE_E, NODE_B, -gbe - gm);
00543   setY (NODE_E, NODE_C, -go);
00544   setY (NODE_E, NODE_E, gbe + go + gm);
00545   setY (NODE_E, NODE_S, 0);
00546   setY (NODE_S, NODE_B, 0);
00547   setY (NODE_S, NODE_C, 0);
00548   setY (NODE_S, NODE_E, 0);
00549   setY (NODE_S, NODE_S, 0);
00550 #endif
00551 }
00552 
00553 void bjt::saveOperatingPoints (void) {
00554   nr_double_t Vbe, Vbc;
00555   Vbe = real (getV (NODE_B) - getV (NODE_E)) * pol;
00556   Vbc = real (getV (NODE_B) - getV (NODE_C)) * pol;
00557   Ucs = real (getV (NODE_S) - getV (NODE_C)) * pol;
00558   setOperatingPoint ("Vbe", Vbe);
00559   setOperatingPoint ("Vbc", Vbc);
00560   setOperatingPoint ("Vce", Vbe - Vbc);
00561   setOperatingPoint ("Vcs", Ucs);
00562   if (deviceEnabled (cbcx)) {
00563     Ubx = real (cbcx->getV (NODE_1) - cbcx->getV (NODE_2)) * pol;
00564     setOperatingPoint ("Vbx", Ubx);
00565   }
00566 }
00567 
00568 void bjt::loadOperatingPoints (void) {
00569   Ube = getOperatingPoint ("Vbe");
00570   Ubc = getOperatingPoint ("Vbc");
00571   Uce = getOperatingPoint ("Vce");
00572   Ucs = getOperatingPoint ("Vcs");
00573 }
00574 
00575 void bjt::calcOperatingPoints (void) {
00576 
00577   // fetch device model parameters
00578   nr_double_t Cje0 = getScaledProperty ("Cje");
00579   nr_double_t Vje  = getScaledProperty ("Vje");
00580   nr_double_t Mje  = getPropertyDouble ("Mje");
00581   nr_double_t Cjc0 = getScaledProperty ("Cjc");
00582   nr_double_t Vjc  = getScaledProperty ("Vjc");
00583   nr_double_t Mjc  = getPropertyDouble ("Mjc");
00584   nr_double_t Xcjc = getPropertyDouble ("Xcjc");
00585   nr_double_t Cjs0 = getScaledProperty ("Cjs");
00586   nr_double_t Vjs  = getScaledProperty ("Vjs");
00587   nr_double_t Mjs  = getPropertyDouble ("Mjs");
00588   nr_double_t Fc   = getPropertyDouble ("Fc");
00589   nr_double_t Vtf  = getPropertyDouble ("Vtf");
00590   nr_double_t Tf   = getPropertyDouble ("Tf");
00591   nr_double_t Xtf  = getPropertyDouble ("Xtf");
00592   nr_double_t Itf  = getScaledProperty ("Itf");
00593   nr_double_t Tr   = getPropertyDouble ("Tr");
00594 
00595   nr_double_t Cbe, Cbci, Cbcx, Ccs;
00596 
00597   // interpret zero as infinity for that model parameter
00598   Vtf = Vtf > 0 ? 1.0 / Vtf : 0;
00599 
00600   // depletion capacitance of base-emitter diode
00601   Cbe = pnCapacitance (Ube, Cje0, Vje, Mje, Fc);
00602   Qbe = pnCharge (Ube, Cje0, Vje, Mje, Fc);
00603 
00604   // diffusion capacitance of base-emitter diode
00605   if (If != 0.0) {
00606     nr_double_t e, Tff, dTffdUbe, dTffdUbc, a;
00607     a = 1 / (1 + Itf / If);
00608     e = 2 * qucs::exp (std::min (Ubc * Vtf, 709.0));
00609     Tff = Tf * (1 + Xtf * sqr (a) * e);
00610     dTffdUbe = Tf * Xtf * 2 * gif * Itf * cubic (a) / sqr (If) * e;
00611     Cbe += (If * dTffdUbe + Tff * (gif - If / Qb * dQbdUbe)) / Qb;
00612     Qbe += If * Tff / Qb;
00613     dTffdUbc = Tf * Xtf * Vtf * sqr (a) * e;
00614     dQbedUbc = If / Qb * (dTffdUbc - Tff / Qb * dQbdUbc);
00615   }
00616 
00617   // depletion and diffusion capacitance of base-collector diode
00618   Cbci = pnCapacitance (Ubc, Cjc0 * Xcjc, Vjc, Mjc, Fc) + Tr * gir;
00619   Qbci = pnCharge (Ubc, Cjc0 * Xcjc, Vjc, Mjc, Fc) + Tr * Ir;
00620 
00621   // depletion and diffusion capacitance of external base-collector capacitor
00622   Cbcx = pnCapacitance (Ubx, Cjc0 * (1 - Xcjc), Vjc, Mjc, Fc);
00623   Qbcx = pnCharge (Ubx, Cjc0 * (1 - Xcjc), Vjc, Mjc, Fc);
00624 
00625   // depletion capacitance of collector-substrate diode
00626   Ccs = pnCapacitance (Ucs, Cjs0, Vjs, Mjs);
00627   Qcs = pnCharge (Ucs, Cjs0, Vjs, Mjs);
00628 
00629   // finally save the operating points
00630   setOperatingPoint ("Cbe", Cbe);
00631   setOperatingPoint ("Cbci", Cbci);
00632   setOperatingPoint ("Cbcx", Cbcx);
00633   setOperatingPoint ("Ccs", Ccs);
00634   setOperatingPoint ("gmf", gitf);
00635   setOperatingPoint ("gmr", gitr);
00636   setOperatingPoint ("gmu", gbci + gbcn);
00637   setOperatingPoint ("gpi", gbei + gben);
00638   setOperatingPoint ("Rbb", Rbb);
00639   setOperatingPoint ("Ibe", Ibe);
00640   setOperatingPoint ("Ice", It);
00641 }
00642 
00643 void bjt::initSP (void) {
00644   allocMatrixS ();
00645   processCbcx ();
00646   if (deviceEnabled (cbcx)) {
00647     cbcx->initSP ();
00648     cbcx->initNoiseSP ();
00649   }
00650 }
00651 
00652 void bjt::processCbcx (void) {
00653   nr_double_t Xcjc = getPropertyDouble ("Xcjc");
00654   nr_double_t Rbm  = getScaledProperty ("Rbm");
00655   nr_double_t Cjc0 = getScaledProperty ("Cjc");
00656 
00657   /* if necessary then insert external capacitance between internal
00658      collector node and external base node */
00659   if (Rbm != 0.0 && Cjc0 != 0.0 && Xcjc != 1.0) {
00660     if (!deviceEnabled (cbcx)) {
00661       cbcx = splitCapacitor (this, cbcx, "Cbcx", rb->getNode (NODE_1),
00662                              getNode (NODE_C));
00663     }
00664     cbcx->setProperty ("C", getOperatingPoint ("Cbcx"));
00665   }
00666   else {
00667     disableCapacitor (this, cbcx);
00668   }
00669 }
00670 
00671 void bjt::initAC (void) {
00672   allocMatrixMNA ();
00673   processCbcx ();
00674   if (deviceEnabled (cbcx)) {
00675     cbcx->initAC ();
00676     cbcx->initNoiseAC ();
00677   }
00678 }
00679 
00680 void bjt::calcAC (nr_double_t frequency) {
00681   setMatrixY (calcMatrixY (frequency));
00682 }
00683 
00684 void bjt::calcNoiseAC (nr_double_t frequency) {
00685   setMatrixN (calcMatrixCy (frequency));
00686 }
00687 
00688 #define qbeState 0 // base-emitter charge state
00689 #define cbeState 1 // base-emitter current state
00690 #define qbcState 2 // base-collector charge state
00691 #define cbcState 3 // base-collector current state
00692 #define qcsState 4 // collector-substrate charge state
00693 #define ccsState 5 // collector-substrate current state
00694 
00695 #define qbxState 0 // external base-collector charge state
00696 #define cbxState 1 // external base-collector current state
00697 
00698 void bjt::initTR (void) {
00699   setStates (7);
00700   initDC ();
00701   doTR = true;
00702 
00703   // handle external base-collector capacitance appropriately
00704   processCbcx ();
00705   if (deviceEnabled (cbcx)) {
00706     cbcx->initTR ();
00707     cbcx->setProperty ("Controlled", getName ());
00708   }
00709 }
00710 
00711 void bjt::calcTR (nr_double_t t) {
00712   calcDC ();
00713   saveOperatingPoints ();
00714   loadOperatingPoints ();
00715   calcOperatingPoints ();
00716 
00717   nr_double_t Cbe  = getOperatingPoint ("Cbe");
00718   nr_double_t Ccs  = getOperatingPoint ("Ccs");
00719   nr_double_t Cbci = getOperatingPoint ("Cbci");
00720   nr_double_t Cbcx = getOperatingPoint ("Cbcx");
00721 
00722   // handle Rbb and Cbcx appropriately
00723   if (Rbb != 0.0) {
00724     rb->setScaledProperty ("R", Rbb);
00725     rb->calcTR (t);
00726     if (deviceEnabled (cbcx)) {
00727       cbcx->clearI ();
00728       cbcx->clearY ();
00729       cbcx->transientCapacitance (qbxState, NODE_1, NODE_2, Cbcx, Ubx, Qbcx);
00730     }
00731   }
00732 
00733   // usual capacitances
00734   transientCapacitance (qbeState, NODE_B, NODE_E, Cbe, Ube, Qbe);
00735   transientCapacitance (qbcState, NODE_B, NODE_C, Cbci, Ubc, Qbci);
00736   transientCapacitance (qcsState, NODE_S, NODE_C, Ccs, Ucs, Qcs);
00737 
00738   // trans-capacitances
00739   transientCapacitanceC (NODE_B, NODE_E, NODE_B, NODE_C, dQbedUbc, Ubc);
00740 }
00741 
00742 void bjt::excessPhase (int istate, nr_double_t& i, nr_double_t& g) {
00743 
00744   // fetch device properties
00745   nr_double_t Ptf = getPropertyDouble ("Ptf");
00746   nr_double_t Tf = getPropertyDouble ("Tf");
00747   nr_double_t td = deg2rad (Ptf) * Tf;
00748 
00749   // return if nothing todo
00750   if (td == 0.0) return;
00751 
00752   // fill-in current history during initialization
00753   if (getMode () & MODE_INIT) fillState (istate, i);
00754 
00755   // calculate current coefficients C1, C2 and C3
00756   nr_double_t * delta = getDelta ();
00757   nr_double_t c3, c2, c1, dn, ra;
00758   c1 = delta[0] / td;
00759   c2 = 3 * c1;
00760   c1 = c2 * c1;
00761   dn = 1 + c1 + c2;
00762   c1 = c1 / dn;
00763   ra = delta[0] / delta[1];
00764   c2 = (1 + ra + c2) / dn;
00765   c3 = ra / dn;
00766 
00767   // update and save current, update transconductance
00768   i = i * c1 + getState (istate, 1) * c2 - getState (istate, 2) * c3;
00769   setState (istate, i);
00770   g = g * c1;
00771 }
00772 
00773 // properties
00774 PROP_REQ [] = {
00775   { "Is", PROP_REAL, { 1e-16, PROP_NO_STR }, PROP_POS_RANGE },
00776   { "Nf", PROP_REAL, { 1, PROP_NO_STR }, PROP_RNGII (0.1, 100) },
00777   { "Nr", PROP_REAL, { 1, PROP_NO_STR }, PROP_RNGII (0.1, 100) },
00778   { "Ikf", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00779   { "Ikr", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00780   { "Vaf", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00781   { "Var", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00782   { "Ise", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00783   { "Ne", PROP_REAL, { 1.5, PROP_NO_STR }, PROP_RNGII (0.1, 100) },
00784   { "Isc", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00785   { "Nc", PROP_REAL, { 2, PROP_NO_STR }, PROP_RNGII (0.1, 100) },
00786   { "Bf", PROP_REAL, { 100, PROP_NO_STR }, PROP_POS_RANGEX },
00787   { "Br", PROP_REAL, { 1, PROP_NO_STR }, PROP_POS_RANGEX },
00788   { "Rbm", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00789   { "Irb", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00790   { "Cje", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00791   { "Vje", PROP_REAL, { 0.75, PROP_NO_STR }, PROP_RNGXI (0, 10) },
00792   { "Mje", PROP_REAL, { 0.33, PROP_NO_STR }, PROP_RNGII (0, 1) },
00793   { "Cjc", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00794   { "Vjc", PROP_REAL, { 0.75, PROP_NO_STR }, PROP_RNGXI (0, 10) },
00795   { "Mjc", PROP_REAL, { 0.33, PROP_NO_STR }, PROP_RNGII (0, 1) },
00796   { "Xcjc", PROP_REAL, { 1, PROP_NO_STR }, PROP_RNGII (0, 1) },
00797   { "Cjs", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00798   { "Vjs", PROP_REAL, { 0.75, PROP_NO_STR }, PROP_RNGXI (0, 10) },
00799   { "Mjs", PROP_REAL, { 0, PROP_NO_STR }, PROP_RNGII (0, 1) },
00800   { "Fc", PROP_REAL, { 0.5, PROP_NO_STR }, PROP_RNGII (0, 1) },
00801   { "Vtf", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
00802   { "Tf", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00803   { "Xtf", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00804   { "Itf", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00805   { "Tr", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00806   PROP_NO_PROP };
00807 PROP_OPT [] = {
00808   { "Rc", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00809   { "Re", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00810   { "Rb", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00811   { "Kf", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00812   { "Af", PROP_REAL, { 1, PROP_NO_STR }, PROP_POS_RANGE },
00813   { "Ffe", PROP_REAL, { 1, PROP_NO_STR }, PROP_POS_RANGE },
00814   { "Kb", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00815   { "Ab", PROP_REAL, { 1, PROP_NO_STR }, PROP_POS_RANGE },
00816   { "Fb", PROP_REAL, { 1, PROP_NO_STR }, PROP_POS_RANGE },
00817   { "Temp", PROP_REAL, { 26.85, PROP_NO_STR }, PROP_MIN_VAL (K) },
00818   { "Type", PROP_STR, { PROP_NO_VAL, "npn" }, PROP_RNG_BJT },
00819   { "Ptf", PROP_REAL, { 0, PROP_NO_STR }, PROP_RNGII (-180, +180) },
00820   { "Xtb", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
00821   { "Xti", PROP_REAL, { 3, PROP_NO_STR }, PROP_POS_RANGE },
00822   { "Eg", PROP_REAL, { EgSi, PROP_NO_STR }, PROP_POS_RANGE },
00823   { "Tnom", PROP_REAL, { 26.85, PROP_NO_STR }, PROP_MIN_VAL (K) },
00824   { "Area", PROP_REAL, { 1, PROP_NO_STR }, PROP_POS_RANGEX },
00825   PROP_NO_PROP };
00826 struct define_t bjt::cirdef =
00827   { "BJT", 4, PROP_COMPONENT, PROP_NO_SUBSTRATE, PROP_NONLINEAR, PROP_DEF };