-- File: [Cherry]<Thyme>System>C03>Level2Model.mesa
-- Last editted by SChen  February 19, 1984  3:17 PM

DIRECTORY AltoDefs, Real, RealFns, spModelDefs;
Level2Model: 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; -- SqRt

  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; -- initSqRt

  LimitVgs: PROC[vgs, vgsOld, vto: REAL] RETURNS[vgsNew: REAL]= BEGIN
    vtstHigh, vtstLow, vtoPlus35, deltaV: REAL;
    vtstHigh_ 2.0 + 2.0*ABS[vgsOld-vto];
    vtstLow_ vtstHigh/2.0 + 2.0;
    vtoPlus35_ vto + 3.5;
    deltaV_ vgs-vgsOld;
    vgsNew_ vgs;
    IF vgsOld < vto THEN BEGIN
      -- OFF
      IF deltaV > 0 THEN BEGIN
        vTemp: REAL_ vto + 0.5;
        IF vgsNew > vTemp THEN vgsNew_ vTemp; 
      END;
    END ELSE IF vgsOld < vtoPlus35 THEN BEGIN
      -- MIDDLE
      vgsNew_ IF deltaV > 0 THEN MIN[vgsNew, vto + 4.0] -- increasing
              ELSE MAX[vgsNew, vto - 0.5]; -- decreasing
    END ELSE BEGIN
      -- ON
      IF deltaV > 0 THEN BEGIN -- staying on
        IF deltaV > vtstHigh THEN vgsNew_ vgsOld + vtstHigh;
      END ELSE BEGIN -- going off
        IF vgsNew < vtoPlus35 THEN vgsNew_ MAX[vgsNew, vto + 2.0]
        ELSE IF -deltaV > vtstLow THEN vgsNew _ vgsOld - vtstLow;
      END;
    END;
  END; -- LimitVgs

  LimitVds: PROC[vds, vdsOld: REAL] RETURNS[vdsNew: REAL]= BEGIN
    vdsNew_ vds;
    IF vdsOld < 3.5 THEN BEGIN
      vdsNew_ IF vdsNew <= vdsOld THEN MAX[vdsNew, -5]
              ELSE MIN[vdsNew, 4];
    END ELSE BEGIN
      IF vdsNew <= vdsOld THEN 
        {IF vdsNew < 3.5 THEN vdsNew_ MAX[vdsNew, 2]}
      ELSE vdsNew_ MIN[vdsNew, 2.0 + 3.0*vdsOld];
    END;
  END; -- LimitVds

  LimitVbs: PROC[vbs, vbsOld, vt, vcrit: REAL] 
    RETURNS[vbsNew: REAL]= BEGIN
    vbsNew_ vbs;
    IF vbsNew > vcrit THEN BEGIN
      vLimit: REAL_ vt + vt;
      deltaV: REAL_ vbsNew - vbsOld;
      IF ABS[deltaV] > vLimit THEN BEGIN
        IF vbsOld <= 0 THEN vbsNew_ vt*RealFns.Ln[vbsNew/vt]
        ELSE BEGIN
          arg: REAL_ 1.0 + deltaV/vt;
          vbsNew_ IF arg <= 0 THEN vcrit 
                              ELSE vbsOld + vt*RealFns.Ln[arg];
        END;
      END;
    END;
  END; -- LimitVbs

  Baum: PROC[VgsX, Vbin, Eta, Phi, Vbs, Vmax, 
    L, UEffective, gammaD, SArg3: REAL] RETURNS[Vdsat: REAL]= BEGIN

    v1, v2, xv, a1, b1, c1, d1, a3, b3, A, B, C, R, S, 
    r3, s2, P, p0, p2, fi, y3, delta4, xValid: REAL_ 0.0;
    x4, poly4: ARRAY[1..8] OF REAL;
    a4, b4, Sig1, Sig2: ARRAY[1..4] OF REAL;
    Icount, Jcount: CARDINAL_ 0;

    Icount_ Jcount_ 0;
    x4_ poly4_ ALL[0.0];
    a4_ b4_ ALL[0.0];
    Sig1_ [1.0, -1.0, 1.0, -1.0];
    Sig2_ [1.0, 1.0, -1.0, -1.0];
        
    v1_ (VgsX-Vbin)/Eta+Phi-Vbs; 
    v2_ Phi-Vbs;
    xv_ Vmax*L/UEffective;
         
    a1_ gammaD/0.75; b1_ -2.0*(v1+xv); c1_ -2.0*gammaD*xv; 
    d1_ 2.0*v1*(v2+xv)-v2*v2-4.0/3.0*gammaD*SArg3;
        
    A_ -b1; B_ a1*c1-4.0*d1; C_ -d1*(a1*a1-4.0*b1)-c1*c1;
       
    R_ -A*A/3.0+B; r3_ R*R*R;
    S_ 2.0*A*A*A/27.0-A*B/3.0+C; s2_ S*S;
       
    P_ s2/4.0+r3/27.0; p0_ ABS[P]; p2_ SqRt[p0];
       
    -- get y3 value.
    IF P<0.0 THEN BEGIN 
      ro: REAL;
      ro_ SqRt[s2/4.0+p0]; 
      ro_ RealFns.Ln[ro]/3.0; 
      ro_ RealFns.Exp[ro]; 
      fi_ RealFns.ArcTan[-2.0*p2/S, 1]; -- arcTan(-2.0*p2/S); 
      y3_ 2.0*ro*RealFns.Cos[fi/3.0]-A/3.0;
    END ELSE BEGIN -- 550 
      p3, p4: REAL;
      p3_ RealFns.Exp[ RealFns.Ln[ABS[-S/2.0+p2]]/3.0 ];
      p4_ RealFns.Exp[ RealFns.Ln[ABS[-S/2.0-p2]]/3.0 ];
      y3_ p3+p4-A/3.0;
    END;
    
    a3_ SqRt[a1*a1/4.0-b1+y3]; 
    b3_ SqRt[y3*y3/4.0-d1];
    
    -- from a3, b3, Sig1, Sig2, y3 -> a4, b4 -> delta4 -> x4
    FOR I:CARDINAL IN[1..4] DO
      a4[I]_ a1/2.0+Sig1[I]*a3; 
      b4[I]_ y3/2.0+Sig2[I]*b3; 
      delta4_ a4[I]*a4[I]/4.0-b4[I];
      IF delta4>=0.0 THEN BEGIN
        sqrtd4: REAL;
        sqrtd4_ SqRt[delta4];
        Icount_ Icount+1; 
        x4[Icount]_ a4[I]/(-2.0)+sqrtd4; 
        Icount_ Icount+1; 
        x4[Icount]_ a4[I]/(-2.0)-sqrtd4;
      END; 
    ENDLOOP;
    
    FOR J:CARDINAL IN [1..Icount] DO
      tmp: REAL; 
      IF x4[J]<=0.0 AND J<Icount THEN LOOP; 
      
      tmp_ (x4[J]+a1);
      tmp_ (tmp*x4[J]+b1);
      tmp_ (tmp*x4[J]+c1);
      poly4[J]_ (tmp*x4[J]+d1);
      IF ABS[poly4[J]]>1.0E-6 AND J<Icount THEN LOOP; 
      Jcount_ Jcount+1; 
      IF Jcount<=1 THEN xValid_ x4[J];
      IF x4[J]<=xValid THEN xValid_ x4[J];
    ENDLOOP;
    
    IF Jcount>0 THEN Vdsat_ xValid*xValid+Vbs-Phi
    ELSE SIGNAL spModelDefs.Retreat[
      "Error in immediate velocity saturation model"];
  END; -- Baum

  Grove: PROC[gammaD, VgsX, Vbin, Eta, Phi, Vbs: REAL]
    RETURNS[Vdsat: REAL]= BEGIN
    gammaD2, ArgV: REAL;
    gammaD2_ gammaD*gammaD;
    ArgV_ (VgsX-Vbin)/Eta+Phi-Vbs;
    IF ArgV<=0.0 THEN Vdsat_ 0.0
    ELSE BEGIN
      Arg: REAL;
      Arg_ SqRt[1.0+4.0*ArgV/gammaD2];
      Vdsat_ (VgsX-Vbin)/Eta+gammaD2*(1.0-Arg)/2.0;
      IF Vdsat<0.0 THEN Vdsat_0.0;
    END;
  END; -- Grove

  JcnCurrent: PROC[vt, v, isat: REAL] RETURNS[i: REAL]= BEGIN
    t1: REAL_ v/vt;
    IF v > 0 THEN BEGIN
      t1_ MIN[t1, 85.0];
      i_ isat*(RealFns.Exp[t1]- 1.0);
    END ELSE i_ IF t1 <= -1.0 THEN -isat ELSE isat*t1;
  END; -- JcnCurrent

  Pair: TYPE= RECORD[f, df: REAL];

  FandDF: PROC[v, p, sp, sp3: REAL] RETURNS[pair: Pair]= BEGIN
    IF v <= 0 THEN BEGIN
      pair.f_ SqRt[p-v]; 
      pair.df_ -0.5/pair.f;
    END ELSE BEGIN -- 130
      pair.f_ sp/(1.0+0.5*v/p); 
      pair.df_ -0.5*pair.f*pair.f/sp3; 
    END;
  END; -- FandDF
        
  SpiceLevel2: spModelDefs.model= BEGIN

    -- args
    VdX: CARDINAL= 0;
    VgX: CARDINAL= 1;
    VsX: CARDINAL= 2;
    VbX: CARDINAL= 3;

    -- oldArgs
    oldVdX: CARDINAL= 0;
    oldVgX: CARDINAL= 1;
    oldVsX: CARDINAL= 2;
    oldVbX: CARDINAL= 3;

    -- parms
    LX: CARDINAL= 0;
    WX: CARDINAL= 1;
    XdX: CARDINAL= 2;
    XjX: CARDINAL= 3;
    IssatX: CARDINAL= 4;
    IdsatX: CARDINAL= 5;
    BetaX: CARDINAL= 6;
    PhiX: CARDINAL= 7;
    SqRtPhiX: CARDINAL= 8;
    SqRtPhi3X: CARDINAL= 9;
    PbX: CARDINAL= 10;
    SqRtPbX: CARDINAL= 11;
    VbpX: CARDINAL= 12;
    VbiX: CARDINAL= 13;
    GammaX: CARDINAL= 14;
    LambdaX: CARDINAL= 15;
    VtX: CARDINAL= 16;
    UoX: CARDINAL= 17;
    VmaxX: CARDINAL= 18;
    QNfsX: CARDINAL= 19;
    NsubX: CARDINAL= 20;
    UexpX: CARDINAL= 21;
    NeffX: CARDINAL= 22;
    EtaX: CARDINAL= 23;
    FactorX: CARDINAL= 24;
    CoxX: CARDINAL= 25;
    CovGSX: CARDINAL= 26;
    CovGDX: CARDINAL= 27;
    CgbMX: CARDINAL= 28;
    Cgb23rdsX: CARDINAL= 29;
    CbsBottomX: CARDINAL= 30;
    CbdBottomX: CARDINAL= 31;
    CbsSWX: CARDINAL= 32;
    CbdSWX: CARDINAL= 33;
    CbsFwd1X: CARDINAL= 34;
    CbsFwd2X: CARDINAL= 35;
    CbdFwd1X: CARDINAL= 36;
    CbdFwd2X: CARDINAL= 37;
    MjX: CARDINAL= 38;
    MjswX: CARDINAL= 39;
    FcPbX: CARDINAL= 40;
    VonX: CARDINAL= 41;
    minVX: CARDINAL= 42;
    maxVX: CARDINAL= 43;
    SignX: CARDINAL= 44;
    
    reverse: BOOLEAN;
    sign, Vd, Vg, Vs, Vb, oldVd, oldVg, oldVs, oldVb,
    Vds, Vgs, Vbs, Vgd, Vbd, 
    oldVds, oldVgs, oldVbs, oldVgd, oldVbd, 
    Id, Cgb, Cgs, Cgd, Cbs, Cbd, t1,
    Vbin, XdBarg, XdSarg, ArgSS, ArgSD, DBArgS, DBArgD,
    ArgD, ArgS, gamaSD, DBXWD, DBXWS, gammaD, 
    DgdDvb, VgsT, VgdT,
    SArg3, body, UFactor, UEffective, 
    VgsX, bodyS, clFactor, Vdsat, ArgG, beta1, Vcrit: REAL;
    SArg, BArg: Pair;
    
    LimitNodeVoltage: PROC[oldV: REAL]
      RETURNS[newV: REAL]= BEGIN
      newV_ MIN[parms[maxVX], MAX[parms[minVX], oldV]];
    END; -- LimitNodeVoltage

    Id_ 0;
    sign_ parms[SignX];
    Vd_ LimitNodeVoltage[sign*args[VdX]];
    Vg_ LimitNodeVoltage[sign*args[VgX]]; 
    Vs_ LimitNodeVoltage[sign*args[VsX]];
    Vb_ LimitNodeVoltage[sign*args[VbX]];
    oldVd_ LimitNodeVoltage[sign*args[oldVdX]];
    oldVg_ LimitNodeVoltage[sign*args[oldVgX]]; 
    oldVs_ LimitNodeVoltage[sign*args[oldVsX]];
    oldVb_ LimitNodeVoltage[sign*args[oldVbX]];
    
    Vds_ Vd - Vs; Vgs_ Vg - Vs; Vbs_ Vb - Vs; 
    Vgd_ Vg - Vd; Vbd_ Vb - Vd;
    oldVds_ oldVd - oldVs; oldVgs_ oldVg - oldVs; oldVbs_ oldVb - oldVs;
    oldVgd_ oldVg - oldVd; oldVbd_ oldVb - oldVd;

    IF oldVds >= 0 THEN BEGIN
      Vgs_ LimitVgs[Vgs, oldVgs, parms[VonX]];
      Vds_ Vgs - Vgd;
      Vds_ LimitVds[Vds, oldVds];
      Vgd_ Vgs - Vds;
    END ELSE BEGIN
      Vgd_ LimitVgs[Vgd, oldVgd, parms[VonX]];
      Vds_ Vgs - Vgd;
      Vds_ -LimitVds[-Vds, -oldVds];
      Vgs_ Vgd + Vds;
    END;
    IF Vds >= 0 THEN BEGIN
      Vcrit_ parms[VtX]*RealFns.Ln[parms[VtX]/(1.414*parms[IssatX])];
      Vbs_ LimitVbs[Vbs, oldVbs, parms[VtX], Vcrit];
      Vbd_ Vbs - Vds;
      args[VdX]_ sign*(Vs + Vds);
      args[VgX]_ sign*(Vs + Vgs);
      args[VbX]_ sign*(Vs + Vbs);
    END ELSE BEGIN
      Vcrit_ parms[VtX]*RealFns.Ln[parms[VtX]/(1.414*parms[IdsatX])];
      Vbd_ LimitVbs[Vbd, oldVbd, parms[VtX], Vcrit];
      Vbs_ Vbd + Vds;
      args[VsX]_ sign*(Vd - Vds);
      args[VgX]_ sign*(Vd + Vgd);
      args[VbX]_ sign*(Vd + Vbd);
    END;

    -- get Ibs and Ibd
    results[6]_ sign*JcnCurrent[parms[VtX], Vbs, parms[IssatX]];
    results[7]_ sign*JcnCurrent[parms[VtX], Vbd, parms[IdsatX]];
      
    IF Vds >= 0 THEN reverse_ FALSE
    ELSE {reverse_ TRUE; Vds_ -Vds;
          t1_ Vbs; Vbs_ Vbd; Vbd_ t1;
          t1_ Vgs; Vgs_ Vgd; Vgd_ t1;
          };

    -- 100
    SArg_ FandDF[Vbs, parms[PhiX], parms[SqRtPhiX], parms[SqRtPhi3X]];
    BArg_ FandDF[Vbd, parms[PhiX], parms[SqRtPhiX], parms[SqRtPhi3X]];
    
    -- 200, Calculate parms[VonX], narrow channel effect
    Vbin_ parms[VbiX]+parms[FactorX]*(parms[PhiX]-Vbs);
    XdBarg_ parms[XdX]*BArg.f; 
    XdSarg_ parms[XdX]*SArg.f; 

    -- short channel effect with Vds # 0
    ArgSS_ ArgSD_ DBArgS_ DBArgD_ 0.0;
    IF parms[XjX] > 0.0 THEN BEGIN
      ArgXS, ArgXD: REAL;
      ArgXS_ 1.0+2.0*XdSarg/parms[XjX];
      ArgS_ SqRt[ArgXS]; 
      ArgSS_ 0.5*parms[XjX]/parms[LX]*(ArgS-1.0);
      ArgXD_ 1.0+2.0*XdBarg/parms[XjX];
      ArgD_ SqRt[ArgXD]; 
      ArgSD_ 0.5*parms[XjX]/parms[LX]*(ArgD-1.0);
    END;
    
    --(205)
    gamaSD_ parms[GammaX]*(1.0-ArgSS-ArgSD);
    DBXWD_ parms[XdX]*BArg.df; 
    DBXWS_ parms[XdX]*SArg.df; 
    IF parms[XjX] > 0.0 THEN BEGIN
      DBArgS_ 0.5/parms[LX]*DBXWS/ArgS;
      DBArgD_ 0.5/parms[LX]*DBXWD/ArgD;
    END;
    
    DgdDvb_ -parms[GammaX]*(DBArgS+DBArgD); 
    
    --220 
    parms[VonX]_ Vbin+gamaSD*SArg.f;
    Vdsat_ 0.0;

    --225
    BEGIN
      IF parms[QNfsX] = 0.0 THEN BEGIN -- 230
        VgsT_ Vgs-parms[VonX];
        IF VgsT<=0.0 THEN GOTO done; -- (Id=0) 1050 
      END ELSE BEGIN
        CdOnCo, Xn: REAL;
        CdOnCo_ -(gamaSD*SArg.df+DgdDvb*SArg.f)+parms[FactorX];
        Xn_ 1.0+parms[QNfsX]/parms[CoxX]+CdOnCo; 
        parms[VonX]_ parms[VonX]+parms[VtX]*Xn; 
        ArgG_ 1.0/(parms[VtX]*Xn);
        VgsT_ Vgs-parms[VonX];
      END; 
        
      -- 300
      SArg3_ SArg.f*SArg.f*SArg.f;
      gammaD_ gamaSD; 
      body_ BArg.f*BArg.f*BArg.f-SArg3; 
      
      -- 400, EVALUATE EFFECTIVE MOBILITY AND ITS DERIVATIVES
      IF VgsT > parms[VbpX] THEN BEGIN
        UFactor_ RealFns.Exp[ parms[UexpX]*RealFns.Ln[parms[VbpX]/VgsT] ]; 
        UEffective_ parms[UoX]*UFactor; 
      END ELSE BEGIN  -- 410
        UFactor_ 1.0; 
        UEffective_ parms[UoX]; 
      END;
      
      -- 500 EVALUATE SATURATION VOLTAGE AND ITS DERIVATIVES
      -- ACCORDING TO GROVE-FROHMAN EQUATION
      VgsX_ Vgs;
      gammaD_ gamaSD/parms[EtaX]; 
      IF parms[QNfsX]#0.0 THEN VgsX_ MAX[Vgs, parms[VonX]];
      Vdsat_ IF gammaD <= 0.0 THEN MAX[0.0, (VgsX-Vbin)/parms[EtaX]]
      ELSE Grove[gammaD, VgsX, Vbin, parms[EtaX], parms[PhiX], Vbs]; 
      
      -- 545
      IF parms[VmaxX] > 0.0 THEN Vdsat_ Baum[
        VgsX, Vbin, parms[EtaX], parms[PhiX], Vbs, parms[VmaxX], 
        parms[LX], UEffective, gammaD, SArg3];
        
      -- 600, EVALUATE EFFECTIVE CHANNEL LENGTH AND ITS DERIVATIVES
      IF Vds # 0.0 THEN BEGIN
        lFactor, ArgV, xdv, lv, ls, SArgV, BSArg: REAL;
         
        gammaD_ gamaSD; 
        BSArg_ IF Vbs <= Vdsat THEN SqRt[Vdsat-Vbs+parms[PhiX]]
               ELSE parms[SqRtPhiX]/(1.0+0.5*(Vbs-Vdsat)/parms[PhiX]);
          
        bodyS_ BSArg*BSArg*BSArg-SArg3; 
        IF parms[VmaxX] <=0.0 THEN BEGIN
          IF parms[NsubX] # 0.0 AND parms[LambdaX] <= 0.0 THEN BEGIN
            Arg: REAL;
            ArgV_ (Vds-Vdsat)/4.0;
            SArgV_ SqRt[1.0+ArgV*ArgV];
            Arg_ SqRt[ArgV+SArgV]; 
            lFactor_ parms[XdX]/(parms[LX]*Vds);
            parms[LambdaX]_ lFactor*Arg;
          END;
        END ELSE BEGIN -- 603 
          ArgV_ (VgsX-Vbin)/parms[EtaX]-Vdsat;
          xdv_ parms[XdX]/SqRt[parms[NeffX]]; 
          lv_ parms[VmaxX]*xdv/(2.0*UEffective); 
          IF parms[NsubX] # 0.0 AND parms[LambdaX] <= 0.0 THEN BEGIN
            ArgV_ MAX[Vds-Vdsat, 0.0]; 
            ls_ SqRt[lv*lv+ArgV]; 
            lFactor_ xdv/(parms[LX]*Vds); 
            parms[LambdaX]_ lFactor*(ls-lv);
          END;
        END; 
      END; -- block 600
      
        -- 620, LIMIT CHANNEL SHORTENING AT PUNCH-THROUGH 
      BEGIN
        leffective, deltaL, xwb, ld: REAL;
        xwb_ parms[XdX]*parms[SqRtPbX];
        ld_ parms[LX]-xwb;
        clFactor_ 1.0-parms[LambdaX]*Vds; 
        leffective_ parms[LX]*clFactor; 
        deltaL_ parms[LambdaX]*Vds*parms[LX];
        IF parms[NsubX]=0.0 THEN xwb_ 0.25E-6; 
        IF leffective<xwb THEN BEGIN  
          leffective_ xwb/(1.0+(deltaL-ld)/xwb);
          clFactor_ leffective/parms[LX]; 
        END;
      END;

      -- 700 EVALUATE EFFECTIVE beta (EFFECTIVE KP) 
      beta1_ parms[BetaX]*UFactor/clFactor;
 
      -- TEST FOR MODE OF OPERATION AND BRANCH APPROPRIATELY
      gammaD_ gamaSD; 
      IF Vds <= 1.0E-8 THEN GOTO done; -- (Id=0) 1050
      -- 730
      IF VgsT <= 0.0 THEN BEGIN -- else goto 900 
        -- SUBTHRESHOLD REGION, IF Vdsat<=0.0 THEN Id_ 0.0
        IF  Vdsat > 0.0 THEN BEGIN -- 830
          VdsOn, IdsOn, ArgGVgsT: REAL;
          VdsOn_ MIN[Vdsat,Vds];
          ArgGVgsT_ ArgG*VgsT;
          IF Vds > Vdsat THEN body_ bodyS;
          IdsOn_ beta1*( (parms[VonX] - Vbin -
                   parms[EtaX]*VdsOn*0.5)*VdsOn - gammaD*body/1.5 );
          -- if ArgGVgsT < -ln(IdsOn) - 34.539 then Id < 1E-15 ~ 0
          Id_ IF ArgGVgsT < (- 34.539 - RealFns.Ln[IdsOn]) THEN 0
              ELSE sign*IdsOn*RealFns.Exp[ArgGVgsT];
        END;
      
      END ELSE 
        Id_ IF Vds <= Vdsat THEN -- 900 LINEAR REGION
              sign*beta1*((Vgs-Vbin-parms[EtaX]*Vds/2.0)*Vds-gammaD*body/1.5) 
            ELSE -- 1000 SATURATION REGION
              sign*beta1*((Vgs-Vbin-parms[EtaX]*Vdsat/2.0)*Vdsat-gammaD*bodyS/1.5);
    EXITS
      done=> NULL;
    END;
      
    -- capacitances, and results
    VgdT_ Vgd - parms[VonX];

    -- Cgb:
    Cgb_ IF VgsT <= 0 THEN parms[CgbMX] ELSE 0;
                       
    -- Cgs and Cgd:

    t1_ VgsT+VgdT;
    t1_ t1*t1; -- (Vgs + Vgd - 2Von)^2
    Cgs_ IF VgsT <= 0.0 THEN parms[CovGSX]
       ELSE IF VgdT <= 0.0 THEN parms[Cgb23rdsX] + parms[CovGSX]
       ELSE parms[CovGSX] + parms[Cgb23rdsX]*(1.0 - VgdT*VgdT/t1);
    Cgd_ IF VgdT <= 0.0 THEN parms[CovGDX]
       ELSE parms[CovGDX] + parms[Cgb23rdsX]*(1.0 - VgsT*VgsT/t1);

    -- Cbs and Cbd

    IF Vbs >= parms[FcPbX] THEN -- 520
      Cbs_ parms[CbsFwd1X] + Vbs*parms[CbsFwd2X]
    ELSE BEGIN
      t1_ RealFns.Ln[1.0-Vbs/parms[PbX]];
      Cbs_ parms[CbsBottomX]*RealFns.Exp[-parms[MjX]*t1] + parms[CbsSWX]*RealFns.Exp[-parms[MjswX]*t1]
    END;
    
    IF Vbd >= parms[FcPbX] THEN -- 520
      Cbd_ parms[CbdFwd1X] + Vbd*parms[CbdFwd2X]
    ELSE BEGIN
      t1_ RealFns.Ln[1.0-Vbd/parms[PbX]];
      Cbd_ parms[CbdBottomX]*RealFns.Exp[-parms[MjX]*t1] + parms[CbdSWX]*RealFns.Exp[-parms[MjswX]*t1]
    END;
    
    -- results order: Id, Cgb, Cgs, Cgd, Cbs, Cbd, Ibs, Ibd.
    IF reverse THEN BEGIN
      results[0]_ -Id;
      results[1]_ Cgb;
      results[2]_ Cgd;
      results[3]_ Cgs;
      results[4]_ Cbd;
      results[5]_ Cbs;
    END ELSE BEGIN
      results[0]_ Id;
      results[1]_ Cgb;
      results[2]_ Cgs;
      results[3]_ Cgd;
      results[4]_ Cbs;
      results[5]_ Cbd;
    END;
  END; -- SpiceLevel2  
  
  initSqRt[];
  
  spModelDefs.EnterModels["SpiceLevel2", SpiceLevel2, 4, 45, 8];
END.
2/12/84:-
  created to support Spice level2 model.
2/19/84:-
  completed addition of Voltage limiting stmts; modified Poly4 calculation.
2/22/84:- 

  added protection for week inversion current calculation to insure that Exp(ArgG*VgsT) doesn't blow up.