bjt2load.c

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

                if(ckt->CKTmode & MODEINITPRED) {
                    xfact = ckt->CKTdelta/ckt->CKTdeltaOld[1];
                    *(ckt->CKTstate0 + here->BJT2vbe) = 
                            *(ckt->CKTstate1 + here->BJT2vbe);
                    vbe = (1+xfact)**(ckt->CKTstate1 + here->BJT2vbe)-
                            xfact* *(ckt->CKTstate2 + here->BJT2vbe);
                    *(ckt->CKTstate0 + here->BJT2vbc) = 
                            *(ckt->CKTstate1 + here->BJT2vbc);
                    vbc = (1+xfact)**(ckt->CKTstate1 + here->BJT2vbc)-
                            xfact* *(ckt->CKTstate2 + here->BJT2vbc);
                    *(ckt->CKTstate0 + here->BJT2cc) = 
                            *(ckt->CKTstate1 + here->BJT2cc);
                    *(ckt->CKTstate0 + here->BJT2cb) = 
                            *(ckt->CKTstate1 + here->BJT2cb);
                    *(ckt->CKTstate0 + here->BJT2gpi) = 
                            *(ckt->CKTstate1 + here->BJT2gpi);
                    *(ckt->CKTstate0 + here->BJT2gmu) = 
                            *(ckt->CKTstate1 + here->BJT2gmu);
                    *(ckt->CKTstate0 + here->BJT2gm) = 
                            *(ckt->CKTstate1 + here->BJT2gm);
                    *(ckt->CKTstate0 + here->BJT2go) = 
                            *(ckt->CKTstate1 + here->BJT2go);
                    *(ckt->CKTstate0 + here->BJT2gx) = 
                            *(ckt->CKTstate1 + here->BJT2gx);
                    *(ckt->CKTstate0 + here->BJT2vsub) = 
                            *(ckt->CKTstate1 + here->BJT2vsub);
                    vsub = (1+xfact)**(ckt->CKTstate1 + here->BJT2vsub)-
                            xfact* *(ckt->CKTstate2 + here->BJT2vsub);
                } else {
#endif /* PREDICTOR */
                    /*
                     *   compute new nonlinear branch voltages
                     */
                    vbe=model->BJT2type*(
                        *(ckt->CKTrhsOld+here->BJT2basePrimeNode)-
                        *(ckt->CKTrhsOld+here->BJT2emitPrimeNode));
                    vbc=model->BJT2type*(
                        *(ckt->CKTrhsOld+here->BJT2basePrimeNode)-
                        *(ckt->CKTrhsOld+here->BJT2colPrimeNode));
                    vsub=model->BJT2type*model->BJT2subs*(
                        *(ckt->CKTrhsOld+here->BJT2substNode)-
                        *(ckt->CKTrhsOld+here->BJT2substConNode));
#ifndef PREDICTOR
                }
#endif /* PREDICTOR */
                delvbe=vbe- *(ckt->CKTstate0 + here->BJT2vbe);
                delvbc=vbc- *(ckt->CKTstate0 + here->BJT2vbc);
                vbx=model->BJT2type*(
                    *(ckt->CKTrhsOld+here->BJT2baseNode)-
                    *(ckt->CKTrhsOld+here->BJT2colPrimeNode));
                vsub=model->BJT2type*model->BJT2subs*(
                    *(ckt->CKTrhsOld+here->BJT2substNode)-
                    *(ckt->CKTrhsOld+here->BJT2substConNode));
                cchat= *(ckt->CKTstate0 + here->BJT2cc)+(*(ckt->CKTstate0 + 
                        here->BJT2gm)+ *(ckt->CKTstate0 + here->BJT2go))*delvbe-
                        (*(ckt->CKTstate0 + here->BJT2go)+*(ckt->CKTstate0 +
                        here->BJT2gmu))*delvbc;
                cbhat= *(ckt->CKTstate0 + here->BJT2cb)+ *(ckt->CKTstate0 + 
                        here->BJT2gpi)*delvbe+ *(ckt->CKTstate0 + here->BJT2gmu)*
                        delvbc;
#ifndef NOBYPASS
                /*
                 *    bypass if solution has not changed
                 */
                /* the following collections of if's would be just one
                 * if the average compiler could handle it, but many
                 * find the expression too complicated, thus the split.
                 */
                if( (ckt->CKTbypass) &&
                        (!(ckt->CKTmode & MODEINITPRED)) &&
                        (fabs(delvbe) < (ckt->CKTreltol*MAX(fabs(vbe),
                            fabs(*(ckt->CKTstate0 + here->BJT2vbe)))+
                            ckt->CKTvoltTol)) )
                    if( (fabs(delvbc) < ckt->CKTreltol*MAX(fabs(vbc),
                            fabs(*(ckt->CKTstate0 + here->BJT2vbc)))+
                            ckt->CKTvoltTol) )
                    if( (fabs(cchat-*(ckt->CKTstate0 + here->BJT2cc)) < 
                            ckt->CKTreltol* MAX(fabs(cchat),
                            fabs(*(ckt->CKTstate0 + here->BJT2cc)))+
                            ckt->CKTabstol) )
                    if( (fabs(cbhat-*(ckt->CKTstate0 + here->BJT2cb)) < 
                            ckt->CKTreltol* MAX(fabs(cbhat),
                            fabs(*(ckt->CKTstate0 + here->BJT2cb)))+
                            ckt->CKTabstol) ) {
                    /*
                     * bypassing....
                     */
                    vbe = *(ckt->CKTstate0 + here->BJT2vbe);
                    vbc = *(ckt->CKTstate0 + here->BJT2vbc);
                    cc = *(ckt->CKTstate0 + here->BJT2cc);
                    cb = *(ckt->CKTstate0 + here->BJT2cb);
                    gpi = *(ckt->CKTstate0 + here->BJT2gpi);
                    gmu = *(ckt->CKTstate0 + here->BJT2gmu);
                    gm = *(ckt->CKTstate0 + here->BJT2gm);
                    go = *(ckt->CKTstate0 + here->BJT2go);
                    gx = *(ckt->CKTstate0 + here->BJT2gx);
                    geqcb = *(ckt->CKTstate0 + here->BJT2geqcb);
                    gcsub = *(ckt->CKTstate0 + here->BJT2gcsub);
                    geqbx = *(ckt->CKTstate0 + here->BJT2geqbx);
                    vsub = *(ckt->CKTstate0 + here->BJT2vsub);
                    gdsub = *(ckt->CKTstate0 + here->BJT2gdsub);
                    cdsub = *(ckt->CKTstate0 + here->BJT2cdsub);
                    goto load;
                }
#endif /*NOBYPASS*/
                /*
                 *   limit nonlinear branch voltages
                 */
                ichk1=1;
                vbe = DEVpnjlim(vbe,*(ckt->CKTstate0 + here->BJT2vbe),vt,
                        here->BJT2tVcrit,&icheck);
                vbc = DEVpnjlim(vbc,*(ckt->CKTstate0 + here->BJT2vbc),vt,
                        here->BJT2tVcrit,&ichk1);
                if (ichk1 == 1) icheck=1;
                vsub = DEVpnjlim(vsub,*(ckt->CKTstate0 + here->BJT2vsub),vt,
                        here->BJT2tSubVcrit,&ichk1);
                if (ichk1 == 1) icheck=1;
            }
            /*
             *   determine dc current and derivitives
             */
next1:      vtn=vt*model->BJT2emissionCoeffF;
            if(vbe >= -3*vtn){
                evbe=exp(vbe/vtn);
                cbe=csat*(evbe-1);
                gbe=csat*evbe/vtn;
            } else {
                arg=3*vtn/(vbe*CONSTe);
                arg = arg * arg * arg;
                cbe = -csat*(1+arg);
                gbe = csat*3*arg/vbe;
            }
            if (c2 == 0) {
                cben=0;
                gben=0;
            } else {
                if(vbe >= -3*vte){
                    evben=exp(vbe/vte);
                    cben=c2*(evben-1);
                    gben=c2*evben/vte;
                } else {
                    arg=3*vte/(vbe*CONSTe);
                    arg = arg * arg * arg;
                    cben = -c2*(1+arg);
                    gben = c2*3*arg/vbe;
                }
            }
            gben+=ckt->CKTgmin;
            cben+=ckt->CKTgmin*vbe;
            vtn=vt*model->BJT2emissionCoeffR;
            if(vbc >= -3*vtn) {
                evbc=exp(vbc/vtn);
                cbc=csat*(evbc-1);
                gbc=csat*evbc/vtn;
            } else {
                arg=3*vtn/(vbc*CONSTe);
                arg = arg * arg * arg;
                cbc = -csat*(1+arg);
                gbc = csat*3*arg/vbc;
            }
            if (c4 == 0) {
                cbcn=0;
                gbcn=0;
            } else {
                if(vbc >= -3*vtc) {
                    evbcn=exp(vbc/vtc);
                    cbcn=c4*(evbcn-1);
                    gbcn=c4*evbcn/vtc;
                } else {
                    arg=3*vtc/(vbc*CONSTe);
                    arg = arg * arg * arg;
                    cbcn = -c4*(1+arg);
                    gbcn = c4*3*arg/vbc;
                }
            }
            gbcn+=ckt->CKTgmin;
            cbcn+=ckt->CKTgmin*vbc;
            if(vsub <= -3*vt) {
                arg=3*vt/(vsub*CONSTe);
                arg = arg * arg * arg;
                gdsub = csubsat*3*arg/vsub+ckt->CKTgmin;
                cdsub = -csubsat*(1+arg)+ckt->CKTgmin*vsub;
            } else {
                evsub = exp(MIN(MAX_EXP_ARG,vsub/vt));
                gdsub = csubsat*evsub/vt + ckt->CKTgmin;
                cdsub = csubsat*(evsub-1) + ckt->CKTgmin*vsub;
            }
            /*
             *   determine base charge terms
             */
            q1=1/(1-model->BJT2invEarlyVoltF*vbc-model->BJT2invEarlyVoltR*vbe);
            if(oik == 0 && oikr == 0) {
                qb=q1;
                dqbdve=q1*qb*model->BJT2invEarlyVoltR;
                dqbdvc=q1*qb*model->BJT2invEarlyVoltF;
            } else {
                q2=oik*cbe+oikr*cbc;
                arg=MAX(0,1+4*q2);
                sqarg=1;
                if(arg != 0) sqarg=sqrt(arg);
                qb=q1*(1+sqarg)/2;
                dqbdve=q1*(qb*model->BJT2invEarlyVoltR+oik*gbe/sqarg);
                dqbdvc=q1*(qb*model->BJT2invEarlyVoltF+oikr*gbc/sqarg);
            }
            /*
             *   weil's approx. for excess phase applied with backward-
             *   euler integration
             */
            cc=0;
            cex=cbe;
            gex=gbe;
            if(ckt->CKTmode & (MODETRAN | MODEAC) && td != 0) {
                arg1=ckt->CKTdelta/td;
                arg2=3*arg1;
                arg1=arg2*arg1;
                denom=1+arg1+arg2;
                arg3=arg1/denom;
                if(ckt->CKTmode & MODEINITTRAN) {
                    *(ckt->CKTstate1 + here->BJT2cexbc)=cbe/qb;
                    *(ckt->CKTstate2 + here->BJT2cexbc)=
                            *(ckt->CKTstate1 + here->BJT2cexbc);
                }
                cc=(*(ckt->CKTstate1 + here->BJT2cexbc)*(1+ckt->CKTdelta/
                        ckt->CKTdeltaOld[1]+arg2)-
                        *(ckt->CKTstate2 + here->BJT2cexbc)*ckt->CKTdelta/
                        ckt->CKTdeltaOld[1])/denom;
                cex=cbe*arg3;
                gex=gbe*arg3;
                *(ckt->CKTstate0 + here->BJT2cexbc)=cc+cex/qb;
            }
            /*
             *   determine dc incremental conductances
             */
            cc=cc+(cex-cbc)/qb-cbc/here->BJT2tBetaR-cbcn;
            cb=cbe/here->BJT2tBetaF+cben+cbc/here->BJT2tBetaR+cbcn;
            gx=rbpr+rbpi/qb;
            if(xjrb != 0) {
                arg1=MAX(cb/xjrb,1e-9);
                arg2=(-1+sqrt(1+14.59025*arg1))/2.4317/sqrt(arg1);
                arg1=tan(arg2);
                gx=rbpr+3*rbpi*(arg1-arg2)/arg2/arg1/arg1;
            }
            if(gx != 0) gx=1/gx;
            gpi=gbe/here->BJT2tBetaF+gben;
            gmu=gbc/here->BJT2tBetaR+gbcn;
            go=(gbc+(cex-cbc)*dqbdvc/qb)/qb;
            gm=(gex-(cex-cbc)*dqbdve/qb)/qb-go;
            if( (ckt->CKTmode & (MODETRAN | MODEAC)) ||
                    ((ckt->CKTmode & MODETRANOP) && (ckt->CKTmode & MODEUIC)) ||
                    (ckt->CKTmode & MODEINITSMSIG)) {
                /*
                 *   charge storage elements
                 */
                tf=model->BJT2transitTimeF;
                tr=model->BJT2transitTimeR;
                czbe=here->BJT2tBEcap*here->BJT2area;
                pe=here->BJT2tBEpot;
                xme=model->BJT2junctionExpBE;
                cdis=model->BJT2baseFractionBCcap;
		
		if (model->BJT2subs == VERTICAL)
		    ctot=here->BJT2tBCcap*here->BJT2areab;
		else    
                    ctot=here->BJT2tBCcap*here->BJT2areac;
		    
		czbc=ctot*cdis;
                czbx=ctot-czbc;
                pc=here->BJT2tBCpot;
                xmc=model->BJT2junctionExpBC;
                fcpe=here->BJT2tDepCap;
		
		if (model->BJT2subs == VERTICAL)