psmodel.c

ngspice又一个电子CAD仿真软件代码.功能更全

/*
 	Parker-Skellern MESFET model
 	
	Copyright (C) 1994, 1995, 1996  Macquarie University                    
	All Rights Reserved 
	Author: Anthony Parker
	Date:	2  Feb 1994  created
	        9  Feb 1994  correct NaN problem in strong cut-off region
	        20 MAR 1994  corrected capacitance initialization  
	        24 MAR 1994  added parameter MVST  
	        28 MAR 1994  reorganized declaration scopes 
            19 APR 1994  added new parameters: PS_HFETA, PS_HFE1, PS_HFE2,
                             PS_HFG1, and PS_HFG2
            18 May 1994  corrected 1/0 error when PS_VSUB=0
            15 Jul 1994  corrected errors in acload routine
            10 Aug 1995  added PS_VSUB to gds += gm*PS_VSUB*mvst*(vgt-vgst*(a..
			12 Sep 1995  changed _XXX to PS_XXX to aid portability
            13 Sep 1995  change to give arg=1-1/subfac; 
			                           if(vst!=0) gds+=gm*PS_VSUB..;
                                        gm *= arg;
			10 Feb 1996  change to names to match MicroSim code.
			5  Jul 1996  corrected diode eq (change Gmin*vgs to Gmin*vgd).
 
*****************************************************************************/

/*-----------
| functions defined in this file are:
    PSids()          returns dc drain source current and assigns other
                     current and branch conductances
    qgg()            static function that returns gate charge
    PScharge()       returns gate-source and gate-drain charge and capacitance
    PSacload()       returns small-signal conductance elements
    PSinstanceinit() initializes model parameters
 */

#define PSMODEL_C      /* activate local definitions in psmesfet.h */
#include "psmodel.h"

/*-----------
| dc current and conductance calculation */
double
PSids(ckt, model, here, vgs, vgd, igs, igd, ggs, ggd, Gm, Gds)
cref *ckt;
modl *model;
inst *here;
double vgs;
double vgd;
double *igs;
double *igd;
double *ggs;
double *ggd;
double *Gm;
double *Gds;
{
#define FX  -10.0 /* not too small else fatal rounding error in (rpt-a_rpt) */
#define MX  40.0  /* maximum exponential argument */
#define EMX 2.353852668370199842e17  /* exp(MX) */

   double idrain, arg;
   double area = AREA;

   { /* gate junction diodes */
      double zz;
      { /* gate-junction forward conduction */
         double Gmin   = GMIN;
         double Vt     = NVT;
         double isat   = IS   * area;
         if ((arg=vgs/Vt) > FX) {
            if(arg < MX) {
               *ggs=(zz=isat*exp(arg))/Vt+Gmin; *igs= zz -isat +Gmin*vgs;
            } else {
               *ggs=(zz=isat*EMX)/Vt+Gmin;     *igs=zz*(arg-MX+1)-isat+Gmin*vgs;
            }
         } else {
            *ggs = Gmin;                       *igs = -isat + Gmin * vgs;
         }
         if ((arg=vgd/Vt) > FX) {
            if(arg < MX) {
               *ggd=(zz=isat*exp(arg))/Vt+Gmin; *igd= zz -isat +Gmin*vgd;
            } else {
               *ggd=(zz=isat*EMX)/Vt+Gmin;     *igd=zz*(arg-MX+1)-isat+Gmin*vgd;
            }
         } else {
            *ggd = Gmin;                       *igd = -isat + Gmin * vgd;
         }
      }
      { /* gate-junction reverse 'breakdown' conduction */
         double Vbd    = VBD;
         double ibd    = IBD  * area;
         if ((arg=-vgs/Vbd) > FX) {
            if(arg < MX) {
               *ggs += (zz=ibd*exp(arg))/Vbd;  *igs -= zz-ibd;
            } else {
               *ggs += (zz=ibd*EMX)/Vbd;       *igs -= zz*((arg-MX)+1) - ibd;
            }
         } else                                *igs += ibd;
         if ((arg=-vgd/Vbd) > FX) {
            if(arg < MX) {
               *ggd += (zz=ibd*exp(arg))/Vbd;  *igd -= zz-ibd;
            } else {
               *ggd += (zz=ibd*EMX)/Vbd;       *igd -= zz*((arg-MX)+1) - ibd;
            }
         } else                                *igd += ibd;
      }
   }
    
   { /*  compute drain current and derivitives */
      double gm, gds;
      double vdst = vgs - vgd;
      double stepofour = STEP * FOURTH;
      { /* Include rate dependent threshold modulation */
         double vgst, dvgd, dvgs, h, vgdtrap, vgstrap, eta, gam;
         double vto = VTO;
         double LFg = LFGAM,   LFg1 = LFG1,  LFg2 = LFG2;
         double HFg = HFGAM,   HFg1 = HFG1,  HFg2 = HFG2;
         double HFe = HFETA,   HFe1 = HFE1,  HFe2 = HFE2;
         if(TRAN_ANAL) {
            double taug = TAUG;
            h = taug/(taug + stepofour);  h*=h; h*=h; /*4th power*/
            VGDTRAP_NOW = vgdtrap = h*VGDTRAP_BEFORE + (1-h) * vgd;
            VGSTRAP_NOW = vgstrap = h*VGSTRAP_BEFORE + (1-h) * vgs;
         } else {
            h = 0;
            VGDTRAP_NOW = vgdtrap = vgd;
            VGSTRAP_NOW = vgstrap = vgs;
         }
         vgst = vgs - vto;
         vgst -= (      LFg - LFg1*vgstrap + LFg2*vgdtrap)*vgdtrap;
         vgst += (eta = HFe - HFe1*vgdtrap + HFe2*vgstrap)*(dvgs = vgstrap-vgs);
         vgst += (gam = HFg - HFg1*vgstrap + HFg2*vgdtrap)*(dvgd = vgdtrap-vgd);
         { /* Exponential Subthreshold effect ids(vgst,vdst) */
            double vgt, subfac;
            double mvst = MVST;
            double vst = VSUB * (1 + mvst*vdst);
            if (vgst > FX*vst) {
               if (vgst > (arg=MX*vst)) { /* numerically large */
                  vgt = (EMX/(subfac = EMX+1))*(vgst-arg) + arg;
               } else  /* limit gate bias exponentially */
                  vgt = vst * log( subfac=(1 + exp(vgst/vst)) );
               { /* Dual Power-law ids(vgt,vdst) */
                  double mQ  = Q;
                  double PmQ = P - mQ;
                  double dvpd_dvdst=(double)D3*pow(vgt,PmQ);
                  double vdp = vdst*dvpd_dvdst; /*D3=P/Q/((VBI-vto)^PmQ)*/
                  { /* Early saturation effect ids(vgt,vdp) */
                     double za  = (double)ZA; /* sqrt(1 + Z)/2 */
                     double mxi = MXI;
                     double vsatFac = vgt/(mxi*vgt  + (double)XI_WOO);
                     double vsat=vgt/(1 + vsatFac);
                     double  aa = za*vdp+vsat/2.0;
                     double a_aa = aa-vsat;
                     double  rpt = sqrt( aa * aa + (arg=vsat*vsat*Z/4.0));
                     double a_rpt = sqrt(a_aa * a_aa + arg);
                     double vdt = (rpt - a_rpt);
                     double dvdt_dvdp = za * (aa/rpt - a_aa/a_rpt);
                     double dvdt_dvgt = (vdt - vdp*dvdt_dvdp)
                           *(1 + mxi*vsatFac*vsatFac)/(1 + vsatFac)/vgt;
                     { /* Intrinsic Q-law FET equation ids(vgt,vdt) */
                        gds=pow(vgt-vdt,mQ-1);
                        idrain = vdt*gds + vgt*(gm=pow(vgt,mQ-1)-gds);
                        gds *= mQ;
                        gm *= mQ;
                     }
                     gm += gds*dvdt_dvgt;
                     gds *= dvdt_dvdp;
                  }
                  gm += gds*PmQ*vdp/vgt;
                  gds *= dvpd_dvdst;
               }
               arg = 1 - 1/subfac;
               if(vst != 0) gds += gm*VSUB*mvst*(vgt-vgst*arg)/vst;
               gm *= arg;
            } else { /* in extreme cut-off (numerically) */
               idrain = gm = gds = 0.0e0;
            }
         }
         gds += gm*(arg = h*gam +
                   (1-h)*(HFe1*dvgs-HFg2*dvgd+2*LFg2*vgdtrap-LFg1*vgstrap+LFg));
         gm *= 1 - h*eta + (1-h)*(HFe2*dvgs -HFg1*dvgd + LFg1*vgdtrap) - arg;
      }
      { /* apply channel length modulation and beta scaling */
         double lambda = LAM;
         double beta   = BETA  * area;
         gm *= (arg = beta*(1 + lambda*vdst));
         gds = beta*lambda*idrain + gds*arg;
         idrain *= arg;
      }
        
      { /* apply thermal reduction of drain current */
        double h, pfac, pAverage;
        double delta = DELT / area;
        if(TRAN_ANAL) {
           double taud = TAUD;
           h = taud/(taud + stepofour);    h*=h; h*=h;
           POWR_NOW = pAverage = h*POWR_BEFORE + (1-h)*vdst*idrain;
        } else {
           POWR_NOW = POWR_BEFORE = pAverage = vdst*idrain;  h = 0;
        }
        idrain /= (pfac = 1+pAverage*delta);
        *Gm  = gm * (arg = (h*delta*POWR_BEFORE + 1)/pfac/pfac);
        *Gds = gds * arg - (1-h) * delta*idrain*idrain;
      }
   }
   return(idrain);
}

/*-----------
| code based on Statz et. al. capacitance model,  IEEE Tran ED Feb 87 */
static double
qgg(vgs, vgd, gamma, pb, alpha, vto, vmax, xc, cgso, cgdo, cgs, cgd)
double vgs, vgd, gamma, pb, alpha, vto, vmax, xc, cgso, cgdo, *cgs, *cgd;
{
   double qrt, ext, Cgso, cpm, cplus, cminus;
   double vds   = vgs-vgd;
   double d1_xc = 1-xc;
   double vert  = sqrt( vds * vds + alpha );
   double veff  = 0.5*(vgs + vgd + vert) + gamma*vds;
   double vnr   = d1_xc*(veff-vto);
   double vnrt  = sqrt( vnr*vnr + 0.04 );
   double vnew  = veff + 0.5*(vnrt - vnr);
   if ( vnew < vmax ) {
      ext  = 0;
      qrt  = sqrt(1 - vnew/pb);
      Cgso = 0.5*cgso/qrt*(1+xc + d1_xc*vnr/vnrt);
   } else {
      double vx  = 0.5*(vnew-vmax);
      double par = 1+vx/(pb-vmax);
      qrt  = sqrt(1 - vmax/pb);
      ext  = vx*(1 + par)/qrt;
      Cgso = 0.5*cgso/qrt*(1+xc + d1_xc*vnr/vnrt) * par;
   }
   cplus = 0.5*(1 + (cpm = vds/vert)); cminus = cplus - cpm;
   *cgs = Cgso*(cplus +gamma) + cgdo*(cminus+gamma);
   *cgd = Cgso*(cminus-gamma) + cgdo*(cplus -gamma);
   return(cgso*((pb+pb)*(1-qrt) + ext) + cgdo*(veff - vert));
}

/*-----------
| call during ac analysis initialisation or during transient analysis */
void
PScharge(ckt, model, here, vgs, vgd, capgs, capgd)
cref *ckt;
modl *model;
inst *here;
double vgs;
double vgd;
double *capgs;
double *capgd;
{
#define QGG(a,b,c,d) qgg(a,b,gac,phib,alpha,vto,vmax,xc,czgs,czgd,c,d)
/*  double qgg(); */
      
   double czgs = CGS * AREA;
   double czgd = CGD * AREA;
   double vto   = VTO;
   double alpha = (double)ALPHA;  /* (XI*woo/(XI+1)/2)^2 */
   double xc    = XC;
   double vmax  = VMAX;
   double phib  = VBI;
   double gac   = ACGAM;
   
   if(/*TRAN_INIT ||*/ !TRAN_ANAL)
       QGS_NOW = QGD_NOW = QGS_BEFORE = QGD_BEFORE
          = QGG(vgs,vgd,capgs,capgd);
   else {
      double cgsna,cgsnc;
      double cgdna,cgdnb, a_cap;
      double vgs1  = VGS1;
      double vgd1  = VGD1;
      double qgga=QGG(vgs ,vgd ,&cgsna,&cgdna);
      double qggb=QGG(vgs1,vgd ,&a_cap,&cgdnb);
      double qggc=QGG(vgs ,vgd1,&cgsnc,&a_cap);
      double qggd=QGG(vgs1,vgd1,&a_cap,&a_cap);
      QGS_NOW = QGS_BEFORE + 0.5 * (qgga-qggb + qggc-qggd);
      QGD_NOW = QGD_BEFORE + 0.5 * (qgga-qggc + qggb-qggd);
      *capgs = 0.5 * (cgsna + cgsnc);
      *capgd = 0.5 * (cgdna + cgdnb);
   }
}


/*-----------
| call for each frequency in ac analysis */
void
PSacload(ckt, model, here, vgs, vgd, ids, omega, Gm, xGm, Gds, xGds)
cref *ckt;
modl *model;
inst *here;
double vgs;
double vgd;
double ids;
double omega;
double *Gm;
double *xGm;
double *Gds;
double *xGds;
{
    double arg;
    double vds = vgs - vgd;
    double LFgam = LFGAM;
    double LFg1 = LFG1;
    double LFg2 = LFG2*vgd;
    double HFg1 = HFG1;
    double HFg2 = HFG2*vgd;
    double HFeta  = HFETA;
    double HFe1 = HFE1;
    double HFe2 = HFE2*vgs;
    double hfgam= HFGAM - HFg1*vgs + HFg2;
    double eta  = HFeta - HFe1*vgd + HFe2;
    double lfga = LFgam - LFg1*vgs + LFg2 + LFg2;
    double gmo  = *Gm/(1 - lfga + LFg1*vgd);
    
    double wtg = TAUG * omega;
    double wtgdet = 1 + wtg*wtg;
    double gwtgdet = gmo/wtgdet;

    double gdsi = (arg=hfgam - lfga)*gwtgdet;
    double gdsr = arg*gmo - gdsi;
    double gmi  = (eta + LFg1*vgd)*gwtgdet + gdsi;
 
    double xgds = wtg*gdsi;
    double  gds = *Gds + gdsr;
    double  xgm = -wtg*gmi;
    double   gm = gmi + gmo*(1 -eta - hfgam);
                
    double delta = DELT / AREA;
    double wtd = TAUD * omega ;
    double wtddet = 1 + wtd * wtd;
    double fac = delta * ids;
    double del = 1/(1 - fac * vds);
    double dd = (del-1) / wtddet;
    double dr = del - dd;
    double di = wtd * dd;
                
    double cdsqr = fac * ids * del * wtd/wtddet;

    *Gm   = dr*gm  - di*xgm;
    *xGm  = di*gm  + dr*xgm;
 
    *Gds  = dr*gds - di*xgds + cdsqr*wtd;
    *xGds = di*gds + dr*xgds + cdsqr;
}


void    /* call when temperature changes */
PSinstanceinit(model, here)
modl *model;
inst *here;
{
#ifndef PARAM_CAST     /* allow casting to parameter type */
#define PARAM_CAST     /* if not specified then don't cast */
#endif
    
    double woo = (VBI - VTO);
    XI_WOO = PARAM_CAST (XI * woo);
    ZA     = PARAM_CAST (sqrt(1 + Z)/2);
    ALPHA  = PARAM_CAST (XI_WOO*XI_WOO/(XI+1)/(XI+1)/ 4);
    D3     = PARAM_CAST (P/Q/pow(woo,(P - Q)));
}