1 /*
2 Parker-Skellern MESFET model
3
4 Copyright (C) 1994, 1995, 1996 Macquarie University
5 All Rights Reserved
6 Author: Anthony Parker
7 Date: 2 Feb 1994 created
8 9 Feb 1994 correct NaN problem in strong cut-off region
9 20 MAR 1994 corrected capacitance initialization
10 24 MAR 1994 added parameter MVST
11 28 MAR 1994 reorganized declaration scopes
12 19 APR 1994 added new parameters: PS_HFETA, PS_HFE1, PS_HFE2,
13 PS_HFG1, and PS_HFG2
14 18 May 1994 corrected 1/0 error when PS_VSUB=0
15 15 Jul 1994 corrected errors in acload routine
16 10 Aug 1995 added PS_VSUB to gds += gm*PS_VSUB*mvst*(vgt-vgst*(a..
17 12 Sep 1995 changed _XXX to PS_XXX to aid portability
18 13 Sep 1995 change to give arg=1-1/subfac;
19 if(vst!=0) gds+=gm*PS_VSUB..;
20 gm *= arg;
21 10 Feb 1996 change to names to match MicroSim code.
22 5 Jul 1996 corrected diode eq (change Gmin*vgs to Gmin*vgd).
23
24 *****************************************************************************/
25
26 /*-----------
27 | functions defined in this file are:
28 PSids() returns dc drain source current and assigns other
29 current and branch conductances
30 qgg() static function that returns gate charge
31 PScharge() returns gate-source and gate-drain charge and capacitance
32 PSacload() returns small-signal conductance elements
33 PSinstanceinit() initializes model parameters
34 */
35
36 #define PSMODEL_C /* activate local definitions in psmesfet.h */
37 #include "psmodel.h"
38
39 /*-----------
40 | dc current and conductance calculation */
41 double
PSids(cref * ckt,modl * model,inst * here,double vgs,double vgd,double * igs,double * igd,double * ggs,double * ggd,double * Gm,double * Gds)42 PSids(
43 cref *ckt,
44 modl *model,
45 inst *here,
46 double vgs,
47 double vgd,
48 double *igs,
49 double *igd,
50 double *ggs,
51 double *ggd,
52 double *Gm,
53 double *Gds)
54 {
55 #define FX -10.0 /* not too small else fatal rounding error in (rpt-a_rpt) */
56 #define MX 40.0 /* maximum exponential argument */
57 #define EMX 2.353852668370199842e17 /* exp(MX) */
58
59 double idrain, arg;
60 double area = AREA;
61
62 { /* gate junction diodes */
63 double zz;
64 { /* gate-junction forward conduction */
65 double Gmin = GMIN;
66 double Vt = NVT;
67 double isat = IS * area;
68 if ((arg=vgs/Vt) > FX) {
69 if(arg < MX) {
70 *ggs=(zz=isat*exp(arg))/Vt+Gmin; *igs= zz -isat +Gmin*vgs;
71 } else {
72 *ggs=(zz=isat*EMX)/Vt+Gmin; *igs=zz*(arg-MX+1)-isat+Gmin*vgs;
73 }
74 } else {
75 *ggs = Gmin; *igs = -isat + Gmin * vgs;
76 }
77 if ((arg=vgd/Vt) > FX) {
78 if(arg < MX) {
79 *ggd=(zz=isat*exp(arg))/Vt+Gmin; *igd= zz -isat +Gmin*vgd;
80 } else {
81 *ggd=(zz=isat*EMX)/Vt+Gmin; *igd=zz*(arg-MX+1)-isat+Gmin*vgd;
82 }
83 } else {
84 *ggd = Gmin; *igd = -isat + Gmin * vgd;
85 }
86 }
87 { /* gate-junction reverse 'breakdown' conduction */
88 double Vbd = VBD;
89 double ibd = IBD * area;
90 if ((arg=-vgs/Vbd) > FX) {
91 if(arg < MX) {
92 *ggs += (zz=ibd*exp(arg))/Vbd; *igs -= zz-ibd;
93 } else {
94 *ggs += (zz=ibd*EMX)/Vbd; *igs -= zz*((arg-MX)+1) - ibd;
95 }
96 } else *igs += ibd;
97 if ((arg=-vgd/Vbd) > FX) {
98 if(arg < MX) {
99 *ggd += (zz=ibd*exp(arg))/Vbd; *igd -= zz-ibd;
100 } else {
101 *ggd += (zz=ibd*EMX)/Vbd; *igd -= zz*((arg-MX)+1) - ibd;
102 }
103 } else *igd += ibd;
104 }
105 }
106
107 { /* compute drain current and derivitives */
108 double gm, gds;
109 double vdst = vgs - vgd;
110 double stepofour = STEP * FOURTH;
111 { /* Include rate dependent threshold modulation */
112 double vgst, dvgd, dvgs, h, vgdtrap, vgstrap, eta, gam;
113 double vto = VTO;
114 double LFg = LFGAM, LFg1 = LFG1, LFg2 = LFG2;
115 double HFg = HFGAM, HFg1 = HFG1, HFg2 = HFG2;
116 double HFe = HFETA, HFe1 = HFE1, HFe2 = HFE2;
117 if(TRAN_ANAL) {
118 double taug = TAUG;
119 h = taug/(taug + stepofour); h*=h; h*=h; /*4th power*/
120 VGDTRAP_NOW = vgdtrap = h*VGDTRAP_BEFORE + (1-h) * vgd;
121 VGSTRAP_NOW = vgstrap = h*VGSTRAP_BEFORE + (1-h) * vgs;
122 } else {
123 h = 0;
124 VGDTRAP_NOW = vgdtrap = vgd;
125 VGSTRAP_NOW = vgstrap = vgs;
126 }
127 vgst = vgs - vto;
128 vgst -= ( LFg - LFg1*vgstrap + LFg2*vgdtrap)*vgdtrap;
129 vgst += (eta = HFe - HFe1*vgdtrap + HFe2*vgstrap)*(dvgs = vgstrap-vgs);
130 vgst += (gam = HFg - HFg1*vgstrap + HFg2*vgdtrap)*(dvgd = vgdtrap-vgd);
131 { /* Exponential Subthreshold effect ids(vgst,vdst) */
132 double vgt, subfac;
133 double mvst = MVST;
134 double vst = VSUB * (1 + mvst*vdst);
135 if (vgst > FX*vst) {
136 if (vgst > (arg=MX*vst)) { /* numerically large */
137 vgt = (EMX/(subfac = EMX+1))*(vgst-arg) + arg;
138 } else /* limit gate bias exponentially */
139 vgt = vst * log( subfac=(1 + exp(vgst/vst)) );
140 { /* Dual Power-law ids(vgt,vdst) */
141 double mQ = Q;
142 double PmQ = P - mQ;
143 double dvpd_dvdst=(double)D3*pow(vgt,PmQ);
144 double vdp = vdst*dvpd_dvdst; /*D3=P/Q/((VBI-vto)^PmQ)*/
145 { /* Early saturation effect ids(vgt,vdp) */
146 double za = (double)ZA; /* sqrt(1 + Z)/2 */
147 double mxi = MXI;
148 double vsatFac = vgt/(mxi*vgt + (double)XI_WOO);
149 double vsat=vgt/(1 + vsatFac);
150 double aa = za*vdp+vsat/2.0;
151 double a_aa = aa-vsat;
152 double rpt = sqrt( aa * aa + (arg=vsat*vsat*Z/4.0));
153 double a_rpt = sqrt(a_aa * a_aa + arg);
154 double vdt = (rpt - a_rpt);
155 double dvdt_dvdp = za * (aa/rpt - a_aa/a_rpt);
156 double dvdt_dvgt = (vdt - vdp*dvdt_dvdp)
157 *(1 + mxi*vsatFac*vsatFac)/(1 + vsatFac)/vgt;
158 { /* Intrinsic Q-law FET equation ids(vgt,vdt) */
159 gds=pow(vgt-vdt,mQ-1);
160 idrain = vdt*gds + vgt*(gm=pow(vgt,mQ-1)-gds);
161 gds *= mQ;
162 gm *= mQ;
163 }
164 gm += gds*dvdt_dvgt;
165 gds *= dvdt_dvdp;
166 }
167 gm += gds*PmQ*vdp/vgt;
168 gds *= dvpd_dvdst;
169 }
170 arg = 1 - 1/subfac;
171 if(vst != 0) gds += gm*VSUB*mvst*(vgt-vgst*arg)/vst;
172 gm *= arg;
173 } else { /* in extreme cut-off (numerically) */
174 idrain = gm = gds = 0.0e0;
175 }
176 }
177 gds += gm*(arg = h*gam +
178 (1-h)*(HFe1*dvgs-HFg2*dvgd+2*LFg2*vgdtrap-LFg1*vgstrap+LFg));
179 gm *= 1 - h*eta + (1-h)*(HFe2*dvgs -HFg1*dvgd + LFg1*vgdtrap) - arg;
180 }
181 { /* apply channel length modulation and beta scaling */
182 double lambda = LAM;
183 double beta = BETA * area;
184 gm *= (arg = beta*(1 + lambda*vdst));
185 gds = beta*lambda*idrain + gds*arg;
186 idrain *= arg;
187 }
188
189 { /* apply thermal reduction of drain current */
190 double h, pfac, pAverage;
191 double delta = DELT / area;
192 if(TRAN_ANAL) {
193 double taud = TAUD;
194 h = taud/(taud + stepofour); h*=h; h*=h;
195 POWR_NOW = pAverage = h*POWR_BEFORE + (1-h)*vdst*idrain;
196 } else {
197 POWR_NOW = POWR_BEFORE = pAverage = vdst*idrain; h = 0;
198 }
199 idrain /= (pfac = 1+pAverage*delta);
200 *Gm = gm * (arg = (h*delta*POWR_BEFORE + 1)/pfac/pfac);
201 *Gds = gds * arg - (1-h) * delta*idrain*idrain;
202 }
203 }
204 return(idrain);
205 }
206
207 /*-----------
208 | code based on Statz et. al. capacitance model, IEEE Tran ED Feb 87 */
209 static double
qgg(double vgs,double vgd,double gamma,double pb,double alpha,double vto,double vmax,double xc,double cgso,double cgdo,double * cgs,double * cgd)210 qgg(
211 double vgs,
212 double vgd,
213 double gamma,
214 double pb,
215 double alpha,
216 double vto,
217 double vmax,
218 double xc,
219 double cgso,
220 double cgdo,
221 double *cgs,
222 double *cgd
223 )
224 /*vgs, vgd, gamma, pb, alpha, vto, vmax, xc, cgso, cgdo, *cgs, *cgd;*/
225 {
226 double qrt, ext, Cgso, cpm, cplus, cminus;
227 double vds = vgs-vgd;
228 double d1_xc = 1-xc;
229 double vert = sqrt( vds * vds + alpha );
230 double veff = 0.5*(vgs + vgd + vert) + gamma*vds;
231 double vnr = d1_xc*(veff-vto);
232 double vnrt = sqrt( vnr*vnr + 0.04 );
233 double vnew = veff + 0.5*(vnrt - vnr);
234 if ( vnew < vmax ) {
235 ext = 0;
236 qrt = sqrt(1 - vnew/pb);
237 Cgso = 0.5*cgso/qrt*(1+xc + d1_xc*vnr/vnrt);
238 } else {
239 double vx = 0.5*(vnew-vmax);
240 double par = 1+vx/(pb-vmax);
241 qrt = sqrt(1 - vmax/pb);
242 ext = vx*(1 + par)/qrt;
243 Cgso = 0.5*cgso/qrt*(1+xc + d1_xc*vnr/vnrt) * par;
244 }
245 cplus = 0.5*(1 + (cpm = vds/vert)); cminus = cplus - cpm;
246 *cgs = Cgso*(cplus +gamma) + cgdo*(cminus+gamma);
247 *cgd = Cgso*(cminus-gamma) + cgdo*(cplus -gamma);
248 return(cgso*((pb+pb)*(1-qrt) + ext) + cgdo*(veff - vert));
249 }
250
251 /*-----------
252 | call during ac analysis initialisation or during transient analysis */
253 void
PScharge(cref * ckt,modl * model,inst * here,double vgs,double vgd,double * capgs,double * capgd)254 PScharge(
255 cref *ckt,
256 modl *model,
257 inst *here,
258 double vgs,
259 double vgd,
260 double *capgs,
261 double *capgd
262 )
263 {
264 #define QGG(a,b,c,d) qgg(a,b,gac,phib,alpha,vto,vmax,xc,czgs,czgd,c,d)
265 /* double qgg(); */
266
267 double czgs = CGS * AREA;
268 double czgd = CGD * AREA;
269 double vto = VTO;
270 double alpha = (double)ALPHA; /* (XI*woo/(XI+1)/2)^2 */
271 double xc = XC;
272 double vmax = VMAX;
273 double phib = VBI;
274 double gac = ACGAM;
275
276 if(/*TRAN_INIT ||*/ !TRAN_ANAL)
277 QGS_NOW = QGD_NOW = QGS_BEFORE = QGD_BEFORE
278 = QGG(vgs,vgd,capgs,capgd);
279 else {
280 double cgsna,cgsnc;
281 double cgdna,cgdnb, a_cap;
282 double vgs1 = VGS1;
283 double vgd1 = VGD1;
284 double qgga=QGG(vgs ,vgd ,&cgsna,&cgdna);
285 double qggb=QGG(vgs1,vgd ,&a_cap,&cgdnb);
286 double qggc=QGG(vgs ,vgd1,&cgsnc,&a_cap);
287 double qggd=QGG(vgs1,vgd1,&a_cap,&a_cap);
288 QGS_NOW = QGS_BEFORE + 0.5 * (qgga-qggb + qggc-qggd);
289 QGD_NOW = QGD_BEFORE + 0.5 * (qgga-qggc + qggb-qggd);
290 *capgs = 0.5 * (cgsna + cgsnc);
291 *capgd = 0.5 * (cgdna + cgdnb);
292 }
293 }
294
295
296 /*-----------
297 | call for each frequency in ac analysis */
298 void
PSacload(cref * ckt,modl * model,inst * here,double vgs,double vgd,double ids,double omega,double * Gm,double * xGm,double * Gds,double * xGds)299 PSacload(
300 cref *ckt,
301 modl *model,
302 inst *here,
303 double vgs,
304 double vgd,
305 double ids,
306 double omega,
307 double *Gm,
308 double *xGm,
309 double *Gds,
310 double *xGds
311 )
312 {
313 double arg;
314 double vds = vgs - vgd;
315 double LFgam = LFGAM;
316 double LFg1 = LFG1;
317 double LFg2 = LFG2*vgd;
318 double HFg1 = HFG1;
319 double HFg2 = HFG2*vgd;
320 double HFeta = HFETA;
321 double HFe1 = HFE1;
322 double HFe2 = HFE2*vgs;
323 double hfgam= HFGAM - HFg1*vgs + HFg2;
324 double eta = HFeta - HFe1*vgd + HFe2;
325 double lfga = LFgam - LFg1*vgs + LFg2 + LFg2;
326 double gmo = *Gm/(1 - lfga + LFg1*vgd);
327
328 double wtg = TAUG * omega;
329 double wtgdet = 1 + wtg*wtg;
330 double gwtgdet = gmo/wtgdet;
331
332 double gdsi = (arg=hfgam - lfga)*gwtgdet;
333 double gdsr = arg*gmo - gdsi;
334 double gmi = (eta + LFg1*vgd)*gwtgdet + gdsi;
335
336 double xgds = wtg*gdsi;
337 double gds = *Gds + gdsr;
338 double xgm = -wtg*gmi;
339 double gm = gmi + gmo*(1 -eta - hfgam);
340
341 double delta = DELT / AREA;
342 double wtd = TAUD * omega ;
343 double wtddet = 1 + wtd * wtd;
344 double fac = delta * ids;
345 double del = 1/(1 - fac * vds);
346 double dd = (del-1) / wtddet;
347 double dr = del - dd;
348 double di = wtd * dd;
349
350 double cdsqr = fac * ids * del * wtd/wtddet;
351
352 NG_IGNORE(ckt);
353
354 *Gm = dr*gm - di*xgm;
355 *xGm = di*gm + dr*xgm;
356
357 *Gds = dr*gds - di*xgds + cdsqr*wtd;
358 *xGds = di*gds + dr*xgds + cdsqr;
359 }
360
361
362 void /* call when temperature changes */
PSinstanceinit(modl * model,inst * here)363 PSinstanceinit(
364 modl *model,
365 inst *here
366 )
367 {
368 #ifndef PARAM_CAST /* allow casting to parameter type */
369 #define PARAM_CAST /* if not specified then don't cast */
370 #endif
371
372 double woo = (VBI - VTO);
373 XI_WOO = PARAM_CAST (XI * woo);
374 ZA = PARAM_CAST (sqrt(1 + Z)/2);
375 ALPHA = PARAM_CAST (XI_WOO*XI_WOO/(XI+1)/(XI+1)/ 4);
376 D3 = PARAM_CAST (P/Q/pow(woo,(P - Q)));
377 }
378