1 /**********
2 Copyright 1990 Regents of the University of California.  All rights reserved.
3 Author: 1988 Jaijeet S Roychowdhury
4 **********/
5 
6 #include "ngspice/ngspice.h"
7 #include "ngspice/cktdefs.h"
8 #include "bjtdefs.h"
9 #include "ngspice/const.h"
10 #include "ngspice/distodef.h"
11 #include "ngspice/sperror.h"
12 #include "ngspice/devdefs.h"
13 #include "ngspice/suffix.h"
14 
15 /*
16  * This function initialises the Taylor coeffs for the
17  * BJT's in the circuit
18  */
19 
20 /* actually load the current resistance value into the sparse matrix
21  * previously provided */
22 
BJTdSetup(GENmodel * inModel,CKTcircuit * ckt)23 int BJTdSetup(GENmodel *inModel, CKTcircuit *ckt)
24 {
25     BJTmodel *model = (BJTmodel*)inModel;
26     BJTinstance *here;
27     double arg;
28     double c2;
29     double c4;
30     double lcapbe1, lcapbe2, lcapbe3;
31     double lcapbx1, lcapbx2, lcapbx3;
32     double cb = 0.0;
33     double cbc;
34     double cbcn;
35     double cbe;
36     double cben;
37     double cdis;
38     double csat;
39     double ctot;
40     double czbc;
41     double czbcf2;
42     double czbe;
43     double czbef2;
44     double czbx;
45     double czbxf2;
46     double czcs;
47     double evbc;
48     double evbcn;
49     double evbe;
50     double evben;
51     double f1;
52     double f2;
53     double f3;
54     double fcpc;
55     double fcpe;
56     double gbb1 = 0.0;
57     double gbc;
58     double gbcn;
59     double gbe;
60     double gbe2, gbe3;
61     double gbc2,gbc3;
62     double gben2,gben3;
63     double gbcn2,gbcn3;
64     double gben;
65     double gbb2 = 0.0, gbb3 = 0.0;
66     double oik;
67     double oikr;
68     double ovtf;
69     double pc;
70     double pe;
71     double ps;
72     double q1;
73     double q2;
74     double qb;
75     double rbpi;
76     double rbpr;
77     double sarg;
78     double sqarg;
79     double tf;
80     double tr;
81     double vbc;
82     double vbe;
83     double vbx;
84     double vsc;
85     double vt;
86     double vtc;
87     double vte;
88     double vtn;
89     double xjrb;
90     double xjtf;
91     double xmc;
92     double xme;
93     double xms;
94     double xtf;
95     double vbed;
96     double vbb;
97 
98     double lcapbc1 = 0.0;
99     double lcapbc2 = 0.0;
100     double lcapbc3 = 0.0;
101 
102     double lcapsc1 = 0.0;
103     double lcapsc2 = 0.0;
104     double lcapsc3 = 0.0;
105     double ic;
106     double dummy;
107     Dderivs d_p, d_q, d_r;
108     Dderivs d_dummy, d_q1, d_qb, d_dummy2;
109     Dderivs d_arg, d_sqarg, d_ic, d_q2;
110     Dderivs d_z, d_tanz, d_vbb, d_ibb, d_rbb;
111     Dderivs d_ib, d_cbe, d_tff, d_qbe;
112 
113     d_p.value = 0.0;
114     d_p.d1_p = 1.0;
115     d_p.d1_q = 0.0;
116     d_p.d1_r = 0.0;
117     d_p.d2_p2 = 0.0;
118     d_p.d2_q2 = 0.0;
119     d_p.d2_r2 = 0.0;
120     d_p.d2_pq = 0.0;
121     d_p.d2_qr = 0.0;
122     d_p.d2_pr = 0.0;
123     d_p.d3_p3 = 0.0;
124     d_p.d3_q3 = 0.0;
125     d_p.d3_r3 = 0.0;
126     d_p.d3_p2q = 0.0;
127     d_p.d3_p2r = 0.0;
128     d_p.d3_pq2 = 0.0;
129     d_p.d3_q2r = 0.0;
130     d_p.d3_pr2 = 0.0;
131     d_p.d3_qr2 = 0.0;
132     d_p.d3_pqr = 0.0;
133 
134     EqualDeriv(&d_q, &d_p);
135     d_q.d1_q = 1.0;
136     d_q.d1_p = 0.0;
137 
138     EqualDeriv(&d_r, &d_p);
139     d_r.d1_r = 1.0;
140     d_r.d1_p = 0.0;
141 
142 
143     /*  loop through all the models */
144     for( ; model != NULL; model = BJTnextModel(model)) {
145 
146         /* loop through all the instances of the model */
147         for (here = BJTinstances(model); here != NULL ;
148              here=BJTnextInstance(here)) {
149 
150             vt = here->BJTtemp * CONSTKoverQ;
151 
152 
153             /*
154              *   dc model paramters
155              */
156             csat=here->BJTtSatCur*here->BJTarea * here->BJTm;
157             rbpr=here->BJTtminBaseResist/(here->BJTarea * here->BJTm);
158             rbpi=here->BJTtbaseResist/(here->BJTarea * here->BJTm)-rbpr;
159             oik=here->BJTtinvRollOffF/(here->BJTarea * here->BJTm);
160             c2=here->BJTtBEleakCur*here->BJTarea * here->BJTm;
161             vte=here->BJTtleakBEemissionCoeff*vt;
162             oikr=here->BJTtinvRollOffR/(here->BJTarea * here->BJTm);
163             c4=here->BJTtBCleakCur * here->BJTm;
164             if (model->BJTsubs == VERTICAL)
165                 c4 *= here->BJTareab;
166             else
167                 c4 *= here->BJTareac; /* lateral transistor */
168             vtc=here->BJTtleakBCemissionCoeff*vt;
169             xjrb=here->BJTtbaseCurrentHalfResist*here->BJTarea * here->BJTm;
170 
171 
172             /*
173              *   initialization
174              */
175             vbe= model->BJTtype*(*(ckt->CKTrhsOld + here->BJTbasePrimeNode) -
176                                  *(ckt->CKTrhsOld + here->BJTemitPrimeNode));
177             vbc= model->BJTtype*(*(ckt->CKTrhsOld + here->BJTbaseNode) -
178                                  *(ckt->CKTrhsOld + here->BJTcolPrimeNode));
179             vbx=model->BJTtype*(
180                 *(ckt->CKTrhsOld+here->BJTbaseNode)-
181                 *(ckt->CKTrhsOld+here->BJTcolPrimeNode));
182             vsc=model->BJTtype*(
183                 *(ckt->CKTrhsOld+here->BJTsubstNode)-
184                 *(ckt->CKTrhsOld+here->BJTcolPrimeNode));
185 
186             vbb=model->BJTtype*(
187                 *(ckt->CKTrhsOld+here->BJTbaseNode) -
188                 *(ckt->CKTrhsOld+here->BJTbasePrimeNode));
189 
190 
191             vbed = vbe; /* this is just a dummy variable
192                          * it is the delayed vbe to be
193                          * used in the delayed gm generator
194                          */
195 
196 
197             /* ic = f1(vbe,vbc,vbed) + f2(vbc) + f3(vbc)
198              *
199              * we shall calculate the taylor coeffs of ic wrt vbe,
200              * vbed, and vbc and store them away.  the equations f1 f2
201              * and f3 are given elsewhere; we shall start off with f1,
202              * compute derivs. upto third order and then do f2 and f3
203              * and add their derivatives.
204              *
205              * Since f1 above is a function of three variables, it
206              * will be convenient to use derivative structures to
207              * compute the derivatives of f1. For this computation,
208              * p=vbe, q=vbc, r=vbed.
209              *
210              * ib = f1(vbe) + f2(vbc) (not the same f's as above, in
211              * case you are wondering!)  the gbe's gbc's gben's and
212              * gbcn's are convenient subsidiary variables.
213              *
214              * irb = f(vbe, vbc, vbb) - the vbe & vbc dependencies
215              * arise from the qb term.
216              *
217              * qbe = f1(vbe,vbc) + f2(vbe)
218              *
219              * derivative structures will be used again in the above
220              * two equations. p=vbe, q=vbc, r=vbb.
221              *
222              * qbc = f(vbc) ; qbx = f(vbx)
223              *
224              * qss = f(vsc) */
225             /*
226              *   determine dc current and derivitives
227              */
228             vtn=vt*here->BJTtemissionCoeffF;
229             if(vbe > -5*vtn){
230                 evbe=exp(vbe/vtn);
231                 cbe=csat*(evbe-1)+ckt->CKTgmin*vbe;
232                 gbe=csat*evbe/vtn+ckt->CKTgmin;
233                 gbe2 = csat*evbe/vtn/vtn;
234                 gbe3 = gbe2/vtn;
235 
236                 /* note - these are actually derivs, not Taylor
237                  * coeffs. - not divided by 2! and 3!
238                  */
239                 if (c2 == 0) {
240                     cben=0;
241                     gben=gben2=gben3=0;
242                 } else {
243                     evben=exp(vbe/vte);
244                     cben=c2*(evben-1);
245                     gben=c2*evben/vte;
246                     gben2=gben/vte;
247                     gben3=gben2/vte;
248                 }
249             } else {
250                 gbe = -csat/vbe+ckt->CKTgmin;
251                 gbe2=gbe3=gben2=gben3=0;
252                 cbe=gbe*vbe;
253                 gben = -c2/vbe;
254                 cben=gben*vbe;
255             }
256             vtn=vt*here->BJTtemissionCoeffR;
257             if(vbc > -5*vtn) {
258                 evbc=exp(vbc/vtn);
259                 cbc=csat*(evbc-1)+ckt->CKTgmin*vbc;
260                 gbc=csat*evbc/vtn+ckt->CKTgmin;
261                 gbc2=csat*evbc/vtn/vtn;
262                 gbc3=gbc2/vtn;
263                 if (c4 == 0) {
264                     cbcn=0;
265                     gbcn=0;
266                     gbcn2=gbcn3=0;
267                 } else {
268                     evbcn=exp(vbc/vtc);
269                     cbcn=c4*(evbcn-1);
270                     gbcn=c4*evbcn/vtc;
271                     gbcn2=gbcn/vtc;
272                     gbcn3=gbcn2/vtc;
273                 }
274             } else {
275                 gbc = -csat/vbc+ckt->CKTgmin;
276                 gbc2=gbc3=0;
277                 cbc = gbc*vbc;
278                 gbcn = -c4/vbc;
279                 gbcn2=gbcn3=0;
280                 cbcn=gbcn*vbc;
281             }
282             /*
283              *   determine base charge terms
284              */
285             /* q1 is a function of 2 variables p=vbe and q=vbc. r=
286              * anything
287              */
288             q1=1/(1-here->BJTtinvEarlyVoltF*vbc-here->BJTtinvEarlyVoltR*vbe);
289             dummy = (1-here->BJTtinvEarlyVoltF*vbc-
290                      here->BJTtinvEarlyVoltR*vbe);
291             EqualDeriv(&d_dummy, &d_p);
292             d_dummy.value = dummy;
293             d_dummy.d1_p = - here->BJTtinvEarlyVoltR;
294             d_dummy.d1_q = - here->BJTtinvEarlyVoltF;
295 
296             /* q1 = 1/dummy */
297             InvDeriv(&d_q1, &d_dummy);
298 
299             /* now q1 and its derivatives are set up */
300             if(oik == 0 && oikr == 0) {
301                 qb=q1;
302                 EqualDeriv(&d_qb, &d_q1);
303             } else {
304                 q2=oik*cbe+oikr*cbc;
305                 EqualDeriv(&d_q2, &d_p);
306                 d_q2.value = q2;
307                 d_q2.d1_p = oik*gbe;
308                 d_q2.d1_q = oikr*gbc;
309                 d_q2.d2_p2 = oik*gbe2;
310                 d_q2.d2_q2 = oikr*gbc2;
311                 d_q2.d3_p3 = oik*gbe3;
312                 d_q2.d3_q3 = oikr*gbc3;
313                 arg=MAX(0,1+4*q2);
314                 if (arg == 0.)
315                 {
316                     EqualDeriv(&d_arg,&d_p);
317                     d_arg.d1_p = 0.0;
318                 }
319                 else
320                 {
321                     TimesDeriv(&d_arg,&d_q2,4.0);
322                     d_arg.value += 1.;
323                 }
324                 sqarg=1;
325                 EqualDeriv(&d_sqarg,&d_p);
326                 d_sqarg.value = 1.0;
327                 d_sqarg.d1_p = 0.0;
328                 if(arg != 0){
329                     sqarg=sqrt(arg);
330                     SqrtDeriv(&d_sqarg, &d_arg);
331                 }
332 
333                 qb=q1*(1+sqarg)/2;
334                 dummy = 1 + sqarg;
335                 EqualDeriv(&d_dummy, &d_sqarg);
336                 d_dummy.value += 1.0;
337                 MultDeriv(&d_qb, &d_q1, &d_dummy);
338                 TimesDeriv(&d_qb, &d_qb, 0.5);
339             }
340 
341             ic = (cbe - cbc)/qb;
342             /* cbe is a fn of vbed only; cbc of vbc; and qb of vbe and vbc */
343             /* p=vbe, q=vbc, r=vbed; now dummy = cbe - cbc */
344             EqualDeriv(&d_dummy, &d_p);
345             d_dummy.d1_p = 0.0;
346             d_dummy.value = cbe-cbc;
347             d_dummy.d1_r = gbe;
348             d_dummy.d2_r2 = gbe2;
349             d_dummy.d3_r3 = gbe3;
350             d_dummy.d1_q = -gbc;
351             d_dummy.d2_q2 = -gbc2;
352             d_dummy.d3_q3 = -gbc3;
353 
354             DivDeriv(&d_ic, &d_dummy, &d_qb);
355 
356 
357             d_ic.value -= cbc/here->BJTtBetaR + cbcn;
358             d_ic.d1_q -= gbc/here->BJTtBetaR + gbcn;
359             d_ic.d2_q2 -= gbc2/here->BJTtBetaR + gbcn2;
360             d_ic.d3_q3 -= gbc3/here->BJTtBetaR + gbcn3;
361 
362             /* check this point: where is the f2(vbe) contribution to ic ? */
363 
364             /* ic derivatives all set up now */
365             /* base spread resistance */
366 
367             if ( !((rbpr == 0.0) && (rbpi == 0.0)))
368             {
369                 cb=cbe/here->BJTtBetaF+cben+cbc/here->BJTtBetaR+cbcn;
370                 /* we are calculating derivatives w.r.t cb itself */
371                 /*
372                   gx=rbpr+rbpi/qb;
373                 */
374 
375                 if (cb != 0.0) {
376                     if((xjrb != 0.0) && (rbpi != 0.0)) {
377                         /* p = ib, q, r = anything */
378                         dummy=MAX(cb/xjrb,1e-9);
379                         EqualDeriv(&d_dummy, &d_p);
380                         d_dummy.value = dummy;
381                         d_dummy.d1_p = 1/xjrb;
382                         SqrtDeriv(&d_dummy, &d_dummy);
383                         TimesDeriv(&d_dummy, &d_dummy, 2.4317);
384 
385                         /*
386                           dummy2=(-1+sqrt(1+14.59025*MAX(cb/xjrb,1e-9)));
387                         */
388                         EqualDeriv(&d_dummy2, &d_p);
389                         d_dummy2.value = 1+14.59025*MAX(cb/xjrb,1e-9);
390                         d_dummy2.d1_p = 14.59025/xjrb;
391                         SqrtDeriv(&d_dummy2, &d_dummy2);
392                         d_dummy2.value -= 1.0;
393 
394                         DivDeriv(&d_z, &d_dummy2, &d_dummy);
395                         TanDeriv(&d_tanz, &d_z);
396 
397                         /* now using dummy = tanz - z and dummy2 =
398                            z*tanz*tanz */
399                         TimesDeriv(&d_dummy, &d_z, -1.0);
400                         PlusDeriv(&d_dummy, &d_dummy, &d_tanz);
401 
402                         MultDeriv(&d_dummy2, &d_tanz, &d_tanz);
403                         MultDeriv(&d_dummy2, &d_dummy2, &d_z);
404 
405                         DivDeriv(&d_rbb , &d_dummy, &d_dummy2);
406                         TimesDeriv(&d_rbb,&d_rbb, 3.0*rbpi);
407                         d_rbb.value += rbpr;
408 
409                         MultDeriv(&d_vbb, &d_rbb, &d_p);
410 
411                         /* power series inversion to get the
412                            conductance derivatives */
413 
414                         if (d_vbb.d1_p != 0) {
415                             gbb1 = 1/d_vbb.d1_p;
416                             gbb2 = -(d_vbb.d2_p2*0.5)*gbb1*gbb1;
417                             gbb3 = gbb1*gbb1*gbb1*gbb1*(-(d_vbb.d3_p3/6.0)
418                                                         + 2*(d_vbb.d2_p2*0.5)*(d_vbb.d2_p2*0.5)*gbb1);
419                         }
420                         else
421                             printf("\nd_vbb.d1_p = 0 in base spread resistance calculations\n");
422 
423 
424                         /* r = vbb */
425                         EqualDeriv(&d_ibb, &d_r);
426                         d_ibb.value = cb;
427                         d_ibb.d1_r = gbb1;
428                         d_ibb.d2_r2 = 2*gbb2;
429                         d_ibb.d3_r3 = 6.0*gbb3;
430                     }
431                     else
432                     {
433                         /*
434                           rbb = rbpr + rbpi/qb;
435                           ibb = vbb /rbb; = f(vbe, vbc, vbb)
436                         */
437 
438                         EqualDeriv(&d_rbb,&d_p);
439                         d_rbb.d1_p = 0.0;
440                         if (rbpi != 0.0) {
441                             InvDeriv(&d_rbb, &d_qb);
442                             TimesDeriv(&d_rbb, &d_rbb,rbpi);
443                         }
444                         d_rbb.value += rbpr;
445 
446                         EqualDeriv(&d_ibb,&d_r);
447                         d_ibb.value = vbb;
448                         DivDeriv(&d_ibb,&d_ibb,&d_rbb);
449 
450                     }
451                 }
452                 else
453                 {
454                     EqualDeriv(&d_ibb,&d_r);
455                     if (rbpr != 0.0)
456                         d_ibb.d1_r = 1/rbpr;
457                 }
458             }
459             else
460             {
461                 EqualDeriv(&d_ibb,&d_p);
462                 d_ibb.d1_p = 0.0;
463             }
464 
465             /* formulae for base spread resistance over! */
466 
467             /* ib term */
468 
469             EqualDeriv(&d_ib, &d_p);
470             d_ib.value = cb;
471             d_ib.d1_p = gbe/here->BJTtBetaF + gben;
472             d_ib.d2_p2 = gbe2/here->BJTtBetaF + gben2;
473             d_ib.d3_p3 = gbe3/here->BJTtBetaF + gben3;
474 
475             d_ib.d1_q = gbc/here->BJTtBetaR + gbcn;
476             d_ib.d2_q2 = gbc2/here->BJTtBetaR + gbcn2;
477             d_ib.d3_q3 = gbc3/here->BJTtBetaR + gbcn3;
478 
479             /* ib term over */
480             /*
481              *   charge storage elements
482              */
483             tf=here->BJTttransitTimeF;
484             tr=here->BJTttransitTimeR;
485             czbe=here->BJTtBEcap*here->BJTarea * here->BJTm;
486             pe=here->BJTtBEpot;
487             xme=here->BJTtjunctionExpBE;
488             cdis=model->BJTbaseFractionBCcap;
489             ctot=here->BJTtBCcap * here->BJTm;
490             if (model->BJTsubs == VERTICAL)
491                 ctot *= here->BJTareab;
492             else
493                 ctot *= here->BJTareac;
494             czbc=ctot*cdis;
495             czbx=ctot-czbc;
496             pc=here->BJTtBCpot;
497             xmc=here->BJTtjunctionExpBC;
498             fcpe=here->BJTtDepCap;
499             czcs=model->BJTcapSub * here->BJTm;
500             if (model->BJTsubs == VERTICAL)
501                 czcs *= here->BJTareac;
502             else
503                 czcs *= here->BJTareab;
504             ps=model->BJTpotentialSubstrate;
505             xms=model->BJTexponentialSubstrate;
506             xtf=model->BJTtransitTimeBiasCoeffF;
507             ovtf=model->BJTtransitTimeVBCFactor;
508             xjtf=here->BJTttransitTimeHighCurrentF*here->BJTarea * here->BJTm;
509             if(tf != 0 && vbe >0) {
510                 EqualDeriv(&d_cbe, &d_p);
511                 d_cbe.value = cbe;
512                 d_cbe.d1_p = gbe;
513                 d_cbe.d2_p2 = gbe2;
514                 d_cbe.d3_p3 = gbe3;
515                 if(xtf != 0){
516                     if(ovtf != 0) {
517                         /* dummy = exp ( vbc*ovtf) */
518                         EqualDeriv(&d_dummy, &d_q);
519                         d_dummy.value = vbc*ovtf;
520                         d_dummy.d1_q = ovtf;
521                         ExpDeriv(&d_dummy, &d_dummy);
522                     }
523                     else
524                     {
525                         EqualDeriv(&d_dummy,&d_p);
526                         d_dummy.value = 1.0;
527                         d_dummy.d1_p = 0.0;
528                     }
529                     if(xjtf != 0) {
530                         EqualDeriv(&d_dummy2, &d_cbe);
531                         d_dummy2.value += xjtf;
532                         DivDeriv(&d_dummy2, &d_cbe, &d_dummy2);
533                         MultDeriv (&d_dummy2, &d_dummy2, &d_dummy2);
534                     }
535                     else
536                     {
537                         EqualDeriv(&d_dummy2,&d_p);
538                         d_dummy2.value = 1.0;
539                         d_dummy2.d1_p = 0.0;
540                     }
541 
542                     MultDeriv(&d_tff, &d_dummy, &d_dummy2);
543                     TimesDeriv(&d_tff, &d_tff, tf*xtf);
544                     d_tff.value += tf;
545                 }
546                 else
547                 {
548                     EqualDeriv(&d_tff,&d_p);
549                     d_tff.value = tf;
550                     d_tff.d1_p = 0.0;
551                 }
552 
553 
554 
555 
556                 /* qbe = tff/qb*cbe */
557 
558                 /*
559                   dummy = tff/qb;
560                 */
561                 /* these are the cbe coeffs */
562                 DivDeriv(&d_dummy, &d_tff, &d_qb);
563                 MultDeriv(&d_qbe, &d_dummy, &d_cbe);
564 
565             }
566             else
567             {
568                 EqualDeriv(&d_qbe, &d_p);
569                 d_qbe.value = 0.0;
570                 d_qbe.d1_p = 0.0;
571             }
572             if (vbe < fcpe) {
573                 arg=1-vbe/pe;
574                 sarg=exp(-xme*log(arg));
575                 lcapbe1 = czbe*sarg;
576                 lcapbe2 =
577                     0.5*czbe*xme*sarg/(arg*pe);
578                 lcapbe3 =
579                     czbe*xme*(xme+1)*sarg/(arg*arg*pe*pe*6);
580             } else {
581                 f1=here->BJTtf1;
582                 f2=here->BJTtf2;
583                 f3=here->BJTtf3;
584                 czbef2=czbe/f2;
585                 lcapbe1 = czbef2*(f3+xme*vbe/pe);
586                 lcapbe2 = 0.5*xme*czbef2/pe;
587                 lcapbe3 = 0.0;
588             }
589             d_qbe.d1_p += lcapbe1;
590             d_qbe.d2_p2 += lcapbe2*2.;
591             d_qbe.d3_p3 += lcapbe3*6.;
592 
593 
594             fcpc=here->BJTtf4;
595             f1=here->BJTtf5;
596             f2=here->BJTtf6;
597             f3=here->BJTtf7;
598             if (vbc < fcpc) {
599                 arg=1-vbc/pc;
600                 sarg=exp(-xmc*log(arg));
601                 lcapbc1 = czbc*sarg;
602                 lcapbc2 =
603                     0.5*czbc*xmc*sarg/(arg*pc);
604                 lcapbc3 =
605                     czbc*xmc*(xmc+1)*sarg/(arg*arg*pc*pc*6);
606             } else {
607                 czbcf2=czbc/f2;
608                 lcapbc1 = czbcf2*(f3+xmc*vbc/pc);
609                 lcapbc2 = 0.5*xmc*czbcf2/pc;
610                 lcapbc3 = 0;
611             }
612             if(vbx < fcpc) {
613                 arg=1-vbx/pc;
614                 sarg=exp(-xmc*log(arg));
615                 lcapbx1 = czbx*sarg;
616                 lcapbx2 =
617                     0.5*czbx*xmc*sarg/(arg*pc);
618                 lcapbx3 =
619                     czbx*xmc*(xmc+1)*sarg/(arg*arg*pc*pc*6);
620             } else {
621                 czbxf2=czbx/f2;
622                 lcapbx1 = czbxf2*(f3+xmc*vbx/pc);
623                 lcapbx2 = 0.5*xmc*czbxf2/pc;
624                 lcapbx3 = 0;
625             }
626             if(vsc < 0){
627                 arg=1-vsc/ps;
628                 sarg=exp(-xms*log(arg));
629                 lcapsc1 = czcs*sarg;
630                 lcapsc2 =
631                     0.5*czcs*xms*sarg/(arg*ps);
632                 lcapsc3 =
633                     czcs*xms*(xms+1)*sarg/(arg*arg*ps*ps*6);
634             } else {
635                 lcapsc1 = czcs*(1+xms*vsc/ps);
636                 lcapsc2 = czcs*0.5*xms/ps;
637                 lcapsc3 = 0;
638             }
639 
640             /*
641              *   store small-signal parameters
642              */
643             here->ic_x = d_ic.d1_p;
644             here->ic_y = d_ic.d1_q;
645             here->ic_xd = d_ic.d1_r;
646             here->ic_x2 = 0.5*model->BJTtype*d_ic.d2_p2;
647             here->ic_y2 = 0.5*model->BJTtype*d_ic.d2_q2;
648             here->ic_w2 = 0.5*model->BJTtype*d_ic.d2_r2;
649             here->ic_xy = model->BJTtype*d_ic.d2_pq;
650             here->ic_yw = model->BJTtype*d_ic.d2_qr;
651             here->ic_xw = model->BJTtype*d_ic.d2_pr;
652             here->ic_x3 = d_ic.d3_p3/6.;
653             here->ic_y3 = d_ic.d3_q3/6.;
654             here->ic_w3 = d_ic.d3_r3/6.;
655             here->ic_x2w = 0.5*d_ic.d3_p2r;
656             here->ic_x2y = 0.5*d_ic.d3_p2q;
657             here->ic_y2w = 0.5*d_ic.d3_q2r;
658             here->ic_xy2 = 0.5*d_ic.d3_pq2;
659             here->ic_xw2 = 0.5*d_ic.d3_pr2;
660             here->ic_yw2 = 0.5*d_ic.d3_qr2;
661             here->ic_xyw = d_ic.d3_pqr;
662 
663             here->ib_x = d_ib.d1_p;
664             here->ib_y = d_ib.d1_q;
665             here->ib_x2 = 0.5*model->BJTtype*d_ib.d2_p2;
666             here->ib_y2 = 0.5*model->BJTtype*d_ib.d2_q2;
667             here->ib_xy = model->BJTtype*d_ib.d2_pq;
668             here->ib_x3 = d_ib.d3_p3/6.;
669             here->ib_y3 = d_ib.d3_q3/6.;
670             here->ib_x2y = 0.5*d_ib.d3_p2q;
671             here->ib_xy2 = 0.5*d_ib.d3_pq2;
672 
673             here->ibb_x = d_ibb.d1_p;
674             here->ibb_y = d_ibb.d1_q;
675             here->ibb_z = d_ibb.d1_r;
676             here->ibb_x2 = 0.5*model->BJTtype*d_ibb.d2_p2;
677             here->ibb_y2 = 0.5*model->BJTtype*d_ibb.d2_q2;
678             here->ibb_z2 = 0.5*model->BJTtype*d_ibb.d2_r2;
679             here->ibb_xy = model->BJTtype*d_ibb.d2_pq;
680             here->ibb_yz = model->BJTtype*d_ibb.d2_qr;
681             here->ibb_xz = model->BJTtype*d_ibb.d2_pr;
682             here->ibb_x3 = d_ibb.d3_p3/6.;
683             here->ibb_y3 = d_ibb.d3_q3/6.;
684             here->ibb_z3 = d_ibb.d3_r3/6.;
685             here->ibb_x2z = 0.5*d_ibb.d3_p2r;
686             here->ibb_x2y = 0.5*d_ibb.d3_p2q;
687             here->ibb_y2z = 0.5*d_ibb.d3_q2r;
688             here->ibb_xy2 = 0.5*d_ibb.d3_pq2;
689             here->ibb_xz2 = 0.5*d_ibb.d3_pr2;
690             here->ibb_yz2 = 0.5*d_ibb.d3_qr2;
691             here->ibb_xyz = d_ibb.d3_pqr;
692 
693             here->qbe_x = d_qbe.d1_p;
694             here->qbe_y = d_qbe.d1_q;
695             here->qbe_x2 = 0.5*model->BJTtype*d_qbe.d2_p2;
696             here->qbe_y2 = 0.5*model->BJTtype*d_qbe.d2_q2;
697             here->qbe_xy = model->BJTtype*d_qbe.d2_pq;
698             here->qbe_x3 = d_qbe.d3_p3/6.;
699             here->qbe_y3 = d_qbe.d3_q3/6.;
700             here->qbe_x2y = 0.5*d_qbe.d3_p2q;
701             here->qbe_xy2 = 0.5*d_qbe.d3_pq2;
702 
703             here->capbc1 = lcapbc1;
704             here->capbc2 = lcapbc2;
705             here->capbc3 = lcapbc3;
706 
707             here->capbx1 = lcapbx1;
708             here->capbx2 = lcapbx2;
709             here->capbx3 = lcapbx3;
710 
711             here->capsc1 = lcapsc1;
712             here->capsc2 = lcapsc2;
713             here->capsc3 = lcapsc3;
714         }
715     }
716     return(OK);
717 }
718