home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 13 / 13.iso / p / p064 / 3.ddi / MOS2.C < prev    next >
Encoding:
C/C++ Source or Header  |  1991-07-01  |  21.5 KB  |  771 lines
  1. /****************************************************************************/
  2. /*        Copyright (C) 1986, 1987, 1988, 1989, 1990, 1991                  */
  3. /*          By MicroSim Corporation, All Rights Reserved                    */
  4. /****************************************************************************/
  5. /* mos2.c
  6.  *   $Revision:   1.16  $
  7.  *   $Author:   pwt  $
  8.  *   $Date:   26 Jun 1991 17:24:12  $ */
  9.  
  10. #include "option1.h"
  11. #include "mos.fd"
  12.  
  13. /* these are set in MOS2eval and used there, and later in MOS2cap */
  14.  
  15. static double
  16.     Cox,
  17.     Gamasd,
  18.     Qspof,
  19.     Vbin,
  20.     Vth;
  21.  
  22. #define PARAM(a)     ((double)Model->a)
  23. #define DELTA_MOS    PARAM(M_delta)
  24. #define GAMMA        PARAM(M_gamma)
  25. #define KP        PARAM(M_kp)
  26. #define LAMBDA        PARAM(M_lambda)
  27. #define NEFF        PARAM(M_neff)
  28. #define NFS        PARAM(M_nfs)
  29. #define NSUB        PARAM(M_nsub)
  30. #define PB        PARAM(M_pb)
  31. #define PHI        PARAM(M_phi)
  32. #define TOX        PARAM(M_tox)
  33. #define TYPE        (Model->M_type)
  34. #define UCRIT        PARAM(M_ucrit)
  35. #define UEXP        PARAM(M_uexp)
  36. #define UO        PARAM(M_uo)
  37. #define VBI        PARAM(M_vbi)
  38. #define VBP        UCRIT        /* = bypass voltage (level 2) */
  39. #define VMAX        PARAM(M_vmax)
  40. #define XD        PARAM(M_xd)
  41. #define XJ        PARAM(M_xj)
  42.  
  43. #ifndef cc_noline
  44. #line 4 "MOS2eval"
  45. #endif
  46.  
  47. void MOS2eval(Model,Len,Wid,Mdev,Vgs,Vds,Vbs,Von,Vdsat,Id,Gds,Gm,Gmbs,/*TITLE*/
  48.           Charge,Cggb,Cgdb,Cgsb,Cbgb,Cbdb,Cbsb,Qgate,Qchan,Qbulk)
  49.  
  50. struct M_ *Model;          /* (pointer to) model */
  51. double    Len, Wid, Mdev,      /* effective device size */
  52.         Vgs, Vds, Vbs,      /* terminal voltages */
  53.         *Von, *Vdsat,            /* (pointer to) device voltages (returned) */
  54.         *Id,              /* (pointer to) drain current (returned) */
  55.         *Gds, *Gm, *Gmbs;      /* (pointer to) device conductances (returned) */
  56. int        Charge;          /* do charge calculations (YES/NO) */
  57. double      *Cggb, *Cgdb, *Cgsb,      /* (pointer to) device capacitances (returned) */
  58.       *Cbgb, *Cbdb, *Cbsb,
  59.       *Qgate, *Qchan, *Qbulk; /* (pointer to) device charges (returned) */
  60.  
  61. /*****************************************************************************
  62. * Purpose
  63. *   evaluate MOS level 2 (analytical, one-dimensional) model
  64. *
  65. * History
  66. * 86-08-22 pwt    - creation
  67. * 87-10-07 pwt    - re-work as part of BSIM addition
  68. * 88-01-05 pwt    - parameter M_tox now oxide thickness, not Cox
  69. * 89-01-25 swg    - limit v2 and v1 in quartic eqn. solution, even more
  70. *          than they were, to prevent VAX floating point overflow.
  71. *          Also limit sarg3.
  72. *        - Init Gmbs in one case where it wasn't
  73. * 89-04-15 pwt    - add device multiplier, to fix narrow-channel calculation
  74. * 89-06-12 whjb - Made use of EXP() macro
  75. * 90-02-12 sv    - add local temperatures to models.
  76. * whjb 04 Mar 91 - Corrected bug in weak inversion and in mobility modulation
  77. * whjb 08 Jun 91 - Corrected bug in derivative of "body" in linear and sat.
  78. * pwt  26 Jun 91 - remove divide-by-zero problem, analytically, in sub-thres.
  79. *
  80. *****************************************************************************/
  81. {
  82. /* Local variables */
  83.   int    ivmflg;            /* mystery counter */
  84.  
  85.   double
  86.     temp,
  87.     lambda,
  88.     vgst, vgsx, vpof, vdsat1,
  89.     eta, factor,
  90.     beta1,            /* effective beta (kp') */
  91.     ufact, ueff,        /* mobility */
  92.     dudvgs, dudvbs, dudvds,
  93.     xn, dxndvb, dxndvd,    /* weak inversion: Von = Vth + xn*VT */
  94.     dodvbs,    dodvds,
  95.     argg,            /* = 1/(xn*VT) */
  96.     body, dbdvds,
  97.     sphi, sphi3,        /* phi^(1/2) & phi^(3/2) */
  98.     sarg, barg,        /* charge effects */
  99.     sarg3,
  100.     dsrgdb, d2sdb2,        /* 1st & 2nd derivatives */
  101.     dbrgdb, d2bdb2,
  102.     gammad,
  103.     dgddvb, dgdvds, dgdvbs,
  104.     dgddb2,
  105.     dsdvgs, dsdvbs,
  106.     bsarg, dbsrdb,        /* effective channel length args. */
  107.     bodys, gdbdvs,
  108.     dldvgs, dldvds, dldvbs,
  109.     clfact,
  110.     gdbdv;
  111. #define    vt    VT
  112.  
  113.   static double        /* arrays for scattering limited velocity evaluation */
  114.     sig1[] = { 0., 1., -1.,  1., -1. },
  115.     sig2[] = { 0., 1.,  1., -1., -1. };
  116.  
  117. /* Code */
  118.   Cox  = Wid*Len*EPSOX/TOX;
  119.  
  120. /* compute some useful quantities */
  121.  
  122.   sphi   = sqrt(PHI);
  123.   sphi3  = sphi*PHI;
  124.   lambda = LAMBDA;            /* gets re-evaluated (sometimes) */
  125.  
  126.   if ( Vbs > 0. ) {
  127.     sarg   = sphi/(1 + .5*Vbs/PHI);
  128.     temp   = -sarg/sphi3;
  129.     dsrgdb = .5*temp*sarg;
  130.     d2sdb2 = dsrgdb*temp;
  131.     }
  132.   else {
  133.     temp   = PHI-Vbs;
  134.     sarg   = sqrt(temp);
  135.     dsrgdb = -.5/sarg;
  136.     d2sdb2 =  .5*dsrgdb/temp;
  137.     }
  138.  
  139.   if ( Vds-Vbs < 0. ) {
  140.     barg   = sphi/(1 + .5*(Vbs-Vds)/PHI);
  141.     temp   = -barg/sphi3;
  142.     dbrgdb = .5*temp*barg;
  143.     d2bdb2 = dbrgdb*temp;
  144.     }
  145.   else {
  146.     temp   = PHI+Vds-Vbs;
  147.     barg   = sqrt(temp);
  148.     dbrgdb = -.5/barg;
  149.     d2bdb2 =  .5*dbrgdb/temp;
  150.     }
  151.  
  152. /* calculate threshold voltage */
  153.  
  154.   if ( GAMMA <= 0. || NSUB  <= 0. ) {    /* no short-channel effects */
  155.     Gamasd = gammad = GAMMA;
  156.     dgddvb = dgdvds = dgddb2 = 0.;
  157.     }
  158.   else {                /* short channel effect w/Vds != 0. */
  159.     if ( XJ > 0.) {
  160.       double argxs, argss, args;
  161.       double argxd, argsd, argd;
  162.       double dbxws, dbxwd;
  163.       double dbargs, dbargd;
  164.       double xws, xwd;
  165.  
  166.       xws    = XD*sarg;            /* depletion charge effect: source */
  167.       argxs  = 1+2*xws/XJ;
  168.       args   = sqrt(argxs);
  169.       argss  = (args-1)*.5*XJ/Len;
  170.  
  171.       xwd    = XD*barg;            /* depletion charge effect: drain */
  172.       argxd  = 1+2*xwd/XJ;
  173.       argd   = sqrt(argxd);
  174.       argsd  = (argd-1)*.5*XJ/Len;
  175.       Gamasd = (1-argsd-argss)*GAMMA;
  176.  
  177.       dbxws  = XD*dsrgdb;        /* derivatives of corrections */
  178.       dbargs = .5*dbxws/(Len*args);
  179.  
  180.       dbxwd  = XD*dbrgdb;
  181.       dbargd = .5*dbxwd/(Len*argd);
  182.       dgdvds = dbargd*GAMMA;
  183.  
  184.       dgddb2 = .5*GAMMA*XD/Len*(
  185.                  ( d2sdb2 + dsrgdb*dsrgdb*XD/(XJ*argxs) )/args
  186.                  +
  187.                  ( d2bdb2 + dbrgdb*dbrgdb*XD/(XJ*argxd) )/argd
  188.                  );
  189.  
  190.       dgddvb = -GAMMA*(dbargs+dbargd);
  191.       }
  192.     else {
  193.       dgdvds = dgddb2 = dgddvb = 0.;
  194.       Gamasd = gammad = GAMMA;
  195.     } }
  196.  
  197. /*** this section only useful if DELTA > 0 ***/
  198. /*factor = .125*DELTA_MOS*TWOPI*EPSSIL/Cox*Len; <- confusing, but correct */
  199.   factor = (TWOPI*DELTA_MOS*EPSSIL*TOX)/(8*EPSOX*Wid/Mdev);
  200.     /* narrow channel effect */
  201.   eta    = 1+factor;
  202.   Vbin   = TYPE*VBI + factor*(PHI-Vbs);
  203.   *Von   = Vbin + Gamasd*sarg;
  204.   Vth    = *Von;
  205.   *Vdsat = 0.;
  206.  
  207.   if ( NFS == 0. ||
  208.        Cox    == 0. ) {
  209.     vgst = Vgs - *Von;
  210.     if ( vgst <= 0. ) {        /* cut-off region */
  211.       *Gds = 0.;
  212.       goto L1050;
  213.     } }
  214.   else {
  215.     double cdonco = factor - Gamasd*dsrgdb - dgddvb*sarg;
  216.  
  217. /*    xn = 1 + cdonco + CHARGE*NFS/Cox*Wid*Len;*/
  218. /*    xn    = 1 + cdonco + CHARGE*NFS/COX;*/
  219.     xn    = 1 + cdonco + CHARGE*NFS*TOX/EPSOX;
  220.     temp  = xn*vt;
  221.     *Von += temp;
  222.     vgst  = Vgs - *Von;
  223.     argg  = 1/temp;
  224.     }
  225.  
  226. /* compute some more useful quantities */
  227.  
  228.   body   = barg*barg*barg - ( sarg3 = sarg*sarg*sarg );
  229.   dbdvds = 2*barg*barg*(-dbrgdb);    /* dbrgdb = -dbdvbs */
  230.   gammad = Gamasd;
  231.   dgdvbs = dgddvb;
  232.   gdbdv  = 2*gammad*( barg*barg*dbrgdb - sarg*sarg*dsrgdb );
  233.   dodvbs = -factor + dgdvbs*sarg + gammad*dsrgdb;
  234.   if ( NFS != 0. && Cox != 0. ) {
  235.     dxndvb = 2*dgdvbs*dsrgdb + gammad*d2sdb2 + dgddb2*sarg;
  236.     dodvbs = dxndvb*vt + dodvbs;
  237.     dxndvd = dgdvds*dsrgdb;
  238.     dodvds = dxndvd*vt + dgdvds*sarg;
  239.     }
  240.  
  241. /* evaluate effective mobility and its derivatives */
  242.  
  243.   if ( Cox != 0. && vgst > VBP ) {
  244.     ufact  = pow( VBP/vgst, UEXP);
  245.     ueff   = ufact*UO;
  246.     temp   = ufact*UEXP/vgst;
  247.     dudvgs = -temp;
  248.     dudvbs = temp*dodvbs;
  249.     dudvds = 0.;
  250.     }
  251.   else {
  252.     ufact  = 1.;
  253.     ueff   = UO;
  254.     dudvgs = dudvbs = dudvds = 0.;
  255.     }
  256.  
  257. /* evaluate Vd(sat) & derivatives according to Grove-Frohman equation */
  258.  
  259. /*** re-visit (? re-code) this section ***/
  260.  
  261.   vgsx   = Vgs;
  262.   gammad = Gamasd/eta;
  263.   dgdvbs = dgddvb;
  264.  
  265.   if ( NFS != 0. && Cox != 0. ) vgsx = MAX(Vgs,*Von);
  266.  
  267.   if ( gammad <= 0. ) {
  268.     vpof   = MAX( eta*Vds + Vbin,  0.);
  269.     *Vdsat = MAX( (vgsx-Vbin)/eta, 0.);
  270.     vdsat1 = MAX( (vpof-Vbin)/eta, 0.);
  271.     dsdvgs = 1.;
  272.     dsdvbs = 0.;
  273.     }
  274.   else {
  275.     double
  276.       gammd2 = gammad*gammad,
  277.       argv   = (vgsx-Vbin)/eta + PHI - Vbs;
  278.  
  279.     if ( argv <= 0. ) {
  280.       vpof  = Vth;
  281.       *Vdsat = vdsat1 = dsdvgs = dsdvbs = 0.;
  282.       }
  283.     else {
  284.       double arg = sqrt( 1 + 4*argv/gammd2 );
  285.  
  286.       *Vdsat = (vgsx-Vbin)/eta + .5*gammd2*(1.-arg);
  287.       *Vdsat = MAX(*Vdsat, 0.);
  288.       if ( Charge == YES ) {
  289.         double
  290.           arg1  = gammd2/(eta*eta),
  291.           arg2  = Vds - .5*arg1,
  292.           argsq = (arg2 + .5*arg1 + PHI - Vbs)*arg1, /* <- isn't this
  293.           argsq = (    Vds        + PHI - Vbs)*arg1,    <- THIS ? */
  294.           argv1;
  295.  
  296.         if ( argsq < 0.) vpof = Vth;
  297.         else             vpof = Vbin + eta*( arg2 + .5*arg1 + sqrt(argsq));
  298.  
  299.         argv1 = (vpof-Vbin)/eta + PHI - Vbs;
  300.         if ( argv1 > 0.) {
  301.           arg1   = sqrt( 1 + 4*argv1/gammd2 );
  302.           vdsat1 = (vpof-Vbin)/eta + .5*gammd2*(1-arg1);
  303.           vdsat1 = MAX( vdsat1, 0. );
  304.           }
  305.         else
  306.           vdsat1 = 0.;
  307.         }
  308.       dsdvgs = ( 1- 1/arg )/eta;
  309.       dsdvbs = ( gammad*(1-arg) + 2*argv/(gammad*arg) )/eta*dgdvbs + 1/arg + factor*dsdvgs;
  310.     } }
  311.  
  312. /* store Vd(sat) as above in vpof (pinch-off voltage ) */
  313.  
  314. /* evaluate Vd(sat) & derivatives
  315.    according to Baum's theory of scattering velocity saturation */
  316.  
  317. /*** re-visit (? re-code) this section: solution to quartic eqn.
  318.     ? separate routine ***/
  319.  
  320.   if ( VMAX > 0.) {
  321.     int i, iknt, j, jknt;
  322.     double xvalid, a3, b3, y3, delta4, a4[5], b4[5], x4[9], poly4[9];
  323.     double xv, v2, v1, a1, b1, c1, d1, a, b, c, r, s, s2, p, p0, p2, sarg3lim;
  324.  
  325.       xv = VMAX*Len/ueff;
  326.       v2 = PHI - Vbs;
  327.       /* limit v2, v1, and sarg3 so the following calculations do not exceed
  328.      the vax floating point range (1.7e38) */
  329.       if (fabs(v2) > 5e2) v2 = (v2 > 0. ? 5e2 : -5e2);
  330.       v1 = (vgsx-Vbin)/eta + v2;
  331.       if (fabs(v1) > 5e2) v1 = (v1 > 0. ? 5e2 : -5e2);
  332.  
  333.       a1 = gammad/.75;
  334.       b1 = -2*(v1+xv);
  335.       c1 = -2*gammad*xv;
  336.       if (fabs(sarg3) > 1e8) sarg3lim = (sarg3 > 0. ? 1e8 : -1e8);
  337.       else             sarg3lim = sarg3;
  338.       d1 =  2*v1*(v2+xv) - v2*v2 - a1*sarg3lim;
  339.  
  340.       a  = -b1;
  341.       b  =  a1*c1 - 4*d1;
  342.       c  = -d1*( a1*a1 - 4*b1 ) - c1*c1;
  343.  
  344.       r  = -a*a/3 + b;
  345.       s  = 2*a*a*a/27 - a*b/3 + c;
  346.  
  347.       s2 = s*s;
  348.       p  = s2/4 + r*r*(r/27);
  349.       p0 = fabs(p);
  350.       p2 = sqrt(p0);
  351.  
  352.     if ( p < 0. ) {
  353.       double
  354.         ro = pow( sqrt( s2/4 + p0 ), 1./3.),
  355.         fi = atan(-2*p2/s);
  356.  
  357.       y3 = 2*ro*cos(fi/3) - a/3;
  358.       }
  359.     else {
  360.       double
  361.         p3 = pow( fabs( -s/2 + p2 ), 1./3.),
  362.         p4 = pow( fabs( -s/2 - p2 ), 1./3.);
  363.  
  364.       y3 = p3 + p4 - a/3;
  365.       }
  366.  
  367.     iknt = 0;
  368.     a3 = sqrt( a1*a1/4 - b1 + y3 );
  369.     b3 = sqrt( y3*y3/4 - d1 );
  370.  
  371.     for ( i = 1; i <= 4; i++ ) {    /* start index @ 1 (even in C) */
  372.  
  373.       a4[i]  = a1/2 + sig1[i]*a3;
  374.       b4[i]  = y3/2 + sig2[i]*b3;
  375.  
  376.       delta4 = a4[i]*a4[i]/4 - b4[i];
  377.  
  378.       if ( delta4 < 0.) continue;
  379.  
  380.       x4[++iknt] = -a4[i]/2 + sqrt(delta4);
  381.       x4[++iknt] = -a4[i]/2 - sqrt(delta4);
  382.       }
  383.  
  384.     jknt = 0;
  385.     for ( j = 1; j <= iknt; j++) {
  386.  
  387.       if ( x4[j] <= 0. ) continue;
  388.  
  389.       poly4[j] = POLY4(x4[j],d1,c1,b1,a1,1);
  390.  
  391.       if ( fabs(poly4[j]) > 1E-6 ) continue;
  392.  
  393.       xvalid = ++jknt <= 1 ? x4[j] : MIN( xvalid, x4[j] );
  394.       }
  395.  
  396.     if ( jknt > 0 ) *Vdsat = xvalid*xvalid + Vbs - PHI;
  397.     else            ivmflg++;
  398.     }
  399.  
  400. /* evaluate effective channel length & derivatives */
  401.  
  402.   if ( Vds == 0.) {
  403.  
  404. L610:
  405.     dldvgs =
  406.     dldvds =
  407.     dldvbs = 0.;
  408.  
  409.     }
  410.   else {
  411.     double dldsat;
  412.  
  413.     gammad = Gamasd;
  414.  
  415.     if ( Vbs - *Vdsat > 0. ) {
  416.       bsarg  = sphi/(1+ .5*(Vbs - *Vdsat)/PHI);
  417.       dbsrdb = -.5*bsarg*bsarg/sphi3;
  418.       }
  419.     else {
  420.       bsarg  = sqrt(*Vdsat-Vbs+PHI);
  421.       dbsrdb = -.5/bsarg;
  422.       }
  423.  
  424.     bodys  = bsarg*bsarg*bsarg - sarg3;
  425.     gdbdvs = 2*gammad*( bsarg*bsarg*dbsrdb - sarg*sarg*dsrgdb );
  426.  
  427.     if ( VMAX > 0.) {
  428.       double xls, xlfact,
  429.         xdv    = XD/sqrt(NEFF),
  430.         xlv    = xdv*VMAX/(2*ueff),
  431.         argv   = (vgsx-Vbin)/eta - *Vdsat,
  432.         vqchan = argv - gammad*bsarg,
  433.         dqdsat = gammad*dbsrdb - 1,
  434.         vl     = VMAX*Len,
  435.         dfunds = vl*dqdsat - ueff*vqchan,
  436.         dfundg = (vl - ueff*(*Vdsat))/eta,
  437.         dfundb = -vl*( 1 + dqdsat - factor/eta )
  438.                  + ueff*( gdbdvs - dgdvbs*bodys/1.5 )/eta;
  439.  
  440.       dsdvgs = -dfundg/dfunds;
  441.       dsdvbs = -dfundb/dfunds;
  442.  
  443.       if ( NSUB==0.|| LAMBDA > 0.) goto L610;
  444.  
  445.       argv   = Vds - *Vdsat;
  446.       if (argv >= 0.0) {
  447.         xls    = sqrt( argv + xlv*xlv );
  448.         xlfact = (xdv/Len)/Vds;
  449.         lambda = xlfact*( xls - xlv );
  450.         dldsat = xdv/(2*Len*xls);
  451.         }
  452.       else {
  453.         lambda = dldsat = 0.0;
  454.         }
  455.       }
  456.     else {
  457.       double arg, argv, sargv, xlfact;
  458.  
  459.       if ( NSUB==0.|| LAMBDA > 0.) goto L610;
  460.  
  461.       argv   = (Vds - *Vdsat)/4;
  462.       sargv  = sqrt( 1 + argv*argv );
  463.       arg    = sqrt( argv + sargv );
  464.       xlfact = (XD/Len)/Vds;
  465.       lambda = xlfact*arg;
  466.       dldsat = lambda*Vds/(8*sargv);
  467.       }
  468.  
  469.     dldvgs = dldsat*dsdvgs;
  470.     dldvds = dldsat - lambda;
  471.     dldvbs = dldsat*dsdvbs;
  472.     }
  473.  
  474. /* limit channel shortening at punch-through */
  475.  
  476.   {
  477. /*  xwb    = XD*sbiarg;        moved from above */
  478.  
  479.     double xleff,
  480.       xwb    = XD*sqrt(PB),
  481.       xld    = Len-xwb;
  482.  
  483.     clfact = 1 - lambda*Vds;
  484.     xleff  = clfact*Len;
  485.     dldvds = -lambda-dldvds;
  486.  
  487.     if ( NSUB == 0.) xwb = .25E-6;
  488.  
  489.     if ( xleff < xwb ) {
  490.       double dfact,
  491.         deltal = lambda*Vds*Len;
  492.  
  493.       xleff  = xwb/( 1 + (deltal-xld)/xwb );
  494.       clfact = xleff/Len;
  495.  
  496.       dfact   = xleff*xleff/(xwb*xwb);
  497.       dldvgs *= dfact;
  498.       dldvds *= dfact;
  499.       dldvbs *= dfact;
  500.     } }
  501.  
  502. /* evaluate effective beta (effective k') */
  503.  
  504.   beta1 = (ufact/clfact)*KP*Wid/Len;
  505.  
  506. /* test for mode of operation and branch appropriately */
  507.  
  508.   gammad = Gamasd;
  509.   dgdvbs = dgddvb;
  510.  
  511. /*if (vds.gt.1.0d-8) goto 730 */
  512.   if ( Vds > 0. ) {            /* forward operation */
  513.  
  514.     if ( Vgs <= *Von ) {        /* sub-threshold region */
  515.  
  516.       if ( *Vdsat <= 0. ) {
  517.         *Gds = 0.;
  518.         if ( Vgs <= Vth ) goto L1050;
  519.         *Id = *Gm = *Gmbs = 0.;
  520.         }
  521.       else {
  522.         double cdson, gdson, gbson, didvds, gmw, argt, expg, arg;
  523.         double vdson = MIN(*Vdsat,Vds);
  524.  
  525.         if ( Vds > *Vdsat ) {
  526.           barg   = bsarg;
  527.           dbrgdb = dbsrdb;
  528.           body   = bodys;
  529.           gdbdv  = gdbdvs;
  530.           }
  531.         cdson  = beta1*( (*Von - Vbin - .5*eta*vdson)*vdson - gammad*body/1.5);
  532.         didvds = beta1*(*Von - Vbin - eta*vdson - gammad*barg);
  533.  
  534.         gdson  = -cdson*dldvds/clfact - beta1*dgdvds*body/1.5;
  535.         if ( Vds < *Vdsat ) gdson = gdson + didvds;
  536.  
  537.         gbson  = -cdson*dldvbs/clfact + beta1*( dodvbs*vdson + factor*vdson - dgdvbs*body/1.5 - gdbdv );
  538.         if ( Vds > *Vdsat ) gbson = gbson + didvds*dsdvbs;
  539.  
  540.         argt = argg*(Vgs - *Von);
  541.         expg = EXP(argg*(Vgs - *Von));
  542.     /* check expg value to avoid underflow or overflow in later use */
  543.     if (expg < MINREAL) expg = 0.0;
  544.         arg = cdson*( dudvgs/ufact - dldvgs/clfact );
  545.         *Gm = arg + beta1*(*Vdsat);
  546.         if (Vds > *Vdsat)
  547.           *Gm += beta1*(*Von - Vbin - eta*(*Vdsat) - gammad*bsarg)*dsdvgs;
  548.  
  549.         *Id = cdson*expg;
  550.         gmw = *Id*argg;
  551.         *Gm = gmw;
  552.         if ( Vds > *Vdsat )
  553.           *Gm = gmw + didvds*dsdvgs*expg;
  554.  
  555.         *Gds  = gdson*expg - *Gm*dodvds - gmw*(Vgs - *Von)*dxndvd/xn;
  556.         *Gmbs = gbson*expg - *Gm*dodvbs - gmw*(Vgs - *Von)*dxndvb/xn;
  557.       } }
  558.     else if ( Vds > *Vdsat ) {        /* saturation region */
  559.       double arg;
  560.  
  561.       *Id   = beta1*( ( Vgs - Vbin - eta*(*Vdsat)/2 )*(*Vdsat) - gammad*bodys/1.5 );
  562.       arg   = *Id*( dudvgs/ufact - dldvgs/clfact );
  563.       *Gm   = arg + beta1*(*Vdsat) + beta1*( Vgs - Vbin - eta*(*Vdsat) - gammad*bsarg)*dsdvgs;
  564.       *Gds  = -*Id*dldvds/clfact - beta1*dgdvds*bodys/1.5;
  565.       arg   = *Id*( dudvbs/ufact - dldvbs/clfact );
  566.       *Gmbs = arg - beta1*( gdbdvs + dgdvbs*bodys/1.5 - factor*(*Vdsat) ) + beta1*( Vgs - Vbin - eta*(*Vdsat) - gammad*bsarg )*dsdvbs;
  567.       }
  568.     else {                /* linear region */
  569.       double arg;
  570.  
  571.       *Id   = beta1*( ( Vgs - Vbin - eta*Vds/2 )*Vds - gammad*body/1.5 );
  572.       arg   = *Id*( dudvgs/ufact - dldvgs/clfact );
  573.       *Gm   = arg + beta1*Vds;
  574.       arg   = *Id*( dudvds/ufact - dldvds/clfact );
  575.       *Gds  = arg + beta1*( Vgs - Vbin - eta*Vds - gammad*dbdvds - dgdvds*body/1.5 );
  576.       arg   = *Id*( dudvbs/ufact - dldvbs/clfact );
  577.       *Gmbs = arg - beta1*( gdbdv + dgdvbs*body/1.5 - factor*Vds );
  578.       }
  579.  
  580. /* compute charges for "on" region */
  581.  
  582.     if ( Charge == YES ) {
  583.       if ( Vgs > Vth ) {
  584.         Mqspof(Model,Vds,Vbs,Vgs,vpof,*Von,*Vdsat,vdsat1,Cggb,Cgdb,Cgsb,Cbgb,Cbdb,Cbsb,Qgate,Qchan,Qbulk);
  585.         }
  586.       else {
  587.         MOS2cap(Model,Vds,Vbs,Vgs,*Vdsat,Cggb,Cgdb,Cgsb,Cbgb,Cbdb,Cbsb,Qgate,Qchan,Qbulk);
  588.         Qspof = 0.;
  589.       } }
  590.     else
  591.       *Qgate = *Qbulk = *Qchan = Qspof = 0.;
  592.     }
  593.   else {                /* reverse operation */
  594.     if ( Vgs > *Von )
  595.       *Gds = beta1*( Vgs - Vbin - gammad*sarg );
  596.  
  597.     else if ( NFS == 0. || Cox == 0. )
  598.       *Gds = 0.;
  599.  
  600.     else {
  601.       double argt = argg*(Vgs - *Von);
  602.  
  603.       argt = EXP(argt);
  604.  
  605.       *Gds = beta1*(*Von - Vbin - gammad*sarg)*argt;
  606.       }
  607. L1050:
  608.     *Id = *Gm = *Gmbs = 0.;
  609.     if ( Charge == YES ) {
  610.       MOS2cap(Model,Vds,Vbs,Vgs,*Vdsat,Cggb,Cgdb,Cgsb,Cbgb,Cbdb,Cbsb,Qgate,Qchan,Qbulk);
  611.       Qspof = 0.;
  612.       }
  613.     else
  614.       *Qgate = *Qbulk = *Qchan = Qspof = 0.;
  615.     }
  616.  
  617.   } /* end MOS2eval */
  618.  
  619. #ifndef cc_noline
  620. #line 4 "MOS2cap"
  621. #endif
  622.  
  623. void MOS2cap(Model, Vds, Vbs, Vgs, Vdsat,                 /*TITLE*/
  624.          Cggb, Cgdb, Cgsb,
  625.          Cbgb, Cbdb, Cbsb,
  626.          Qgate, Qchan, Qbulk)
  627.  
  628. struct M_  *Model;
  629. double     Vds, Vbs, Vgs, Vdsat,
  630.            *Cggb, *Cgdb, *Cgsb,
  631.            *Cbgb, *Cbdb, *Cbsb,
  632.            *Qgate, *Qchan, *Qbulk;
  633.  
  634. /*****************************************************************************
  635. * Purpose
  636. *   calculate MOS level 2 capacitance and charge
  637. *
  638. * Returned value
  639. *   none
  640. *
  641. * Discussion
  642. *   none
  643. *
  644. * Author
  645. *   PWT - 29 Aug 86 - creation                        
  646. *   PWT   10 Oct 87 - re-work as part of BSIM addition
  647. *
  648. *****************************************************************************/
  649. {
  650. /* Local variables */
  651.   double vgb, vg;
  652.  
  653. /* Code */
  654. /* initialize charges */
  655.  
  656.   *Qgate = *Qbulk = 0.;
  657.  
  658. /* reference voltages for charge computation */
  659.  
  660.   vgb = Vgs - Vbs;
  661.   vg  = vgb - Vbin + PHI;
  662.  
  663. /* determine operating region */
  664.  
  665.   if ( Vgs > Vth ) {            /* "on" region */
  666.     double
  667.       vbd, vd, vsat,
  668.       vs, vs2, vsp5, vs1p5,
  669.       ve, ve2, vep5, ve1p5,
  670.       term0, term1, term2, term3, term4, term5, term6, term7,
  671.       term10, term11, term12,
  672.       term20, term21, term22,
  673.       argn, argd,
  674.       dvedvd, 
  675.       dgndve, dddve,
  676.       dgndvs, dddvs;
  677.  
  678.     vbd = Vbs - Vds;
  679.     vd  = MAX( PHI - vbd, 1e-8);
  680.  
  681.     vs    = MAX( PHI - Vbs, 1e-8);
  682.     vs2   = vs*vs;
  683. /*  vs3   = vs2*vs;    never used */
  684. /*  vs5   = vs3*vs2;    never used */
  685.     vsp5  = sqrt(vs);
  686.     vs1p5 = vsp5*vs;
  687. /*  vs2p5 = vs1p5*vs;    never used */
  688.  
  689.     vsat = Vdsat + vs;
  690.  
  691.     if ( vd < vsat ) {
  692.       ve     = vd;
  693.       dvedvd = 1.;
  694. /*    dvedvg = 0.;    never used */
  695.       }
  696.     else {
  697.       ve     = vsat;
  698.       dvedvd = 0.;
  699. /*    dvedvg = 0.;    never used */
  700.       }
  701.     ve2   = ve*ve;
  702. /*  ve3   = ve2*ve;    never used */
  703. /*  ve5   = ve2*ve3;    never used */
  704.     vep5  = sqrt(ve);
  705.     ve1p5 = vep5*ve;
  706. /*  ve2p5 = ve1p5*ve;    never used */
  707.  
  708.     term0  = ve + vs;
  709.     term1  = vep5 + vsp5;
  710.     term2  = vep5*vsp5;
  711.     term3  = ve2 + vs2;
  712.     term4  = ve*vs;
  713.     term5  = term0*term1;
  714.     term6  = (term3 + term4) + term2*term0;
  715.     term7  = (term3 + term4)*term1;
  716.     term10 = vep5 + .5*vsp5;
  717.     term11 = 1.5*ve + vsp5*term10;
  718.     term12 = 2*ve1p5 + vsp5*term11;
  719.     term20 = .5*vep5 + vsp5;
  720.     term21 = 1.5*vs + vep5*term20;
  721.     term22 = 2*vs1p5 + vep5*term21;
  722.     argn   = .5*vg*term5 - .4*Gamasd*term6 - term7/3;
  723.     argd   = vg*term1 - Gamasd*( term0 + term2 )/1.5 - .5*term1*term0;
  724. /*  argd2  = argd*argd;
  725.  *** never used */
  726.  
  727.     *Qgate = Cox*( vg - argn/argd );
  728.     dgndve = .5*vg*term11 - .4*Gamasd*term12 - (2.5*ve2 + vsp5*term12)/3;
  729.     dddve  = .5*vg - Gamasd*term10/1.5 - .5*term11;
  730. /*  dqgdve = -Cox/argd*( dgndve - ( vg - *Qgate/Cox )*dddve );
  731.  *** never used */
  732.     dgndvs = .5*vg*term21 - .4*Gamasd*term22 - ( 2.5*vs2 + vep5*term22 )/3;
  733.     dddvs  = .5*vg - Gamasd*term20/1.5 - .5*term21;
  734.  
  735.     *Cgdb  = -Cox/(argd*vep5)*( dgndve - ( vg - *Qgate/Cox )*dddve )*dvedvd;
  736.     *Cgsb  = -Cox/(argd*vsp5)*( dgndvs - ( vg - *Qgate/Cox )*dddvs );
  737.     *Cggb  = Cox*( 1 - term1/argd*( .5*term0 - vg + *Qgate/Cox ) );
  738.  
  739.     argn   = vg*( term0 + term2 )/1.5 - .5*Gamasd*term5 - .4*term6;
  740.     dgndve = vg*term10/1.5 - .5*Gamasd*term11 - .4*term12;
  741.     dgndvs = vg*term20/1.5 - .5*Gamasd*term21 - .4*term22;
  742.  
  743.     *Qbulk = -Gamasd*Cox*argn/argd;
  744.  
  745.     *Cbdb  = -Cox/(vep5*argd)*( *Qbulk/Cox*dddve + Gamasd*dgndve )*dvedvd;
  746.     *Cbsb  = -Cox/(vsp5*argd)*( *Qbulk/Cox*dddvs + Gamasd*dgndvs );
  747.     *Cbgb  = -Cox/argd*( Gamasd*( term0 + term2 )/1.5 + *Qbulk/Cox*term1 );
  748.     }
  749.   else {                /* "off" region */
  750.  
  751.     if ( vg > 0. ) {
  752.       double
  753.         gamma2 = .5*Gamasd,
  754.         sqarg  = sqrt( gamma2*gamma2 + vg );
  755.  
  756.       *Qgate = Gamasd*Cox*( sqarg - gamma2 );
  757.       *Cggb  = .5*Cox*Gamasd/sqarg;
  758.       }
  759.     else {
  760.       *Qgate = Cox*vg;
  761.       *Cggb  = Cox;
  762.       }
  763.     *Qbulk = -*Qgate;
  764.     *Cbgb  = -*Cggb;
  765.     *Cgdb  =  *Cgsb = *Cbdb = *Cbsb = 0.;
  766.     }
  767.  
  768.   *Qchan = -( *Qgate + *Qbulk );
  769.  
  770.   } /* end MOS2cap */
  771.