-- File: [Cherry]<Thyme>System>03>spModels.mesa
-- Last editted:
-- SChen April 29, 1983  12:03 PM
-- Wilhelm April 12, 1982  9:28 AM, reformated by Barth and stored under
--   [Cherry]<Barth>Thyme>1.97> .
-- Details at end of file.
DIRECTORY AltoDefs, Real, RealFns, spModelDefs;
spModels:  PROGRAM
  IMPORTS RealFns, spModelDefs =
  BEGIN

    modelError:  SIGNAL[s:  STRING] = CODE;

-- Temporary square root kluge until RealFns is fixed.

    SingleReal:  TYPE = RECORD[m2:    CARDINAL, -- Backwards!!!!
                               sign:  BOOLEAN,
                               exp:   [0..256),
                               m1:    [0..128)];

    SquareReal:  TYPE = RECORD[m2:     CARDINAL, -- Backwards!!!!
                               sign:   BOOLEAN,
                               expD2:  [0..128),
                               index:  [0..16),
                               m1:     [0..16)];

    guesses:  ARRAY[0..16) OF REAL;

    SqRt:  PROCEDURE[x: REAL] RETURNS[y:  REAL] =
      BEGIN
        xFmt:  SquareReal;
        yFmt:  SingleReal;

        IF x < 0.0 THEN SIGNAL modelError["SqRt of non-positive"];
        xFmt ← LOOPHOLE[x, SquareReal];
        yFmt ← LOOPHOLE[guesses[xFmt.index], SingleReal];
        yFmt.exp ← yFmt.exp + xFmt.expD2 - 63;
        y ← LOOPHOLE[yFmt, REAL];
        y ← LOOPHOLE[LOOPHOLE[y + x/y, LONG CARDINAL] - 40000000B, REAL];
        y ← LOOPHOLE[LOOPHOLE[y + x/y, LONG CARDINAL] - 40000000B, REAL];
      END;

    initSqRt:  PROCEDURE =
      BEGIN
        i:  CARDINAL;
        x1, x2:  REAL;
        xFmt:  SquareReal;

        FOR i IN [0..16) DO
          xFmt ← [0, FALSE, 63, i, 0];
          x1 ← LOOPHOLE[xFmt, REAL];
          IF i < 15 THEN xFmt ← [0, FALSE, 63, i + 1, 0]
          ELSE xFmt ← [0, FALSE, 64, 0, 0];
          x2 ← LOOPHOLE[xFmt, REAL];
          guesses[i] ← 2.0*(RealFns.SqRt[x1]*x2 -
                            RealFns.SqRt[x2]*x1)/(x2-x1)
        ENDLOOP
      END;


    NFET:  spModelDefs.model =
      BEGIN
        VgX:  CARDINAL = 0;
        VsX:  CARDINAL = 1;
        VdX:  CARDINAL = 2;
        VbX:  CARDINAL = 3;

        VtX:     CARDINAL = 0;
        hKp:     CARDINAL = 1;
        Phi2:    CARDINAL = 2;
        Phi2fb:  CARDINAL = 3;
        Kb:      CARDINAL = 4;
        ftKb:    CARDINAL = 5;
        hKbsq:   CARDINAL = 6;
        qKbsq:   CARDINAL = 7;
        Vfb:     CARDINAL = 8;
        Cox:     CARDINAL = 9;
        Cgsd0:   CARDINAL = 10;
        Cgsdp:   CARDINAL = 11;
        Ioff:    CARDINAL = 12;
        Io:      CARDINAL = 13;
        kTq:     CARDINAL = 14;
        expMax:  CARDINAL = 15;

        IdX:   CARDINAL = 0;
        CgbX:  CARDINAL = 1;
        CgsX:  CARDINAL = 2;
        CgdX:  CARDINAL = 3;

        Vg, Vs, Vd, Vb, Id:  REAL;
        Cgs, Cgd:  REAL;
        Vt, Vgst, Vgdt, VdSAT, Vsbphi2, Vdx, SqrtV:  REAL;
        t1, t2:  REAL;
        reverse:  BOOLEAN ← FALSE;

        Vg ← args[VgX];
        Vb ← args[VbX];
        Vd ← args[VdX];
        Vs ← args[VsX];

        t1 ← Vb - Vs;
        IF t1 > 0.0 THEN
          BEGIN
            t1 ← t1/parms[kTq];
            IF t1 >= parms[expMax] THEN
              SIGNAL spModelDefs.Retreat["NFET -- Source."];
            results[4] ← parms[Io]*(RealFns.Exp[t1] - 1.0)
          END
        ELSE results[4] ← -parms[Io];

        t2 ← Vb - Vd;
        IF t2 > 0.0 THEN
          BEGIN
            t2 ← t2/parms[kTq];
            IF t2 >= parms[expMax] THEN
              SIGNAL spModelDefs.Retreat["NFET -- Drain."];
            results[5] ← parms[Io]*(RealFns.Exp[t2] - 1.0)
          END
        ELSE results[5] ← -parms[Io];

        IF Vs > Vd THEN
          BEGIN
            reverse ← TRUE;
            t1 ← Vs;
            Vs ← Vd;
            Vd ← t1
          END;

        Vsbphi2 ← MAX[Vs - Vb + parms[Phi2], 0.0];
        SqrtV ← SqRt[MAX[0.0, parms[qKbsq] + Vg - parms[Vfb] - Vb]];
        Vt ← IF Vs = Vb THEN parms[VtX]
             ELSE parms[Phi2fb] + parms[Kb]*SqRt[Vsbphi2];
        Vgst ← Vg - Vs - Vt;
        IF Vgst > 0 THEN
          BEGIN
            VdSAT ← Vg - parms[Phi2fb] + parms[hKbsq] - parms[Kb]*SqrtV;
            IF VdSAT < Vs THEN SIGNAL modelError["Negative VdSAT"];
            Vdx ← IF Vd > VdSAT THEN VdSAT ELSE Vd;
            t1 ← Vsbphi2*SqRt[Vsbphi2];
            t2 ← Vdx - Vb + parms[Phi2];
            t2 ← t2*SqRt[t2];
            Id ← parms[hKp]*((2.0*(Vg - parms[Phi2fb]) - Vs - Vdx)*
                             (Vdx - Vs) + parms[ftKb]*(t1 - t2))
          END
        ELSE Id ← parms[Ioff];

        IF Id < 0.0 THEN
          IF Id < -parms[Ioff]THEN
            SIGNAL modelError["Negative drain current"]
          ELSE Id ← parms[Ioff];

        t1 ← Vg - Vb - parms[Vfb];
        IF t1 <= 0.0 THEN results[CgbX] ← parms[Cox]
        ELSE results[CgbX] ← 0.5*parms[Kb]*parms[Cox]/SqrtV;

        Vgdt ← Vg - Vd - Vt;
        t1 ← Vgst + Vgdt;
        t1 ← t1*t1;
        IF Vgst <= 0.0 THEN Cgs ← parms[Cgsdp]
        ELSE
          IF Vgdt < 0.0 THEN Cgs ← parms[Cgsd0] + parms[Cgsdp]
          ELSE Cgs ← parms[Cgsdp] + parms[Cgsd0]*(1.0 - Vgdt*Vgdt/t1);
        IF Vgdt <= 0.0 THEN Cgd ← parms[Cgsdp]
        ELSE Cgd ← parms[Cgsdp] + parms[Cgsd0]*(1.0 - Vgst*Vgst/t1);

        IF reverse THEN
          BEGIN
            results[IdX]  ← -Id;
            results[CgsX] ← Cgd;
            results[CgdX] ← Cgs
          END
        ELSE
          BEGIN
            results[IdX]  ← Id;
            results[CgsX] ← Cgs;
            results[CgdX] ← Cgd
          END
      END;

    PFET:  spModelDefs.model =
      BEGIN
        VgX:  CARDINAL = 0;
        VsX:  CARDINAL = 1;
        VdX:  CARDINAL = 2;
        VbX:  CARDINAL = 3;

        VtX:     CARDINAL = 0;
        hKp:     CARDINAL = 1;
        Phi2:    CARDINAL = 2;
        Phi2fb:  CARDINAL = 3;
        Kb
:      CARDINAL = 4;
        ftKb:    CARDINAL = 5;
        hKbsq:   CARDINAL = 6;
        qKbsq:   CARDINAL = 7;
        Vfb:     CARDINAL = 8;
        Cox:     CARDINAL = 9;
        Cgsd0:   CARDINAL = 10;
        Cgsdp:   CARDINAL = 11;
        Ioff:    CARDINAL = 12;
        Io:      CARDINAL = 13;
        kTq:     CARDINAL = 14;
        expMax:  CARDINAL = 15;

        IdX:   CARDINAL = 0;
        CgbX:  CARDINAL = 1;
        CgsX:  CARDINAL = 2;
        CgdX:  CARDINAL = 3;

        Vg, Vs, Vd, Vb, Id:  REAL;
        Cgs, Cgd:  REAL;
        Vt, Vgst, Vgdt, VdSAT, Vsbphi2, Vdx, SqrtV:  REAL;
        t1, t2:  REAL;
        reverse:  BOOLEAN ← FALSE;

        Vg ← args[VgX];
        Vb ← args[VbX];
        Vd ← args[VdX];
        Vs ← args[VsX];

        t1 ← Vs - Vb;
        IF t1 > 0.0 THEN
          BEGIN
            t1 ← t1/parms[kTq];
            IF t1 >= parms[expMax] THEN
              SIGNAL spModelDefs.Retreat["PFET -- Source."];
            results[4] ← parms[Io]*(RealFns.Exp[t1] - 1.0)
          END
        ELSE results[4] ← -parms[Io];

        t2 ← Vd - Vb;
        IF t2 > 0.0 THEN
          BEGIN
            t2 ← t2/parms[kTq];
            IF t2 >= parms[expMax] THEN
              SIGNAL spModelDefs.Retreat["PFET -- Drain."];
            results[5] ← parms[Io]*(RealFns.Exp[t2] - 1.0)
          END
        ELSE results[5] ← -parms[Io];

        IF Vs < Vd THEN
          BEGIN
            reverse ← TRUE;
            t1 ← Vs;
            Vs ← Vd;
            Vd ← t1
          END;

        Vsbphi2 ← MIN[Vs - Vb + parms[Phi2], 0.0];
        SqrtV ← SqRt[-MIN[0.0, -parms[qKbsq] + Vg - parms[Vfb] - Vb]];
        Vt ← IF Vs = Vb THEN parms[VtX]
             ELSE parms[Phi2fb] - parms[Kb]*SqRt[-Vsbphi2];
        Vgst ← Vg - Vs - Vt;
        IF Vgst < 0 THEN
          BEGIN
            VdSAT ← Vg - parms[Phi2fb] - parms[hKbsq] + parms[Kb]*SqrtV;
            IF VdSAT > Vs THEN SIGNAL modelError["Positive VdSAT"];
            Vdx ← IF Vd < VdSAT THEN VdSAT ELSE Vd;
            t1 ← Vsbphi2*SqRt[-Vsbphi2];
            t2 ← Vdx - Vb + parms[Phi2];
            t2 ← t2*SqRt[-t2];
            Id ← parms[hKp]*((2.0*(Vg - parms[Phi2fb]) - Vs - Vdx)*
                             (Vdx - Vs) - parms[ftKb]*(t1 - t2))
          END
        ELSE Id ← parms[Ioff];

        IF Id < 0.0 THEN
          IF Id < -parms[Ioff]THEN
            SIGNAL modelError["Negative drain current"]
          ELSE Id ← parms[Ioff];

        t1 ← Vg - Vb - parms[Vfb];
        IF t1 >= 0.0 THEN results[CgbX] ← parms[Cox]
        ELSE results[CgbX] ← 0.5*parms[Kb]*parms[Cox]/SqrtV;

        Vgdt ← Vg - Vd - Vt;
        t1 ← Vgst + Vgdt;
        t1 ← t1*t1;
        IF Vgst >= 0.0 THEN Cgs ← parms[Cgsdp]
        ELSE
          IF Vgdt > 0.0 THEN Cgs ← parms[Cgsd0] + parms[Cgsdp]
          ELSE Cgs ← parms[Cgsdp] + parms[Cgsd0]*(1.0 - Vgdt*Vgdt/t1);
        IF Vgdt >= 0.0 THEN Cgd ← parms[Cgsdp]
        ELSE Cgd ← parms[Cgsdp] + parms[Cgsd0]*(1.0 - Vgst*Vgst/t1);

        IF reverse THEN
          BEGIN
            results[IdX]  ← -Id;
            results[CgsX] ← Cgd;
            results[CgdX] ← Cgs
          END
        ELSE
          BEGIN
            results[IdX]  ← Id;
            results[CgsX] ← Cgs;
            results[CgdX] ← Cgd
          END
      END;

    XMOS:  spModelDefs.model =
      BEGIN
        VgX:  CARDINAL = 0;
        VsX:  CARDINAL = 1;
        VdX:  CARDINAL = 2;
        VbX:  CARDINAL = 3;

        Beta0:   CARDINAL = 0;
        Theta0:  CARDINAL = 1;
        Vth0:    CARDINAL = 2;
        v1:      CARDINAL = 3;
        v2:      CARDINAL = 4;
        TwoPhi:  CARDINAL = 5;
        sq2Phi:  CARDINAL = 6;
        A:       CARDINAL = 7;
        Az:      CARDINAL = 8;
        B:       CARDINAL = 9;
        Bz:      CARDINAL = 10;
        C0:      CARDINAL = 11;
        C0z:     CARDINAL = 12;
        C1:      CARDINAL = 13;
        C2:      CARDINAL = 14;
        Cox:     CARDINAL = 15;
        Cgxp:    CARDINAL = 16;
        Cgx0:    CARDINAL = 17;
        Cdep:    CARDINAL = 18;
        Vfb:     CARDINAL = 19;
        Io:      CARDINAL = 20;
        kTq:     CARDINAL = 21;
        expMax:  CARDINAL = 22;
        leak:    CARDINAL = 23;

        IdX:    CARDINAL = 0;
        CgbX:   CARDINAL = 1;
        CgsX:   CARDINAL = 2;
        CgdX:   CARDINAL = 3;
        Sjunc:  CARDINAL = 4;
        Djunc:  CARDINAL = 5;

        reverse:  BOOLEAN ← FALSE;
        Vg, Vs, Vd, Vb, Id:  REAL;
        Vds, Vth, Vgsth, VdsSAT, beta, IdSAT, At, As, z:  REAL;
        Cgs, Cgd, RCdep, V1, V2, Fs, Fd:  REAL;
        t1, t2:  REAL;

        Vg ← args[VgX];
        Vb ← args[VbX];
        Vd ← args[VdX];
        Vs ← args[VsX];

        t1 ← Vb - Vs;
        IF t1 >= 0.0 THEN
          BEGIN
            t1 ← t1/parms[kTq];
            IF t1 >= parms[expMax] THEN
              SIGNAL spModelDefs.Retreat["XMOS -- Source."];
            results[Sjunc] ← parms[Io]*(RealFns.Exp[t1] - 1.0)
          END
        ELSE results[Sjunc] ← -parms[Io];

        t2 ← Vb - Vd;
        IF t2 >= 0.0 THEN
          BEGIN
            t2 ← t2/parms[kTq];
            IF t2 >= parms[expMax] THEN
              SIGNAL spModelDefs.Retreat["XMOS -- Drain."];
            results[Djunc] ← parms[Io]*(RealFns.Exp[t2] - 1.0)
          END
        ELSE results[Djunc] ← -parms[Io];

        IF Vs > Vd THEN
          BEGIN
            reverse ← TRUE;
            t1 ← Vs;
            Vs ← Vd;
            Vd ← t1
          END;
        Vds ← Vd - Vs;

        z     ← SqRt[Vs - Vb + TwoPhi];
        t1    ← z - parms[sq2Phi];
        Vth   ← parms[Vth0] + (parms[v1] + parms[v2]*t1)*t1;
        Vgsth ← Vg - Vs - Vth;

        V1 ← Vgsth + 0.5*(Vth - parms[Vfb]);
        IF Vgsth > 0.0 THEN
          BEGIN
            At ← parms[A] + parms[Az]/z + (parms[B] + parms[Bz]/z)*Vgsth;
            VdsSAT ← Vgsth/At;
            beta ← parms[Beta0]/(1.0 + parms[Theta0]*Vgsth);
            IF Vds > VdsSAT THEN
              BEGIN
                IdSAT ← 0.5*beta*At*VdsSAT*VdsSAT;
                As    ← 1.0/SqRt[parms[C0] + parms[C0z]/z +
                                 (parms[C1] + parms[C2]*IdSAT)*IdSAT];
                Id    ← IdSAT/(1.0 - As*SqRt[Vds - VdsSAT])
              END
            ELSE Id ← beta*Vds*(Vgsth - 0.5*At*Vds);
            IF V1 > Vds THEN
              BEGIN
                t1 ← 2.0*V1 - Vds;
                t1 ← t1*t1;
                Fs ← V1*(3.0*V1 - 2.0*Vds)/t1;
                t2 ← V1 - Vds;
                Fd ← 3.0*t2*t2/t1;
                V2 ← Vth - parms[Vfb] + 0.5*Vds + Vs - Vb
              END
            ELSE
              BEGIN
                Fs ← 1.0;
                Fd ← 0.0;
                V2 ← Vth - parms[Vfb] + 0.5*V1 + Vs - Vb
              END
          END
        ELSE
          BEGIN
            Fs ← IF V1 < 0.0 THEN 0.0
                 ELSE 2.0*V1/(Vth - parms[Vfb]);
            Fd ← 0.0;
            V2 ← IF V1 > 0.0 THEN
                   IF V1 > Vds THEN
                      V1 + 0.5*(Vth - parms[Vfb] + Vds) + Vs - Vb
                   ELSE
                      Vds + 0.5*(Vth - parms[Vfb] + Vds) + Vs - Vb
                 ELSE
                   IF Vg - Vb < parms[Vfb] THEN 0.0
                   ELSE Vg - Vb - parms[Vfb];
            Id ← parms[Io] + parms[leak]*Vds
          END;

        Cgs ← Fs*parms[Cgx0];
        Cgd ← Fd*parms[CgdX];
        t1 ← Cox - Cgs - Cgd;
        IF t1 <= 0.0 THEN modelError["Bulk capacitance -- XMOS"];
        RCdep ← SqRt[V2]/parms[Cdep];
        results[CgbX] ← t1/(1.0 + RCdep*t1);

        IF reverse THEN
          BEGIN
            results[IdX]  ← -Id;
            results[CgsX] ← Cgd;
            results[CgdX] ← Cgs
          END
        ELSE
          BEGIN
            results[IdX]  ← Id;
            results[CgsX] ← Cgs;
            results[CgdX] ← Cgd
          END
      END;

    dcDiode:  spModelDefs.model =
      BEGIN
        Vd, e:  REAL;

        Va:  CARDINAL = 0;
        Vc:  CARDINAL = 1;

        Io:      CARDINAL = 0;
        kTq:     CARDINAL = 1;
        expMax:  CARDINAL = 2;

        Vd ← args[Va] - args[Vc];
        IF Vd + parms[kTq] <= 0.0 THEN results[0] ← -parms[Io]
        ELSE
          BEGIN
            e ← Vd/parms[kTq];
            IF e >= parms[expMax] THEN
              SIGNAL spModelDefs.Retreat["dcDiode."];
            results[0] ← parms[Io]*(RealFns.Exp[e] - 1.0)
          END
      END;

    CSIM: spModelDefs.model= BEGIN
    -- args
      VgX: CARDINAL= 0;
      VsX: CARDINAL= 1;
      VdX: CARDINAL= 2;
      VbX: CARDINAL= 3;
    -- parms
      vfb: CARDINAL= 0;
      twoPhiF: CARDINAL= 1;
      k1: CARDINAL= 2;
      k2: CARDINAL= 3;
      eta: CARDINAL= 4;
      beta0: CARDINAL= 5;
      u0: CARDINAL= 6;
      u1: CARDINAL= 7;
      -- above are the 8 electrical parameters of CSIM model
      -- cf. Scharfetter's memo of 12/30/82.
      MCgbo: CARDINAL= 8;
      Cgbo23rds: CARDINAL= 9;
      Cov: CARDINAL= 10;
      Cbso: CARDINAL= 11;
      Cbdo: CARDINAL= 12;
      phiB: CARDINAL= 13;
      -- above 6 parameters are used for capacitances calculation.
      -- cf. MOSAID model by Dick Foss
      Ios: CARDINAL= 14;
      Iod: CARDINAL= 15;
      kTq: CARDINAL= 16;
      expMax: CARDINAL= 17;
      NoCap: CARDINAL= 18;
      nChannel: CARDINAL= 19; -- -1: pType; else nType;
   -- results       
      IdX : CARDINAL= 0;
      CgbX: CARDINAL= 1;
      CgsX: CARDINAL= 2;
      CgdX: CARDINAL= 3;
      CbsX: CARDINAL= 4;
      CbdX: CARDINAL= 5;
      Ibs: CARDINAL= 6;
      Ibd: CARDINAL= 7;
    	
      Vg, Vs, Vd, Vb, Vds, Vgs, Vbs, Vbd: REAL;
      Id, Cgs, Cgd, Cbs, Cbd: REAL;
      reverse: BOOLEAN← FALSE;
      t1, TwoPhiFMVbs, Sqrt2PhiFMVbs, EtaVds,
        VgsMVth, VgdMVth, Vth: REAL;
    	
      IF parms[nChannel]= -1 THEN 
        {Vg← -args[VgX]; Vb← -args[VbX]; 
         Vd← -args[VdX]; Vs← -args[VsX]}
      ELSE  
        {parms[nChannel]← 1;
         Vg← args[VgX]; Vb← args[VbX]; Vd← args[VdX]; Vs← args[VsX]};
        
      -- get Ibs
      IF Vb > Vs THEN BEGIN
        t1← (Vb-Vs)/parms[kTq];
        IF t1 >= parms[expMax] THEN SIGNAL spModelDefs.Retreat[
          "CSIM -> Vbs forward biased too much"];
        results[Ibs]← parms[nChannel]*parms[Ios]*(RealFns.Exp[t1]- 1.0);
      END ELSE results[Ibs]← -parms[nChannel]*parms[Ios];
        
      -- get Ibd  
      IF Vb > Vd THEN BEGIN
        t1← (Vb-Vd)/parms[kTq];
        IF t1 >= parms[expMax] THEN SIGNAL spModelDefs.Retreat[
          "CSIM -> Vbd forward biased too much"];
        results[Ibd]← parms[nChannel]*parms[Iod]*(RealFns.Exp[t1]- 1.0);
      END ELSE results[Ibd]← -parms[nChannel]*parms[Iod];
          
      IF Vs > Vd THEN {reverse← TRUE; t1← Vs; Vs← Vd; Vd← t1};
          
      Vds← Vd - Vs; Vgs← Vg - Vs; Vbs← Vb - Vs; Vbd← Vb - Vd;
      
      -- if it ever reaches here, Vbs must be less than 2phif.
      -- cf. expMax calculation in ThymeBasics.thy
      TwoPhiFMVbs← parms[twoPhiF]-Vbs;
      Sqrt2PhiFMVbs← RealFns.SqRt[TwoPhiFMVbs];
      EtaVds← parms[eta]*Vds;
      
      Vth← parms[vfb]+parms[twoPhiF]+parms[k1]*Sqrt2PhiFMVbs-EtaVds
           -parms[k2]*TwoPhiFMVbs;
      VgsMVth← Vgs - Vth;
        
      IF VgsMVth <= 0.0 THEN Id← 0
      ELSE BEGIN -- Vgs>Vth:
             
        g: REAL← 1.0-1.0/(1.744+0.8364*TwoPhiFMVbs);
        a: REAL← 1.0 + 0.5*g*parms[k1]/Sqrt2PhiFMVbs;
        alpha: REAL← a*(1.0 + parms[u1]*VgsMVth); -- a(1+u1*(Vgs-Vth))
        VdSat: REAL← VgsMVth/alpha;
        BetaP: REAL← parms[nChannel]*parms[beta0]
                     /(1.0+parms[u0]*VgsMVth); -- beta0/(1+u0(Vgs-Vth))
         
        Id← IF Vds<VdSat THEN BetaP*Vds*(VgsMVth - 0.5*alpha*Vds)
                         ELSE BetaP*VgsMVth*VgsMVth/(2.0*alpha);
            
      END; -- end of Vgs>Vth
      
      IF parms[NoCap]#0 THEN {
        results[CgbX]← results[CgsX]← results[CgdX]← 
          results[CbsX]← results[CbdX]← 0; 
        results[IdX]← IF reverse THEN -Id ELSE Id;
        RETURN};
      
      VgdMVth← Vg - Vd - Vth;

      -- Cgb:
      results[CgbX]← IF VgsMVth <= 0 THEN parms[MCgbo] ELSE 0;
    		     
      -- Cgs and Cgd:
      t1← VgsMVth+VgdMVth;
      t1← t1*t1; -- (Vgs + Vgd - 2Vth)↑2
      Cgs← IF VgsMVth <= 0.0 THEN parms[Cov]
           ELSE IF VgdMVth <= 0.0 THEN parms[Cgbo23rds] + parms[Cov]
           ELSE parms[Cov] + 
                parms[Cgbo23rds]*(1.0 - VgdMVth*VgdMVth/t1);
      Cgd← IF VgdMVth <= 0.0 THEN parms[Cov]
           ELSE parms[Cov] + 
                parms[Cgbo23rds]*(1.0 - VgsMVth*VgsMVth/t1);
           
      -- Cbs and Cbd
      -- as phiB>2phif, the following should never blow up.
      Cbs← parms[Cbso]/SqRt[1.0 - Vbs/parms[phiB]];
      Cbd← parms[Cbdo]/SqRt[1.0 - Vbd/parms[phiB]];
           
      IF reverse THEN BEGIN
        results[IdX] ← -Id;
        results[CgsX]← Cgd;
        results[CgdX]← Cgs;
        results[CbsX]← Cbd;
        results[CbdX]← Cbs;
      END ELSE BEGIN
        results[IdX] ← Id;
        results[CgsX]← Cgs;
        results[CgdX]← Cgd;
        results[CbsX]← Cbs;
        results[CbdX]← Cbd;
      END;
        
    END; -- CSIM

    acDiode: spModelDefs.model= BEGIN
      Vc: CARDINAL= 0; -- cathode voltage index
      Va: CARDINAL= 1; -- anode voltage index
  
      Co: CARDINAL= 0;
      Pb: CARDINAL= 1;
      Vt: CARDINAL= 2;
      Io: CARDINAL= 3;
      expMax: CARDINAL= 4;
  
      C: CARDINAL= 0;
      I: CARDINAL= 1;
  
      t1: REAL;
      V: REAL← args[Va] - args[Vc];
      IF V > 0 THEN BEGIN -- forward biased
        t1← V/parms[Vt];
        IF t1 >= parms[expMax] THEN SIGNAL spModelDefs.Retreat[
          "acDiode forward biased too much"];
        results[I]← -parms[Io]*(RealFns.Exp[t1]- 1.0);
      END ELSE results[I]← parms[Io];
      results[C]← parms[Co]/SqRt[1.0 - V/parms[Pb]];
      END; -- acDiode
          
    initSqRt[];
    
    spModelDefs.EnterModels["CSIM",       CSIM,       4, 20, 8];
    spModelDefs.EnterModels["acDiode",    acDiode,    2,  5, 2];
    spModelDefs.EnterModels["dcDiode",    dcDiode,    2,  3, 1];
    spModelDefs.EnterModels["NFET",       NFET,       4, 16, 6];
    spModelDefs.EnterModels["PFET",       PFET,       4, 16, 6];
    spModelDefs.EnterModels["XMOS",       XMOS,       4, 24, 6];
  END.
4/29/83:- 
  original: [Cherry]<Barth>Thyme>1.97>spModels.mesa
  changes: added and "enter"ed CSIM and acDiode model.