home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 13 / 13.iso / p / p064 / 3.ddi / BJT.C < prev    next >
Encoding:
C/C++ Source or Header  |  1991-07-01  |  20.7 KB  |  752 lines
  1. /****************************************************************************/
  2. /*        Copyright (C) 1986, 1987, 1988, 1989, 1990, 1991                  */
  3. /*          By MicroSim Corporation, All Rights Reserved                    */
  4. /****************************************************************************/
  5. /* bjt.c
  6.  *   $Revision:   1.17  $
  7.  *   $Author:   pwt  $
  8.  *   $Date:   14 May 1991 17:59:34  $ */
  9.  
  10. /******************* USERS OF DEVICE EQUATIONS OPTION ***********************/
  11. /******** The procedure for changing the bipolar model parameters  **********/
  12. /******** and device equations is the same as for the MOSFET.  See **********/
  13. /******** the comments in the files mos.c and m.h for details.     **********/
  14.  
  15. #define Q_DEVICE
  16.  
  17. #include "option1.h"
  18.  
  19. #ifdef USE_DECL_ARGS
  20. double qsf(double);    /* this declaration will go into some header file */
  21. int Bjt(struct q_ *,int ,int ,int );
  22. #else
  23. double qsf();
  24. int Bjt();
  25. #endif
  26.  
  27. int Bjt(    /* Process bipolar transistor for DC and TRANSIENT analyses */
  28.   Instance,    /* a device to evaluate */
  29.   ModeFl,    /* mode control flag:   see "tran.h" */
  30.   InitFl,    /* initialization flag: see "tran.h" */
  31.   LoadFl    /* NO: bypass MatPrm load if solution unchanged */
  32.   )        /* return: did not converge (YES/NO) */
  33.   struct q_ *Instance;
  34.   int ModeFl, InitFl, LoadFl;
  35.  
  36. /*
  37. pwt    19 Aug 86    creation
  38. pwt    16 Mar 88    add shared data area to Q.H, change source to match
  39. whjb    04 Mar 89    Added substrate current to bypass check
  40. whjb    06 Apr 89    LoadFl==YES forces full recompute
  41. pwt    15 Apr 89    set T=0 cap. currents
  42. pwt    15 Oct 89    enhancements for quasi-saturation effect
  43. pwt    13 Dec 89    add high-injection parameter NK effect
  44. sv      12 Feb 90    add local temperatures to models.
  45. pwt    04 May 90    minor speed-up in some calculations
  46. pwt    08 Aug 90    apply LoadFl condition to parasitic resistors
  47. */
  48. { struct Q_ *model;    /* device model */
  49.   double
  50.     isat,        /* ideal junction saturation current */
  51.     vbc, vbe, vbx, vbn,    /* voltages */
  52.     vce, vjs,        /* vjs is c-to-S for NPN, b-to-S for PNP */
  53.     vtc, vte, vtn,
  54.     gpi, gmu,        /* conductances: active */
  55.     gm, go, gx,
  56.     gbc, gbe, gjs,
  57.     gnbn, gnbc,
  58.     gbcn, gben,
  59.     geq_cb,        /* equiv. cond. of transit time as funct. of Vbc */
  60.     geq_bx,        /* equiv. cond. of b-x capacitance */
  61.     geq_bn,        /* equiv. cond. of b-n capacitance */
  62.     gex, iex,        /* excess phase contributions */
  63.     gxn, gxc,        /* quasi-saturation value */
  64.     kqb,        /* base charge factors */
  65.     dkqb_dvbe,
  66.     dkqb_dvbc,
  67.     ic, ib, in,        /* currents */
  68.     isc, ise, iss,
  69.     ibc, ibe, ijs,
  70.     ibcn, iben,
  71.     ic_hat, ib_hat,    /* predicted currents */
  72.     in_hat,
  73.     type;        /* NPN=+1, PNP=-1 */
  74. #ifdef INC_IJS
  75.   double is_hat;
  76. #endif
  77.  
  78. #define    vt    VT
  79.  
  80.   int
  81.     qsat,
  82.     jlim_fl = YES,    /* junction limiting flag: YES, NO */
  83.     pred_fl  = NO,    /* pred. calculated  flag: YES, NO */
  84.     nonconv  = NO;    /* non-convergence   flag: YES, NO */
  85.  
  86. /* defines & declarations for state-vector access */
  87.  
  88.   struct qsv_def *sv[MSTVCT];
  89.  
  90.   EVAL_SV(sv,Instance,q_sda.q_sv);
  91.  
  92. /* defines for model parameter access: simplify code appearance */
  93.  
  94. #define PARAM(a) ((double)model->a)
  95. #define LPNP    (model->Q_lpnp)
  96. #define AF    PARAM(Q_af)
  97. #define BF    PARAM(Q_bf)
  98. #define BR    PARAM(Q_br)
  99. #define CJC    PARAM(Q_cjc)
  100. #define CJE    PARAM(Q_cje)
  101. #define CJS    PARAM(Q_cjs)
  102. #define EG    PARAM(Q_eg)
  103. #define FC    PARAM(Q_fc)
  104. #define GAMMA    PARAM(Q_gamma)
  105. #define IKF    PARAM(Q_ikf)
  106. #define IKR    PARAM(Q_ikr)
  107. #define IRB    PARAM(Q_irb)
  108. #define IS    PARAM(Q_is)
  109. #define ISC    PARAM(Q_isc)
  110. #define ISE    PARAM(Q_ise)
  111. #define ISS    PARAM(Q_iss)
  112. #define ITF    PARAM(Q_itf)
  113. #define KF    PARAM(Q_kf)
  114. #define MJC    PARAM(Q_mjc)
  115. #define MJE    PARAM(Q_mje)
  116. #define MJS    PARAM(Q_mjs)
  117. #define NC    PARAM(Q_nc)
  118. #define NE    PARAM(Q_ne)
  119. #define NF    PARAM(Q_nf)
  120. #define NK    PARAM(Q_nk)
  121. #define NR    PARAM(Q_nr)
  122. #define NS    PARAM(Q_ns)
  123. #define PTF    PARAM(Q_ptf)
  124. #define QCO    PARAM(Q_qco)
  125. #define RB    PARAM(Q_rb)
  126. #define RBM    PARAM(Q_rbm)
  127. #define RC    PARAM(Q_rc)
  128. #define RCO    PARAM(Q_rco)
  129. #define RE    PARAM(Q_re)
  130. #define TF    PARAM(Q_tf)
  131. #define TR    PARAM(Q_tr)
  132. #define VAF    PARAM(Q_vaf)
  133. #define VAR    PARAM(Q_var)
  134. #define VJC    PARAM(Q_vjc)
  135. #define VJE    PARAM(Q_vje)
  136. #define VJS    PARAM(Q_vjs)
  137. #define VO    PARAM(Q_vo)
  138. #define VTF    PARAM(Q_vtf)
  139. #define XCJC    PARAM(Q_xcjc)
  140. #define XTB    PARAM(Q_xtb)
  141. #define XTF    PARAM(Q_xtf)
  142. #define XTI    PARAM(Q_xti)
  143.  
  144. #define FCPE    PARAM(Q_fcpe)
  145. #define F1BE    PARAM(Q_f1be)
  146. #define F2BE    PARAM(Q_f2be)
  147. #define F3BE    PARAM(Q_f3be)
  148. #define FCPC    PARAM(Q_fcpc)
  149. #define F1BC    PARAM(Q_f1bc)
  150. #define F2BC    PARAM(Q_f2bc)
  151. #define F3BC    PARAM(Q_f3bc)
  152. #define VCRIT    PARAM(Q_vcrit)
  153.  
  154. #define TD    PTF    /* = TF*PF (in radians) */
  155.  
  156. #define IN0    (Instance->qcv_in)
  157. #define IC0    (Instance->qcv_ic)
  158. #define IB0    (Instance->qcv_ib)
  159. #define IJS0    (Instance->qcv_ijs)
  160.  
  161. #define GPI0    (Instance->qcv_gpi)
  162. #define GMU0    (Instance->qcv_gmu)
  163. #define GM0    (Instance->qcv_gm)
  164. #define GO0    (Instance->qcv_go)
  165. #define GX0    (Instance->qcv_gx)
  166. #define GJS0    (Instance->qcv_gjs)
  167. #define GEQBN0    (Instance->qcv_gnbeq)    /* change to ~gbneq for consistency */
  168. #define GEQCB0    (Instance->qcv_gcbeq)
  169. /*#define GEQBX0    (Instance->qcv_gbxeq)    not used */
  170.  
  171. #define GNBN0    (Instance->qcv_gnbn)
  172. #define GNBC0    (Instance->qcv_gnbc)
  173.  
  174. #define CBN0    (Instance->q_sda.q_ac.qac_cbn)
  175. #define CBC0    (Instance->q_sda.q_ac.qac_cbc)
  176. #define CBE0    (Instance->q_sda.q_ac.qac_cbe)
  177. #define CJS0    (Instance->q_sda.q_ac.qac_cjs)
  178. #define CBX0    (Instance->q_sda.q_ac.qac_cBc)
  179.  
  180. #define DevOff    (Instance->q_off)
  181. #define AREA    (Instance->q_area)
  182.  
  183.   model = Instance->q_model;        /* find the model */
  184.  
  185.   type  = model->Q_type > 0 ? 1.: -1.;
  186.  
  187. /* initialization */
  188.  
  189.   geq_bx = geq_cb = geq_bn = 0.;
  190.  
  191.   qsat = RCO!=0.? YES : NO;
  192.  
  193.   if      ( InitFl==INSTV0 ) {        /* use prev. iteration */
  194.     double
  195.       vb = VOLTAGE(q_b),
  196.       vc = VOLTAGE(q_c);
  197.  
  198.     vbe = type*( vb - VOLTAGE(q_e) );
  199.     vbc = type*( vb - vc );
  200.     vbn = type*( vb - VOLTAGE(q_n) );
  201.     vbx = type*( VOLTAGE(q_B) - vc );
  202.     vjs = LPNP
  203.         ? type*( vb - VOLTAGE(q_S) )
  204.         : type*( VOLTAGE(q_S) - vc );
  205.     }
  206.   else if ( InitFl==INTRAN ) {        /* use prev. time-step */
  207.     vbe = Q_VBE(1);
  208.     vbc = Q_VBC(1);
  209.     vbn = Q_VBN(1);
  210.     vjs = Q_VJS(1);
  211.     if ( ModeFl==MDTRAN && NOSOLV==YES ) {
  212.       vbx = type*( Instance->q_vbe - Instance->q_vce );
  213.       }
  214.     else {
  215.       vbx = type*( VOLTAGE(q_B) - VOLTAGE(q_c) );
  216.       }
  217.     }
  218.   else if ( InitFl==INOFF &&        /* set "off" devices */
  219.             DevOff==YES)
  220.     vbe = vbc = vbn = vjs = vbx = 0.;
  221.  
  222.   else if ( InitFl==ININIT ) {        /* use IC= values */
  223.  
  224.     if ( ModeFl==MDBPTR && NOSOLV==YES ) {
  225.       vce = type*Instance->q_vce;
  226.       vbe = type*Instance->q_vbe;
  227.       vbx = vbc = vbn = vbe-vce;
  228.       vjs = 0.;
  229.       }
  230.     else {
  231.       vbe = ( DevOff ? 0. : VCRIT );
  232.       vbc = vbn = vbx = vjs = 0.;
  233.       }
  234.     }
  235.   else {
  236.     double
  237.       del_vbe, del_vbc, del_vbn, del_vjs, del_in, del_ic, del_ib;
  238. #ifdef INC_IJS
  239.     double del_is;
  240. #endif
  241.  
  242.     if ( InitFl==INPRDCT ) {        /* extrapolate value */
  243.       double arg = DELTA/DELOLD[1];
  244.  
  245.       vbe = arg*( Q_VBE(1) - Q_VBE(2) ) + Q_VBE(1);
  246.       vbc = arg*( Q_VBC(1) - Q_VBC(2) ) + Q_VBC(1);
  247.       vbn = arg*( Q_VBN(1) - Q_VBN(2) ) + Q_VBN(1);
  248.       vjs = arg*( Q_VJS(1) - Q_VJS(2) ) + Q_VJS(1);
  249.  
  250.       *sv[0] = *sv[1];
  251.  
  252.       vbx = type*( VOLTAGE(q_B) - VOLTAGE(q_c) );
  253.       }
  254.     else {                /* use current value */
  255.       double
  256.         vb = VOLTAGE(q_b),
  257.         vc = VOLTAGE(q_c);
  258.  
  259.       vbe = type*( vb - VOLTAGE(q_e) );
  260.       vbc = type*( vb - vc );
  261.       vbn = type*( vb - VOLTAGE(q_n) );
  262.       vbx = type*( VOLTAGE(q_B) - vc );
  263.       vjs = LPNP
  264.           ? type*( vb - VOLTAGE(q_S) )
  265.           : type*( VOLTAGE(q_S) - vc );
  266.       }
  267.  
  268. /* compute new non-linear branch voltage */
  269.  
  270.     del_vbe = vbe - Q_VBE(0);
  271.     del_vbc = vbc - Q_VBC(0);
  272.     del_vbn = vbn - Q_VBN(0);
  273.     del_vjs = vjs - Q_VJS(0);
  274.  
  275.     del_ic =   del_vbe*(GM0+GO0)
  276.              - del_vbc*(GMU0+GNBC0+GO0)
  277.              - del_vbn*GNBN0
  278.              - (   LPNP ? 0. :  del_vjs*GJS0 );
  279. /*           - ( MDTRAN ? del_vbx*GEQBX0 : 0.);    note: del_vbx not available */
  280.  
  281.     ic_hat = del_ic + IC0;
  282.  
  283.     del_ib =   del_vbe*GPI0
  284.              + del_vbc*(GMU0+GEQCB0)
  285.              + (  LPNP  ? del_vjs*GJS0   : 0.)
  286.              + ( MDTRAN ? del_vbn*GEQBN0 : 0.);
  287.  
  288.     ib_hat = del_ib + IB0;
  289.  
  290.     if ( qsat ) {
  291.       del_in =   del_vbc*GNBC0
  292.                + del_vbn*( GNBN0 - ( MDTRAN ? GEQBN0 : 0.));
  293.       in_hat = del_in + IN0;
  294.       }
  295.     else
  296.       in_hat = del_in = 0.;
  297.  
  298. #ifdef INC_IJS
  299.     del_is = del_vjs*GJS0;
  300.     is_hat = del_is + IJS0;
  301. #endif
  302.  
  303. /* bypass if solution not changed */
  304.  
  305.     if ( LoadFl ) ;
  306.     else if ( InitFl==INPRDCT && DELTA!=DELOLD[1] ) ;
  307.     else if ( qsat &&
  308.               fabs(del_in)  >= TOL(in_hat,     IN0,RELTOL,ABSTOL) ) ;
  309.     else if ( fabs(del_ic)  >= TOL(ic_hat,     IC0,RELTOL,ABSTOL) ) ;
  310.     else if ( fabs(del_vbc) >= TOL(   vbc,Q_VBC(0),RELTOL, VNTOL) ) ;
  311.     else if ( fabs(del_vbe) >= TOL(   vbe,Q_VBE(0),RELTOL, VNTOL) ) ;
  312.     else if ( qsat &&
  313.               fabs(del_vbn) >= TOL(   vbn,Q_VBN(0),RELTOL, VNTOL) ) ;
  314.     else if ( fabs(del_vjs) >= TOL(   vjs,Q_VJS(0),RELTOL, VNTOL) ) ;
  315.     else if ( fabs(del_ib)  >= TOL(ib_hat,     IB0,RELTOL,ABSTOL) ) ;
  316. #ifdef INC_IJS
  317.     else if ( fabs(del_is)  >= TOL(is_hat,    IJS0,RELTOL,ABSTOL) ) ;
  318. #endif
  319.     else {
  320.       goto done;
  321.       }
  322.     pred_fl = YES;
  323.  
  324. /* solution changed: limit new junction voltages */
  325.  
  326.     jlim_fl  = PNJLIM(vbe,Q_VBE(0),vt,VCRIT);
  327.     jlim_fl |= PNJLIM(vbc,Q_VBC(0),vt,VCRIT);
  328.  
  329.     if (  qsat  ) jlim_fl |= PNJLIM(vbn,Q_VBN(0),vt,VCRIT);
  330.     if ( ISS!=0.) jlim_fl |= PNJLIM(vjs,Q_VJS(0),vt,VCRIT);
  331.       
  332.     } /* end of initialization */
  333.  
  334. /* compute DC current & derivatives */
  335.  
  336.   isat = AREA*IS;
  337.   ise  = AREA*ISE;
  338.   isc  = AREA*ISC;
  339.   iss  = AREA*ISS;
  340.   vte  = NE*vt;
  341.   vtc  = NC*vt;
  342.   vtn  = NF*vt;
  343.  
  344.   if ( vbe > -10*vtn ) {        /* b-e "active" */
  345.     double evbe = EXP(vbe/vtn);
  346.  
  347.     gbe = (evbe/vtn)*isat + GMIN;
  348.     ibe = (evbe-1)*isat + vbe*GMIN;
  349.  
  350.     if ( ise==0. )
  351.       gben = iben = 0.;
  352.     else {
  353.       double evben = EXP(vbe/vte);
  354.  
  355.       gben = (evben/vte)*ise;
  356.       iben = (evben-1)*ise;
  357.     } }
  358.   else {                /* b-e reverse bias */
  359.     gbe  = GMIN;
  360.     ibe  = vbe*GMIN - isat;
  361.     gben = 0.;
  362.     iben = ( ise==0. ? 0. : -ise );
  363.     }
  364.  
  365.   vtn = NR*vt;
  366.   if ( vbc > -10*vtn ) {        /* b-c "active" */
  367.     double evbc = EXP(vbc/vtn);
  368.  
  369.     gbc = (evbc/vtn)*isat + GMIN;
  370.     ibc = (evbc-1)*isat + vbc*GMIN;
  371.  
  372.     if ( isc==0. )
  373.       gbcn = ibcn = 0.;
  374.     else {
  375.       double evbcn = EXP(vbc/vtc);
  376.  
  377.       gbcn = (evbcn/vtc)*isc;
  378.       ibcn = (evbcn-1)*isc;
  379.     } }
  380.   else {                /* b-c reverse bias */
  381.     gbc  = GMIN;
  382.     ibc  = vbc*GMIN - isat;
  383.     gbcn = 0.;
  384.     ibcn = ( isc==0. ? 0. : -isc );
  385.     }
  386.  
  387.   if ( iss==0. ) {            /* no diode (SPICE compatible) */
  388.     gjs = ijs = 0.;
  389.     }
  390.   else {
  391.     vtn = NS*vt;
  392.     if ( vjs > -10*vtn ) {        /* substrate junction "active" */
  393.       double evjs = EXP(vjs/vtn);
  394.  
  395.       gjs = (evjs/vtn)*iss + GMIN;
  396.       ijs = (evjs-1)*iss + vjs*GMIN;
  397.       }
  398.     else {                /* substrate junction reverse bias */
  399.       gjs = GMIN;
  400.       ijs = vjs*GMIN - iss;
  401.     } }
  402.  
  403.   { double                /* determine base charge terms */
  404.       kq1 = 1/(1 - vbc*VAF - vbe*VAR);
  405.  
  406.     kqb = kq1;
  407.     dkqb_dvbc = kq1*kqb*VAF;
  408.     dkqb_dvbe = kq1*kqb*VAR;
  409.  
  410.     if ( IKF!=0. || IKR!=0.) {
  411.       double arg, darg,            /* due to MODCHK: */
  412.         oikf = IKF/AREA,        /* IKF is really 1/ikf */
  413.         oikr = IKR/AREA,        /* IKR is really 1/ikr */
  414.         kq2  = ibe*oikf + ibc*oikr;
  415.  
  416.       arg = 1+4*kq2;
  417.       if ( arg<=0.) darg = arg = 1.;    /* erronious situation */
  418.       else {
  419.         if ( NK==.5 )
  420.           darg = ( arg = sqrt( darg = arg) )/darg;
  421.         else
  422.           darg = ( arg = pow(  darg = arg, NK) )/darg;
  423.         }
  424.  
  425.       kqb       = kq1*(1+arg)*.5;
  426.       dkqb_dvbc = kq1*(kqb*VAF + oikr*gbc*darg);
  427.       dkqb_dvbe = kq1*(kqb*VAR + oikf*gbe*darg);
  428.     } }
  429.  
  430.   if ( ModeFl==MDTRAN && TD!=0.) {    /* Weil's approximation for */
  431.     double                /* excess phase applied with */
  432.       arg1  = DELTA/TD,            /* backward Euler integration */
  433.       arg2  = 3*arg1,
  434.       denom = ( arg1 *= arg2 ) + arg2 + 1,
  435.       arg3  = arg1/denom;
  436.  
  437.     gex = arg3*gbe;
  438.     iex = arg3*ibe;
  439.  
  440.     if ( InitFl==INTRAN )        /* set-up initial current */
  441.       Q_IEXBC(2) = Q_IEXBC(1) = ibe/kqb;
  442.  
  443.     ic = ( Q_IEXBC(1)*( 1 + DELTA/DELOLD[1] + arg2 ) -
  444.            Q_IEXBC(2)*DELTA/DELOLD[1] )/denom;
  445.  
  446.     Q_IEXBC(0) = ic + iex/kqb;
  447.     }
  448.   else {
  449.     ic  = 0.;
  450.     gex = gbe;
  451.     iex = ibe;
  452.     }
  453.  
  454.   { double                /* combine currents and determine */
  455.       arg1 = ibc/BR + ibcn,        /* DC incremental conductances */
  456.       arg2 = (iex-ibc)/kqb;
  457.  
  458.     ic += arg2 - arg1;
  459.     ib  = ibe/BF + iben + arg1;
  460.  
  461.     gpi = gbe/BF + gben;
  462.     gmu = gbc/BR + gbcn;
  463.     go  = (gbc + arg2*dkqb_dvbc)/kqb;
  464.     gm  = (gex - arg2*dkqb_dvbe)/kqb - go;
  465.     }
  466.  
  467.   if ( RB==0.) gx = 0.;            /* determine base resistance */
  468.   else {
  469.     if ( IRB==0.) gx = RBM + (RB-RBM)/kqb;
  470.     else {
  471.       double
  472.         arg1 = ib/(AREA*IRB),
  473.         arg2;
  474.  
  475.       if ( arg1<=1E-16 ) gx = RB;    /* too small re. 1.0 */
  476.       else {
  477.         arg1 = arg1*14.59025044449664;    /* (ib/ibr)*(12/pi)^2 */
  478.  
  479.         arg2 = .25*TWOPI*(sqrt(1+arg1)-1)/sqrt(arg1);
  480.  
  481.         arg1 = tan(arg2);
  482.  
  483.         gx   = RBM + (RB-RBM)*(arg1-arg2)*3/(arg2*arg1*arg1);
  484.       } }
  485.     gx = AREA/gx;
  486.     }
  487.  
  488.   if ( qsat ) {                /* quasi-saturation */
  489.     double
  490.       rco = RCO/AREA;
  491.  
  492.     gxc = GAMMA*exp(vbc/vt);
  493.     gxn = GAMMA*exp(vbn/vt);
  494.  
  495.     if ( vbc==vbn ) {
  496.       in   = 0.;
  497.       gnbn = -(
  498.       gnbc = (gxn/(2*(qsf(gxn)+2)) + 1)/rco );
  499.       }
  500.     else {
  501.       double
  502.         vnc  = vbc-vbn,
  503.         tmp  = fabs(vnc)+VO,
  504.         rtmp = rco*tmp,
  505.         qsbn = qsf(gxn),
  506.         qsbc = qsf(gxc),
  507.         sign = vnc>0.? 1.: -1.;
  508.  
  509.       in   = ( vt*( (qsbc-qsbn) + log((qsbn+2)/(qsbc+2)) ) + vnc )*VO/rtmp;
  510.  
  511.       gnbc = (gxc/(2*(qsbc+2)) + 1)*VO/rtmp - sign*in/tmp;
  512.  
  513.       gnbn = sign*in/tmp - (gxn/(2*(qsbn+2)) + 1)*VO/rtmp;
  514.     } }
  515.   else
  516.     in = gnbc = gnbn = 0.;
  517.  
  518. /* charge calculations */
  519.  
  520.   if ( ModeFl==MDTRAN ||
  521.        InitFl==INSTV0 ||
  522.      ( ModeFl==MDBPTR && NOSOLV==YES ) ) {
  523.     double
  524.       cbn, cbc, cbe, cbx, cjs,
  525.       cbeo = AREA*CJE,        /* scale 0-bias cap. */
  526.       ctot = AREA*CJC,
  527.       cbco = ctot*XCJC,
  528.       cbxo = ctot-cbco,
  529.       cjso = AREA*CJS,
  530.       itf  = AREA*ITF;
  531.  
  532.     if ( TF!=0. && vbe > 0.) {
  533.       double argtf = 0., arg2 = 0., arg3 = 0.;
  534.  
  535.       if ( XTF!=0.) {
  536.         argtf = XTF;
  537.         if ( VTF!=0.) argtf *= exp(vbc*VTF);
  538.         arg2 = argtf;
  539.         if ( itf!=0. ) {
  540.           double temp = ibe/(ibe+itf);
  541.  
  542.           argtf = argtf*temp*temp;
  543.           arg2  = argtf*(3-temp-temp);
  544.           }
  545.         arg3 = ibe*argtf*VTF;
  546.         }
  547.       ibe    = ibe*(1+argtf)/kqb;
  548.       gbe    = (gbe*(1+arg2) - ibe*dkqb_dvbe)/kqb;
  549.       geq_cb = (arg3-ibe*dkqb_dvbc)*TF/kqb;
  550.       }
  551.  
  552.     PNJCHG(Q_QCBE(0),cbe,cbeo,vbe,VJE,FCPE,MJE,F1BE,F2BE,F3BE);
  553.     Q_QCBE(0) += ibe*TF;
  554.     cbe       += gbe*TF;
  555.  
  556.     PNJCHG(Q_QCBC(0),cbc,cbco,vbc,VJC,FCPC,MJC,F1BC,F2BC,F3BC);
  557.     Q_QCBC(0) += ibc*TR;
  558.     cbc       += gbc*TR;
  559.  
  560.     PNJCHG(Q_QCBX(0),cbx,cbxo,vbx,VJC,FCPC,MJC,F1BC,F2BC,F3BC);
  561.  
  562.     PNJCHG0(Q_QCJS(0),cjs,cjso,vjs,VJS,MJS);
  563.  
  564.     if ( qsat ) {            /* quasi-saturation */
  565.       double
  566.         qco = AREA*QCO,
  567.         tmp = .5*qco/vt;
  568.  
  569.       Q_QCBN(0)  = qco*(qsf(gxn)-.5*GAMMA);
  570.       Q_QCBC(0) += qco*(qsf(gxc)-.5*GAMMA);
  571.  
  572.       cbn  = gxn*tmp/sqrt(1+gxn);
  573.       cbc += gxc*tmp/sqrt(1+gxc);
  574.       }
  575.     else
  576.       cbn = Q_QCBN(0) = 0.;
  577.  
  578.     Q_ICBN(0) =        /* reset, so T=0 results (e.g. Probe) are correct */
  579.     Q_ICBC(0) =
  580.     Q_ICBE(0) =
  581.     Q_ICJS(0) =
  582.     Q_ICBX(0) = 0.;
  583.  
  584.     if ( InitFl==INSTV0) {        /* store small-signal parameters */
  585.       CBN0 = cbn;            /* and update the state vector */            
  586.       CBC0 = cbc;
  587.       CBE0 = cbe;
  588.       CJS0 = cjs;
  589.       CBX0 = cbx;
  590.  
  591.       goto load;
  592.       }
  593.  
  594.     if ( !(ModeFl==MDBPTR && NOSOLV==NO) ) {    /* transient analysis */
  595.       double g_eq, i_eq;
  596.  
  597.       if ( InitFl==INTRAN ) {            /* for 1st transient step: */
  598.         Q_QCBN(2) = Q_QCBN(1) = Q_QCBN(0);    /* set "old" charges so */
  599.         Q_QCBC(2) = Q_QCBC(1) = Q_QCBC(0);    /* integration works */
  600.         Q_QCBE(2) = Q_QCBE(1) = Q_QCBE(0);
  601.         Q_QCJS(2) = Q_QCJS(1) = Q_QCJS(0);
  602.         Q_QCBX(2) = Q_QCBX(1) = Q_QCBX(0);
  603.         }
  604.       geq_cb *= AG1;
  605.  
  606.       INTEGR8(geq_bn,i_eq,cbn,Q_QIBN(0),Q_QIBN(1),Q_QIBN(2));
  607.  
  608.       INTEGR8(g_eq,  i_eq,cbc,Q_QIBC(0),Q_QIBC(1),Q_QIBC(2));
  609.       gmu +=  g_eq;
  610.  
  611.       INTEGR8(g_eq,  i_eq,cbe,Q_QIBE(0),Q_QIBE(1),Q_QIBE(2));
  612.       gpi +=  g_eq;
  613.  
  614.       INTEGR8(g_eq,  i_eq,cjs,Q_QIJS(0),Q_QIJS(1),Q_QIJS(2));
  615.       gjs +=  g_eq;
  616.       ijs +=  Q_ICJS(0);
  617.  
  618.       INTEGR8(geq_bx,i_eq,cbx,Q_QIBX(0),Q_QIBX(1),Q_QIBX(2));
  619.  
  620.       ib += Q_ICBC(0) + Q_ICBE(0);
  621.       ic -= Q_ICBC(0);
  622.  
  623.       if ( InitFl==INTRAN ) {            /* for 1st transient step: */
  624.         Q_ICBN(1) = Q_ICBN(0);            /* set "old" currents so */
  625.         Q_ICBC(1) = Q_ICBC(0);            /* time-step calc. works */
  626.         Q_ICBE(1) = Q_ICBE(0);
  627.         Q_ICJS(1) = Q_ICJS(0);
  628.         Q_ICBX(1) = Q_ICBX(0);
  629.         }
  630.     } } /* end of charge calculations */
  631.  
  632. /* check convergence */
  633.  
  634.   if ( pred_fl && (
  635.          jlim_fl ||
  636.  
  637.          qsat && fabs(in_hat-in) > TOL(in_hat,in,RELTOL,ABSTOL) ||
  638.  
  639.          fabs(ic_hat-ic) > TOL(ic_hat,ic,RELTOL,ABSTOL) ||
  640.          fabs(ib_hat-ib) > TOL(ib_hat,ib,RELTOL,ABSTOL)
  641.          )
  642.        ) nonconv = YES;
  643.  
  644. load:
  645.   { double
  646.       gcpr, gepr,
  647.       ieq_nc = qsat ? type*( in - vbc*gnbc - vbn*gnbn ) : 0.,
  648.       ieq_cb = type*( ic - vbe*(gm+go) + vbc*(gmu+go) ),
  649.       ieq_be = type*( ic + ib - vbe*(gm+go+gpi) + vbc*(go-geq_cb) ),
  650.       ieq_js = type*( ijs - vjs*gjs )*( LPNP ? -1 : 1 ),
  651.       ieq_bx = ModeFl==MDTRAN ? type*( Q_ICBX(0) - vbx*geq_bx ) : 0.,
  652.       ieq_bn = ModeFl==MDTRAN ? type*( Q_ICBN(0) - vbn*geq_bn ) : 0.;
  653.  
  654.     /* load current vector */
  655.  
  656.     Y_MATRIX(q_iB) = -ieq_bx;
  657.     Y_MATRIX(q_ic) = ( LPNP ? 0.: ieq_js ) - ieq_cb + ieq_nc + ieq_bx;
  658.     Y_MATRIX(q_ib) = ( LPNP ? ieq_js : 0.) + ieq_cb - ieq_be - ieq_bn;
  659.     Y_MATRIX(q_ie) =  ieq_be;
  660.     Y_MATRIX(q_iS) = -ieq_js;
  661.     Y_MATRIX(q_in) =  ieq_bn - ieq_nc;
  662.  
  663.     /* load conductance matrix */
  664.  
  665.     Y_MATRIX(q_Bb) = Y_MATRIX(q_bB) = -gx;
  666.     Y_MATRIX(q_Bc) = Y_MATRIX(q_cB) = -geq_bx;
  667.     Y_MATRIX(q_BB) = gx + geq_bx;
  668.  
  669.     Y_MATRIX(q_bn) = -geq_bn;
  670.     Y_MATRIX(q_bc) = -gmu - geq_cb;
  671.     Y_MATRIX(q_be) = -gpi;
  672.     Y_MATRIX(q_bb) = ( LPNP ? gjs : 0.) + gx + gmu + gpi + geq_cb + geq_bn;
  673.  
  674.     Y_MATRIX(q_cn) =  gnbn;
  675.     Y_MATRIX(q_nc) = -gnbc;
  676.     Y_MATRIX(q_nb) =  gnbn - geq_bn + gnbc;
  677.     Y_MATRIX(q_nn) = -gnbn + geq_bn + ( gcpr = AREA*RC );
  678.  
  679.     Y_MATRIX(q_cb) =  gm - gmu - gnbn - gnbc;
  680.     Y_MATRIX(q_ce) = -gm - go;
  681.     Y_MATRIX(q_cc) = ( LPNP ? 0.: gjs ) + gmu + go + geq_bx + gnbc;
  682.  
  683.     Y_MATRIX(q_eb) = -gm - gpi - geq_cb;
  684.     Y_MATRIX(q_ec) = -go + geq_cb;
  685.     Y_MATRIX(q_ee) =  gm + go + gpi + ( gepr = AREA*RE );
  686.  
  687.     Y_MATRIX(q_Sj) = Y_MATRIX(q_jS) = -( Y_MATRIX(q_SS) = gjs );
  688.  
  689.     if ( LoadFl ) {
  690.       Y_MATRIX(q_Cn) = Y_MATRIX(q_nC) = -( Y_MATRIX(q_CC) = gcpr );
  691.       Y_MATRIX(q_Ee) = Y_MATRIX(q_eE) = -( Y_MATRIX(q_EE) = gepr );
  692.     } }
  693.  
  694. #ifdef DEBUG
  695. printf("\n*** %s at time = %g",Instance->q_name,TIME);
  696. printf("\n  q_Bb = %g",Y_MATRIX(q_Bb));
  697. printf("\n  q_bB = %g",Y_MATRIX(q_bB));
  698. printf("\n  q_BB = %g",Y_MATRIX(q_BB));
  699. printf("\n  q_bb = %g",Y_MATRIX(q_bb));
  700. printf("\n  q_Bc = %g",Y_MATRIX(q_Bc));
  701. printf("\n  q_bc = %g",Y_MATRIX(q_bc));
  702. printf("\n  q_be = %g",Y_MATRIX(q_be));
  703. printf("\n  q_cB = %g",Y_MATRIX(q_cB));
  704. printf("\n  q_cb = %g",Y_MATRIX(q_cb));
  705. printf("\n  q_Cc = %g",Y_MATRIX(q_Cc));
  706. printf("\n  q_cC = %g",Y_MATRIX(q_cC));
  707. printf("\n  q_CC = %g",Y_MATRIX(q_CC));
  708. printf("\n  q_cc = %g",Y_MATRIX(q_cc));
  709. printf("\n  q_ce = %g",Y_MATRIX(q_ce));
  710. printf("\n  q_eb = %g",Y_MATRIX(q_eb));
  711. printf("\n  q_ec = %g",Y_MATRIX(q_ec));
  712. printf("\n  q_Ee = %g",Y_MATRIX(q_Ee));
  713. printf("\n  q_eE = %g",Y_MATRIX(q_eE));
  714. printf("\n  q_EE = %g",Y_MATRIX(q_EE));
  715. printf("\n  q_ee = %g",Y_MATRIX(q_ee));
  716. printf("\n  q_iB = %g",Y_MATRIX(q_iB));
  717. printf("\n  q_ib = %g",Y_MATRIX(q_ib));
  718. printf("\n  q_ic = %g",Y_MATRIX(q_ic));
  719. printf("\n  q_ie = %g",Y_MATRIX(q_ie));
  720. printf("\n  q_iS = %g",Y_MATRIX(q_iS));
  721. printf("\n  q_jS = %g",Y_MATRIX(q_jS));
  722. printf("\n  q_Sj = %g",Y_MATRIX(q_Sj));
  723. printf("\n  q_SS = %g",Y_MATRIX(q_SS));
  724. #endif
  725.  
  726.   Q_VBE(0) = vbe;
  727.   Q_VBC(0) = vbc;
  728.   Q_VJS(0) = vjs;
  729.   Q_VBN(0) = vbn;
  730.  
  731.   IN0  = in;
  732.   IC0  = ic;
  733.   IB0  = ib;
  734.   IJS0 = ijs;
  735.  
  736.   GPI0 = gpi;
  737.   GMU0 = gmu;
  738.   GM0  = gm;
  739.   GO0  = go;
  740.   GX0  = gx;
  741.   GJS0 = gjs;
  742.  
  743.   GEQBN0 = geq_bn;
  744.   GEQCB0 = geq_cb;
  745.  
  746.   GNBN0  = gnbn;
  747.   GNBC0  = gnbc;
  748.  
  749. done:
  750.   return nonconv;
  751.   }
  752.