dune_core/
bsim2.rs

1//! BSIM2 — level 5 in the 3f5 (`dev/bsim2/b2ld.c`, `b2eval.c`,
2//! `b2temp.c`, `b2set.c`, `b2mpar.c`; the load shell is BSIM1's
3//! with the size-dependent parameters in a per-geometry struct).
4//! The drain current blends a subthreshold exponential and the
5//! strong-inversion law through a cubic spline on `Vgeff` between
6//! `VGLOW` and `VGHIGH`, with vertical-field mobility (`UA`/`UB`),
7//! velocity saturation (`U1`), the tanh Beta interpolation up to
8//! `VDD`, and the `AI`/`BI` impact-ionization factor. ngspice-47
9//! carries this code unchanged.
10
11use crate::consts::{ROOT2, VT0};
12use crate::diode::pnjlim;
13use crate::matrix::{AcCtx, AnKind, InitMode, LoadCtx, MatrixBuilder, Slot, TRASH, Var};
14use crate::mos1::{MosGeom, fetlim, limvds};
15
16/// `.model … nmos(level=5 …)` parameters, 3f5 names (`b2.c`).
17/// [`Bsim2Model::new`] carries the `b2set.c` defaults; junction
18/// potentials are clamped to ≥ 0.1 like BSIM1's.
19#[derive(Clone, Debug, Default)]
20pub struct Bsim2Model {
21    pub nmos: bool,
22    /// `VFB` and its L/W dependences.
23    pub vfb: f64,
24    pub lvfb: f64,
25    pub wvfb: f64,
26    /// `PHI` and its L/W dependences.
27    pub phi: f64,
28    pub lphi: f64,
29    pub wphi: f64,
30    /// `K1` and its L/W dependences.
31    pub k1: f64,
32    pub lk1: f64,
33    pub wk1: f64,
34    /// `K2` and its L/W dependences.
35    pub k2: f64,
36    pub lk2: f64,
37    pub wk2: f64,
38    /// `ETA0` and its L/W dependences.
39    pub eta0: f64,
40    pub leta0: f64,
41    pub weta0: f64,
42    /// `ETAB` and its L/W dependences.
43    pub etab: f64,
44    pub letab: f64,
45    pub wetab: f64,
46    /// `MU0B` and its L/W dependences.
47    pub mu0b: f64,
48    pub lmu0b: f64,
49    pub wmu0b: f64,
50    /// `MUS0` and its L/W dependences.
51    pub mus0: f64,
52    pub lmus0: f64,
53    pub wmus0: f64,
54    /// `MUSB` and its L/W dependences.
55    pub musb: f64,
56    pub lmusb: f64,
57    pub wmusb: f64,
58    /// `MU20` and its L/W dependences.
59    pub mu20: f64,
60    pub lmu20: f64,
61    pub wmu20: f64,
62    /// `MU2B` and its L/W dependences.
63    pub mu2b: f64,
64    pub lmu2b: f64,
65    pub wmu2b: f64,
66    /// `MU2G` and its L/W dependences.
67    pub mu2g: f64,
68    pub lmu2g: f64,
69    pub wmu2g: f64,
70    /// `MU30` and its L/W dependences.
71    pub mu30: f64,
72    pub lmu30: f64,
73    pub wmu30: f64,
74    /// `MU3B` and its L/W dependences.
75    pub mu3b: f64,
76    pub lmu3b: f64,
77    pub wmu3b: f64,
78    /// `MU3G` and its L/W dependences.
79    pub mu3g: f64,
80    pub lmu3g: f64,
81    pub wmu3g: f64,
82    /// `MU40` and its L/W dependences.
83    pub mu40: f64,
84    pub lmu40: f64,
85    pub wmu40: f64,
86    /// `MU4B` and its L/W dependences.
87    pub mu4b: f64,
88    pub lmu4b: f64,
89    pub wmu4b: f64,
90    /// `MU4G` and its L/W dependences.
91    pub mu4g: f64,
92    pub lmu4g: f64,
93    pub wmu4g: f64,
94    /// `UA0` and its L/W dependences.
95    pub ua0: f64,
96    pub lua0: f64,
97    pub wua0: f64,
98    /// `UAB` and its L/W dependences.
99    pub uab: f64,
100    pub luab: f64,
101    pub wuab: f64,
102    /// `UB0` and its L/W dependences.
103    pub ub0: f64,
104    pub lub0: f64,
105    pub wub0: f64,
106    /// `UBB` and its L/W dependences.
107    pub ubb: f64,
108    pub lubb: f64,
109    pub wubb: f64,
110    /// `U10` and its L/W dependences.
111    pub u10: f64,
112    pub lu10: f64,
113    pub wu10: f64,
114    /// `U1B` and its L/W dependences.
115    pub u1b: f64,
116    pub lu1b: f64,
117    pub wu1b: f64,
118    /// `U1D` and its L/W dependences.
119    pub u1d: f64,
120    pub lu1d: f64,
121    pub wu1d: f64,
122    /// `N0` and its L/W dependences.
123    pub n0: f64,
124    pub ln0: f64,
125    pub wn0: f64,
126    /// `NB` and its L/W dependences.
127    pub nb: f64,
128    pub lnb: f64,
129    pub wnb: f64,
130    /// `ND` and its L/W dependences.
131    pub nd: f64,
132    pub lnd: f64,
133    pub wnd: f64,
134    /// `VOF0` and its L/W dependences.
135    pub vof0: f64,
136    pub lvof0: f64,
137    pub wvof0: f64,
138    /// `VOFB` and its L/W dependences.
139    pub vofb: f64,
140    pub lvofb: f64,
141    pub wvofb: f64,
142    /// `VOFD` and its L/W dependences.
143    pub vofd: f64,
144    pub lvofd: f64,
145    pub wvofd: f64,
146    /// `AI0` and its L/W dependences.
147    pub ai0: f64,
148    pub lai0: f64,
149    pub wai0: f64,
150    /// `AIB` and its L/W dependences.
151    pub aib: f64,
152    pub laib: f64,
153    pub waib: f64,
154    /// `BI0` and its L/W dependences.
155    pub bi0: f64,
156    pub lbi0: f64,
157    pub wbi0: f64,
158    /// `BIB` and its L/W dependences.
159    pub bib: f64,
160    pub lbib: f64,
161    pub wbib: f64,
162    /// `VGHIGH` and its L/W dependences.
163    pub vghigh: f64,
164    pub lvghigh: f64,
165    pub wvghigh: f64,
166    /// `VGLOW` and its L/W dependences.
167    pub vglow: f64,
168    pub lvglow: f64,
169    pub wvglow: f64,
170    /// `MU0` (no L/W triple of its own).
171    pub mu0: f64,
172    /// `DL`/`DW` in µm; `TOX` in µm.
173    pub dl: f64,
174    pub dw: f64,
175    pub tox: f64,
176    /// `TEMP` in °C (BSIM2 uses it for its own Vtm).
177    pub temp: f64,
178    /// `VDD`/`VGG`/`VBB`: measurement ranges.
179    pub vdd: f64,
180    pub vgg: f64,
181    pub vbb: f64,
182    /// Overlap capacitances per meter.
183    pub cgso: f64,
184    pub cgdo: f64,
185    pub cgbo: f64,
186    /// `XPART` (`channelChargePartitionFlag`).
187    pub xpart: f64,
188    pub rsh: f64,
189    pub js: f64,
190    pub pb: f64,
191    pub mj: f64,
192    pub pbsw: f64,
193    pub mjsw: f64,
194    pub cj: f64,
195    pub cjsw: f64,
196    /// `WDF`/`DELL` (accepted, unused by the load).
197    pub wdf: f64,
198    pub dell: f64,
199    /// `KF`/`AF` (noise only; ngspice-47's `b2noi.c` — the 3f5's
200    /// BSIM2 has no noise and no such parameters). Defaults 0 and 1;
201    /// `new` sets `AF`, the derived `Default` leaves it 0.
202    pub kf: f64,
203    pub af: f64,
204}
205
206impl Bsim2Model {
207    /// The `b2set.c` defaults (everything not listed is 0).
208    pub fn new(nmos: bool) -> Self {
209        Bsim2Model {
210            nmos,
211            vfb: -1.0,
212            phi: 0.75,
213            k1: 0.8,
214            ua0: 0.2,
215            u10: 0.1,
216            mu0: 400.0,
217            mus0: 500.0,
218            mu20: 1.5,
219            mu30: 10.0,
220            n0: 1.4,
221            nb: 0.5,
222            vof0: 1.8,
223            vghigh: 0.2,
224            vglow: -0.15,
225            tox: 0.03,
226            temp: 27.0,
227            vdd: 5.0,
228            vgg: 5.0,
229            // Positive in the 3f5 too (`b2set.c`), though it is a
230            // lower Vbs clamp.
231            vbb: 5.0,
232            wdf: 10.0,
233            af: 1.0,
234            ..Default::default()
235        }
236    }
237
238    pub(crate) fn needs_prime(&self, squares: f64) -> bool {
239        self.rsh != 0.0 && squares != 0.0
240    }
241}
242
243/// Per-geometry parameters (`bsim2SizeDependParam`, `b2temp.c`):
244/// `p0 + pL·(1e-6/EffL) + pW·(1e-6/EffW)` plus the derived
245/// constants the load needs.
246#[derive(Debug)]
247pub(crate) struct Bsim2Size {
248    pub type_: f64,
249    pub vfb: f64,
250    pub phi: f64,
251    pub k1: f64,
252    pub k2: f64,
253    pub eta0: f64,
254    pub eta_b: f64,
255    pub beta0: f64,
256    pub beta0_b: f64,
257    pub betas0: f64,
258    pub betas_b: f64,
259    pub beta20: f64,
260    pub beta2_b: f64,
261    pub beta2_g: f64,
262    pub beta30: f64,
263    pub beta3_b: f64,
264    pub beta3_g: f64,
265    pub beta40: f64,
266    pub beta4_b: f64,
267    pub beta4_g: f64,
268    pub ua0: f64,
269    pub ua_b: f64,
270    pub ub0: f64,
271    pub ub_b: f64,
272    pub u10: f64,
273    pub u1_b: f64,
274    pub u1_d: f64,
275    pub n0: f64,
276    pub n_b: f64,
277    pub n_d: f64,
278    pub vof0: f64,
279    pub vof_b: f64,
280    pub vof_d: f64,
281    pub ai0: f64,
282    pub ai_b: f64,
283    pub bi0: f64,
284    pub bi_b: f64,
285    pub vghigh: f64,
286    pub vglow: f64,
287    pub cox_wl: f64,
288    pub one_third_cox_wl: f64,
289    pub two_third_cox_wl: f64,
290    pub gs_overlap: f64,
291    pub gd_overlap: f64,
292    pub gb_overlap: f64,
293    /// The overlap capacitances `B2load` recomputes per call: drawn
294    /// W for gate-source/drain and effective L for gate-bulk (the
295    /// `pParam` ones above, on effective W, feed only `B2acLoad`).
296    pub ld_gs_overlap: f64,
297    pub ld_gd_overlap: f64,
298    pub ld_gb_overlap: f64,
299    pub phis3: f64,
300    /// `pParam->Arg`: dBeta1/dVbs.
301    pub arg: f64,
302    pub vt0: f64,
303    pub drain_conductance: f64,
304    pub source_conductance: f64,
305    // Model-level copies.
306    pub vdd: f64,
307    pub vdd2: f64,
308    pub vgg2: f64,
309    pub vbb2: f64,
310    /// `B2Vtm = 8.625e-5·(temp°C + 273)`.
311    pub vtm: f64,
312    pub js: f64,
313    pub pb: f64,
314    pub pbsw: f64,
315    pub mj: f64,
316    pub mjsw: f64,
317    pub cj: f64,
318    pub cjsw: f64,
319    pub xpart: f64,
320    /// Noise: `KF`, `AF`, `B2Cox` (F/cm²) and `DW`/`DL` (µm).
321    pub kf: f64,
322    pub af: f64,
323    pub cox: f64,
324    pub delta_w: f64,
325    pub delta_l: f64,
326}
327
328impl Bsim2Size {
329    pub fn new(name: &str, m: &Bsim2Model, g: &MosGeom) -> Result<Self, String> {
330        let typ = if m.nmos { 1.0 } else { -1.0 };
331        let pb = f64::max(m.pb, 0.1);
332        let pbsw = f64::max(m.pbsw, 0.1);
333        if m.tox <= 0.0 {
334            return Err(format!("{name}: BSIM2 model needs TOX > 0"));
335        }
336        let cox = 3.453e-13 / (m.tox * 1.0e-4); // F/cm²
337
338        let eff_l = g.l - m.dl * 1e-6;
339        let eff_w = g.w - m.dw * 1e-6;
340        if eff_l <= 0.0 {
341            return Err(format!("{name}: effective channel length <= 0"));
342        }
343        if eff_w <= 0.0 {
344            return Err(format!("{name}: effective channel width <= 0"));
345        }
346        let inv_l = 1.0e-6 / eff_l;
347        let inv_w = 1.0e-6 / eff_w;
348        let p = |p0: f64, pl: f64, pw: f64| p0 + pl * inv_l + pw * inv_w;
349
350        let vfb = p(m.vfb, m.lvfb, m.wvfb);
351        let phi = p(m.phi, m.lphi, m.wphi);
352        let k1 = p(m.k1, m.lk1, m.wk1);
353        let k2 = p(m.k2, m.lk2, m.wk2);
354        let mut beta0 = m.mu0;
355        let beta0_b = p(m.mu0b, m.lmu0b, m.wmu0b);
356        let mut betas0 = p(m.mus0, m.lmus0, m.wmus0);
357        if betas0 < 1.01 * beta0 {
358            betas0 = 1.01 * beta0;
359        }
360        let mut betas_b = p(m.musb, m.lmusb, m.wmusb);
361        let tmp = betas0 - beta0 - beta0_b * m.vbb;
362        if -betas_b * m.vbb > tmp {
363            betas_b = -tmp / m.vbb;
364        }
365        let beta20 = p(m.mu20, m.lmu20, m.wmu20);
366        let beta2_b = p(m.mu2b, m.lmu2b, m.wmu2b);
367        let beta2_g = p(m.mu2g, m.lmu2g, m.wmu2g);
368        let mut beta30 = p(m.mu30, m.lmu30, m.wmu30);
369        let mut beta3_b = p(m.mu3b, m.lmu3b, m.wmu3b);
370        let beta3_g = p(m.mu3g, m.lmu3g, m.wmu3g);
371        let mut beta40 = p(m.mu40, m.lmu40, m.wmu40);
372        let mut beta4_b = p(m.mu4b, m.lmu4b, m.wmu4b);
373        let beta4_g = p(m.mu4g, m.lmu4g, m.wmu4g);
374
375        let cox_w_over_l = cox * eff_w / eff_l;
376        beta0 *= cox_w_over_l;
377        let beta0_b = beta0_b * cox_w_over_l;
378        betas0 *= cox_w_over_l;
379        betas_b *= cox_w_over_l;
380        beta30 *= cox_w_over_l;
381        beta3_b *= cox_w_over_l;
382        let beta3_g = beta3_g * cox_w_over_l;
383        beta40 *= cox_w_over_l;
384        beta4_b *= cox_w_over_l;
385        let beta4_g = beta4_g * cox_w_over_l;
386
387        let mut n0 = p(m.n0, m.ln0, m.wn0);
388        if n0 < 0.0 {
389            n0 = 0.0;
390        }
391
392        let sqrt_phi = phi.sqrt();
393        let phis3 = sqrt_phi * phi;
394        let cox_wl = cox * eff_l * eff_w * 1.0e4;
395        let series = |squares: f64| {
396            let r = m.rsh * squares;
397            if r != 0.0 { 1.0 / r } else { 0.0 }
398        };
399        let vt0 = vfb + phi + k1 * sqrt_phi - k2 * phi;
400
401        Ok(Bsim2Size {
402            type_: typ,
403            vfb,
404            phi,
405            k1,
406            k2,
407            eta0: p(m.eta0, m.leta0, m.weta0),
408            eta_b: p(m.etab, m.letab, m.wetab),
409            beta0,
410            beta0_b,
411            betas0,
412            betas_b,
413            beta20,
414            beta2_b,
415            beta2_g,
416            beta30,
417            beta3_b,
418            beta3_g,
419            beta40,
420            beta4_b,
421            beta4_g,
422            ua0: p(m.ua0, m.lua0, m.wua0),
423            ua_b: p(m.uab, m.luab, m.wuab),
424            ub0: p(m.ub0, m.lub0, m.wub0),
425            ub_b: p(m.ubb, m.lubb, m.wubb),
426            u10: p(m.u10, m.lu10, m.wu10),
427            u1_b: p(m.u1b, m.lu1b, m.wu1b),
428            u1_d: p(m.u1d, m.lu1d, m.wu1d),
429            n0,
430            n_b: p(m.nb, m.lnb, m.wnb),
431            n_d: p(m.nd, m.lnd, m.wnd),
432            vof0: p(m.vof0, m.lvof0, m.wvof0),
433            vof_b: p(m.vofb, m.lvofb, m.wvofb),
434            vof_d: p(m.vofd, m.lvofd, m.wvofd),
435            ai0: p(m.ai0, m.lai0, m.wai0),
436            ai_b: p(m.aib, m.laib, m.waib),
437            bi0: p(m.bi0, m.lbi0, m.wbi0),
438            bi_b: p(m.bib, m.lbib, m.wbib),
439            vghigh: p(m.vghigh, m.lvghigh, m.wvghigh),
440            vglow: p(m.vglow, m.lvglow, m.wvglow),
441            cox_wl,
442            one_third_cox_wl: cox_wl / 3.0,
443            two_third_cox_wl: 2.0 * cox_wl / 3.0,
444            gs_overlap: m.cgso * eff_w,
445            gd_overlap: m.cgdo * eff_w,
446            gb_overlap: m.cgbo * eff_l,
447            ld_gs_overlap: m.cgso * g.w,
448            ld_gd_overlap: m.cgdo * g.w,
449            ld_gb_overlap: m.cgbo * eff_l,
450            phis3,
451            arg: betas_b - beta0_b - m.vdd * (beta3_b - m.vdd * beta4_b),
452            vt0,
453            drain_conductance: series(g.nrd),
454            source_conductance: series(g.nrs),
455            vdd: m.vdd,
456            vdd2: 2.0 * m.vdd,
457            vgg2: 2.0 * m.vgg,
458            vbb2: 2.0 * m.vbb,
459            vtm: 8.625e-5 * (m.temp + 273.0),
460            js: m.js,
461            pb,
462            pbsw,
463            mj: m.mj,
464            mjsw: m.mjsw,
465            cj: m.cj,
466            cjsw: m.cjsw,
467            xpart: m.xpart,
468            kf: m.kf,
469            af: m.af,
470            cox,
471            delta_w: m.dw,
472            delta_l: m.dl,
473        })
474    }
475}
476
477// State offsets (`bsim2def.h` — identical to BSIM1's layout).
478const VBD: usize = 0;
479const VBS: usize = 1;
480const VGS: usize = 2;
481const VDS: usize = 3;
482const CD: usize = 4;
483const CBS: usize = 5;
484const CBD: usize = 6;
485const GM: usize = 7;
486const GDS: usize = 8;
487const GMBS: usize = 9;
488const GBD: usize = 10;
489const GBS: usize = 11;
490const QB: usize = 12;
491const IQB: usize = 13;
492const QG: usize = 14;
493const IQG: usize = 15;
494const QD: usize = 16;
495const IQD: usize = 17;
496const CGGB: usize = 18;
497const CGDB: usize = 19;
498const CGSB: usize = 20;
499const CBGB: usize = 21;
500const CBDB: usize = 22;
501const CBSB: usize = 23;
502const CAPBD: usize = 24;
503const CAPBS: usize = 26;
504const CDGB: usize = 28;
505const CDDB: usize = 29;
506const CDSB: usize = 30;
507const QBS: usize = 33;
508const QBD: usize = 34;
509pub(crate) const BSIM2_NUM_STATES: usize = 35;
510
511// Slot order (`b2set.c` — same as BSIM1/MOS1).
512const D_D: usize = 0;
513const G_G: usize = 1;
514const S_S: usize = 2;
515const B_B: usize = 3;
516const DP_DP: usize = 4;
517const SP_SP: usize = 5;
518const D_DP: usize = 6;
519const G_B: usize = 7;
520const G_DP: usize = 8;
521const G_SP: usize = 9;
522const S_SP: usize = 10;
523const B_DP: usize = 11;
524const B_SP: usize = 12;
525const DP_SP: usize = 13;
526const DP_D: usize = 14;
527const B_G: usize = 15;
528const DP_G: usize = 16;
529const SP_G: usize = 17;
530const SP_S: usize = 18;
531const DP_B: usize = 19;
532const SP_B: usize = 20;
533const SP_DP: usize = 21;
534
535/// Everything `B2evaluate` hands back.
536#[derive(Clone, Copy, Default)]
537struct Eval {
538    gm: f64,
539    gds: f64,
540    gmbs: f64,
541    qg: f64,
542    qb: f64,
543    qd: f64,
544    cggb: f64,
545    cgdb: f64,
546    cgsb: f64,
547    cbgb: f64,
548    cbdb: f64,
549    cbsb: f64,
550    cdgb: f64,
551    cddb: f64,
552    cdsb: f64,
553    cdrain: f64,
554    von: f64,
555    vdsat: f64,
556}
557
558pub(crate) struct Bsim2I {
559    drain: Var,
560    gate: Var,
561    source: Var,
562    bulk: Var,
563    d_p: Var,
564    s_p: Var,
565    t: Bsim2Size,
566    geom: MosGeom,
567    off: bool,
568    ic_vds: f64,
569    ic_vgs: f64,
570    ic_vbs: f64,
571    ic_given: crate::devices::IcGiven,
572    st: usize,
573    s: [Slot; 22],
574    mode: f64,
575    von: f64,
576    vdsat: f64,
577}
578
579impl Bsim2I {
580    #[allow(clippy::too_many_arguments)]
581    pub fn new(
582        nodes: [Var; 4],
583        primes: [Var; 2],
584        t: Bsim2Size,
585        geom: MosGeom,
586        off: bool,
587        ic: (f64, f64, f64),
588        st: usize,
589    ) -> Self {
590        let von = t.vt0;
591        Bsim2I {
592            drain: nodes[0],
593            gate: nodes[1],
594            source: nodes[2],
595            bulk: nodes[3],
596            d_p: primes[0],
597            s_p: primes[1],
598            t,
599            geom,
600            off,
601            ic_vds: ic.0,
602            ic_vgs: ic.1,
603            ic_vbs: ic.2,
604            ic_given: crate::devices::IcGiven::ALL,
605            st,
606            s: [TRASH; 22],
607            mode: 1.0,
608            von,
609            vdsat: 0.0,
610        }
611    }
612}
613
614impl Bsim2I {
615    /// `B2evaluate`: drain current, derivatives and terminal
616    /// charges at `vds ≥ 0` (the caller swaps terminals for the
617    /// inverse mode). The diagnostic `valuetypeflag` modes of the C
618    /// (negative XPART values turning Ids into a probe of gm, qg,
619    /// …) are not carried; every flag gets the plain
620    /// `Ids = max(Ids, 1e-50)` clamp.
621    #[allow(clippy::too_many_lines)]
622    fn evaluate(&self, vds_in: f64, vbs_in: f64, vgs_in: f64, charges: bool) -> Eval {
623        let p = &self.t;
624        let mut e = Eval::default();
625
626        // Range clamps (`b2eval.c` top).
627        let vbs = f64::max(vbs_in, p.vbb2);
628        let vgs = f64::min(vgs_in, p.vgg2);
629        let vds = f64::min(vds_in, p.vdd2);
630
631        // Threshold voltage.
632        let (phisb, d_phisb_dvb, t1s, d_t1s_dvb);
633        if vbs <= 0.0 {
634            phisb = p.phi - vbs;
635            d_phisb_dvb = -1.0;
636            t1s = phisb.sqrt();
637            d_t1s_dvb = -0.5 / t1s;
638        } else {
639            let tmp = p.phi / (p.phi + vbs);
640            phisb = p.phi * tmp;
641            d_phisb_dvb = -tmp * tmp;
642            t1s = p.phis3 / (p.phi + 0.5 * vbs);
643            d_t1s_dvb = -0.5 * t1s * t1s / p.phis3;
644        }
645        let eta = p.eta0 + p.eta_b * vbs;
646        let ua = p.ua0 + p.ua_b * vbs;
647        let ub = p.ub0 + p.ub_b * vbs;
648        let u1s = p.u10 + p.u1_b * vbs;
649        let vth = p.vfb + p.phi + p.k1 * t1s - p.k2 * phisb - eta * vds;
650        e.von = vth;
651        let d_vth_dvd = -eta;
652        let d_vth_dvb = p.k1 * d_t1s_dvb + p.k2 - p.eta_b * vds;
653        let vgst = vgs - vth;
654        let tmp = 1.0 / (1.744 + 0.8364 * phisb);
655        let gg = 1.0 - tmp;
656        let d_gg_dvb = 0.8364 * tmp * tmp * d_phisb_dvb;
657        let t0 = gg / t1s;
658        let tmp1 = 0.5 * t0 * p.k1;
659        let aa = 1.0 + tmp1;
660        let d_aa_dvb = (aa - 1.0) * (d_gg_dvb / gg - d_t1s_dvb / t1s);
661        let inv_aa = 1.0 / aa;
662        let vghigh = p.vghigh;
663        let vglow = p.vglow;
664
665        // Effective gate drive: direct, exponential subthreshold,
666        // or the cubic spline between VGLOW and VGHIGH.
667        let (vgeff, d_vgeff_dvg, d_vgeff_dvd, d_vgeff_dvb);
668        let mut exp0 = 0.0;
669        let mut exp1 = 0.0;
670        let mut n = 0.0;
671        if vgst >= vghigh || p.n0 == 0.0 {
672            vgeff = vgst;
673            d_vgeff_dvg = 1.0;
674            d_vgeff_dvd = -d_vth_dvd;
675            d_vgeff_dvb = -d_vth_dvb;
676        } else {
677            let vof = p.vof0 + p.vof_b * vbs + p.vof_d * vds;
678            n = p.n0 + p.n_b / t1s + p.n_d * vds;
679            let tmp = 0.5 / (n * p.vtm);
680            let exparg1 = f64::max(-vds / p.vtm, -30.0);
681            exp1 = exparg1.exp();
682            let tmp1 = f64::max(1.0 - exp1, 1.0e-18);
683            let tmp2 = 2.0 * aa * tmp1;
684            if vgst <= vglow {
685                let exparg = f64::max(vgst * tmp, -30.0);
686                exp0 = (0.5 * vof + exparg).exp();
687                vgeff = tmp2.sqrt() * p.vtm * exp0;
688                let t0 = n * p.vtm;
689                d_vgeff_dvg = vgeff * tmp;
690                d_vgeff_dvd =
691                    d_vgeff_dvg * (n / tmp1 * exp1 - d_vth_dvd - vgst * p.n_d / n + t0 * p.vof_d);
692                d_vgeff_dvb = d_vgeff_dvg
693                    * (p.vof_b * t0 - d_vth_dvb
694                        + p.n_b * vgst / (n * t1s * t1s) * d_t1s_dvb
695                        + t0 * inv_aa * d_aa_dvb);
696            } else {
697                let exparg = f64::max(vglow * tmp, -30.0);
698                exp0 = (0.5 * vof + exparg).exp();
699                let con3 = (2.0 * aa * (1.0 - exp1)).sqrt() * p.vtm * exp0;
700                let con1 = vghigh;
701                let con4 = con3 * tmp;
702                let sqr_high = vghigh * vghigh;
703                let sqr_low = vglow * vglow;
704                let cub_high = vghigh * sqr_high;
705                let cub_low = vglow * sqr_low;
706                let t0 = 2.0 * vghigh;
707                let t1 = 2.0 * vglow;
708                let t2 = 3.0 * sqr_high;
709                let t3 = 3.0 * sqr_low;
710                let t4 = vghigh - vglow;
711                let t5 = sqr_high - sqr_low;
712                let t6 = cub_high - cub_low;
713                let t7 = con1 - con3;
714                let delta = 1.0 / ((t1 - t0) * t6 + (t2 - t3) * t5 + (t0 * t3 - t1 * t2) * t4);
715                let coeffb =
716                    (t1 - con4 * t0) * t6 + (con4 * t2 - t3) * t5 + (t0 * t3 - t1 * t2) * t7;
717                let coeffc = (con4 - 1.0) * t6 + (t2 - t3) * t7 + (t3 - con4 * t2) * t4;
718                let coeffd = (t1 - t0) * t7 + (1.0 - con4) * t5 + (con4 * t0 - t1) * t4;
719                let coeffa = sqr_high * (coeffc + coeffd * t0);
720                vgeff = (coeffa + vgst * (coeffb + vgst * (coeffc + vgst * coeffd))) * delta;
721                d_vgeff_dvg = (coeffb + vgst * (2.0 * coeffc + 3.0 * vgst * coeffd)) * delta;
722                let t7b = con3 * tmp;
723                let t8 = d_t1s_dvb * p.n_b / (t1s * t1s * n);
724                let t9 = n * p.vtm;
725                let d_con3_dvd = t7b * (n * exp1 / tmp1 - vglow * p.n_d / n + t9 * p.vof_d);
726                let d_con3_dvb = t7b * (t9 * inv_aa * d_aa_dvb + vglow * t8 + t9 * p.vof_b);
727                let d_con4_dvd = tmp * d_con3_dvd - t7b * p.n_d / n;
728                let d_con4_dvb = tmp * d_con3_dvb + t7b * t8;
729                let d_coeffb_dvd =
730                    d_con4_dvd * (t2 * t5 - t0 * t6) + d_con3_dvd * (t1 * t2 - t0 * t3);
731                let d_coeffc_dvd = d_con4_dvd * (t6 - t2 * t4) + d_con3_dvd * (t3 - t2);
732                let d_coeffd_dvd = d_con4_dvd * (t0 * t4 - t5) + d_con3_dvd * (t0 - t1);
733                let d_coeffa_dvd = sqr_high * (d_coeffc_dvd + d_coeffd_dvd * t0);
734                d_vgeff_dvd = -d_vgeff_dvg * d_vth_dvd
735                    + (d_coeffa_dvd
736                        + vgst * (d_coeffb_dvd + vgst * (d_coeffc_dvd + vgst * d_coeffd_dvd)))
737                        * delta;
738                let d_coeffb_dvb =
739                    d_con4_dvb * (t2 * t5 - t0 * t6) + d_con3_dvb * (t1 * t2 - t0 * t3);
740                let d_coeffc_dvb = d_con4_dvb * (t6 - t2 * t4) + d_con3_dvb * (t3 - t2);
741                let d_coeffd_dvb = d_con4_dvb * (t0 * t4 - t5) + d_con3_dvb * (t0 - t1);
742                let d_coeffa_dvb = sqr_high * (d_coeffc_dvb + d_coeffd_dvb * t0);
743                d_vgeff_dvb = -d_vgeff_dvg * d_vth_dvb
744                    + (d_coeffa_dvb
745                        + vgst * (d_coeffb_dvb + vgst * (d_coeffc_dvb + vgst * d_coeffd_dvb)))
746                        * delta;
747            }
748        }
749
750        let mut ids = 0.0;
751        let mut gm = 0.0;
752        let mut gds = 0.0;
753        let mut gmb = 0.0;
754        let mut vdsat = 0.0;
755        // The charge block reuses these when `Vgst >= VGHIGH` (the C
756        // keeps the DC section's values in the same locals).
757        let (mut dc_dvdsat_dvd, mut dc_dvdsat_dvg, mut dc_dvdsat_dvb) = (0.0, 0.0, 0.0);
758
759        if vgeff > 0.0 {
760            let mut uvert = 1.0 + vgeff * (ua + vgeff * ub);
761            uvert = f64::max(uvert, 0.2);
762            let inv_uvert = 1.0 / uvert;
763            let t8 = ua + 2.0 * ub * vgeff;
764            let d_uvert_dvg = t8 * d_vgeff_dvg;
765            let d_uvert_dvd = t8 * d_vgeff_dvd;
766            let d_uvert_dvb = t8 * d_vgeff_dvb + vgeff * (p.ua_b + vgeff * p.ub_b);
767            let t8 = u1s * inv_aa * inv_uvert;
768            let vc = t8 * vgeff;
769            let t9 = vc * inv_uvert;
770            let d_vc_dvg = t8 * d_vgeff_dvg - t9 * d_uvert_dvg;
771            let d_vc_dvd = t8 * d_vgeff_dvd - t9 * d_uvert_dvd;
772            let d_vc_dvb = t8 * d_vgeff_dvb + p.u1_b * vgeff * inv_aa * inv_uvert
773                - vc * inv_aa * d_aa_dvb
774                - t9 * d_uvert_dvb;
775            let tmp2 = (1.0 + 2.0 * vc).sqrt();
776            let kk = 0.5 * (1.0 + vc + tmp2);
777            let inv_kk = 1.0 / kk;
778            let d_kk_dvc = 0.5 + 0.5 / tmp2;
779            let sqrt_kk = kk.sqrt();
780            let t8 = inv_aa / sqrt_kk;
781            vdsat = f64::max(vgeff * t8, 1.0e-18);
782            let inv_vdsat = 1.0 / vdsat;
783            let t9 = 0.5 * vdsat * inv_kk * d_kk_dvc;
784            let d_vdsat_dvd = t8 * d_vgeff_dvd - t9 * d_vc_dvd;
785            let d_vdsat_dvg = t8 * d_vgeff_dvg - t9 * d_vc_dvg;
786            let d_vdsat_dvb = t8 * d_vgeff_dvb - t9 * d_vc_dvb - vdsat * inv_aa * d_aa_dvb;
787            dc_dvdsat_dvd = d_vdsat_dvd;
788            dc_dvdsat_dvg = d_vdsat_dvg;
789            dc_dvdsat_dvb = d_vdsat_dvb;
790
791            let beta0 = p.beta0 + p.beta0_b * vbs;
792            let betas = p.betas0 + p.betas_b * vbs;
793            let beta2 = p.beta20 + p.beta2_b * vbs + p.beta2_g * vgs;
794            let beta3 = p.beta30 + p.beta3_b * vbs + p.beta3_g * vgs;
795            let beta4 = p.beta40 + p.beta4_b * vbs + p.beta4_g * vgs;
796            let beta1 = betas - (beta0 + p.vdd * (beta3 - p.vdd * beta4));
797            let t0 = f64::min(vds * beta2 * inv_vdsat, 30.0);
798            let t1 = t0.exp();
799            let t2 = t1 * t1;
800            let t3 = t2 + 1.0;
801            let tanh = (t2 - 1.0) / t3;
802            let sqrsech = 4.0 * t2 / (t3 * t3);
803            let beta = beta0 + beta1 * tanh + vds * (beta3 - beta4 * vds);
804            let t4 = beta1 * sqrsech * inv_vdsat;
805            let t5 = p.vdd * tanh;
806            let d_beta_dvd = beta3 - 2.0 * beta4 * vds + t4 * (beta2 - t0 * d_vdsat_dvd);
807            let d_beta_dvg = t4 * (p.beta2_g * vds - t0 * d_vdsat_dvg) + p.beta3_g * (vds - t5)
808                - p.beta4_g * (vds * vds - p.vdd * t5);
809            let d_beta1_dvb = p.arg;
810            let d_beta_dvb = p.beta0_b
811                + d_beta1_dvb * tanh
812                + vds * (p.beta3_b - vds * p.beta4_b)
813                + t4 * (p.beta2_b * vds - t0 * d_vdsat_dvb);
814
815            if vgst > vglow {
816                if vds <= vdsat {
817                    // Triode region.
818                    let t3 = vds * inv_vdsat;
819                    let t4 = t3 - 1.0;
820                    let t2 = 1.0 - p.u1_d * t4 * t4;
821                    let u1 = u1s * t2;
822                    let mut utot = uvert + u1 * vds;
823                    utot = f64::max(utot, 0.5);
824                    let inv_utot = 1.0 / utot;
825                    let t5 = 2.0 * u1s * p.u1_d * inv_vdsat * t4;
826                    let d_u1_dvd = t5 * (t3 * d_vdsat_dvd - 1.0);
827                    let d_u1_dvg = t5 * t3 * d_vdsat_dvg;
828                    let d_u1_dvb = t5 * t3 * d_vdsat_dvb + p.u1_b * t2;
829                    let d_utot_dvd = d_uvert_dvd + u1 + vds * d_u1_dvd;
830                    let d_utot_dvg = d_uvert_dvg + vds * d_u1_dvg;
831                    let d_utot_dvb = d_uvert_dvb + vds * d_u1_dvb;
832                    let tmp1 = vgeff - 0.5 * aa * vds;
833                    let tmp3 = tmp1 * vds;
834                    let betaeff = beta * inv_utot;
835                    ids = betaeff * tmp3;
836                    let t6 = ids / betaeff * inv_utot;
837                    gds = t6 * (d_beta_dvd - betaeff * d_utot_dvd)
838                        + betaeff * (tmp1 + (d_vgeff_dvd - 0.5 * aa) * vds);
839                    gm = t6 * (d_beta_dvg - betaeff * d_utot_dvg) + betaeff * vds * d_vgeff_dvg;
840                    gmb = t6 * (d_beta_dvb - betaeff * d_utot_dvb)
841                        + betaeff * vds * (d_vgeff_dvb - 0.5 * vds * d_aa_dvb);
842                } else {
843                    // Saturation.
844                    let tmp1 = vgeff * inv_aa * inv_kk;
845                    let tmp3 = 0.5 * vgeff * tmp1;
846                    let betaeff = beta * inv_uvert;
847                    ids = betaeff * tmp3;
848                    let t0 = ids / betaeff * inv_uvert;
849                    let t1 = betaeff * vgeff * inv_aa * inv_kk;
850                    let t2 = ids * inv_kk * d_kk_dvc;
851                    if p.ai0 != 0.0 {
852                        let ai = p.ai0 + p.ai_b * vbs;
853                        let bi = p.bi0 + p.bi_b * vbs;
854                        let t5 = f64::min(bi / (vds - vdsat), 30.0);
855                        let t6 = (-t5).exp();
856                        let fr = 1.0 + ai * t6;
857                        let t7 = t5 / (vds - vdsat);
858                        let t8 = (1.0 - fr) * t7;
859                        let d_fr_dvd = t8 * (d_vdsat_dvd - 1.0);
860                        let d_fr_dvg = t8 * d_vdsat_dvg;
861                        let d_fr_dvb =
862                            t8 * d_vdsat_dvb + t6 * (p.ai_b - ai * p.bi_b / (vds - vdsat));
863                        gds = (t0 * (d_beta_dvd - betaeff * d_uvert_dvd) + t1 * d_vgeff_dvd
864                            - t2 * d_vc_dvd)
865                            * fr
866                            + ids * d_fr_dvd;
867                        gm = (t0 * (d_beta_dvg - betaeff * d_uvert_dvg) + t1 * d_vgeff_dvg
868                            - t2 * d_vc_dvg)
869                            * fr
870                            + ids * d_fr_dvg;
871                        gmb = (t0 * (d_beta_dvb - betaeff * d_uvert_dvb) + t1 * d_vgeff_dvb
872                            - t2 * d_vc_dvb
873                            - ids * inv_aa * d_aa_dvb)
874                            * fr
875                            + ids * d_fr_dvb;
876                        ids *= fr;
877                    } else {
878                        gds = t0 * (d_beta_dvd - betaeff * d_uvert_dvd) + t1 * d_vgeff_dvd
879                            - t2 * d_vc_dvd;
880                        gm = t0 * (d_beta_dvg - betaeff * d_uvert_dvg) + t1 * d_vgeff_dvg
881                            - t2 * d_vc_dvg;
882                        gmb = t0 * (d_beta_dvb - betaeff * d_uvert_dvb) + t1 * d_vgeff_dvb
883                            - t2 * d_vc_dvb
884                            - ids * inv_aa * d_aa_dvb;
885                    }
886                }
887            } else {
888                // Subthreshold current (Vgst <= Vglow).
889                let t0 = exp0 * exp0;
890                let t1 = exp1;
891                ids = beta * p.vtm * p.vtm * t0 * (1.0 - t1);
892                let t2 = ids / beta;
893                let t4 = n * p.vtm;
894                let t3 = ids / t4;
895                let (fr, d_fr_dvd, d_fr_dvg, d_fr_dvb);
896                if vds > vdsat && p.ai0 != 0.0 {
897                    let ai = p.ai0 + p.ai_b * vbs;
898                    let bi = p.bi0 + p.bi_b * vbs;
899                    let t5 = f64::min(bi / (vds - vdsat), 30.0);
900                    let t6 = (-t5).exp();
901                    fr = 1.0 + ai * t6;
902                    let t7 = t5 / (vds - vdsat);
903                    let t8 = (1.0 - fr) * t7;
904                    d_fr_dvd = t8 * (d_vdsat_dvd - 1.0);
905                    d_fr_dvg = t8 * d_vdsat_dvg;
906                    d_fr_dvb = t8 * d_vdsat_dvb + t6 * (p.ai_b - ai * p.bi_b / (vds - vdsat));
907                } else {
908                    fr = 1.0;
909                    d_fr_dvd = 0.0;
910                    d_fr_dvg = 0.0;
911                    d_fr_dvb = 0.0;
912                }
913                gds = (t2 * d_beta_dvd
914                    + t3 * (p.vof_d * t4 - d_vth_dvd - p.n_d * vgst / n)
915                    + beta * p.vtm * t0 * t1)
916                    * fr
917                    + ids * d_fr_dvd;
918                gm = (t2 * d_beta_dvg + t3) * fr + ids * d_fr_dvg;
919                gmb = (t2 * d_beta_dvb
920                    + t3 * (p.vof_b * t4 - d_vth_dvb + p.n_b * vgst / (n * t1s * t1s) * d_t1s_dvb))
921                    * fr
922                    + ids * d_fr_dvb;
923                ids *= fr;
924            }
925        }
926
927        // `gds = MAX(gds, 1e-20)`.
928        gds = f64::max(gds, 1.0e-20);
929
930        // Charges: zeroed outside transient needs (XPART > 1 would
931        // also zero them, as in the C).
932        if charges && p.xpart <= 1.0 {
933            let (vbseff, d_vbseff_dvb);
934            if vbs < 0.0 {
935                vbseff = vbs;
936                d_vbseff_dvb = 1.0;
937            } else {
938                vbseff = p.phi - phisb;
939                d_vbseff_dvb = -d_phisb_dvb;
940            }
941            let arg1 = vgs - vbseff - p.vfb;
942            let arg2 = arg1 - vgst;
943            let qbulk = p.one_third_cox_wl * arg2;
944            let d_qbulk_dvb = p.one_third_cox_wl * (d_vth_dvb - d_vbseff_dvb);
945            let d_qbulk_dvd = p.one_third_cox_wl * d_vth_dvd;
946            if arg1 <= 0.0 {
947                // Accumulation.
948                e.qg = p.cox_wl * arg1;
949                e.qb = -e.qg;
950                e.cggb = p.cox_wl;
951                e.cgsb = -e.cggb * (1.0 - d_vbseff_dvb);
952                e.cbgb = -p.cox_wl;
953                e.cbsb = -e.cgsb;
954            } else if vgst <= 0.0 {
955                // Subthreshold.
956                let t2 = arg1 / arg2;
957                let t3 = t2 * t2 * (p.cox_wl - p.two_third_cox_wl * t2);
958                e.qg = p.cox_wl * arg1 * (1.0 - t2 * (1.0 - t2 / 3.0));
959                e.qb = -e.qg;
960                e.cggb = p.cox_wl * (1.0 - t2 * (2.0 - t2));
961                let tmp = t3 * d_vth_dvb - (e.cggb + t3) * d_vbseff_dvb;
962                e.cgdb = t3 * d_vth_dvd;
963                e.cgsb = -(e.cggb + e.cgdb + tmp);
964                e.cbgb = -e.cggb;
965                e.cbdb = -e.cgdb;
966                e.cbsb = -e.cgsb;
967            } else {
968                // Strong inversion: recompute Vdsat without Vgeff
969                // smoothing below VGHIGH (the C's inner block).
970                let mut vdsat_c = vdsat;
971                let mut d_vdsat_dvd = dc_dvdsat_dvd;
972                let mut d_vdsat_dvg = dc_dvdsat_dvg;
973                let mut d_vdsat_dvb = dc_dvdsat_dvb;
974                if vgst < p.vghigh {
975                    let mut uvert = 1.0 + vgst * (ua + vgst * ub);
976                    uvert = f64::max(uvert, 0.2);
977                    let inv_uvert = 1.0 / uvert;
978                    let d_uvert_dvg = ua + 2.0 * ub * vgst;
979                    let d_uvert_dvd = -d_uvert_dvg * d_vth_dvd;
980                    let d_uvert_dvb = -d_uvert_dvg * d_vth_dvb + vgst * (p.ua_b + vgst * p.ub_b);
981                    let t8 = u1s * inv_aa * inv_uvert;
982                    let vc = t8 * vgst;
983                    let t9 = vc * inv_uvert;
984                    let d_vc_dvg = t8 - t9 * d_uvert_dvg;
985                    let d_vc_dvd = -t8 * d_vth_dvd - t9 * d_uvert_dvd;
986                    let d_vc_dvb = -t8 * d_vth_dvb + p.u1_b * vgst * inv_aa * inv_uvert
987                        - vc * inv_aa * d_aa_dvb
988                        - t9 * d_uvert_dvb;
989                    let tmp2 = (1.0 + 2.0 * vc).sqrt();
990                    let kk = 0.5 * (1.0 + vc + tmp2);
991                    let inv_kk = 1.0 / kk;
992                    let d_kk_dvc = 0.5 + 0.5 / tmp2;
993                    let t8 = inv_aa / kk.sqrt();
994                    vdsat_c = vgst * t8;
995                    // `*vdsat` reports the charge block's Vdsat.
996                    vdsat = vdsat_c;
997                    let t9 = 0.5 * vdsat_c * inv_kk * d_kk_dvc;
998                    d_vdsat_dvd = -t8 * d_vth_dvd - t9 * d_vc_dvd;
999                    d_vdsat_dvg = t8 - t9 * d_vc_dvg;
1000                    d_vdsat_dvb = -t8 * d_vth_dvb - t9 * d_vc_dvb - vdsat_c * inv_aa * d_aa_dvb;
1001                }
1002                if vds >= vdsat_c {
1003                    // Saturation region.
1004                    e.cggb = p.two_third_cox_wl;
1005                    e.cgdb = -e.cggb * d_vth_dvd + d_qbulk_dvd;
1006                    let tmp = -e.cggb * d_vth_dvb + d_qbulk_dvb;
1007                    e.cgsb = -(e.cggb + e.cgdb + tmp);
1008                    e.cbgb = 0.0;
1009                    e.cbdb = -d_qbulk_dvd;
1010                    e.cbsb = d_qbulk_dvd + d_qbulk_dvb;
1011                    e.cdgb = -0.4 * e.cggb;
1012                    let tmp = -e.cdgb * d_vth_dvb;
1013                    e.cddb = -e.cdgb * d_vth_dvd;
1014                    e.cdsb = -(e.cdgb + e.cddb + tmp);
1015                    e.qb = -qbulk;
1016                    e.qg = p.two_third_cox_wl * vgst + qbulk;
1017                    e.qd = e.cdgb * vgst;
1018                } else {
1019                    // Linear region.
1020                    let t7 = vds / vdsat_c;
1021                    let t8 = vgst / vdsat_c;
1022                    let t6 = t7 * t8;
1023                    let t9 = 1.0 - t7;
1024                    let vgdt = vgst * t9;
1025                    let t0 = vgst / (vgst + vgdt);
1026                    let t1 = vgdt / (vgst + vgdt);
1027                    let t5 = t0 * t1;
1028                    let t2 = 1.0 - t1 + t5;
1029                    let t3 = 1.0 - t0 + t5;
1030                    let d_vgdt_dvg = t9 + t6 * d_vdsat_dvg;
1031                    let d_vgdt_dvd = t6 * d_vdsat_dvd - t8 - t9 * d_vth_dvd;
1032                    let d_vgdt_dvb = t6 * d_vdsat_dvb - t9 * d_vth_dvb;
1033                    e.qg = p.two_third_cox_wl * (vgst + vgdt - vgdt * t0) + qbulk;
1034                    e.qb = -qbulk;
1035                    e.qd = -p.one_third_cox_wl
1036                        * (0.2 * vgdt + 0.8 * vgst + vgdt * t1 + 0.2 * t5 * (vgdt - vgst));
1037                    e.cggb = p.two_third_cox_wl * (t2 + t3 * d_vgdt_dvg);
1038                    let tmp = d_qbulk_dvb + p.two_third_cox_wl * (t3 * d_vgdt_dvb - t2 * d_vth_dvb);
1039                    e.cgdb = p.two_third_cox_wl * (t3 * d_vgdt_dvd - t2 * d_vth_dvd) + d_qbulk_dvd;
1040                    e.cgsb = -(e.cggb + e.cgdb + tmp);
1041                    let t2b = 0.8 - 0.4 * t1 * (2.0 * t1 + t0 + t0 * (t1 - t0));
1042                    let t3b = 0.2 + t1 + t0 * (1.0 - 0.4 * t0 * (t1 + 3.0 * t0));
1043                    e.cdgb = -p.one_third_cox_wl * (t2b + t3b * d_vgdt_dvg);
1044                    let tmp = p.one_third_cox_wl * (t2b * d_vth_dvb - t3b * d_vgdt_dvb);
1045                    e.cddb = p.one_third_cox_wl * (t2b * d_vth_dvd - t3b * d_vgdt_dvd);
1046                    e.cdsb = -(e.cdgb + tmp + e.cddb);
1047                    e.cbgb = 0.0;
1048                    e.cbdb = -d_qbulk_dvd;
1049                    e.cbsb = d_qbulk_dvd + d_qbulk_dvb;
1050                }
1051            }
1052        }
1053
1054        e.gm = gm;
1055        e.gds = gds;
1056        e.gmbs = gmb;
1057        e.cdrain = f64::max(ids, 1.0e-50);
1058        e.vdsat = vdsat;
1059        e
1060    }
1061}
1062
1063impl crate::devices::Device for Bsim2I {
1064    fn par_load(&self) -> bool {
1065        true
1066    }
1067
1068    /// `B2getic` (`b2getic.c`): VBS, VDS, VGS not given
1069    /// on the card from the external nodes.
1070    fn set_ic(&mut self, rhs: &[f64]) {
1071        use crate::devices::getic;
1072        let g = self.ic_given;
1073        getic(g, 2, &mut self.ic_vbs, rhs, self.bulk, self.source);
1074        getic(g, 0, &mut self.ic_vds, rhs, self.drain, self.source);
1075        getic(g, 1, &mut self.ic_vgs, rhs, self.gate, self.source);
1076    }
1077
1078    fn set_ic_given(&mut self, given: crate::devices::IcGiven) {
1079        self.ic_given = given;
1080    }
1081
1082    fn alloc(&mut self, mb: &mut MatrixBuilder) {
1083        let (d, g, s, b) = (self.drain, self.gate, self.source, self.bulk);
1084        let (dp, sp) = (self.d_p, self.s_p);
1085        self.s = [
1086            mb.slot(d, d),
1087            mb.slot(g, g),
1088            mb.slot(s, s),
1089            mb.slot(b, b),
1090            mb.slot(dp, dp),
1091            mb.slot(sp, sp),
1092            mb.slot(d, dp),
1093            mb.slot(g, b),
1094            mb.slot(g, dp),
1095            mb.slot(g, sp),
1096            mb.slot(s, sp),
1097            mb.slot(b, dp),
1098            mb.slot(b, sp),
1099            mb.slot(dp, sp),
1100            mb.slot(dp, d),
1101            mb.slot(b, g),
1102            mb.slot(dp, g),
1103            mb.slot(sp, g),
1104            mb.slot(sp, s),
1105            mb.slot(dp, b),
1106            mb.slot(sp, b),
1107            mb.slot(sp, dp),
1108        ];
1109    }
1110
1111    /// `B2convTest` (`b2cvtest.c`).
1112    fn conv_test(&self, ctx: &crate::matrix::ConvCtx) -> bool {
1113        let typ = self.t.type_;
1114        let st = self.st;
1115        let s0 = |o: usize| ctx.states.get(0, st + o);
1116        let vbs = typ * (ctx.rhs[self.bulk] - ctx.rhs[self.s_p]);
1117        let vgs = typ * (ctx.rhs[self.gate] - ctx.rhs[self.s_p]);
1118        let vds = typ * (ctx.rhs[self.d_p] - ctx.rhs[self.s_p]);
1119        let vbd = vbs - vds;
1120        let vgd = vgs - vds;
1121        let vgdo = s0(VGS) - s0(VDS);
1122        let delvbs = vbs - s0(VBS);
1123        let delvbd = vbd - s0(VBD);
1124        let delvgs = vgs - s0(VGS);
1125        let delvds = vds - s0(VDS);
1126        let delvgd = vgd - vgdo;
1127        let cdhat = if self.mode >= 0.0 {
1128            s0(CD) - s0(GBD) * delvbd + s0(GMBS) * delvbs + s0(GM) * delvgs + s0(GDS) * delvds
1129        } else {
1130            s0(CD) - (s0(GBD) - s0(GMBS)) * delvbd - s0(GM) * delvgd + s0(GDS) * delvds
1131        };
1132        let cbhat = s0(CBS) + s0(CBD) + s0(GBD) * delvbd + s0(GBS) * delvbs;
1133        let cd = s0(CD);
1134        let cb = s0(CBS) + s0(CBD);
1135        let tol = ctx.reltol * f64::max(cdhat.abs(), cd.abs()) + ctx.abstol;
1136        if (cdhat - cd).abs() >= tol {
1137            return false;
1138        }
1139        let tol = ctx.reltol * f64::max(cbhat.abs(), cb.abs()) + ctx.abstol;
1140        (cbhat - cb).abs() <= tol
1141    }
1142
1143    #[allow(clippy::too_many_lines)]
1144    fn load(&mut self, ctx: &mut LoadCtx) {
1145        let typ = self.t.type_;
1146        let st = self.st;
1147
1148        let drain_sat_cur = f64::max(self.geom.ad * self.t.js, 1.0e-15);
1149        let source_sat_cur = f64::max(self.geom.as_ * self.t.js, 1.0e-15);
1150        let von_inst = typ * self.von;
1151        let vt0_inst = typ * self.t.vt0;
1152
1153        let mut check = true;
1154        let (mut vbs, mut vgs, mut vds);
1155        if ctx.mode == InitMode::Smsig {
1156            vbs = ctx.states.get(0, st + VBS);
1157            vgs = ctx.states.get(0, st + VGS);
1158            vds = ctx.states.get(0, st + VDS);
1159        } else if ctx.mode == InitMode::Tran {
1160            vbs = ctx.states.get(1, st + VBS);
1161            vgs = ctx.states.get(1, st + VGS);
1162            vds = ctx.states.get(1, st + VDS);
1163        } else if ctx.mode == InitMode::Jct && !self.off {
1164            vds = typ * self.ic_vds;
1165            vgs = typ * self.ic_vgs;
1166            vbs = typ * self.ic_vbs;
1167            if vds == 0.0 && vgs == 0.0 && vbs == 0.0 && !(ctx.kind == AnKind::TranOp && ctx.uic) {
1168                vbs = -1.0;
1169                vgs = vt0_inst;
1170                vds = 0.0;
1171            }
1172        } else if matches!(ctx.mode, InitMode::Jct | InitMode::Fix) && self.off {
1173            vbs = 0.0;
1174            vgs = 0.0;
1175            vds = 0.0;
1176        } else {
1177            if ctx.mode == InitMode::Pred {
1178                let xfact = ctx.pred_xfact;
1179                let mut pred = |o: usize| {
1180                    let s1 = ctx.states.get(1, st + o);
1181                    ctx.states.set(0, st + o, s1);
1182                    (1.0 + xfact) * s1 - xfact * ctx.states.get(2, st + o)
1183                };
1184                vbs = pred(VBS);
1185                vgs = pred(VGS);
1186                vds = pred(VDS);
1187                let vbd0 = ctx.states.get(0, st + VBS) - ctx.states.get(0, st + VDS);
1188                ctx.states.set(0, st + VBD, vbd0);
1189                for o in [
1190                    CD, CBS, CBD, GM, GDS, GMBS, GBD, GBS, CGGB, CBGB, CBSB, CGDB, CGSB, CBDB,
1191                    CDGB, CDDB, CDSB,
1192                ] {
1193                    let v = ctx.states.get(1, st + o);
1194                    ctx.states.set(0, st + o, v);
1195                }
1196            } else {
1197                vbs = typ * (ctx.sol[self.bulk] - ctx.sol[self.s_p]);
1198                vgs = typ * (ctx.sol[self.gate] - ctx.sol[self.s_p]);
1199                vds = typ * (ctx.sol[self.d_p] - ctx.sol[self.s_p]);
1200            }
1201            let mut vbd = vbs - vds;
1202            let mut vgd = vgs - vds;
1203            let vgdo = ctx.states.get(0, st + VGS) - ctx.states.get(0, st + VDS);
1204
1205            if ctx.states.get(0, st + VDS) >= 0.0 {
1206                vgs = fetlim(vgs, ctx.states.get(0, st + VGS), von_inst);
1207                vds = vgs - vgd;
1208                vds = limvds(vds, ctx.states.get(0, st + VDS));
1209            } else {
1210                vgd = fetlim(vgd, vgdo, von_inst);
1211                vds = vgs - vgd;
1212                vds = -limvds(-vds, -ctx.states.get(0, st + VDS));
1213                vgs = vgd + vds;
1214            }
1215            if vds >= 0.0 {
1216                let vcrit = VT0 * f64::ln(VT0 / (ROOT2 * source_sat_cur));
1217                vbs = pnjlim(vbs, ctx.states.get(0, st + VBS), VT0, vcrit, &mut check);
1218            } else {
1219                let vcrit = VT0 * f64::ln(VT0 / (ROOT2 * drain_sat_cur));
1220                vbd = pnjlim(vbd, ctx.states.get(0, st + VBD), VT0, vcrit, &mut check);
1221                vbs = vbd + vds;
1222            }
1223        }
1224
1225        let vbd = vbs - vds;
1226        let vgd = vgs - vds;
1227        let vgb = vgs - vbs;
1228
1229        // Bulk junctions on `CONSTvt0` (not the device temperature).
1230        let (gbs, cbs);
1231        if vbs <= 0.0 {
1232            gbs = source_sat_cur / VT0 + ctx.gmin;
1233            cbs = gbs * vbs;
1234        } else {
1235            let evbs = f64::exp(vbs / VT0);
1236            gbs = source_sat_cur * evbs / VT0 + ctx.gmin;
1237            cbs = source_sat_cur * (evbs - 1.0) + ctx.gmin * vbs;
1238        }
1239        let (gbd, cbd);
1240        if vbd <= 0.0 {
1241            gbd = drain_sat_cur / VT0 + ctx.gmin;
1242            cbd = gbd * vbd;
1243        } else {
1244            let evbd = f64::exp(vbd / VT0);
1245            gbd = drain_sat_cur * evbd / VT0 + ctx.gmin;
1246            cbd = drain_sat_cur * (evbd - 1.0) + ctx.gmin * vbd;
1247        }
1248
1249        self.mode = if vds >= 0.0 { 1.0 } else { -1.0 };
1250        let mode = self.mode;
1251
1252        let charges = ctx.kind == AnKind::Tran
1253            || ctx.mode == InitMode::Smsig
1254            || (ctx.kind == AnKind::TranOp && ctx.uic);
1255
1256        // `B2evaluate`, with the terminal swap for the inverse mode:
1257        // the drain-referred charge/capacitance rows come back on
1258        // the source side when `vds < 0`.
1259        let ev = if vds >= 0.0 {
1260            self.evaluate(vds, vbs, vgs, charges)
1261        } else {
1262            self.evaluate(-vds, vbd, vgd, charges)
1263        };
1264        let (cggb, cgdb, cgsb, cbgb, cbdb, cbsb);
1265        let (cdgb, cddb, cdsb, csgb, cssb, csdb);
1266        let (qgate, qbulk, qdrn, qsrc);
1267        if vds >= 0.0 {
1268            cggb = ev.cggb;
1269            cgdb = ev.cgdb;
1270            cgsb = ev.cgsb;
1271            cbgb = ev.cbgb;
1272            cbdb = ev.cbdb;
1273            cbsb = ev.cbsb;
1274            cdgb = ev.cdgb;
1275            cddb = ev.cddb;
1276            cdsb = ev.cdsb;
1277            csgb = 0.0;
1278            cssb = 0.0;
1279            csdb = 0.0;
1280            qgate = ev.qg;
1281            qbulk = ev.qb;
1282            qdrn = ev.qd;
1283            qsrc = 0.0;
1284        } else {
1285            cggb = ev.cggb;
1286            cgsb = ev.cgdb;
1287            cgdb = ev.cgsb;
1288            cbgb = ev.cbgb;
1289            cbsb = ev.cbdb;
1290            cbdb = ev.cbsb;
1291            csgb = ev.cdgb;
1292            cssb = ev.cddb;
1293            csdb = ev.cdsb;
1294            cdgb = 0.0;
1295            cddb = 0.0;
1296            cdsb = 0.0;
1297            qgate = ev.qg;
1298            qbulk = ev.qb;
1299            qsrc = ev.qd;
1300            qdrn = 0.0;
1301        }
1302        let (gm, gds, gmbs, cdrain) = (ev.gm, ev.gds, ev.gmbs, ev.cdrain);
1303        self.von = typ * ev.von;
1304        self.vdsat = typ * ev.vdsat;
1305
1306        let cd = mode * cdrain - cbd;
1307
1308        // Depletion capacitances.
1309        let (mut capbs, mut capbd) = (0.0, 0.0);
1310        if charges {
1311            let t = &self.t;
1312            let czbd = t.cj * self.geom.ad;
1313            let czbs = t.cj * self.geom.as_;
1314            let czbdsw = t.cjsw * self.geom.pd;
1315            let czbssw = t.cjsw * self.geom.ps;
1316            let (phi_b, phi_bsw, mj, mjsw) = (t.pb, t.pbsw, t.mj, t.mjsw);
1317            if vbs < 0.0 {
1318                let arg = 1.0 - vbs / phi_b;
1319                let argsw = 1.0 - vbs / phi_bsw;
1320                let sarg = f64::exp(-mj * arg.ln());
1321                let sargsw = f64::exp(-mjsw * argsw.ln());
1322                ctx.states.set(
1323                    0,
1324                    st + QBS,
1325                    phi_b * czbs * (1.0 - arg * sarg) / (1.0 - mj)
1326                        + phi_bsw * czbssw * (1.0 - argsw * sargsw) / (1.0 - mjsw),
1327                );
1328                capbs = czbs * sarg + czbssw * sargsw;
1329            } else {
1330                ctx.states.set(
1331                    0,
1332                    st + QBS,
1333                    vbs * (czbs + czbssw)
1334                        + vbs * vbs * (czbs * mj * 0.5 / phi_b + czbssw * mjsw * 0.5 / phi_bsw),
1335                );
1336                capbs = czbs + czbssw + vbs * (czbs * mj / phi_b + czbssw * mjsw / phi_bsw);
1337            }
1338            if vbd < 0.0 {
1339                let arg = 1.0 - vbd / phi_b;
1340                let argsw = 1.0 - vbd / phi_bsw;
1341                let sarg = f64::exp(-mj * arg.ln());
1342                let sargsw = f64::exp(-mjsw * argsw.ln());
1343                ctx.states.set(
1344                    0,
1345                    st + QBD,
1346                    phi_b * czbd * (1.0 - arg * sarg) / (1.0 - mj)
1347                        + phi_bsw * czbdsw * (1.0 - argsw * sargsw) / (1.0 - mjsw),
1348                );
1349                capbd = czbd * sarg + czbdsw * sargsw;
1350            } else {
1351                ctx.states.set(
1352                    0,
1353                    st + QBD,
1354                    vbd * (czbd + czbdsw)
1355                        + vbd * vbd * (czbd * mj * 0.5 / phi_b + czbdsw * mjsw * 0.5 / phi_bsw),
1356                );
1357                capbd = czbd + czbdsw + vbd * (czbd * mj / phi_b + czbdsw * mjsw / phi_bsw);
1358            }
1359        }
1360
1361        // Convergence.
1362        if (!self.off || ctx.mode != InitMode::Fix) && check {
1363            ctx.noncon = true;
1364        }
1365
1366        ctx.states.set(0, st + VBS, vbs);
1367        ctx.states.set(0, st + VBD, vbd);
1368        ctx.states.set(0, st + VGS, vgs);
1369        ctx.states.set(0, st + VDS, vds);
1370        ctx.states.set(0, st + CD, cd);
1371        ctx.states.set(0, st + CBS, cbs);
1372        ctx.states.set(0, st + CBD, cbd);
1373        ctx.states.set(0, st + GM, gm);
1374        ctx.states.set(0, st + GDS, gds);
1375        ctx.states.set(0, st + GMBS, gmbs);
1376        ctx.states.set(0, st + GBD, gbd);
1377        ctx.states.set(0, st + GBS, gbs);
1378        ctx.states.set(0, st + CGGB, cggb);
1379        ctx.states.set(0, st + CGDB, cgdb);
1380        ctx.states.set(0, st + CGSB, cgsb);
1381        ctx.states.set(0, st + CBGB, cbgb);
1382        ctx.states.set(0, st + CBDB, cbdb);
1383        ctx.states.set(0, st + CBSB, cbsb);
1384        ctx.states.set(0, st + CDGB, cdgb);
1385        ctx.states.set(0, st + CDDB, cddb);
1386        ctx.states.set(0, st + CDSB, cdsb);
1387        ctx.states.set(0, st + CAPBS, capbs);
1388        ctx.states.set(0, st + CAPBD, capbd);
1389
1390        // Gate charges and companion conductances (`B2mosCap`).
1391        let (mut ceqqg, mut ceqqb, mut ceqqd) = (0.0, 0.0, 0.0);
1392        let mut gc = [0.0f64; 12]; // gcggb gcgdb gcgsb gcbgb gcbdb gcbsb gcdgb gcddb gcdsb gcsgb gcsdb gcssb
1393        let t = &self.t;
1394        if charges {
1395            let ag0 = ctx.ag[0];
1396            let (args, c_b, c_d): ([f64; 8], [f64; 3], [f64; 3]);
1397            if mode > 0.0 {
1398                args = [
1399                    t.ld_gd_overlap,
1400                    t.ld_gs_overlap,
1401                    t.ld_gb_overlap,
1402                    capbd,
1403                    capbs,
1404                    cggb,
1405                    cgdb,
1406                    cgsb,
1407                ];
1408                c_b = [cbgb, cbdb, cbsb];
1409                c_d = [cdgb, cddb, cdsb];
1410            } else {
1411                args = [
1412                    t.ld_gs_overlap,
1413                    t.ld_gd_overlap,
1414                    t.ld_gb_overlap,
1415                    capbs,
1416                    capbd,
1417                    cggb,
1418                    cgsb,
1419                    cgdb,
1420                ];
1421                c_b = [cbgb, cbsb, cbdb];
1422                c_d = [csgb, cssb, csdb];
1423            }
1424            // `B2mosCap`: equivalent conductances…
1425            let gcdgb_l = (c_d[0] - args[0]) * ag0;
1426            let gcddb_l = (c_d[1] + args[3] + args[0]) * ag0;
1427            let gcdsb_l = c_d[2] * ag0;
1428            let gcsgb_l = -(args[5] + c_b[0] + c_d[0] + args[1]) * ag0;
1429            let gcsdb_l = -(args[6] + c_b[1] + c_d[1]) * ag0;
1430            let gcssb_l = (args[4] + args[1] - (args[7] + c_b[2] + c_d[2])) * ag0;
1431            let gcggb_l = (args[5] + args[0] + args[1] + args[2]) * ag0;
1432            let gcgdb_l = (args[6] - args[0]) * ag0;
1433            let gcgsb_l = (args[7] - args[1]) * ag0;
1434            let gcbgb_l = (c_b[0] - args[2]) * ag0;
1435            let gcbdb_l = (c_b[1] - args[3]) * ag0;
1436            let gcbsb_l = (c_b[2] - args[4]) * ag0;
1437            // …and total terminal charges.
1438            let (vg1, vg2) = if mode > 0.0 { (vgd, vgs) } else { (vgs, vgd) };
1439            let qgd = args[0] * vg1;
1440            let qgs = args[1] * vg2;
1441            let qgb = args[2] * vgb;
1442            let qgate = qgate + qgd + qgs + qgb;
1443            let qbulk = qbulk - qgb;
1444            let (mut qdrn, mut qsrc) = (qdrn, qsrc);
1445            if mode > 0.0 {
1446                qdrn -= qgd;
1447                qsrc = -(qgate + qbulk + qdrn);
1448            } else {
1449                qsrc -= qgs;
1450                qdrn = -(qgate + qbulk + qsrc);
1451            }
1452            let _ = qsrc;
1453
1454            if mode > 0.0 {
1455                gc = [
1456                    gcggb_l, gcgdb_l, gcgsb_l, gcbgb_l, gcbdb_l, gcbsb_l, gcdgb_l, gcddb_l,
1457                    gcdsb_l, gcsgb_l, gcsdb_l, gcssb_l,
1458                ];
1459            } else {
1460                // The reverse call swaps d and s in the outputs.
1461                gc = [
1462                    gcggb_l, gcgsb_l, gcgdb_l, gcbgb_l, gcbsb_l, gcbdb_l, gcsgb_l, gcssb_l,
1463                    gcsdb_l, gcdgb_l, gcdsb_l, gcddb_l,
1464                ];
1465                // and the gc arrays above are already in (g,b,d,s)
1466                // row order for the swapped call, matching the C's
1467                // pointer permutation.
1468            }
1469
1470            ctx.states.set(0, st + QG, qgate);
1471            ctx.states
1472                .set(0, st + QD, qdrn - ctx.states.get(0, st + QBD));
1473            ctx.states.set(
1474                0,
1475                st + QB,
1476                qbulk + ctx.states.get(0, st + QBD) + ctx.states.get(0, st + QBS),
1477            );
1478
1479            if ctx.mode == InitMode::Smsig {
1480                return;
1481            }
1482            if ctx.mode == InitMode::Tran {
1483                for o in [QB, QG, QD] {
1484                    let v = ctx.states.get(0, st + o);
1485                    ctx.states.set(1, st + o, v);
1486                }
1487            }
1488            ctx.integrate(0.0, st + QB);
1489            ctx.integrate(0.0, st + QG);
1490            ctx.integrate(0.0, st + QD);
1491
1492            // Equivalent charge currents.
1493            let cqgate = ctx.states.get(0, st + IQG);
1494            let cqbulk = ctx.states.get(0, st + IQB);
1495            let cqdrn = ctx.states.get(0, st + IQD);
1496            ceqqg = cqgate - gc[0] * vgb + gc[1] * vbd + gc[2] * vbs;
1497            ceqqb = cqbulk - gc[3] * vgb + gc[4] * vbd + gc[5] * vbs;
1498            ceqqd = cqdrn - gc[6] * vgb + gc[7] * vbd + gc[8] * vbs;
1499            if ctx.mode == InitMode::Tran {
1500                for o in [IQB, IQG, IQD] {
1501                    let v = ctx.states.get(0, st + o);
1502                    ctx.states.set(1, st + o, v);
1503                }
1504            }
1505        }
1506
1507        // Load the current vector.
1508        let ceqbs = typ * (cbs - (gbs - ctx.gmin) * vbs);
1509        let ceqbd = typ * (cbd - (gbd - ctx.gmin) * vbd);
1510        let ceqqg = typ * ceqqg;
1511        let ceqqb = typ * ceqqb;
1512        let ceqqd = typ * ceqqd;
1513        let (xnrm, xrev, cdreq);
1514        if mode >= 0.0 {
1515            xnrm = 1.0;
1516            xrev = 0.0;
1517            cdreq = typ * (cdrain - gds * vds - gm * vgs - gmbs * vbs);
1518        } else {
1519            xnrm = 0.0;
1520            xrev = 1.0;
1521            cdreq = -typ * (cdrain + gds * vds - gm * vgd - gmbs * vbd);
1522        }
1523        ctx.add_rhs(self.gate, -ceqqg);
1524        ctx.add_rhs(self.bulk, -(ceqbs + ceqbd + ceqqb));
1525        ctx.add_rhs(self.d_p, ceqbd - cdreq - ceqqd);
1526        ctx.add_rhs(self.s_p, cdreq + ceqbs + ceqqg + ceqqb + ceqqd);
1527
1528        // Load the Y matrix.
1529        let gdpr = t.drain_conductance;
1530        let gspr = t.source_conductance;
1531        ctx.add(self.s[D_D], gdpr);
1532        ctx.add(self.s[G_G], gc[0]);
1533        ctx.add(self.s[S_S], gspr);
1534        ctx.add(self.s[B_B], gbd + gbs - gc[3] - gc[4] - gc[5]);
1535        ctx.add(self.s[DP_DP], gdpr + gds + gbd + xrev * (gm + gmbs) + gc[7]);
1536        ctx.add(
1537            self.s[SP_SP],
1538            gspr + gds + gbs + xnrm * (gm + gmbs) + gc[11],
1539        );
1540        ctx.add(self.s[D_DP], -gdpr);
1541        ctx.add(self.s[G_B], -gc[0] - gc[1] - gc[2]);
1542        ctx.add(self.s[G_DP], gc[1]);
1543        ctx.add(self.s[G_SP], gc[2]);
1544        ctx.add(self.s[S_SP], -gspr);
1545        ctx.add(self.s[B_G], gc[3]);
1546        ctx.add(self.s[B_DP], -gbd + gc[4]);
1547        ctx.add(self.s[B_SP], -gbs + gc[5]);
1548        ctx.add(self.s[DP_D], -gdpr);
1549        ctx.add(self.s[DP_G], (xnrm - xrev) * gm + gc[6]);
1550        ctx.add(
1551            self.s[DP_B],
1552            -gbd + (xnrm - xrev) * gmbs - gc[6] - gc[7] - gc[8],
1553        );
1554        ctx.add(self.s[DP_SP], -gds - xnrm * (gm + gmbs) + gc[8]);
1555        ctx.add(self.s[SP_G], -(xnrm - xrev) * gm + gc[9]);
1556        ctx.add(self.s[SP_S], -gspr);
1557        ctx.add(
1558            self.s[SP_B],
1559            -gbs - (xnrm - xrev) * gmbs - gc[9] - gc[10] - gc[11],
1560        );
1561        ctx.add(self.s[SP_DP], -gds - xrev * (gm + gmbs) + gc[10]);
1562    }
1563
1564    /// `B2noise` (ngspice-47 `b2noi.c`): see [`crate::bsim1::b12_noise`].
1565    fn noise_suffixes(&self) -> Option<&'static [&'static str]> {
1566        Some(crate::bsim1::B12_NOISE_SUFFIXES)
1567    }
1568
1569    fn noise_dens(&self, ctx: &crate::matrix::NoiseCtx) -> Vec<f64> {
1570        let t = &self.t;
1571        let n = crate::bsim1::B12Noise {
1572            drain: self.drain,
1573            source: self.source,
1574            d_p: self.d_p,
1575            s_p: self.s_p,
1576            gdpr: t.drain_conductance,
1577            gspr: t.source_conductance,
1578            gm: ctx.states.get(0, self.st + GM),
1579            cd: ctx.states.get(0, self.st + CD),
1580            kf: t.kf,
1581            af: t.af,
1582            flicker_area: (self.geom.w - t.delta_w * 1e-6) * (self.geom.l - t.delta_l * 1e-6),
1583            cox: t.cox,
1584        };
1585        crate::bsim1::b12_noise(ctx, &n)
1586    }
1587
1588    fn trunc(&self, t: &crate::integ::TerrCtx, time_step: &mut f64) {
1589        crate::integ::terr(self.st + QB, t, time_step);
1590        crate::integ::terr(self.st + QG, t, time_step);
1591        crate::integ::terr(self.st + QD, t, time_step);
1592    }
1593
1594    /// `B2acLoad`.
1595    fn load_ac(&mut self, ctx: &mut AcCtx) {
1596        let st = self.st;
1597        let s0 = |o: usize| ctx.states.get(0, st + o);
1598        let (xnrm, xrev) = if self.mode < 0.0 {
1599            (0.0, 1.0)
1600        } else {
1601            (1.0, 0.0)
1602        };
1603        let gdpr = self.t.drain_conductance;
1604        let gspr = self.t.source_conductance;
1605        let gm = s0(GM);
1606        let gds = s0(GDS);
1607        let gmbs = s0(GMBS);
1608        let gbd = s0(GBD);
1609        let gbs = s0(GBS);
1610        let capbd = s0(CAPBD);
1611        let capbs = s0(CAPBS);
1612        let cggb = s0(CGGB);
1613        let cgsb = s0(CGSB);
1614        let cgdb = s0(CGDB);
1615        let cbgb = s0(CBGB);
1616        let cbsb = s0(CBSB);
1617        let cbdb = s0(CBDB);
1618        let cdgb = s0(CDGB);
1619        let cdsb = s0(CDSB);
1620        let cddb = s0(CDDB);
1621        let o = ctx.omega;
1622        let xcdgb = (cdgb - self.t.gd_overlap) * o;
1623        let xcddb = (cddb + capbd + self.t.gd_overlap) * o;
1624        let xcdsb = cdsb * o;
1625        let xcsgb = -(cggb + cbgb + cdgb + self.t.gs_overlap) * o;
1626        let xcsdb = -(cgdb + cbdb + cddb) * o;
1627        let xcssb = (capbs + self.t.gs_overlap - (cgsb + cbsb + cdsb)) * o;
1628        let xcggb = (cggb + self.t.gd_overlap + self.t.gs_overlap + self.t.gb_overlap) * o;
1629        let xcgdb = (cgdb - self.t.gd_overlap) * o;
1630        let xcgsb = (cgsb - self.t.gs_overlap) * o;
1631        let xcbgb = (cbgb - self.t.gb_overlap) * o;
1632        let xcbdb = (cbdb - capbd) * o;
1633        let xcbsb = (cbsb - capbs) * o;
1634
1635        ctx.add_b(self.s[G_G], xcggb);
1636        ctx.add_b(self.s[B_B], -xcbgb - xcbdb - xcbsb);
1637        ctx.add_b(self.s[DP_DP], xcddb);
1638        ctx.add_b(self.s[SP_SP], xcssb);
1639        ctx.add_b(self.s[G_B], -xcggb - xcgdb - xcgsb);
1640        ctx.add_b(self.s[G_DP], xcgdb);
1641        ctx.add_b(self.s[G_SP], xcgsb);
1642        ctx.add_b(self.s[B_G], xcbgb);
1643        ctx.add_b(self.s[B_DP], xcbdb);
1644        ctx.add_b(self.s[B_SP], xcbsb);
1645        ctx.add_b(self.s[DP_G], xcdgb);
1646        ctx.add_b(self.s[DP_B], -xcdgb - xcddb - xcdsb);
1647        ctx.add_b(self.s[DP_SP], xcdsb);
1648        ctx.add_b(self.s[SP_G], xcsgb);
1649        ctx.add_b(self.s[SP_B], -xcsgb - xcsdb - xcssb);
1650        ctx.add_b(self.s[SP_DP], xcsdb);
1651        ctx.add_g(self.s[D_D], gdpr);
1652        ctx.add_g(self.s[S_S], gspr);
1653        ctx.add_g(self.s[B_B], gbd + gbs);
1654        ctx.add_g(self.s[DP_DP], gdpr + gds + gbd + xrev * (gm + gmbs));
1655        ctx.add_g(self.s[SP_SP], gspr + gds + gbs + xnrm * (gm + gmbs));
1656        ctx.add_g(self.s[D_DP], -gdpr);
1657        ctx.add_g(self.s[S_SP], -gspr);
1658        ctx.add_g(self.s[B_DP], -gbd);
1659        ctx.add_g(self.s[B_SP], -gbs);
1660        ctx.add_g(self.s[DP_D], -gdpr);
1661        ctx.add_g(self.s[DP_G], (xnrm - xrev) * gm);
1662        ctx.add_g(self.s[DP_B], -gbd + (xnrm - xrev) * gmbs);
1663        ctx.add_g(self.s[DP_SP], -gds - xnrm * (gm + gmbs));
1664        ctx.add_g(self.s[SP_G], -(xnrm - xrev) * gm);
1665        ctx.add_g(self.s[SP_S], -gspr);
1666        ctx.add_g(self.s[SP_B], -gbs - (xnrm - xrev) * gmbs);
1667        ctx.add_g(self.s[SP_DP], -gds - xrev * (gm + gmbs));
1668    }
1669
1670    /// `B2pzLoad` (`b2pzld.c`, `b1pzld.c` with the overlaps read from
1671    /// `pParam`): see [`crate::bsim1::bsim_pz_load`].
1672    fn load_pz(&mut self, ctx: &mut crate::matrix::PzCtx) -> Result<(), &'static str> {
1673        let ov = crate::bsim1::BsimOverlaps {
1674            gd: self.t.gd_overlap,
1675            gs: self.t.gs_overlap,
1676            gb: self.t.gb_overlap,
1677        };
1678        let (gdpr, gspr) = (self.t.drain_conductance, self.t.source_conductance);
1679        crate::bsim1::bsim_pz_load(ctx, &self.s, self.st, self.mode, gdpr, gspr, &ov);
1680        Ok(())
1681    }
1682}
dune_core/bsim2.rs

dune_core/
bsim2.rs
