-- File: [Thyme]<Thyme>System>Spice>Level2Model.mesa
-- Last editted by SChen  8-Feb-84 20:21:27

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

  SpiceLevel2: spModelDefs.model= BEGIN
    -- args
    Vd: REAL_ args[0];
    Vg: REAL_ args[1];
    Vs: REAL_ args[2];
    Vb: REAL_ args[3];

    -- parms
    l: REAL= parms[0];
    w: REAL= parms[1];
    xd: REAL= parms[2];
    xj: REAL= parms[3];
    Issat: REAL= parms[4];
    Idsat: REAL= parms[5];
    beta: REAL= parms[6];
    phi: REAL= parms[7];
    phiB: REAL= parms[8];
    vbp: REAL= parms[9];
    vbi: REAL= parms[10];
    gamma: REAL= parms[11];
    lambda: REAL_ parms[12];
    vt: REAL= parms[13];
    uo: REAL= parms[14];
    vmax: REAL= parms[15];
    QNfs: REAL= parms[16];
    nsub: REAL= parms[17];
    uexp: REAL= parms[18];
    neff: REAL= parms[19];
    eta: REAL= parms[20];
    factor: REAL= parms[21];

    Cox: REAL= parms[22];
    CovGS: REAL= parms[23];
    CovDS: REAL= parms[24];
    CgbM: REAL= parms[25];
    Cgb23rds: REAL= parms[26];

    CbsBottom: REAL= parms[27];
    CbdBottom: REAL= parms[28];
    CbsSW: REAL= parms[29];
    CbdSW: REAL= parms[30];
    CbsFwd1: REAL= parms[31];
    CbsFwd2: REAL= parms[32];
    CbdFwd1: REAL= parms[33];
    CbdFwd2: REAL= parms[34];
    FcByPb: REAL= parms[35];
    sign: REAL= parms[36];

    Vds, Vgs, Vbs, Vbd, Id, Ibs, Ibd, 
    Cgb, Cgs, Cgd, Cbs, Cbd, t1: REAL_ 0.0;
    reverse: BOOLEAN_ FALSE;
    SArg, DsargDb, SqrtPhi, SqrtPhi3, BArg, DbargDb,
    Vbin, XdBarg, XdSarg, ArgSS, ArgSD, DBArgS, DBArgD,
    ArgD, ArgS, ArgXS, ArgXD, gamaSD, DBXWD, DBXWS, gammaD, 
    DgdDvb, CdOnCo, Xn, VgsT, ArgG, SArg3, SqrtPhib, body,
    UFactor, UEffective, VgsX, ArgV, Arg, BSArg, DbsargDb, bodyS,
    clFactor, von, vdsat, beta1: REAL_ 0.0;
    
    IF sign= -1 THEN {Vd_ -Vd; Vg_ -Vg; Vs_ -Vs; Vb_ -Vb};
    
    IF Vs > Vd THEN {reverse_ TRUE; t1_ Vs; Vs_ Vd; Vd_ t1};
      
    Vds_ Vd - Vs; Vgs_ Vg - Vs; Vbs_ Vb - Vs; Vbd_ Vb - Vd;
    
    -- get Ibs
    t1_ Vbs/vt;
    IF Vbs > 0 THEN BEGIN
      t1_ MIN[t1, 85.0];
      Ibs_ sign*Issat*(RealFns.Exp[t1]- 1.0);
    END ELSE Ibs_ sign*Issat*t1;
    
    -- get Ibd  
    t1_ Vbd/vt;
    IF Vbd > 0 THEN BEGIN
      t1_ MIN[t1, 85.0];
      Ibd_ sign*Idsat*(RealFns.Exp[t1]- 1.0);
    END ELSE Ibd_ sign*Idsat*t1;
      
    -- 100 
    IF Vbs <= 0.0 THEN BEGIN 
      SArg_ SqRt[phi-Vbs]; 
      DsargDb_ -0.5/SArg;
    END ELSE BEGIN -- ( 110 )
      SqrtPhi_ SqRt[phi]; 
      SqrtPhi3_ phi*SqrtPhi;
      SArg_ SqrtPhi/(1.0+0.5*Vbs/phi); 
      DsargDb_ -0.5*SArg*SArg/SqrtPhi3; 
    END;
  
    -- 120
    IF Vds >= Vbs THEN BEGIN 
      BArg_ SqRt[phi+Vds-Vbs]; 
      DbargDb_ -0.5/BArg;
    END ELSE BEGIN -- 130
      BArg_ SqrtPhi/(1.0+0.5*(Vbs-Vds)/phi); 
      DbargDb_ -0.5*BArg*BArg/SqrtPhi3; 
    END;
    
    -- 200, Calculate von, narrow channel effect
    Vbin_ vbi+factor*(phi-Vbs);
    IF (gamma<=0.0) OR (nsub<=0.0) THEN BEGIN -- 215
      gamaSD_ gammaD_ gamma;
      DgdDvb_ 0.0;
    END ELSE BEGIN 
      XdBarg_ xd*BArg; 
      XdSarg_ xd*SArg; 
   
      -- short channel effect with Vds # 0
      ArgSS_ ArgSD_ DBArgS_ DBArgD_ 0.0;
      IF xj > 0.0 THEN BEGIN
        ArgXS_ 1.0+2.0*XdSarg/xj;
        ArgS_ SqRt[ArgXS]; 
        ArgSS_ 0.5*xj/l*(ArgS-1.0);
        ArgXD_ 1.0+2.0*XdBarg/xj;
        ArgD_ SqRt[ArgXD]; 
        ArgSD_ 0.5*xj/l*(ArgD-1.0);
      END;
      
      --(205)
      gamaSD_ gamma*(1.0-ArgSS-ArgSD);
      DBXWD_ xd*DbargDb; 
      DBXWS_ xd*DsargDb; 
      IF xj>0.0 THEN BEGIN
        DBArgS_ 0.5/l*DBXWS/ArgS;
        DBArgD_ 0.5/l*DBXWD/ArgD;
      END;
      
      DgdDvb_ -gamma*(DBArgS+DBArgD); 
    END;
    
    --220 
    von_ Vbin+gamaSD*SArg;
    vdsat_ 0.0;

    --225
    BEGIN
      IF QNfs = 0.0 THEN BEGIN -- 230
        VgsT_ Vgs-von;
        IF VgsT<=0.0 THEN GOTO done; -- (Id=0) 1050 
      END ELSE BEGIN
        CdOnCo_ -(gamaSD*DsargDb+DgdDvb*SArg)+factor;
        Xn_ 1.0+QNfs/Cox+CdOnCo; 
        von_ von+vt*Xn; 
        ArgG_ 1.0/(vt*Xn);
        VgsT_ Vgs-von;
      END; 
        
      -- 300
      SArg3_ SArg*SArg*SArg;
      SqrtPhib_ SqRt[phiB];
      gammaD_ gamaSD; 
      body_ BArg*BArg*BArg-SArg3; 
      
      -- 400, EVALUATE EFFECTIVE MOBILITY AND ITS DERIVATIVES
      IF VgsT > vbp THEN BEGIN
        UFactor_ RealFns.Exp[ uexp*RealFns.Ln[vbp/VgsT] ]; 
        UEffective_ uo*UFactor; 
      END ELSE BEGIN  -- 410
        UFactor_ 1.0; 
        UEffective_ uo; 
      END;
      
      -- 500 EVALUATE SATURATION VOLTAGE AND ITS DERIVATIVES ACCORDING TO
      -- GROVE-FROHMAN EQUATION
      VgsX_ Vgs;
      gammaD_ gamaSD/eta; 
      IF QNfs#0.0 THEN VgsX_ MAX[Vgs, von];
      IF gammaD<=0.0 THEN vdsat_ MAX[0.0, (VgsX-Vbin)/eta]
      ELSE BEGIN
        gammaD2: REAL_ gammaD*gammaD;
        ArgV_ (VgsX-Vbin)/eta+phi-Vbs;
        IF ArgV<=0.0 THEN vdsat_ 0.0
        ELSE BEGIN
          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; -- vdsat_ DMAX1(vdsat,0.0);
        END;
      END; 
      
      -- 545
      IF vmax > 0.0 THEN BEGIN
      -- EVALUATE SATURATION VOLTAGE AND ITS DERIVATIVES ACCORDING TO
      -- BAUM'S THEORY OF SCATTERING VELOCITY SATURATION 
        v1, v2, xv, a1, b1, c1, d1, a3, b3, A, B, C, R, S, r3, s2, P, p0, p2, fi, 
          y3, delta4, xValid: REAL;
        x4, poly4, a4, b4, Sig1, Sig2: ARRAY[1..4] OF REAL;
        Icount, Jcount: CARDINAL_ 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_ SqRt[s2/4.0+p0]; 
          ro_ RealFns.Ln[ro]/3.0; 
          ro_ RealFns.Exp[ro]; 
          fi_ RealFns.ArcTan[-2.0*p2, S]; -- arcTan(-2.0*p2/S); 
          y3_ 2.0*ro*RealFns.Cos[fi/3.0]-A/3.0;
        END ELSE BEGIN -- 550 
          p3: REAL_ RealFns.Exp[ RealFns.Ln[ABS[-S/2.0+p2]]/3.0 ];
          p4: REAL_ 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
            Icount_ Icount+1; 
            x4[Icount]_ -a4[I]/2.0+SqRt[delta4]; 
            Icount_ Icount+1; 
            x4[Icount]_ -a4[I]/2.0-SqRt[delta4];
          END; 
        ENDLOOP;
        
        FOR J:CARDINAL IN [1..Icount] DO 
          IF x4[J]<=0.0 THEN LOOP; 
          
          poly4[J]_ x4[J]*x4[J]*x4[J]*x4[J]+a1*x4[J]*x4[J]*x4[J]; 
          poly4[J]_ poly4[J]+b1*x4[J]*x4[J]+c1*x4[J]+d1;
          
          IF ABS[poly4[J]]>1.0E-6 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; -- Block 545
        
      -- 600, EVALUATE EFFECTIVE CHANNEL LENGTH AND ITS DERIVATIVES
      IF Vds # 0.0 THEN BEGIN
        lFactor, ArgV, xdv, lv, ls, SArgV: REAL;
         
        gammaD_ gamaSD; 
        BSArg_ IF Vbs <= vdsat THEN SqRt[vdsat-Vbs+phi];
               ELSE SqrtPhi/(1.0+0.5*(Vbs-vdsat)/phi);
          
        bodyS_ BSArg*BSArg*BSArg-SArg3; 
        IF vmax<=0.0 THEN BEGIN
          IF nsub # 0.0 AND lambda <= 0.0 THEN BEGIN
            ArgV_ (Vds-vdsat)/4.0;
            SArgV_ SqRt[1.0+ArgV*ArgV];
            Arg_ SqRt[ArgV+SArgV]; 
            lFactor_ xd/(l*Vds);
            lambda_ lFactor*Arg;
	END;
        END ELSE BEGIN -- 603 
          ArgV_ (VgsX-Vbin)/eta-vdsat;
          xdv_ xd/SqRt[neff]; 
          lv_ vmax*xdv/(2.0*UEffective); 
          IF nsub # 0.0 AND lambda <= 0.0 THEN BEGIN
            ArgV_ MAX[Vds-vdsat, 0.0]; 
            ls_ SqRt[lv*lv+ArgV]; 
            lFactor_ xdv/(l*Vds); 
            lambda_ lFactor*(ls-lv);
	END;
        END;
      END; -- block 600
      
      -- 620, LIMIT CHANNEL SHORTENING AT PUNCH-THROUGH 
      BEGIN
        leffective, deltaL: REAL;
        xwb: REAL_ xd*SqrtPhib;
        ld: REAL_ l-xwb;
        clFactor_ 1.0-lambda*Vds; 
        leffective_ l*clFactor; 
        deltaL_ lambda*Vds*l;
        IF nsub=0.0 THEN xwb_ 0.25E-6; 
        IF leffective<xwb THEN BEGIN  
          leffective_ xwb/(1.0+(deltaL-ld)/xwb);
          clFactor_ leffective/l; 
        END;
      END;
  
      -- 700 EVALUATE EFFECTIVE beta (EFFECTIVE KP) 
      beta1_ beta*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 
        VdsOn, IdsOn: REAL;
        
        -- SUBTHRESHOLD REGION, IF vdsat<=0.0 THEN Id_ 0.0
        IF  vdsat > 0.0 THEN BEGIN -- 830
          VdsOn_ MIN[vdsat,Vds];
          IF Vds>vdsat THEN body_ bodyS;
          -- 850 
          IdsOn_ beta1*((von-Vbin-eta*VdsOn*0.5)*VdsOn-gammaD*body/1.5);
          Id_ sign*IdsOn*RealFns.Exp[ArgG*VgsT];
        END;
        
      END ELSE 
        Id_ IF Vds <= vdsat THEN -- 900 LINEAR REGION
              sign*beta1*((Vgs-Vbin-eta*Vds/2.0)*Vds-gammaD*body/1.5) 
            ELSE -- 1000 SATURATION REGION
              sign*beta1*((Vgs-Vbin-eta*vdsat/2.0)*vdsat-gammaD*bodyS/1.5);
    EXITS
      done=> NULL;
    END;
    
    -- capacitances, and results
    VgdT: REAL_ Vg - Vd - von;

    -- Cgb:
    Cgb_ IF VgsT <= 0 THEN CgbM ELSE 0;
  		     
    -- Cgs and Cgd:

    t1_ VgsT+VgdT;
    t1_ t1*t1; -- (Vgs + Vgd - 2von)^2
    Cgs_ IF VgsT <= 0.0 THEN CovGS
       ELSE IF VgdT <= 0.0 THEN Cgb23rds + CovGS
       ELSE CovGS + Cgb23rds*(1.0 - VgdT*VgdT/t1);
    Cgd_ IF VgdT <= 0.0 THEN CovDS
       ELSE CovDS + Cgb23rds*(1.0 - VgsT*VgsT/t1);

    -- Cbs and Cbd

    IF Vbs >= FcByPb THEN -- 520
      Cbs_ CbsFwd1 + Vbs*CbsFwd2
    ELSE BEGIN
      t1_ RealFns.Ln[1.0-Vbs/PhiB];
      Cbs_ CbsBottom*RealFns.Exp[-Mj*t1] + CbsSW*RealFns.Exp[-Mjsw*t1]
    END;
    
    IF Vbd >= FcByPb THEN -- 520
      Cbd_ CbdFwd1 + Vbd*CbdFwd2
    ELSE BEGIN
      t1_ RealFns.Ln[1.0-Vbd/PhiB];
      Cbd_ CbdBottom*RealFns.Exp[-Mj*t1] + CbdSW*RealFns.Exp[-Mjsw*t1]
    END;
    
    -- results order: Id, Cgb, Cgs, Cgd, Cbs, Cbd, Ibs, Ibd.
    results_ IF reverse THEN [-Id, Cgb, Cgd, Cgs, Cbd, Cbs, Ibd, Ibs]
                        ELSE [ Id, Cgb, Cgs, Cgd, Cbs, Cbd, Ibs, Ibd];
    END;
  END; -- SpiceLevel2  
  
  initSqRt[];
  
  spModelDefs.EnterModels["SpiceLevel2", SpiceLevel2, 4, 37, 8];

END.