home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 13 / 13.iso / p / p064 / 3.ddi / GASFET.C < prev    next >
Encoding:
C/C++ Source or Header  |  1991-07-01  |  17.1 KB  |  573 lines
  1. /****************************************************************************/
  2. /*        Copyright (C) 1986, 1987, 1988, 1989, 1990, 1991                  */
  3. /*          By MicroSim Corporation, All Rights Reserved                    */
  4. /****************************************************************************/
  5. /* gasfet.c
  6.  *   $Revision:   1.18  $
  7.  *   $Author:   pwt  $
  8.  *   $Date:   07 Nov 1990 10:35:18  $ */
  9.  
  10. /******************* USERS OF DEVICE EQUATIONS OPTION ***********************/
  11. /******** The procedure for changing the GaAs FET 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 B_DEVICE
  16.  
  17. #include "option1.h"
  18.  
  19. #ifdef USE_DECL_ARGS
  20. int GaAsfet(struct b_ *,int ,int ,int );
  21. static double
  22.   GaAs2cap(struct B_ *,double,double,double,double,double,double *,double *);
  23. #else
  24. int GaAsfet();
  25. static double GaAs2cap();
  26. #endif
  27.  
  28. int GaAsfet(    /* Process GaAsfet for DC and TRANSIENT analyses */
  29.   Instance,    /* a device to evaluate */
  30.   ModeFl,    /* mode control flag:   see "tran.h" */
  31.   InitFl,    /* initialization flag: see "tran.h" */
  32.   LoadFl    /* NO: bypass MatPrm load if solution unchanged */
  33.   )        /* return: did not converge (YES/NO) */
  34.   struct b_ *Instance;
  35.   int ModeFl, InitFl, LoadFl;
  36.  
  37. /*
  38. pwt    15 Aug 86    creation
  39. pwt    ?? Jul 87    add Raytheon model (level==2)
  40. pwt    02 Nov 87    correct usage of TAU (transit time)
  41. pwt    16 Mar 88    add shared data area to J.H, change source to match
  42. pwt    07 Mar 89    correct sign on reverse gm for Raytheon model
  43. whjb    06 Apr 89    LoadFl==YES forces full recalculation
  44. pwt    15 Apr 89    set T=0 cap. current
  45. sv    12 Feb 90    add local temperatures to models.
  46. pwt    12 Apr 90    add TriQuint model (level==3)
  47. pwt    04 May 90    minor speed-up in some calculations
  48. pwt    08 Aug 90    correct IC= for device type
  49. pwt    07 Nov 90    add VMAX and VDELTA for TriQuint
  50. */
  51. { struct B_ *model;    /* device model */
  52.   double
  53.     beta,        /* device transconductance */
  54.     isat,        /* junction saturation current */
  55.     vds, vgs, vgd,    /* terminal voltages */
  56.     vte,        /* effective thermal voltage */
  57.     gm,            /* forward transconductance */
  58.     gds, ggs, ggd,    /* conductances */
  59.     ig, id, igd,    /* currents */
  60.     idrain,        /* drain current */
  61.     ig_hat, id_hat,    /* predicted currents */
  62.     cgs, cgd, cds,    /* device capacitances */
  63.     type;        /* +1 (so far) */
  64.  
  65. #define vt    VT
  66.  
  67.   int
  68.     level,
  69.     jlim_fl = YES,    /* junction limiting flag: YES, NO */
  70.     pred_fl = NO,    /* pred. calculated  flag: YES, NO */
  71.     nonconv = NO;    /* non-convergence   flag: YES, NO */
  72.  
  73. /* defines & declarations for state-vector access */
  74.  
  75.   struct bsv_def *sv[MSTVCT];
  76.  
  77.   EVAL_SV(sv,Instance,b_sda.b_sv);
  78.  
  79. /* defines for model parameter access: simplify code appearance */
  80.  
  81. #define PARAM(a) ((double)model->a)
  82. #define AF    PARAM(B_af)
  83. #define ALPHA    PARAM(B_alpha)
  84. #define B    PARAM(B_b)
  85. #define BETA    PARAM(B_beta)
  86. #define CDS    PARAM(B_cds)
  87. #define CGD    PARAM(B_cgd)
  88. #define CGS    PARAM(B_cgs)
  89. #define DELT    PARAM(B_delta)    /* TriQuint "delta" */
  90. #define FC    PARAM(B_fc)
  91. #define FCPB    FC        /* = fc*vbi */
  92. #define F1    PARAM(B_f1)
  93. #define F2    PARAM(B_f2)
  94. #define F3    PARAM(B_f3)
  95. #define GAMMA    PARAM(B_gamma)
  96. #define IS    PARAM(B_is)
  97. #define KF    PARAM(B_kf)
  98. #define LAMBDA    PARAM(B_lambda)
  99. #define LEVEL    PARAM(B_level)
  100. #define M    PARAM(B_m)
  101. #define N    PARAM(B_n)
  102. #define Q    PARAM(B_q)
  103. #define RD    PARAM(B_rd)
  104. #define RG    PARAM(B_rg)
  105. #define RS    PARAM(B_rs)
  106. #define TAU    PARAM(B_tau)
  107. #define VBI    PARAM(B_vbi)
  108. #define VCRIT    PARAM(B_vcrit)
  109. #define VDELTA    PARAM(B_vdelta)    /* Statz "delta" */
  110. #define VMAX    PARAM(B_vmax)
  111. #define VTO    PARAM(B_vto)
  112.  
  113. #define IG0    (Instance->bcv_ig)
  114. #define ID0    (Instance->bcv_id)
  115. #define IGD0    (Instance->bcv_igd)
  116. #define GM0    (Instance->bcv_gm)
  117. #define GDS0    (Instance->bcv_gds)
  118. #define GGS0    (Instance->bcv_ggs)
  119. #define GGD0    (Instance->bcv_ggd)
  120.  
  121. #define CGS0    (Instance->b_sda.b_ac.bac_cgs)
  122. #define CGD0    (Instance->b_sda.b_ac.bac_cgd)
  123. #define CDS0    (Instance->b_sda.b_ac.bac_cds)
  124.  
  125. #define DevOff    (Instance->b_off)
  126. #define AREA    (Instance->b_area)
  127.  
  128.   model = Instance->b_model;        /* find the model */
  129.  
  130.   type  = model->B_type > 0 ? 1.: -1.;
  131.  
  132. /* initialization */
  133.  
  134.   if      ( InitFl==INSTV0 ) {        /* use prev. iteration */
  135.     vgs = ( VOLTAGE(b_g) - VOLTAGE(b_s) )*type;
  136.     vgd = ( VOLTAGE(b_g) - VOLTAGE(b_d) )*type;
  137.     }
  138.   else if ( InitFl==INTRAN ) {        /* use prev. time-step */
  139.     vgs = B_VGS(1);
  140.     vgd = B_VGD(1);
  141.     }
  142.   else if ( InitFl==INOFF &&        /* set "off" devices */
  143.             DevOff==YES )
  144.     vgs = vgd = 0.;
  145.  
  146.   else if ( InitFl==ININIT ) {        /* use IC= values */
  147.  
  148.     if ( ModeFl==MDBPTR && NOSOLV==YES ) {
  149.       vgs = (    Instance->b_vgs)*type;
  150.       vgd = (vgs-Instance->b_vds)*type;
  151.       }
  152.     else
  153.       vgd = vgs = ( DevOff ? 0.: -1.);
  154.     }
  155.   else {
  156.     double del_vgs, del_vgd, del_vds, del_id, del_ig;
  157.  
  158.     if ( InitFl==INPRDCT ) {        /* extrapolate value */
  159.       double arg = DELTA/DELOLD[1];
  160.  
  161.       vgs = arg*( B_VGS(1) - B_VGS(2) ) + B_VGS(1);
  162.       vgd = arg*( B_VGD(1) - B_VGD(2) ) + B_VGD(1);
  163.       *sv[0] = *sv[1];
  164.       }
  165.     else {                /* use current value */
  166.       vgs = ( VOLTAGE(b_g) - VOLTAGE(b_s) )*type;
  167.       vgd = ( VOLTAGE(b_g) - VOLTAGE(b_d) )*type;
  168.       }
  169.  
  170. /* compute new non-linear branch voltage */
  171.  
  172.     del_vgd = vgd - B_VGD(0);
  173.     del_vgs = vgs - B_VGS(0);
  174.     del_vds = del_vgs - del_vgd;
  175.  
  176.     del_id = del_vds*GDS0 - del_vgd*GGD0 +
  177.              ( B_VGS(0) > B_VGD(0) ? del_vgs : del_vgd )*GM0;
  178.     id_hat = del_id + ID0;
  179.  
  180.     del_ig = del_vgs*GGS0 + del_vgd*GGD0;
  181.     ig_hat = del_ig + IG0;
  182.  
  183. /* bypass if solution not changed */
  184.  
  185.     if ( InitFl==INPRDCT && DELTA!=DELOLD[1] || LoadFl==YES ) ;
  186.     else if ( fabs(del_id)  >= TOL(id_hat,     ID0,RELTOL,ABSTOL) ) ;
  187.     else if ( fabs(del_vgd) >= TOL(   vgd,B_VGD(0),RELTOL, VNTOL) ) ;
  188.     else if ( fabs(del_vgs) >= TOL(   vgs,B_VGS(0),RELTOL, VNTOL) ) ;
  189.     else if ( fabs(del_ig)  >= TOL(ig_hat,     IG0,RELTOL,ABSTOL) ) ;
  190.     else {
  191.       goto done;
  192.       }
  193.     pred_fl = YES;
  194.  
  195. /* solution changed: limit non-linear branch voltages */
  196.  
  197.     jlim_fl  = PNJLIM(vgs,B_VGS(0),vt,VCRIT);
  198.     jlim_fl |= PNJLIM(vgd,B_VGD(0),vt,VCRIT);
  199.  
  200.     FETlim(&vgs,B_VGS(0),VTO);
  201.     FETlim(&vgd,B_VGD(0),VTO);
  202.  
  203.     } /* end of initialization */
  204.  
  205. /* compute DC current and derivatives */
  206.  
  207.   level = NINT(LEVEL);
  208.   beta  = AREA*BETA;
  209.   isat  = AREA*IS;
  210.   vte   = N*vt;
  211.   vds   = vgs-vgd;
  212.  
  213.   if ( vgs > -10*vte ) {        /* currents re. vgs */
  214.     double evgs = EXP(vgs/vte);
  215.  
  216.     ggs = (evgs/vte)*isat + GMIN;
  217.     ig  = (evgs - 1)*isat + vgs*GMIN;
  218.     }
  219.   else {
  220.     ggs = GMIN;
  221.     ig  = ggs*vgs - isat;
  222.     }
  223.  
  224.   if ( vgd > -10*vte ) {        /* currents re. vgd */
  225.     double evgd = EXP(vgd/vte);
  226.  
  227.     ggd = (evgd/vte)*isat + GMIN;
  228.     igd = (evgd - 1)*isat + vgd*GMIN;
  229.     }
  230.   else {
  231.     ggd = GMIN;
  232.     igd = ggd*vgd - isat;
  233.     }
  234.  
  235.   ig += igd;
  236.  
  237. /* compute drain current & derivatives */
  238.  
  239.   if ( vds >= 0.) {            /* forward mode */
  240.     double vgst = vgs-VTO;
  241.  
  242.     if ( level==3 ) vgst += GAMMA*vds;    /* TriQuit correction */
  243.  
  244.     if ( vgst <= 0.)            /* cut-off */
  245.       idrain = gm = gds = 0.;
  246.     else {
  247.       double
  248.         vgstsq = vgst*vgst,
  249.         prod   = 1 + vds*LAMBDA,
  250.         betap  = prod*beta;
  251.  
  252.       if ( level==1 ) {            /* Curtice: linear & saturated */
  253.         double arg = tanh(vds*ALPHA);
  254.  
  255.         idrain = betap*vgstsq*arg;
  256.         gm     = 2*betap*vgst*arg;
  257.         gds    = beta*vgstsq*arg*LAMBDA + betap*vgstsq*(1-arg*arg)*ALPHA;
  258.         }
  259.       else if ( level==2 ) {        /* Raytheon model */
  260.         double arg = 1 + vgst*B;
  261.  
  262.         if ( vds < 3/ALPHA ) {        /* Raytheon: linear */
  263.           double
  264.             afact = 1 - vds*ALPHA/3,
  265.             lfact = 1 - afact*afact*afact;
  266.  
  267.           idrain = betap*vgstsq*lfact/arg;
  268.           gm     = betap*vgst*(1+arg)*lfact/(arg*arg);
  269.           gds    = beta*vgstsq*(afact*afact*prod*ALPHA + lfact*LAMBDA)/arg;
  270.           }
  271.         else {                /* Raytheon: saturated */
  272.           idrain = betap*vgstsq/arg;
  273.           gm     = betap*vgst*(1+arg)/(arg*arg);
  274.           gds    = beta*vgstsq*LAMBDA/arg;
  275.         } }
  276.       else {                /* TriQuint model */
  277.  
  278.         if ( vds < 3/ALPHA ) {        /* TriQuint: linear */
  279.           double
  280.             afact   = 1 - vds*ALPHA/3,
  281.             lfact   = 1 - afact*afact*afact,
  282.             lnvgst  = log(vgst),
  283.             izero   = beta*exp(Q*lnvgst)*lfact,
  284.             gmzero  = Q*beta*exp((Q-1)*lnvgst)*lfact,
  285.             gdszero = beta*exp((Q-1)*lnvgst)*(lfact*Q*GAMMA
  286.                       + ALPHA*afact*afact*vgst),
  287.             pfact   = 1/(1+DELT*vds*izero);
  288.  
  289.           idrain = izero*pfact;
  290.           gm     = gmzero*pfact*pfact;
  291.           gds    = gdszero*pfact*pfact - DELT*idrain*idrain;
  292.           }
  293.         else {                /* TriQuint: saturated */
  294.           double
  295.             lnvgst  = log(vgst),
  296.             izero   = beta*exp(Q*lnvgst),
  297.             gmzero  = Q*beta*exp((Q-1)*lnvgst),
  298.             gdszero = Q*GAMMA*beta*exp((Q-1)*lnvgst),
  299.             pfact   = 1/(1+DELT*vds*izero);
  300.  
  301.           idrain = izero*pfact;
  302.           gm     = gmzero*pfact*pfact;
  303.           gds    = gdszero*pfact*pfact - DELT*idrain*idrain;
  304.         } }
  305.     } }
  306.   else {                /* reverse mode */
  307.     double vgdt = vgd-VTO;
  308.  
  309.     if ( level==3 ) vgdt -= GAMMA*vds;    /* TriQuit correction */
  310.  
  311.     if ( vgdt <= 0. )            /* cut-off */
  312.       idrain = gm = gds = 0.;
  313.     else {
  314.       double
  315.         vgdtsq = vgdt*vgdt,
  316.         prod   = 1 - vds*LAMBDA,
  317.         betap  = prod*beta;
  318.  
  319.       if ( level==1 ) {            /* Curtice: linear & saturated */
  320.         double arg = tanh(-vds*ALPHA);
  321.  
  322.         idrain = -betap*vgdtsq*arg;
  323.         gm     = 2*betap*vgdt*arg;
  324.         gds    = beta*vgdtsq*arg*LAMBDA + betap*vgdtsq*(1-arg*arg)*ALPHA;
  325.         }
  326.       else if ( level==2 ) {        /* Raytheon model */
  327.         double arg = 1 + vgdt*B;
  328.  
  329.         if ( -vds < 3/ALPHA ) {        /* Raytheon: linear */
  330.           double
  331.             afact = 1 + vds*ALPHA/3,
  332.             lfact = 1 - afact*afact*afact;
  333.  
  334.           idrain = -betap*vgdtsq*lfact/arg;
  335.           gm     = betap*vgdt*(1+arg)*lfact/(arg*arg);
  336.           gds    = beta*vgdtsq*(afact*afact*prod*ALPHA + lfact*LAMBDA)/arg;
  337.           }
  338.         else {                /* Raytheon: saturated */
  339.           idrain = -betap*vgdtsq/arg;
  340.           gm     = betap*vgdt*(1+arg)/(arg*arg);
  341.           gds    = beta*vgdtsq*LAMBDA/arg;
  342.         } }
  343.       else {                /* TriQuint model */
  344.  
  345.         if ( -vds < 3/ALPHA ) {        /* TriQuint: linear */
  346.           double
  347.             afact   = 1 + vds*ALPHA/3,
  348.             lfact   = 1 - afact*afact*afact,
  349.             lnvgdt  = log(vgdt),
  350.             izero   = beta*exp(Q*lnvgdt)*lfact,
  351.             gmzero  = Q*beta*exp((Q-1)*lnvgdt)*lfact,
  352.             gdszero = beta*exp((Q-1)*lnvgdt)*(lfact*Q*GAMMA
  353.                       + ALPHA*afact*afact*vgdt),
  354.             pfact   = 1/(1-DELT*vds*izero);       /* power correction */
  355.  
  356.           idrain = -izero*pfact;
  357.           gm     = gmzero*pfact*pfact;
  358.           gds    = gdszero*pfact*pfact - DELT*idrain*idrain;
  359.           }
  360.         else {                /* TriQuint: saturated */
  361.           double
  362.             lnvgdt  = log(vgdt),
  363.             izero   = beta*exp(Q*lnvgdt),
  364.             gmzero  = Q*beta*exp((Q-1)*lnvgdt),
  365.             gdszero = Q*GAMMA*beta*exp((Q-1)*lnvgdt),
  366.             pfact   = 1/(1-DELT*vds*izero);
  367.  
  368.           idrain = -izero*pfact;
  369.           gm     = gmzero*pfact*pfact;
  370.           gds    = gdszero*pfact*pfact - DELT*idrain*idrain;
  371.         } }
  372.     } }
  373.  
  374.   if ( ModeFl==MDTRAN && TAU!=0.) {    /* Weil's approximation for */
  375.     double                /* time delay applied with */
  376.       arg1  = DELTA/TAU,        /* backward Euler integration */
  377.       arg2  = 3*arg1,
  378.       denom = ( arg1 *= arg2 ) + arg2 + 1,
  379.       arg3  = arg1/denom;
  380.  
  381.     if ( InitFl==INTRAN )        /* set-up initial current */
  382.       B_IEXDS(2) = B_IEXDS(1) = idrain;
  383.  
  384.     idrain *= arg3;
  385.     gm     *= arg3;
  386.  
  387.     id = ( B_IEXDS(1)*( 1 + DELTA/DELOLD[1] + arg2 ) -
  388.            B_IEXDS(2)*DELTA/DELOLD[1] )/denom;
  389.  
  390.     B_IEXDS(0) = id + idrain;
  391.     }
  392.   else
  393.     id = 0.;
  394.  
  395.   id += idrain-igd;            /* equiv. drain current */
  396.  
  397. /* charge calculations */
  398.  
  399.   if ( ModeFl==MDTRAN ||
  400.        InitFl==INSTV0 ||
  401.      ( ModeFl==MDBPTR && NOSOLV==YES ) ) {
  402.     double
  403.       cgso = CGS*AREA,        /* scale 0-bias cap. */
  404.       cgdo = CGD*AREA,
  405.       cdso = CDS*AREA;
  406.  
  407.     if ( level==1 ) {
  408.       PNJCHG(B_QCGS(0),cgs,cgso,vgs,VBI,FCPB,M,F1,F2,F3);
  409.       PNJCHG(B_QCGD(0),cgd,cgdo,vgd,VBI,FCPB,M,F1,F2,F3);
  410.       }
  411.     else {
  412.       double
  413.         vcap = 1/ALPHA,
  414.         qgga = GaAs2cap(model,     vgs,     vgd,vcap,cgso,cgdo,&cgs,&cgd),
  415.         qggb = GaAs2cap(model,B_VGS(1),     vgd,vcap,cgso,cgdo,NULL,NULL),
  416.         qggc = GaAs2cap(model,     vgs,B_VGD(1),vcap,cgso,cgdo,NULL,NULL),
  417.         qggd = GaAs2cap(model,B_VGS(1),B_VGD(1),vcap,cgso,cgdo,NULL,NULL);
  418.  
  419.       if ( ModeFl==INTRAN )    /* set-up for initial charge equation */
  420.         B_QCGS(1) = B_QCGD(1) = qgga;
  421.  
  422.       B_QCGS(0) = .5*(qgga-qggb + qggc-qggd) + B_QCGS(1);
  423.       B_QCGD(0) = .5*(qgga-qggc + qggb-qggd) + B_QCGD(1);
  424.       }
  425.  
  426.     B_QCDS(0) = cdso*vds;
  427.     cds = cdso;
  428.  
  429.     B_ICGS(0) =        /* reset, so T=0 results (e.g. Probe) are correct */
  430.     B_ICGD(0) =
  431.     B_ICDS(0) = 0.;
  432.  
  433.     if ( InitFl==INSTV0 ) {        /* store small-signal parameters */
  434.       CGS0 = cgs;            /* and update the state vector */
  435.       CGD0 = cgd;
  436.       CDS0 = cds;
  437.  
  438.       goto load;
  439.       }
  440.  
  441.     if ( !(ModeFl==MDBPTR && NOSOLV==NO) ) {    /* transient analysis */
  442.       double g_eq, i_eq;
  443.  
  444.       if ( InitFl==INTRAN ) {            /* for 1st transient step: */
  445.         B_QCGS(2) = B_QCGS(1) = B_QCGS(0);    /* set "old" charges so */
  446.         B_QCGD(2) = B_QCGD(1) = B_QCGD(0);    /* integration works */
  447.         B_QCDS(2) = B_QCDS(1) = B_QCDS(0);
  448.         }
  449.  
  450.       INTEGR8(g_eq,i_eq,cgs,B_QIGS(0),B_QIGS(1),B_QIGS(2));
  451.       ggs  += g_eq;
  452.  
  453.       INTEGR8(g_eq,i_eq,cgd,B_QIGD(0),B_QIGD(1),B_QIGD(2));
  454.       ggd  += g_eq;
  455.  
  456.       INTEGR8(g_eq,i_eq,cds,B_QIDS(0),B_QIDS(1),B_QIDS(2));
  457.       gds  += g_eq;
  458.  
  459.       ig  += B_ICGD(0) + B_ICGS(0);
  460.       igd += B_ICGD(0);
  461.       id  += B_ICDS(0) - B_ICGD(0);
  462.  
  463.       if ( InitFl==INTRAN ) {            /* for 1st transient step: */
  464.         B_ICGS(1) = B_ICGS(0);            /* set "old" currents so */
  465.         B_ICGD(1) = B_ICGD(0);            /* time-step calc. works */
  466.         B_ICDS(1) = B_ICDS(0);
  467.         }
  468.  
  469.     } } /* end of charge calculations */
  470.  
  471. /* check convergence */
  472.  
  473.   if ( pred_fl && (
  474.          jlim_fl ||
  475.          fabs(id_hat-id) > TOL(id_hat,id,RELTOL,ABSTOL) ||
  476.          fabs(ig_hat-ig) > TOL(ig_hat,ig,RELTOL,ABSTOL)
  477.          )
  478.        ) nonconv = YES;
  479.  
  480. load:
  481.   { int xfwd = ( vds > 0. ? YES : NO );
  482.     double
  483.       gdpr, gspr,
  484.       ieq_gd =      igd - ggd*vgd,
  485.       ieq_gs = ig - igd - ggs*vgs,
  486.       ieq_dr = id + igd - gds*vds - gm*( xfwd ? vgs : -vgd );
  487.  
  488.     /* load current vector */
  489.  
  490.     Y_MATRIX(b_ig) = -type*(ieq_gs + ieq_gd);
  491.     Y_MATRIX(b_id) =  type*(ieq_gd - ieq_dr);
  492.     Y_MATRIX(b_is) =  type*(ieq_gs + ieq_dr);
  493.  
  494.     /* load conductance matrix */
  495.  
  496.  
  497.     Y_MATRIX(b_ds) = ( xfwd ? -gm :  0. ) - gds;
  498.     Y_MATRIX(b_dg) = ( xfwd ?  gm : -gm ) - ggd;
  499.     Y_MATRIX(b_dd) = ( xfwd ?  0. :  gm ) + gds + ggd + ( gdpr = AREA*RD );
  500.  
  501.     Y_MATRIX(b_sd) = ( xfwd ?  0. : -gm ) - gds;
  502.     Y_MATRIX(b_sg) = ( xfwd ? -gm :  gm ) - ggs;
  503.     Y_MATRIX(b_ss) = ( xfwd ?  gm :  0. ) + gds + ggs + ( gspr = AREA*RS );
  504.  
  505.     Y_MATRIX(b_gd) = -ggd;
  506.     Y_MATRIX(b_gs) = -ggs;
  507.     Y_MATRIX(b_gg) =  ggd + ggs + RG;
  508.  
  509.     if ( LoadFl ) {
  510.       Y_MATRIX(b_dD) = Y_MATRIX(b_Dd) = -( Y_MATRIX(b_DD) = gdpr );
  511.       Y_MATRIX(b_Ss) = Y_MATRIX(b_sS) = -( Y_MATRIX(b_SS) = gspr );
  512.       Y_MATRIX(b_Gg) = Y_MATRIX(b_gG) = -( Y_MATRIX(b_GG) = RG );
  513.     } }
  514.  
  515.   B_VGS(0) = vgs;
  516.   B_VGD(0) = vgd;
  517.  
  518.   IG0 = ig ;
  519.   ID0 = id ;
  520.  
  521.   GM0  = gm ;
  522.   GDS0 = gds;
  523.   GGS0 = ggs;
  524.   GGD0 = ggd;
  525.  
  526. done:
  527.   return nonconv;
  528.   }
  529.  
  530. static double        /* calculate channel charge for Raytheon model */
  531.   GaAs2cap(model,vgs,vgd,vcap,cgso,cgdo,p_cgs,p_cgd)
  532.  
  533. struct B_ *model;
  534. double   vgs,vgd,vcap,cgso,cgdo,*p_cgs,*p_cgd;
  535.  
  536. /* Author
  537.   PWT   87-07-00        creation (after SPICE3A7)
  538.   PWT   87-10-29        correct Cgd formula, re-work code for speed
  539.   pwt    07 Nov 90    pass in *model instead of model parameters
  540. */
  541.  
  542. { double ext, qroot, qggval,
  543.     veroot = sqrt( (vgs - vgd)*(vgs - vgd) + vcap*vcap ),
  544.     veff1  = .5*(veroot + vgs + vgd),
  545.     veff2  = veff1 - veroot,
  546.     vnroot = sqrt( (veff1 - VTO)*(veff1 - VTO) + VDELTA*VDELTA ),
  547.     vnew1  = .5*(vnroot + veff1 + VTO),
  548.     vnew3  = vnew1;
  549.  
  550.   if ( vnew1 < VMAX )
  551.     ext = 0.;
  552.   else {
  553.     vnew1 = VMAX;
  554.     ext   = (vnew3-VMAX)/sqrt(1-VMAX/VBI);
  555.     }
  556.   qroot  = sqrt(1 - vnew1/VBI),
  557.   qggval = cgso*(2*VBI*(1-qroot) + ext) + cgdo*veff2;
  558.  
  559.   if ( p_cgs!=NULL && p_cgd!=NULL ) {    /* calculate capacitances */
  560.     double
  561.       par1   = .5*(1 + (veff1-VTO)/vnroot),
  562.       cgsx   = par1*cgso/qroot,
  563.       cfact  = (vgs-vgd)/veroot,
  564.       cplus  = .5*(1 + cfact),
  565.       cminus = cplus - cfact;
  566.  
  567.     *p_cgs = cgsx*cplus  + cgdo*cminus;
  568.     *p_cgd = cgsx*cminus + cgdo*cplus;
  569.     }
  570.  
  571.   return qggval;
  572.   }
  573.