module thinfilm_mos(nd, ngf, ns, ngb);
 inout nd,ngf,ns,ngb;
 electrical nd, ngf, ns, ngb;
 parameter real l=10u from (0:inf),
                w=10u from (0:inf);
 parameter real toxf   = 20n from (0:inf),      // m
                toxb   = 0.5u from (0:inf),     // m
                nsub   = 1e14 from (0:inf),     // 1/cm^3
                ngate  = 1e19 from (0:inf),     // 1/cm^3
                nbody  = 5e15 from (0:inf),     // 1/cm^3
                tb     = 0.1u from (0:inf),     // m
                u0     = 700 from (0:inf),      // cm-2/(v*S)
                lambda = 0.05;                  //
 parameter integer type=1 from [-1:1] exclude 0;

 real phi,      // body potential.
      vfbf,     // flat-band voltage - front channel.
      vfbb,     // flat-band voltage - back channel.
      vtfa,     // threshold voltage - back channel accumulated.
      vgba,     // vgb for accumulation at back surface.
      vgbi,     // vgb for inversion at back surface.
      vtff,     // threshold voltage.
      wkf,      // work-function, front-channel.
      wkb,      // work-function, back-channel.
      alpha,    // capacitance ratio.
      cob,      // capacitance back-gate to body.
      cof,      // capacitance front-gate to body.
      cb,       // body intrinsic capacitance.
      cbb,      // series body / back-gate capacitance.
      cfb,      // series front-gate / body capacitance.
      cfbb,     // series front-gate / body / back-gate cap.
      qb,       // fixed depleted body charge.
      kp,       // K-prime.
      qgf,      // front-gate charge.
      qgb,      // back-gate charge.
      qn,       // channel charge.
      qd,       // drain component of channel charge.
      qs;       // source component of channel charge.

 analog
 begin
   @(initial_step)
   begin : my
     real Vt,ni,epsox,epsil;

     ni = 9.6e9; // 1/cm^3
     epsox = 3.9*`P_EPS0;
     epsil = 11.7*`P_EPS0;
     Vt = `P_K*$temperature()/`P_Q;
     phi = 2*Vt*ln(nbody/ni);
     wkf = Vt*ln(ngate/ni) - phi/2;
     wkb = Vt*ln(nsub/ni) - phi/2;
     vfbf = wkf; // front-channel fixed charge assumed zero.
     vfbb = wkb; // back-channel fixed charge assumed zero.
     qb = `P_Q*nbody*1e6*tb;
     cob = epsox/toxb;
     cof = epsox/toxf;
     cb = epsil/tb;
     cbb = cob*cb/(cob + cb);
     cfb = cof*cb/(cof + cb);
     cfbb = cfb*cob/(cfb + cob);
     alpha = cbb/cof;
     vtfa = vfbf + (1 + cb/cof)*phi - qb/(2*cof);
     vgba = type*vfbb - phi*cb/cob - qb/(2*cob);
     vgbi = type*vfbb + phi - qb/(2*cob);
     kp = w*u0*1e-4*cof/l;
     back_surf = 0;
   end

   vgfs = type*V(ngf,ns);
   vds  = type*V(nd,ns);
   vgbs = type*V(ngb,ns);

// calc. threshold and saturation voltages.
   vtff = vtfa - (vgbs - vgba)*cbb/cof;
   vdsat = (vgfs - vtff)/(1 + alpha);

// drain current calculations.
   if (vgfs < vtff)
   begin

     // front-channel in accumulation / cutoff region(s).
     id = 0;
     qn = 0;
     qd = 0;
     qs = 0;
     qgf = w*l*cfbb*(vgfs - wkf - qb/(2*cbb) -
                     (vgbs - vfbb + qb/(2*cob)));
     qgb = - (qgf + w*l*qb);
   end
   else
     if (vds < vdsat)
     begin : LinearRegion

       // front-channel in linear region.
       real tmp1;
       id = kp*((vgfs - vtff)*vds - 0.5*(1 + cbb/cof)*vds*vds);
       id = id*(1 + lambda*vds);
       tmp1 = (1 + alpha)*vds;
       qn = -w*l*cof*(vgfs - vtff - tmp1/2 +
                           tmp1*tmp1/(12*(vgfs - vtff - tmp1/2)));
       qd = 0.4*qn;
       qs = 0.6*qn;
       qgf = w*l*cof*(vgfs - wkf - phi - vds/2 +
                           tmp1*vds/(12*(vgfs - vtff - tmp1/2)));
       qgb = - (qgf + qn + w*l*qb);
     end
     else
     begin

       // front-channel in saturation.
       id = 0.5*kp*(pow((vgfs - vtff), 2))/(1 + cbb/cof);
       id = id*(1 + lambda*vds);
       qn = -w*l*cof*(2/3)*(vgfs - vtff);
       qd = 0.4*qn;
       qs = 0.6*qn;
       qgf = w*l*cof*(vgfs - wkf - phi -
                           ((vgfs - vtff)/(3*(1 + alpha))));
       qgb = - (qgf + qn + w*l*qb);
     end

   // intrinsic device.
   I(nd,ns)   <+ type*id;
   I(nd,ngb)  <+ type*ddt(qd);
   I(ns,ngb)  <+ type*ddt(qs);
   I(ngf,ngb) <+ type*ddt(qgf);

`ifdef CHECK_BACK_SURFACE
   // check back-surface contraints. save the state in back_surf
   // at final, use back_surf to print out any violations.
   // if back_surf != 0, already have a warning, so continue.

   if (!back_surf)
     if(vgbs > vgbi)    back_surf = 1;
     else
       if (vgbs < vgba) back_surf = 2;
`endif


`ifdef CHECK_BACK_SURFACE
   @(final_step)
     case (back_surf)
       1 : $strobe("Back-surface went into inversion during analysis.\n");
       2 : $strobe("Back-surface went into accumulation during analysis.\n");
     endcase
`endif

 end
endmodule