soi3load.c

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

             *     which is a bit complex mostly for reasons of continuity
             *     and robustness - c'est la vie.
             *     
             */

            /* the following variables are local to this code block until
             * it is obvious that they can be made global
             */
            double vg;
            double eta_s; /* 1+Citf/Cof */
            double gamma,sigma;
            double Egy = 0.0,vgy,Sgy;
            double psi_ss;
            double AL,A0,LL,L0;
            double ALx,A0x,EAL = 0.0,EA0 = 0.0;
            double logterm0,logtermL;
            double Vgconst0,VgconstL;
            double Egx0 = 0.0,EgxL = 0.0;
            double vgx0,vgxL;
            double psi_si0,psi_siL,psi_st0;
            double Esi0 = 0.0,EsiL = 0.0,Ess0 = 0.0;
            double theta2;
            double TVF; /* thermal velocity-saturation factor */
            double Mmob;
            double sqt0,sqt0one;
            double alpha,delta,Bf;
            double vGBF,vGT,vGBT;
            double PSI;
            double S;
            double psi_sLsat;
            double Esd0 = 0.0,EsdL = 0.0;

            double psi_s0,psi_sL;
            double fL,f0;
            double DfL_Dvgfb,DfL_Dvdb,DfL_Dvsb,DfL_DdeltaT;
            double Df0_Dvgfb,Df0_Dvdb,Df0_Dvsb,Df0_DdeltaT;
            double ld,Dld_Dvgfb,Dld_Dvdb,Dld_Dvsb,Dld_DdeltaT;

            double ich0;

            double pDpsi_si0_Dvgx0,pDpsi_siL_DvgxL;
            double Dvgx0_Dvgfb,pDvgx0_Dvg,Dvgx0_Dvsb,Dvgx0_Dvdb,Dvgx0_DdeltaT;
            double DvgxL_Dvgfb,pDvgxL_Dvg,DvgxL_Dvsb,DvgxL_Dvdb,DvgxL_DdeltaT;
            double pDpsi_si0_Dvsb,pDpsi_siL_Dvdb;
            double Dpsi_si0_Dvgfb,Dpsi_si0_Dvsb,Dpsi_si0_Dvdb;
            double Dpsi_siL_Dvgfb,Dpsi_siL_Dvsb,Dpsi_siL_Dvdb;
            double Dpsi_ss_Dvgfb,Dpsi_ss_Dvsb,Dpsi_ss_Dvdb;
            double Dpsi_st0_Dvgfb,Dpsi_st0_Dvsb,Dpsi_st0_Dvdb;
            double pDdelta_Dpsi_st0;
            double Dalpha_Dvgfb,Dalpha_Dvdb,Dalpha_Dvsb,Dalpha_DdeltaT;
            double DPSI_Dvgfb,DPSI_Dvdb,DPSI_Dvsb,DPSI_DdeltaT;
            double pDBf_Dpsi_st0,DvGT_Dvgfb,DvGT_Dvsb,DvGT_Dvdb,DvGT_DdeltaT;
            double D,DS_Dvgfb,DS_Dvsb,DS_Dvdb,DS_DdeltaT;
            double Dpsi_sLsat_Dvgfb,Dpsi_sLsat_Dvsb,Dpsi_sLsat_Dvdb;
            double gmg,gmd,gms,gmt;
            double Dpsi_si0_DdeltaT,Dpsi_siL_DdeltaT;
            double Dpsi_ss_DdeltaT,Dpsi_st0_DdeltaT;
            double Dpsi_sLsat_DdeltaT;
            double Dpsi_sL_Dvgfb,Dpsi_s0_Dvgfb,Dpsi_sL_Dvdb,Dpsi_s0_Dvdb;
            double Dpsi_sL_Dvsb,Dpsi_s0_Dvsb,Dpsi_sL_DdeltaT,Dpsi_s0_DdeltaT;

            double vdsat;
            double Dvdsat_Dvgfb,Dvdsat_Dvsb,Dvdsat_Dvdb,Dvdsat_DdeltaT;
            int    m,i;
            double vds2m,vdsat2m,Em,vdslim;
            double Dvdslim_Dvgfb,Dvdslim_Dvdb,Dvdslim_Dvsb,Dvdslim_DdeltaT;
            double Vmx;
            double DVmx_Dvgfb,DVmx_Dvdb,DVmx_Dvsb,DVmx_DdeltaT;
            double Vm1,Em1,Vm1x;
            double DVm1x_Dvgfb,DVm1x_Dvdb,DVm1x_Dvsb,DVm1x_DdeltaT;
            double vgeff,lm,lmeff,Elm,Dlm_Dvgfb,Dlm_Dvdb,Dlm_Dvsb,Dlm_DdeltaT;

            double Mminus1=0;
            double EM,betaM;
            double gMg,gMd,gMs;
            double Fm; /* mobility degradation factor */
            double Y; /* channel length modulation factor */
            double tmp,tmp1; /* temporary var to aid pre-calculation */
            double TMF; /* thermal mobility factor */
            int vgx0trans,vgxLtrans; /* flags to indicate if vg transform performed */
            int psi_s0_trans,psi_sL_trans=0;
            int Ess0_trans,Esd0_trans,EsdL_trans;
            int A0trans,ALtrans;

            double vT,EchiD,EchiD1;
            double BetaBJTeff;

/* Now we use a nasty trick - if device is A over T, must "swap" drain and source
   potentials.  we do a literal switch and change it back for the outside world. */
            if (here->SOI3mode == -1) {
              tmp = vsb;
              vsb = vdb;
              vdb = tmp;
            }

/* Intrinsic Electrical Bit - has a bit of thermal due to TTC */            
            sigma = model->SOI3sigma/(EffectiveLength);
            eta_s = 1 + (model->SOI3C_ssf/model->SOI3frontOxideCapFactor);
            vg = vgfb - here->SOI3tVfbF * model->SOI3type
                      + sigma*(here->SOI3mode*vds)
                      - model->SOI3chiFB*deltaT;
            gamma = model->SOI3gamma * (1 + (model->SOI3deltaW)/(here->SOI3w))
                                     * (1 - (model->SOI3deltaL)/EffectiveLength);
/* must introduce some limiting to prevent exp overflow */
            if (vg > (vt*MAX_EXP_ARG)) {
              vgy = vg;
            } else {
              Egy = exp(vg/vt);
              vgy = vt*log(1+Egy);
            }
            Sgy = sqrt((vgy/eta_s) + (gamma*gamma)/(4*eta_s*eta_s));
            psi_ss = (Sgy-0.5*gamma/eta_s)*(Sgy-0.5*gamma/eta_s);
            A0 = here->SOI3tPhi + vsb;
            AL = here->SOI3tPhi + vdb;
            if (A0 > (vt*MAX_EXP_ARG)) {
              A0x = A0;
              A0trans = 0;
            } else {
              EA0 = exp(A0/vt);
              A0x = vt * log(1+EA0);
              A0trans = 1;
            }
            if (AL > (vt*MAX_EXP_ARG)) {
              ALx = AL;
              ALtrans = 0;
            } else {
              EAL = exp(AL/vt);
              ALx = vt * log(1+EAL);
              ALtrans = 1;
            }
/* if vg is very large then don't need to transform
   vg into vgx.
*/
            Vgconst0 = eta_s*A0x - vt*(eta_s - 1) + gamma*sqrt(A0x);
            if ((vg-Vgconst0) < (5 * vt * MAX_EXP_ARG)) {
              vgx0trans = 1;
              Egx0 = exp((vg - Vgconst0)/(5*vt));
              vgx0 = (5*vt) * log(1 + Egx0) + Vgconst0;
            } else {
              vgx0trans = 0; /* no transform performed */
              vgx0 = vg;
            }


            VgconstL = eta_s*ALx - vt*(eta_s - 1) + gamma*sqrt(ALx);
            if ((vg-VgconstL) < (5 * vt * MAX_EXP_ARG)) {
              vgxLtrans = 1;
              EgxL = exp((vg - VgconstL)/(5*vt));
              vgxL = (5*vt) * log(1 + EgxL) + VgconstL;
            } else {
              vgxLtrans = 0; /* no transform performed */
              vgxL = vg;
            }

            L0 = vgx0 - eta_s*(A0 - vt);
            LL = vgxL - eta_s*(AL - vt);
            logterm0 = ((L0/gamma)*(L0/gamma) - A0)/vt;
            logtermL = ((LL/gamma)*(LL/gamma) - AL)/vt;
            if (logterm0<=0) { /* can only happen due to numerical problems with sqrt*/
              psi_si0 = A0 + model->SOI3chiPHI*deltaT + vt*log(vt/(gamma*gamma));
            } else {
              psi_si0 = A0 + model->SOI3chiPHI*deltaT + vt*log(logterm0);
            }
            if (logtermL<=0) {
              psi_siL = AL + model->SOI3chiPHI*deltaT + vt*log(vt/(gamma*gamma));
            } else {
              psi_siL = AL + model->SOI3chiPHI*deltaT + vt*log(logtermL);
            }

            if ((psi_si0-psi_ss) < vt*MAX_EXP_ARG) {
              Ess0 = exp((psi_si0-psi_ss)/vt);
              Ess0_trans = 1;
              if (psi_si0 > (vt * MAX_EXP_ARG)) { /* psi_si0 is BIG */
                psi_st0 = psi_si0 - vt*log(1+Ess0);
                psi_s0_trans = 0;
              } else {
                Esi0 = exp(psi_si0/vt);
                psi_st0 = vt*log(1+(Esi0/(1+Ess0)));
                psi_s0_trans = 1;
              }
            } else {
              psi_st0 = psi_ss;
              Ess0_trans = 0;
              psi_s0_trans = 0;
            }

/* now for psi_sLsat - has thermal influence ! (if vel sat is included) */
            TMF = exp(-model->SOI3k*log(1+(deltaT/here->SOI3temp)));
            sqt0 = sqrt(psi_st0+1E-25);
            sqt0one = sqrt(1+psi_st0);
/* delta is formulated this way so we can change it here and nowhere else */
            delta = 0.5/sqt0one;
            pDdelta_Dpsi_st0 = -0.25/((1+psi_st0)*sqt0one);
            alpha = eta_s+gamma*delta;
            vGBF = vg - psi_st0;
            Bf = sqt0 - delta*psi_st0;
            vGT = vg - gamma*Bf;
            vGBT =  vGT + alpha*vt;

            if(model->SOI3vsat == 0) {
              theta2 = 0;
              TVF = 0.0; /* could be any nominal value */
            } else {
              TVF  = 0.8*exp((here->SOI3temp+deltaT)/600);
              theta2 = (here->SOI3tSurfMob/(model->SOI3vsat))*1e-2 /*cm to m*/
                       *TMF*(1+TVF)/(EffectiveLength*(1+model->SOI3TVF0));
                        /* theta2=1/(L.Ec) */
            }
            Mmob = theta2 - 0.5*model->SOI3theta;

            PSI = MAX(0,(vGT/alpha) - psi_st0);
            S = 0.5*(1+sqrt(1 + (2*Mmob*PSI)/(1+model->SOI3theta*vGBF)));
            psi_sLsat = psi_st0 + PSI/S;

            if ((psi_si0 - psi_sLsat)<(vt*MAX_EXP_ARG)) {
              Esd0 = exp((psi_si0 - psi_sLsat)/vt);
              Esd0_trans = 1;
              if (psi_s0_trans) {
                psi_s0 = vt*log(1+Esi0/(1+Esd0));
              } else {
                psi_s0 = psi_si0 - vt*log(1+Esd0);
              }
            } else {
              psi_s0 = psi_sLsat;
              Esd0_trans = 0;
            }

            if ((psi_siL - psi_sLsat)<(vt*MAX_EXP_ARG)) {
              EsdL = exp((psi_siL - psi_sLsat)/vt);
              EsdL_trans = 1;
              if (psi_siL < (vt * MAX_EXP_ARG)) {
                EsiL = exp(psi_siL/vt);
                psi_sL = vt*log(1+EsiL/(1+EsdL));
                psi_sL_trans = 1;
              } else {
                psi_sL = psi_siL - vt*log(1+EsdL);
                psi_sL_trans = 0;
              }
            } else {
              psi_sL = psi_sLsat;
              EsdL_trans = 0;
            }
/* if, after all that, we have impossible situation due to numerical limiting schemes ... */

            if (psi_s0>psi_sL) {
              psi_sL = psi_s0;
            }

/* now we can finally work out surface potential for source and drain end */

/* now surface potential is ready */

            f0 = (vGBT - 0.5*alpha*psi_s0)*psi_s0;
            fL = (vGBT - 0.5*alpha*psi_sL)*psi_sL;

            ich0 = Beta * (fL - f0); /* This is "intrinsic" bit */

/* now for derivatives - they're a bit of a nightmare   *
 * notation  pD... means PARTIAL derivative whilst D... *
 * means FULL derivative                                */

/* first sub-threshold */

            if (vg > (vt*MAX_EXP_ARG)) {
              Dpsi_ss_Dvgfb = (1-0.5*gamma/(eta_s*Sgy))/eta_s;
            } else {
              Dpsi_ss_Dvgfb = (1-0.5*gamma/(eta_s*Sgy))*(Egy/(1+Egy))/eta_s;
            }
            Dpsi_ss_Dvsb  = Dpsi_ss_Dvgfb*(-sigma);
            Dpsi_ss_Dvdb  = Dpsi_ss_Dvgfb*(sigma);
            Dpsi_ss_DdeltaT = Dpsi_ss_Dvgfb*(-model->SOI3chiFB);

/* now for strong inversion */

            pDpsi_si0_Dvgx0 = 2*vt*L0/(L0*L0-gamma*gamma*A0);
            pDpsi_siL_DvgxL = 2*vt*LL/(LL*LL-gamma*gamma*AL);

  /* if vg is transformed, must have deriv accordingly, else it is 1 */

            if (vgx0trans) {
              Dvgx0_Dvgfb = pDvgx0_Dvg  = Egx0/(1+Egx0);
              if (A0trans) {
                /* JimB - 27/8/98 */
                /* Get divide by zero errors when A0x is equal to zero */
                /* Temporary fix - add small constant to denominator */
/*                Dvgx0_Dvsb  = pDvgx0_Dvg * (-sigma) +
                               (1-pDvgx0_Dvg)*(eta_s+0.5*gamma/sqrt(A0x))*
                               (EA0/(1+EA0)); */
                Dvgx0_Dvsb  = pDvgx0_Dvg * (-sigma) +
                               (1-pDvgx0_Dvg)*(eta_s+0.5*gamma/(sqrt(A0x)+1e-25))*
                               (EA0/(1+EA0));
                /* End JimB */
              } else {
                /* JimB - 27/8/97 */
                /* Temporary fix - add small constant to denominator */
                /*Dvgx0_Dvsb  = pDvgx0_Dvg * (-sigma) +
                               (1-pDvgx0_Dvg)*(eta_s+0.5*gamma/sqrt(A0)); */
                Dvgx0_Dvsb  = pDvgx0_Dvg * (-sigma) +
                               (1-pDvgx0_Dvg)*(eta_s+0.5*gamma/(sqrt(A0)+1e-25));