1/* $Id: d_bjt.model,v 26.136 2009/12/08 02:03:49 al Exp $ -*- C++ -*-
2 * Copyright (C) 2002 Albert Davis
3 * Author: Albert Davis <aldavis@gnu.org>
4 *
5 * This file is part of "Gnucap", the Gnu Circuit Analysis Package
6 *
7 * This program is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 3, or (at your option)
10 * any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program; if not, write to the Free Software
19 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
20 * 02110-1301, USA.
21 *------------------------------------------------------------------
22 * Berkeley BJT model
23 * Derived from Spice code, both from 2g6 and 3f4
24 * Recoded for Gnucap model compiler, Al Davis, 2002
25 *------------------------------------------------------------------
26 * data structures and defaults for bjt model.
27 *
28 * netlist syntax:
29 * device:  qxxxx c b e s mname <device args>
30 * model:   .model mname NPN <args>
31 *	or  .model mname PNP <args>
32 *
33 * known BUGs
34 * 1. excess phase partially implemented, disabled, as if PTF == 0.
35 */
36/*--------------------------------------------------------------------------*/
37h_headers {
38#include "d_diode.h"
39}
40cc_headers {
41#include "u_limit.h"
42}
43/*--------------------------------------------------------------------------*/
44device BUILT_IN_BJT {
45  parse_name bjt;
46  model_type BUILT_IN_BJT;
47  id_letter Q;
48  circuit {
49    sync;
50    ports {c b e} {s};
51    local_nodes {
52      ic  short_to=c   short_if="!OPT::rstray || m->rc == 0.";
53      ib short_to=b
54	short_if="!OPT::rstray || (m->rb == 0. && m->rbm == 0.)";
55      ie short_to=e  short_if="!OPT::rstray || m->re == 0.";
56    }
57    // basics
58    cpoly_g Ice {ic,ie  ib,ie} state=ccexxx; // cce, go, gm
59    cpoly_g Ipi {ib ie} state=cpixxx; // cpi, gpi
60    cpoly_g Imu {ib ic}  state=cmuxxx; // cmu, gmu
61    // junction charge effects
62    fpoly_cap Cbx {b ic} state=qbx; // qbx, cqbx
63    fpoly_cap Cbc {ib ic} state=qbc; // qbc, cqbc
64    fpoly_cap Ccs {s ic} state=qcs // qcs, cqcs
65	omit="!(_n[n_s].n_())";
66    fpoly_cap Cbe {ib ie ib ic} state=qbe; // qbe, cqbe, geqbc
67    // parasitics
68    resistor Rc {c ic} value="m->rc / c->area"
69	omit="!OPT::rstray || m->rc == 0."; // gcpr
70    resistor Re {e ie} value="m->re / c->area"
71	omit="!OPT::rstray || m->re == 0."; // gepr
72    cpoly_g Yb {b ib} state=ixxxx // gx.  should be "eval"
73	omit="!OPT::rstray || (m->rb == 0. && m->rbm == 0.)";
74    capacitor Cbcp {b c} value="m->cbcp * c->area"
75	omit="!OPT::cstray || m->cbcp == 0.";
76    capacitor Cbep {b e} value="m->cbep * c->area"
77	omit="!OPT::cstray || m->cbep == 0.";
78    capacitor Cbs {b s} value="(m->cbsp + m->ccsp) * c->area"
79	omit="!OPT::cstray || ((m->cbsp + m->ccsp) == 0.)";
80  }
81  tr_probe {
82    v = "@n_c[V] - @n_e[V]";
83    "vbei{nt}" = "vbe";
84    "vbci{nt}" = "vbc";
85    "vbxi{nt}" = "vbx";
86    "vcsi{nt}" = "vcs";
87    vbs = "@n_b[V] - @n_s[V]";
88    vbe = "@n_b[V] - @n_e[V]";
89    vbc = "@n_b[V] - @n_c[V]";
90    vbx = "@n_b[V] - @n_ib[V]";
91    vcs = "@n_c[V] - @n_s[V]";
92    vcb = "@n_c[V] - @n_b[V]";
93    vce = "@n_c[V] - @n_e[V]";
94    ves = "@n_e[V] - @n_s[V]";
95    veb = "@n_e[V] - @n_b[V]";
96    vec = "@n_e[V] - @n_c[V]";
97    vb = "@n_b[V]";
98    vc = "@n_c[V]";
99    ve = "@n_e[V]";
100    vs = "@n_s[V]";
101    vbi = "@n_ib[V]";
102    vci = "@n_ic[V]";
103    vei = "@n_ie[V]";
104
105    "i" = "cce";
106    "ice" = "cce";
107    "iceo{ffset}" = "ccexxx";
108    hoe = "go";
109    "ro{e}" = "(go != 0.) ? 1/go : BIGBIG";
110
111    ipi = "cpi";
112    "ipio{ffset}" = "cpixxx";
113    rpi = "(gpi != 0.) ? 1/gpi : BIGBIG";
114    hie = "(gpi != 0.) ? 1/gpi : BIGBIG";
115
116    imu = "cmu";
117    "imuo{ffset}" = "cmuxxx";
118    rmu = "(gmu != 0.) ? 1/gmu : BIGBIG";
119
120    ib = "cpi + cmu";
121    rx = "(gx != NA) ? 1/gx : 0.";
122
123    ic = "cce - cmu";
124    ie = "-cce -cpi";
125
126    cbx = "cqbx";
127    cbc = "cqbc";
128    cmu = "cqbc";
129    ccs = "cqcs";
130    cbe = "cqbe";
131    cpi = "cqbe";
132
133    p = "@Ice[P] + @Ipi[P] + @Imu[P] + @Rc[P] + @Re[P] + @Yb[P]
134       + @Cbx[P] + @Cbc[P] + @Ccs[P] + @Cbe[P]";
135    pd = "@Ice[PD] + @Ipi[PD] + @Imu[PD] + @Rc[PD] + @Re[PD] + @Yb[PD]
136       + @Cbx[PD] + @Cbc[PD] + @Ccs[PD] + @Cbe[PD]";
137    ps = "@Ice[PS] + @Ipi[PS] + @Imu[PS] + @Rc[PS] + @Re[PS] + @Yb[PS]
138       + @Cbx[PS] + @Cbc[PS] + @Ccs[PS] + @Cbe[PS]";
139
140    status = "static_cast<double>(converged() * 2)";
141  }
142  device {
143    calculated_parameters {
144      double vbe "B-E voltage"; // gather
145      double vbc "B-C voltage";
146      double vbx "B-C voltage (ext base)";
147      double vcs "C-S voltage";
148
149      double cce "collector-emitter current";
150      double ccexxx;
151      double go "Small signal output conductance";
152      double gm "Small signal transconductance";
153
154      double cpi "emitter-base current";
155      double cpixxx;
156      double gpi "Small signal input conductance - pi";
157
158      double cmu "collector-base current";
159      double cmuxxx;
160      double gmu "Small signal conductance - mu";
161
162      double ixxxx "Current offset at base node, constant" default=0.;
163      double gx "dix/dvbb Conductance from base to internal base";
164
165      double qbx "Charge storage B-X junction";
166      double cqbx "Cap. due to charge storage in B-X jct."; // cbx, capbx
167
168      double qbc "Charge storage B-C junction";
169      double cqbc "Cap. due to charge storage in B-C jct."; // cmu, capbc
170
171      double qcs "Charge storage C-S junction";
172      double cqcs "Cap. due to charge storage in C-S jct."; // ccs, capcs
173
174      double qbe "Charge storage B-E junction";
175      double cqbe "Cap. due to charge storage in B-E jct."; // cpi, capbe
176      double geqcb "d(Ieb)/d(Vcb)";
177
178      double cexbc_0 "Total Capacitance in B-X junction"; // ??
179      double cexbc_1 "";
180      double cexbc_2 "";
181      double _dt_0 "time step";
182      double _dt_1 "";
183    }
184  }
185  common {
186    unnamed area;
187    raw_parameters {
188      double area "area factor" name=Area positive default=1.;
189      bool off "device initially off" name=OFF default=false print_test=off;
190      double icvbe "Initial B-E voltage" name="ICVBE" default=NA;
191      double icvce "Initial C-E voltage" name="ICVCE" default=NA;
192      double temp_c "instance temperature" name="TEMP" default=NA;
193    }
194    calculated_parameters {
195      double oik  "" calculate="m->invrollofff / c->area";
196      double oikr "" calculate="m->invrolloffr / c->area";
197    }
198  }
199  tr_eval {
200    int foo=3;
201  }
202}
203/*--------------------------------------------------------------------------*/
204model BUILT_IN_BJT {
205  level 1;
206  dev_type BUILT_IN_BJT;
207  hide_base;
208  inherit BUILT_IN_DIODE;
209  public_keys {
210    npn polarity=pN;
211    pnp polarity=pP;
212    npn1 polarity=pN;
213    pnp1 polarity=pP;
214  }
215  independent {
216    override {
217      double kf "Flicker Noise Coefficient" name=KF default=0.;
218      double af "Flicker Noise Exponent" name=AF default=1.;
219    }
220    raw_parameters {
221      int level "dummy" name=LEVEL default=1 print_test=false;
222      // basic
223      double bf "Ideal forward beta" name=BF alt_name=BFM default=100.;
224      double br "Ideal reverse beta" name=BR alt_name=BRM default=1.;
225
226      double ibc "bc Saturation Current"
227	name=IBC default=NA print_test="ibe != ibc"
228	final_default="((has_hard_value(i_s)) ? i_s : 1e-16)";
229      double ibe "be Saturation Current"
230	name=IBE default=NA print_test="ibe != ibc"
231	final_default="((has_hard_value(i_s)) ? i_s : 1e-16)";
232      double i_s "Saturation Current"
233	name=IS default=NA print_test="ibe == ibc"
234	final_default="((ibe == ibc) ? ibe : NA)";
235
236      //double iss "bulk to collector or base sat current (ignored)"
237      //name=ISS default=0;
238
239      //double expli "current explosion factor (ignored)"
240      //name=EXPLI default=1e15;
241
242      double nf "Forward emission coefficient" name=NF default=1.;
243      double nr "Reverse emission coefficient" name=NR default=1.;
244      //double ns "substrate leakage emission coefficient (ignored)"
245      //name=NS default=1;
246
247      //int subs "geometry, +1=vertical, -1=lateral (ignored)", name=SUBS
248      //default="(polarity==pP) ? -1 : 1";
249
250      // base width modulation
251      double vaf "Forward Early voltage" name=VAf alt_name=VA\0VBF default=NA;
252      double var "Reverse Early voltage" name=VAR alt_name=VB default=NA;
253
254      // low current beta degeneration
255      double isc "B-C leakage saturation current"
256	name=ISC default=NA final_default="(c4*ibc)";
257      double c4 "obsolete, don't use" name=C4 alt_name=JLC default=0.;
258      double nc "B-C leakage emission coefficient" name=NC default=2.;
259
260      double ise "B-E leakage saturation current"
261	name=ISE default=NA final_default="(c2*ibe)";
262      double c2 "obsolete, don't use" name=C2 alt_name=JLE default=0.;
263      double ne "B-E leakage emission coefficient" name=NE default=1.5;
264
265      // high current beta degeneration
266      double ikf "Forward beta roll-off corner current"
267	name=IKf alt_name=IK\0JBF default=NA;
268      double ikr "reverse beta roll-off corner current"
269	name=IKR alt_name=JBR default=NA;
270
271      // parasitic resistance
272      double irb "Current for base resistance=(rb+rbm)/2"
273	name=IRB alt_name=JRB\0IOB default=NA;
274      double rb "Zero bias base resistance" name=RB default=0.;
275      double rbm "Minimum base resistance"
276	name=RBM default=NA final_default=rb;
277      double re "Emitter resistance" name=RE default=0.;
278      double rc "Collector resistance" name=RC default=0.;
279
280      // parasitic capacitance
281      double cbcp "external BC capacitance"
282	name=CBCP print_test="cbcp!=0." default=0.;
283      double cbep "external BE capacitance"
284	name=CBEP print_test="cbep!=0." default=0.;
285      double cbsp "external BS capacitance (lateral)"
286	name=CBSP print_test="cbsp!=0." default=0.;
287      double ccsp "external BS capacitance (vertical)"
288	name=CCSP print_test="ccsp!=0." default=0.;
289
290      // junction capacitance
291      double cjc "Zero bias B-C depletion capacitance" name=CJC default=0.;
292      double cje "Zero bias B-E depletion capacitance" name=CJE default=0.;
293      double cjs "Zero bias C-S capacitance"
294	name=CJS alt_name=CCS\0CSUB default=0.;
295
296      double fc "Forward bias junction fit parameter"
297	name=FC default=NA final_default=.5 quiet_max=.9999;
298
299      double mjc "B-C junction grading coefficient"
300	name=MJC alt_name=MC\0MJ default=.33;
301      double mje "B-E junction grading coefficient"
302	name=MJE alt_name=ME default=.33;
303      double mjs "Substrate junction grading coefficient"
304	name=MJS alt_name=MSub\0MS\0ESUB default=0.; // alt name ESUB??
305
306      double vjc "B-C built in potential" name=VJC alt_name=PC default=.75;
307      double vje "B-E built in potential" name=VJE alt_name=PE default=.75;
308      double vjs "Substrate junction built in potential"
309	name=VJS alt_name=PSub\0PS default=.75;
310
311      double xcjc "Fraction of B-C cap to internal base"
312	name=XCJC alt_name=CDIS default=1.;
313
314      // transit time
315      double itf "High current dependence of TF"
316	name=ITF alt_name=JTF default=0.;
317      double ptf "Excess phase" name=PTf default=0.;
318
319      double tf "Ideal forward transit time" name=TF default=0.;
320      double tr "Ideal reverse transit time" name=TR default=0.;
321
322      double vtf "Voltage giving VBC dependence of TF" name=VTF default=NA;
323      double xtf "Coefficient for bias dependence of TF" name=XTF default=0.;
324
325      // temperature effects
326      double xtb "Forward and reverse beta temp. exp."
327	name=XTB alt_name=TB default=0.; // alt name TCB
328      double xti "Temp. exponent for IS" name=XTI default=3.;
329      double eg "Energy gap for IS temp. dependency" name=EG default=1.11;
330    }
331    calculated_parameters {
332      double tnom_k "nominal temperature, kelvin"
333	calculate="_tnom_c + P_CELSIUS0";
334      polarity_t polarity "" default=pN;
335      double invearlyvoltf "Inverse early voltage:forward" name=INVEARLYVOLTF
336	calculate="(has_nz_value(vaf)) ? 1./vaf : 0.";
337      double invearlyvoltr "Inverse early voltage:reverse" name=INVEARLYVOLTR
338	calculate="(has_nz_value(var)) ? 1./var : 0.";
339      double invrollofff "Inverse roll off - forward" name=INVROLLOFFF
340	calculate="(has_nz_value(ikf)) ? 1./ikf : 0.";
341      double invrolloffr "Inverse roll off - reverse" name=INVROLLOFFR
342	calculate="(has_nz_value(ikr)) ? 1./ikr : 0.";
343      double transtimevbcfact "Transit time VBC factor" name=TRANSTIMEVBCFACT
344	calculate="(has_nz_value(vtf)) ? 1./(vtf*1.44) : 0.";
345      double excessphasefactor "Excess phase fact." name=EXCESSPHASEFACTOR
346        calculate="(ptf * DTOR) * tf";
347      double xfc "" calculate="log(1 - fc)";
348      double f2 "" calculate="exp((1 + mje) * xfc)";
349      double f3 "" calculate="1 - fc * (1 + mje)";
350      double f6 "" calculate="exp((1 + mjc) * xfc)";
351      double f7 "" calculate="1 - fc * (1 + mjc)";
352    }
353  }
354  temperature_dependent {
355    calculated_parameters {
356      double vt "";
357      //double is "BJTtSatCur" calculate = "m->is * factor";
358      double ibc "BJTtSatCur" calculate = "m->ibc * factor";
359      double ibe "BJTtSatCur" calculate = "m->ibe * factor";
360      double BetaF "BJTtBetaF" calculate = "m->bf * bfactor";
361      double BetaR "BJTtBetaR" calculate = "m->br * bfactor";
362      double BEleakCur "BJTtBEleakCur"
363	calculate = "m->ise * exp(factlog/m->ne) / bfactor";
364      double BCleakCur "BJTtBCleakCur"
365	calculate = "m->isc * exp(factlog/m->nc) / bfactor";
366      double BEpot "BJTtBEpot";
367      double BEcap "BJTtBEcap";
368      double DepCap "";
369      double f1 "";
370      double BCpot "BJTtBCpot";
371      double BCcap "BJTtBCcap";
372      double f4 "";
373      double f5 "";
374      double Vcrit "" calculate="vt * log(vt / (M_SQRT2 * m->ibe))";
375    }
376    code_pre {
377      const double reftemp = 300.15;
378      double temp_k =
379	((has_hard_value(c->temp_c)) ? c->temp_c : d->_sim->_temp_c) + P_CELSIUS0;
380      double fact1 = m->tnom_k / reftemp;
381      double fact2 = temp_k / reftemp;
382      double tempratio  = temp_k / m->tnom_k; // fact2/fact1
383      double kt = temp_k * P_K;
384      vt = temp_k * P_K_Q;
385      double egap = 1.16 - (7.02e-4*temp_k*temp_k) / (temp_k+1108.); // egfet
386      // double arg = (m->eg*tempratio - egap) / (2*kt);
387      double arg = -egap/(2 * kt) + 1.1150877 / (P_K * (reftemp+reftemp));
388      double pbfact = -2 * vt * (1.5 * log(fact2) + P_Q * arg);
389      double ratlog = log(tempratio);
390      double ratio1 = tempratio - 1;
391      double factlog = ratio1 * m->eg / vt + m->xti * ratlog;
392      double factor = exp(factlog);
393      double bfactor = exp(ratlog * m->xtb);
394
395    }
396    code_post {
397      {
398	double pbo = (m->vje - pbfact) / fact1;
399	BEpot = fact2 * pbo + pbfact;
400	double gmaold = (m->vje - pbo) / pbo;
401	double gmanew = (BEpot - pbo) / pbo;
402	BEcap = (m->cje / (1 + m->mje * (4e-4*(m->tnom_k-reftemp)-gmaold)))
403	  * (1 + m->mje * (4e-4*(temp_k-reftemp)-gmanew));
404	DepCap = m->fc * BEpot;
405	f1 = BEpot * (1 - exp((1 - m->mje) * m->xfc)) / (1 - m->mje);
406      }
407      {
408	double pbo = (m->vjc - pbfact) / fact1;
409	BCpot = fact2 * pbo + pbfact;
410	double gmaold = (m->vjc - pbo) / pbo;
411	double gmanew = (BCpot - pbo) / pbo;
412	BCcap = (m->cjc / (1 + m->mjc
413			   * (4e-4*(m->tnom_k-reftemp)-gmaold)))
414	  * (1 + m->mjc * (4e-4*(temp_k-reftemp)-gmanew));
415	f4 = m->fc * BCpot;
416	f5 = BCpot * (1 - exp((1 - m->mjc) * m->xfc)) / (1 - m->mjc);
417      }
418    }
419  }
420  tr_eval {
421    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
422    trace0("--------------------------");
423    trace1(d->long_label().c_str(), d->evaliter());
424    trace4("", d->vbe, d->vbc, d->vbx, d->vcs);
425    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
426    double cbe, gbe;
427    { // d->cpi, d->gpi, d->cpixxx
428      // uses: d->vbe
429      double cben, gben;
430      double vtn = t->vt * m->nf;
431      double csat = t->ibe * c->area;
432      double C2 = t->BEleakCur * c->area;
433      if (d->vbe > -5 * vtn) {
434	double evbe = exp(d->vbe / vtn);
435	cbe = csat * (evbe-1) + OPT::gmin * d->vbe;
436	gbe = csat * evbe/vtn + OPT::gmin;
437	if (C2 == 0.) {
438	  cben = 0.;
439	  gben = 0.;
440	}else{
441	  double vte = m->ne * t->vt;
442	  double evben = exp(d->vbe / vte);
443	  cben = C2 * (evben-1);
444	  gben = C2 * evben/vte;
445	}
446	trace4("vbe on", cbe, gbe, cben, gben);
447      }else{
448	gbe = -csat/d->vbe + OPT::gmin;
449	cbe = gbe * d->vbe;
450	gben = -C2 / d->vbe;
451	cben = gben * d->vbe;
452	trace4("vbe off", cbe, gbe, cben, gben);
453      }
454      d->cpi = cbe / t->BetaF + cben;
455      d->gpi = gbe / t->BetaF + gben;
456      d->cpixxx = d->cpi - d->vbe * d->gpi;
457      trace3("", t->BetaF, d->cpi, d->gpi);
458    }
459    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
460    double cbc, gbc, cbcn;
461    { // d->cmu, d->gmu, d->cmuxxx
462      // uses: d->vbc
463      double gbcn;
464      double vtn = t->vt * m->nr;
465      double csat = t->ibc * c->area;
466      double C4 = t->BCleakCur * c->area;
467      if (d->vbc > -5 * vtn) {
468	double evbc = exp(d->vbc / vtn);
469	cbc = csat * (evbc-1) + OPT::gmin * d->vbc;
470	gbc = csat * evbc/vtn + OPT::gmin;
471	if (C4 == 0.) {
472	  cbcn = 0.;
473	  gbcn = 0.;
474	}else{
475	  double vtc = m->nc * t->vt;
476	  double evbcn = exp(d->vbc / vtc);
477	  cbcn = C4 * (evbcn-1);
478	  gbcn = C4 * evbcn/vtc;
479	}
480	trace4("vbc on", cbc, gbc, cbcn, gbcn);
481      }else{
482	gbc = -csat/d->vbc + OPT::gmin;
483	cbc = gbc * d->vbc;
484	gbcn = -C4 / d->vbc;
485	cbcn = gbcn * d->vbc;
486	trace4("vbc off", cbc, gbc, cbcn, gbcn);
487      }
488      d->cmu = cbc / t->BetaR + cbcn;
489      d->gmu = gbc / t->BetaR + gbcn;
490      d->cmuxxx = d->cmu - d->vbc * d->gmu;
491      trace3("", t->BetaR, d->cmu, d->gmu);
492    }
493    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
494    //   determine base charge terms
495    double qb, dqbdve, dqbdvc;
496    {
497      double q1 = 1 / (1 - m->invearlyvoltf*d->vbc - m->invearlyvoltr*d->vbe);
498      if(c->oik == 0. && c->oikr == 0.) {
499	qb = q1;
500	dqbdve = q1 * qb * m->invearlyvoltr;
501	dqbdvc = q1 * qb * m->invearlyvoltf;
502	trace4("!oik", q1, qb, dqbdve, dqbdvc);
503      }else{
504	double q2 = c->oik * cbe + c->oikr * cbc;
505	double arg = std::max(0., 1+4*q2);
506	double sqarg = (arg == 0.) ? 1 : sqrt(arg);
507	qb = q1 * (1+sqarg) / 2;
508	dqbdve = q1 * (qb * m->invearlyvoltr + c->oik  * gbe / sqarg);
509	dqbdvc = q1 * (qb * m->invearlyvoltf + c->oikr * gbc / sqarg);
510	trace2("", c->oik, c->oikr);
511	trace4("oik", q1, qb, dqbdve, dqbdvc);
512      }
513    }
514    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
515    // weil's approx. for excess phase applied with backward-euler integration
516    {
517      double cc = 0.;
518      double cex = cbe;
519      double gex = gbe;
520      if (0 && m->excessphasefactor != 0.) {
521	unreachable();
522	incomplete(); // doesn't save old values of cexbc, so disabled
523	double arg1 = d->_dt_0 / m->excessphasefactor;
524	double arg2 = 3 * arg1;
525	arg1 *= arg2;
526	double denom = 1 + arg1 + arg2;
527	double arg3 = arg1 / denom;
528	if (_sim->is_initial_step()) {
529	  d->cexbc_2 = d->cexbc_1 = cbe / qb;
530	}else{
531	}
532	cc = (d->cexbc_1 * (1 + d->_dt_0/d->_dt_1 + arg2)
533	      - d->cexbc_2 * d->_dt_0/d->_dt_1) / denom;
534	cex *= arg3;
535	gex *= arg3;
536      }
537      d->cexbc_0 = cc + cex / qb;
538
539      d->cce = cc + (cex-cbc)/qb - cbc/t->BetaR - cbcn;
540      d->go = (gbc + (cex-cbc)*dqbdvc / qb) / qb;
541      d->gm = (gex - (cex-cbc)*dqbdve / qb) / qb - d->go;
542      d->ccexxx = d->cce - ((d->vbe - d->vbc) * d->go + d->vbe * d->gm);
543      trace4("", d->cce, d->go, d->gm, d->cce/t->vt);
544    }
545    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
546    // d->gx
547    // should be moved to a private eval
548    // may use d->cpi, d->cmu, qb
549    {
550      if (!OPT::rstray || (!has_nz_value(m->rb) && !has_nz_value(m->rbm))) {
551	trace3("", m->rb, m->irb, d->gx);
552	assert(d->gx == NA);
553      }else{
554	double rx = NA;
555	double rbpr = m->rbm / c->area;
556	double rbpi = m->rb / c->area - rbpr;
557	if (has_nz_value(m->irb)) {itested();//554
558	  // base resistance lowering at high current
559	  double cb = d->cpi + d->cmu;
560	  double xjrb = m->irb * c->area;
561	  double arg1 = std::max(cb/xjrb, 1e-9);
562	  double arg2 = (-1 + sqrt(1+14.59025*arg1)) / 2.4317 / sqrt(arg1);
563	  arg1 = tan(arg2);
564	  rx = rbpr + 3 * rbpi * (arg1-arg2) / arg2 / arg1 / arg1;
565	}else{
566	  rx = rbpr + rbpi / qb;
567	}
568	trace3("", m->rb, m->irb, rx);
569	assert(rx != NA);
570	assert(rx != 0.);
571	d->gx = 1 / rx;
572	trace1("", d->gx);
573      }
574    }
575    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
576    const bool charge_computation_needed = OPT::cstray;
577    if (charge_computation_needed) {
578      if (has_nz_value(m->tf) && d->vbe > 0.) {
579	double argtf = NA;
580	double arg2  = NA;
581	double arg3  = NA;
582	if (has_nz_value(m->xtf)) {itested();//579
583	  if (m->transtimevbcfact != 0.) {
584	    argtf = m->xtf * exp(d->vbc * m->transtimevbcfact);
585	  }else{untested();
586	    argtf = m->xtf;
587	  }
588	  arg2 = argtf;
589	  if(has_nz_value(m->itf)) {
590	    double temp = cbe / (cbe + m->itf * c->area);
591	    argtf *= temp*temp;
592	    arg2  *= (3-temp-temp);
593	  }else{
594	  }
595	  arg3 = cbe * argtf * m->transtimevbcfact;
596	}else{
597	  arg3 = arg2 = argtf = 0.;
598	}
599	assert(argtf != NA);
600	assert(arg2  != NA);
601	assert(arg3  != NA);
602	cbe *= (1+argtf) / qb;
603	gbe = (gbe * (1+arg2) - cbe * dqbdve) / qb;
604	d->geqcb=m->tf*(arg3-cbe*dqbdvc)/qb;
605      }else{
606	d->geqcb=0.;
607      }
608      { // d->qbe, d->cqbe
609	// uses: d->vbe, cbe, gbe
610	double czbe = t->BEcap * c->area;
611	if (d->vbe < t->DepCap) {
612	  double arg = 1 - d->vbe / t->BEpot;
613	  double sarg = pow(arg, -m->mje);
614	  d->qbe = m->tf * cbe + t->BEpot * czbe * (1-arg*sarg) / (1 - m->mje);
615	  d->cqbe = m->tf * gbe + czbe * sarg;
616	}else{
617	  double czbef2 = czbe / m->f2;
618	  d->qbe = m->tf * cbe + czbe * t->f1
619	    + czbef2 * (m->f3 * (d->vbe - t->DepCap)
620			+ (m->mje / (2. * t->BEpot))
621			* (d->vbe * d->vbe - t->DepCap * t->DepCap));
622	  d->cqbe = m->tf * gbe + czbef2 * (m->f3 + m->mje*d->vbe / t->BEpot);
623	}
624      }
625      { // d->qbc, d->cqbc
626	// uses: d->vbc, cbc, gbc
627	double czbc = t->BCcap * c->area * m->xcjc;
628	if (d->vbc < t->f4) {
629	  double arg = 1 - d->vbc / t->BCpot;
630	  double sarg = pow(arg, -m->mjc);
631	  d->qbc = m->tr *cbc + t->BCpot *czbc * (1 - arg*sarg) / (1 - m->mjc);
632	  d->cqbc = m->tr * gbc + czbc * sarg;
633	}else{
634	  double czbcf2 = czbc / m->f6;
635	  d->qbc = m->tr * cbc + czbc * t->f5
636	    + czbcf2 * (m->f7 * (d->vbc-t->f4)
637		+ (m->mjc/(t->BCpot+t->BCpot)) * (d->vbc*d->vbc-t->f4*t->f4));
638	  d->cqbc = m->tr * gbc + czbcf2 * (m->f7 + m->mjc * d->vbc/t->BCpot);
639	}
640      }
641      { // d->qbx, d->cqbx
642	// uses: d->vbx
643	double czbx = t->BCcap * c->area * (1 - m->xcjc);
644	if (d->vbx < t->f4) {
645	  double arg = 1 - d->vbx / t->BCpot;
646	  double sarg = pow(arg, -m->mjc);
647	  d->qbx = t->BCpot * czbx * (1 - arg*sarg) / (1 - m->mjc);
648	  d->cqbx = czbx * sarg;
649	}else{
650	  double czbxf2 = czbx / m->f6;
651	  d->qbx = czbx * t->f5 + czbxf2
652	    * (m->f7 * (d->vbx-t->f4)
653	      + (m->mjc / (t->BCpot+t->BCpot)) * (d->vbx*d->vbx-t->f4*t->f4));
654	  d->cqbx = czbxf2 * (m->f7 + m->mjc * d->vbx / t->BCpot);
655	}
656      }
657      { // d->qcs, d->cqcs
658	// uses: d->vcs
659	double czcs = m->cjs * c->area;
660	if (d->vcs < 0.) {
661	  double arg = 1 - d->vcs / m->vjs;
662	  double sarg = pow(arg, -m->mjs);
663	  d->qcs = m->vjs * czcs * (1 - arg*sarg) / (1 - m->mjs);
664	  d->cqcs = czcs * sarg;
665	}else{
666	  d->qcs = d->vcs * czcs * (1 + m->mjs * d->vcs / (2 * m->vjs));
667	  d->cqcs = czcs * (1 + m->mjs * d->vcs / m->vjs);
668	}
669      }
670    }
671  }
672}
673/*--------------------------------------------------------------------------*/
674cc_direct {
675/*--------------------------------------------------------------------------*/
676bool DEV_BUILT_IN_BJT::tr_needs_eval()const
677{
678  if (is_q_for_eval()) {
679    untested();
680    return false;
681  }else if (!converged()) {
682    return true;
683  }else{
684    const COMMON_BUILT_IN_BJT* c = prechecked_cast<const COMMON_BUILT_IN_BJT*>(common());
685    assert(c);
686    const MODEL_BUILT_IN_BJT* m=prechecked_cast<const MODEL_BUILT_IN_BJT*>(c->model());
687    assert(m);
688    polarity_t polarity = m->polarity;
689    return !(conchk(vbc, polarity*(_n[n_ib].v0()-_n[n_ic].v0()),
690		    OPT::vntol)
691	     && conchk(vbe, polarity*(_n[n_ib].v0()-_n[n_ie].v0()),
692		       OPT::vntol)
693	     && conchk(vcs, polarity*(_n[n_ic].v0()-_n[n_s].v0()),
694		       OPT::vntol));
695  }
696}
697/*--------------------------------------------------------------------------*/
698bool DEV_BUILT_IN_BJT::do_tr()
699{
700  const COMMON_BUILT_IN_BJT* c = prechecked_cast<const COMMON_BUILT_IN_BJT*>(common());
701  assert(c);
702  const MODEL_BUILT_IN_BJT* m = prechecked_cast<const MODEL_BUILT_IN_BJT*>(c->model());
703  assert(m);
704  const TDP_BUILT_IN_BJT T(this);
705  const TDP_BUILT_IN_BJT* t = &T;
706
707  if(_sim->is_initial_step()) {	// initial guess
708    if (c->off) {
709      vbe = 0.;
710    }else{
711      double vt = (_sim->_temp_c + P_CELSIUS0) * P_K_Q;
712      vbe = vt * log(vt / (M_SQRT2 * m->ibe));
713    }
714    vbc = 0.;
715    /* ERROR:  need to initialize VCS, VBX here */
716    vcs = vbx = 0.;
717  }else{				// normal gather
718    vbe = pnj_limit((m->polarity * volts_limited(_n[n_ib], _n[n_ie])),
719		    vbe, t->vt, t->Vcrit);
720    vbc = pnj_limit((m->polarity * volts_limited(_n[n_ib], _n[n_ic])),
721		    vbc, t->vt, t->Vcrit);
722    vbx = m->polarity * volts_limited(_n[n_b], _n[n_ic]);
723    vcs = m->polarity * volts_limited(_n[n_s], _n[n_ic]);
724  }
725
726  if (_sim->uic_now()) {itested();//736
727    if (has_good_value(c->icvbe)) {untested();//737
728      vbe = m->polarity * c->icvbe;
729    }else{itested();//739
730    }
731    if (has_good_value(c->icvce)) {untested();//741
732      vbc = vbe - m->polarity * c->icvce;
733      vbx = vbc;
734    }else{itested();//744
735    }
736  }else{
737  }
738
739  m->tr_eval(this);
740  switch (m->polarity) {
741  case pP:
742    cce = -cce;
743    ccexxx = -ccexxx;
744    cpi = -cpi;
745    cpixxx = -cpixxx;
746    cmu = -cmu;
747    cmuxxx = -cmuxxx;
748    assert(ixxxx == 0.);
749    qbx = -qbx;
750    qbc = -qbc;
751    qcs = -qcs;
752    qbe = -qbe;
753    break;
754  case pN:
755    // leave it as is
756    break;
757  }
758
759  assert(subckt());
760  set_converged(subckt()->do_tr());
761  return converged();
762}
763}
764/*--------------------------------------------------------------------------*/
765/*--------------------------------------------------------------------------*/
766