1 /**********
2 Copyright 1999 Regents of the University of California.  All rights reserved.
3 Author: Weidong Liu and Pin Su         Feb 1999
4 Author: 1998 Samuel Fung, Dennis Sinitsky and Stephen Tang
5 Modified by Pin Su, Wei Jin 99/9/27
6 Modified by Paolo Nenzi 2002
7 File: b3soifdld.c          98/5/01
8 **********/
9 
10 /*
11  * Revision 2.1  99/9/27 Pin Su
12  * BSIMFD2.1 release
13  */
14 
15 #include "ngspice/ngspice.h"
16 #include "ngspice/cktdefs.h"
17 #include "b3soifddef.h"
18 #include "ngspice/trandefs.h"
19 #include "ngspice/const.h"
20 #include "ngspice/sperror.h"
21 #include "ngspice/devdefs.h"
22 #include "ngspice/suffix.h"
23 
24 
25 #define MAX_EXP 5.834617425e14
26 #define MIN_EXP 1.713908431e-15
27 #define EXP_THRESHOLD 34.0
28 #define EPSOX 3.453133e-11
29 #define EPSSI 1.03594e-10
30 #define Charge_q 1.60219e-19
31 #define KboQ 8.617087e-5  /*  Kb / q   */
32 #define Eg300 1.115   /*  energy gap at 300K  */
33 #define DELTA_1 0.02
34 #define DELTA_2 0.02
35 #define DELTA_3 0.02
36 #define DELTA_4 0.02
37 #define DELT_Vbs0eff 0.02
38 #define DELT_Vbsmos  0.005
39 #define DELT_Vbseff  0.005
40 #define DELT_Xcsat   0.2
41 #define DELT_Vbs0dio 1e-7
42 #define DELTA_VFB  0.02
43 #define DELTA_Vcscv  0.0004
44 #define DELT_Vbsdio 0.01
45 #define CONST_2OV3 0.6666666666
46 #define OFF_Vbsdio  2e-2
47 #define OFF_Vbs0_dio 2.02e-2
48 #define QEX_FACT  20
49 
50 
51     /* B3SOIFDSmartVbs(Vbs, Old, here, check)
52      *  Smart Vbs guess.
53      */
54 
55 static double
B3SOIFDSmartVbs(double New,double Old,B3SOIFDinstance * here,CKTcircuit * ckt,int * check)56 B3SOIFDSmartVbs(double New, double Old, B3SOIFDinstance *here,
57                 CKTcircuit *ckt, int *check)
58 {
59    NG_IGNORE(Old);
60    NG_IGNORE(check);
61 
62    /* only do it for floating body and DC */
63    if (here->B3SOIFDfloat && (ckt->CKTmode & (MODEDC | MODEDCOP)))
64    {
65       /* Vbs cannot be negative in DC */
66       if (New < 0.0)  New = 0.0;
67    }
68    return(New);
69 }
70 
71 
72     /* B3SOIFDlimit(vnew,vold)
73      *  limits the per-iteration change of any absolute voltage value
74      */
75 
76 static double
B3SOIFDlimit(double vnew,double vold,double limit,int * check)77 B3SOIFDlimit(double vnew, double vold, double limit, int *check)
78 {
79     double T0, T1;
80 
81     if (isnan (vnew) || isnan (vold))
82     {
83 	fprintf(stderr, "Alberto says:  YOU TURKEY!  The limiting function received NaN.\n");
84 	fprintf(stderr, "New prediction returns to 0.0!\n");
85         vnew = 0.0;
86         *check = 1;
87     }
88     T0 = vnew - vold;
89     T1 = fabs(T0);
90     if (T1 > limit) {
91         if (T0 > 0.0)
92             vnew = vold + limit;
93         else
94             vnew = vold - limit;
95 	*check = 1;
96     }
97     return vnew;
98 }
99 
100 
101 
102 int
B3SOIFDload(GENmodel * inModel,CKTcircuit * ckt)103 B3SOIFDload(GENmodel *inModel, CKTcircuit *ckt)
104 {
105 B3SOIFDmodel *model = (B3SOIFDmodel*)inModel;
106 B3SOIFDinstance *here;
107 int selfheat;
108 
109 double ag0, qgd, qgs, von, cbhat, VgstNVt, ExpVgst = 0.0;
110 double cdhat, cdreq, ceqbd, ceqbs, ceqqb, ceqqd, ceqqg, ceq, geq;
111 double delvbd, delvbs, delvds, delvgd, delvgs;
112 double Vfbeff, dVfbeff_dVd, dVfbeff_dVb, V3, V4;
113 double gcgdb, gcggb, gcgsb, gcgeb, gcgT;
114 double gcsdb, gcsgb, gcssb, gcseb, gcsT;
115 double gcddb, gcdgb, gcdsb, gcdeb, gcdT;
116 double gcbdb, gcbgb, gcbsb, gcbeb, gcbT;
117 double gcedb, gcegb, gcesb, gceeb, gceT;
118 double gcTt, gTtg, gTtb, gTte, gTtdp, gTtt, gTtsp;
119 double vbd, vbs, vds, vgb, vgd, vgs, vgdo;
120 #ifndef PREDICTOR
121 double xfact;
122 #endif
123 double vg, vd, vs, vp, ve, vb;
124 double Vds, Vgs, Vbs, Gmbs, FwdSum, RevSum;
125 double Vgs_eff, Vfb, dVfb_dVb, dVfb_dVd, dVfb_dT;
126 double Phis, dPhis_dVb, sqrtPhis, dsqrtPhis_dVb, Vth = 0.0, dVth_dVb;
127 double dVth_dVd, dVth_dT;
128 double Vgst, dVgs_eff_dVg;
129 double n, dn_dVb, Vtm;
130 double ExpArg, V0;
131 double ueff = 0.0, dueff_dVg, dueff_dVd, dueff_dVb, dueff_dT;
132 double Esat, Vdsat = 0.0;
133 double EsatL, dEsatL_dVg, dEsatL_dVd, dEsatL_dVb, dEsatL_dT;
134 double dVdsat_dVg, dVdsat_dVb, dVdsat_dVd, dVdsat_dT, Vasat;
135 double dVasat_dVg, dVasat_dVb, dVasat_dVd, dVasat_dT;
136 double Va, dVa_dVd, dVa_dVg, dVa_dVb, dVa_dT;
137 double Vbseff, dVbseff_dVb;
138 double One_Third_CoxWL, Two_Third_CoxWL, CoxWL;
139 double T0, dT0_dVg, dT0_dVd, dT0_dVb, dT0_dVc, dT0_dVe, dT0_dT;
140 double T1, dT1_dVg, dT1_dVd, dT1_dVb, dT1_dVc, dT1_dVe, dT1_dT;
141 double T2, dT2_dVg, dT2_dVd, dT2_dVb, dT2_dVc, dT2_dVe, dT2_dT;
142 double T3, dT3_dVg, dT3_dVd, dT3_dVb, dT3_dVc, dT3_dVe, dT3_dT;
143 double T4, dT4_dVg, dT4_dVd, dT4_dVb, dT4_dVe, dT4_dT;
144 double T5, dT5_dVe;
145 double T6, dT6_dVe, dT6_dT;
146 double T7;
147 double T8;
148 double T9;
149 double T10;
150 double T11;
151 double Abulk, dAbulk_dVb, Abulk0, dAbulk0_dVb;
152 double VACLM, dVACLM_dVg, dVACLM_dVd, dVACLM_dVb, dVACLM_dT;
153 double VADIBL, dVADIBL_dVg, dVADIBL_dVd, dVADIBL_dVb, dVADIBL_dT;
154 double Xdep, dXdep_dVb, lt1, dlt1_dVb, ltw, dltw_dVb;
155 double Delt_vth, dDelt_vth_dVb, dDelt_vth_dT;
156 double Theta0, dTheta0_dVb;
157 double TempRatio, tmp1, tmp2, tmp3, tmp4;
158 double DIBL_Sft, dDIBL_Sft_dVd, Lambda, dLambda_dVg;
159 double a1;
160 
161 double Vgsteff = 0.0, dVgsteff_dVg, dVgsteff_dVd, dVgsteff_dVb;
162 double dVgsteff_dVe, dVgsteff_dT;
163 double Vdseff = 0.0, dVdseff_dVg, dVdseff_dVd, dVdseff_dVb, dVdseff_dT;
164 double VdseffCV, dVdseffCV_dVg, dVdseffCV_dVd, dVdseffCV_dVb;
165 double diffVds;
166 double dAbulk_dVg, dn_dVd ;
167 double beta, dbeta_dVg, dbeta_dVd, dbeta_dVb, dbeta_dT;
168 double gche, dgche_dVg, dgche_dVd, dgche_dVb, dgche_dT;
169 double fgche1, dfgche1_dVg, dfgche1_dVd, dfgche1_dVb, dfgche1_dT;
170 double fgche2, dfgche2_dVg, dfgche2_dVd, dfgche2_dVb, dfgche2_dT;
171 double Idl, dIdl_dVg, dIdl_dVd, dIdl_dVb, dIdl_dT;
172 double Ids = 0.0, Gm, Gds = 0.0, Gmb;
173 double CoxWovL;
174 double Rds, dRds_dVg, dRds_dVb, dRds_dT, WVCox, WVCoxRds;
175 double Vgst2Vtm, dVgst2Vtm_dT, VdsatCV, dVdsatCV_dVg, dVdsatCV_dVb;
176 double Leff, Weff, dWeff_dVg, dWeff_dVb;
177 double AbulkCV, dAbulkCV_dVb;
178 double qgdo, qgso, cgdo, cgso;
179 
180 double dxpart, sxpart;
181 
182 struct b3soifdSizeDependParam *pParam;
183 int ByPass, Check, ChargeComputationNeeded = 0, error;
184 
185 double gbbsp, gbbdp, gbbg, gbbb, gbbe, gbbp, gbbT;
186 double gddpsp, gddpdp, gddpg, gddpb, gddpe, gddpT;
187 double gsspsp, gsspdp, gsspg, gsspb, gsspe, gsspT;
188 double Gbpbs, Gbpgs, Gbpds, Gbpes, Gbpps, GbpT;
189 double ves, ved, veb, vge = 0.0, delves, vedo, delved;
190 double vps, vpd, Vps, delvps;
191 double Vbd, Ves, Vesfb, sqrtXdep, DeltVthtemp, dDeltVthtemp_dT;
192 double Vbp, dVbp_dVp, dVbp_dVb, dVbp_dVg, dVbp_dVd, dVbp_dVe, dVbp_dT;
193 double Vpsdio, dVpsdio_dVg, dVpsdio_dVd, dVpsdio_dVe, dVpsdio_dVp, dVpsdio_dT;
194 double DeltVthw, dDeltVthw_dVb, dDeltVthw_dT;
195 double dVbseff_dVd, dVbseff_dVe, dVbseff_dT;
196 double dVdsat_dVc, dVasat_dVc, dVACLM_dVc, dVADIBL_dVc, dVa_dVc;
197 double dfgche1_dVc, dfgche2_dVc, dgche_dVc, dVdseff_dVc, dIdl_dVc;
198 double Gm0, Gds0, Gmb0, GmT0, Gmc, Gme, GmT, dVbseff_dVg;
199 double dDIBL_Sft_dVb;
200 double diffVdsii  ;
201 double Idgidl = 0.0, Gdgidld, Gdgidlg, Isgidl = 0.0, Gsgidlg;
202 double Gjsd, Gjsb, GjsT, Gjdd, Gjdb, GjdT;
203 double Ibp = 0.0, Iii = 0.0, Giid, Giig, Giib, Giie, GiiT, Gcd, Gcb, GcT;
204 double ceqbody, ceqbodcon = 0.0;
205 double gppg = 0.0, gppdp = 0.0, gppb = 0.0, gppe = 0.0;
206 double gppp = 0.0, gppsp = 0.0, gppT;
207 double delTemp, deldelTemp, Temp;
208 double ceqth, ceqqth;
209 double K1;
210 double qjs = 0.0, gcjsbs, gcjsT;
211 double qjd = 0.0, gcjdbs, gcjdds, gcjdT;
212 double qge;
213 double ceqqe;
214 double ni, Eg, Cbox, Nfb, CboxWL;
215 double dVfbeff_dVrg, Cbe = 0.0;
216 double qinv = 0.0, qgate = 0.0, qbody = 0.0, qdrn = 0.0, qsrc, qsub = 0.0;
217 double cqgate, cqbody = 0.0, cqdrn = 0.0, cqsub, cqtemp;
218 double Cgg, Cgd, Cgb, Cge;
219 double Csg, Csd, Csb, Cse, Cbg = 0.0, Cbd = 0.0, Cbb = 0.0;
220 double Cgg1, Cgb1, Cgd1, Csg1, Csd1, Csb1;
221 double Vbs0t = 0.0, dVbs0t_dT ;
222 double Vbs0 = 0.0,dVbs0_dVe, dVbs0_dT;
223 double Vbs0eff = 0.0 ,dVbs0eff_dVg ,dVbs0eff_dVd ,dVbs0eff_dVe, dVbs0eff_dT;
224 double Vbs0teff = 0.0, dVbs0teff_dVg, dVbs0teff_dVd;
225 double dVbs0teff_dVe, dVbs0teff_dT;
226 double dVbsdio_dVg, dVbsdio_dVd, dVbsdio_dVe;
227 double dVbsdio_dVb, dVbsdio_dT;
228 double Vthfd = 0.0,dVthfd_dVd ,dVthfd_dVe, dVthfd_dT;
229 double Vbs0mos = 0.0 ,dVbs0mos_dVe, dVbs0mos_dT;
230 double Vbsmos ,dVbsmos_dVg ,dVbsmos_dVb ,dVbsmos_dVd, dVbsmos_dVe, dVbsmos_dT;
231 double Abeff ,dAbeff_dVg ,dAbeff_dVb, dAbeff_dVc;
232 double Vcs ,dVcs_dVg ,dVcs_dVb ,dVcs_dVd ,dVcs_dVe, dVcs_dT;
233 double Xcsat = 0.0 ,dXcsat_dVg , dXcsat_dVc;
234 double Vdsatii ,dVdsatii_dVg ,dVdsatii_dVd, dVdsatii_dVb, dVdsatii_dT;
235 double Vdseffii ,dVdseffii_dVg ,dVdseffii_dVd, dVdseffii_dVb, dVdseffii_dT;
236 double VcsCV = 0.0;
237 double VdsCV = 0.0, dVdsCV_dVg = 0.0, dVdsCV_dVb = 0.0;
238 double dVdsCV_dVd = 0.0, dVdsCV_dVc = 0.0;
239 double Phisd ,dPhisd_dVg ,dPhisd_dVb ,dPhisd_dVd,  dPhisd_dVc;
240 double sqrtPhisd;
241 double Xc = 0.0;
242 double Ic = 0.0;
243 double Ibs = 0.0;
244 double Ibd = 0.0;
245 double Denomi ,dDenomi_dVg ,dDenomi_dVd ,dDenomi_dVb , dDenomi_dT;
246 double Qbf = 0.0;
247 double Qsubs1 = 0.0;
248 double Qsubs2 = 0.0;
249 double Qsub0 = 0.0  ,dQsub0_dVg   ,dQsub0_dVb  ,dQsub0_dVd ;
250 double Qac0 = 0.0 ,dQac0_dVb   ,dQac0_dVd;
251 double Qdep0 ,dQdep0_dVb;
252 double Qe1 = 0.0;
253 double Ce1g ,Ce1b ,Ce1d ,Ce1e, Ce1T;
254 double Ce2g ,Ce2b ,Ce2d ,Ce2e, Ce2T;
255 double Qe2 = 0.0;
256 double dQac0_dVrg, Vbsdio = 0.0, dQsub0_dVrg;
257 
258 /*  for self-heating  */
259 double vbi, vfbb, phi, sqrtPhi, Xdep0, jbjt, jdif, jrec, jtun, u0temp, vsattemp;
260 double rds0, ua, ub, uc;
261 double dvbi_dT, dvfbb_dT, djbjt_dT, djdif_dT, djrec_dT, djtun_dT, du0temp_dT;
262 double dvsattemp_dT, drds0_dT, dua_dT, dub_dT, duc_dT, dni_dT, dVtm_dT;
263 double dVfbeff_dT, dQac0_dT, dQsub0_dT, dVdsCV_dT = 0.0, dPhisd_dT;
264 double CbT, CsT, CgT;
265 
266 double Qex, dQex_dVg, dQex_dVb, dQex_dVd, dQex_dVe, dQex_dT;
267 
268 /* clean up last */
269 FILE *fpdebug = NULL;
270 /* end clean up */
271 int nandetect;
272 static int nanfound = 0;
273 char nanmessage [12];
274 
275 double m;
276 
277 for (; model != NULL; model = B3SOIFDnextModel(model))
278 {    for (here = B3SOIFDinstances(model); here != NULL;
279           here = B3SOIFDnextInstance(here))
280      {
281           Check = 0;
282           ByPass = 0;
283           selfheat = (model->B3SOIFDshMod == 1) && (here->B3SOIFDrth0 != 0.0);
284 	  pParam = here->pParam;
285 
286           if (here->B3SOIFDdebugMod > 3)
287           {
288              if (model->B3SOIFDtype > 0)
289                 fpdebug = fopen("b3soifdn.log", "a");
290              else
291                 fpdebug = fopen("b3soifdp.log", "a");
292 
293              fprintf(fpdebug, "******* Time : %.5e ******* Device:  %s  Iteration:  %d\n",
294                      ckt->CKTtime, here->B3SOIFDname, here->B3SOIFDiterations);
295           }
296 
297           if ((ckt->CKTmode & MODEINITSMSIG))
298 	  {   vbs = *(ckt->CKTstate0 + here->B3SOIFDvbs);
299               vgs = *(ckt->CKTstate0 + here->B3SOIFDvgs);
300               ves = *(ckt->CKTstate0 + here->B3SOIFDves);
301               vps = *(ckt->CKTstate0 + here->B3SOIFDvps);
302               vds = *(ckt->CKTstate0 + here->B3SOIFDvds);
303               delTemp = *(ckt->CKTstate0 + here->B3SOIFDdeltemp);
304 
305               vg = *(ckt->CKTrhsOld + here->B3SOIFDgNode);
306               vd = *(ckt->CKTrhsOld + here->B3SOIFDdNodePrime);
307               vs = *(ckt->CKTrhsOld + here->B3SOIFDsNodePrime);
308               vp = *(ckt->CKTrhsOld + here->B3SOIFDpNode);
309               ve = *(ckt->CKTrhsOld + here->B3SOIFDeNode);
310               vb = *(ckt->CKTrhsOld + here->B3SOIFDbNode);
311 
312               if (here->B3SOIFDdebugMod > 2)
313               {
314                   fprintf(fpdebug, "... INIT SMSIG ...\n");
315               }
316               if (here->B3SOIFDdebugMod > 0)
317               {
318                  fprintf(stderr,
319                          "DC op. point converge with %d iterations\n",
320                          here->B3SOIFDiterations);
321               }
322           }
323 	  else if ((ckt->CKTmode & MODEINITTRAN))
324 	  {   vbs = *(ckt->CKTstate1 + here->B3SOIFDvbs);
325               vgs = *(ckt->CKTstate1 + here->B3SOIFDvgs);
326               ves = *(ckt->CKTstate1 + here->B3SOIFDves);
327               vps = *(ckt->CKTstate1 + here->B3SOIFDvps);
328               vds = *(ckt->CKTstate1 + here->B3SOIFDvds);
329               delTemp = *(ckt->CKTstate1 + here->B3SOIFDdeltemp);
330 
331               vg = *(ckt->CKTrhsOld + here->B3SOIFDgNode);
332               vd = *(ckt->CKTrhsOld + here->B3SOIFDdNodePrime);
333               vs = *(ckt->CKTrhsOld + here->B3SOIFDsNodePrime);
334               vp = *(ckt->CKTrhsOld + here->B3SOIFDpNode);
335               ve = *(ckt->CKTrhsOld + here->B3SOIFDeNode);
336               vb = *(ckt->CKTrhsOld + here->B3SOIFDbNode);
337 
338               if (here->B3SOIFDdebugMod > 2)
339               {
340                  fprintf(fpdebug, "... Init Transient ....\n");
341               }
342               if (here->B3SOIFDdebugMod > 0)
343               {
344                  fprintf(stderr, "Transient operation point converge with %d iterations\n",
345 here->B3SOIFDiterations);
346               }
347               here->B3SOIFDiterations = 0;
348           }
349 	  else if ((ckt->CKTmode & MODEINITJCT) && !here->B3SOIFDoff)
350 	  {   vds = model->B3SOIFDtype * here->B3SOIFDicVDS;
351               vgs = model->B3SOIFDtype * here->B3SOIFDicVGS;
352               ves = model->B3SOIFDtype * here->B3SOIFDicVES;
353               vbs = model->B3SOIFDtype * here->B3SOIFDicVBS;
354               vps = model->B3SOIFDtype * here->B3SOIFDicVPS;
355 
356 	      vg = vd = vs = vp = ve = 0.0;
357 
358               here->B3SOIFDiterations = 0;  /*  initialize iteration number  */
359 
360               delTemp = 0.0;
361 	      here->B3SOIFDphi = pParam->B3SOIFDphi;
362 
363 
364               if (here->B3SOIFDdebugMod > 2)
365 		fprintf(fpdebug, "... INIT JCT ...\n");
366 
367 	      if ((vds == 0.0) && (vgs == 0.0) && (vbs == 0.0) &&
368 	         ((ckt->CKTmode & (MODETRAN | MODEAC|MODEDCOP |
369 		   MODEDCTRANCURVE)) || (!(ckt->CKTmode & MODEUIC))))
370 	      {   vbs = 0.0;
371 		  vgs = model->B3SOIFDtype*0.1 + pParam->B3SOIFDvth0;
372 		  vds = 0.0;
373 		  ves = 0.0;
374 		  vps = 0.0;
375 	      }
376 	  }
377 	  else if ((ckt->CKTmode & (MODEINITJCT | MODEINITFIX)) &&
378 		  (here->B3SOIFDoff))
379 	  {    delTemp = vps = vbs = vgs = vds = ves = 0.0;
380                vg = vd = vs = vp = ve = 0.0;
381                here->B3SOIFDiterations = 0;  /*  initialize iteration number  */
382 	  }
383 	  else
384 	  {
385 #ifndef PREDICTOR
386 	       if ((ckt->CKTmode & MODEINITPRED))
387 	       {   xfact = ckt->CKTdelta / ckt->CKTdeltaOld[1];
388 		   *(ckt->CKTstate0 + here->B3SOIFDvbs) =
389 			 *(ckt->CKTstate1 + here->B3SOIFDvbs);
390 		   vbs = (1.0 + xfact)* (*(ckt->CKTstate1 + here->B3SOIFDvbs))
391 			 - (xfact * (*(ckt->CKTstate2 + here->B3SOIFDvbs)));
392 		   *(ckt->CKTstate0 + here->B3SOIFDvgs) =
393 			 *(ckt->CKTstate1 + here->B3SOIFDvgs);
394 		   vgs = (1.0 + xfact)* (*(ckt->CKTstate1 + here->B3SOIFDvgs))
395 			 - (xfact * (*(ckt->CKTstate2 + here->B3SOIFDvgs)));
396 		   *(ckt->CKTstate0 + here->B3SOIFDves) =
397 			 *(ckt->CKTstate1 + here->B3SOIFDves);
398 		   ves = (1.0 + xfact)* (*(ckt->CKTstate1 + here->B3SOIFDves))
399 			 - (xfact * (*(ckt->CKTstate2 + here->B3SOIFDves)));
400 		   *(ckt->CKTstate0 + here->B3SOIFDvps) =
401 			 *(ckt->CKTstate1 + here->B3SOIFDvps);
402 		   vps = (1.0 + xfact)* (*(ckt->CKTstate1 + here->B3SOIFDvps))
403 			 - (xfact * (*(ckt->CKTstate2 + here->B3SOIFDvps)));
404 		   *(ckt->CKTstate0 + here->B3SOIFDvds) =
405 			 *(ckt->CKTstate1 + here->B3SOIFDvds);
406 		   vds = (1.0 + xfact)* (*(ckt->CKTstate1 + here->B3SOIFDvds))
407 			 - (xfact * (*(ckt->CKTstate2 + here->B3SOIFDvds)));
408 		   *(ckt->CKTstate0 + here->B3SOIFDvbd) =
409 			 *(ckt->CKTstate0 + here->B3SOIFDvbs)
410 			 - *(ckt->CKTstate0 + here->B3SOIFDvds);
411 
412                    *(ckt->CKTstate0 + here->B3SOIFDvg) = *(ckt->CKTstate1 + here->B3SOIFDvg);
413                    *(ckt->CKTstate0 + here->B3SOIFDvd) = *(ckt->CKTstate1 + here->B3SOIFDvd);
414                    *(ckt->CKTstate0 + here->B3SOIFDvs) = *(ckt->CKTstate1 + here->B3SOIFDvs);
415                    *(ckt->CKTstate0 + here->B3SOIFDvp) = *(ckt->CKTstate1 + here->B3SOIFDvp);
416                    *(ckt->CKTstate0 + here->B3SOIFDve) = *(ckt->CKTstate1 + here->B3SOIFDve);
417 
418                    /* Only predict ve */
419                    ve = (1.0 + xfact)* (*(ckt->CKTstate1 + here->B3SOIFDve))
420 
421                         - (xfact * (*(ckt->CKTstate2 + here->B3SOIFDve)));
422                    /* Then update vg, vs, vb, vd, vp base on ve */
423                    vs = ve - model->B3SOIFDtype * ves;
424                    vg = model->B3SOIFDtype * vgs + vs;
425                    vd = model->B3SOIFDtype * vds + vs;
426                    vb = model->B3SOIFDtype * vbs + vs;
427                    vp = model->B3SOIFDtype * vps + vs;
428 
429 		   delTemp = (1.0 + xfact)* (*(ckt->CKTstate1 +
430 			 here->B3SOIFDdeltemp))-(xfact * (*(ckt->CKTstate2 +
431 			 here->B3SOIFDdeltemp)));
432 
433 		   if (selfheat)
434 		   {
435 		       here->B3SOIFDphi = 2.0 * here->B3SOIFDvtm
436 					* log (pParam->B3SOIFDnpeak /
437 					       here->B3SOIFDni);
438 		   }
439 
440 		   if (here->B3SOIFDdebugMod > 0)
441 		   {
442                       fprintf(stderr, "Time = %.6e converge with %d iterations\n", ckt->CKTtime, here->B3SOIFDiterations);
443                    }
444 		   if (here->B3SOIFDdebugMod > 2)
445 		   {
446 		      fprintf(fpdebug, "... PREDICTOR calculation ....\n");
447 		   }
448                    here->B3SOIFDiterations = 0;
449 	       }
450 	       else
451 	       {
452 #endif /* PREDICTOR */
453 
454                    vg = B3SOIFDlimit(*(ckt->CKTrhsOld + here->B3SOIFDgNode),
455                                  *(ckt->CKTstate0 + here->B3SOIFDvg), 3.0, &Check);
456                    vd = B3SOIFDlimit(*(ckt->CKTrhsOld + here->B3SOIFDdNodePrime),
457                                  *(ckt->CKTstate0 + here->B3SOIFDvd), 3.0, &Check);
458                    vs = B3SOIFDlimit(*(ckt->CKTrhsOld + here->B3SOIFDsNodePrime),
459                                  *(ckt->CKTstate0 + here->B3SOIFDvs), 3.0, &Check);
460                    vp = B3SOIFDlimit(*(ckt->CKTrhsOld + here->B3SOIFDpNode),
461                                  *(ckt->CKTstate0 + here->B3SOIFDvp), 3.0, &Check);
462                    ve = B3SOIFDlimit(*(ckt->CKTrhsOld + here->B3SOIFDeNode),
463                                  *(ckt->CKTstate0 + here->B3SOIFDve), 3.0, &Check);
464                    delTemp = *(ckt->CKTrhsOld + here->B3SOIFDtempNode);
465 
466 		   vbs = model->B3SOIFDtype * (*(ckt->CKTrhsOld+here->B3SOIFDbNode)
467                                 - *(ckt->CKTrhsOld+here->B3SOIFDsNodePrime));
468 
469 		   vps = model->B3SOIFDtype * (vp - vs);
470 		   vgs = model->B3SOIFDtype * (vg - vs);
471 		   ves = model->B3SOIFDtype * (ve - vs);
472 		   vds = model->B3SOIFDtype * (vd - vs);
473 
474 		   if (here->B3SOIFDdebugMod > 2)
475 		   {
476 		      fprintf(fpdebug, "... DC calculation ....\n");
477 fprintf(fpdebug, "Vg = %.10f; Vb = %.10f; Vs = %.10f\n",
478 			 *(ckt->CKTrhsOld + here->B3SOIFDgNode),
479 			 *(ckt->CKTrhsOld + here->B3SOIFDbNode),
480 			 *(ckt->CKTrhsOld + here->B3SOIFDsNode));
481 fprintf(fpdebug, "Vd = %.10f; Vsp = %.10f; Vdp = %.10f\n",
482 			 *(ckt->CKTrhsOld + here->B3SOIFDdNode),
483 			 *(ckt->CKTrhsOld + here->B3SOIFDsNodePrime),
484 			 *(ckt->CKTrhsOld + here->B3SOIFDdNodePrime));
485 fprintf(fpdebug, "Ve = %.10f; Vp = %.10f; delTemp = %.10f\n",
486 			 *(ckt->CKTrhsOld + here->B3SOIFDeNode),
487 			 *(ckt->CKTrhsOld + here->B3SOIFDpNode),
488 			 *(ckt->CKTrhsOld + here->B3SOIFDtempNode));
489 
490 		   }
491 
492 #ifndef PREDICTOR
493 	       }
494 #endif /* PREDICTOR */
495 
496 	       vbd = vbs - vds;
497 	       vgd = vgs - vds;
498                ved = ves - vds;
499 	       vgdo = *(ckt->CKTstate0 + here->B3SOIFDvgs)
500 		    - *(ckt->CKTstate0 + here->B3SOIFDvds);
501 	       vedo = *(ckt->CKTstate0 + here->B3SOIFDves)
502 		    - *(ckt->CKTstate0 + here->B3SOIFDvds);
503 	       delvbs = vbs - *(ckt->CKTstate0 + here->B3SOIFDvbs);
504 	       delvbd = vbd - *(ckt->CKTstate0 + here->B3SOIFDvbd);
505 	       delvgs = vgs - *(ckt->CKTstate0 + here->B3SOIFDvgs);
506 	       delves = ves - *(ckt->CKTstate0 + here->B3SOIFDves);
507 	       delvps = vps - *(ckt->CKTstate0 + here->B3SOIFDvps);
508 	       deldelTemp = delTemp - *(ckt->CKTstate0 + here->B3SOIFDdeltemp);
509 	       delvds = vds - *(ckt->CKTstate0 + here->B3SOIFDvds);
510 	       delvgd = vgd - vgdo;
511                delved = ved - vedo;
512 
513 	       if (here->B3SOIFDmode >= 0)
514 	       {
515                    cdhat = here->B3SOIFDcd + (here->B3SOIFDgm-here->B3SOIFDgjdg) * delvgs
516                          + (here->B3SOIFDgds - here->B3SOIFDgjdd) * delvds
517                          + (here->B3SOIFDgmbs - here->B3SOIFDgjdb) * delvbs
518                          + (here->B3SOIFDgme - here->B3SOIFDgjde) * delves
519 			 + (here->B3SOIFDgmT - here->B3SOIFDgjdT) * deldelTemp;
520 	       }
521 	       else
522 	       {
523                    cdhat = here->B3SOIFDcd + (here->B3SOIFDgm-here->B3SOIFDgjdg) * delvgd
524                          - (here->B3SOIFDgds - here->B3SOIFDgjdd) * delvds
525                          + (here->B3SOIFDgmbs - here->B3SOIFDgjdb) * delvbd
526                          + (here->B3SOIFDgme - here->B3SOIFDgjde) * delved
527                          + (here->B3SOIFDgmT - here->B3SOIFDgjdT) * deldelTemp;
528 
529 	       }
530 	       cbhat = here->B3SOIFDcb + here->B3SOIFDgbgs * delvgs
531 		     + here->B3SOIFDgbbs * delvbs + here->B3SOIFDgbds * delvds
532                      + here->B3SOIFDgbes * delves + here->B3SOIFDgbps * delvps
533 		     + here->B3SOIFDgbT * deldelTemp;
534 
535 #ifndef NOBYPASS
536 	   /* following should be one big if connected by && all over
537 	    * the place, but some C compilers can't handle that, so
538 	    * we split it up here to let them digest it in stages
539 	    */
540 
541 	       if (here->B3SOIFDdebugMod > 3)
542                {
543 fprintf(fpdebug, "Convergent Criteria : vbs %d  vds %d  vgs %d  ves %d  vps %d  temp %d\n",
544 	((fabs(delvbs) < (ckt->CKTreltol * MAX(fabs(vbs),
545 	fabs(*(ckt->CKTstate0+here->B3SOIFDvbs))) + ckt->CKTvoltTol))) ? 1 : 0,
546 	((fabs(delvds) < (ckt->CKTreltol * MAX(fabs(vds),
547 	fabs(*(ckt->CKTstate0+here->B3SOIFDvds))) + ckt->CKTvoltTol))) ? 1 : 0,
548 	((fabs(delvgs) < (ckt->CKTreltol * MAX(fabs(vgs),
549 	fabs(*(ckt->CKTstate0+here->B3SOIFDvgs))) + ckt->CKTvoltTol))) ? 1 : 0,
550 	((fabs(delves) < (ckt->CKTreltol * MAX(fabs(ves),
551 	fabs(*(ckt->CKTstate0+here->B3SOIFDves))) + ckt->CKTvoltTol))) ? 1 : 0,
552 	((fabs(delvps) < (ckt->CKTreltol * MAX(fabs(vps),
553 	fabs(*(ckt->CKTstate0+here->B3SOIFDvps))) + ckt->CKTvoltTol))) ? 1 : 0,
554 	((fabs(deldelTemp) < (ckt->CKTreltol * MAX(fabs(delTemp),
555 	fabs(*(ckt->CKTstate0+here->B3SOIFDdeltemp))) + ckt->CKTvoltTol*1e4))) ? 1 : 0);
556 fprintf(fpdebug, "delCd %.4e, delCb %.4e\n",  fabs(cdhat - here->B3SOIFDcd) ,
557         fabs(cbhat - here->B3SOIFDcb));
558 
559                }
560 	       if ((!(ckt->CKTmode & MODEINITPRED)) && (ckt->CKTbypass) && Check == 0)
561 	       if ((fabs(delvbs) < (ckt->CKTreltol * MAX(fabs(vbs),
562 		   fabs(*(ckt->CKTstate0+here->B3SOIFDvbs))) + ckt->CKTvoltTol))  )
563 	       if ((fabs(delvbd) < (ckt->CKTreltol * MAX(fabs(vbd),
564 		   fabs(*(ckt->CKTstate0+here->B3SOIFDvbd))) + ckt->CKTvoltTol))  )
565 	       if ((fabs(delvgs) < (ckt->CKTreltol * MAX(fabs(vgs),
566 		   fabs(*(ckt->CKTstate0+here->B3SOIFDvgs))) + ckt->CKTvoltTol)))
567 	       if ((fabs(delves) < (ckt->CKTreltol * MAX(fabs(ves),
568 		   fabs(*(ckt->CKTstate0+here->B3SOIFDves))) + ckt->CKTvoltTol)))
569 	       if ( (here->B3SOIFDbodyMod == 0) || (here->B3SOIFDbodyMod == 2) ||
570                   (fabs(delvps) < (ckt->CKTreltol * MAX(fabs(vps),
571 		   fabs(*(ckt->CKTstate0+here->B3SOIFDvps))) + ckt->CKTvoltTol)) )
572 	       if ( (here->B3SOIFDtempNode == 0)  ||
573                   (fabs(deldelTemp) < (ckt->CKTreltol * MAX(fabs(delTemp),
574 		   fabs(*(ckt->CKTstate0+here->B3SOIFDdeltemp)))
575 		   + ckt->CKTvoltTol*1e4)))
576 	       if ((fabs(delvds) < (ckt->CKTreltol * MAX(fabs(vds),
577 		   fabs(*(ckt->CKTstate0+here->B3SOIFDvds))) + ckt->CKTvoltTol)))
578 	       if ((fabs(cdhat - here->B3SOIFDcd) < ckt->CKTreltol
579 		   * MAX(fabs(cdhat),fabs(here->B3SOIFDcd)) + ckt->CKTabstol))
580 	       if ((fabs(cbhat - here->B3SOIFDcb) < ckt->CKTreltol
581 		   * MAX(fabs(cbhat),fabs(here->B3SOIFDcb)) + ckt->CKTabstol) )
582 	       {   /* bypass code */
583 	           vbs = *(ckt->CKTstate0 + here->B3SOIFDvbs);
584 	           vbd = *(ckt->CKTstate0 + here->B3SOIFDvbd);
585 	           vgs = *(ckt->CKTstate0 + here->B3SOIFDvgs);
586 	           ves = *(ckt->CKTstate0 + here->B3SOIFDves);
587 	           vps = *(ckt->CKTstate0 + here->B3SOIFDvps);
588 	           vds = *(ckt->CKTstate0 + here->B3SOIFDvds);
589 	           delTemp = *(ckt->CKTstate0 + here->B3SOIFDdeltemp);
590 
591 		   /*  calculate Vds for temperature conductance calculation
592 		       in bypass (used later when filling Temp node matrix)  */
593 		   Vds = here->B3SOIFDmode > 0 ? vds : -vds;
594 
595 	           vgd = vgs - vds;
596 	           vgb = vgs - vbs;
597                    veb = ves - vbs;
598 
599 	           if (here->B3SOIFDdebugMod > 2)
600 	           {
601 fprintf(stderr, "Bypass for %s...\n", here->B3SOIFDname);
602 	    	      fprintf(fpdebug, "... By pass  ....\n");
603 		      fprintf(fpdebug, "vgs=%.4f, vds=%.4f, vbs=%.4f, ",
604 			   vgs, vds, vbs);
605 	    	      fprintf(fpdebug, "ves=%.4f, vps=%.4f\n", ves, vps);
606 		   }
607 		   if ((ckt->CKTmode & (MODETRAN | MODEAC)) ||
608 		      ((ckt->CKTmode & MODETRANOP) && (ckt->CKTmode & MODEUIC)))
609 		   {   ByPass = 1;
610 		       goto line755;
611 		   }
612 		   else
613 		   {   goto line850;
614 		   }
615 	       }
616 
617 
618 #endif /*NOBYPASS*/
619 		       von = here->B3SOIFDvon;
620 
621 			if ((here->B3SOIFDdebugMod > 1) || (here->B3SOIFDdebugMod == -1))
622 		       {
623 			   here->B3SOIFDdum1 = here->B3SOIFDdum2 = here->B3SOIFDdum3 = 0.0;
624                            here->B3SOIFDdum4 = here->B3SOIFDdum5 = 0.0;
625 			   Qac0 = Qsub0 = Qsubs1 = Qsubs2 = Qbf = Qe1 = Qe2 = 0.0;
626 			   qjs = qjd = Cbg = Cbb = Cbd = Cbe = Xc = qdrn = qgate = 0.0;
627 			   qbody = qsub = 0.0;
628 		       }
629 
630 		       if (here->B3SOIFDdebugMod > 2) {
631 			      fprintf(fpdebug, "Limited : vgs = %.8f\n", vgs);
632 			      fprintf(fpdebug, "Limited : vds = %.8f\n", vds);
633 		       }
634 
635                        if (*(ckt->CKTstate0 + here->B3SOIFDvds) >= 0.0)
636                           T0 = *(ckt->CKTstate0 + here->B3SOIFDvbs);
637                        else
638                           T0 = *(ckt->CKTstate0 + here->B3SOIFDvbd);
639 
640 		       if (here->B3SOIFDdebugMod > 2)
641 			  fprintf(fpdebug, "Before lim : vbs = %.8f, after = ", T0);
642 
643 		       if (vds >= 0.0)
644 		       {
645 		           vbs = B3SOIFDlimit(vbs, T0, 0.2, &Check);
646                            vbs = B3SOIFDSmartVbs(vbs, T0, here, ckt, &Check);
647 			   vbd = vbs - vds;
648                            vb = model->B3SOIFDtype * vbs + vs;
649 		           if (here->B3SOIFDdebugMod > 2)
650 			      fprintf(fpdebug, "%.8f\n", vbs);
651 		       } else
652 		       {
653 		           vbd = B3SOIFDlimit(vbd, T0, 0.2, &Check);
654                            vbd = B3SOIFDSmartVbs(vbd, T0, here, ckt, &Check);
655 			   vbs = vbd + vds;
656                            vb = model->B3SOIFDtype * vbs + vd;
657 		           if (here->B3SOIFDdebugMod > 2)
658 			      fprintf(fpdebug, "%.8f\n", vbd);
659 		       }
660 
661 		       delTemp =B3SOIFDlimit(delTemp, *(ckt->CKTstate0 + here->B3SOIFDdeltemp),5.0,&Check);
662 
663                   }
664 
665 /*  Calculate temperature dependent values for self-heating effect  */
666 		  Temp = delTemp + ckt->CKTtemp;
667 /* for debugging
668   Temp = ckt->CKTtemp;
669   selfheat = 1;
670   if (here->B3SOIFDname[1] == '2')
671   {
672      Temp += 0.01;
673   } */
674 		  TempRatio = Temp / model->B3SOIFDtnom;
675 
676 		  if (selfheat) {
677 		      Vtm = KboQ * Temp;
678 
679                       T0 = 1108.0 + Temp;
680 		      T5 = Temp * Temp;
681 		      Eg = 1.16 - 7.02e-4 * T5 / T0;
682 		      T1 = ((7.02e-4 * T5) - T0 * (14.04e-4 * Temp)) / T0 / T0;
683                       /*  T1 = dEg / dT  */
684 
685                       T2 = 1.9230584e-4;  /*  T2 = 1 / 300.15^(3/2)  */
686 		      T5 = sqrt(Temp);
687                       T3 = 1.45e10 * Temp * T5 * T2;
688                       T4 = exp(21.5565981 - Eg / (2.0 * Vtm));
689 		      ni = T3 * T4;
690                       dni_dT = 2.175e10 * T2 * T5 * T4 + T3 * T4 *
691                                (-Vtm * T1 + Eg * KboQ) / (2.0 * Vtm * Vtm);
692 
693                       T0 = log(1.0e20 * pParam->B3SOIFDnpeak / (ni * ni));
694 		      vbi = Vtm * T0;
695                       dvbi_dT = KboQ * T0 + Vtm * (-2.0 * dni_dT / ni);
696 
697 		      if (pParam->B3SOIFDnsub > 0) {
698                          T0 = log(pParam->B3SOIFDnpeak / pParam->B3SOIFDnsub);
699 		         vfbb = -model->B3SOIFDtype * Vtm*T0;
700                          dvfbb_dT = -model->B3SOIFDtype * KboQ*T0;
701                       }
702 		      else {
703                          T0 = log(-pParam->B3SOIFDnpeak*pParam->B3SOIFDnsub/ni/ni);
704 		         vfbb = -model->B3SOIFDtype * Vtm*T0;
705                          dvfbb_dT = -model->B3SOIFDtype *
706                                    (KboQ * T0 + Vtm * 2.0 * dni_dT / ni);
707                       }
708 
709 /*		      phi = 2.0 * Vtm * log(pParam->B3SOIFDnpeak / ni);  */
710 		      phi = here->B3SOIFDphi;
711 		      sqrtPhi = sqrt(phi);
712 		      Xdep0 = sqrt(2.0 * EPSSI / (Charge_q
713 				         * pParam->B3SOIFDnpeak * 1.0e6))
714 				         * sqrtPhi;
715 		      /*  Save the values below for phi calculation in B3SOIFDaccept()  */
716 		      here->B3SOIFDvtm = Vtm;
717 		      here->B3SOIFDni = ni;
718 
719                       /*  Use dTx_dVe variables to act as dTx_dT variables  */
720 
721                       T8 = 1 / model->B3SOIFDtnom;
722                       T7 = model->B3SOIFDxbjt / pParam->B3SOIFDndiode;
723 		      T0 = pow(TempRatio, T7);
724                       dT0_dVe = T7 * pow(TempRatio, T7 - 1.0) * T8;
725 
726                       T7 = model->B3SOIFDxdif / pParam->B3SOIFDndiode;
727 		      T1 = pow(TempRatio, T7);
728                       dT1_dVe = T7 * pow(TempRatio, T7 - 1.0) * T8;
729 
730                       T7 = model->B3SOIFDxrec / pParam->B3SOIFDndiode / 2.0;
731 		      T2 = pow(TempRatio, T7);
732                       dT2_dVe = T7 * pow(TempRatio, T7 - 1.0) * T8;
733 
734 		      T3 = TempRatio - 1.0;
735 		      T4 = Eg300 / pParam->B3SOIFDndiode / Vtm * T3;
736                       dT4_dVe = Eg300 / pParam->B3SOIFDndiode / Vtm / Vtm *
737                                 (Vtm * T8 - T3 * KboQ);
738 		      T5 = exp(T4);
739                       dT5_dVe = dT4_dVe * T5;
740 		      T6 = sqrt(T5);
741                       dT6_dVe = 0.5 / T6 * dT5_dVe;
742 
743 		      jbjt = pParam->B3SOIFDisbjt * T0 * T5;
744 		      jdif = pParam->B3SOIFDisdif * T1 * T5;
745 		      jrec = pParam->B3SOIFDisrec * T2 * T6;
746                       djbjt_dT = pParam->B3SOIFDisbjt * (T0 * dT5_dVe + T5 * dT0_dVe);
747                       djdif_dT = pParam->B3SOIFDisdif * (T1 * dT5_dVe + T5 * dT1_dVe);
748                       djrec_dT = pParam->B3SOIFDisrec * (T2 * dT6_dVe + T6 * dT2_dVe);
749 
750                       T7 = model->B3SOIFDxtun / pParam->B3SOIFDntun;
751 		      T0 = pow(TempRatio, T7);
752 		      jtun = pParam->B3SOIFDistun * T0;
753                       djtun_dT = pParam->B3SOIFDistun * T7 * pow(TempRatio, T7 - 1.0) * T8;
754 
755 		      u0temp = pParam->B3SOIFDu0 * pow(TempRatio, pParam->B3SOIFDute);
756                       du0temp_dT = pParam->B3SOIFDu0 * pParam->B3SOIFDute *
757                                    pow(TempRatio, pParam->B3SOIFDute - 1.0) * T8;
758 
759 		      vsattemp = pParam->B3SOIFDvsat - pParam->B3SOIFDat * T3;
760                       dvsattemp_dT = -pParam->B3SOIFDat * T8;
761 
762 		      rds0 = (pParam->B3SOIFDrdsw + pParam->B3SOIFDprt
763 		          * T3) / pParam->B3SOIFDrds0denom;
764                       drds0_dT = pParam->B3SOIFDprt / pParam->B3SOIFDrds0denom * T8;
765 
766 		      ua = pParam->B3SOIFDuatemp + pParam->B3SOIFDua1 * T3;
767 		      ub = pParam->B3SOIFDubtemp + pParam->B3SOIFDub1 * T3;
768 		      uc = pParam->B3SOIFDuctemp + pParam->B3SOIFDuc1 * T3;
769                       dua_dT = pParam->B3SOIFDua1 * T8;
770                       dub_dT = pParam->B3SOIFDub1 * T8;
771                       duc_dT = pParam->B3SOIFDuc1 * T8;
772 		  }
773 		  else {
774                       vbi = pParam->B3SOIFDvbi;
775                       vfbb = pParam->B3SOIFDvfbb;
776                       phi = pParam->B3SOIFDphi;
777                       sqrtPhi = pParam->B3SOIFDsqrtPhi;
778                       Xdep0 = pParam->B3SOIFDXdep0;
779                       jbjt = pParam->B3SOIFDjbjt;
780                       jdif = pParam->B3SOIFDjdif;
781                       jrec = pParam->B3SOIFDjrec;
782                       jtun = pParam->B3SOIFDjtun;
783                       u0temp = pParam->B3SOIFDu0temp;
784                       vsattemp = pParam->B3SOIFDvsattemp;
785                       rds0 = pParam->B3SOIFDrds0;
786                       ua = pParam->B3SOIFDua;
787                       ub = pParam->B3SOIFDub;
788                       uc = pParam->B3SOIFDuc;
789                       dni_dT = dvbi_dT = dvfbb_dT = djbjt_dT = djdif_dT = 0.0;
790                       djrec_dT = djtun_dT = du0temp_dT = dvsattemp_dT = 0.0;
791                       drds0_dT = dua_dT = dub_dT = duc_dT = 0.0;
792 		  }
793 
794 		  /* TempRatio used for Vth and mobility */
795 		  if (selfheat) {
796 		      TempRatio = Temp / model->B3SOIFDtnom - 1.0;
797 		  }
798 		  else {
799 		      TempRatio =  ckt->CKTtemp / model->B3SOIFDtnom - 1.0;
800 		  }
801 
802 		  /* determine DC current and derivatives */
803 		  vbd = vbs - vds;
804 		  vgd = vgs - vds;
805 		  vgb = vgs - vbs;
806 		  ved = ves - vds;
807 		  veb = ves - vbs;
808 		  vge = vgs - ves;
809 		  vpd = vps - vds;
810 
811 
812 		  if (vds >= 0.0)
813 		  {   /* normal mode */
814 		      here->B3SOIFDmode = 1;
815 		      Vds = vds;
816 		      Vgs = vgs;
817 		      Vbs = vbs;
818 		      Vbd = vbd;
819 		      Ves = ves;
820 		      Vps = vps;
821 		  }
822 		  else
823 		  {   /* inverse mode */
824 		      here->B3SOIFDmode = -1;
825 		      Vds = -vds;
826 		      Vgs = vgd;
827 		      Vbs = vbd;
828 		      Vbd = vbs;
829 		      Ves = ved;
830 		      Vps = vpd;
831 		  }
832 
833 
834                   if (here->B3SOIFDdebugMod > 2)
835 		  {
836 		     fprintf(fpdebug, "Vgs=%.4f, Vds=%.4f, Vbs=%.4f, ",
837 			   Vgs, Vds, Vbs);
838 		     fprintf(fpdebug, "Ves=%.4f, Vps=%.4f, Temp=%.1f\n",
839 			   Ves, Vps, Temp);
840 		  }
841 
842 		  Vesfb = Ves - vfbb;
843 		  Cbox = model->B3SOIFDcbox;
844 		  K1 = pParam->B3SOIFDk1;
845 
846 		  ChargeComputationNeeded =
847 			 ((ckt->CKTmode & (MODEAC | MODETRAN | MODEINITSMSIG)) ||
848 			 ((ckt->CKTmode & MODETRANOP) && (ckt->CKTmode & MODEUIC)))
849 			 ? 1 : 0;
850 
851                   if (here->B3SOIFDdebugMod == -1)
852                      ChargeComputationNeeded = 1;
853 
854 
855 
856 
857 /* Poly Gate Si Depletion Effect */
858 		  T0 = pParam->B3SOIFDvfb + phi;
859 		  if ((pParam->B3SOIFDngate > 1.e18) && (pParam->B3SOIFDngate < 1.e25)
860 		       && (Vgs > T0))
861 		  /* added to avoid the problem caused by ngate */
862 		  {   T1 = 1.0e6 * Charge_q * EPSSI * pParam->B3SOIFDngate
863 			 / (model->B3SOIFDcox * model->B3SOIFDcox);
864 		      T4 = sqrt(1.0 + 2.0 * (Vgs - T0) / T1);
865 		      T2 = T1 * (T4 - 1.0);
866 		      T3 = 0.5 * T2 * T2 / T1; /* T3 = Vpoly */
867 		      T7 = 1.12 - T3 - 0.05;
868 		      T6 = sqrt(T7 * T7 + 0.224);
869 		      T5 = 1.12 - 0.5 * (T7 + T6);
870 		      Vgs_eff = Vgs - T5;
871 		      dVgs_eff_dVg = 1.0 - (0.5 - 0.5 / T4) * (1.0 + T7 / T6);
872 		  }
873 		  else
874 		  {   Vgs_eff = Vgs;
875 		      dVgs_eff_dVg = 1.0;
876 		  }
877 
878 
879 		  Leff = pParam->B3SOIFDleff;
880 
881 		  if (selfheat) {
882 		      Vtm = KboQ * Temp;
883                       dVtm_dT = KboQ;
884 		  }
885 		  else {
886 		      Vtm = model->B3SOIFDvtm;
887                       dVtm_dT = 0.0;
888 		  }
889 
890 		  V0 = vbi - phi;
891 
892 /* Prepare Vbs0t */
893 		      T0 = -pParam->B3SOIFDdvbd1 * pParam->B3SOIFDleff / pParam->B3SOIFDlitl;
894 		      T1 = pParam->B3SOIFDdvbd0 * (exp(0.5*T0) + 2*exp(T0));
895 		      T2 = T1 * (vbi - phi);
896 		      T3 = 0.5 * model->B3SOIFDqsi / model->B3SOIFDcsi;
897 		      Vbs0t = phi - T3 + pParam->B3SOIFDvbsa + T2;
898                       if (selfheat)
899                          dVbs0t_dT = T1 * dvbi_dT;
900                       else
901                          dVbs0t_dT = 0.0;
902 
903 /* Prepare Vbs0 */
904 			  T0 = 1 + model->B3SOIFDcsieff / Cbox;
905                           T1 = pParam->B3SOIFDkb1 / T0;
906 			  T2 = T1 * (Vbs0t - Vesfb);
907 
908                           /* T6 is Vbs0 before limiting */
909                           T6 = Vbs0t - T2;
910                           dT6_dVe = T1;
911                           if (selfheat)
912                              dT6_dT = dVbs0t_dT - T1 * (dVbs0t_dT + dvfbb_dT);
913                           else
914                              dT6_dT = 0.0;
915 
916                           /* limit Vbs0 to below phi */
917                           T1 = phi - pParam->B3SOIFDdelp;
918                           T2 = T1 - T6 - DELT_Vbseff;
919                           T3 = sqrt(T2 * T2 + 4.0 * DELT_Vbseff);
920                           Vbs0 = T1 - 0.5 * (T2 + T3);
921                           T4 = 0.5 * (1 + T2/T3);
922                           dVbs0_dVe = T4 * dT6_dVe;
923                           if (selfheat)  dVbs0_dT = T4 * dT6_dT;
924                           else  dVbs0_dT = 0.0;
925 
926 			  T1 = Vbs0t - Vbs0 - DELT_Vbsmos;
927 			  T2 = sqrt(T1 * T1 + DELT_Vbsmos * DELT_Vbsmos);
928 			  T3 = 0.5 * (T1 + T2);
929 			  T4 = T3 * model->B3SOIFDcsieff / model->B3SOIFDqsieff;
930 			  Vbs0mos = Vbs0 - 0.5 * T3 * T4;
931                           T5 = 0.5 * T4 * (1 + T1 / T2);
932 			  dVbs0mos_dVe = dVbs0_dVe * (1 + T5);
933                           if (selfheat)
934 			     dVbs0mos_dT = dVbs0_dT - (dVbs0t_dT - dVbs0_dT) * T5;
935                           else
936 			     dVbs0mos_dT = 0.0;
937 
938 /* Prepare Vthfd - treat Vbs0mos as if it were independent variable Vb */
939 		     Phis = phi - Vbs0mos;
940 		     dPhis_dVb = -1;
941 		     sqrtPhis = sqrt(Phis);
942 		     dsqrtPhis_dVb = -0.5 / sqrtPhis;
943 		     Xdep = Xdep0 * sqrtPhis / sqrtPhi;
944 		     dXdep_dVb = (Xdep0 / sqrtPhi)
945 				 * dsqrtPhis_dVb;
946 		     sqrtXdep = sqrt(Xdep);
947 
948 		     T0 = pParam->B3SOIFDdvt2 * Vbs0mos;
949 		     if (T0 >= - 0.5)
950 		     {   T1 = 1.0 + T0;
951 			 T2 = pParam->B3SOIFDdvt2;
952 		     }
953 		     else /* Added to avoid any discontinuity problems caused by dvt2*/
954 		     {   T4 = 1.0 / (3.0 + 8.0 * T0);
955 			 T1 = (1.0 + 3.0 * T0) * T4;
956 			 T2 = pParam->B3SOIFDdvt2 * T4 * T4;
957 		     }
958 		     lt1 = model->B3SOIFDfactor1 * sqrtXdep * T1;
959 		     dlt1_dVb = model->B3SOIFDfactor1 * (0.5 / sqrtXdep * T1 * dXdep_dVb
960 			      + sqrtXdep * T2);
961 
962 		     T0 = pParam->B3SOIFDdvt2w * Vbs0mos;
963 		     if (T0 >= - 0.5)
964 		     {   T1 = 1.0 + T0;
965 			 T2 = pParam->B3SOIFDdvt2w;
966 		     }
967 		     else /* Added to avoid any discontinuity problems caused by
968 			     dvt2w */
969 		     {   T4 = 1.0 / (3.0 + 8.0 * T0);
970 			 T1 = (1.0 + 3.0 * T0) * T4;
971 			 T2 = pParam->B3SOIFDdvt2w * T4 * T4;
972 		     }
973 		     ltw= model->B3SOIFDfactor1 * sqrtXdep * T1;
974 		     dltw_dVb = model->B3SOIFDfactor1 * (0.5 / sqrtXdep * T1 * dXdep_dVb
975 			      + sqrtXdep * T2);
976 
977 		     T0 = -0.5 * pParam->B3SOIFDdvt1 * Leff / lt1;
978 		     if (T0 > -EXP_THRESHOLD)
979 		     {   T1 = exp(T0);
980 			 dT1_dVb = -T0 / lt1 * T1 * dlt1_dVb;
981 			 Theta0 = T1 * (1.0 + 2.0 * T1);
982 			 dTheta0_dVb = (1.0 + 4.0 * T1) * dT1_dVb;
983 		     }
984 		     else
985 		     {   T1 = MIN_EXP;
986 			 Theta0 = T1 * (1.0 + 2.0 * T1);
987 			 dTheta0_dVb = 0.0;
988 		     }
989 		     here->B3SOIFDthetavth = pParam->B3SOIFDdvt0 * Theta0;
990 		     Delt_vth = here->B3SOIFDthetavth * V0;
991 		     dDelt_vth_dVb = pParam->B3SOIFDdvt0 * dTheta0_dVb * V0;
992                      if (selfheat) dDelt_vth_dT = here->B3SOIFDthetavth * dvbi_dT;
993                      else dDelt_vth_dT = 0.0;
994 
995 		     T0 = -0.5*pParam->B3SOIFDdvt1w * pParam->B3SOIFDweff*Leff/ltw;
996 		     if (T0 > -EXP_THRESHOLD)
997 		     {   T1 = exp(T0);
998 			 T2 = T1 * (1.0 + 2.0 * T1);
999 			 dT1_dVb = -T0 / ltw * T1 * dltw_dVb;
1000 			 dT2_dVb = (1.0 + 4.0 * T1) * dT1_dVb;
1001 		     }
1002 		     else
1003 		     {   T1 = MIN_EXP;
1004 			 T2 = T1 * (1.0 + 2.0 * T1);
1005 			 dT2_dVb = 0.0;
1006 		     }
1007 		     T0 = pParam->B3SOIFDdvt0w * T2;
1008 		     DeltVthw = T0 * V0;
1009 		     dDeltVthw_dVb = pParam->B3SOIFDdvt0w * dT2_dVb * V0;
1010                      if (selfheat)  dDeltVthw_dT = T0 * dvbi_dT;
1011                      else  dDeltVthw_dT = 0.0;
1012 
1013 		     T0 = sqrt(1.0 + pParam->B3SOIFDnlx / Leff);
1014                      T1 = (pParam->B3SOIFDkt1 + pParam->B3SOIFDkt1l / Leff
1015                            + pParam->B3SOIFDkt2 * Vbs0mos);
1016 		     DeltVthtemp = pParam->B3SOIFDk1 * (T0 - 1.0) * sqrtPhi + T1 * TempRatio;
1017                      if (selfheat)   dDeltVthtemp_dT = T1 / model->B3SOIFDtnom;
1018                      else   dDeltVthtemp_dT = 0.0;
1019 
1020 		     tmp2 = model->B3SOIFDtox * phi
1021 			  / (pParam->B3SOIFDweff + pParam->B3SOIFDw0);
1022 
1023 		     T3 = pParam->B3SOIFDeta0 + pParam->B3SOIFDetab * Vbs0mos;
1024 		     if (T3 < 1.0e-4) /* avoid discontinuity problems caused by etab */
1025 		     {   T9 = 1.0 / (3.0 - 2.0e4 * T3);
1026 			 T3 = (2.0e-4 - T3) * T9;
1027 			 T4 = T9 * T9 * pParam->B3SOIFDetab;
1028 			 dT3_dVb = T4;
1029 		     }
1030 		     else
1031 		     {
1032 			 dT3_dVb = pParam->B3SOIFDetab;
1033 		     }
1034 		     DIBL_Sft = T3 * pParam->B3SOIFDtheta0vb0 * Vds;
1035 		     dDIBL_Sft_dVd = T3 * pParam->B3SOIFDtheta0vb0;
1036 		     dDIBL_Sft_dVb = pParam->B3SOIFDtheta0vb0 * Vds * dT3_dVb;
1037 
1038 		     Vthfd = model->B3SOIFDtype * pParam->B3SOIFDvth0 + pParam->B3SOIFDk1
1039 			 * (sqrtPhis - sqrtPhi) - pParam->B3SOIFDk2
1040 			 * Vbs0mos-Delt_vth-DeltVthw +(pParam->B3SOIFDk3 +pParam->B3SOIFDk3b
1041 			 * Vbs0mos) * tmp2 + DeltVthtemp - DIBL_Sft;
1042 
1043 		     T6 = pParam->B3SOIFDk3b * tmp2 - pParam->B3SOIFDk2 +
1044 			  pParam->B3SOIFDkt2 * TempRatio;
1045 		     dVthfd_dVd = -dDIBL_Sft_dVd;
1046 		     T7 = pParam->B3SOIFDk1 * dsqrtPhis_dVb
1047 			  - dDelt_vth_dVb - dDeltVthw_dVb
1048 			  + T6 - dDIBL_Sft_dVb;
1049 		     dVthfd_dVe = T7 * dVbs0mos_dVe;
1050                      if (selfheat)
1051                         dVthfd_dT = dDeltVthtemp_dT - dDelt_vth_dT - dDeltVthw_dT
1052                                   + T7 * dVbs0mos_dT;
1053                      else
1054                         dVthfd_dT = 0.0;
1055 
1056 /* Effective Vbs0 and Vbs0t for all Vgs */
1057 		     T1 = Vthfd - Vgs_eff - DELT_Vbs0eff;
1058 		     T2 = sqrt(T1 * T1 + DELT_Vbs0eff * DELT_Vbs0eff );
1059 
1060 		     Vbs0teff = Vbs0t - 0.5 * (T1 + T2);
1061 		     dVbs0teff_dVg = 0.5  * (1 + T1/T2) * dVgs_eff_dVg;
1062 		     dVbs0teff_dVd = - 0.5 * (1 + T1 / T2) * dVthfd_dVd;
1063 		     dVbs0teff_dVe = - 0.5 * (1 + T1 / T2) * dVthfd_dVe;
1064                      if (selfheat)
1065                         dVbs0teff_dT = dVbs0t_dT - 0.5 * (1 + T1 / T2) * dVthfd_dT;
1066                      else
1067                         dVbs0teff_dT = 0.0;
1068 
1069                      /* Calculate nfb */
1070                      T3 = 1 / (K1 * K1);
1071 			T4 = pParam->B3SOIFDkb3 * Cbox / model->B3SOIFDcox;
1072 			T8 = sqrt(phi - Vbs0mos);
1073 			T5 = sqrt(1 + 4 * T3 * (phi + K1 * T8 - Vbs0mos));
1074 			T6 = 1 + T4 * T5;
1075 			Nfb = model->B3SOIFDnfb = 1 / T6;
1076 			T7 = 2 * T3 * T4 * Nfb * Nfb / T5 * (0.5 * K1 / T8 + 1);
1077 			Vbs0eff = Vbs0 - Nfb * 0.5 * (T1 + T2);
1078 			dVbs0eff_dVg = Nfb * 0.5  * (1 + T1/T2) * dVgs_eff_dVg;
1079 			dVbs0eff_dVd = - Nfb * 0.5 * (1 + T1 / T2) * dVthfd_dVd;
1080 			dVbs0eff_dVe = dVbs0_dVe - Nfb * 0.5 * (1 + T1 / T2)
1081 				     * dVthfd_dVe - T7 * 0.5 * (T1 + T2) * dVbs0mos_dVe;
1082                         if (selfheat)
1083                            dVbs0eff_dT = dVbs0_dT - Nfb * 0.5 * (1 + T1 / T2)
1084                                        * dVthfd_dT - T7 * 0.5 * (T1 + T2) * dVbs0mos_dT;
1085                         else
1086                            dVbs0eff_dT = 0.0;
1087 
1088 /* Simple check of Vbs */
1089 /* Prepare Vbsdio */
1090                         Vbs = Vbsdio = Vbs0eff;
1091                         dVbsdio_dVg = dVbs0eff_dVg;
1092                         dVbsdio_dVd = dVbs0eff_dVd;
1093                         dVbsdio_dVe = dVbs0eff_dVe;
1094                         dVbsdio_dT = dVbs0eff_dT;
1095                         dVbsdio_dVb = 0.0;
1096 
1097 /* Prepare Vbseff */
1098 			T1 = Vbs0teff - Vbsdio - DELT_Vbsmos;
1099 			T2 = sqrt(T1 * T1 + DELT_Vbsmos * DELT_Vbsmos);
1100 			T3 = 0.5 * (T1 + T2);
1101 			T5 = 0.5 * (1 + T1/T2);
1102 			dT3_dVg = T5 * (dVbs0teff_dVg - dVbsdio_dVg);
1103 			dT3_dVd = T5 * (dVbs0teff_dVd - dVbsdio_dVd);
1104 			dT3_dVb = - T5 * dVbsdio_dVb;
1105 			dT3_dVe = T5 * (dVbs0teff_dVe - dVbsdio_dVe);
1106                         if (selfheat)  dT3_dT = T5 * (dVbs0teff_dT - dVbsdio_dT);
1107                         else  dT3_dT = 0.0;
1108 			T4 = T3 * model->B3SOIFDcsieff / model->B3SOIFDqsieff;
1109 
1110 			Vbsmos = Vbsdio - 0.5 * T3 * T4;
1111 			dVbsmos_dVg = dVbsdio_dVg - T4 * dT3_dVg;
1112 			dVbsmos_dVd = dVbsdio_dVd - T4 * dT3_dVd;
1113 			dVbsmos_dVb = dVbsdio_dVb - T4 * dT3_dVb;
1114 			dVbsmos_dVe = dVbsdio_dVe - T4 * dT3_dVe;
1115                         if (selfheat)  dVbsmos_dT = dVbsdio_dT - T4 * dT3_dT;
1116                         else  dVbsmos_dT = 0.0;
1117 
1118 /* Prepare Vcs */
1119 		     Vcs = Vbsdio - Vbs0eff;
1120 		     dVcs_dVb = dVbsdio_dVb;
1121 		     dVcs_dVg = dVbsdio_dVg - dVbs0eff_dVg;
1122 		     dVcs_dVd = dVbsdio_dVd - dVbs0eff_dVd;
1123 		     dVcs_dVe = dVbsdio_dVe - dVbs0eff_dVe;
1124                      dVcs_dT = dVbsdio_dT - dVbs0eff_dT;
1125 
1126 /* Check Vps */
1127                      /* Note : if Vps is less Vbs0eff => non-physical */
1128                      T1 = Vps - Vbs0eff + DELT_Vbs0dio;
1129                      T2 = sqrt(T1 * T1 + DELT_Vbs0dio * DELT_Vbs0dio);
1130                      T3 = 0.5 * (1 + T1/T2);
1131                      Vpsdio = Vbs0eff + 0.5 * (T1 + T2);
1132                      dVpsdio_dVg = (1 - T3) * dVbs0eff_dVg;
1133                      dVpsdio_dVd = (1 - T3) * dVbs0eff_dVd;
1134                      dVpsdio_dVe = (1 - T3) * dVbs0eff_dVe;
1135                      if (selfheat)  dVpsdio_dT = (1 - T3) * dVbs0eff_dT;
1136                      else  dVpsdio_dT = 0.0;
1137                      dVpsdio_dVp = T3;
1138                      Vbp = Vbsdio - Vpsdio;
1139                      dVbp_dVb = dVbsdio_dVb;
1140                      dVbp_dVg = dVbsdio_dVg - dVpsdio_dVg;
1141                      dVbp_dVd = dVbsdio_dVd - dVpsdio_dVd;
1142                      dVbp_dVe = dVbsdio_dVe - dVpsdio_dVe;
1143                      dVbp_dT = dVbsdio_dT - dVpsdio_dT;
1144                      dVbp_dVp = - dVpsdio_dVp;
1145 
1146                   here->B3SOIFDvbsdio = Vbsdio;
1147                   here->B3SOIFDvbs0eff = Vbs0eff;
1148 
1149 		  T1 = phi - pParam->B3SOIFDdelp;
1150 		  T2 = T1 - Vbsmos - DELT_Vbseff;
1151 		  T3 = sqrt(T2 * T2 + 4.0 * DELT_Vbseff * T1);
1152 		  Vbseff = T1 - 0.5 * (T2 + T3);
1153 		  T4 = 0.5 * (1 + T2/T3);
1154 		  dVbseff_dVg = T4 * dVbsmos_dVg;
1155 		  dVbseff_dVd = T4 * dVbsmos_dVd;
1156 		  dVbseff_dVb = T4 * dVbsmos_dVb;
1157 		  dVbseff_dVe = T4 * dVbsmos_dVe;
1158                   if (selfheat)  dVbseff_dT = T4 * dVbsmos_dT;
1159                   else  dVbseff_dT = 0.0;
1160 
1161                   here->B3SOIFDvbseff = Vbseff;
1162 
1163 		  Phis = phi - Vbseff;
1164 		  dPhis_dVb = -1;
1165 		  sqrtPhis = sqrt(Phis);
1166 		  dsqrtPhis_dVb = -0.5 / sqrtPhis ;
1167 
1168 		  Xdep = Xdep0 * sqrtPhis / sqrtPhi;
1169 		  dXdep_dVb = (Xdep0 / sqrtPhi)
1170 			    * dsqrtPhis_dVb;
1171 
1172 /* Vth Calculation */
1173 		  T3 = sqrt(Xdep);
1174 
1175 		  T0 = pParam->B3SOIFDdvt2 * Vbseff;
1176 		  if (T0 >= - 0.5)
1177 		  {   T1 = 1.0 + T0;
1178 		      T2 = pParam->B3SOIFDdvt2 ;
1179 		  }
1180 		  else /* Added to avoid any discontinuity problems caused by dvt2 */
1181 		  {   T4 = 1.0 / (3.0 + 8.0 * T0);
1182 		      T1 = (1.0 + 3.0 * T0) * T4;
1183 		      T2 = pParam->B3SOIFDdvt2 * T4 * T4 ;
1184 		  }
1185 		  lt1 = model->B3SOIFDfactor1 * T3 * T1;
1186 		  dlt1_dVb =model->B3SOIFDfactor1 * (0.5 / T3 * T1 * dXdep_dVb + T3 * T2);
1187 
1188 		  T0 = pParam->B3SOIFDdvt2w * Vbseff;
1189 		  if (T0 >= - 0.5)
1190 		  {   T1 = 1.0 + T0;
1191 		      T2 = pParam->B3SOIFDdvt2w ;
1192 		  }
1193 		  else /* Added to avoid any discontinuity problems caused by dvt2w */
1194 		  {   T4 = 1.0 / (3.0 + 8.0 * T0);
1195 		      T1 = (1.0 + 3.0 * T0) * T4;
1196 		      T2 = pParam->B3SOIFDdvt2w * T4 * T4 ;
1197 		  }
1198 		  ltw= model->B3SOIFDfactor1 * T3 * T1;
1199 		  dltw_dVb=model->B3SOIFDfactor1*(0.5 / T3 * T1 * dXdep_dVb + T3 * T2);
1200 
1201 		  T0 = -0.5 * pParam->B3SOIFDdvt1 * Leff / lt1;
1202 		  if (T0 > -EXP_THRESHOLD)
1203 		  {   T1 = exp(T0);
1204 		      Theta0 = T1 * (1.0 + 2.0 * T1);
1205 		      dT1_dVb = -T0 / lt1 * T1 * dlt1_dVb;
1206 		      dTheta0_dVb = (1.0 + 4.0 * T1) * dT1_dVb;
1207 		  }
1208 		  else
1209 		  {   T1 = MIN_EXP;
1210 		      Theta0 = T1 * (1.0 + 2.0 * T1);
1211 		      dTheta0_dVb = 0.0;
1212 		  }
1213 
1214 		  here->B3SOIFDthetavth = pParam->B3SOIFDdvt0 * Theta0;
1215 		  Delt_vth = here->B3SOIFDthetavth * V0;
1216 		  dDelt_vth_dVb = pParam->B3SOIFDdvt0 * dTheta0_dVb * V0;
1217                   if (selfheat)  dDelt_vth_dT = here->B3SOIFDthetavth * dvbi_dT;
1218                   else  dDelt_vth_dT = 0.0;
1219 
1220 		  T0 = -0.5 * pParam->B3SOIFDdvt1w * pParam->B3SOIFDweff * Leff / ltw;
1221 		  if (T0 > -EXP_THRESHOLD)
1222 		  {   T1 = exp(T0);
1223 		      T2 = T1 * (1.0 + 2.0 * T1);
1224 		      dT1_dVb = -T0 / ltw * T1 * dltw_dVb;
1225 		      dT2_dVb = (1.0 + 4.0 * T1) * dT1_dVb;
1226 		  }
1227 		  else
1228 		  {   T1 = MIN_EXP;
1229 		      T2 = T1 * (1.0 + 2.0 * T1);
1230 		      dT2_dVb = 0.0;
1231 		  }
1232 
1233 		  T0 = pParam->B3SOIFDdvt0w * T2;
1234 		  DeltVthw = T0 * V0;
1235 		  dDeltVthw_dVb = pParam->B3SOIFDdvt0w * dT2_dVb * V0;
1236                   if (selfheat)   dDeltVthw_dT = T0 * dvbi_dT;
1237                   else   dDeltVthw_dT = 0.0;
1238 
1239 		  T0 = sqrt(1.0 + pParam->B3SOIFDnlx / Leff);
1240                   T1 = (pParam->B3SOIFDkt1 + pParam->B3SOIFDkt1l / Leff
1241                         + pParam->B3SOIFDkt2 * Vbseff);
1242 		  DeltVthtemp = pParam->B3SOIFDk1 * (T0 - 1.0) * sqrtPhi + T1 * TempRatio;
1243                   if (selfheat)
1244                      dDeltVthtemp_dT = T1 / model->B3SOIFDtnom;
1245                   else
1246                      dDeltVthtemp_dT = 0.0;
1247 
1248 		  tmp2 = model->B3SOIFDtox * phi
1249 		       / (pParam->B3SOIFDweff + pParam->B3SOIFDw0);
1250 
1251 		  T3 = pParam->B3SOIFDeta0 + pParam->B3SOIFDetab * Vbseff;
1252 		  if (T3 < 1.0e-4) /* avoid  discontinuity problems caused by etab */
1253 		  {   T9 = 1.0 / (3.0 - 2.0e4 * T3);
1254 		      T3 = (2.0e-4 - T3) * T9;
1255 		      T4 = T9 * T9 * pParam->B3SOIFDetab;
1256 		      dT3_dVb = T4 ;
1257 		  }
1258 		  else
1259 		  {
1260 		      dT3_dVb = pParam->B3SOIFDetab ;
1261 		  }
1262 		  DIBL_Sft = T3 * pParam->B3SOIFDtheta0vb0 * Vds;
1263 		  dDIBL_Sft_dVd = pParam->B3SOIFDtheta0vb0 * T3;
1264 		  dDIBL_Sft_dVb = pParam->B3SOIFDtheta0vb0 * Vds * dT3_dVb;
1265 
1266 		  Vth = model->B3SOIFDtype * pParam->B3SOIFDvth0 + pParam->B3SOIFDk1
1267 		      * (sqrtPhis - sqrtPhi) - pParam->B3SOIFDk2
1268 		      * Vbseff- Delt_vth - DeltVthw +(pParam->B3SOIFDk3 + pParam->B3SOIFDk3b
1269 		      * Vbseff) * tmp2 + DeltVthtemp - DIBL_Sft;
1270 
1271 		  here->B3SOIFDvon = Vth;
1272 
1273 		  T6 = pParam->B3SOIFDk3b * tmp2 - pParam->B3SOIFDk2
1274 		       + pParam->B3SOIFDkt2 * TempRatio;
1275 		  dVth_dVb = pParam->B3SOIFDk1 * dsqrtPhis_dVb
1276 			   - dDelt_vth_dVb - dDeltVthw_dVb
1277 			   + T6 - dDIBL_Sft_dVb;  /*  this is actually dVth_dVbseff  */
1278 		  dVth_dVd = -dDIBL_Sft_dVd;
1279                   if (selfheat)  dVth_dT = dDeltVthtemp_dT - dDelt_vth_dT - dDeltVthw_dT;
1280                   else  dVth_dT = 0.0;
1281 
1282 /* Calculate n */
1283 		  T2 = pParam->B3SOIFDnfactor * EPSSI / Xdep;
1284 		  dT2_dVb = - T2 / Xdep * dXdep_dVb;
1285 
1286 		  T3 = pParam->B3SOIFDcdsc + pParam->B3SOIFDcdscb * Vbseff
1287 		       + pParam->B3SOIFDcdscd * Vds;
1288 		  dT3_dVb = pParam->B3SOIFDcdscb;
1289 		  dT3_dVd = pParam->B3SOIFDcdscd;
1290 
1291 		  T4 = (T2 + T3 * Theta0 + pParam->B3SOIFDcit) / model->B3SOIFDcox;
1292 		  dT4_dVb = (dT2_dVb + Theta0 * dT3_dVb + dTheta0_dVb * T3)
1293                             / model->B3SOIFDcox;
1294 		  dT4_dVd = Theta0 * dT3_dVd / model->B3SOIFDcox;
1295 
1296 		  if (T4 >= -0.5)
1297 		  {   n = 1.0 + T4;
1298 		      dn_dVb = dT4_dVb;
1299 		      dn_dVd = dT4_dVd;
1300 		  }
1301 		  else
1302 		   /* avoid  discontinuity problems caused by T4 */
1303 		  {   T0 = 1.0 / (3.0 + 8.0 * T4);
1304 		      n = (1.0 + 3.0 * T4) * T0;
1305 		      T0 *= T0;
1306 		      dn_dVb = T0 * dT4_dVb;
1307 		      dn_dVd = T0 * dT4_dVd;
1308 		  }
1309 
1310 /* Effective Vgst (Vgsteff) Calculation */
1311 
1312 		  Vgst = Vgs_eff - Vth;
1313 
1314 		  T10 = 2.0 * n * Vtm;
1315 		  VgstNVt = Vgst / T10;
1316 		  ExpArg = (2.0 * pParam->B3SOIFDvoff - Vgst) / T10;
1317 
1318 		  /* MCJ: Very small Vgst */
1319 		  if (VgstNVt > EXP_THRESHOLD)
1320 		  {   Vgsteff = Vgst;
1321                       /* T0 is dVgsteff_dVbseff */
1322                       T0 = -dVth_dVb;
1323 		      dVgsteff_dVg = dVgs_eff_dVg + T0 * dVbseff_dVg;
1324 		      dVgsteff_dVd = -dVth_dVd + T0 * dVbseff_dVd;
1325 		      dVgsteff_dVb = T0 * dVbseff_dVb;
1326                       dVgsteff_dVe = T0 * dVbseff_dVe;
1327                       if (selfheat)
1328                          dVgsteff_dT  = -dVth_dT + T0 * dVbseff_dT;
1329                       else
1330                          dVgsteff_dT = 0.0;
1331 		  }
1332 		  else if (ExpArg > EXP_THRESHOLD)
1333 		  {   T0 = (Vgst - pParam->B3SOIFDvoff) / (n * Vtm);
1334 		      ExpVgst = exp(T0);
1335 		      Vgsteff = Vtm * pParam->B3SOIFDcdep0 / model->B3SOIFDcox * ExpVgst;
1336 		      T3 = Vgsteff / (n * Vtm) ;
1337                       /* T1 is dVgsteff_dVbseff */
1338 		      T1  = -T3 * (dVth_dVb + T0 * Vtm * dn_dVb);
1339 		      dVgsteff_dVg = T3 * dVgs_eff_dVg + T1 * dVbseff_dVg;
1340 		      dVgsteff_dVd = -T3 * (dVth_dVd + T0 * Vtm * dn_dVd) + T1 * dVbseff_dVd;
1341                       dVgsteff_dVe = T1 * dVbseff_dVe;
1342                       dVgsteff_dVb = T1 * dVbseff_dVb;
1343                       if (selfheat)
1344                          dVgsteff_dT = -T3 * (dVth_dT + T0 * dVtm_dT * n)
1345                                      + Vgsteff / Temp + T1 * dVbseff_dT;
1346                       else
1347                          dVgsteff_dT = 0.0;
1348 		  }
1349 		  else
1350 		  {   ExpVgst = exp(VgstNVt);
1351 		      T1 = T10 * log(1.0 + ExpVgst);
1352 		      dT1_dVg = ExpVgst / (1.0 + ExpVgst);
1353 		      dT1_dVb = -dT1_dVg * (dVth_dVb + Vgst / n * dn_dVb)
1354 			      + T1 / n * dn_dVb;
1355 		      dT1_dVd = -dT1_dVg * (dVth_dVd + Vgst / n * dn_dVd)
1356 			      + T1 / n * dn_dVd;
1357                       T3 = (1.0 / Temp);
1358                       if (selfheat)
1359                          dT1_dT = -dT1_dVg * (dVth_dT + Vgst * T3) + T1 * T3;
1360                       else
1361                          dT1_dT = 0.0;
1362 
1363 		      dT2_dVg = -model->B3SOIFDcox / (Vtm * pParam->B3SOIFDcdep0)
1364 			      * exp(ExpArg);
1365 		      T2 = 1.0 - T10 * dT2_dVg;
1366 		      dT2_dVd = -dT2_dVg * (dVth_dVd - 2.0 * Vtm * ExpArg * dn_dVd)
1367 			      + (T2 - 1.0) / n * dn_dVd;
1368 		      dT2_dVb = -dT2_dVg * (dVth_dVb - 2.0 * Vtm * ExpArg * dn_dVb)
1369 			      + (T2 - 1.0) / n * dn_dVb;
1370                       if (selfheat)
1371                          dT2_dT = -dT2_dVg * (dVth_dT - ExpArg * T10 * T3);
1372                       else
1373                          dT2_dT = 0.0;
1374 
1375 		      Vgsteff = T1 / T2;
1376 		      T3 = T2 * T2;
1377                       /*  T4 is dVgsteff_dVbseff  */
1378 		      T4 = (T2 * dT1_dVb - T1 * dT2_dVb) / T3;
1379                       dVgsteff_dVb = T4 * dVbseff_dVb;
1380                       dVgsteff_dVe = T4 * dVbseff_dVe;
1381 		      dVgsteff_dVg = (T2 * dT1_dVg - T1 * dT2_dVg) / T3 * dVgs_eff_dVg
1382                                      + T4 * dVbseff_dVg;
1383 		      dVgsteff_dVd = (T2 * dT1_dVd - T1 * dT2_dVd) / T3 + T4 * dVbseff_dVd;
1384                       if (selfheat)
1385                          dVgsteff_dT = (T2 * dT1_dT - T1 * dT2_dT) / T3 + T4 * dVbseff_dT;
1386                       else
1387                          dVgsteff_dT = 0.0;
1388 		  }
1389 		  Vgst2Vtm = Vgsteff + 2.0 * Vtm;
1390                   if (selfheat)  dVgst2Vtm_dT = 2.0 * dVtm_dT;
1391                   else  dVgst2Vtm_dT = 0.0;
1392 
1393 /* Calculate Effective Channel Geometry */
1394 		  T9 = sqrtPhis - sqrtPhi;
1395 		  Weff = pParam->B3SOIFDweff - 2.0 * (pParam->B3SOIFDdwg * Vgsteff
1396 		       + pParam->B3SOIFDdwb * T9);
1397 		  dWeff_dVg = -2.0 * pParam->B3SOIFDdwg;
1398 		  dWeff_dVb = -2.0 * pParam->B3SOIFDdwb * dsqrtPhis_dVb;
1399 
1400 		  if (Weff < 2.0e-8) /* to avoid the discontinuity problem due to Weff*/
1401 		  {   T0 = 1.0 / (6.0e-8 - 2.0 * Weff);
1402 		      Weff = 2.0e-8 * (4.0e-8 - Weff) * T0;
1403 		      T0 *= T0 * 4.0e-16;
1404 		      dWeff_dVg *= T0;
1405 		      dWeff_dVb *= T0;
1406 		  }
1407 
1408 		  T0 = pParam->B3SOIFDprwg * Vgsteff + pParam->B3SOIFDprwb * T9;
1409 		  if (T0 >= -0.9)
1410 		  {   Rds = rds0 * (1.0 + T0);
1411 		      dRds_dVg = rds0 * pParam->B3SOIFDprwg;
1412 		      dRds_dVb = rds0 * pParam->B3SOIFDprwb * dsqrtPhis_dVb;
1413                       if (selfheat)  dRds_dT = (1.0 + T0) * drds0_dT;
1414                       else  dRds_dT = 0.0;
1415 		  }
1416 		  else
1417 		   /* to avoid the discontinuity problem due to prwg and prwb*/
1418 		  {   T1 = 1.0 / (17.0 + 20.0 * T0);
1419 		      Rds = rds0 * (0.8 + T0) * T1;
1420 		      T1 *= T1;
1421 		      dRds_dVg = rds0 * pParam->B3SOIFDprwg * T1;
1422 		      dRds_dVb = rds0 * pParam->B3SOIFDprwb * dsqrtPhis_dVb
1423 			       * T1;
1424                       if (selfheat)  dRds_dT = (0.8 + T0) * T1 * drds0_dT;
1425                       else  dRds_dT = 0.0;
1426 		  }
1427 
1428 /* Calculate Abulk */
1429                   if (pParam->B3SOIFDa0 == 0.0)
1430                   {
1431                      Abulk0 = Abulk = dAbulk0_dVb = dAbulk_dVg = dAbulk_dVb = 0.0;
1432                   }
1433                   else
1434                   {
1435 		     T1 = 0.5 * pParam->B3SOIFDk1 / sqrtPhi;
1436 		     T9 = sqrt(model->B3SOIFDxj * Xdep);
1437 		     tmp1 = Leff + 2.0 * T9;
1438 		     T5 = Leff / tmp1;
1439 		     tmp2 = pParam->B3SOIFDa0 * T5;
1440 		     tmp3 = pParam->B3SOIFDweff + pParam->B3SOIFDb1;
1441 		     tmp4 = pParam->B3SOIFDb0 / tmp3;
1442 		     T2 = tmp2 + tmp4;
1443 		     dT2_dVb = -T9 * tmp2 / tmp1 / Xdep * dXdep_dVb;
1444 		     T6 = T5 * T5;
1445 		     T7 = T5 * T6;
1446 
1447 		     Abulk0 = T1 * T2;
1448 		     dAbulk0_dVb = T1 * dT2_dVb;
1449 
1450 		     T8 = pParam->B3SOIFDags * pParam->B3SOIFDa0 * T7;
1451 		     dAbulk_dVg = -T1 * T8;
1452 		     Abulk = Abulk0 + dAbulk_dVg * Vgsteff;
1453 
1454 		     dAbulk_dVb = dAbulk0_dVb - T8 * Vgsteff * 3.0 * T1 * dT2_dVb
1455                                              / tmp2;
1456                   }
1457 
1458                   if (Abulk0 < 0.01)
1459                   {
1460                      T9 = 1.0 / (3.0 - 200.0 * Abulk0);
1461                      Abulk0 = (0.02 - Abulk0) * T9;
1462                      dAbulk0_dVb *= T9 * T9;
1463                   }
1464 
1465                   if (Abulk < 0.01)
1466                   {
1467                      T9 = 1.0 / (3.0 - 200.0 * Abulk);
1468                      Abulk = (0.02 - Abulk) * T9;
1469                      dAbulk_dVb *= T9 * T9;
1470                   }
1471 
1472 		  T2 = pParam->B3SOIFDketa * Vbseff;
1473 		  if (T2 >= -0.9)
1474 		  {   T0 = 1.0 / (1.0 + T2);
1475 		      dT0_dVb = -pParam->B3SOIFDketa * T0 * T0 ;
1476 		  }
1477 		  else
1478 		  /* added to avoid the problems caused by Keta */
1479 		  {   T1 = 1.0 / (0.8 + T2);
1480 		      T0 = (17.0 + 20.0 * T2) * T1;
1481 		      dT0_dVb = -pParam->B3SOIFDketa * T1 * T1 ;
1482 		  }
1483 		  dAbulk_dVg *= T0;
1484 		  dAbulk_dVb = dAbulk_dVb * T0 + Abulk * dT0_dVb;
1485 		  dAbulk0_dVb = dAbulk0_dVb * T0 + Abulk0 * dT0_dVb;
1486 		  Abulk *= T0;
1487 		  Abulk0 *= T0;
1488 
1489 		  Abulk += 1;
1490 		  Abulk0 += 1;
1491 
1492 /* Prepare Abeff */
1493 		      T0 = pParam->B3SOIFDabp * Vgst2Vtm;
1494 		      T1 = 1 - Vcs / T0 - DELT_Xcsat;
1495 		      T2 = sqrt(T1 * T1 + DELT_Xcsat * DELT_Xcsat);
1496 		      T3 = 1 - 0.5 * (T1 + T2);
1497 		      T5 = -0.5 * (1 + T1 / T2);
1498 		      dT1_dVg = Vcs / Vgst2Vtm / T0;
1499 		      dT3_dVg = T5 * dT1_dVg;
1500 		      dT1_dVc = - 1 / T0;
1501 		      dT3_dVc = T5 * dT1_dVc;
1502 
1503 		      Xcsat = pParam->B3SOIFDmxc * T3 * T3 + (1 - pParam->B3SOIFDmxc)*T3;
1504 		      T4 = 2 * pParam->B3SOIFDmxc * T3 + (1 - pParam->B3SOIFDmxc);
1505 		      dXcsat_dVg = T4 * dT3_dVg;
1506 		      dXcsat_dVc = T4 * dT3_dVc;
1507 
1508 		      Abeff = Xcsat * Abulk + (1 - Xcsat) * model->B3SOIFDadice;
1509 		      T0 = Xcsat * dAbulk_dVg + Abulk * dXcsat_dVg;
1510 		      dAbeff_dVg = T0 - model->B3SOIFDadice * dXcsat_dVg;
1511 		      dAbeff_dVb = Xcsat * dAbulk_dVb;
1512 		      dAbeff_dVc = (Abulk - model->B3SOIFDadice) * dXcsat_dVc;
1513                  here->B3SOIFDabeff = Abeff;
1514 
1515 /* Mobility calculation */
1516 		  if (model->B3SOIFDmobMod == 1)
1517 		  {   T0 = Vgsteff + Vth + Vth;
1518 		      T2 = ua + uc * Vbseff;
1519 		      T3 = T0 / model->B3SOIFDtox;
1520 		      T5 = T3 * (T2 + ub * T3);
1521 		      dDenomi_dVg = (T2 + 2.0 * ub * T3) / model->B3SOIFDtox;
1522 		      dDenomi_dVd = dDenomi_dVg * 2 * dVth_dVd;
1523 		      dDenomi_dVb = dDenomi_dVg * 2 * dVth_dVb + uc * T3 ;
1524                       if (selfheat)
1525                          dDenomi_dT = dDenomi_dVg * 2 * dVth_dT
1526                                     + (dua_dT + Vbseff * duc_dT
1527                                     + dub_dT * T3 ) * T3;
1528                       else
1529                          dDenomi_dT = 0.0;
1530 		  }
1531 		  else if (model->B3SOIFDmobMod == 2)
1532 		  {   T5 = Vgsteff / model->B3SOIFDtox * (ua
1533 			 + uc * Vbseff + ub * Vgsteff
1534 			 / model->B3SOIFDtox);
1535 		      dDenomi_dVg = (ua + uc * Vbseff
1536 				  + 2.0 * ub * Vgsteff / model->B3SOIFDtox)
1537 				  / model->B3SOIFDtox;
1538 		      dDenomi_dVd = 0.0;
1539 		      dDenomi_dVb = Vgsteff * uc / model->B3SOIFDtox ;
1540                       if (selfheat)
1541                          dDenomi_dT = Vgsteff / model->B3SOIFDtox
1542                                     * (dua_dT + Vbseff * duc_dT + dub_dT
1543                                     * Vgsteff / model->B3SOIFDtox);
1544                       else
1545                          dDenomi_dT = 0.0;
1546 		  }
1547 		  else  /*  mobMod == 3  */
1548 		  {   T0 = Vgsteff + Vth + Vth;
1549 		      T2 = 1.0 + uc * Vbseff;
1550 		      T3 = T0 / model->B3SOIFDtox;
1551 		      T4 = T3 * (ua + ub * T3);
1552 		      T5 = T4 * T2;
1553 		      dDenomi_dVg = (ua + 2.0 * ub * T3) * T2
1554 				  / model->B3SOIFDtox;
1555 		      dDenomi_dVd = dDenomi_dVg * 2.0 * dVth_dVd;
1556 		      dDenomi_dVb = dDenomi_dVg * 2.0 * dVth_dVb
1557 				  + uc * T4 ;
1558                       if (selfheat)
1559                          dDenomi_dT = dDenomi_dVg * 2.0 * dVth_dT
1560                                     + (dua_dT + dub_dT * T3) * T3 * T2
1561                                     + T4 * Vbseff * duc_dT;
1562                       else
1563                          dDenomi_dT = 0.0;
1564 		  }
1565 
1566 		  if (T5 >= -0.8)
1567 		  {   Denomi = 1.0 + T5;
1568 		  }
1569 		  else /* Added to avoid the discontinuity problem caused by ua and ub*/
1570 		  {   T9 = 1.0 / (7.0 + 10.0 * T5);
1571 		      Denomi = (0.6 + T5) * T9;
1572 		      T9 *= T9;
1573 		      dDenomi_dVg *= T9;
1574 		      dDenomi_dVd *= T9;
1575 		      dDenomi_dVb *= T9;
1576                       if (selfheat)  dDenomi_dT *= T9;
1577                       else   dDenomi_dT = 0.0;
1578 		  }
1579 
1580 		  here->B3SOIFDueff = ueff = u0temp / Denomi;
1581 		  T9 = -ueff / Denomi;
1582 		  dueff_dVg = T9 * dDenomi_dVg;
1583 		  dueff_dVd = T9 * dDenomi_dVd;
1584 		  dueff_dVb = T9 * dDenomi_dVb;
1585                   if (selfheat)  dueff_dT = T9 * dDenomi_dT + du0temp_dT / Denomi;
1586                   else  dueff_dT = 0.0;
1587 
1588 /* Saturation Drain Voltage  Vdsat */
1589 		  WVCox = Weff * vsattemp * model->B3SOIFDcox;
1590 		  WVCoxRds = WVCox * Rds;
1591 
1592 /*                  dWVCoxRds_dT = WVCox * dRds_dT
1593                                  + Weff * model->B3SOIFDcox * Rds * dvsattemp_dT; */
1594 
1595 		  Esat = 2.0 * vsattemp / ueff;
1596 		  EsatL = Esat * Leff;
1597 		  T0 = -EsatL /ueff;
1598 		  dEsatL_dVg = T0 * dueff_dVg;
1599 		  dEsatL_dVd = T0 * dueff_dVd;
1600 		  dEsatL_dVb = T0 * dueff_dVb;
1601                   if (selfheat)
1602                      dEsatL_dT = T0 * dueff_dT + EsatL / vsattemp * dvsattemp_dT;
1603                   else
1604                      dEsatL_dT = 0.0;
1605 
1606 		  /* Sqrt() */
1607 		  a1 = pParam->B3SOIFDa1;
1608 		  if (a1 == 0.0)
1609 		  {   Lambda = pParam->B3SOIFDa2;
1610 		      dLambda_dVg = 0.0;
1611 		  }
1612 		  else if (a1 > 0.0)
1613 /* Added to avoid the discontinuity problem caused by a1 and a2 (Lambda) */
1614 		  {   T0 = 1.0 - pParam->B3SOIFDa2;
1615 		      T1 = T0 - pParam->B3SOIFDa1 * Vgsteff - 0.0001;
1616 		      T2 = sqrt(T1 * T1 + 0.0004 * T0);
1617 		      Lambda = pParam->B3SOIFDa2 + T0 - 0.5 * (T1 + T2);
1618 		      dLambda_dVg = 0.5 * pParam->B3SOIFDa1 * (1.0 + T1 / T2);
1619 		  }
1620 		  else
1621 		  {   T1 = pParam->B3SOIFDa2 + pParam->B3SOIFDa1 * Vgsteff - 0.0001;
1622 		      T2 = sqrt(T1 * T1 + 0.0004 * pParam->B3SOIFDa2);
1623 		      Lambda = 0.5 * (T1 + T2);
1624 		      dLambda_dVg = 0.5 * pParam->B3SOIFDa1 * (1.0 + T1 / T2);
1625 		  }
1626 
1627 		  if (Rds > 0)
1628 		  {   tmp2 = dRds_dVg / Rds + dWeff_dVg / Weff;
1629 		      tmp3 = dRds_dVb / Rds + dWeff_dVb / Weff;
1630 		  }
1631 		  else
1632 		  {   tmp2 = dWeff_dVg / Weff;
1633 		      tmp3 = dWeff_dVb / Weff;
1634 		  }
1635 		  if ((Rds == 0.0) && (Lambda == 1.0))
1636 		  {   T0 = 1.0 / (Abeff * EsatL + Vgst2Vtm);
1637 		      tmp1 = 0.0;
1638 		      T1 = T0 * T0;
1639 		      T2 = Vgst2Vtm * T0;
1640 		      T3 = EsatL * Vgst2Vtm;
1641 		      Vdsat = T3 * T0;
1642 
1643 		      dT0_dVg = -(Abeff * dEsatL_dVg + EsatL * dAbeff_dVg + 1.0) * T1;
1644 		      dT0_dVd = -(Abeff * dEsatL_dVd) * T1;
1645 		      dT0_dVb = -(Abeff * dEsatL_dVb + EsatL * dAbeff_dVb) * T1;
1646                       dT0_dVc = 0.0;
1647                       if (selfheat)
1648 		         dT0_dT  = -(Abeff * dEsatL_dT + dVgst2Vtm_dT) * T1;
1649                       else dT0_dT  = 0.0;
1650 
1651 		      dVdsat_dVg = T3 * dT0_dVg + T2 * dEsatL_dVg + EsatL * T0;
1652 		      dVdsat_dVd = T3 * dT0_dVd + T2 * dEsatL_dVd;
1653 		      dVdsat_dVb = T3 * dT0_dVb + T2 * dEsatL_dVb;
1654 		      dVdsat_dVc = 0.0;
1655                       if (selfheat)
1656 		         dVdsat_dT  = T3 * dT0_dT  + T2 * dEsatL_dT
1657 				    + EsatL * T0 * dVgst2Vtm_dT;
1658                       else dVdsat_dT  = 0.0;
1659 		  }
1660 		  else
1661 		  {   tmp1 = dLambda_dVg / (Lambda * Lambda);
1662 		      T9 = Abeff * WVCoxRds;
1663 		      T8 = Abeff * T9;
1664 		      T7 = Vgst2Vtm * T9;
1665 		      T6 = Vgst2Vtm * WVCoxRds;
1666 		      T0 = 2.0 * Abeff * (T9 - 1.0 + 1.0 / Lambda);
1667 		      dT0_dVg = 2.0 * (T8 * tmp2 - Abeff * tmp1
1668 			      + (2.0 * T9 + 1.0 / Lambda - 1.0) * dAbeff_dVg);
1669 /*		      dT0_dVb = 2.0 * (T8 * tmp3  this is equivalent to one below, but simpler
1670 			      + (2.0 * T9 + 1.0 / Lambda - 1.0) * dAbeff_dVg);  */
1671 		      dT0_dVb = 2.0 * (T8 * (2.0 / Abeff * dAbeff_dVb + tmp3)
1672 			      + (1.0 / Lambda - 1.0) * dAbeff_dVb);
1673 		      dT0_dVd = 0.0;
1674 		      dT0_dVc = 0.0;
1675 
1676                       if (selfheat)
1677                       {
1678 		         tmp4 = dRds_dT / Rds + dvsattemp_dT / vsattemp;
1679 		         dT0_dT  = 2.0 * T8 * tmp4;
1680                       } else tmp4 = dT0_dT = 0.0;
1681 
1682 		      T1 = Vgst2Vtm * (2.0 / Lambda - 1.0) + Abeff * EsatL + 3.0 * T7;
1683 
1684 		      dT1_dVg = (2.0 / Lambda - 1.0) - 2.0 * Vgst2Vtm * tmp1
1685 			      + Abeff * dEsatL_dVg + EsatL * dAbeff_dVg + 3.0 * (T9
1686 			      + T7 * tmp2 + T6 * dAbeff_dVg);
1687 		      dT1_dVb = Abeff * dEsatL_dVb + EsatL * dAbeff_dVb
1688 			      + 3.0 * (T6 * dAbeff_dVb + T7 * tmp3);
1689 		      dT1_dVd = Abeff * dEsatL_dVd;
1690 		      dT1_dVc = 0.0;
1691 
1692 
1693                       if (selfheat)
1694                       {
1695 		         tmp4 += dVgst2Vtm_dT / Vgst2Vtm;
1696 		         dT1_dT  = (2.0 / Lambda - 1.0) * dVgst2Vtm_dT
1697 				 + Abeff * dEsatL_dT + 3.0 * T7 * tmp4;
1698                       } else dT1_dT = 0.0;
1699 
1700 		      T2 = Vgst2Vtm * (EsatL + 2.0 * T6);
1701 		      dT2_dVg = EsatL + Vgst2Vtm * dEsatL_dVg
1702 			      + T6 * (4.0 + 2.0 * Vgst2Vtm * tmp2);
1703 		      dT2_dVb = Vgst2Vtm * (dEsatL_dVb + 2.0 * T6 * tmp3);
1704 		      dT2_dVd = Vgst2Vtm * dEsatL_dVd;
1705                       if (selfheat)
1706 		         dT2_dT  = Vgst2Vtm * dEsatL_dT + EsatL * dVgst2Vtm_dT
1707 				 + 2.0 * T6 * (dVgst2Vtm_dT + Vgst2Vtm * tmp4);
1708                       else
1709 		         dT2_dT  = 0.0;
1710 
1711 		      T3 = sqrt(T1 * T1 - 2.0 * T0 * T2);
1712 		      Vdsat = (T1 - T3) / T0;
1713 
1714 		      dVdsat_dVg = (dT1_dVg - (T1 * dT1_dVg - dT0_dVg * T2
1715 				 - T0 * dT2_dVg) / T3 - Vdsat * dT0_dVg) / T0;
1716 		      dVdsat_dVb = (dT1_dVb - (T1 * dT1_dVb - dT0_dVb * T2
1717 				 - T0 * dT2_dVb) / T3 - Vdsat * dT0_dVb) / T0;
1718 		      dVdsat_dVd = (dT1_dVd - (T1 * dT1_dVd - T0 * dT2_dVd) / T3) / T0;
1719                       dVdsat_dVc  = 0.0;
1720                       if (selfheat)
1721 		         dVdsat_dT  = (dT1_dT - (T1 * dT1_dT - dT0_dT * T2
1722 				    - T0 * dT2_dT) / T3 - Vdsat * dT0_dT) / T0;
1723                       else dVdsat_dT  = 0.0;
1724 		  }
1725 		  here->B3SOIFDvdsat = Vdsat;
1726 
1727 /* Vdsatii for impact ionization */
1728                   if (pParam->B3SOIFDaii > 0.0)
1729                   {
1730                      if (pParam->B3SOIFDcii != 0.0)
1731                      {
1732                         T0 = pParam->B3SOIFDcii / sqrt(3.0) + pParam->B3SOIFDdii;
1733                         /* Hard limit Vds to T0 => T4  i.e. limit T0 to 3.0 */
1734                         T1 = Vds - T0 - 0.1;
1735                         T2 = sqrt(T1 * T1 + 0.4);
1736                         T3 = T0 + 0.5 * (T1 + T2);
1737                         dT3_dVd = 0.5 * (1.0 + T1/T2);
1738 
1739                         T4 = T3 - pParam->B3SOIFDdii;
1740                         T5 = pParam->B3SOIFDcii / T4;
1741                         T0 = T5 * T5;
1742                         dT0_dVd = - 2 * T0 / T4 * dT3_dVd;
1743                      } else
1744                      {
1745                         T0 = dT0_dVd = 0.0;
1746                      }
1747                      T0 += 1.0;
1748 
1749                      T3 = pParam->B3SOIFDaii + pParam->B3SOIFDbii / Leff;
1750                      T4 = 1.0 / (T0 * Vgsteff + T3 * EsatL);
1751                      T5 = -T4 * T4;
1752                      T6 = Vgsteff * T4;
1753                      T7 = EsatL * Vgsteff;
1754                      Vdsatii = T7 * T4;
1755 
1756                      dT4_dVg = T5 * (T0 + T3 * dEsatL_dVg);
1757                      dT4_dVb = T5 * T3 * dEsatL_dVb;
1758                      dT4_dVd = T5 * (Vgsteff * dT0_dVd + T3 * dEsatL_dVd);
1759 
1760                      if (selfheat) dT4_dT = T5 * (T3 * dEsatL_dT);
1761                      else          dT4_dT = 0.0;
1762 
1763                      T8 = T4 * Vgsteff;
1764                      dVdsatii_dVg = T7 * dT4_dVg + T4 * (EsatL + Vgsteff * dEsatL_dVg);
1765                      dVdsatii_dVb = T7 * dT4_dVb + T8 * dEsatL_dVb;
1766                      dVdsatii_dVd = T7 * dT4_dVd + T8 * dEsatL_dVd;
1767                      if (selfheat) dVdsatii_dT  = T7 * dT4_dT  + T8 * dEsatL_dT;
1768                      else          dVdsatii_dT  = 0.0;
1769                   } else
1770                   {
1771                     Vdsatii = Vdsat;
1772                     dVdsatii_dVg = dVdsat_dVg;
1773                     dVdsatii_dVb = dVdsat_dVb;
1774                     dVdsatii_dVd = dVdsat_dVd;
1775                     dVdsatii_dT  = dVdsat_dT;
1776                   }
1777 
1778 /* Effective Vds (Vdseff) Calculation */
1779 		  T1 = Vdsat - Vds - pParam->B3SOIFDdelta;
1780 		  dT1_dVg = dVdsat_dVg;
1781 		  dT1_dVd = dVdsat_dVd - 1.0;
1782 		  dT1_dVb = dVdsat_dVb;
1783 		  dT1_dVc = dVdsat_dVc;
1784 		  dT1_dT  = dVdsat_dT;
1785 
1786 		  T2 = sqrt(T1 * T1 + 4.0 * pParam->B3SOIFDdelta * Vdsat);
1787 		  T0 = T1 / T2;
1788 		  T3 = 2.0 * pParam->B3SOIFDdelta / T2;
1789 		  dT2_dVg = T0 * dT1_dVg + T3 * dVdsat_dVg;
1790 		  dT2_dVd = T0 * dT1_dVd + T3 * dVdsat_dVd;
1791 		  dT2_dVb = T0 * dT1_dVb + T3 * dVdsat_dVb;
1792                   dT2_dVc = 0.0;
1793                   if (selfheat)
1794 		     dT2_dT  = T0 * dT1_dT  + T3 * dVdsat_dT;
1795                   else dT2_dT  = 0.0;
1796 
1797 		  Vdseff = Vdsat - 0.5 * (T1 + T2);
1798 		  dVdseff_dVg = dVdsat_dVg - 0.5 * (dT1_dVg + dT2_dVg);
1799 		  dVdseff_dVd = dVdsat_dVd - 0.5 * (dT1_dVd + dT2_dVd);
1800 		  dVdseff_dVb = dVdsat_dVb - 0.5 * (dT1_dVb + dT2_dVb);
1801                   dVdseff_dVc = 0.0;
1802                   if (selfheat)
1803 		     dVdseff_dT  = dVdsat_dT  - 0.5 * (dT1_dT  + dT2_dT);
1804                   else dVdseff_dT  = 0.0;
1805 
1806 		  if (Vdseff > Vds)
1807 		      Vdseff = Vds; /* This code is added to fixed the problem
1808 				       caused by computer precision when
1809 				       Vds is very close to Vdseff. */
1810 		  diffVds = Vds - Vdseff;
1811 
1812 /* Effective Vdsii for Iii calculation */
1813 		  T1 = Vdsatii - Vds - pParam->B3SOIFDdelta;
1814 
1815 		  T2 = sqrt(T1 * T1 + 4.0 * pParam->B3SOIFDdelta * Vdsatii);
1816 		  T0 = T1 / T2;
1817 		  T3 = 2.0 * pParam->B3SOIFDdelta / T2;
1818                   T4 = T0 + T3;
1819 		  dT2_dVg = T4 * dVdsatii_dVg;
1820 		  dT2_dVd = T4 * dVdsatii_dVd - T0;
1821                   dT2_dVb = T4 * dVdsatii_dVb;
1822                   if (selfheat) dT2_dT = T4*dVdsatii_dT;
1823                   else dT2_dT  = 0.0;
1824 
1825 		  Vdseffii = Vdsatii - 0.5 * (T1 + T2);
1826 		  dVdseffii_dVg = 0.5 * (dVdsatii_dVg - dT2_dVg);
1827 		  dVdseffii_dVd = 0.5 * (dVdsatii_dVd - dT2_dVd + 1.0);
1828 		  dVdseffii_dVb = 0.5 * (dVdsatii_dVb - dT2_dVb);
1829                   if (selfheat)
1830 		     dVdseffii_dT  = 0.5 * (dVdsatii_dT - dT2_dT);
1831                   else dVdseffii_dT  = 0.0;
1832 		  diffVdsii = Vds - Vdseffii;
1833 
1834 /* Calculate VAsat */
1835 		  tmp4 = 1.0 - 0.5 * Abeff * Vdsat / Vgst2Vtm;
1836 		  T9 = WVCoxRds * Vgsteff;
1837 		  T8 = T9 / Vgst2Vtm;
1838 		  T0 = EsatL + Vdsat + 2.0 * T9 * tmp4;
1839 
1840 		  T7 = 2.0 * WVCoxRds * tmp4;
1841 		  dT0_dVg = dEsatL_dVg + dVdsat_dVg + T7 * (1.0 + tmp2 * Vgsteff)
1842                           - T8 * (Abeff * dVdsat_dVg - Abeff * Vdsat / Vgst2Vtm
1843 			  + Vdsat * dAbeff_dVg);
1844 
1845 		  dT0_dVb = dEsatL_dVb + dVdsat_dVb + T7 * tmp3 * Vgsteff
1846 			  - T8 * (dAbeff_dVb * Vdsat + Abeff * dVdsat_dVb);
1847 		  dT0_dVd = dEsatL_dVd + dVdsat_dVd - T8 * Abeff * dVdsat_dVd;
1848                   dT0_dVc = 0.0;
1849 
1850                   if (selfheat)
1851                   {
1852 		     tmp4 = dRds_dT / Rds + dvsattemp_dT / vsattemp;
1853 		     dT0_dT  = dEsatL_dT + dVdsat_dT + T7 * tmp4 * Vgsteff
1854 			     - T8 * (Abeff * dVdsat_dT - Abeff * Vdsat * dVgst2Vtm_dT
1855 			     / Vgst2Vtm);
1856                   } else
1857                      dT0_dT = 0.0;
1858 
1859 		  T9 = WVCoxRds * Abeff;
1860 		  T1 = 2.0 / Lambda - 1.0 + T9;
1861 		  dT1_dVg = -2.0 * tmp1 +  WVCoxRds * (Abeff * tmp2 + dAbeff_dVg);
1862 		  dT1_dVb = dAbeff_dVb * WVCoxRds + T9 * tmp3;
1863                   dT1_dVc = 0.0;
1864                   if (selfheat)
1865 		     dT1_dT  = T9 * tmp4;
1866                   else
1867 		     dT1_dT  = 0.0;
1868 
1869 		  Vasat = T0 / T1;
1870 		  dVasat_dVg = (dT0_dVg - Vasat * dT1_dVg) / T1;
1871 		  dVasat_dVb = (dT0_dVb - Vasat * dT1_dVb) / T1;
1872 		  dVasat_dVd = dT0_dVd / T1;
1873                   dVasat_dVc = 0.0;
1874                   if (selfheat) dVasat_dT  = (dT0_dT  - Vasat * dT1_dT)  / T1;
1875                   else dVasat_dT  = 0.0;
1876 
1877 /* Calculate VACLM */
1878 		  if ((pParam->B3SOIFDpclm > 0.0) && (diffVds > 1.0e-10))
1879 		  {   T0 = 1.0 / (pParam->B3SOIFDpclm * Abeff * pParam->B3SOIFDlitl);
1880 		      dT0_dVb = -T0 / Abeff * dAbeff_dVb;
1881 		      dT0_dVg = -T0 / Abeff * dAbeff_dVg;
1882                       dT0_dVc = 0.0;
1883 
1884 		      T2 = Vgsteff / EsatL;
1885 		      T1 = Leff * (Abeff + T2);
1886 		      dT1_dVg = Leff * ((1.0 - T2 * dEsatL_dVg) / EsatL + dAbeff_dVg);
1887 		      dT1_dVb = Leff * (dAbeff_dVb - T2 * dEsatL_dVb / EsatL);
1888 		      dT1_dVd = -T2 * dEsatL_dVd / Esat;
1889 		      dT1_dVc = 0.0;
1890                       if (selfheat) dT1_dT  = -T2 * dEsatL_dT / Esat;
1891                       else dT1_dT  = 0.0;
1892 
1893 		      T9 = T0 * T1;
1894 		      VACLM = T9 * diffVds;
1895 		      dVACLM_dVg = T0 * dT1_dVg * diffVds - T9 * dVdseff_dVg
1896 				 + T1 * diffVds * dT0_dVg;
1897 		      dVACLM_dVb = (dT0_dVb * T1 + T0 * dT1_dVb) * diffVds
1898 				 - T9 * dVdseff_dVb;
1899 		      dVACLM_dVd = T0 * dT1_dVd * diffVds + T9 * (1.0 - dVdseff_dVd);
1900                       dVACLM_dVc = 0.0;
1901                       if (selfheat)
1902 		         dVACLM_dT  = T0 * dT1_dT * diffVds - T9 * dVdseff_dT;
1903                       else dVACLM_dT  = 0.0;
1904 
1905 		  }
1906 		  else
1907 		  {   VACLM = MAX_EXP;
1908 		      dVACLM_dVd = dVACLM_dVg = dVACLM_dVb = dVACLM_dVc = dVACLM_dT = 0.0;
1909 		  }
1910 
1911 
1912 /* Calculate VADIBL */
1913 		  if (pParam->B3SOIFDthetaRout > 0.0)
1914 		  {   T8 = Abeff * Vdsat;
1915 		      T0 = Vgst2Vtm * T8;
1916 		      T1 = Vgst2Vtm + T8;
1917 		      dT0_dVg = Vgst2Vtm * Abeff * dVdsat_dVg + T8
1918 			      + Vgst2Vtm * Vdsat * dAbeff_dVg;
1919 		      dT1_dVg = 1.0 + Abeff * dVdsat_dVg + Vdsat * dAbeff_dVg;
1920 		      dT1_dVb = dAbeff_dVb * Vdsat + Abeff * dVdsat_dVb;
1921 		      dT0_dVb = Vgst2Vtm * dT1_dVb;
1922 		      dT1_dVd = Abeff * dVdsat_dVd;
1923 		      dT0_dVd = Vgst2Vtm * dT1_dVd;
1924                       dT0_dVc = dT1_dVc = 0.0;
1925                       if (selfheat)
1926                       {
1927 		         dT0_dT  = dVgst2Vtm_dT * T8 + Abeff * Vgst2Vtm * dVdsat_dT;
1928 		         dT1_dT  = dVgst2Vtm_dT + Abeff * dVdsat_dT;
1929                       } else
1930                          dT0_dT = dT1_dT = 0.0;
1931 
1932 		      T9 = T1 * T1;
1933 		      T2 = pParam->B3SOIFDthetaRout;
1934 		      VADIBL = (Vgst2Vtm - T0 / T1) / T2;
1935 		      dVADIBL_dVg = (1.0 - dT0_dVg / T1 + T0 * dT1_dVg / T9) / T2;
1936 		      dVADIBL_dVb = (-dT0_dVb / T1 + T0 * dT1_dVb / T9) / T2;
1937 		      dVADIBL_dVd = (-dT0_dVd / T1 + T0 * dT1_dVd / T9) / T2;
1938                       dVADIBL_dVc = 0.0;
1939                       if (selfheat)
1940 		         dVADIBL_dT = (dVgst2Vtm_dT - dT0_dT/T1 + T0*dT1_dT/T9) / T2;
1941                       else dVADIBL_dT = 0.0;
1942 
1943 		      T7 = pParam->B3SOIFDpdiblb * Vbseff;
1944 		      if (T7 >= -0.9)
1945 		      {   T3 = 1.0 / (1.0 + T7);
1946 			  VADIBL *= T3;
1947 			  dVADIBL_dVg *= T3;
1948 			  dVADIBL_dVb = (dVADIBL_dVb - VADIBL * pParam->B3SOIFDpdiblb)
1949 				      * T3;
1950 			  dVADIBL_dVd *= T3;
1951 			  dVADIBL_dVc *= T3;
1952 			  if (selfheat)  dVADIBL_dT  *= T3;
1953 			  else  dVADIBL_dT  = 0.0;
1954 		      }
1955 		      else
1956 /* Added to avoid the discontinuity problem caused by pdiblcb */
1957 		      {   T4 = 1.0 / (0.8 + T7);
1958 			  T3 = (17.0 + 20.0 * T7) * T4;
1959 			  dVADIBL_dVg *= T3;
1960 			  dVADIBL_dVb = dVADIBL_dVb * T3
1961 				      - VADIBL * pParam->B3SOIFDpdiblb * T4 * T4;
1962 			  dVADIBL_dVd *= T3;
1963 			  dVADIBL_dVc *= T3;
1964 			  if (selfheat)  dVADIBL_dT  *= T3;
1965 			  else  dVADIBL_dT  = 0.0;
1966 			  VADIBL *= T3;
1967 		      }
1968 		  }
1969 		  else
1970 		  {   VADIBL = MAX_EXP;
1971 		      dVADIBL_dVd = dVADIBL_dVg = dVADIBL_dVb = dVADIBL_dVc
1972 			= dVADIBL_dT = 0.0;
1973 		  }
1974 
1975 /* Calculate VA */
1976 
1977 		  T8 = pParam->B3SOIFDpvag / EsatL;
1978 		  T9 = T8 * Vgsteff;
1979 		  if (T9 > -0.9)
1980 		  {   T0 = 1.0 + T9;
1981 		      dT0_dVg = T8 * (1.0 - Vgsteff * dEsatL_dVg / EsatL);
1982 		      dT0_dVb = -T9 * dEsatL_dVb / EsatL;
1983 		      dT0_dVd = -T9 * dEsatL_dVd / EsatL;
1984                       if (selfheat)
1985 		         dT0_dT  = -T9 * dEsatL_dT / EsatL;
1986                       else
1987 		         dT0_dT  = 0.0;
1988 		  }
1989 		  else /* Added to avoid the discontinuity problems caused by pvag */
1990 		  {   T1 = 1.0 / (17.0 + 20.0 * T9);
1991 		      T0 = (0.8 + T9) * T1;
1992 		      T1 *= T1;
1993 		      dT0_dVg = T8 * (1.0 - Vgsteff * dEsatL_dVg / EsatL) * T1;
1994 
1995 		      T9 *= T1 / EsatL;
1996 		      dT0_dVb = -T9 * dEsatL_dVb;
1997 		      dT0_dVd = -T9 * dEsatL_dVd;
1998                       if (selfheat)
1999 		         dT0_dT  = -T9 * dEsatL_dT;
2000                       else
2001 		         dT0_dT  = 0.0;
2002 		  }
2003 
2004 		  tmp1 = VACLM * VACLM;
2005 		  tmp2 = VADIBL * VADIBL;
2006 		  tmp3 = VACLM + VADIBL;
2007 
2008 		  T1 = VACLM * VADIBL / tmp3;
2009 		  tmp3 *= tmp3;
2010 		  dT1_dVg = (tmp1 * dVADIBL_dVg + tmp2 * dVACLM_dVg) / tmp3;
2011 		  dT1_dVd = (tmp1 * dVADIBL_dVd + tmp2 * dVACLM_dVd) / tmp3;
2012 		  dT1_dVb = (tmp1 * dVADIBL_dVb + tmp2 * dVACLM_dVb) / tmp3;
2013                   dT1_dVc = 0.0;
2014                   if (selfheat)
2015 		     dT1_dT  = (tmp1 * dVADIBL_dT  + tmp2 * dVACLM_dT ) / tmp3;
2016                   else dT1_dT  = 0.0;
2017 
2018 		  Va = Vasat + T0 * T1;
2019 		  dVa_dVg = dVasat_dVg + T1 * dT0_dVg + T0 * dT1_dVg;
2020 		  dVa_dVd = dVasat_dVd + T1 * dT0_dVd + T0 * dT1_dVd;
2021 		  dVa_dVb = dVasat_dVb + T1 * dT0_dVb + T0 * dT1_dVb;
2022                   dVa_dVc = 0.0;
2023                   if (selfheat)
2024 		     dVa_dT  = dVasat_dT  + T1 * dT0_dT  + T0 * dT1_dT;
2025                   else dVa_dT  = 0.0;
2026 
2027 /* Calculate Ids */
2028 		  CoxWovL = model->B3SOIFDcox * Weff / Leff;
2029 		  beta = ueff * CoxWovL;
2030 		  dbeta_dVg = CoxWovL * dueff_dVg + beta * dWeff_dVg / Weff;
2031 		  dbeta_dVd = CoxWovL * dueff_dVd;
2032 		  dbeta_dVb = CoxWovL * dueff_dVb + beta * dWeff_dVb / Weff;
2033 		  if (selfheat)  dbeta_dT  = CoxWovL * dueff_dT;
2034 		  else  dbeta_dT  = 0.0;
2035 
2036 		  T0 = 1.0 - 0.5 * Abeff * Vdseff / Vgst2Vtm;
2037 		  dT0_dVg = -0.5 * (Abeff * dVdseff_dVg
2038 			  - Abeff * Vdseff / Vgst2Vtm + Vdseff * dAbeff_dVg) / Vgst2Vtm;
2039 		  dT0_dVd = -0.5 * Abeff * dVdseff_dVd / Vgst2Vtm;
2040 		  dT0_dVb = -0.5 * (Abeff * dVdseff_dVb + dAbeff_dVb * Vdseff)
2041 			  / Vgst2Vtm;
2042 		  dT0_dVc = 0.0;
2043 		  if (selfheat)
2044                      dT0_dT  = -0.5 * (Abeff * dVdseff_dT
2045                              - Abeff * Vdseff / Vgst2Vtm * dVgst2Vtm_dT)
2046 			     / Vgst2Vtm;
2047 		  else dT0_dT = 0.0;
2048 
2049 		  fgche1 = Vgsteff * T0;
2050 		  dfgche1_dVg = Vgsteff * dT0_dVg + T0;
2051 		  dfgche1_dVd = Vgsteff * dT0_dVd;
2052 		  dfgche1_dVb = Vgsteff * dT0_dVb;
2053 		  dfgche1_dVc = 0.0;
2054 		  if (selfheat)  dfgche1_dT  = Vgsteff * dT0_dT;
2055 		  else  dfgche1_dT  = 0.0;
2056 
2057 		  T9 = Vdseff / EsatL;
2058 		  fgche2 = 1.0 + T9;
2059 		  dfgche2_dVg = (dVdseff_dVg - T9 * dEsatL_dVg) / EsatL;
2060 		  dfgche2_dVd = (dVdseff_dVd - T9 * dEsatL_dVd) / EsatL;
2061 		  dfgche2_dVb = (dVdseff_dVb - T9 * dEsatL_dVb) / EsatL;
2062 		  dfgche2_dVc = 0.0;
2063 		  if (selfheat)  dfgche2_dT  = (dVdseff_dT  - T9 * dEsatL_dT)  / EsatL;
2064 		  else  dfgche2_dT  = 0.0;
2065 
2066 		  gche = beta * fgche1 / fgche2;
2067 		  dgche_dVg = (beta * dfgche1_dVg + fgche1 * dbeta_dVg
2068 			    - gche * dfgche2_dVg) / fgche2;
2069 		  dgche_dVd = (beta * dfgche1_dVd + fgche1 * dbeta_dVd
2070 			    - gche * dfgche2_dVd) / fgche2;
2071 		  dgche_dVb = (beta * dfgche1_dVb + fgche1 * dbeta_dVb
2072 			    - gche * dfgche2_dVb) / fgche2;
2073 		  dgche_dVc = 0.0;
2074 		  if (selfheat)
2075 		     dgche_dT  = (beta * dfgche1_dT  + fgche1 * dbeta_dT
2076 			       - gche * dfgche2_dT)  / fgche2;
2077 		  else dgche_dT  = 0.0;
2078 
2079 		  T0 = 1.0 + gche * Rds;
2080 		  T9 = Vdseff / T0;
2081 		  Idl = gche * T9;
2082 
2083 /*  Whoa, these formulas for the derivatives of Idl are convoluted, but I
2084     verified them to be correct  */
2085 
2086 		  dIdl_dVg = (gche * dVdseff_dVg + T9 * dgche_dVg) / T0
2087 			   - Idl * gche / T0 * dRds_dVg ;
2088 		  dIdl_dVd = (gche * dVdseff_dVd + T9 * dgche_dVd) / T0;
2089 		  dIdl_dVb = (gche * dVdseff_dVb + T9 * dgche_dVb
2090 			   - Idl * dRds_dVb * gche) / T0;
2091 		  dIdl_dVc  = 0.0;
2092 		  if (selfheat)
2093 		     dIdl_dT  = (gche * dVdseff_dT + T9 * dgche_dT
2094 			      - Idl * dRds_dT * gche) / T0;
2095 		  else dIdl_dT  = 0.0;
2096 
2097 		  T9 =  diffVds / Va;
2098 		  T0 =  1.0 + T9;
2099 		  here->B3SOIFDids = Ids = Idl * T0;
2100 
2101 		  Gm0 = T0 * dIdl_dVg - Idl * (dVdseff_dVg + T9 * dVa_dVg) / Va;
2102 		  Gds0 = T0 * dIdl_dVd + Idl * (1.0 - dVdseff_dVd
2103 			    - T9 * dVa_dVd) / Va;
2104 		  Gmb0 = T0 * dIdl_dVb - Idl * (dVdseff_dVb + T9 * dVa_dVb) / Va;
2105                   Gmc = 0.0;
2106                   if (selfheat)
2107 		     GmT0 = T0 * dIdl_dT - Idl * (dVdseff_dT + T9 * dVa_dT) / Va;
2108                   else GmT0 = 0.0;
2109 
2110 /* This includes all dependencies from Vgsteff, Vbseff, Vcs */
2111 
2112 		  Gm = Gm0 * dVgsteff_dVg + Gmb0 * dVbseff_dVg + Gmc * dVcs_dVg;
2113 		  Gmb = Gm0 * dVgsteff_dVb + Gmb0 * dVbseff_dVb + Gmc * dVcs_dVb;
2114 		  Gds = Gm0 * dVgsteff_dVd + Gmb0 * dVbseff_dVd + Gmc * dVcs_dVd + Gds0;
2115 		  Gme = Gm0 * dVgsteff_dVe + Gmb0 * dVbseff_dVe + Gmc * dVcs_dVe;
2116 		  if (selfheat)
2117 		     GmT = Gm0 * dVgsteff_dT + Gmb0 * dVbseff_dT  + Gmc * dVcs_dT + GmT0;
2118 		  else GmT = 0.0;
2119 
2120 
2121 /* calculate substrate current Iii */
2122 
2123 		  Giig = Giib = Giid = Giie = GiiT = 0.0;
2124 		  here->B3SOIFDiii = Iii = 0.0;
2125 
2126 		  Idgidl = Gdgidld = Gdgidlg = 0.0;
2127 		  Isgidl = Gsgidlg = 0;
2128 
2129 		  here->B3SOIFDibs = Ibs = 0.0;
2130 		  here->B3SOIFDibd = Ibd = 0.0;
2131 		  here->B3SOIFDic = Ic = 0.0;
2132 
2133                   Gjsb = Gjsd = GjsT = 0.0;
2134                   Gjdb = Gjdd = GjdT = 0.0;
2135                   Gcd = Gcb = GcT = 0.0;
2136 
2137 		  here->B3SOIFDibp = Ibp = 0.0;
2138 		  here->B3SOIFDgbpbs = here->B3SOIFDgbpgs = here->B3SOIFDgbpds = 0.0;
2139 		  here->B3SOIFDgbpes = here->B3SOIFDgbpps = here->B3SOIFDgbpT = here->B3SOIFDcbodcon = 0.0;
2140                   Gbpbs = Gbpgs = Gbpds = Gbpes = Gbpps = GbpT = 0.0;
2141 
2142 		  /*  Current going out of drainprime node into the drain of device  */
2143 		  /*  "node" means the SPICE circuit node  */
2144 
2145                   here->B3SOIFDcdrain = Ids + Ic;
2146                   here->B3SOIFDcd = Ids + Ic - Ibd + Iii + Idgidl;
2147                   here->B3SOIFDcb = Ibs + Ibd + Ibp - Iii - Idgidl - Isgidl;
2148 
2149    		  here->B3SOIFDgds = Gds + Gcd;
2150    		  here->B3SOIFDgm = Gm;
2151    		  here->B3SOIFDgmbs = Gmb + Gcb;
2152 		  here->B3SOIFDgme = Gme;
2153 		  if (selfheat)
2154 		      here->B3SOIFDgmT = GmT + GcT;
2155 		  else
2156 		      here->B3SOIFDgmT = 0.0;
2157 
2158                   /*  note that sign is switched because power flows out
2159                       of device into the temperature node.
2160                       Currently ommit self-heating due to bipolar current
2161                       because it can cause convergence problem*/
2162 
2163                   here->B3SOIFDgtempg = -Gm * Vds;
2164                   here->B3SOIFDgtempb = -Gmb * Vds;
2165                   here->B3SOIFDgtempe = -Gme * Vds;
2166                   here->B3SOIFDgtempT = -GmT * Vds;
2167                   here->B3SOIFDgtempd = -Gds * Vds - Ids;
2168 		  here->B3SOIFDcth = - Ids * Vds - model->B3SOIFDtype *
2169                                   (here->B3SOIFDgtempg * Vgs + here->B3SOIFDgtempb * Vbs
2170                                  + here->B3SOIFDgtempe * Ves + here->B3SOIFDgtempd * Vds)
2171                                  - here->B3SOIFDgtempT * delTemp;
2172 
2173                   /*  Body current which flows into drainprime node from the drain of device  */
2174 
2175 		  here->B3SOIFDgjdb = Gjdb - Giib;
2176 		  here->B3SOIFDgjdd = Gjdd - (Giid + Gdgidld);
2177 		  here->B3SOIFDgjdg = - (Giig + Gdgidlg);
2178 		  here->B3SOIFDgjde = - Giie;
2179 		  if (selfheat) here->B3SOIFDgjdT = GjdT - GiiT;
2180 		  else here->B3SOIFDgjdT = 0.0;
2181 		  here->B3SOIFDcjd = Ibd - Iii - Idgidl - here->B3SOIFDminIsub/2
2182 				 - (here->B3SOIFDgjdb * Vbs + here->B3SOIFDgjdd * Vds
2183 				 +  here->B3SOIFDgjdg * Vgs + here->B3SOIFDgjde * Ves
2184 				 +  here->B3SOIFDgjdT * delTemp);
2185 
2186 		  /*  Body current which flows into sourceprime node from the source of device  */
2187 
2188    		  here->B3SOIFDgjsb = Gjsb;
2189    		  here->B3SOIFDgjsd = Gjsd;
2190 		  here->B3SOIFDgjsg = - Gsgidlg;
2191 		  if (selfheat) here->B3SOIFDgjsT = GjsT;
2192 		  else here->B3SOIFDgjsT = 0.0;
2193 		  here->B3SOIFDcjs = Ibs - Isgidl - here->B3SOIFDminIsub/2
2194 				  - (here->B3SOIFDgjsb * Vbs + here->B3SOIFDgjsd * Vds
2195 				  +  here->B3SOIFDgjsg * Vgs + here->B3SOIFDgjsT * delTemp);
2196 
2197 		  /*  Current flowing into body node  */
2198 
2199 		  here->B3SOIFDgbbs = Giib - Gjsb - Gjdb - Gbpbs;
2200 		  here->B3SOIFDgbgs = Giig + Gdgidlg + Gsgidlg - Gbpgs;
2201    		  here->B3SOIFDgbds = Giid + Gdgidld - Gjsd - Gjdd - Gbpds;
2202 		  here->B3SOIFDgbes = Giie - Gbpes;
2203 		  here->B3SOIFDgbps = - Gbpps;
2204 		  if (selfheat) here->B3SOIFDgbT = GiiT - GjsT - GjdT - GbpT;
2205 		  else here->B3SOIFDgbT = 0.0;
2206 		  here->B3SOIFDcbody = Iii + Idgidl + Isgidl - Ibs - Ibd - Ibp + here->B3SOIFDminIsub
2207 				   - (here->B3SOIFDgbbs * Vbs + here->B3SOIFDgbgs * Vgs
2208 				   + here->B3SOIFDgbds * Vds + here->B3SOIFDgbps * Vps
2209 				   + here->B3SOIFDgbes * Ves + here->B3SOIFDgbT * delTemp);
2210 
2211 	/* Calculate Qinv for Noise analysis */
2212 
2213 		  T1 = Vgsteff * (1.0 - 0.5 * Abeff * Vdseff / Vgst2Vtm);
2214 		  here->B3SOIFDqinv = -model->B3SOIFDcox * pParam->B3SOIFDweff * Leff * T1;
2215 
2216 	/*  Begin CV (charge) model  */
2217 
2218 		  if ((model->B3SOIFDxpart < 0) || (!ChargeComputationNeeded))
2219 		  {   qgate  = qdrn = qsrc = qbody = 0.0;
2220 		      here->B3SOIFDcggb = here->B3SOIFDcgsb = here->B3SOIFDcgdb = 0.0;
2221 		      here->B3SOIFDcdgb = here->B3SOIFDcdsb = here->B3SOIFDcddb = 0.0;
2222 		      here->B3SOIFDcbgb = here->B3SOIFDcbsb = here->B3SOIFDcbdb = 0.0;
2223 		      goto finished;
2224 		  }
2225 		  else
2226 		  {
2227 		       CoxWL  = model->B3SOIFDcox * pParam->B3SOIFDweffCV
2228 			      * pParam->B3SOIFDleffCV;
2229 
2230                        /* By using this Vgsteff,cv, discontinuity in moderate
2231                           inversion charges can be avoid.  However, in capMod=3,
2232                           Vdsat from IV is used.  The dVdsat_dVg is referred to
2233                           the IV Vgsteff and therefore induces error in the charges
2234                           derivatives.  Fortunately, Vgsteff,iv and Vgsteff,cv are
2235                           different only in subthreshold where Qsubs is neglectible.
2236                           So the errors in derivatives is not a serious problem */
2237 
2238 		       if ((VgstNVt > -EXP_THRESHOLD) && (VgstNVt < EXP_THRESHOLD))
2239 		       {   ExpVgst *= ExpVgst;
2240 			   Vgsteff = n * Vtm * log(1.0 + ExpVgst);
2241 			   T0 = ExpVgst / (1.0 + ExpVgst);
2242 			   T1 = -T0 * (dVth_dVb + Vgst / n * dn_dVb) + Vgsteff / n * dn_dVb;
2243 			   dVgsteff_dVd = -T0 * (dVth_dVd + Vgst / n * dn_dVd)
2244 					  + Vgsteff / n * dn_dVd + T1 * dVbseff_dVd;
2245 			   dVgsteff_dVg = T0 * dVgs_eff_dVg + T1 * dVbseff_dVg;
2246                            dVgsteff_dVb = T1 * dVbseff_dVb;
2247                            dVgsteff_dVe = T1 * dVbseff_dVe;
2248                            if (selfheat)
2249 			      dVgsteff_dT  = -T0 * (dVth_dT + Vgst / Temp) + Vgsteff / Temp
2250 					  + T1 * dVbseff_dT;
2251                            else dVgsteff_dT  = 0.0;
2252 		       }
2253 
2254 		       Vfb = Vth - phi - pParam->B3SOIFDk1 * sqrtPhis;
2255 
2256 		       dVfb_dVb = dVth_dVb - pParam->B3SOIFDk1 * dsqrtPhis_dVb;
2257 		       dVfb_dVd = dVth_dVd;
2258 		       dVfb_dT  = dVth_dT;
2259 
2260 		      if ((model->B3SOIFDcapMod == 2) || (model->B3SOIFDcapMod == 3))
2261 		      {
2262 			  /* Necessary because charge behaviour very strange at
2263 			     Vgsteff = 0 */
2264 			  Vgsteff += 1e-4;
2265 
2266 			  /* Something common in capMod 2 and 3 */
2267 			  V3 = Vfb - Vgs_eff + Vbseff - DELTA_3;
2268 			  if (Vfb <= 0.0)
2269 			  {   T0 = sqrt(V3 * V3 - 4.0 * DELTA_3 * Vfb);
2270 			      T2 = -DELTA_3 / T0;
2271 			  }
2272 			  else
2273 			  {   T0 = sqrt(V3 * V3 + 4.0 * DELTA_3 * Vfb);
2274 			      T2 = DELTA_3 / T0;
2275 			  }
2276 
2277 			  T1 = 0.5 * (1.0 + V3 / T0);
2278 			  Vfbeff = Vfb - 0.5 * (V3 + T0);
2279 			  dVfbeff_dVd = (1.0 - T1 - T2) * dVfb_dVd;
2280 			  dVfbeff_dVb = (1.0 - T1 - T2) * dVfb_dVb - T1;
2281 			  dVfbeff_dVrg = T1 * dVgs_eff_dVg;
2282 			  if (selfheat) dVfbeff_dT  = (1.0 - T1 - T2) * dVfb_dT;
2283                           else  dVfbeff_dT = 0.0;
2284 
2285 			  Qac0 = -CoxWL * (Vfbeff - Vfb);
2286 			  dQac0_dVrg = -CoxWL * dVfbeff_dVrg;
2287 			  dQac0_dVd = -CoxWL * (dVfbeff_dVd - dVfb_dVd);
2288 			  dQac0_dVb = -CoxWL * (dVfbeff_dVb - dVfb_dVb);
2289 			  if (selfheat) dQac0_dT  = -CoxWL * (dVfbeff_dT - dVfb_dT);
2290                           else  dQac0_dT = 0.0;
2291 
2292 			  T0 = 0.5 * K1;
2293 			  T3 = Vgs_eff - Vfbeff - Vbseff - Vgsteff;
2294 			  if (pParam->B3SOIFDk1 == 0.0)
2295 			  {   T1 = 0.0;
2296 			      T2 = 0.0;
2297 			  }
2298 			  else if (T3 < 0.0)
2299 			  {   T1 = T0 + T3 / pParam->B3SOIFDk1;
2300 			      T2 = CoxWL;
2301 			  }
2302 			  else
2303 			  {   T1 = sqrt(T0 * T0 + T3);
2304 			      T2 = CoxWL * T0 / T1;
2305 			  }
2306 
2307 			  Qsub0 = CoxWL * K1 * (T0 - T1);
2308 
2309 			  dQsub0_dVrg = T2 * (dVfbeff_dVrg - dVgs_eff_dVg);
2310 			  dQsub0_dVg = T2;
2311 			  dQsub0_dVd = T2 * dVfbeff_dVd;
2312 			  dQsub0_dVb = T2 * (dVfbeff_dVb + 1);
2313 			  if (selfheat) dQsub0_dT  = T2 * dVfbeff_dT;
2314                           else  dQsub0_dT = 0.0;
2315 
2316 			  One_Third_CoxWL = CoxWL / 3.0;
2317 			  Two_Third_CoxWL = 2.0 * One_Third_CoxWL;
2318 			  AbulkCV = Abulk0 * pParam->B3SOIFDabulkCVfactor;
2319 			  dAbulkCV_dVb = pParam->B3SOIFDabulkCVfactor * dAbulk0_dVb;
2320 
2321 			  /*  This is actually capMod=2 calculation  */
2322 			  VdsatCV = Vgsteff / AbulkCV;
2323 			  dVdsatCV_dVg = 1.0 / AbulkCV;
2324 			  dVdsatCV_dVb = -VdsatCV * dAbulkCV_dVb / AbulkCV;
2325                           VdsatCV += 1e-5;
2326 
2327 			  V4 = VdsatCV - Vds - DELTA_4;
2328 			  T0 = sqrt(V4 * V4 + 4.0 * DELTA_4 * VdsatCV);
2329 			  VdseffCV = VdsatCV - 0.5 * (V4 + T0);
2330 			  T1 = 0.5 * (1.0 + V4 / T0);
2331 			  T2 = DELTA_4 / T0;
2332 			  T3 = (1.0 - T1 - T2) / AbulkCV;
2333 			  dVdseffCV_dVg = T3;
2334 			  dVdseffCV_dVd = T1;
2335 			  dVdseffCV_dVb = -T3 * VdsatCV * dAbulkCV_dVb;
2336 
2337 			  if (model->B3SOIFDcapMod == 2)
2338 			  {
2339 			     /* VdsCV Make it compatible with capMod 3 */
2340 			     VdsCV = VdseffCV;
2341 			     dVdsCV_dVg = dVdseffCV_dVg;
2342 			     dVdsCV_dVd = dVdseffCV_dVd;
2343 			     dVdsCV_dVb = dVdseffCV_dVb;
2344 			 }
2345 			 else if (model->B3SOIFDcapMod == 3)
2346 			 {
2347 			    /* Front gate strong inversion depletion charge */
2348 			    /* VdssatCV calculation */
2349 
2350 			     T1 = Vgsteff + K1*sqrtPhis + 0.5*K1*K1;
2351 			     T2 = Vgsteff + K1*sqrtPhis + Phis + 0.25*K1*K1;
2352 
2353 			     dT1_dVb = K1*dsqrtPhis_dVb;
2354 			     dT2_dVb = dT1_dVb + dPhis_dVb;
2355 			     dT1_dVg = dT2_dVg = 1;
2356 
2357                              /* Note VdsatCV is redefined in capMod = 3 */
2358 			     VdsatCV = T1 - K1*sqrt(T2);
2359 
2360 			     dVdsatCV_dVb = dT1_dVb - K1/2/sqrt(T2)*dT2_dVb;
2361 			     dVdsatCV_dVg = dT1_dVg - K1/2/sqrt(T2)*dT2_dVg;
2362 
2363 			     T1 = VdsatCV - Vdsat;
2364 			     dT1_dVg = dVdsatCV_dVg - dVdsat_dVg;
2365 			     dT1_dVb = dVdsatCV_dVb - dVdsat_dVb;
2366 			     dT1_dVd = - dVdsat_dVd;
2367 			     dT1_dVc = - dVdsat_dVc;
2368 			     dT1_dT  = - dVdsat_dT;
2369 
2370 			     if (!(T1 == 0.0))
2371 			     {  T3 = -0.5 * Vdsat / T1;  /* Vdsmax */
2372 				T2 = T3 * Vdsat;
2373 				T4 = T2 + T1 * T3 * T3;     /* fmax */
2374 				if ((Vdseff > T2) && (T1 < 0))
2375 				{
2376 				   VdsCV = T4;
2377 				   T5 = -0.5 / (T1 * T1);
2378 				   dT3_dVg = T5 * (T1 * dVdsat_dVg - Vdsat * dT1_dVg);
2379 				   dT3_dVb = T5 * (T1 * dVdsat_dVb - Vdsat * dT1_dVb);
2380 				   dT3_dVd = T5 * (T1 * dVdsat_dVd - Vdsat * dT1_dVd);
2381                                    dT3_dVc=0.0;
2382                                    if (selfheat)
2383 				      dT3_dT=T5 * (T1 * dVdsat_dT  - Vdsat * dT1_dT);
2384                                    else  dT3_dT=0.0;
2385 
2386 				   dVdsCV_dVd = T3 * dVdsat_dVd + Vdsat * dT3_dVd
2387 					      + T3 * (2 * T1 * dT3_dVd + T3 * dT1_dVd);
2388 				   dVdsCV_dVg = T3 * dVdsat_dVg + Vdsat * dT3_dVg
2389 					      + T3 * (2 * T1 * dT3_dVg + T3 * dT1_dVg);
2390 				   dVdsCV_dVb = T3 * dVdsat_dVb + Vdsat * dT3_dVb
2391 					      + T3 * (2 * T1 * dT3_dVb + T3 * dT1_dVb);
2392                                    dVdsCV_dVc = 0.0;
2393                                    if (selfheat)
2394 				      dVdsCV_dT  = T3 * dVdsat_dT  + Vdsat * dT3_dT
2395 				   	         + T3 * (2 * T1 * dT3_dT  + T3 * dT1_dT );
2396                                    else  dVdsCV_dT = 0.0;
2397 				} else
2398 				{  T5 = Vdseff / Vdsat;
2399 				   T6 = T5 * T5;
2400 				   T7 = 2 * T1 * T5 / Vdsat;
2401 				   T8 = T7 / Vdsat;
2402 				   VdsCV = Vdseff + T1 * T6;
2403 				   dVdsCV_dVd = dVdseff_dVd + T8 *
2404 					      ( Vdsat * dVdseff_dVd
2405 					      - Vdseff * dVdsat_dVd)
2406 					      + T6 * dT1_dVd;
2407 				   dVdsCV_dVb = dVdseff_dVb + T8 * ( Vdsat *
2408 						dVdseff_dVb - Vdseff * dVdsat_dVb)
2409 						+ T6 * dT1_dVb;
2410 				   dVdsCV_dVg = dVdseff_dVg + T8 *
2411 					      ( Vdsat * dVdseff_dVg
2412 					      - Vdseff * dVdsat_dVg)
2413 					      + T6 * dT1_dVg;
2414                                    dVdsCV_dVc = 0.0;
2415 				   if (selfheat)
2416 				      dVdsCV_dT  = dVdseff_dT  + T8 *
2417 					         ( Vdsat * dVdseff_dT
2418 					         - Vdseff * dVdsat_dT )
2419 					         + T6 * dT1_dT ;
2420 				   else dVdsCV_dT  = 0.0;
2421 				}
2422 			     } else
2423 			     {  VdsCV = Vdseff;
2424 				dVdsCV_dVb = dVdseff_dVb;
2425 				dVdsCV_dVd = dVdseff_dVd;
2426 				dVdsCV_dVg = dVdseff_dVg;
2427 				dVdsCV_dVc = dVdseff_dVc;
2428 				dVdsCV_dT  = dVdseff_dT;
2429 			     }
2430 			     if (VdsCV < 0.0) VdsCV = 0.0;
2431 			     VdsCV += 1e-4;
2432 
2433 			     if (VdsCV > (VdsatCV - 1e-7))
2434 			     {
2435 				VdsCV = VdsatCV - 1e-7;
2436 			     }
2437 			     Phisd = Phis + VdsCV;
2438 			     dPhisd_dVb = dPhis_dVb + dVdsCV_dVb;
2439 			     dPhisd_dVd = dVdsCV_dVd;
2440 			     dPhisd_dVg = dVdsCV_dVg;
2441 			     dPhisd_dVc = dVdsCV_dVc;
2442 			     dPhisd_dT  = dVdsCV_dT;
2443 			     sqrtPhisd = sqrt(Phisd);
2444 
2445 			     /* Qdep0 - Depletion charge at Vgs=Vth */
2446 			     T10 = CoxWL * K1;
2447 			     Qdep0 = T10 * sqrtPhis;
2448 			     dQdep0_dVb = T10 * dsqrtPhis_dVb;
2449 			  }    /* End of if capMod == 3 */
2450 
2451 		         /* Something common in both capMod 2 or 3 */
2452 
2453 		         /* Backgate charge */
2454 		         CboxWL = pParam->B3SOIFDkb3 * Cbox * pParam->B3SOIFDweffCV
2455                                 * pParam->B3SOIFDleffCV;
2456 			 Cbg = Cbb = Cbd = Cbe = CbT = 0.0;
2457 			 Ce2g = Ce2b = Ce2d = Ce2e = Ce2T = 0.0;
2458                          Qbf = Qe2 = 0.0;
2459 
2460 			 T2 = - 0.5 * model->B3SOIFDcboxt * pParam->B3SOIFDweffCV
2461 			    * pParam->B3SOIFDleffCV;
2462 			 Qe1 = T2 * VdsCV - CboxWL * (Vbs0eff - Vesfb);
2463 			 Ce1g = T2 * (dVdsCV_dVg * dVgsteff_dVg + dVdsCV_dVb * dVbseff_dVg
2464                                 + dVdsCV_dVc * dVcs_dVg) - CboxWL * dVbs0eff_dVg;
2465 			 Ce1d = T2 * (dVdsCV_dVg * dVgsteff_dVd + dVdsCV_dVb * dVbseff_dVd
2466                                 + dVdsCV_dVc * dVcs_dVd + dVdsCV_dVd)
2467                                 - CboxWL * dVbs0eff_dVd;
2468                          Ce1b = 0.0;
2469 			 Ce1e = T2 * (dVdsCV_dVg * dVgsteff_dVe + dVdsCV_dVb * dVbseff_dVe
2470                                 + dVdsCV_dVc * dVcs_dVe)
2471                                 - CboxWL * (dVbs0eff_dVe - 1.0);
2472                          if (selfheat)
2473                             Ce1T = T2 * (dVdsCV_dVg  * dVgsteff_dT  + dVdsCV_dVb * dVbseff_dT
2474                                    + dVdsCV_dVc * dVcs_dT + dVdsCV_dT)
2475                                    - CboxWL * (dVbs0eff_dT + dvfbb_dT);
2476                          else Ce1T = 0.0;
2477 
2478 			  /* Total inversion charge */
2479 			  T0 = AbulkCV * VdseffCV;
2480 			  T1 = 12.0 * (Vgsteff - 0.5 * T0 + 1e-20);
2481 			  T2 = VdseffCV / T1;
2482 			  T3 = T0 * T2;
2483 
2484 			  T4 = (1.0 - 12.0 * T2 * T2 * AbulkCV);
2485 			  T5 = (6.0 * T0 * (4.0 * Vgsteff - T0) / (T1 * T1) - 0.5);
2486 			  T6 = 12.0 * T2 * T2 * Vgsteff;
2487 
2488 			  qinv = CoxWL * (Vgsteff - 0.5 * VdseffCV + T3);
2489 			  Cgg1 = CoxWL * (T4 + T5 * dVdseffCV_dVg);
2490 			  Cgd1 = CoxWL * T5 * dVdseffCV_dVd;
2491 			  Cgb1 = CoxWL * (T5 * dVdseffCV_dVb + T6 * dAbulkCV_dVb);
2492 
2493 			  /* Inversion charge partitioning into S / D */
2494 			  if (model->B3SOIFDxpart > 0.5)
2495 			  {   /* 0/100 Charge partition model */
2496 			     T1 = T1 + T1;
2497 			     qsrc = -CoxWL * (0.5 * Vgsteff + 0.25 * T0
2498 				  - T0 * T0 / T1);
2499 			     T7 = (4.0 * Vgsteff - T0) / (T1 * T1);
2500 			     T4 = -(0.5 + 24.0 * T0 * T0 / (T1 * T1));
2501 			     T5 = -(0.25 * AbulkCV - 12.0 * AbulkCV * T0 * T7);
2502 			     T6 = -(0.25 * VdseffCV - 12.0 * T0 * VdseffCV * T7);
2503 			     Csg1 = CoxWL * (T4 + T5 * dVdseffCV_dVg);
2504 			     Csd1 = CoxWL * T5 * dVdseffCV_dVd;
2505 			     Csb1 = CoxWL * (T5 * dVdseffCV_dVb + T6 * dAbulkCV_dVb);
2506 
2507 			  }
2508 			  else if (model->B3SOIFDxpart < 0.5)
2509 			  {   /* 40/60 Charge partition model */
2510 			     T1 = T1 / 12.0;
2511 			     T2 = 0.5 * CoxWL / (T1 * T1);
2512 			     T3 = Vgsteff * (2.0 * T0 * T0 / 3.0 + Vgsteff
2513 				* (Vgsteff - 4.0 * T0 / 3.0))
2514 				- 2.0 * T0 * T0 * T0 / 15.0;
2515 			     qsrc = -T2 * T3;
2516 			     T7 = 4.0 / 3.0 * Vgsteff * (Vgsteff - T0)
2517 				+ 0.4 * T0 * T0;
2518 			     T4 = -2.0 * qsrc / T1 - T2 * (Vgsteff * (3.0
2519 				* Vgsteff - 8.0 * T0 / 3.0)
2520 				   + 2.0 * T0 * T0 / 3.0);
2521 			     T5 = (qsrc / T1 + T2 * T7) * AbulkCV;
2522 			     T6 = (qsrc / T1 * VdseffCV + T2 * T7 * VdseffCV);
2523 			     Csg1 = T4 + T5 * dVdseffCV_dVg;
2524 			     Csd1 = T5 * dVdseffCV_dVd;
2525 			     Csb1 = T5 * dVdseffCV_dVb + T6 * dAbulkCV_dVb;
2526 			  }
2527 			  else
2528 			  {   /* 50/50 Charge partition model */
2529 			     qsrc = - 0.5 * qinv;
2530 			     Csg1 = - 0.5 * Cgg1;
2531 			     Csb1 = - 0.5 * Cgb1;
2532 			     Csd1 = - 0.5 * Cgd1;
2533 			  }
2534 
2535 			  Csg = Csg1 * dVgsteff_dVg + Csb1 * dVbseff_dVg;
2536 			  Csd = Csd1 + Csg1 * dVgsteff_dVd + Csb1 * dVbseff_dVd;
2537    			  Csb = Csg1 * dVgsteff_dVb + Csb1 * dVbseff_dVb;
2538 			  Cse = Csg1 * dVgsteff_dVe + Csb1 * dVbseff_dVe;
2539                           if (selfheat)
2540                              CsT = Csg1 * dVgsteff_dT  + Csb1 * dVbseff_dT;
2541                           else  CsT = 0.0;
2542 
2543                           Qex=dQex_dVg=dQex_dVb=dQex_dVd=dQex_dVe=dQex_dT=0.0;
2544 
2545 			  qgate = qinv - (Qbf + Qe2);
2546    			  qbody = (Qbf - Qe1 + Qex);
2547 			  qsub = Qe1 + Qe2 - Qex;
2548 			  qdrn = -(qinv + qsrc);
2549 
2550 			  Cgg = (Cgg1 * dVgsteff_dVg + Cgb1 * dVbseff_dVg) - Cbg ;
2551 			  Cgd = (Cgd1  + Cgg1 * dVgsteff_dVd + Cgb1 * dVbseff_dVd)-Cbd;
2552 			  Cgb = (Cgb1 * dVbseff_dVb + Cgg1 * dVgsteff_dVb) - Cbb;
2553 			  Cge = (Cgg1 * dVgsteff_dVe + Cgb1 * dVbseff_dVe) - Cbe;
2554                           if (selfheat)
2555                              CgT = (Cgg1 * dVgsteff_dT  + Cgb1 * dVbseff_dT ) - CbT;
2556                           else  CgT = 0.0;
2557 
2558 			  here->B3SOIFDcggb = Cgg - Ce2g;
2559 			  here->B3SOIFDcgsb = - (Cgg  + Cgd  + Cgb  + Cge)
2560 					  + (Ce2g + Ce2d + Ce2b + Ce2e);
2561 			  here->B3SOIFDcgdb = Cgd - Ce2d;
2562 			  here->B3SOIFDcgeb = Cge - Ce2e;
2563                           here->B3SOIFDcgT = CgT - Ce2T;
2564 
2565 			  here->B3SOIFDcbgb = Cbg - Ce1g + dQex_dVg;
2566 			  here->B3SOIFDcbsb = -(Cbg  + Cbd  + Cbb  + Cbe)
2567 					  + (Ce1g + Ce1d + Ce1b + Ce1e)
2568                                           - (dQex_dVg + dQex_dVd + dQex_dVb + dQex_dVe);
2569 			  here->B3SOIFDcbdb = Cbd - Ce1d + dQex_dVd;
2570 			  here->B3SOIFDcbeb = Cbe - Ce1e + dQex_dVe;
2571                           here->B3SOIFDcbT = CbT - Ce1T + dQex_dT;
2572 
2573 			  here->B3SOIFDcegb = Ce1g + Ce2g - dQex_dVg;
2574 			  here->B3SOIFDcesb = -(Ce1g + Ce1d + Ce1b + Ce1e)
2575 					  -(Ce2g + Ce2d + Ce2b + Ce2e)
2576                                           +(dQex_dVg + dQex_dVd + dQex_dVb + dQex_dVe);
2577 			  here->B3SOIFDcedb = Ce1d + Ce2d - dQex_dVd;
2578 			  here->B3SOIFDceeb = Ce1e + Ce2e - dQex_dVe;
2579                           here->B3SOIFDceT = Ce1T + Ce2T - dQex_dT;
2580 
2581 			  here->B3SOIFDcdgb = -(Cgg + Cbg + Csg);
2582 			  here->B3SOIFDcddb = -(Cgd + Cbd + Csd);
2583 			  here->B3SOIFDcdeb = -(Cge + Cbe + Cse);
2584                           here->B3SOIFDcdT = -(CgT + CbT + CsT);
2585 			  here->B3SOIFDcdsb = (Cgg + Cgd + Cgb + Cge
2586 					  + Cbg + Cbd + Cbb + Cbe
2587 					  + Csg + Csd + Csb + Cse);
2588 
2589 		      } /* End of if capMod == 2 or capMod ==3 */
2590 		  }
2591 
2592 	finished: /* returning Values to Calling Routine */
2593 		  /*
2594 		   *  COMPUTE EQUIVALENT DRAIN CURRENT SOURCE
2595 		   */
2596 		  if (ChargeComputationNeeded)
2597 		  {
2598                       qjs = qjd = 0.0;
2599                       gcjdds = gcjdbs = gcjdT = 0.0;
2600                       gcjsbs = gcjsT = 0.0;
2601 
2602 		      qdrn -= qjd;
2603 		      qbody += (qjs + qjd);
2604 		      qsrc = -(qgate + qbody + qdrn + qsub);
2605 
2606 		      /* Update the conductance */
2607 		      here->B3SOIFDcddb -= gcjdds;
2608                       here->B3SOIFDcdT -= gcjdT;
2609 		      here->B3SOIFDcdsb += gcjdds + gcjdbs;
2610 
2611 		      here->B3SOIFDcbdb += (gcjdds);
2612                       here->B3SOIFDcbT += (gcjdT + gcjsT);
2613 		      here->B3SOIFDcbsb -= (gcjdds + gcjdbs + gcjsbs);
2614 
2615 		      /* Extrinsic Bottom S/D to substrate charge */
2616 		      T10 = -model->B3SOIFDtype * ves;
2617 		      /* T10 is vse without type conversion */
2618 		      if ( ((pParam->B3SOIFDnsub > 0) && (model->B3SOIFDtype > 0)) ||
2619 		           ((pParam->B3SOIFDnsub < 0) && (model->B3SOIFDtype < 0)) )
2620 		      {
2621 		         if (T10 < pParam->B3SOIFDvsdfb)
2622 		         {  here->B3SOIFDqse = here->B3SOIFDcsbox * (T10 - pParam->B3SOIFDvsdfb);
2623 			    here->B3SOIFDgcse = here->B3SOIFDcsbox;
2624 		         }
2625 		         else if (T10 < pParam->B3SOIFDsdt1)
2626 		         {  T0 = T10 - pParam->B3SOIFDvsdfb;
2627 			    T1 = T0 * T0;
2628 			    here->B3SOIFDqse = T0 * (here->B3SOIFDcsbox -
2629                                              pParam->B3SOIFDst2 / 3 * T1) ;
2630 			    here->B3SOIFDgcse = here->B3SOIFDcsbox - pParam->B3SOIFDst2 * T1;
2631 		         }
2632 		         else if (T10 < pParam->B3SOIFDvsdth)
2633 		         {  T0 = T10 - pParam->B3SOIFDvsdth;
2634 			    T1 = T0 * T0;
2635 			    here->B3SOIFDqse = here->B3SOIFDcsmin * T10 + here->B3SOIFDst4 +
2636                                              pParam->B3SOIFDst3 / 3 * T0 * T1;
2637    			 here->B3SOIFDgcse = here->B3SOIFDcsmin + pParam->B3SOIFDst3 * T1;
2638 		         }
2639 		         else
2640 		         {  here->B3SOIFDqse = here->B3SOIFDcsmin * T10 + here->B3SOIFDst4;
2641 			    here->B3SOIFDgcse = here->B3SOIFDcsmin;
2642 		         }
2643 		      } else
2644 		      {
2645 		         if (T10 < pParam->B3SOIFDvsdth)
2646 		         {  here->B3SOIFDqse = here->B3SOIFDcsmin * (T10 - pParam->B3SOIFDvsdth);
2647 			    here->B3SOIFDgcse = here->B3SOIFDcsmin;
2648 		         }
2649 		         else if (T10 < pParam->B3SOIFDsdt1)
2650 		         {  T0 = T10 - pParam->B3SOIFDvsdth;
2651    			    T1 = T0 * T0;
2652    			    here->B3SOIFDqse = T0 * (here->B3SOIFDcsmin - pParam->B3SOIFDst2 / 3 * T1) ;
2653 			    here->B3SOIFDgcse = here->B3SOIFDcsmin - pParam->B3SOIFDst2 * T1;
2654 		         }
2655 		         else if (T10 < pParam->B3SOIFDvsdfb)
2656 		         {  T0 = T10 - pParam->B3SOIFDvsdfb;
2657 			    T1 = T0 * T0;
2658 			    here->B3SOIFDqse = here->B3SOIFDcsbox * T10 + here->B3SOIFDst4 +
2659                                              pParam->B3SOIFDst3 / 3 * T0 * T1;
2660 			    here->B3SOIFDgcse = here->B3SOIFDcsbox + pParam->B3SOIFDst3 * T1;
2661 		         }
2662 		         else
2663 		         {  here->B3SOIFDqse = here->B3SOIFDcsbox * T10 + here->B3SOIFDst4;
2664 			    here->B3SOIFDgcse = here->B3SOIFDcsbox;
2665 		         }
2666 		      }
2667 
2668 		      /* T11 is vde without type conversion */
2669 		      T11 = model->B3SOIFDtype * (vds - ves);
2670 		      if ( ((pParam->B3SOIFDnsub > 0) && (model->B3SOIFDtype > 0)) ||
2671 		           ((pParam->B3SOIFDnsub < 0) && (model->B3SOIFDtype < 0)) )
2672 		      {
2673 		         if (T11 < pParam->B3SOIFDvsdfb)
2674 		         {  here->B3SOIFDqde = here->B3SOIFDcdbox * (T11 - pParam->B3SOIFDvsdfb);
2675 			    here->B3SOIFDgcde = here->B3SOIFDcdbox;
2676 		         }
2677 		         else if (T11 < pParam->B3SOIFDsdt1)
2678 		         {  T0 = T11 - pParam->B3SOIFDvsdfb;
2679    			    T1 = T0 * T0;
2680    			    here->B3SOIFDqde = T0 * (here->B3SOIFDcdbox - pParam->B3SOIFDdt2 / 3 * T1) ;
2681 			    here->B3SOIFDgcde = here->B3SOIFDcdbox - pParam->B3SOIFDdt2 * T1;
2682 		         }
2683 		         else if (T11 < pParam->B3SOIFDvsdth)
2684 		         {  T0 = T11 - pParam->B3SOIFDvsdth;
2685 			    T1 = T0 * T0;
2686 			    here->B3SOIFDqde = here->B3SOIFDcdmin * T11 + here->B3SOIFDdt4 +
2687                                              pParam->B3SOIFDdt3 / 3 * T0 * T1;
2688 			    here->B3SOIFDgcde = here->B3SOIFDcdmin + pParam->B3SOIFDdt3 * T1;
2689 		         }
2690 		         else
2691 		         {  here->B3SOIFDqde = here->B3SOIFDcdmin * T11 + here->B3SOIFDdt4;
2692 			    here->B3SOIFDgcde = here->B3SOIFDcdmin;
2693 		         }
2694 		      } else
2695 		      {
2696 		         if (T11 < pParam->B3SOIFDvsdth)
2697 		         {  here->B3SOIFDqde = here->B3SOIFDcdmin * (T11 - pParam->B3SOIFDvsdth);
2698 			    here->B3SOIFDgcde = here->B3SOIFDcdmin;
2699 		         }
2700 		         else if (T11 < pParam->B3SOIFDsdt1)
2701 		         {  T0 = T11 - pParam->B3SOIFDvsdth;
2702    			    T1 = T0 * T0;
2703    			    here->B3SOIFDqde = T0 * (here->B3SOIFDcdmin - pParam->B3SOIFDdt2 / 3 * T1) ;
2704 			    here->B3SOIFDgcde = here->B3SOIFDcdmin - pParam->B3SOIFDdt2 * T1;
2705 		         }
2706 		         else if (T11 < pParam->B3SOIFDvsdfb)
2707 		         {  T0 = T11 - pParam->B3SOIFDvsdfb;
2708 			    T1 = T0 * T0;
2709 			    here->B3SOIFDqde = here->B3SOIFDcdbox * T11 + here->B3SOIFDdt4 +
2710                                              pParam->B3SOIFDdt3 / 3 * T0 * T1;
2711 			    here->B3SOIFDgcde = here->B3SOIFDcdbox + pParam->B3SOIFDdt3 * T1;
2712 		         }
2713 		         else
2714 		         {  here->B3SOIFDqde = here->B3SOIFDcdbox * T11 + here->B3SOIFDdt4;
2715 			    here->B3SOIFDgcde = here->B3SOIFDcdbox;
2716 		         }
2717 		      }
2718 
2719 		      /* Extrinsic : Sidewall fringing S/D charge */
2720 		      here->B3SOIFDqse += pParam->B3SOIFDcsesw * T10;
2721 		      here->B3SOIFDgcse += pParam->B3SOIFDcsesw;
2722 		      here->B3SOIFDqde += pParam->B3SOIFDcdesw * T11;
2723 		      here->B3SOIFDgcde += pParam->B3SOIFDcdesw;
2724 
2725 		      /* All charge are mutliplied with type at the end, but qse and qde
2726 			 have true polarity => so pre-mutliplied with type */
2727 		      here->B3SOIFDqse *= model->B3SOIFDtype;
2728 		      here->B3SOIFDqde *= model->B3SOIFDtype;
2729 		  }
2730 
2731 		  here->B3SOIFDxc = Xc;
2732 		  here->B3SOIFDcbb = Cbb;
2733 		  here->B3SOIFDcbd = Cbd;
2734 		  here->B3SOIFDcbg = Cbg;
2735 		  here->B3SOIFDqbf = Qbf;
2736 		  here->B3SOIFDqjs = qjs;
2737 		  here->B3SOIFDqjd = qjd;
2738 
2739                   if (here->B3SOIFDdebugMod == -1)
2740                      ChargeComputationNeeded = 0;
2741 
2742 		  /*
2743 		   *  check convergence
2744 		   */
2745 		  if ((here->B3SOIFDoff == 0) || (!(ckt->CKTmode & MODEINITFIX)))
2746 		  {   if (Check == 1)
2747 		      {   ckt->CKTnoncon++;
2748 if (here->B3SOIFDdebugMod > 2)
2749    fprintf(fpdebug, "Check is on, noncon=%d\n", ckt->CKTnoncon++);
2750 
2751 		      }
2752 		  }
2753 
2754                   *(ckt->CKTstate0 + here->B3SOIFDvg) = vg;
2755                   *(ckt->CKTstate0 + here->B3SOIFDvd) = vd;
2756                   *(ckt->CKTstate0 + here->B3SOIFDvs) = vs;
2757                   *(ckt->CKTstate0 + here->B3SOIFDvp) = vp;
2758                   *(ckt->CKTstate0 + here->B3SOIFDve) = ve;
2759 
2760 		  *(ckt->CKTstate0 + here->B3SOIFDvbs) = vbs;
2761 		  *(ckt->CKTstate0 + here->B3SOIFDvbd) = vbd;
2762 		  *(ckt->CKTstate0 + here->B3SOIFDvgs) = vgs;
2763 		  *(ckt->CKTstate0 + here->B3SOIFDvds) = vds;
2764 		  *(ckt->CKTstate0 + here->B3SOIFDves) = ves;
2765 		  *(ckt->CKTstate0 + here->B3SOIFDvps) = vps;
2766 		  *(ckt->CKTstate0 + here->B3SOIFDdeltemp) = delTemp;
2767 
2768 		  /* bulk and channel charge plus overlaps */
2769 
2770 		  if (!ChargeComputationNeeded)
2771 		      goto line850;
2772 #ifndef NOBYPASS
2773 	line755:
2774 #endif
2775 		  ag0 = ckt->CKTag[0];
2776 
2777 		  T0 = vgd + DELTA_1;
2778 		  T1 = sqrt(T0 * T0 + 4.0 * DELTA_1);
2779 		  T2 = 0.5 * (T0 - T1);
2780 
2781 		  T3 = pParam->B3SOIFDweffCV * pParam->B3SOIFDcgdl;
2782 		  T4 = sqrt(1.0 - 4.0 * T2 / pParam->B3SOIFDckappa);
2783 		  cgdo = pParam->B3SOIFDcgdo + T3 - T3 * (1.0 - 1.0 / T4)
2784 			 * (0.5 - 0.5 * T0 / T1);
2785 		  qgdo = (pParam->B3SOIFDcgdo + T3) * vgd - T3 * (T2
2786 			 + 0.5 * pParam->B3SOIFDckappa * (T4 - 1.0));
2787 
2788 		  T0 = vgs + DELTA_1;
2789 		  T1 = sqrt(T0 * T0 + 4.0 * DELTA_1);
2790 		  T2 = 0.5 * (T0 - T1);
2791 		  T3 = pParam->B3SOIFDweffCV * pParam->B3SOIFDcgsl;
2792 		  T4 = sqrt(1.0 - 4.0 * T2 / pParam->B3SOIFDckappa);
2793 		  cgso = pParam->B3SOIFDcgso + T3 - T3 * (1.0 - 1.0 / T4)
2794 			 * (0.5 - 0.5 * T0 / T1);
2795 		  qgso = (pParam->B3SOIFDcgso + T3) * vgs - T3 * (T2
2796 			 + 0.5 * pParam->B3SOIFDckappa * (T4 - 1.0));
2797 
2798 
2799 		  if (here->B3SOIFDmode > 0)
2800 		  {   gcdgb = (here->B3SOIFDcdgb - cgdo) * ag0;
2801 		      gcddb = (here->B3SOIFDcddb + cgdo + here->B3SOIFDgcde) * ag0;
2802 		      gcdsb = here->B3SOIFDcdsb * ag0;
2803 		      gcdeb = (here->B3SOIFDcdeb - here->B3SOIFDgcde) * ag0;
2804                       gcdT = model->B3SOIFDtype * here->B3SOIFDcdT * ag0;
2805 
2806 		      gcsgb = -(here->B3SOIFDcggb + here->B3SOIFDcbgb + here->B3SOIFDcdgb
2807 			    + here->B3SOIFDcegb + cgso) * ag0;
2808 		      gcsdb = -(here->B3SOIFDcgdb + here->B3SOIFDcbdb + here->B3SOIFDcddb
2809 			    + here->B3SOIFDcedb) * ag0;
2810 		      gcssb = (cgso + here->B3SOIFDgcse - (here->B3SOIFDcgsb + here->B3SOIFDcbsb
2811 			    + here->B3SOIFDcdsb + here->B3SOIFDcesb)) * ag0;
2812 		      gcseb = -(here->B3SOIFDgcse + here->B3SOIFDcgeb + here->B3SOIFDcbeb + here->B3SOIFDcdeb
2813 			    + here->B3SOIFDceeb) * ag0;
2814                       gcsT = - model->B3SOIFDtype * (here->B3SOIFDcgT + here->B3SOIFDcbT + here->B3SOIFDcdT
2815                             + here->B3SOIFDceT) * ag0;
2816 
2817 		      gcggb = (here->B3SOIFDcggb + cgdo + cgso + pParam->B3SOIFDcgeo) * ag0;
2818 		      gcgdb = (here->B3SOIFDcgdb - cgdo) * ag0;
2819 		      gcgsb = (here->B3SOIFDcgsb - cgso) * ag0;
2820 		      gcgeb = (here->B3SOIFDcgeb - pParam->B3SOIFDcgeo) * ag0;
2821                       gcgT = model->B3SOIFDtype * here->B3SOIFDcgT * ag0;
2822 
2823 		      gcbgb = here->B3SOIFDcbgb * ag0;
2824 		      gcbdb = here->B3SOIFDcbdb * ag0;
2825 		      gcbsb = here->B3SOIFDcbsb * ag0;
2826    		      gcbeb = here->B3SOIFDcbeb * ag0;
2827                       gcbT = model->B3SOIFDtype * here->B3SOIFDcbT * ag0;
2828 
2829 		      gcegb = (here->B3SOIFDcegb - pParam->B3SOIFDcgeo) * ag0;
2830 		      gcedb = (here->B3SOIFDcedb - here->B3SOIFDgcde) * ag0;
2831 		      gcesb = (here->B3SOIFDcesb - here->B3SOIFDgcse) * ag0;
2832 		      gceeb = (here->B3SOIFDgcse + here->B3SOIFDgcde +
2833                                here->B3SOIFDceeb + pParam->B3SOIFDcgeo) * ag0;
2834 
2835                       gceT = model->B3SOIFDtype * here->B3SOIFDceT * ag0;
2836 
2837                       gcTt = pParam->B3SOIFDcth * ag0;
2838 
2839 		      sxpart = 0.6;
2840 		      dxpart = 0.4;
2841 
2842 		      /* Lump the overlap capacitance and S/D parasitics */
2843 		      qgd = qgdo;
2844 		      qgs = qgso;
2845 		      qge = pParam->B3SOIFDcgeo * vge;
2846 		      qgate += qgd + qgs + qge;
2847 		      qdrn += here->B3SOIFDqde - qgd;
2848 		      qsub -= qge + here->B3SOIFDqse + here->B3SOIFDqde;
2849 		      qsrc = -(qgate + qbody + qdrn + qsub);
2850 		  }
2851 		  else
2852 		  {   gcsgb = (here->B3SOIFDcdgb - cgso) * ag0;
2853 		      gcssb = (here->B3SOIFDcddb + cgso + here->B3SOIFDgcse) * ag0;
2854 		      gcsdb = here->B3SOIFDcdsb * ag0;
2855 		      gcseb = (here->B3SOIFDcdeb - here->B3SOIFDgcse) * ag0;
2856                       gcsT = model->B3SOIFDtype * here->B3SOIFDcdT * ag0;
2857 
2858 		      gcdgb = -(here->B3SOIFDcggb + here->B3SOIFDcbgb + here->B3SOIFDcdgb
2859 			    + here->B3SOIFDcegb + cgdo) * ag0;
2860 		      gcdsb = -(here->B3SOIFDcgdb + here->B3SOIFDcbdb + here->B3SOIFDcddb
2861 			    + here->B3SOIFDcedb) * ag0;
2862 		      gcddb = (cgdo + here->B3SOIFDgcde - (here->B3SOIFDcgsb + here->B3SOIFDcbsb
2863 			    + here->B3SOIFDcdsb + here->B3SOIFDcesb)) * ag0;
2864 		      gcdeb = -(here->B3SOIFDgcde + here->B3SOIFDcgeb + here->B3SOIFDcbeb + here->B3SOIFDcdeb
2865 			    + here->B3SOIFDceeb) * ag0;
2866                       gcdT = - model->B3SOIFDtype * (here->B3SOIFDcgT + here->B3SOIFDcbT
2867                              + here->B3SOIFDcdT + here->B3SOIFDceT) * ag0;
2868 
2869 		      gcggb = (here->B3SOIFDcggb + cgdo + cgso + pParam->B3SOIFDcgeo) * ag0;
2870 		      gcgsb = (here->B3SOIFDcgdb - cgso) * ag0;
2871 		      gcgdb = (here->B3SOIFDcgsb - cgdo) * ag0;
2872 		      gcgeb = (here->B3SOIFDcgeb - pParam->B3SOIFDcgeo) * ag0;
2873                       gcgT = model->B3SOIFDtype * here->B3SOIFDcgT * ag0;
2874 
2875 		      gcbgb = here->B3SOIFDcbgb * ag0;
2876 		      gcbsb = here->B3SOIFDcbdb * ag0;
2877 		      gcbdb = here->B3SOIFDcbsb * ag0;
2878 		      gcbeb = here->B3SOIFDcbeb * ag0;
2879                       gcbT = model->B3SOIFDtype * here->B3SOIFDcbT * ag0;
2880 
2881 		      gcegb = (here->B3SOIFDcegb - pParam->B3SOIFDcgeo) * ag0;
2882 		      gcesb = (here->B3SOIFDcedb - here->B3SOIFDgcse) * ag0;
2883 		      gcedb = (here->B3SOIFDcesb - here->B3SOIFDgcde) * ag0;
2884 		      gceeb = (here->B3SOIFDceeb + pParam->B3SOIFDcgeo +
2885 			       here->B3SOIFDgcse + here->B3SOIFDgcde) * ag0;
2886                       gceT = model->B3SOIFDtype * here->B3SOIFDceT * ag0;
2887 
2888                       gcTt = pParam->B3SOIFDcth * ag0;
2889 
2890 		      dxpart = 0.6;
2891 		      sxpart = 0.4;
2892 
2893 		      /* Lump the overlap capacitance */
2894 		      qgd = qgdo;
2895 		      qgs = qgso;
2896 		      qge = pParam->B3SOIFDcgeo * vge;
2897 		      qgate += qgd + qgs + qge;
2898 		      qsrc = qdrn - qgs + here->B3SOIFDqse;
2899 		      qsub -= qge + here->B3SOIFDqse + here->B3SOIFDqde;
2900 		      qdrn = -(qgate + qbody + qsrc + qsub);
2901 		  }
2902 
2903 		  here->B3SOIFDcgdo = cgdo;
2904 		  here->B3SOIFDcgso = cgso;
2905 
2906                   if (ByPass) goto line860;
2907 
2908 		  *(ckt->CKTstate0 + here->B3SOIFDqe) = qsub;
2909 		  *(ckt->CKTstate0 + here->B3SOIFDqg) = qgate;
2910 		  *(ckt->CKTstate0 + here->B3SOIFDqd) = qdrn;
2911                   *(ckt->CKTstate0 + here->B3SOIFDqb) = qbody;
2912 		  if ((model->B3SOIFDshMod == 1) && (here->B3SOIFDrth0!=0.0))
2913                      *(ckt->CKTstate0 + here->B3SOIFDqth) = pParam->B3SOIFDcth * delTemp;
2914 
2915 
2916 		  /* store small signal parameters */
2917 		  if (ckt->CKTmode & MODEINITSMSIG)
2918 		  {   goto line1000;
2919 		  }
2920 		  if (!ChargeComputationNeeded)
2921 		      goto line850;
2922 
2923 
2924 		  if (ckt->CKTmode & MODEINITTRAN)
2925 		  {   *(ckt->CKTstate1 + here->B3SOIFDqb) =
2926 			    *(ckt->CKTstate0 + here->B3SOIFDqb);
2927 		      *(ckt->CKTstate1 + here->B3SOIFDqg) =
2928 			    *(ckt->CKTstate0 + here->B3SOIFDqg);
2929 		      *(ckt->CKTstate1 + here->B3SOIFDqd) =
2930 			    *(ckt->CKTstate0 + here->B3SOIFDqd);
2931 		      *(ckt->CKTstate1 + here->B3SOIFDqe) =
2932 			    *(ckt->CKTstate0 + here->B3SOIFDqe);
2933                       *(ckt->CKTstate1 + here->B3SOIFDqth) =
2934                             *(ckt->CKTstate0 + here->B3SOIFDqth);
2935 		  }
2936 
2937                   error = NIintegrate(ckt, &geq, &ceq,0.0,here->B3SOIFDqb);
2938                   if (error) return(error);
2939                   error = NIintegrate(ckt, &geq, &ceq, 0.0, here->B3SOIFDqg);
2940                   if (error) return(error);
2941                   error = NIintegrate(ckt,&geq, &ceq, 0.0, here->B3SOIFDqd);
2942                   if (error) return(error);
2943                   error = NIintegrate(ckt,&geq, &ceq, 0.0, here->B3SOIFDqe);
2944                   if (error) return(error);
2945 		  if ((model->B3SOIFDshMod == 1) && (here->B3SOIFDrth0!=0.0))
2946                   {
2947                       error = NIintegrate(ckt, &geq, &ceq, 0.0, here->B3SOIFDqth);
2948                       if (error) return (error);
2949                   }
2950 
2951 		  goto line860;
2952 
2953 	line850:
2954 		  /* initialize to zero charge conductance and current */
2955 		  ceqqe = ceqqg = ceqqb = ceqqd = ceqqth= 0.0;
2956 
2957 		  gcdgb = gcddb = gcdsb = gcdeb = gcdT = 0.0;
2958 		  gcsgb = gcsdb = gcssb = gcseb = gcsT = 0.0;
2959 		  gcggb = gcgdb = gcgsb = gcgeb = gcgT = 0.0;
2960 		  gcbgb = gcbdb = gcbsb = gcbeb = gcbT = 0.0;
2961 		  gcegb = gcedb = gceeb = gcesb = gceT = 0.0;
2962                   gcTt = 0.0;
2963 
2964 		  sxpart = (1.0 - (dxpart = (here->B3SOIFDmode > 0) ? 0.4 : 0.6));
2965 
2966 		  goto line900;
2967 
2968 	line860:
2969 		  /* evaluate equivalent charge current */
2970 
2971                   cqgate = *(ckt->CKTstate0 + here->B3SOIFDcqg);
2972                   cqbody = *(ckt->CKTstate0 + here->B3SOIFDcqb);
2973                   cqdrn = *(ckt->CKTstate0 + here->B3SOIFDcqd);
2974                   cqsub = *(ckt->CKTstate0 + here->B3SOIFDcqe);
2975                   cqtemp = *(ckt->CKTstate0 + here->B3SOIFDcqth);
2976 
2977                   here->B3SOIFDcb += cqbody;
2978                   here->B3SOIFDcd += cqdrn;
2979 
2980 		  ceqqg = cqgate - gcggb * vgb + gcgdb * vbd + gcgsb * vbs
2981                           - gcgeb * veb - gcgT * delTemp;
2982 		  ceqqb = cqbody - gcbgb * vgb + gcbdb * vbd + gcbsb * vbs
2983 			  - gcbeb * veb - gcbT * delTemp;
2984 		  ceqqd = cqdrn - gcdgb * vgb + gcddb * vbd + gcdsb * vbs
2985                           - gcdeb * veb - gcdT * delTemp;
2986 		  ceqqe = cqsub - gcegb * vgb + gcedb * vbd + gcesb * vbs
2987 			  - gceeb * veb - gceT * delTemp;;
2988                   ceqqth = cqtemp - gcTt * delTemp;
2989 
2990 		  if (ckt->CKTmode & MODEINITTRAN)
2991 		  {   *(ckt->CKTstate1 + here->B3SOIFDcqe) =
2992 			    *(ckt->CKTstate0 + here->B3SOIFDcqe);
2993 		      *(ckt->CKTstate1 + here->B3SOIFDcqb) =
2994 			    *(ckt->CKTstate0 + here->B3SOIFDcqb);
2995 		      *(ckt->CKTstate1 + here->B3SOIFDcqg) =
2996 			    *(ckt->CKTstate0 + here->B3SOIFDcqg);
2997 		      *(ckt->CKTstate1 + here->B3SOIFDcqd) =
2998 			    *(ckt->CKTstate0 + here->B3SOIFDcqd);
2999 		      *(ckt->CKTstate1 + here->B3SOIFDcqth) =
3000 			    *(ckt->CKTstate0 + here->B3SOIFDcqth);
3001 		  }
3002 
3003 		  /*
3004 		   *  load current vector
3005 		   */
3006 	line900:
3007 
3008                   m = here->B3SOIFDm;
3009 
3010 		  if (here->B3SOIFDmode >= 0)
3011 		  {   Gm = here->B3SOIFDgm;
3012 		      Gmbs = here->B3SOIFDgmbs;
3013 		      Gme = here->B3SOIFDgme;
3014 		      GmT = model->B3SOIFDtype * here->B3SOIFDgmT;
3015 		      FwdSum = Gm + Gmbs + Gme;
3016 		      RevSum = 0.0;
3017 		      cdreq = model->B3SOIFDtype * (here->B3SOIFDcdrain - here->B3SOIFDgds * vds
3018 			    - Gm * vgs - Gmbs * vbs - Gme * ves - GmT * delTemp);
3019 		      /* ceqbs now is compatible with cdreq, ie. going in is +ve */
3020 		      /* Equivalent current source from the diode */
3021 		      ceqbs = here->B3SOIFDcjs;
3022 		      ceqbd = here->B3SOIFDcjd;
3023 		      /* Current going in is +ve */
3024 		      ceqbody = -here->B3SOIFDcbody;
3025 		      ceqth = here->B3SOIFDcth;
3026 		      ceqbodcon = here->B3SOIFDcbodcon;
3027 
3028 		      gbbg  = -here->B3SOIFDgbgs;
3029 		      gbbdp = -here->B3SOIFDgbds;
3030 		      gbbb  = -here->B3SOIFDgbbs;
3031 		      gbbe  = -here->B3SOIFDgbes;
3032 		      gbbp  = -here->B3SOIFDgbps;
3033 		      gbbT  = -model->B3SOIFDtype * here->B3SOIFDgbT;
3034 		      gbbsp = - ( gbbg + gbbdp + gbbb + gbbe + gbbp);
3035 
3036 		      gddpg  = -here->B3SOIFDgjdg;
3037 		      gddpdp = -here->B3SOIFDgjdd;
3038 		      gddpb  = -here->B3SOIFDgjdb;
3039 		      gddpe  = -here->B3SOIFDgjde;
3040 		      gddpT  = -model->B3SOIFDtype * here->B3SOIFDgjdT;
3041 		      gddpsp = - ( gddpg + gddpdp + gddpb + gddpe);
3042 
3043 		      gsspg  = -here->B3SOIFDgjsg;
3044 		      gsspdp = -here->B3SOIFDgjsd;
3045 		      gsspb  = -here->B3SOIFDgjsb;
3046 		      gsspe  = 0.0;
3047 		      gsspT  = -model->B3SOIFDtype * here->B3SOIFDgjsT;
3048 		      gsspsp = - (gsspg + gsspdp + gsspb + gsspe);
3049 
3050 		      gppg = -here->B3SOIFDgbpgs;
3051 		      gppdp = -here->B3SOIFDgbpds;
3052 		      gppb = -here->B3SOIFDgbpbs;
3053 		      gppe = -here->B3SOIFDgbpes;
3054 		      gppp = -here->B3SOIFDgbpps;
3055 		      gppT = -model->B3SOIFDtype * here->B3SOIFDgbpT;
3056 		      gppsp = - (gppg + gppdp + gppb + gppe + gppp);
3057 
3058                       gTtg  = here->B3SOIFDgtempg;
3059                       gTtb  = here->B3SOIFDgtempb;
3060                       gTte  = here->B3SOIFDgtempe;
3061                       gTtdp = here->B3SOIFDgtempd;
3062                       gTtt  = here->B3SOIFDgtempT;
3063                       gTtsp = - (gTtg + gTtb + gTte + gTtdp);
3064 		  }
3065 		  else
3066 		  {   Gm = -here->B3SOIFDgm;
3067 		      Gmbs = -here->B3SOIFDgmbs;
3068 		      Gme = -here->B3SOIFDgme;
3069 		      GmT = -model->B3SOIFDtype * here->B3SOIFDgmT;
3070 		      FwdSum = 0.0;
3071 		      RevSum = -(Gm + Gmbs + Gme);
3072 		      cdreq = -model->B3SOIFDtype * (here->B3SOIFDcdrain + here->B3SOIFDgds*vds
3073 			    + Gm * vgd + Gmbs * vbd + Gme * (ves - vds) + GmT * delTemp);
3074 		      ceqbs = here->B3SOIFDcjd;
3075 		      ceqbd = here->B3SOIFDcjs;
3076 		      /* Current going in is +ve */
3077 		      ceqbody = -here->B3SOIFDcbody;
3078 		      ceqth = here->B3SOIFDcth;
3079 		      ceqbodcon = here->B3SOIFDcbodcon;
3080 
3081 		      gbbg  = -here->B3SOIFDgbgs;
3082 		      gbbb  = -here->B3SOIFDgbbs;
3083 		      gbbe  = -here->B3SOIFDgbes;
3084 		      gbbp  = -here->B3SOIFDgbps;
3085 		      gbbsp = -here->B3SOIFDgbds;
3086 		      gbbT  = -model->B3SOIFDtype * here->B3SOIFDgbT;
3087 		      gbbdp = - ( gbbg + gbbsp + gbbb + gbbe + gbbp);
3088 
3089 		      gddpg  = -here->B3SOIFDgjsg;
3090 		      gddpsp = -here->B3SOIFDgjsd;
3091 		      gddpb  = -here->B3SOIFDgjsb;
3092 		      gddpe  = 0.0;
3093 		      gddpT  = -model->B3SOIFDtype * here->B3SOIFDgjsT;
3094 		      gddpdp = - (gddpg + gddpsp + gddpb + gddpe);
3095 
3096 		      gsspg  = -here->B3SOIFDgjdg;
3097 		      gsspsp = -here->B3SOIFDgjdd;
3098 		      gsspb  = -here->B3SOIFDgjdb;
3099 		      gsspe  = -here->B3SOIFDgjde;
3100 		      gsspT  = -model->B3SOIFDtype * here->B3SOIFDgjdT;
3101 		      gsspdp = - ( gsspg + gsspsp + gsspb + gsspe);
3102 
3103 		      gppg = -here->B3SOIFDgbpgs;
3104 		      gppsp = -here->B3SOIFDgbpds;
3105 		      gppb = -here->B3SOIFDgbpbs;
3106 		      gppe = -here->B3SOIFDgbpes;
3107 		      gppp = -here->B3SOIFDgbpps;
3108 		      gppT = -model->B3SOIFDtype * here->B3SOIFDgbpT;
3109 		      gppdp = - (gppg + gppsp + gppb + gppe + gppp);
3110 
3111                       gTtg  = here->B3SOIFDgtempg;
3112                       gTtb  = here->B3SOIFDgtempb;
3113                       gTte  = here->B3SOIFDgtempe;
3114                       gTtsp = here->B3SOIFDgtempd;
3115                       gTtt  = here->B3SOIFDgtempT;
3116                       gTtdp = - (gTtg + gTtb + gTte + gTtsp);
3117 		  }
3118 
3119 		   if (model->B3SOIFDtype < 0)
3120 		   {
3121 		       ceqbodcon = -ceqbodcon;
3122 		       ceqbody = -ceqbody;
3123 		       ceqbs = -ceqbs;
3124 		       ceqbd = -ceqbd;
3125 		       ceqqg = -ceqqg;
3126 		       ceqqb = -ceqqb;
3127 		       ceqqd = -ceqqd;
3128 		       ceqqe = -ceqqe;
3129 		   }
3130 
3131 		   (*(ckt->CKTrhs + here->B3SOIFDgNode) -= m * ceqqg);
3132 		   (*(ckt->CKTrhs + here->B3SOIFDdNodePrime) += m * (ceqbd - cdreq - ceqqd));
3133 		   (*(ckt->CKTrhs + here->B3SOIFDsNodePrime) += m * (cdreq + ceqbs + ceqqg
3134 							  + ceqqb + ceqqd + ceqqe));
3135 		   (*(ckt->CKTrhs + here->B3SOIFDeNode) -= m * ceqqe);
3136 
3137                    if (here->B3SOIFDbodyMod == 1) {
3138 		       (*(ckt->CKTrhs + here->B3SOIFDpNode) += m * ceqbodcon);
3139                    }
3140 
3141 		   if (selfheat) {
3142 		       (*(ckt->CKTrhs + here->B3SOIFDtempNode) -= m * (ceqth + ceqqth));
3143                    }
3144 
3145 
3146 
3147  if ((here->B3SOIFDdebugMod > 1) || (here->B3SOIFDdebugMod == -1))
3148 		   {
3149 	              *(ckt->CKTrhs + here->B3SOIFDvbsNode) = here->B3SOIFDvbsdio;
3150 		      *(ckt->CKTrhs + here->B3SOIFDidsNode) = here->B3SOIFDids;
3151 		      *(ckt->CKTrhs + here->B3SOIFDicNode) = here->B3SOIFDic;
3152 		      *(ckt->CKTrhs + here->B3SOIFDibsNode) = here->B3SOIFDibs;
3153 		      *(ckt->CKTrhs + here->B3SOIFDibdNode) = here->B3SOIFDibd;
3154 		      *(ckt->CKTrhs + here->B3SOIFDiiiNode) = here->B3SOIFDiii;
3155 		      *(ckt->CKTrhs + here->B3SOIFDigidlNode) = here->B3SOIFDigidl;
3156 		      *(ckt->CKTrhs + here->B3SOIFDitunNode) = here->B3SOIFDitun;
3157 		      *(ckt->CKTrhs + here->B3SOIFDibpNode) = here->B3SOIFDibp;
3158 		      *(ckt->CKTrhs + here->B3SOIFDabeffNode) = here->B3SOIFDabeff;
3159 		      *(ckt->CKTrhs + here->B3SOIFDvbs0effNode) = here->B3SOIFDvbs0eff;
3160 		      *(ckt->CKTrhs + here->B3SOIFDvbseffNode) = here->B3SOIFDvbseff;
3161 		      *(ckt->CKTrhs + here->B3SOIFDxcNode) = here->B3SOIFDxc;
3162 		      *(ckt->CKTrhs + here->B3SOIFDcbbNode) = here->B3SOIFDcbb;
3163 		      *(ckt->CKTrhs + here->B3SOIFDcbdNode) = here->B3SOIFDcbd;
3164 		      *(ckt->CKTrhs + here->B3SOIFDcbgNode) = here->B3SOIFDcbg;
3165 		      *(ckt->CKTrhs + here->B3SOIFDqbfNode) = here->B3SOIFDqbf;
3166 		      *(ckt->CKTrhs + here->B3SOIFDqjsNode) = here->B3SOIFDqjs;
3167 		      *(ckt->CKTrhs + here->B3SOIFDqjdNode) = here->B3SOIFDqjd;
3168 
3169                       /* clean up last */
3170 		      *(ckt->CKTrhs + here->B3SOIFDgmNode) = Gm;
3171 		      *(ckt->CKTrhs + here->B3SOIFDgmbsNode) = Gmbs;
3172 		      *(ckt->CKTrhs + here->B3SOIFDgdsNode) = Gds;
3173 		      *(ckt->CKTrhs + here->B3SOIFDgmeNode) = Gme;
3174 		      *(ckt->CKTrhs + here->B3SOIFDqdNode) = qdrn;
3175 		      *(ckt->CKTrhs + here->B3SOIFDcbeNode) = Cbe;
3176 		      *(ckt->CKTrhs + here->B3SOIFDvbs0teffNode) = Vbs0teff;
3177 		      *(ckt->CKTrhs + here->B3SOIFDvthNode) = here->B3SOIFDvon;
3178 		      *(ckt->CKTrhs + here->B3SOIFDvgsteffNode) = Vgsteff;
3179 		      *(ckt->CKTrhs + here->B3SOIFDxcsatNode) = Xcsat;
3180 		      *(ckt->CKTrhs + here->B3SOIFDqaccNode) = -Qac0;
3181 		      *(ckt->CKTrhs + here->B3SOIFDqsub0Node) = Qsub0;
3182 		      *(ckt->CKTrhs + here->B3SOIFDqsubs1Node) = Qsubs1;
3183 		      *(ckt->CKTrhs + here->B3SOIFDqsubs2Node) = Qsubs2;
3184 		      *(ckt->CKTrhs + here->B3SOIFDvdscvNode) = VdsCV;
3185 		      *(ckt->CKTrhs + here->B3SOIFDvcscvNode) = VcsCV;
3186 		      *(ckt->CKTrhs + here->B3SOIFDqgNode) = qgate;
3187 		      *(ckt->CKTrhs + here->B3SOIFDqbNode) = qbody;
3188 		      *(ckt->CKTrhs + here->B3SOIFDqeNode) = qsub;
3189 		      *(ckt->CKTrhs + here->B3SOIFDdum1Node) = here->B3SOIFDdum1;
3190 		      *(ckt->CKTrhs + here->B3SOIFDdum2Node) = here->B3SOIFDdum2;
3191 		      *(ckt->CKTrhs + here->B3SOIFDdum3Node) = here->B3SOIFDdum3;
3192 		      *(ckt->CKTrhs + here->B3SOIFDdum4Node) = here->B3SOIFDdum4;
3193 		      *(ckt->CKTrhs + here->B3SOIFDdum5Node) = here->B3SOIFDdum5;
3194                       /* end clean up last */
3195 		   }
3196 
3197 
3198 		   /*
3199 		    *  load y matrix
3200 		    */
3201 		      (*(here->B3SOIFDEgPtr) += m * gcegb);
3202    		      (*(here->B3SOIFDEdpPtr) += m * gcedb);
3203    		      (*(here->B3SOIFDEspPtr) += m * gcesb);
3204 		      (*(here->B3SOIFDGePtr) += m * gcgeb);
3205 		      (*(here->B3SOIFDDPePtr) += m * (Gme + gddpe + gcdeb));
3206 		      (*(here->B3SOIFDSPePtr) += m * (gsspe - Gme + gcseb));
3207 
3208 		   (*(here->B3SOIFDEePtr) += m * gceeb);
3209 
3210 		   (*(here->B3SOIFDGgPtr) += m * (gcggb + ckt->CKTgmin));
3211 		   (*(here->B3SOIFDGdpPtr) += m * (gcgdb - ckt->CKTgmin));
3212 		   (*(here->B3SOIFDGspPtr) += m * gcgsb );
3213 
3214 		   (*(here->B3SOIFDDPgPtr) += m * ((Gm + gcdgb) + gddpg - ckt->CKTgmin));
3215 		   (*(here->B3SOIFDDPdpPtr) += m * ((here->B3SOIFDdrainConductance
3216 					 + here->B3SOIFDgds + gddpdp
3217 					 + RevSum + gcddb) + ckt->CKTgmin));
3218 		   (*(here->B3SOIFDDPspPtr) -= m * (-gddpsp + here->B3SOIFDgds + FwdSum - gcdsb));
3219 
3220 		   (*(here->B3SOIFDDPdPtr) -= m * here->B3SOIFDdrainConductance);
3221 
3222 		   (*(here->B3SOIFDSPgPtr) += m * (gcsgb - Gm + gsspg));
3223 		   (*(here->B3SOIFDSPdpPtr) -= m * (here->B3SOIFDgds - gsspdp + RevSum - gcsdb));
3224 		   (*(here->B3SOIFDSPspPtr) += m * (here->B3SOIFDsourceConductance
3225 					 + here->B3SOIFDgds + gsspsp
3226 					 + FwdSum + gcssb));
3227 		   (*(here->B3SOIFDSPsPtr) -= m * here->B3SOIFDsourceConductance);
3228 
3229 
3230 		   (*(here->B3SOIFDDdPtr) += m * here->B3SOIFDdrainConductance);
3231 		   (*(here->B3SOIFDDdpPtr) -= m * here->B3SOIFDdrainConductance);
3232 
3233 
3234 		   (*(here->B3SOIFDSsPtr) += m * here->B3SOIFDsourceConductance);
3235 		   (*(here->B3SOIFDSspPtr) -= m * here->B3SOIFDsourceConductance);
3236 
3237 		   if (here->B3SOIFDbodyMod == 1) {
3238 		      (*(here->B3SOIFDBpPtr) -= m * gppp);
3239 		      (*(here->B3SOIFDPbPtr) += m * gppb);
3240 		      (*(here->B3SOIFDPpPtr) += m * gppp);
3241 			(*(here->B3SOIFDPgPtr) += m * gppg);
3242 			(*(here->B3SOIFDPdpPtr) += m * gppdp);
3243 			(*(here->B3SOIFDPspPtr) += m * gppsp);
3244 			(*(here->B3SOIFDPePtr) += m * gppe);
3245 		   }
3246 
3247 		   if (selfheat)
3248                    {
3249 		      (*(here->B3SOIFDDPtempPtr) += m * (GmT + gddpT + gcdT));
3250 		      (*(here->B3SOIFDSPtempPtr) += m * (-GmT + gsspT + gcsT));
3251 		      (*(here->B3SOIFDBtempPtr) += m * (gbbT + gcbT));
3252                       (*(here->B3SOIFDEtempPtr) += m * gceT);
3253                       (*(here->B3SOIFDGtempPtr) += m * gcgT);
3254 		      if (here->B3SOIFDbodyMod == 1) {
3255 			  (*(here->B3SOIFDPtempPtr) += m * gppT);
3256 		      }
3257 		      (*(here->B3SOIFDTemptempPtr) += m * (gTtt  + 1/pParam->B3SOIFDrth + gcTt));
3258                       (*(here->B3SOIFDTempgPtr) += m * gTtg);
3259                       (*(here->B3SOIFDTempbPtr) += m * gTtb);
3260                       (*(here->B3SOIFDTempePtr) += m * gTte);
3261                       (*(here->B3SOIFDTempdpPtr) += m * gTtdp);
3262                       (*(here->B3SOIFDTempspPtr) += m * gTtsp);
3263 		   }
3264 
3265 		   if ((here->B3SOIFDdebugMod > 1) || (here->B3SOIFDdebugMod == -1))
3266 		   {
3267 		      *(here->B3SOIFDVbsPtr) += m * 1;
3268 		      *(here->B3SOIFDIdsPtr) += m * 1;
3269 		      *(here->B3SOIFDIcPtr) += m * 1;
3270 		      *(here->B3SOIFDIbsPtr) += m * 1;
3271 		      *(here->B3SOIFDIbdPtr) += m * 1;
3272 		      *(here->B3SOIFDIiiPtr) += m * 1;
3273 		      *(here->B3SOIFDIgidlPtr) += m * 1;
3274 		      *(here->B3SOIFDItunPtr) += m * 1;
3275 		      *(here->B3SOIFDIbpPtr) += m * 1;
3276 		      *(here->B3SOIFDAbeffPtr) += m * 1;
3277 		      *(here->B3SOIFDVbs0effPtr) += m * 1;
3278 		      *(here->B3SOIFDVbseffPtr) += m * 1;
3279 		      *(here->B3SOIFDXcPtr) += m * 1;
3280 		      *(here->B3SOIFDCbgPtr) += m * 1;
3281 		      *(here->B3SOIFDCbbPtr) += m * 1;
3282 		      *(here->B3SOIFDCbdPtr) += m * 1;
3283 		      *(here->B3SOIFDqbPtr) += m * 1;
3284 		      *(here->B3SOIFDQbfPtr) += m * 1;
3285 		      *(here->B3SOIFDQjsPtr) += m * 1;
3286 		      *(here->B3SOIFDQjdPtr) += m * 1;
3287 
3288                       /* clean up last */
3289 		      *(here->B3SOIFDGmPtr) += m * 1;
3290 		      *(here->B3SOIFDGmbsPtr) += m * 1;
3291 		      *(here->B3SOIFDGdsPtr) += m * 1;
3292 		      *(here->B3SOIFDGmePtr) += m * 1;
3293 		      *(here->B3SOIFDVbs0teffPtr) += m * 1;
3294 		      *(here->B3SOIFDVgsteffPtr) += m * 1;
3295 		      *(here->B3SOIFDCbePtr) += m * 1;
3296 		      *(here->B3SOIFDVthPtr) += m * 1;
3297 		      *(here->B3SOIFDXcsatPtr) += m * 1;
3298 		      *(here->B3SOIFDVdscvPtr) += m * 1;
3299 		      *(here->B3SOIFDVcscvPtr) += m * 1;
3300 		      *(here->B3SOIFDQaccPtr) += m * 1;
3301 		      *(here->B3SOIFDQsub0Ptr) += m * 1;
3302 		      *(here->B3SOIFDQsubs1Ptr) += m * 1;
3303 		      *(here->B3SOIFDQsubs2Ptr) += m * 1;
3304 		      *(here->B3SOIFDqgPtr) += m * 1;
3305 		      *(here->B3SOIFDqdPtr) += m * 1;
3306 		      *(here->B3SOIFDqePtr) += m * 1;
3307 		      *(here->B3SOIFDDum1Ptr) += m * 1;
3308 		      *(here->B3SOIFDDum2Ptr) += m * 1;
3309 		      *(here->B3SOIFDDum3Ptr) += m * 1;
3310 		      *(here->B3SOIFDDum4Ptr) += m * 1;
3311 		      *(here->B3SOIFDDum5Ptr) += m * 1;
3312                       /* end clean up last */
3313 		   }
3314 
3315 	line1000:  ;
3316 
3317 /*  Here NaN will be detected in any conductance or equivalent current.  Note
3318     that nandetect is initialized within the "if" statements  */
3319 
3320                    if ((nandetect = isnan (*(here->B3SOIFDGgPtr))) != 0)
3321                    { strcpy (nanmessage, "GgPtr"); }
3322                    else if ((nandetect = isnan (*(here->B3SOIFDGdpPtr))) != 0)
3323                    { strcpy (nanmessage, "GdpPtr"); }
3324                    else if ((nandetect = isnan (*(here->B3SOIFDGspPtr))) != 0)
3325                    { strcpy (nanmessage, "GspPtr"); }
3326                    else if ((nandetect = isnan (*(here->B3SOIFDDPgPtr))) != 0)
3327                    { strcpy (nanmessage, "DPgPtr"); }
3328                    else if ((nandetect = isnan (*(here->B3SOIFDDPdpPtr))) != 0)
3329                    { strcpy (nanmessage, "DPdpPtr"); }
3330                    else if ((nandetect = isnan (*(here->B3SOIFDDPspPtr))) != 0)
3331                    { strcpy (nanmessage, "DPspPtr"); }
3332                    else if ((nandetect = isnan (*(here->B3SOIFDSPgPtr))) != 0)
3333                    { strcpy (nanmessage, "SPgPtr"); }
3334                    else if ((nandetect = isnan (*(here->B3SOIFDSPdpPtr))) != 0)
3335                    { strcpy (nanmessage, "SPdpPtr"); }
3336                    else if ((nandetect = isnan (*(here->B3SOIFDSPspPtr))) != 0)
3337                    { strcpy (nanmessage, "SPspPtr"); }
3338                    else if ((nandetect = isnan (*(here->B3SOIFDEePtr))) != 0)
3339                    { strcpy (nanmessage, "EePtr"); }
3340 
3341                    /*  At this point, nandetect = 0 if none of the
3342                        conductances checked so far are NaN  */
3343 
3344                    if (nandetect == 0)
3345                    {
3346                      if ((nandetect = isnan (*(here->B3SOIFDEgPtr))) != 0)
3347                       { strcpy (nanmessage, "EgPtr"); }
3348                      else if ((nandetect = isnan (*(here->B3SOIFDEdpPtr))) != 0)
3349                       { strcpy (nanmessage, "EdpPtr"); }
3350                      else if ((nandetect = isnan (*(here->B3SOIFDEspPtr))) != 0)
3351                       { strcpy (nanmessage, "EspPtr"); }
3352                      else if ((nandetect = isnan (*(here->B3SOIFDGePtr))) != 0)
3353                       { strcpy (nanmessage, "GePtr"); }
3354                      else if ((nandetect = isnan (*(here->B3SOIFDDPePtr))) != 0)
3355                       { strcpy (nanmessage, "DPePtr"); }
3356                      else if ((nandetect = isnan (*(here->B3SOIFDSPePtr))) != 0)
3357                       { strcpy (nanmessage, "SPePtr"); } }
3358 
3359                    /*  Now check if self-heating caused NaN if nothing else
3360                        has so far (check tempnode current also)  */
3361 
3362                    if (selfheat && nandetect == 0)
3363                    {
3364                      if ((nandetect = isnan (*(here->B3SOIFDTemptempPtr))) != 0)
3365                       { strcpy (nanmessage, "TemptempPtr"); }
3366                      else if ((nandetect = isnan (*(here->B3SOIFDTempgPtr))) != 0)
3367                       { strcpy (nanmessage, "TempgPtr"); }
3368                      else if ((nandetect = isnan (*(here->B3SOIFDTempbPtr))) != 0)
3369                       { strcpy (nanmessage, "TempbPtr"); }
3370                      else if ((nandetect = isnan (*(here->B3SOIFDTempePtr))) != 0)
3371                       { strcpy (nanmessage, "TempePtr"); }
3372                      else if ((nandetect = isnan (*(here->B3SOIFDTempdpPtr))) != 0)
3373                       { strcpy (nanmessage, "TempdpPtr"); }
3374                      else if ((nandetect = isnan (*(here->B3SOIFDTempspPtr))) != 0)
3375                       { strcpy (nanmessage, "TempspPtr"); }
3376                      else if ((nandetect = isnan (*(here->B3SOIFDGtempPtr))) != 0)
3377                       { strcpy (nanmessage, "GtempPtr"); }
3378                      else if ((nandetect = isnan (*(here->B3SOIFDDPtempPtr))) != 0)
3379                       { strcpy (nanmessage, "DPtempPtr"); }
3380                      else if ((nandetect = isnan (*(here->B3SOIFDSPtempPtr))) != 0)
3381                       { strcpy (nanmessage, "SPtempPtr"); }
3382                      else if ((nandetect = isnan (*(here->B3SOIFDEtempPtr))) != 0)
3383                       { strcpy (nanmessage, "EtempPtr"); }
3384                      else if ((nandetect = isnan (*(here->B3SOIFDBtempPtr))) != 0)
3385                       { strcpy (nanmessage, "BtempPtr"); }
3386                      else if ((nandetect = isnan (*(ckt->CKTrhs + here->B3SOIFDtempNode))) != 0)
3387                       { strcpy (nanmessage, "tempNode"); }
3388                    }
3389 
3390                    /*  Lastly, check all equivalent currents (tempnode is
3391                        checked above  */
3392 
3393                    if (nandetect == 0)
3394                    {
3395                       if ((nandetect = isnan (*(ckt->CKTrhs
3396                                                 + here->B3SOIFDgNode))) != 0)
3397                       { strcpy (nanmessage, "gNode"); }
3398                       else if ((nandetect = isnan (*(ckt->CKTrhs
3399                                                      + here->B3SOIFDbNode))) != 0)
3400                       { strcpy (nanmessage, "bNode"); }
3401                       else if ((nandetect = isnan (*(ckt->CKTrhs
3402                                                      + here->B3SOIFDdNodePrime))) != 0)
3403                       { strcpy (nanmessage, "dpNode"); }
3404                       else if ((nandetect = isnan (*(ckt->CKTrhs
3405                                                   + here->B3SOIFDsNodePrime))) != 0)
3406                       { strcpy (nanmessage, "spNode"); }
3407                       else if ((nandetect = isnan (*(ckt->CKTrhs
3408                                                      + here->B3SOIFDeNode))) != 0)
3409                       { strcpy (nanmessage, "eNode"); }
3410                    }
3411 
3412                    /*  Now print error message if NaN detected.  Note that
3413                        error will only be printed once (the first time it is
3414                        encountered) each time SPICE is run since nanfound is
3415                        static variable  */
3416 
3417                    if (nanfound == 0 && nandetect)
3418                    {
3419                       fprintf(stderr, "Alberto says:  YOU TURKEY!  %s is NaN for instance %s at time %g!\n", nanmessage, here->B3SOIFDname, ckt->CKTtime);
3420                       nanfound = nandetect;
3421 		      fprintf(stderr, " The program exit!\n");
3422 		      controlled_exit(EXIT_FAILURE);
3423                    }
3424 
3425 		   if (here->B3SOIFDdebugMod > 2)
3426 		   {
3427       fprintf(fpdebug, "Ids = %.4e, Ic = %.4e, cqdrn = %.4e, gmin=%.3e\n",
3428                         Ids, Ic, cqdrn, ckt->CKTgmin);
3429 		      fprintf(fpdebug, "Iii = %.4e, Idgidl = %.4e, Ibs = %.14e\n",
3430 			      Iii, Idgidl, Ibs);
3431 		      fprintf(fpdebug, "Ibd = %.4e, Ibp = %.4e\n", Ibd, Ibp);
3432 		      fprintf(fpdebug, "qbody = %.5e, qbf = %.5e, qbe = %.5e\n",
3433 			      qbody, Qbf, -(Qe1+Qe2));
3434 		      fprintf(fpdebug, "qbs = %.5e, qbd = %.5e\n", qjs, qjd);
3435 		      fprintf(fpdebug, "qdrn = %.5e, qinv = %.5e\n", qdrn, qinv);
3436 
3437 
3438 
3439 
3440 /*  I am trying to debug the convergence problems here by printing out
3441     the entire Jacobian and equivalent current matrix  */
3442 
3443   if (here->B3SOIFDdebugMod > 4) {
3444   fprintf(fpdebug, "Ibtot = %.6e;\t Cbtot = %.6e;\n", Ibs+Ibp+Ibd-Iii-Idgidl-Isgidl, cqbody);
3445   fprintf(fpdebug, "ceqg = %.6e;\t ceqb = %.6e;\t ceqdp = %.6e;\t ceqsp = %.6e;\n",
3446                     *(ckt->CKTrhs + here->B3SOIFDgNode),
3447                     *(ckt->CKTrhs + here->B3SOIFDbNode),
3448                     *(ckt->CKTrhs + here->B3SOIFDdNodePrime),
3449                     *(ckt->CKTrhs + here->B3SOIFDsNodePrime));
3450   fprintf(fpdebug, "ceqe = %.6e;\t ceqp = %.6e;\t ceqth = %.6e;\n",
3451                     *(ckt->CKTrhs + here->B3SOIFDeNode),
3452                     *(ckt->CKTrhs + here->B3SOIFDpNode),
3453                     *(ckt->CKTrhs + here->B3SOIFDtempNode));
3454 
3455   fprintf(fpdebug, "Eg = %.5e;\t Edp = %.5e;\t Esp = %.5e;\t Eb = %.5e;\n",
3456 		   *(here->B3SOIFDEgPtr),
3457 		   *(here->B3SOIFDEdpPtr),
3458 		   *(here->B3SOIFDEspPtr),
3459 		   *(here->B3SOIFDEbPtr));
3460   fprintf(fpdebug, "Ee = %.5e;\t Gg = %.5e;\t Gdp = %.5e;\t Gsp = %.5e;\n",
3461 		   *(here->B3SOIFDEePtr),
3462 		   *(here->B3SOIFDGgPtr),
3463 		   *(here->B3SOIFDGdpPtr),
3464 		   *(here->B3SOIFDGspPtr));
3465   fprintf(fpdebug, "Gb = %.5e;\t Ge = %.5e;\t DPg = %.5e;\t DPdp = %.5e;\n",
3466 		   *(here->B3SOIFDGbPtr),
3467 		   *(here->B3SOIFDGePtr),
3468 		   *(here->B3SOIFDDPgPtr),
3469 		   *(here->B3SOIFDDPdpPtr));
3470   fprintf(fpdebug, "DPsp = %.5e;\t DPb = %.5e;\t DPe = %.5e;\t\n",
3471 		   *(here->B3SOIFDDPspPtr),
3472 		   *(here->B3SOIFDDPbPtr),
3473 		   *(here->B3SOIFDDPePtr));
3474   fprintf(fpdebug, "DPd = %.5e;\t SPg = %.5e;\t SPdp = %.5e;\t SPsp = %.5e;\n",
3475 		   *(here->B3SOIFDDPdPtr),
3476 		   *(here->B3SOIFDSPgPtr),
3477 		   *(here->B3SOIFDSPdpPtr),
3478 		   *(here->B3SOIFDSPspPtr));
3479   fprintf(fpdebug, "SPb = %.5e;\t SPe = %.5e;\t SPs = %.5e;\n",
3480 		   *(here->B3SOIFDSPbPtr),
3481 		   *(here->B3SOIFDSPePtr),
3482 		   *(here->B3SOIFDSPsPtr));
3483   fprintf(fpdebug, "Dd = %.5e;\t Ddp = %.5e;\t Ss = %.5e;\t Ssp = %.5e;\n",
3484 		   *(here->B3SOIFDDdPtr),
3485 		   *(here->B3SOIFDDdpPtr),
3486 		   *(here->B3SOIFDSsPtr),
3487 		   *(here->B3SOIFDSspPtr));
3488   fprintf(fpdebug, "Bg = %.5e;\t Bdp = %.5e;\t Bsp = %.5e;\t Bb = %.5e;\n",
3489 		   *(here->B3SOIFDBgPtr),
3490 		   *(here->B3SOIFDBdpPtr),
3491 		   *(here->B3SOIFDBspPtr),
3492 		   *(here->B3SOIFDBbPtr));
3493   fprintf(fpdebug, "Be = %.5e;\t Btot = %.5e;\t DPtot = %.5e;\n",
3494                    *(here->B3SOIFDBePtr),
3495                    *(here->B3SOIFDBgPtr) + *(here->B3SOIFDBdpPtr)
3496                    + *(here->B3SOIFDBspPtr) + *(here->B3SOIFDBbPtr)
3497                    + *(here->B3SOIFDBePtr),
3498 		   *(here->B3SOIFDDPePtr)
3499 		   + *(here->B3SOIFDDPgPtr) + *(here->B3SOIFDDPdpPtr)
3500 		   + *(here->B3SOIFDDPspPtr) + *(here->B3SOIFDDPbPtr));
3501 
3502   if (selfheat) {
3503     fprintf (fpdebug, "DPtemp = %.5e;\t SPtemp = %.5e;\t Btemp = %.5e;\n",
3504                       *(here->B3SOIFDDPtempPtr), *(here->B3SOIFDSPtempPtr),
3505                       *(here->B3SOIFDBtempPtr));
3506     fprintf (fpdebug, "Gtemp = %.5e;\t Etemp = %.5e;\n",
3507                       *(here->B3SOIFDGtempPtr), *(here->B3SOIFDEtempPtr));
3508     fprintf (fpdebug, "Tempg = %.5e;\t Tempdp = %.5e;\t Tempsp = %.5e;\t Tempb = %.5e;\n",
3509 		      *(here->B3SOIFDTempgPtr), *(here->B3SOIFDTempdpPtr),
3510 		      *(here->B3SOIFDTempspPtr), *(here->B3SOIFDTempbPtr));
3511     fprintf (fpdebug, "Tempe = %.5e;\t TempT = %.5e;\t Temptot = %.5e;\n",
3512 		      *(here->B3SOIFDTempePtr), *(here->B3SOIFDTemptempPtr),
3513 		      *(here->B3SOIFDTempgPtr) + *(here->B3SOIFDTempdpPtr)
3514 		      + *(here->B3SOIFDTempspPtr)+ *(here->B3SOIFDTempbPtr)
3515 		      + *(here->B3SOIFDTempePtr));
3516   }
3517 
3518    if (here->B3SOIFDbodyMod == 1)
3519    {
3520       fprintf(fpdebug, "ceqbodcon=%.5e;\t", ceqbodcon);
3521       fprintf(fpdebug, "Bp = %.5e;\t Pb = %.5e;\t Pp = %.5e;\n", -gppp, gppb, gppp);
3522       fprintf(fpdebug, "Pg=%.5e;\t Pdp=%.5e;\t Psp=%.5e;\t Pe=%.5e;\n",
3523                     gppg, gppdp, gppsp, gppe);
3524    }
3525 }
3526 
3527 	if (here->B3SOIFDdebugMod > 3)
3528         {
3529            fprintf(fpdebug, "Vth = %.4f, Vbs0eff = %.8f, Vdsat = %.4f\n",
3530                    Vth, Vbs0eff, Vdsat);
3531            fprintf(fpdebug, "ueff = %g, Vgsteff = %.4f, Vdseff = %.4f\n",
3532                    ueff, Vgsteff, Vdseff);
3533            fprintf(fpdebug, "Vthfd = %.4f, Vbs0mos = %.4f, Vbs0 = %.4f\n",
3534                    Vthfd, Vbs0mos, Vbs0);
3535            fprintf(fpdebug, "Vbs0t = %.4f, Vbsdio = %.8f\n",
3536                    Vbs0t, Vbsdio);
3537         }
3538 
3539               fclose(fpdebug);
3540            }
3541 
3542      here->B3SOIFDiterations++;  /*  increment the iteration counter  */
3543 
3544      }  /* End of Mosfet Instance */
3545 }   /* End of Model Instance */
3546 
3547 
3548 return(OK);
3549 }
3550 
3551