b2eval.c

spice中支持多层次元件模型仿真的可单独运行的插件源码

/**********
Copyright 1990 Regents of the University of California.  All rights reserved.
Author: 1988 Min-Chie Jeng, Hong J. Park, Thomas L. Quarles
**********/

#include "spice.h"
#include <stdio.h>
#include "util.h"
#include "cktdefs.h"
#include "bsim2def.h"
#include "trandefs.h"
#include "const.h"
#include "suffix.h"

/* This routine evaluates the drain current, its derivatives and the
 * charges associated with the gate,bulk and drain terminal
 * using the B2 (Berkeley Short-Channel IGFET Model) Equations.
 */
void
B2evaluate(Vds,Vbs,Vgs,here,model,gm,gds,gmb,qg,qb,qd,cgg,cgd,cgs,
        cbg,cbd,cbs,cdg,cdd,cds,Ids,von,vdsat,ckt)

    register CKTcircuit  *ckt;
    register B2model   *model;
    register B2instance *here;
    double Vds,Vbs,Vgs;
    double *gm,*gds,*gmb,*qg,*qb,*qd,*cgg,*cgd,*cgs,*cbg;
    double *cbd,*cbs,*cdg,*cdd,*cds,*Ids,*von,*vdsat;
 {
    double Vth, Vdsat;
    double Phisb, T1s, Eta, Gg, Aa, Inv_Aa, U1, U1s, Vc, Kk, SqrtKk;
    double dPhisb_dVb, dT1s_dVb, dVth_dVb, dVth_dVd, dAa_dVb, dVc_dVd;
    double dVc_dVg, dVc_dVb, dKk_dVc, dVdsat_dVd, dVdsat_dVg, dVdsat_dVb;
    double dUvert_dVg, dUvert_dVd, dUvert_dVb, dBeta0_dVd, dBeta0_dVb, Inv_Kk;
    double dUtot_dVd, dUtot_dVb, dUtot_dVg, Ai, Bi, Vghigh, Vglow, Vgeff, Vof;
    double Vbseff, Vgst, Vgdt, Qbulk, Utot;
    double T0, T1, T2, T3, T4, T5, Arg1, Arg2, Exp0;
    double tmp, tmp1, tmp2, tmp3, Uvert, Beta1, Beta2, Beta0, dGg_dVb, Exp1;
    double T6, T7, T8, T9, Usat, n, ExpArg, ExpArg1;
    double Beta, dQbulk_dVb, dVgdt_dVg, dVgdt_dVd;
    double dVbseff_dVb, Ua, Ub, dVgdt_dVb, dQbulk_dVd;
    double Con1, Con3, Con4, SqrVghigh, SqrVglow, CubVghigh, CubVglow;
    double delta, Coeffa, Coeffb, Coeffc, Coeffd, Inv_Uvert, Inv_Utot;
    double Inv_Vdsat, tanh, Sqrsech, dBeta1_dVb, dU1_dVd, dU1_dVg, dU1_dVb;
    double Betaeff, FR, dFR_dVd, dFR_dVg, dFR_dVb, Betas, Beta3, Beta4;
    double dBeta_dVd, dBeta_dVg, dBeta_dVb, dVgeff_dVg, dVgeff_dVd, dVgeff_dVb;
    double dCon3_dVd, dCon3_dVb, dCon4_dVd, dCon4_dVb, dCoeffa_dVd, dCoeffa_dVb;
    double dCoeffb_dVd, dCoeffb_dVb, dCoeffc_dVd, dCoeffc_dVb;
    double dCoeffd_dVd, dCoeffd_dVb;
    int ChargeComputationNeeded;
    int	valuetypeflag;			/* added  3/19/90 JSD   */

    if ((ckt->CKTmode & (MODEAC | MODETRAN)) ||
        ((ckt->CKTmode & MODETRANOP) && (ckt->CKTmode & MODEUIC)) ||
        (ckt->CKTmode & MODEINITSMSIG) )
    {   ChargeComputationNeeded  =  1;
    }
    else
    {   ChargeComputationNeeded  =  0;
    }

    if (Vbs < model->B2vbb2) Vbs = model->B2vbb2;
    if (Vgs > model->B2vgg2) Vgs = model->B2vgg2;
    if (Vds > model->B2vdd2) Vds = model->B2vdd2;

/* Threshold Voltage. */
    if (Vbs <= 0.0)
    {   Phisb = here->pParam->B2phi - Vbs;
        dPhisb_dVb = -1.0;
        T1s = sqrt(Phisb);
        dT1s_dVb = -0.5 / T1s; 
    }
    else
    {   tmp = here->pParam->B2phi / (here->pParam->B2phi + Vbs);
        Phisb = here->pParam->B2phi * tmp;
        dPhisb_dVb = -tmp * tmp;
        T1s = here->pParam->Phis3 / (here->pParam->B2phi + 0.5 * Vbs);
        dT1s_dVb = -0.5 * T1s * T1s / here->pParam->Phis3;
    }

    Eta = here->pParam->B2eta0 + here->pParam->B2etaB * Vbs;
    Ua = here->pParam->B2ua0 + here->pParam->B2uaB * Vbs;
    Ub = here->pParam->B2ub0 + here->pParam->B2ubB * Vbs;
    U1s = here->pParam->B2u10 + here->pParam->B2u1B * Vbs;

    Vth = here->pParam->B2vfb + here->pParam->B2phi + here->pParam->B2k1
        * T1s - here->pParam->B2k2 * Phisb - Eta * Vds;
    dVth_dVd = -Eta;
    dVth_dVb = here->pParam->B2k1 * dT1s_dVb + here->pParam->B2k2 
	     - here->pParam->B2etaB * Vds;

    Vgst = Vgs - Vth;

    tmp = 1.0 / (1.744 + 0.8364 * Phisb);
    Gg = 1.0 - tmp;
    dGg_dVb = 0.8364 * tmp * tmp * dPhisb_dVb;
    T0 = Gg / T1s;
    tmp1 = 0.5 * T0 * here->pParam->B2k1;
    Aa = 1.0 + tmp1;
    dAa_dVb = (Aa - 1.0) * (dGg_dVb / Gg - dT1s_dVb / T1s);
    Inv_Aa = 1.0 / Aa;

    Vghigh = here->pParam->B2vghigh;
    Vglow = here->pParam->B2vglow;

    if ((Vgst >= Vghigh) || (here->pParam->B2n0 == 0.0))
    {   Vgeff = Vgst;
	dVgeff_dVg = 1.0;
	dVgeff_dVd = -dVth_dVd;
	dVgeff_dVb = -dVth_dVb;
    }
    else
    {   Vof = here->pParam->B2vof0 + here->pParam->B2vofB * Vbs
	    + here->pParam->B2vofD * Vds;
	n = here->pParam->B2n0 + here->pParam->B2nB / T1s 
	  + here->pParam->B2nD * Vds;
	tmp = 0.5 / (n * model->B2Vtm);

	ExpArg1 = -Vds / model->B2Vtm;
	ExpArg1 = MAX(ExpArg1, -30.0);
	Exp1 = exp(ExpArg1);
	tmp1 = 1.0 - Exp1;
	tmp1 = MAX(tmp1, 1.0e-18);
	tmp2 = 2.0 * Aa * tmp1;

	if (Vgst <= Vglow)
	{
            ExpArg = Vgst * tmp;
	    ExpArg = MAX(ExpArg, -30.0);
	    Exp0 = exp(0.5 * Vof + ExpArg);
	    Vgeff = sqrt(tmp2) * model->B2Vtm * Exp0;
	    T0 = n * model->B2Vtm;
	    dVgeff_dVg = Vgeff * tmp;
	    dVgeff_dVd = dVgeff_dVg * (n / tmp1 * Exp1 - dVth_dVd - Vgst
		       * here->pParam->B2nD / n + T0 * here->pParam->B2vofD);
	    dVgeff_dVb = dVgeff_dVg * (here->pParam->B2vofB * T0
		       - dVth_dVb + here->pParam->B2nB * Vgst
		       / (n * T1s * T1s) * dT1s_dVb + T0 * Inv_Aa * dAa_dVb);
        }
	else
	{
            ExpArg = Vglow * tmp;
	    ExpArg = MAX(ExpArg, -30.0);
	    Exp0 = exp(0.5 * Vof + ExpArg);
	    Vgeff = sqrt(2.0 * Aa * (1.0 - Exp1)) * model->B2Vtm * Exp0;
	    Con1 = Vghigh;
	    Con3 = Vgeff;
	    Con4 = Con3 * tmp;
	    SqrVghigh = Vghigh * Vghigh;
	    SqrVglow = Vglow * Vglow;
	    CubVghigh = Vghigh * SqrVghigh;
	    CubVglow = Vglow * SqrVglow;
	    T0 = 2.0 * Vghigh;
	    T1 = 2.0 * Vglow;
	    T2 = 3.0 * SqrVghigh;
	    T3 = 3.0 * SqrVglow;
	    T4 = Vghigh - Vglow;
            T5 = SqrVghigh - SqrVglow;
            T6 = CubVghigh - CubVglow;
            T7 = Con1 - Con3;
            delta = (T1 - T0) * T6 + (T2 - T3) * T5 + (T0 * T3 - T1 * T2) * T4;
	    delta = 1.0 / delta;
            Coeffb = (T1 - Con4 * T0) * T6 + (Con4 * T2 - T3) * T5 
		   + (T0 * T3 - T1 * T2) * T7;
            Coeffc = (Con4 - 1.0) * T6 + (T2 - T3) * T7 + (T3 - Con4 * T2) * T4;
            Coeffd = (T1 - T0) * T7 + (1.0 - Con4) * T5 + (Con4 * T0 - T1) * T4;
            Coeffa = SqrVghigh * (Coeffc + Coeffd * T0);
            Vgeff = (Coeffa + Vgst * (Coeffb + Vgst * (Coeffc + Vgst * Coeffd)))
		  * delta;
	    dVgeff_dVg = (Coeffb + Vgst * (2.0 * Coeffc + 3.0 * Vgst * Coeffd))
		       * delta;
	    T7 = Con3 * tmp;
	    T8 = dT1s_dVb * here->pParam->B2nB / (T1s * T1s * n);
	    T9 = n * model->B2Vtm;
	    dCon3_dVd = T7 * (n * Exp1 / tmp1 - Vglow * here->pParam->B2nD 
		      / n + T9 * here->pParam->B2vofD);
	    dCon3_dVb = T7 * (T9 * Inv_Aa * dAa_dVb + Vglow * T8
		      + T9 * here->pParam->B2vofB);
	    dCon4_dVd = tmp * dCon3_dVd - T7 * here->pParam->B2nD / n;
	    dCon4_dVb = tmp * dCon3_dVb + T7 * T8;

            dCoeffb_dVd = dCon4_dVd * (T2 * T5 - T0 * T6) + dCon3_dVd 
			* (T1 * T2 - T0 * T3);
            dCoeffc_dVd = dCon4_dVd * (T6 - T2 * T4) + dCon3_dVd * (T3 - T2);
            dCoeffd_dVd = dCon4_dVd * (T0 * T4 - T5) + dCon3_dVd * (T0 - T1);
            dCoeffa_dVd = SqrVghigh * (dCoeffc_dVd + dCoeffd_dVd * T0);

	    dVgeff_dVd = -dVgeff_dVg * dVth_dVd + (dCoeffa_dVd + Vgst 
		       * (dCoeffb_dVd + Vgst * (dCoeffc_dVd + Vgst 
		       * dCoeffd_dVd))) * delta;

            dCoeffb_dVb = dCon4_dVb * (T2 * T5 - T0 * T6) + dCon3_dVb 
			* (T1 * T2 - T0 * T3);
            dCoeffc_dVb = dCon4_dVb * (T6 - T2 * T4) + dCon3_dVb * (T3 - T2);
            dCoeffd_dVb = dCon4_dVb * (T0 * T4 - T5) + dCon3_dVb * (T0 - T1);
            dCoeffa_dVb = SqrVghigh * (dCoeffc_dVb + dCoeffd_dVb * T0);

	    dVgeff_dVb = -dVgeff_dVg * dVth_dVb + (dCoeffa_dVb + Vgst 
		       * (dCoeffb_dVb + Vgst * (dCoeffc_dVb + Vgst 
		       * dCoeffd_dVb))) * delta;
        }
    }

    if (Vgeff > 0.0)
    {
        Uvert = 1.0 + Vgeff * (Ua + Vgeff * Ub);
        Uvert = MAX(Uvert, 0.2);
        Inv_Uvert = 1.0 / Uvert;
        T8 = Ua + 2.0 * Ub * Vgeff;
        dUvert_dVg = T8 * dVgeff_dVg;
        dUvert_dVd = T8 * dVgeff_dVd;
        dUvert_dVb = T8 * dVgeff_dVb + Vgeff * (here->pParam->B2uaB
	           + Vgeff * here->pParam->B2ubB);

        T8 = U1s * Inv_Aa * Inv_Uvert;
        Vc = T8 * Vgeff;
        T9 = Vc * Inv_Uvert;
        dVc_dVg = T8 * dVgeff_dVg - T9 * dUvert_dVg;
        dVc_dVd = T8 * dVgeff_dVd - T9 * dUvert_dVd;
        dVc_dVb = T8 * dVgeff_dVb + here->pParam->B2u1B * Vgeff * Inv_Aa 
	        * Inv_Uvert - Vc * Inv_Aa * dAa_dVb - T9 * dUvert_dVb;


        tmp2 = sqrt(1.0 + 2.0 * Vc);
        Kk = 0.5 * (1.0 + Vc + tmp2);
        Inv_Kk = 1.0 / Kk;
        dKk_dVc = 0.5  + 0.5 / tmp2;
        SqrtKk = sqrt(Kk);

        T8 = Inv_Aa / SqrtKk;
        Vdsat = Vgeff * T8;
        Vdsat = MAX(Vdsat, 1.0e-18);
        Inv_Vdsat = 1.0 / Vdsat;
        T9 = 0.5 * Vdsat * Inv_Kk * dKk_dVc;
        dVdsat_dVd = T8 * dVgeff_dVd - T9 * dVc_dVd;
        dVdsat_dVg = T8 * dVgeff_dVg - T9 * dVc_dVg;
        dVdsat_dVb = T8 * dVgeff_dVb - T9 * dVc_dVb - Vdsat* Inv_Aa * dAa_dVb;

        Beta0 = here->pParam->B2beta0 + here->pParam->B2beta0B * Vbs;
        Betas = here->pParam->B2betas0 + here->pParam->B2betasB * Vbs;
        Beta2 = here->pParam->B2beta20 + here->pParam->B2beta2B * Vbs
	      + here->pParam->B2beta2G * Vgs;
        Beta3 = here->pParam->B2beta30 + here->pParam->B2beta3B * Vbs
	      + here->pParam->B2beta3G * Vgs;
        Beta4 = here->pParam->B2beta40 + here->pParam->B2beta4B * Vbs
	      + here->pParam->B2beta4G * Vgs;
        Beta1 = Betas - (Beta0 + model->B2vdd * (Beta3 - model->B2vdd 
	      * Beta4)); 

        T0 = Vds * Beta2 * Inv_Vdsat;
        T0 = MIN(T0, 30.0);
        T1 = exp(T0);
        T2 = T1 * T1;
        T3 = T2 + 1.0;
        tanh = (T2 - 1.0) / T3;
        Sqrsech = 4.0 * T2 / (T3 * T3);

        Beta = Beta0 + Beta1 * tanh + Vds * (Beta3 - Beta4 * Vds);
        T4 = Beta1 * Sqrsech * Inv_Vdsat;
        T5 = model->B2vdd * tanh;
        dBeta_dVd = Beta3 - 2.0 * Beta4 * Vds + T4 * (Beta2 - T0 * dVdsat_dVd);
        dBeta_dVg = T4 * (here->pParam->B2beta2G * Vds - T0 * dVdsat_dVg) 
	          + here->pParam->B2beta3G * (Vds - T5) 
		  - here->pParam->B2beta4G * (Vds * Vds - model->B2vdd * T5);
        dBeta1_dVb = here->pParam->Arg;
        dBeta_dVb = here->pParam->B2beta0B + dBeta1_dVb * tanh + Vds
	          * (here->pParam->B2beta3B - Vds * here->pParam->B2beta4B)
	          + T4 * (here->pParam->B2beta2B * Vds - T0 * dVdsat_dVb);
	    

        if (Vgst > Vglow)
        {   
	    if (Vds <= Vdsat) /* triode region */
	    {
	        T3 = Vds * Inv_Vdsat;
	        T4 = T3 - 1.0;
	        T2 =  1.0 - here->pParam->B2u1D * T4 * T4;
	        U1 =  U1s * T2;
	        Utot = Uvert + U1 * Vds;
	        Utot = MAX(Utot, 0.5);
	        Inv_Utot = 1.0 / Utot;
	        T5 = 2.0 * U1s * here->pParam->B2u1D * Inv_Vdsat * T4;
	        dU1_dVd = T5 * (T3 * dVdsat_dVd - 1.0);
	        dU1_dVg = T5 * T3 * dVdsat_dVg;
	        dU1_dVb = T5 * T3 * dVdsat_dVb + here->pParam->B2u1B * T2;
	        dUtot_dVd = dUvert_dVd + U1 + Vds * dU1_dVd;
	        dUtot_dVg = dUvert_dVg + Vds * dU1_dVg;
	        dUtot_dVb = dUvert_dVb + Vds * dU1_dVb;

	        tmp1 = (Vgeff - 0.5 * Aa * Vds);
                tmp3 = tmp1 * Vds;
	        Betaeff = Beta * Inv_Utot;
                *Ids = Betaeff * tmp3;
	        T6 = *Ids / Betaeff * Inv_Utot;
	        *gds = T6 * (dBeta_dVd - Betaeff * dUtot_dVd) + Betaeff * (tmp1
	            + (dVgeff_dVd - 0.5 * Aa) * Vds); 
                *gm = T6 * (dBeta_dVg - Betaeff * dUtot_dVg) + Betaeff * Vds 
	           * dVgeff_dVg;
	        *gmb = T6 * (dBeta_dVb - Betaeff * dUtot_dVb) + Betaeff * Vds
	            * (dVgeff_dVb - 0.5 * Vds * dAa_dVb);
           } 
           else  /* Saturation */
           {  tmp1 = Vgeff * Inv_Aa * Inv_Kk;
              tmp3 = 0.5 * Vgeff * tmp1;
	      Betaeff = Beta * Inv_Uvert;
              *Ids = Betaeff * tmp3;
	      T0 = *Ids / Betaeff * Inv_Uvert;
	      T1 = Betaeff * Vgeff * Inv_Aa * Inv_Kk;
	      T2 = *Ids * Inv_Kk * dKk_dVc;

	      if (here->pParam->B2ai0 != 0.0)
	      {
	          Ai = here->pParam->B2ai0 + here->pParam->B2aiB * Vbs;
	          Bi = here->pParam->B2bi0 + here->pParam->B2biB * Vbs;
	          T5 = Bi / (Vds - Vdsat);