bjtdset.c

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

/**********
Copyright 1990 Regents of the University of California.  All rights reserved.
Author: 1988 Jaijeet S Roychowdhury
**********/

#include "ngspice.h"
#include "cktdefs.h"
#include "bjtdefs.h"
#include "const.h"
#include "distodef.h"
#include "sperror.h"
#include "devdefs.h"
#include "suffix.h"

/*
 * This function initialises the Taylor coeffs for the
 * BJT's in the circuit
 */

/* actually load the current resistance value into the sparse matrix
 * previously provided */
void TanDeriv(Dderivs *new, Dderivs *old);
int BJTdSetup(GENmodel *inModel, CKTcircuit *ckt)
{
    BJTmodel *model = (BJTmodel*)inModel;
    BJTinstance *here;
    double arg;
    double c2;
    double c4;
    double lcapbe1, lcapbe2, lcapbe3;
    double lcapbx1, lcapbx2, lcapbx3;
    double cb = 0;
    double cbc;
    double cbcn;
    double cbe;
    double cben;
    double cdis;
    double csat;
    double ctot;
    double czbc;
    double czbcf2;
    double czbe;
    double czbef2;
    double czbx;
    double czbxf2;
    double czcs;
    double evbc;
    double evbcn;
    double evbe;
    double evben;
    double f1;
    double f2;
    double f3;
    double fcpc;
    double fcpe;
    double gbb1 = 0;
    double gbc;
    double gbcn;
    double gbe;
    double gbe2, gbe3;
    double gbc2,gbc3;
    double gben2,gben3;
    double gbcn2,gbcn3;
    double gben;
    double gbb2 = 0, gbb3 = 0;
    double oik;
    double oikr;
    double ovtf;
    double pc;
    double pe;
    double ps;
    double q1;
    double q2;
    double qb;
    double rbpi;
    double rbpr;
    double sarg;
    double sqarg;
    double tf;
    double tr;
    double vbc;
    double vbe;
    double vbx;
    double vsc;
    double vt;
    double vtc;
    double vte;
    double vtn;
    double xjrb;
    double xjtf;
    double xmc;
    double xme;
    double xms;
    double xtf;
    double vbed;
    double vbb;

    double lcapbc1 = 0.0;
    double lcapbc2 = 0.0;
    double lcapbc3 = 0.0;

    double lcapsc1 = 0.0;
    double lcapsc2 = 0.0;
    double lcapsc3 = 0.0;
    double ic;
    double dummy;
    Dderivs d_p, d_q, d_r;
    Dderivs d_dummy, d_q1, d_qb, d_dummy2;
    Dderivs d_arg, d_sqarg, d_ic, d_q2;
    Dderivs d_z, d_tanz, d_vbb, d_ibb, d_rbb;
    Dderivs d_ib, d_cbe, d_tff, d_qbe;

    d_p.value = 0.0;
    d_p.d1_p = 1.0;
    d_p.d1_q = 0.0;
    d_p.d1_r = 0.0;
    d_p.d2_p2 = 0.0;
    d_p.d2_q2 = 0.0;
    d_p.d2_r2 = 0.0;
    d_p.d2_pq = 0.0;
    d_p.d2_qr = 0.0;
    d_p.d2_pr = 0.0;
    d_p.d3_p3 = 0.0;
    d_p.d3_q3 = 0.0;
    d_p.d3_r3 = 0.0;
    d_p.d3_p2q = 0.0;
    d_p.d3_p2r = 0.0;
    d_p.d3_pq2 = 0.0;
    d_p.d3_q2r = 0.0;
    d_p.d3_pr2 = 0.0;
    d_p.d3_qr2 = 0.0;
    d_p.d3_pqr = 0.0;

    EqualDeriv(&d_q, &d_p);
    d_q.d1_q = 1.0;
    d_q.d1_p = 0.0;

    EqualDeriv(&d_r, &d_p);
    d_r.d1_r = 1.0;
    d_r.d1_p = 0.0;


    /*  loop through all the models */
    for( ; model != NULL; model = model->BJTnextModel ) {

        /* loop through all the instances of the model */
        for (here = model->BJTinstances; here != NULL ;
	     here=here->BJTnextInstance) {

            vt = here->BJTtemp * CONSTKoverQ;


            /*
             *   dc model paramters
             */
            csat=here->BJTtSatCur*here->BJTarea * here->BJTm;
            rbpr=model->BJTminBaseResist/(here->BJTarea * here->BJTm);
            rbpi=model->BJTbaseResist/(here->BJTarea * here->BJTm)-rbpr;
            oik=model->BJTinvRollOffF/(here->BJTarea * here->BJTm);
            c2=here->BJTtBEleakCur*here->BJTarea * here->BJTm;
            vte=model->BJTleakBEemissionCoeff*vt;
            oikr=model->BJTinvRollOffR/(here->BJTarea * here->BJTm);
            c4=here->BJTtBCleakCur*here->BJTareab * here->BJTm;
            vtc=model->BJTleakBCemissionCoeff*vt;
            xjrb=model->BJTbaseCurrentHalfResist*here->BJTarea * here->BJTm;


            /*
             *   initialization
             */
	    vbe= model->BJTtype*(*(ckt->CKTrhsOld + here->BJTbasePrimeNode) -
				 *(ckt->CKTrhsOld + here->BJTemitPrimeNode));
	    vbc= model->BJTtype*(*(ckt->CKTrhsOld + here->BJTbaseNode) -
				 *(ckt->CKTrhsOld + here->BJTcolPrimeNode));
	    vbx=model->BJTtype*(
		*(ckt->CKTrhsOld+here->BJTbaseNode)-
		*(ckt->CKTrhsOld+here->BJTcolPrimeNode));
	    vsc=model->BJTtype*(
		*(ckt->CKTrhsOld+here->BJTsubstNode)-
		*(ckt->CKTrhsOld+here->BJTcolPrimeNode));

	    vbb=model->BJTtype*(
		*(ckt->CKTrhsOld+here->BJTbaseNode) -
		*(ckt->CKTrhsOld+here->BJTbasePrimeNode));


	    vbed = vbe; /* this is just a dummy variable
			 * it is the delayed vbe to be
			 * used in the delayed gm generator 
			 */


	    /* ic = f1(vbe,vbc,vbed) + f2(vbc) + f3(vbc)
	     * 
	     * we shall calculate the taylor coeffs of ic wrt vbe,
	     * vbed, and vbc and store them away.  the equations f1 f2
	     * and f3 are given elsewhere; we shall start off with f1,
	     * compute derivs. upto third order and then do f2 and f3
	     * and add their derivatives.
	     *
	     * Since f1 above is a function of three variables, it
	     * will be convenient to use derivative structures to
	     * compute the derivatives of f1. For this computation,
	     * p=vbe, q=vbc, r=vbed.
	     *
	     * ib = f1(vbe) + f2(vbc) (not the same f's as above, in
	     * case you are wondering!)  the gbe's gbc's gben's and
	     * gbcn's are convenient subsidiary variables.
	     *
	     * irb = f(vbe, vbc, vbb) - the vbe & vbc dependencies
	     * arise from the qb term.
	     *
	     * qbe = f1(vbe,vbc) + f2(vbe)
	     *
	     * derivative structures will be used again in the above
	     * two equations. p=vbe, q=vbc, r=vbb.
	     *
	     * qbc = f(vbc) ; qbx = f(vbx)
	     *
	     * qss = f(vsc) */
            /*
             *   determine dc current and derivitives
             */
	    vtn=vt*model->BJTemissionCoeffF;
            if(vbe > -5*vtn){
                evbe=exp(vbe/vtn);
                cbe=csat*(evbe-1)+ckt->CKTgmin*vbe;
                gbe=csat*evbe/vtn+ckt->CKTgmin;
		gbe2 = csat*evbe/vtn/vtn;
		gbe3 = gbe2/vtn;

		/* note - these are actually derivs, not Taylor
		 * coeffs. - not divided by 2! and 3!
		 */
                if (c2 == 0) {
                    cben=0;
                    gben=gben2=gben3=0;
                } else {
                    evben=exp(vbe/vte);
                    cben=c2*(evben-1);
                    gben=c2*evben/vte;
		    gben2=gben/vte;
		    gben3=gben2/vte;
                }
            } else {
                gbe = -csat/vbe+ckt->CKTgmin;
		gbe2=gbe3=gben2=gben3=0;
                cbe=gbe*vbe;
                gben = -c2/vbe;
                cben=gben*vbe;
            }
            vtn=vt*model->BJTemissionCoeffR;
            if(vbc > -5*vtn) {
                evbc=exp(vbc/vtn);
                cbc=csat*(evbc-1)+ckt->CKTgmin*vbc;
                gbc=csat*evbc/vtn+ckt->CKTgmin;
		gbc2=csat*evbc/vtn/vtn;
		gbc3=gbc2/vtn;
                if (c4 == 0) {
                    cbcn=0;
                    gbcn=0;
		    gbcn2=gbcn3=0;
                } else {
                    evbcn=exp(vbc/vtc);
                    cbcn=c4*(evbcn-1);
                    gbcn=c4*evbcn/vtc;
		    gbcn2=gbcn/vtc;
		    gbcn3=gbcn2/vtc;
                }
            } else {
                gbc = -csat/vbc+ckt->CKTgmin;
		gbc2=gbc3=0;
                cbc = gbc*vbc;
                gbcn = -c4/vbc;
		gbcn2=gbcn3=0;
                cbcn=gbcn*vbc;
            }
            /*
             *   determine base charge terms
             */
	    /* q1 is a function of 2 variables p=vbe and q=vbc. r=
	     * anything
	     */
            q1=1/(1-model->BJTinvEarlyVoltF*vbc-model->BJTinvEarlyVoltR*vbe);
	    dummy = (1-model->BJTinvEarlyVoltF*vbc-
		     model->BJTinvEarlyVoltR*vbe);
	    EqualDeriv(&d_dummy, &d_p);
	    d_dummy.value = dummy;
	    d_dummy.d1_p = - model->BJTinvEarlyVoltR;
	    d_dummy.d1_q = -  model->BJTinvEarlyVoltF;

	    /* q1 = 1/dummy */
	    InvDeriv(&d_q1, &d_dummy);

	    /* now q1 and its derivatives are set up */
            if(oik == 0 && oikr == 0) {
                qb=q1;
		EqualDeriv(&d_qb, &d_q1);
            } else {
                q2=oik*cbe+oikr*cbc;
		EqualDeriv(&d_q2, &d_p);
		d_q2.value = q2;
		d_q2.d1_p = oik*gbe;
		d_q2.d1_q = oikr*gbc;
		d_q2.d2_p2 = oik*gbe2;
		d_q2.d2_q2 = oikr*gbc2;
		d_q2.d3_p3 = oik*gbe3;
		d_q2.d3_q3 = oikr*gbc3;
                arg=MAX(0,1+4*q2);
		if (arg == 0.)
		{
		    EqualDeriv(&d_arg,&d_p);
		    d_arg.d1_p = 0.0;
		}
		else
		{
		    TimesDeriv(&d_arg,&d_q2,4.0);
		    d_arg.value += 1.;
		}
                sqarg=1;
		EqualDeriv(&d_sqarg,&d_p);
		d_sqarg.value = 1.0;
		d_sqarg.d1_p = 0.0;
                if(arg != 0){
		    sqarg=sqrt(arg);
		    SqrtDeriv(&d_sqarg, &d_arg);
		}

                qb=q1*(1+sqarg)/2;
		dummy = 1 + sqarg;
		EqualDeriv(&d_dummy, &d_sqarg);
		d_dummy.value += 1.0;
		MultDeriv(&d_qb, &d_q1, &d_dummy);
		TimesDeriv(&d_qb, &d_qb, 0.5);
            }

	    ic = (cbe - cbc)/qb;
	    /* cbe is a fn of vbed only; cbc of vbc; and qb of vbe and vbc */
	    /* p=vbe, q=vbc, r=vbed; now dummy = cbe - cbc */
	    EqualDeriv(&d_dummy, &d_p);
	    d_dummy.d1_p = 0.0;
	    d_dummy.value = cbe-cbc;
	    d_dummy.d1_r = gbe;
	    d_dummy.d2_r2 = gbe2;
	    d_dummy.d3_r3 = gbe3;
	    d_dummy.d1_q = -gbc;
	    d_dummy.d2_q2 = -gbc2;
	    d_dummy.d3_q3 = -gbc3;

	    DivDeriv(&d_ic, &d_dummy, &d_qb);


	    d_ic.value -= cbc/here->BJTtBetaR + cbcn;
	    d_ic.d1_q -= gbc/here->BJTtBetaR + gbcn;