Logo Search packages:      
Sourcecode: gnucap version File versions  Download package

d_mos7.cc

/* $Id: d_mos7.model,v 24.4 2003/04/06 10:36:01 al Exp $ -*- C++ -*-
 * Copyright (C) 2001 Albert Davis
 * Author: Albert Davis <aldavis@ieee.org>
 *
 * This file is part of "Gnucap", the Gnu Circuit Analysis Package
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2, or (at your option)
 * any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
 * 02111-1307, USA.
 *------------------------------------------------------------------
 * Berkeley BSIM3v3.1 model
 * Derived from Spice3f4,Copyright 1990 Regents of the University of California
 * Author: 1991 JianHui Huang and Min-Chie Jeng.
 * Recoded for Gnucap model compiler, Al Davis, 2000
 */
/* This file is automatically generated. DO NOT EDIT */

#include "l_compar.h"
#include "l_denoise.h"
#include "ap.h"
#include "d_mos7.h"
/*--------------------------------------------------------------------------*/
const double NA(NOT_INPUT);
const double INF(BIGBIG);
/*--------------------------------------------------------------------------*/
int MODEL_MOS7::_count = 0;
/*--------------------------------------------------------------------------*/
SDP_MOS7::SDP_MOS7(const COMMON_COMPONENT* cc)
  :SDP_MOS_BASE(cc)
{
  assert(cc);
  const COMMON_MOS* c = prechecked_cast<const COMMON_MOS*>(cc);
  assert(c);
  const MODEL_MOS7* m = prechecked_cast<const MODEL_MOS7*>(c->model());
  assert(m);

      {
      double T0 = pow(c->l_in, m->Lln);
      double T1 = pow(c->w_in, m->Lwn);
      double tmp1 = m->Ll / T0 + m->Lw / T1 + m->Lwl / (T0 * T1);
      dl = m->Lint + tmp1;
      dlc = m->dlc + tmp1;
      }
      {
      double T2 = pow(c->l_in, m->Wln);
      double T3 = pow(c->w_in, m->Wwn);
      double tmp2 = m->Wl / T2 + m->Ww / T3 + m->Wwl / (T2 * T3);
      dw = m->Wint + tmp2;
      dwc = m->dwc + tmp2;
      }

      leff = c->l_in - 2.0 * dl;
      weff = c->w_in - 2.0 * dw;
      leffCV = c->l_in - 2.0 * dlc;
      weffCV = c->w_in - 2.0 * dwc;
      cgate = m->cox * w_eff * l_eff; /* BUG:: not adjusted values?? */
      double L = leff;
      double W = weff;
      if (m->binUnit == 1) {
      L /= MICRON2METER;
      W /= MICRON2METER;
      }
  cdsc = m->cdsc(L, W);
  cdscb = m->cdscb(L, W);
  cdscd = m->cdscd(L, W);
  cit = m->cit(L, W);
  nfactor = m->nfactor(L, W);
  xj = m->xj(L, W);
  vsat = m->vsat(L, W);
  at = m->at(L, W);
  a0 = m->a0(L, W);
  ags = m->ags(L, W);
  a1 = m->a1(L, W);
  a2 = m->a2(L, W);
  keta = m->keta(L, W);
  nsub = m->nsub(L, W);
  npeak = m->npeak(L, W);
  ngate = m->ngate(L, W);
  gamma1 = m->gamma1(L, W);
  gamma2 = m->gamma2(L, W);
  vbx = m->vbx(L, W);
  vbm = m->vbm(L, W);
  xt = m->xt(L, W);
  k1 = m->k1(L, W);
  kt1 = m->kt1(L, W);
  kt1l = m->kt1l(L, W);
  kt2 = m->kt2(L, W);
  k2 = m->k2(L, W);
  k3 = m->k3(L, W);
  k3b = m->k3b(L, W);
  w0 = m->w0(L, W);
  nlx = m->nlx(L, W);
  dvt0 = m->dvt0(L, W);
  dvt1 = m->dvt1(L, W);
  dvt2 = m->dvt2(L, W);
  dvt0w = m->dvt0w(L, W);
  dvt1w = m->dvt1w(L, W);
  dvt2w = m->dvt2w(L, W);
  drout = m->drout(L, W);
  dsub = ((m->dsub.nom() == NA)
    ? drout
    : m->dsub(L, W));
  vth0 = ((m->vth0.nom() == NA)
    ? NA
    : m->vth0(L, W));
  ua1 = m->ua1(L, W);
  ua = m->ua(L, W);
  ub1 = m->ub1(L, W);
  ub = m->ub(L, W);
  uc1 = ((m->uc1.nom() == NA)
    ? ((m->mobMod==3) ? -0.056 : -0.056e-9)
    : m->uc1(L, W));
  uc = ((m->uc.nom() == NA)
    ? ((m->mobMod==3) ? -0.0465 : -0.0465e-9)
    : m->uc(L, W));
  u0 = ((m->u0.nom() == NA)
    ? ((m->polarity == pN) ? 0.067 : 0.025)
    : m->u0(L, W));
  ute = m->ute(L, W);
  voff = m->voff(L, W);
  delta = m->delta(L, W);
  rdsw = m->rdsw(L, W);
  prwg = m->prwg(L, W);
  prwb = m->prwb(L, W);
  prt = m->prt(L, W);
  eta0 = m->eta0(L, W);
  etab = m->etab(L, W);
  pclm = m->pclm(L, W);
  pdibl1 = m->pdibl1(L, W);
  pdibl2 = m->pdibl2(L, W);
  pdiblb = m->pdiblb(L, W);
  pscbe1 = m->pscbe1(L, W);
  pscbe2 = m->pscbe2(L, W);
  pvag = m->pvag(L, W);
  wr = m->wr(L, W);
  dwg = m->dwg(L, W);
  dwb = m->dwb(L, W);
  b0 = m->b0(L, W);
  b1 = m->b1(L, W);
  alpha0 = m->alpha0(L, W);
  beta0 = m->beta0(L, W);
  elm = m->elm(L, W);
  vfbcv = m->vfbcv(L, W);
  cgsl = m->cgsl(L, W);
  cgdl = m->cgdl(L, W);
  ckappa = m->ckappa(L, W);
  cf = ((m->cf.nom() == NA)
    ? 2.0 * E_OX / kPI * log(1.0 + 0.4e-6 / m->tox)
    : m->cf(L, W));
  clc = m->clc(L, W);
  cle = m->cle(L, W);
  abulkCVfactor = 1.0 + pow((clc / leff), cle);
  cgso = (m->cgso + cf) * weffCV;
  cgdo = (m->cgdo + cf) * weffCV;
  cgbo = m->cgbo * leffCV;
  litl = sqrt(3.0 * xj * m->tox);

      if (u0 > 1.0) {
      u0 /= 1.0e4;
      }
      if (m->npeak.nom() == NA) {
      {if (m->gamma1.nom() != NA) {
        double T0 = gamma1 * m->cox;
        npeak = 3.021E22 * T0 * T0;
      }else{
        npeak = 1.7e17;
      }}
      }
      {if (m->k1.nom() != NA && m->k2.nom() != NA) {
      if (m->k1.nom() == NA) {
        k1 = 0.53;
      }
      if (m->k2.nom() == NA) {
        k2 = -0.0186;
      }
      }else{
      vbm = -std::abs(vbm);
      if (m->gamma1.nom() == NA) {
        gamma1 = 5.753e-12 * sqrt(npeak) / m->cox;
      }
      if (m->gamma2.nom() == NA) {
        gamma2 = 5.753e-12 * sqrt(nsub) / m->cox;
      }
      }}
}
/*--------------------------------------------------------------------------*/
TDP_MOS7::TDP_MOS7(const DEV_MOS* d)
  :TDP_MOS_BASE(d)
{
  assert(d);
  const COMMON_MOS* c = prechecked_cast<const COMMON_MOS*>(d->common());
  assert(c);
  const SDP_MOS7* s = prechecked_cast<const SDP_MOS7*>(c->sdp());
  assert(s);
  const MODEL_MOS7* m = prechecked_cast<const MODEL_MOS7*>(c->model());
  assert(m);

      temp = SIM::temp;
      double kt = temp * K;
      double egap = 1.16 - 7.02e-4 * temp * temp / (temp + 1108.0);
  tempratio = temp / m->_tnom;
  tempratio_1 = tempratio - 1;
  vtm = kt / Q;

      double jctTempSatCurDensity;
      double jctSidewallTempSatCurDensity;
      {if (temp != m->_tnom) {
      double T0 = m->egap / m->vt_at_tnom - egap / vtm + m->jctTempExponent
        * log(temp / m->_tnom);
      double T1 = exp(T0 / m->jctEmissionCoeff);
      jctTempSatCurDensity = m->js * T1;
      jctSidewallTempSatCurDensity = m->jctSidewallSatCurDensity * T1;
      }else{
      jctTempSatCurDensity = m->js;
      jctSidewallTempSatCurDensity = m->jctSidewallSatCurDensity;
      }}
      if (jctTempSatCurDensity < 0.0) {
      jctTempSatCurDensity = 0.0;
      }
      if (jctSidewallTempSatCurDensity < 0.0) {
      jctSidewallTempSatCurDensity = 0.0;
      }
      {
      double T0 = (tempratio - 1.0);
      ua = s->ua + s->ua1 * T0;
      ub = s->ub + s->ub1 * T0;
      uc = s->uc + s->uc1 * T0;
      u0temp = s->u0 * pow(tempratio, s->ute); 
      vsattemp = s->vsat - s->at * T0;
      rds0 = (s->rdsw + s->prt * T0) / pow(s->weff * 1E6, s->wr);
      }
      phi = 2.0 * m->vt_at_tnom * log(s->npeak / m->ni);
      sqrtPhi = sqrt(phi);
      phis3 = sqrtPhi * phi;
      Xdep0 = sqrt(2.0 * E_SI / (Q * s->npeak * 1.0e6)) * sqrtPhi;
      vbi = m->vt_at_tnom * log(1.0e20 * s->npeak / (m->ni * m->ni));
      cdep0 = sqrt(Q * E_SI * s->npeak * 1.0e6 / 2.0 / phi);

      {if (m->k1.nom() != NA && m->k2.nom() != NA) {
      k2 = s->k2;
      k1 = s->k1;
      }else{
      double vbx = (m->vbx.nom() == NA)
        ? -std::abs(phi - 7.7348e-4 * s->npeak * s->xt * s->xt)
        : -std::abs(s->vbx);
      double T0 = s->gamma1 - s->gamma2;
      double T1 = sqrt(phi - vbx) - sqrtPhi;
      double T2 = sqrt(phi * (phi - s->vbm)) - phi;
      k2 = T0 * T1 / (2.0 * T2 + s->vbm);
      k1 = s->gamma2 - 2.0 * k2 * sqrt(phi - s->vbm);
      }}
      {if (k2 < 0.) {
      double T0 = 0.5 * k1 / k2;
      vbsc = to_range(-30.0, (0.9 * (phi - T0 * T0)), -3.0);
      }else{
      vbsc = -30.0;
      }}
      vbsc = std::min(vbsc, s->vbm);
      {if (s->vth0 == NA) {
      vfb = -1.0;
      vth0 = m->polarity * (vfb + phi + k1 * sqrtPhi);
      }else{
      vth0 = s->vth0;
      vfb = m->polarity * vth0 - phi - k1 * sqrtPhi;
      }}
      {
      double T1 = sqrt(E_SI / E_OX * m->tox * Xdep0);
      double T0 = exp(-0.5 * s->dsub * s->leff / T1);
      theta0vb0 = (T0 + 2.0 * T0 * T0);

      T0 = exp(-0.5 * s->drout * s->leff / T1);
      double T2 = (T0 + 2.0 * T0 * T0);
      thetaRout = s->pdibl1 * T2 + s->pdibl2;
      }
}
/*--------------------------------------------------------------------------*/
MODEL_MOS7::MODEL_MOS7()
  :MODEL_MOS_BASE(),
   cdsc(2.4e-4),
   cdscb(0.0),
   cdscd(0.0),
   cit(0.0),
   nfactor(1),
   xj(.15e-6),
   vsat(8.0e4),
   at(3.3e4),
   a0(1.0),
   ags(0.0),
   a1(0.0),
   a2(1.0),
   keta(-0.047),
   nsub(6.0e16),
   npeak(NA),
   ngate(0.0),
   gamma1(NA),
   gamma2(NA),
   vbx(NA),
   vbm(-3.0),
   xt(1.55e-7),
   k1(NA),
   kt1(-0.11),
   kt1l(0.0),
   kt2(0.022),
   k2(NA),
   k3(80.0),
   k3b(0.0),
   w0(2.5e-6),
   nlx(1.74e-7),
   dvt0(2.2),
   dvt1(0.53),
   dvt2(-0.032),
   dvt0w(0.0),
   dvt1w(5.3e6),
   dvt2w(-0.032),
   drout(0.56),
   dsub(NA),
   vth0(NA),
   ua1(4.31e-9),
   ua(2.25e-9),
   ub1(-7.61e-18),
   ub(5.87e-19),
   uc1(NA),
   uc(NA),
   u0(NA),
   ute(-1.5),
   voff(-0.08),
   delta(0.01),
   rdsw(0.0),
   prwg(0.0),
   prwb(0.0),
   prt(0.0),
   eta0(0.08),
   etab(-0.07),
   pclm(1.3),
   pdibl1(.39),
   pdibl2(0.0086),
   pdiblb(0.0),
   pscbe1(4.24e8),
   pscbe2(1.0e-5),
   pvag(0.0),
   wr(1.0),
   dwg(0.0),
   dwb(0.0),
   b0(0.0),
   b1(0.0),
   alpha0(0.0),
   beta0(30.0),
   elm(5.0),
   vfbcv(-1.0),
   cgsl(0.0),
   cgdl(0.0),
   ckappa(0.6),
   cf(NA),
   clc(0.1e-6),
   cle(0.6),
   capMod(2),
   nqsMod(0),
   mobMod(1),
   noiMod(1),
   paramChk(0),
   binUnit(1),
   version(3.1),
   tox(150.0e-10),
   xpart(0.0),
   jctSidewallSatCurDensity(0.0),
   mjswg(NA),
   pbswg(NA),
   unitLengthGateSidewallJctCap(NA),
   jctEmissionCoeff(1.0),
   jctTempExponent(3.0),
   Lint(0.0),
   Ll(0.0),
   Lln(1.0),
   Lw(0.0),
   Lwn(1.0),
   Lwl(0.0),
   Wint(0.0),
   Wl(0.0),
   Wln(1.0),
   Ww(0.0),
   Wwn(1.0),
   Wwl(0.0),
   dwc(NA),
   dlc(NA),
   noia(NA),
   noib(NA),
   noic(NA),
   em(4.1e7),
   ef(1.0),
   cox(NA),
   factor1(NA),
   vt_at_tnom(NA),
   ni(NA)
{
  ++_count;
  mjsw = NA;
  pb = NA;
  pbsw = NA;
  cjo = 5.0E-4;
  cgdo = NA;
  cgso = NA;
  cgbo = NA;
  mos_level = LEVEL;
}
/*--------------------------------------------------------------------------*/
bool MODEL_MOS7::parse_front(CS& cmd)
{
  return MODEL_MOS_BASE::parse_front(cmd);
}
/*--------------------------------------------------------------------------*/
bool MODEL_MOS7::parse_params(CS& cmd)
{
  return ONE_OF
    || get(cmd, "DIODElevel", &mos_level)
    || get(cmd, "CDSC", &cdsc)
    || get(cmd, "CDSCB", &cdscb)
    || get(cmd, "CDSCD", &cdscd)
    || get(cmd, "CIT", &cit)
    || get(cmd, "NFACTOR", &nfactor)
    || get(cmd, "XJ", &xj)
    || get(cmd, "VSAT", &vsat)
    || get(cmd, "AT", &at)
    || get(cmd, "A0", &a0)
    || get(cmd, "AGS", &ags)
    || get(cmd, "A1", &a1)
    || get(cmd, "A2", &a2)
    || get(cmd, "KETA", &keta)
    || get(cmd, "NSUB", &nsub)
    || get(cmd, "NCH", &npeak)
    || get(cmd, "NGATE", &ngate)
    || get(cmd, "GAMMA1", &gamma1)
    || get(cmd, "GAMMA2", &gamma2)
    || get(cmd, "VBX", &vbx)
    || get(cmd, "VBM", &vbm)
    || get(cmd, "XT", &xt)
    || get(cmd, "K1", &k1)
    || get(cmd, "KT1", &kt1)
    || get(cmd, "KT1L", &kt1l)
    || get(cmd, "KT2", &kt2)
    || get(cmd, "K2", &k2)
    || get(cmd, "K3", &k3)
    || get(cmd, "K3B", &k3b)
    || get(cmd, "W0", &w0)
    || get(cmd, "NLX", &nlx)
    || get(cmd, "DVT0", &dvt0)
    || get(cmd, "DVT1", &dvt1)
    || get(cmd, "DVT2", &dvt2)
    || get(cmd, "DVT0W", &dvt0w)
    || get(cmd, "DVT1W", &dvt1w)
    || get(cmd, "DVT2W", &dvt2w)
    || get(cmd, "DROUT", &drout)
    || get(cmd, "DSUB", &dsub)
    || get(cmd, "VTH0", &vth0)
    || get(cmd, "UA1", &ua1)
    || get(cmd, "UA", &ua)
    || get(cmd, "UB1", &ub1)
    || get(cmd, "UB", &ub)
    || get(cmd, "UC1", &uc1)
    || get(cmd, "UC", &uc)
    || get(cmd, "U0", &u0)
    || get(cmd, "UTE", &ute)
    || get(cmd, "VOFF", &voff)
    || get(cmd, "DELTA", &delta)
    || get(cmd, "RDSW", &rdsw)
    || get(cmd, "PRWG", &prwg)
    || get(cmd, "PRWB", &prwb)
    || get(cmd, "PRT", &prt)
    || get(cmd, "ETA0", &eta0)
    || get(cmd, "ETAB", &etab)
    || get(cmd, "PCLM", &pclm)
    || get(cmd, "PDIBLC1", &pdibl1)
    || get(cmd, "PDIBLC2", &pdibl2)
    || get(cmd, "PDIBLCB", &pdiblb)
    || get(cmd, "PSCBE1", &pscbe1)
    || get(cmd, "PSCBE2", &pscbe2)
    || get(cmd, "PVAG", &pvag)
    || get(cmd, "WR", &wr)
    || get(cmd, "DWG", &dwg)
    || get(cmd, "DWB", &dwb)
    || get(cmd, "B0", &b0)
    || get(cmd, "B1", &b1)
    || get(cmd, "ALPHA0", &alpha0)
    || get(cmd, "BETA0", &beta0)
    || get(cmd, "ELM", &elm)
    || get(cmd, "VFBCV", &vfbcv)
    || get(cmd, "CGSL", &cgsl)
    || get(cmd, "CGDL", &cgdl)
    || get(cmd, "CKAPPA", &ckappa)
    || get(cmd, "CF", &cf)
    || get(cmd, "CLC", &clc)
    || get(cmd, "CLE", &cle)
    || get(cmd, "CAPMOD", &capMod)
    || get(cmd, "NQSMOD", &nqsMod)
    || get(cmd, "MOBMOD", &mobMod)
    || get(cmd, "NOIMOD", &noiMod)
    || get(cmd, "PARAMCHK", &paramChk)
    || get(cmd, "BINUNIT", &binUnit)
    || get(cmd, "VERSION", &version)
    || get(cmd, "TOX", &tox)
    || get(cmd, "XPART", &xpart)
    || get(cmd, "JSW", &jctSidewallSatCurDensity)
    || get(cmd, "MJSWG", &mjswg)
    || get(cmd, "PBSWG", &pbswg)
    || get(cmd, "CJSWG", &unitLengthGateSidewallJctCap)
    || get(cmd, "NJ", &jctEmissionCoeff)
    || get(cmd, "XTI", &jctTempExponent)
    || get(cmd, "LINT", &Lint)
    || get(cmd, "LL", &Ll)
    || get(cmd, "LLN", &Lln)
    || get(cmd, "LW", &Lw)
    || get(cmd, "LWN", &Lwn)
    || get(cmd, "LWL", &Lwl)
    || get(cmd, "WINT", &Wint)
    || get(cmd, "WL", &Wl)
    || get(cmd, "WLN", &Wln)
    || get(cmd, "WW", &Ww)
    || get(cmd, "WWN", &Wwn)
    || get(cmd, "WWL", &Wwl)
    || get(cmd, "DWC", &dwc)
    || get(cmd, "DLC", &dlc)
    || get(cmd, "NOIA", &noia)
    || get(cmd, "NOIB", &noib)
    || get(cmd, "NOIC", &noic)
    || get(cmd, "EM", &em)
    || get(cmd, "EF", &ef)
    || MODEL_MOS_BASE::parse_params(cmd)
    ;
}
/*--------------------------------------------------------------------------*/
void MODEL_MOS7::parse_finish()
{
  MODEL_MOS_BASE::parse_finish();

      tox = std::max(tox, 1e-20);
      cox = 3.453133e-11 / tox;
  if (mjsw == NA) {
    mjsw = .33;
  }
  if (pb == NA) {
    pb = 1.0;
  }
  pb = std::max(pb, 0.1);
  if (pbsw == NA) {
    pbsw = pb;
  }
  pbsw = std::max(pbsw, 0.1);
  if (cgdo == NA) {
    cgdo = (((dlc != NA) && (dlc > 0.0))
        ? dlc * cox - cgdl.nom()
        : 0.6 * xj.nom() * cox);
  }
  if (cgso == NA) {
    cgso = (((dlc != NA) && (dlc > 0.0))
        ? dlc * cox - cgsl.nom()
        : 0.6 * xj.nom() * cox);
  }
  if (cgbo == NA) {
    cgbo = ((dwc != NA)
        ? 2.0 * dwc  * cox
        : 2.0 * Wint * cox);
  }
  cmodel = ((!cmodel)?1:cmodel);
  needs_isub = ((alpha0.nom()!=NA)&&(beta0.nom()!=NA));
  if (mjswg == NA) {
    mjswg = mjsw;
  }
  if (pbswg == NA) {
    pbswg = pbsw;
  }
  pbswg = std::max(pbswg, 0.1);
  if (unitLengthGateSidewallJctCap == NA) {
    unitLengthGateSidewallJctCap = cjsw;
  }
  if (dwc == NA) {
    dwc = Wint;
  }
  if (dlc == NA) {
    dlc = Lint;
  }
  if (noia == NA) {
    noia = (polarity==pN) ? 1e20 : 9.9e18;
  }
  if (noib == NA) {
    noib = (polarity==pN) ? 5e4 : 2.4e3;
  }
  if (noic == NA) {
    noic = (polarity==pN) ?-1.4e-12 :1.4e-12;
  }
  factor1 = sqrt(tox * E_SI / E_OX);
  vt_at_tnom = _tnom * K/Q;
  ni = (1.45e10 * (_tnom / 300.15) 
      * sqrt(_tnom / 300.15) * exp(21.5565981 - egap / (2.0 * vt_at_tnom)));

      if (npeak.nom() > 1.0e20) {
      npeak.set_nom(npeak.nom() * 1.0e-6);
      }
      if (ngate.nom() > 1.0e23) {
      ngate.set_nom(ngate.nom() * 1.0e-6);
      }
}
/*--------------------------------------------------------------------------*/
SDP_CARD* MODEL_MOS7::new_sdp(const COMMON_COMPONENT* c)const
{
  assert(c);
  {if (dynamic_cast<const COMMON_MOS*>(c)) {
    return new SDP_MOS7(c);
  }else{
    return MODEL_MOS_BASE::new_sdp(c);
  }}
}
/*--------------------------------------------------------------------------*/
void MODEL_MOS7::print_front(OMSTREAM& o)const
{
  MODEL_MOS_BASE::print_front(o);
}
/*--------------------------------------------------------------------------*/
void MODEL_MOS7::print_params(OMSTREAM& o)const
{
  o << "level=7";
  MODEL_MOS_BASE::print_params(o);
  if (mos_level != LEVEL)
    o << "  diodelevel=" << mos_level;
  cdsc.print(o, "cdsc");
  cdscb.print(o, "cdscb");
  cdscd.print(o, "cdscd");
  cit.print(o, "cit");
  nfactor.print(o, "nfactor");
  xj.print(o, "xj");
  vsat.print(o, "vsat");
  at.print(o, "at");
  a0.print(o, "a0");
  ags.print(o, "ags");
  a1.print(o, "a1");
  a2.print(o, "a2");
  keta.print(o, "keta");
  nsub.print(o, "nsub");
  npeak.print(o, "nch");
  ngate.print(o, "ngate");
  gamma1.print(o, "gamma1");
  gamma2.print(o, "gamma2");
  vbx.print(o, "vbx");
  vbm.print(o, "vbm");
  xt.print(o, "xt");
  k1.print(o, "k1");
  kt1.print(o, "kt1");
  kt1l.print(o, "kt1l");
  kt2.print(o, "kt2");
  k2.print(o, "k2");
  k3.print(o, "k3");
  k3b.print(o, "k3b");
  w0.print(o, "w0");
  nlx.print(o, "nlx");
  dvt0.print(o, "dvt0");
  dvt1.print(o, "dvt1");
  dvt2.print(o, "dvt2");
  dvt0w.print(o, "dvt0w");
  dvt1w.print(o, "dvt1w");
  dvt2w.print(o, "dvt2w");
  drout.print(o, "drout");
  dsub.print(o, "dsub");
  vth0.print(o, "vth0");
  ua1.print(o, "ua1");
  ua.print(o, "ua");
  ub1.print(o, "ub1");
  ub.print(o, "ub");
  uc1.print(o, "uc1");
  uc.print(o, "uc");
  u0.print(o, "u0");
  ute.print(o, "ute");
  voff.print(o, "voff");
  delta.print(o, "delta");
  rdsw.print(o, "rdsw");
  prwg.print(o, "prwg");
  prwb.print(o, "prwb");
  prt.print(o, "prt");
  eta0.print(o, "eta0");
  etab.print(o, "etab");
  pclm.print(o, "pclm");
  pdibl1.print(o, "pdiblc1");
  pdibl2.print(o, "pdiblc2");
  pdiblb.print(o, "pdiblcb");
  pscbe1.print(o, "pscbe1");
  pscbe2.print(o, "pscbe2");
  pvag.print(o, "pvag");
  wr.print(o, "wr");
  dwg.print(o, "dwg");
  dwb.print(o, "dwb");
  b0.print(o, "b0");
  b1.print(o, "b1");
  alpha0.print(o, "alpha0");
  beta0.print(o, "beta0");
  elm.print(o, "elm");
  vfbcv.print(o, "vfbcv");
  cgsl.print(o, "cgsl");
  cgdl.print(o, "cgdl");
  ckappa.print(o, "ckappa");
  cf.print(o, "cf");
  clc.print(o, "clc");
  cle.print(o, "cle");
  o << "  capmod=" << capMod;
  o << "  nqsmod=" << nqsMod;
  o << "  mobmod=" << mobMod;
  o << "  noimod=" << noiMod;
  o << "  paramchk=" << paramChk;
  o << "  binunit=" << binUnit;
  o << "  version=" << version;
  o << "  tox=" << tox;
  o << "  xpart=" << xpart;
  o << "  jsw=" << jctSidewallSatCurDensity;
  if (mjswg != NA)
    o << "  mjswg=" << mjswg;
  if (pbswg != NA)
    o << "  pbswg=" << pbswg;
  if (unitLengthGateSidewallJctCap != NA)
    o << "  cjswg=" << unitLengthGateSidewallJctCap;
  o << "  nj=" << jctEmissionCoeff;
  o << "  xti=" << jctTempExponent;
  o << "  lint=" << Lint;
  o << "  ll=" << Ll;
  o << "  lln=" << Lln;
  o << "  lw=" << Lw;
  o << "  lwn=" << Lwn;
  o << "  lwl=" << Lwl;
  o << "  wint=" << Wint;
  o << "  wl=" << Wl;
  o << "  wln=" << Wln;
  o << "  ww=" << Ww;
  o << "  wwn=" << Wwn;
  o << "  wwl=" << Wwl;
  if (dwc != NA)
    o << "  dwc=" << dwc;
  if (dlc != NA)
    o << "  dlc=" << dlc;
  if (noia != NA)
    o << "  noia=" << noia;
  if (noib != NA)
    o << "  noib=" << noib;
  if (noic != NA)
    o << "  noic=" << noic;
  o << "  em=" << em;
  o << "  ef=" << ef;
}
/*--------------------------------------------------------------------------*/
void MODEL_MOS7::print_calculated(OMSTREAM& o)const
{
  MODEL_MOS_BASE::print_calculated(o);
}
/*--------------------------------------------------------------------------*/
bool MODEL_MOS7::is_valid(const COMMON_COMPONENT* cc)const
{
  return MODEL_MOS_BASE::is_valid(cc);
}
/*--------------------------------------------------------------------------*/
void MODEL_MOS7::tr_eval(COMPONENT* brh)const
{
  DEV_MOS* d = prechecked_cast<DEV_MOS*>(brh);
  assert(d);
  const COMMON_MOS* c = prechecked_cast<const COMMON_MOS*>(d->common());
  assert(c);
  const SDP_MOS7* s = prechecked_cast<const SDP_MOS7*>(c->sdp());
  assert(s);
  const MODEL_MOS7* m = this;
  const TDP_MOS7 T(d);
  const TDP_MOS7* t = &T;

    trace3("", d->vds, d->vgs, d->vbs);
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    const double EXP_THRESHOLD = 34.0;
    const double MIN_EXP = 1.713908431e-15;
    const double MAX_EXP = 5.834617425e14;
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    d->reverse_if_needed();
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    double Vbseff, dVbseff_dVb;
    {
      double T0 = d->vbs - t->vbsc - 0.001;
      double T1 = sqrt(T0 * T0 - 0.004 * t->vbsc);
      trace3("", t->vbsc, T0, T1);
      Vbseff = t->vbsc + 0.5 * (T0 + T1);
      dVbseff_dVb = 0.5 * (1.0 + T0 / T1);
      trace2("raw", Vbseff, dVbseff_dVb);

      fixzero(&Vbseff, t->vbsc);
      if (Vbseff < d->vbs) {  // From Spice, to fix numeric problems
      untested();       // inadequately.  Above fixzero should do a 
      Vbseff = d->vbs;  // better job, but I left this in case.
      }
    }
    trace2("fixed", Vbseff, dVbseff_dVb);
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    double Phis, dPhis_dVb, sqrtPhis, dsqrtPhis_dVb;
    {if (Vbseff > 0.0) {
      untested();
      d->sbfwd = true;
      double T0 = t->phi / (t->phi + Vbseff);
      Phis = t->phi * T0;
      dPhis_dVb = -T0 * T0;
      sqrtPhis = t->phis3 / (t->phi + 0.5 * Vbseff);
      dsqrtPhis_dVb = -0.5 * sqrtPhis * sqrtPhis / t->phis3;
      trace0("bs-fwd-bias");
    }else{
      d->sbfwd = false;
      Phis = t->phi - Vbseff;
      dPhis_dVb = -1.0;
      sqrtPhis = sqrt(Phis);
      dsqrtPhis_dVb = -0.5 / sqrtPhis;
      trace0("bs-normal");
    }}
    trace4("", Phis, dPhis_dVb, sqrtPhis, dsqrtPhis_dVb);
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    double Xdep = t->Xdep0 * sqrtPhis / t->sqrtPhi;
    double dXdep_dVb = (t->Xdep0 / t->sqrtPhi) * dsqrtPhis_dVb;
    trace2("", Xdep, dXdep_dVb);
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    double Theta0, dTheta0_dVb;
    {
      double lt1, dlt1_dVb;
      {
      double T3 = sqrt(Xdep);
      double T0 = s->dvt2 * Vbseff;
      double T1, T2;
      {if (T0 >= - 0.5) {
        T1 = 1.0 + T0;
        T2 = s->dvt2;
      }else{
        untested();
        /* Added to avoid any discontinuity problems caused by dvt2 */ 
        double T4 = 1.0 / (3.0 + 8.0 * T0);
        T1 = (1.0 + 3.0 * T0) * T4; 
        T2 = s->dvt2 * T4 * T4;
      }}
      trace4("", T0, T1, T2, T3);
      lt1 = m->factor1 * T3 * T1;
      dlt1_dVb = m->factor1 * (0.5 / T3 * T1 * dXdep_dVb + T3 * T2);
      }
      trace2("", lt1, dlt1_dVb);

      double T0 = -0.5 * s->dvt1 * s->leff / lt1;
      {if (T0 > -EXP_THRESHOLD) {
      double T1 = exp(T0);
      Theta0 = T1 * (1.0 + 2.0 * T1);
      double dT1_dVb = -T0 / lt1 * T1 * dlt1_dVb;
      dTheta0_dVb = (1.0 + 4.0 * T1) * dT1_dVb;
      }else{
      double T1 = MIN_EXP;
      Theta0 = T1 * (1.0 + 2.0 * T1);
      dTheta0_dVb = 0.0;
      }}
    }
    trace2("", Theta0, dTheta0_dVb);
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    double dVth_dVb, dVth_dVd; // d->von
    {
      double V0 = t->vbi - t->phi;
      double T2, dT2_dVb;
      {
      double ltw, dltw_dVb;
      {
        double T3 = sqrt(Xdep);
        double T0 = s->dvt2w * Vbseff;
        double T1, T2;
        {if (T0 >= - 0.5) {
          T1 = 1.0 + T0;
          T2 = s->dvt2w;
        }else{
          untested();
          /* Added to avoid any discontinuity problems caused by dvt2w */ 
          double T4 = 1.0 / (3.0 + 8.0 * T0);
          T1 = (1.0 + 3.0 * T0) * T4; 
          T2 = s->dvt2w * T4 * T4;
        }}
        trace4("", T0, T1, T2, T3);
        ltw = m->factor1 * T3 * T1;
        dltw_dVb = m->factor1 * (0.5 / T3 * T1 * dXdep_dVb + T3 * T2);
      }
      trace2("", ltw, dltw_dVb);
      double T0 = -0.5 * s->dvt1w * s->weff * s->leff / ltw;
      {if (T0 > -EXP_THRESHOLD) {
        double T1 = exp(T0);
        T2 = T1 * (1.0 + 2.0 * T1);
        double dT1_dVb = -T0 / ltw * T1 * dltw_dVb;
        dT2_dVb = (1.0 + 4.0 * T1) * dT1_dVb;
      }else{
        double T1 = MIN_EXP;
        T2 = T1 * (1.0 + 2.0 * T1);
        dT2_dVb = 0.0;
      }}
      T0 = s->dvt0w * T2;
      T2 = T0 * V0;
      dT2_dVb = s->dvt0w * dT2_dVb * V0;
      }
      trace3("", V0, T2, dT2_dVb);
      double T0 = sqrt(1.0 + s->nlx / s->leff);
      double T1 = t->k1 * (T0 - 1.0) * t->sqrtPhi
      + (s->kt1 + s->kt1l / s->leff + s->kt2 * Vbseff) * t->tempratio_1;
      double tmp2 = m->tox * t->phi / (s->weff + s->w0);
      
      double T3 = s->eta0 + s->etab * Vbseff;
      trace4("", T0, T1, tmp2, T3);
      double T4;
      {if (T3 < 1.0e-4) {
      untested();
      /* avoid  discontinuity problems caused by etab */ 
      double T9 = 1.0 / (3.0 - 2.0e4 * T3);
      T3 = (2.0e-4 - T3) * T9;
      T4 = T9 * T9;
      trace3("", T9, T3, T4);
      }else{
      T4 = 1.0;
      trace1("", T4);
      }}
      double thetavth = s->dvt0 * Theta0;
      double Delt_vth = thetavth * V0;
      double dDelt_vth_dVb = s->dvt0 * dTheta0_dVb * V0;
      trace4("", thetavth, t->theta0vb0, Delt_vth, dDelt_vth_dVb);
      double dDIBL_Sft_dVd = T3 * t->theta0vb0;
      double DIBL_Sft = dDIBL_Sft_dVd * d->vds;
      trace2("", dDIBL_Sft_dVd, DIBL_Sft);
      
      trace4("", t->vth0, t->k1, sqrtPhis, t->sqrtPhi);
      trace4("", t->k2, Vbseff, Delt_vth, T2);
      trace4("", s->k3, s->k3b, Vbseff, tmp2);
      trace2("", T1, DIBL_Sft);
      double Vth = m->polarity * t->vth0 + t->k1 * (sqrtPhis - t->sqrtPhi) 
      - t->k2 * Vbseff - Delt_vth - T2 + (s->k3 + s->k3b * Vbseff) * tmp2
      + T1 - DIBL_Sft;
      d->von = Vth; 
      
      dVth_dVb = t->k1 * dsqrtPhis_dVb - t->k2 - dDelt_vth_dVb - dT2_dVb
      + s->k3b * tmp2 - s->etab * d->vds * t->theta0vb0 * T4
      + s->kt2 * t->tempratio_1;
      dVth_dVd = -dDIBL_Sft_dVd; 
      trace3("", Vth, dVth_dVb, dVth_dVd);
    }
    // end trace
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    /* Calculate n */
    double n, dn_dVb, dn_dVd;
    {
      double tmp2 = s->nfactor * E_SI / Xdep;
      double tmp3 = s->cdsc + s->cdscb * Vbseff + s->cdscd * d->vds;
      double tmp4 = (tmp2 + tmp3 * Theta0 + s->cit) / m->cox;
      {if (tmp4 >= -0.5) {
      n = 1.0 + tmp4;
      dn_dVb = (-tmp2 / Xdep * dXdep_dVb + tmp3 * dTheta0_dVb 
              + s->cdscb * Theta0) / m->cox;
      dn_dVd = s->cdscd * Theta0 / m->cox;
      }else{
      /* avoid  discontinuity problems caused by tmp4 */ 
      double T0 = 1.0 / (3.0 + 8.0 * tmp4);
      n = (1.0 + 3.0 * tmp4) * T0;
      T0 *= T0;
      dn_dVb = (-tmp2 / Xdep * dXdep_dVb + tmp3 * dTheta0_dVb
              + s->cdscb * Theta0) / m->cox * T0;
      dn_dVd = s->cdscd * Theta0 / m->cox * T0;
      }}
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    /* Poly Gate Si Depletion Effect */
    double Vgs_eff, dVgs_eff_dVg;
    {
      double T0 = t->vfb + t->phi;
      {if ((s->ngate > 1.e18) && (s->ngate < 1.e25) && (d->vgs > T0)) {
      /* added to avoid the problem caused by ngate */
      double T1 = 1.0e6 * Q * E_SI * s->ngate / (m->cox * m->cox);
      double T4 = sqrt(1.0 + 2.0 * (d->vgs - T0) / T1);
      double T2 = T1 * (T4 - 1.0);
      double T3 = 0.5 * T2 * T2 / T1; /* T3 = Vpoly */
      double T7 = 1.12 - T3 - 0.05;
      double T6 = sqrt(T7 * T7 + 0.224);
      double T5 = 1.12 - 0.5 * (T7 + T6);
      Vgs_eff = d->vgs - T5;
      dVgs_eff_dVg = 1.0 - (0.5 - 0.5 / T4) * (1.0 + T7 / T6); 
      }else{
      Vgs_eff = d->vgs;
      dVgs_eff_dVg = 1.0;
      }}
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    /* Effective Vgst (Vgsteff) Calculation */
    double /*Vgsteff,*/ dVgsteff_dVg, dVgsteff_dVd, dVgsteff_dVb, Vgst2Vtm;
    double VgstNVt, ExpVgst; // d->vgst
    {
      double Vgst = Vgs_eff - d->von;
      double T10 = 2.0 * n * t->vtm;
      VgstNVt = Vgst / T10;
      double ExpArg = (2.0 * s->voff - Vgst) / T10;

      /* MCJ: Very small Vgst */
      {if (VgstNVt > EXP_THRESHOLD) {
      d->vgst = Vgst;
      dVgsteff_dVg = dVgs_eff_dVg;
      dVgsteff_dVd = -dVth_dVd;
      dVgsteff_dVb = -dVth_dVb;
      ExpVgst = NOT_VALID;
      }else if (ExpArg > EXP_THRESHOLD) {
      double T0 = (Vgst - s->voff) / (n * t->vtm);
      ExpVgst = exp(T0);
      d->vgst = t->vtm * t->cdep0 / m->cox * ExpVgst;
      dVgsteff_dVg = d->vgst / (n * t->vtm);
      dVgsteff_dVd = -dVgsteff_dVg * (dVth_dVd + T0 * t->vtm * dn_dVd);
      dVgsteff_dVb = -dVgsteff_dVg * (dVth_dVb + T0 * t->vtm * dn_dVb);
      dVgsteff_dVg *= dVgs_eff_dVg;
      }else{
      ExpVgst = exp(VgstNVt);
      double T1 = T10 * log(1.0 + ExpVgst);
      double dT1_dVg = ExpVgst / (1.0 + ExpVgst);
      double dT1_dVb = -dT1_dVg * (dVth_dVb + Vgst / n * dn_dVb)
                  + T1 / n * dn_dVb; 
      double dT1_dVd = -dT1_dVg * (dVth_dVd + Vgst / n * dn_dVd)
                  + T1 / n * dn_dVd;

      double dT2_dVg = -m->cox / (t->vtm * t->cdep0)
                  * exp(ExpArg);
      double T2 = 1.0 - T10 * dT2_dVg;
      double dT2_dVd = -dT2_dVg * (dVth_dVd - 2.0 * t->vtm * ExpArg * dn_dVd)
                  + (T2 - 1.0) / n * dn_dVd;
      double dT2_dVb = -dT2_dVg * (dVth_dVb - 2.0 * t->vtm * ExpArg * dn_dVb)
                  + (T2 - 1.0) / n * dn_dVb;

      d->vgst = T1 / T2;
      double T3 = T2 * T2;
      dVgsteff_dVg = (T2 * dT1_dVg - T1 * dT2_dVg) / T3 * dVgs_eff_dVg;
      dVgsteff_dVd = (T2 * dT1_dVd - T1 * dT2_dVd) / T3;
      dVgsteff_dVb = (T2 * dT1_dVb - T1 * dT2_dVb) / T3;
      }}
      Vgst2Vtm = d->vgst + 2.0 * t->vtm;
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    /* Calculate Effective Channel Geometry */
    double Weff, dWeff_dVg, dWeff_dVb;
    {
      double T9 = sqrtPhis - t->sqrtPhi;
      Weff = s->weff - 2.0 * (s->dwg * d->vgst + s->dwb * T9); 
      dWeff_dVg = -2.0 * s->dwg;
      dWeff_dVb = -2.0 * s->dwb * dsqrtPhis_dVb;

      if (Weff < 2.0e-8) {
      /* to avoid the discontinuity problem due to Weff*/
      double T0 = 1.0 / (6.0e-8 - 2.0 * Weff);
      Weff = 2.0e-8 * (4.0e-8 - Weff) * T0;
      T0 *= T0 * 4.0e-16;
      dWeff_dVg *= T0;
      dWeff_dVb *= T0;
      }
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    double Rds, dRds_dVg, dRds_dVb;
    {
      double T9 = sqrtPhis - t->sqrtPhi;
      double T0 = s->prwg * d->vgst + s->prwb * T9;
      {if (T0 >= -0.9) {
      Rds = t->rds0 * (1.0 + T0);
      dRds_dVg = t->rds0 * s->prwg;
      dRds_dVb = t->rds0 * s->prwb * dsqrtPhis_dVb;
      }else{
      /* to avoid the discontinuity problem due to prwg and prwb*/
      double T1 = 1.0 / (17.0 + 20.0 * T0);
      Rds = t->rds0 * (0.8 + T0) * T1;
      T1 *= T1;
      dRds_dVg = t->rds0 * s->prwg * T1;
      dRds_dVb = t->rds0 * s->prwb * dsqrtPhis_dVb * T1;
      }}
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    /* Calculate Abulk */
    double Abulk0, dAbulk0_dVb, dAbulk_dVg, Abulk, dAbulk_dVb;
    {
      double T1 = 0.5 * t->k1 / sqrtPhis;
      double dT1_dVb = -T1 / sqrtPhis * dsqrtPhis_dVb;

      double T9 = sqrt(s->xj * Xdep);
      double tmp1 = s->leff + 2.0 * T9;
      double T5 = s->leff / tmp1; 
      double tmp2 = s->a0 * T5;
      double tmp3 = s->weff + s->b1; 
      double tmp4 = s->b0 / tmp3;
      double T2 = tmp2 + tmp4;
      double dT2_dVb = -T9 / tmp1 / Xdep * dXdep_dVb;
      double T6 = T5 * T5;
      double T7 = T5 * T6;

      Abulk0 = 1.0 + T1 * T2; 
      dAbulk0_dVb = T1 * tmp2 * dT2_dVb + T2 * dT1_dVb;

      double T8 = s->ags * s->a0 * T7;
      dAbulk_dVg = -T1 * T8;
      Abulk = Abulk0 + dAbulk_dVg * d->vgst; 
      dAbulk_dVb = dAbulk0_dVb - T8 * d->vgst * (dT1_dVb + 3.0 * T1 * dT2_dVb);

      if (Abulk0 < 0.1) {
      /* added to avoid the problems caused by Abulk0 */
      double T9 = 1.0 / (3.0 - 20.0 * Abulk0);
      Abulk0 = (0.2 - Abulk0) * T9;
      dAbulk0_dVb *= T9 * T9;
      }
      if (Abulk < 0.1) {
      /* added to avoid the problems caused by Abulk */
      double T9 = 1.0 / (3.0 - 20.0 * Abulk);
      Abulk = (0.2 - Abulk) * T9;
      dAbulk_dVb *= T9 * T9;
      }

      double T0, dT0_dVb;
      {
      double T2 = s->keta * Vbseff;
      {if (T2 >= -0.9) {
        T0 = 1.0 / (1.0 + T2);
        dT0_dVb = -s->keta * T0 * T0;
      }else{
        /* added to avoid the problems caused by Keta */
        double T1 = 1.0 / (0.8 + T2);
        T0 = (17.0 + 20.0 * T2) * T1;
        dT0_dVb = -s->keta * T1 * T1;
      }}
      }
      dAbulk_dVg *= T0;
      dAbulk_dVb = dAbulk_dVb * T0 + Abulk * dT0_dVb;
      dAbulk0_dVb = dAbulk0_dVb * T0 + Abulk0 * dT0_dVb;
      Abulk *= T0;
      Abulk0 *= T0;
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    /* Mobility calculation */
    double ueff, dueff_dVg, dueff_dVd, dueff_dVb;
    {
      double Denomi, dDenomi_dVg, dDenomi_dVd, dDenomi_dVb;
      {
      double T5;
      {if (m->mobMod == 1) {
        double T0 = d->vgst + d->von + d->von;
        double T2 = t->ua + t->uc * Vbseff;
        double T3 = T0 / m->tox;
        T5 = T3 * (T2 + t->ub * T3);
        dDenomi_dVg = (T2 + 2.0 * t->ub * T3) / m->tox;
        dDenomi_dVd = dDenomi_dVg * 2.0 * dVth_dVd;
        dDenomi_dVb = dDenomi_dVg * 2.0 * dVth_dVb + t->uc * T3;
      }else if (m->mobMod == 2) {
        T5 = d->vgst / m->tox 
          * (t->ua + t->uc * Vbseff + t->ub * d->vgst / m->tox);
        dDenomi_dVg = (t->ua + t->uc * Vbseff + 2.0 * t->ub * d->vgst / m->tox)
          / m->tox;
        dDenomi_dVd = 0.0;
        dDenomi_dVb = d->vgst * t->uc / m->tox; 
      }else{
        double T0 = d->vgst + d->von + d->von;
        double T2 = 1.0 + t->uc * Vbseff;
        double T3 = T0 / m->tox;
        double T4 = T3 * (t->ua + t->ub * T3);
        T5 = T4 * T2;
        dDenomi_dVg = (t->ua + 2.0 * t->ub * T3) * T2 / m->tox;
        dDenomi_dVd = dDenomi_dVg * 2.0 * dVth_dVd;
        dDenomi_dVb = dDenomi_dVg * 2.0 * dVth_dVb + t->uc * T4;
      }}
      {if (T5 >= -0.8) {
        Denomi = 1.0 + T5;
      }else{
        /* Added to avoid the discontinuity problem caused by ua and ub*/ 
        double T9 = 1.0 / (7.0 + 10.0 * T5);
        Denomi = (0.6 + T5) * T9;
        T9 *= T9;
        dDenomi_dVg *= T9;
        dDenomi_dVd *= T9;
        dDenomi_dVb *= T9;
      }}
      }
      ueff = t->u0temp / Denomi;
      double T9 = -ueff / Denomi;
      dueff_dVg = T9 * dDenomi_dVg;
      dueff_dVd = T9 * dDenomi_dVd;
      dueff_dVb = T9 * dDenomi_dVb;
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    double Esat, EsatL, dEsatL_dVg, dEsatL_dVd, dEsatL_dVb;
    {
      Esat = 2.0 * t->vsattemp / ueff;
      EsatL = Esat * s->leff;
      double T0 = -EsatL /ueff;
      dEsatL_dVg = T0 * dueff_dVg;
      dEsatL_dVd = T0 * dueff_dVd;
      dEsatL_dVb = T0 * dueff_dVb;
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    double Vdsat, dVdsat_dVg, dVdsat_dVd, dVdsat_dVb; // d->vdsat
    double Vasat, dVasat_dVg, dVasat_dVb, dVasat_dVd;
    {
      double WVCoxRds;
      {
      double WVCox = Weff * t->vsattemp * m->cox;
      WVCoxRds = WVCox * Rds;
      }
      double Lambda, dLambda_dVg;
      {
      {if (s->a1 == 0.0) {
        Lambda = s->a2;
        dLambda_dVg = 0.0;
      }else if (s->a1 > 0.0) {
        /* avoid discontinuity problem caused by s->a1 and s->a2 (Lambda) */
        double T0 = 1.0 - s->a2;
        double T1 = T0 - s->a1 * d->vgst - 0.0001;
        double T2 = sqrt(T1 * T1 + 0.0004 * T0);
        Lambda = s->a2 + T0 - 0.5 * (T1 + T2);
        dLambda_dVg = 0.5 * s->a1 * (1.0 + T1 / T2);
      }else{
        double T1 = s->a2 + s->a1 * d->vgst - 0.0001;
        double T2 = sqrt(T1 * T1 + 0.0004 * s->a2);
        Lambda = 0.5 * (T1 + T2);
        dLambda_dVg = 0.5 * s->a1 * (1.0 + T1 / T2);
      }}
      }
      double tmp2, tmp3;
      {if (Rds > 0) {
      tmp2 = dRds_dVg / Rds + dWeff_dVg / Weff;
      tmp3 = dRds_dVb / Rds + dWeff_dVb / Weff;
      }else{
      tmp2 = dWeff_dVg / Weff;
      tmp3 = dWeff_dVb / Weff;
      }}
      double tmp1;
      {
      {if ((Rds == 0.0) && (Lambda == 1.0)) {
        double T0 = 1.0 / (Abulk * EsatL + Vgst2Vtm);
        tmp1 = 0.0;
        double T1 = T0 * T0;
        double T2 = Vgst2Vtm * T0;
        double T3 = EsatL * Vgst2Vtm;
        Vdsat = T3 * T0;
        double dT0_dVg = -(Abulk * dEsatL_dVg + EsatL * dAbulk_dVg + 1.0)*T1;
        double dT0_dVd = -(Abulk * dEsatL_dVd) * T1; 
        double dT0_dVb = -(Abulk * dEsatL_dVb + dAbulk_dVb * EsatL) * T1;   
        dVdsat_dVg = T3 * dT0_dVg + T2 * dEsatL_dVg + EsatL * T0;
        dVdsat_dVd = T3 * dT0_dVd + T2 * dEsatL_dVd;
        dVdsat_dVb = T3 * dT0_dVb + T2 * dEsatL_dVb;   
      }else{
        tmp1 = dLambda_dVg / (Lambda * Lambda);
        double T9 = Abulk * WVCoxRds;
        double T8 = Abulk * T9;
        double T7 = Vgst2Vtm * T9;
        double T6 = Vgst2Vtm * WVCoxRds;
        double T0 = 2.0 * Abulk * (T9 - 1.0 + 1.0 / Lambda); 
        double dT0_dVg = 2.0 * (T8 * tmp2 - Abulk * tmp1
                  + (2.0 * T9 + 1.0 / Lambda - 1.0) * dAbulk_dVg);
        double dT0_dVb = 2.0 * (T8 * (2.0 / Abulk * dAbulk_dVb + tmp3)
                          + (1.0 / Lambda - 1.0) * dAbulk_dVb);
        //double dT0_dVd = 0.0;
        
        double T1 = Vgst2Vtm * (2.0 / Lambda - 1.0) + Abulk * EsatL + 3.0*T7;
        double dT1_dVg = (2.0 / Lambda - 1.0) - 2.0 * Vgst2Vtm * tmp1
          + Abulk * dEsatL_dVg + EsatL * dAbulk_dVg 
          + 3.0 * (T9 + T7 * tmp2 + T6 * dAbulk_dVg);
        double dT1_dVb = Abulk * dEsatL_dVb + EsatL * dAbulk_dVb
          + 3.0 * (T6 * dAbulk_dVb + T7 * tmp3);
        double dT1_dVd = Abulk * dEsatL_dVd;
        
        double T2 = Vgst2Vtm * (EsatL + 2.0 * T6);
        double dT2_dVg = EsatL + Vgst2Vtm * dEsatL_dVg
          + T6 * (4.0 + 2.0 * Vgst2Vtm * tmp2);
        double dT2_dVb = Vgst2Vtm * (dEsatL_dVb + 2.0 * T6 * tmp3);
        double dT2_dVd = Vgst2Vtm * dEsatL_dVd;
        
        double T3 = sqrt(T1 * T1 - 2.0 * T0 * T2);
        Vdsat = (T1 - T3) / T0;
        dVdsat_dVg = (dT1_dVg - (T1 * dT1_dVg - dT0_dVg * T2
                        - T0 * dT2_dVg) / T3 - Vdsat * dT0_dVg) / T0;
        dVdsat_dVb = (dT1_dVb - (T1 * dT1_dVb - dT0_dVb * T2
                        - T0 * dT2_dVb) / T3 - Vdsat * dT0_dVb) / T0;
        dVdsat_dVd = (dT1_dVd - (T1 * dT1_dVd - T0 * dT2_dVd) / T3) / T0;
      }}
      d->vdsat = Vdsat;
      }
      /* Calculate VAsat */
      {
      double tmp4 = 1.0 - 0.5 * Abulk * Vdsat / Vgst2Vtm;
      double T9 = WVCoxRds * d->vgst;
      double T8 = T9 / Vgst2Vtm;
      double T0 = EsatL + Vdsat + 2.0 * T9 * tmp4;
      double T7 = 2.0 * WVCoxRds * tmp4;
      double dT0_dVg = dEsatL_dVg + dVdsat_dVg + T7 * (1.0 + tmp2 * d->vgst)
        - T8 * (Abulk * dVdsat_dVg - Abulk * Vdsat / Vgst2Vtm
              + Vdsat * dAbulk_dVg);   
      double dT0_dVb = dEsatL_dVb + dVdsat_dVb + T7 * tmp3 * d->vgst
        - T8 * (dAbulk_dVb * Vdsat + Abulk * dVdsat_dVb);
      double dT0_dVd = dEsatL_dVd + dVdsat_dVd - T8 * Abulk * dVdsat_dVd;
      T9 = WVCoxRds * Abulk; 
      double T1 = 2.0 / Lambda - 1.0 + T9; 
      double dT1_dVg = -2.0 * tmp1 +  WVCoxRds * (Abulk * tmp2 + dAbulk_dVg);
      double dT1_dVb = dAbulk_dVb * WVCoxRds + T9 * tmp3;
      Vasat = T0 / T1;
      dVasat_dVg = (dT0_dVg - Vasat * dT1_dVg) / T1;
      dVasat_dVb = (dT0_dVb - Vasat * dT1_dVb) / T1;
      dVasat_dVd = dT0_dVd / T1;
      }
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    /* Effective Vds (Vdseff) Calculation */
    double Vdseff, diffVds, dVdseff_dVg, dVdseff_dVd, dVdseff_dVb;
    {
      double T1 = Vdsat - d->vds - s->delta;
      double dT1_dVg = dVdsat_dVg;
      double dT1_dVd = dVdsat_dVd - 1.0;
      double dT1_dVb = dVdsat_dVb;
      trace4("", T1, dT1_dVg, dT1_dVd, dT1_dVb);
      
      double T2 = sqrt(T1 * T1 + 4.0 * s->delta * Vdsat);
      double T0 = T1 / T2;
      double T3 = 2.0 * s->delta / T2;
      trace3("", T2, T0, T3);
      double dT2_dVg = T0 * dT1_dVg + T3 * dVdsat_dVg;
      double dT2_dVd = T0 * dT1_dVd + T3 * dVdsat_dVd;
      double dT2_dVb = T0 * dT1_dVb + T3 * dVdsat_dVb;
      trace3("", dT2_dVg, dT2_dVd, dT2_dVb);
      
      Vdseff      = Vdsat - 0.5 * (T1 + T2);
      dVdseff_dVg = dVdsat_dVg - 0.5 * (dT1_dVg + dT2_dVg); 
      dVdseff_dVd = dVdsat_dVd - 0.5 * (dT1_dVd + dT2_dVd); 
      dVdseff_dVb = dVdsat_dVb - 0.5 * (dT1_dVb + dT2_dVb); 
      trace4("raw", Vdseff, dVdseff_dVg, dVdseff_dVd, dVdseff_dVb);

      fixzero(&Vdseff,      Vdsat);
      fixzero(&dVdseff_dVg, dVdsat_dVg);
      fixzero(&dVdseff_dVd, dVdsat_dVd);
      fixzero(&dVdseff_dVb, dVdsat_dVb);
      if (Vdseff > d->vds) {  // From Spice, to fix numeric problems.
      Vdseff = d->vds;
      }
      trace4("fixed", Vdseff, dVdseff_dVg, dVdseff_dVd, dVdseff_dVb);

      diffVds = d->vds - Vdseff;
      trace2("", Vdseff, diffVds);
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    /* Calculate Ids */
    double Idsa, dIdsa_dVg, dIdsa_dVd, dIdsa_dVb;
    {
      double Va, dVa_dVg, dVa_dVd, dVa_dVb;
      {
      double VACLM, dVACLM_dVg, dVACLM_dVb, dVACLM_dVd;
      {if ((s->pclm > 0.0) && (diffVds > 1.0e-10)) {
        double T0 = 1.0 / (s->pclm * Abulk * s->litl);
        double dT0_dVb = -T0 / Abulk * dAbulk_dVb;
        double dT0_dVg = -T0 / Abulk * dAbulk_dVg; 
        double T2 = d->vgst / EsatL;
        double T1 = s->leff * (Abulk + T2); 
        double dT1_dVg = s->leff * ((1.0-T2*dEsatL_dVg)/EsatL + dAbulk_dVg);
        double dT1_dVb = s->leff * (dAbulk_dVb - T2 * dEsatL_dVb / EsatL);
        double dT1_dVd = -T2 * dEsatL_dVd / Esat;
        double T9 = T0 * T1;
        VACLM = T9 * diffVds;
        dVACLM_dVg = T0 * dT1_dVg * diffVds - T9 * dVdseff_dVg
          + T1 * diffVds * dT0_dVg;
        dVACLM_dVb = (dT0_dVb*T1 + T0*dT1_dVb) * diffVds - T9 * dVdseff_dVb;
        dVACLM_dVd = T0 * dT1_dVd * diffVds + T9 * (1.0 - dVdseff_dVd);
      }else{
        VACLM = MAX_EXP;
        dVACLM_dVd = dVACLM_dVg = dVACLM_dVb = 0.0;
      }}
      double VADIBL, dVADIBL_dVg, dVADIBL_dVb, dVADIBL_dVd;
      {if (t->thetaRout > 0.0) {
        double T8 = Abulk * Vdsat;
        double T0 = Vgst2Vtm * T8;
        double dT0_dVg = Vgst2Vtm * Abulk * dVdsat_dVg + T8
          + Vgst2Vtm * Vdsat * dAbulk_dVg;
        double dT0_dVb = Vgst2Vtm * (dAbulk_dVb*Vdsat + Abulk*dVdsat_dVb);
        double dT0_dVd = Vgst2Vtm * Abulk * dVdsat_dVd;
        double T1 = Vgst2Vtm + T8;
        double dT1_dVg = 1.0 + Abulk * dVdsat_dVg + Vdsat * dAbulk_dVg;
        double dT1_dVb = Abulk * dVdsat_dVb + dAbulk_dVb * Vdsat;
        double dT1_dVd = Abulk * dVdsat_dVd;
        double T9 = T1 * T1;
        double T2 = t->thetaRout;
        VADIBL = (Vgst2Vtm - T0 / T1) / T2;
        dVADIBL_dVg = (1.0 - dT0_dVg / T1 + T0 * dT1_dVg / T9) / T2;
        dVADIBL_dVb = (-dT0_dVb / T1 + T0 * dT1_dVb / T9) / T2;
        dVADIBL_dVd = (-dT0_dVd / T1 + T0 * dT1_dVd / T9) / T2;
        
        double T7 = s->pdiblb * Vbseff;
        {if (T7 >= -0.9) {
          double T3 = 1.0 / (1.0 + T7);
          VADIBL *= T3;
          dVADIBL_dVg *= T3;
          dVADIBL_dVb = (dVADIBL_dVb - VADIBL * s->pdiblb) * T3;
          dVADIBL_dVd *= T3;
        }else{
          /* Added to avoid the discontinuity problem caused by pdiblcb */
          double T4 = 1.0 / (0.8 + T7);
          double T3 = (17.0 + 20.0 * T7) * T4;
          dVADIBL_dVg *= T3;
          dVADIBL_dVb = dVADIBL_dVb * T3 - VADIBL * s->pdiblb * T4 * T4;
          dVADIBL_dVd *= T3;
          VADIBL *= T3;
        }}
      }else{
        VADIBL = MAX_EXP;
        dVADIBL_dVd = dVADIBL_dVg = dVADIBL_dVb = 0.0;
      }}
      double T8 = s->pvag / EsatL;
      double T9 = T8 * d->vgst;
      double T0, dT0_dVg, dT0_dVb, dT0_dVd;
      {if (T9 > -0.9) {
        T0 = 1.0 + T9;
        dT0_dVg = T8 * (1.0 - d->vgst * dEsatL_dVg / EsatL);
        dT0_dVb = -T9 * dEsatL_dVb / EsatL;
        dT0_dVd = -T9 * dEsatL_dVd / EsatL;
      }else{
        /* Added to avoid the discontinuity problems caused by pvag */
        double T1 = 1.0 / (17.0 + 20.0 * T9);
        T0 = (0.8 + T9) * T1;
        T1 *= T1;
        dT0_dVg = T8 * (1.0 - d->vgst * dEsatL_dVg / EsatL) * T1;
        T9 *= T1 / EsatL;
        dT0_dVb = -T9 * dEsatL_dVb;
        dT0_dVd = -T9 * dEsatL_dVd;
      }}
      double tmp1 = VACLM * VACLM;
      double tmp2 = VADIBL * VADIBL;
      double tmp3 = VACLM + VADIBL;
      
      double T1 = VACLM * VADIBL / tmp3;
      tmp3 *= tmp3;
      double dT1_dVg = (tmp1 * dVADIBL_dVg + tmp2 * dVACLM_dVg) / tmp3;
      double dT1_dVd = (tmp1 * dVADIBL_dVd + tmp2 * dVACLM_dVd) / tmp3;
      double dT1_dVb = (tmp1 * dVADIBL_dVb + tmp2 * dVACLM_dVb) / tmp3;
      
      Va = Vasat + T0 * T1;
      dVa_dVg = dVasat_dVg + T1 * dT0_dVg + T0 * dT1_dVg;
      dVa_dVd = dVasat_dVd + T1 * dT0_dVd + T0 * dT1_dVd;
      dVa_dVb = dVasat_dVb + T1 * dT0_dVb + T0 * dT1_dVb;
      }
      double Idl, dIdl_dVg, dIdl_dVd, dIdl_dVb;
      {
      double gche, dgche_dVg, dgche_dVd, dgche_dVb;
      {
        double beta, dbeta_dVg, dbeta_dVd, dbeta_dVb;
        {
          double CoxWovL = m->cox * Weff / s->leff;
          beta = ueff * CoxWovL;
          dbeta_dVg = CoxWovL * dueff_dVg + beta * dWeff_dVg / Weff;
          dbeta_dVd = CoxWovL * dueff_dVd;
          dbeta_dVb = CoxWovL * dueff_dVb + beta * dWeff_dVb / Weff;
        }
        double fgche1, dfgche1_dVg, dfgche1_dVd, dfgche1_dVb;
        {
          double T0 = 1.0 - 0.5 * Abulk * Vdseff / Vgst2Vtm;
          double dT0_dVg = -0.5 * (Abulk * dVdseff_dVg 
            - Abulk * Vdseff / Vgst2Vtm + Vdseff * dAbulk_dVg) / Vgst2Vtm;
          double dT0_dVd = -0.5 * Abulk * dVdseff_dVd / Vgst2Vtm;
          double dT0_dVb = -0.5 * (Abulk*dVdseff_dVb + dAbulk_dVb*Vdseff) 
            / Vgst2Vtm;
          fgche1 = d->vgst * T0;
          dfgche1_dVg = d->vgst * dT0_dVg + T0; 
          dfgche1_dVd = d->vgst * dT0_dVd; 
          dfgche1_dVb = d->vgst * dT0_dVb; 
        }
        double fgche2, dfgche2_dVg, dfgche2_dVd, dfgche2_dVb;
        {
          double T9 = Vdseff / EsatL;
          fgche2 = 1.0 + T9;
          dfgche2_dVg = (dVdseff_dVg - T9 * dEsatL_dVg) / EsatL;
          dfgche2_dVd = (dVdseff_dVd - T9 * dEsatL_dVd) / EsatL;
          dfgche2_dVb = (dVdseff_dVb - T9 * dEsatL_dVb) / EsatL;
        }
        gche = beta * fgche1 / fgche2;
        dgche_dVg = (beta * dfgche1_dVg + fgche1 * dbeta_dVg
                   - gche * dfgche2_dVg) / fgche2;
        dgche_dVd = (beta * dfgche1_dVd + fgche1 * dbeta_dVd
                   - gche * dfgche2_dVd) / fgche2;
        dgche_dVb = (beta * dfgche1_dVb + fgche1 * dbeta_dVb
                   - gche * dfgche2_dVb) / fgche2;
      }
      double T0 = 1.0 + gche * Rds;
      double T9 = Vdseff / T0;
      Idl = gche * T9;
      dIdl_dVg = (gche * dVdseff_dVg + T9 * dgche_dVg) / T0
        - Idl * gche / T0 * dRds_dVg; 
      dIdl_dVd = (gche * dVdseff_dVd + T9 * dgche_dVd) / T0; 
      dIdl_dVb = (gche*dVdseff_dVb + T9*dgche_dVb - Idl*dRds_dVb*gche) / T0; 
      }
      double T9 =  diffVds / Va;
      double T0 =  1.0 + T9;
      Idsa = Idl * T0;
      dIdsa_dVg = T0 * dIdl_dVg - Idl * (dVdseff_dVg + T9 * dVa_dVg) / Va;
      dIdsa_dVd = T0 * dIdl_dVd + Idl * (1.0 - dVdseff_dVd - T9*dVa_dVd) / Va;
      dIdsa_dVb = T0 * dIdl_dVb - Idl * (dVdseff_dVb + T9 * dVa_dVb) / Va;
      trace4("", Idsa, dIdsa_dVg, dIdsa_dVd, dIdsa_dVb);
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    // d->ids, d->gds, d->gmf, d->gmbf
    {
      double VASCBE, dVASCBE_dVg, dVASCBE_dVd, dVASCBE_dVb;
      {if (s->pscbe2 > 0.0) {
      {if (diffVds > s->pscbe1 * s->litl / EXP_THRESHOLD) {
        double T0 =  s->pscbe1 * s->litl / diffVds;
        VASCBE = s->leff * exp(T0) / s->pscbe2;
        double T1 = T0 * VASCBE / diffVds;
        dVASCBE_dVg = T1 * dVdseff_dVg;
        dVASCBE_dVd = -T1 * (1.0 - dVdseff_dVd);
        dVASCBE_dVb = T1 * dVdseff_dVb;
      }else{
        VASCBE = MAX_EXP * s->leff/s->pscbe2;
        dVASCBE_dVg = dVASCBE_dVd = dVASCBE_dVb = 0.0;
      }}
      }else{
      VASCBE = MAX_EXP;
      dVASCBE_dVg = dVASCBE_dVd = dVASCBE_dVb = 0.0;
      }}
      double T9 = diffVds / VASCBE;
      double T0 = 1.0 + T9;
      double Ids = Idsa * T0;
      double Gm = T0*dIdsa_dVg - Idsa*(dVdseff_dVg + T9*dVASCBE_dVg) / VASCBE;
      double Gds = T0 * dIdsa_dVd 
      + Idsa * (1.0 - dVdseff_dVd - T9 * dVASCBE_dVd) / VASCBE;
      double Gmb = T0 * dIdsa_dVb
      - Idsa * (dVdseff_dVb + T9 * dVASCBE_dVb) / VASCBE;
      trace3("", T0, dIdsa_dVb, (T0 * dIdsa_dVb));
      trace4("", dVdseff_dVb, T9, dVASCBE_dVb, (dVdseff_dVb + T9*dVASCBE_dVb));
      trace3("", Idsa, VASCBE, (Idsa*(dVdseff_dVb+T9*dVASCBE_dVb)/VASCBE));

      Gds += Gm * dVgsteff_dVd;
      Gmb += Gm * dVgsteff_dVb;
      Gm *= dVgsteff_dVg;
      Gmb *= dVbseff_dVb;
      trace4("", Ids, Gm, Gds, Gmb);
      trace0("=========================");

      d->gds = Gds;
      {if (d->reversed){
      d->ids  = -Ids;
      d->gmr  = Gm;
      d->gmbr = Gmb;
      d->gmf = d->gmbf = 0;
      }else{
      d->ids  = Ids;
      d->gmf  = Gm;
      d->gmbf = Gmb;
      d->gmr = d->gmbr = 0.;
      }}
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    // d->isub, d->gbbs, d->gbgs, d->gbds
    {
      /* calculate substrate current Isub */
      double Isub, Gbd, Gbb, Gbg;
      {if ((s->alpha0 <= 0.0) || (s->beta0 <= 0.0)) {
      Isub = Gbd = Gbb = Gbg = 0.0;
      trace4("no-isub", Isub, Gbd, Gbb, Gbg);
      }else{
      double T2 = s->alpha0 / s->leff;
      double T1, dT1_dVg, dT1_dVd, dT1_dVb;
      {if (diffVds > s->beta0 / EXP_THRESHOLD) {
        double T0 = -s->beta0 / diffVds;
        T1 = T2 * diffVds * exp(T0);
        double T3 = T1 / diffVds * (T0 - 1.0);
        trace3("", T0, T2, T3);
        dT1_dVg = T3 * dVdseff_dVg;
        dT1_dVd = -T3 * (1.0 - dVdseff_dVd);
        dT1_dVb = T3 * dVdseff_dVb;
        trace4("vds > ?", T1, dT1_dVg, dT1_dVd, dT1_dVb);
      }else{
        double T3 = T2 * MIN_EXP;
        trace2("", T2, T3);
        T1 = T3 * diffVds;
        dT1_dVg = -T3 * dVdseff_dVg;
        dT1_dVd = T3 * (1.0 - dVdseff_dVd);
        dT1_dVb = -T3 * dVdseff_dVb;
        trace4("vds < ?", T1, dT1_dVg, dT1_dVd, dT1_dVb);
      }}
      Isub = T1 * Idsa;
      Gbg = T1 * dIdsa_dVg + Idsa * dT1_dVg;
      Gbd = T1 * dIdsa_dVd + Idsa * dT1_dVd;
      Gbb = T1 * dIdsa_dVb + Idsa * dT1_dVb;
      trace4("raw", Isub, Gbd, Gbb, Gbg);
      
      Gbd += Gbg * dVgsteff_dVd;
      Gbb += Gbg * dVgsteff_dVb;
      Gbg *= dVgsteff_dVg;
      Gbb *= dVbseff_dVb; /* bug fixing */
      }}
      trace4("", Isub, Gbd, Gbb, Gbg);
      d->isub = Isub;
      d->gbbs = Gbb;
      d->gbgs = Gbg;
      d->gbds = Gbd;
      //double d__csub = Isub - (Gbb * Vbseff + Gbd * d->vds + Gbg * d->vgs);
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    /* Calculate Qinv for Noise analysis */
    {
      //double T1 = d->vgst * (1.0 - 0.5 * Abulk * Vdseff / Vgst2Vtm);
      //double d__qinv = -m->cox * Weff * s->leff * T1;
    }
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    // ends line 2020 (finished)
    {
      const bool ChargeComputationNeeded = true;
      {if ((m->xpart < 0) || (!ChargeComputationNeeded)) {
      d->qgate = d->qdrn = d->qbulk = 0.0;
      d->cggb = d->cgsb = d->cgdb = 0.0;
      d->cdgb = d->cdsb = d->cddb = 0.0;
      d->cbgb = d->cbsb = d->cbdb = 0.0;
      //goto finished;
      }else if (m->capMod == 0) {
      // block ends 1621 this 1309
      {if (Vbseff < 0.0) {  // redefinition
        Vbseff = d->vbs;
        dVbseff_dVb = 1.0;
      }else{
        Vbseff = t->phi - Phis;
        dVbseff_dVb = -dPhis_dVb;
      }}
      double Vfb = s->vfbcv; // possible improper redefinition later
      double Vth = Vfb + t->phi + t->k1 * sqrtPhis; 
      dVth_dVb = t->k1 * dsqrtPhis_dVb; // redefinition
      double Vgst = Vgs_eff - Vth;
      //double dVgst_dVb = -dVth_dVb;
      //double dVgst_dVg = dVgs_eff_dVg; 
      double CoxWL = m->cox * s->weffCV * s->leffCV;
      double Arg1 = Vgs_eff - Vbseff - Vfb;
      // ends 1618 this 1328
      {if (Arg1 <= 0.0) {
        d->qgate = CoxWL * Arg1;
        d->cggb = CoxWL * dVgs_eff_dVg;
        d->cgdb = 0.0;
        d->cgsb = CoxWL * (dVbseff_dVb - dVgs_eff_dVg);
        
        d->qbulk = -d->qgate;
        d->cdgb = 0.0;
        d->cddb = 0.0;
        d->cdsb = 0.0;
        
        d->qdrn = 0.0;
        d->cbgb = -CoxWL * dVgs_eff_dVg;
        d->cbdb = 0.0;
        d->cbsb = -d->cgsb;
      }else if (Vgst <= 0.0) {
        double T1 = 0.5 * t->k1;
        double T2 = sqrt(T1 * T1 + Arg1);
        double T0 = CoxWL * T1 / T2;
        d->qgate = CoxWL * t->k1 * (T2 - T1);
        d->cggb = T0 * dVgs_eff_dVg;
        d->cgdb = 0.0;
        d->cgsb = T0 * (dVbseff_dVb - dVgs_eff_dVg);
        
        d->qbulk = -d->qgate;
        d->cbgb = -d->cggb;
        d->cbdb = 0.0;
        d->cbsb = -d->cgsb;
   
        d->qdrn = 0.0;
        d->cdgb = 0.0;
        d->cddb = 0.0;
        d->cdsb = 0.0;
      }else{
        double One_Third_CoxWL = CoxWL / 3.0;
        double Two_Third_CoxWL = 2.0 * One_Third_CoxWL;
        // redefine Vdsat, dVdsat_dVg, dVdsat_dVb
        {
          double AbulkCV = Abulk0 * s->abulkCVfactor;
          double dAbulkCV_dVb = s->abulkCVfactor * dAbulk0_dVb;
          Vdsat = Vgst / AbulkCV;
          dVdsat_dVg = dVgs_eff_dVg / AbulkCV;
          dVdsat_dVb = - (Vdsat * dAbulkCV_dVb + dVth_dVb)/ AbulkCV; 
        }
        {if (m->xpart > 0.5) {
          /* 0/100 Charge petition model */
          {if (d->vds >= Vdsat) {
            /* saturation region */
            double T1 = Vdsat / 3.0;
            double T2 = -One_Third_CoxWL * dVdsat_dVb;
            d->qgate = CoxWL * (Vgs_eff - Vfb - t->phi - T1);
            d->cggb = One_Third_CoxWL * (3.0 - dVdsat_dVg) * dVgs_eff_dVg;
            d->cgsb = -(d->cggb + T2);
            d->cgdb = 0.0;
            
            double T2a = -Two_Third_CoxWL * Vgst;
            double T3 = -(T2 + Two_Third_CoxWL * dVth_dVb);
            d->qbulk = -(d->qgate + T2a);
            d->cbgb = -(d->cggb - Two_Third_CoxWL * dVgs_eff_dVg);
            d->cbsb = -(d->cbgb + T3);
            d->cbdb = 0.0;

            d->qdrn = 0.0;
            d->cdgb = 0.0;
            d->cddb = 0.0;
            d->cdsb = 0.0;
          }else{
            /* linear region */
            double Alphaz = Vgst / Vdsat;
            double T1 = 2.0 * Vdsat - d->vds;
            double T2 = d->vds / (3.0 * T1);
            double T3 = T2 * d->vds;
            double T9 = 0.25 * CoxWL;
            double T4 = T9 * Alphaz;
            double T7 = 2.0 * d->vds - T1 - 3.0 * T3;
            double T8 = T3 - T1 - 2.0 * d->vds;
            d->qgate = CoxWL * (Vgs_eff - Vfb - t->phi - 0.5 * (d->vds - T3));
            double T10 = T4 * T8;
            d->qdrn = T4 * T7;
            d->qbulk = -(d->qgate + d->qdrn + T10);
            
            double T5 = T3 / T1;
            d->cggb = CoxWL * (1.0 - T5 * dVdsat_dVg) * dVgs_eff_dVg;
            double T11 = -CoxWL * T5 * dVdsat_dVb;
            d->cgdb = CoxWL * (T2 - 0.5 + 0.5 * T5);
            d->cgsb = -(d->cggb + T11 + d->cgdb);

            double T6 = 1.0 / Vdsat;
            double dAlphaz_dVg = T6 * (1.0 - Alphaz * dVdsat_dVg);
            double dAlphaz_dVb = -T6 * (dVth_dVb + Alphaz * dVdsat_dVb);
            T7 = T9 * T7;
            T8 = T9 * T8;
            T9 = 2.0 * T4 * (1.0 - 3.0 * T5);
            d->cdgb = (T7 * dAlphaz_dVg - T9 * dVdsat_dVg) * dVgs_eff_dVg;
            double T12 = T7 * dAlphaz_dVb - T9 * dVdsat_dVb;
            d->cddb = T4 * (3.0 - 6.0 * T2 - 3.0 * T5);
            d->cdsb = -(d->cdgb + T12 + d->cddb);

            T9 = 2.0 * T4 * (1.0 + T5);
            T10 = (T8 * dAlphaz_dVg - T9 * dVdsat_dVg) * dVgs_eff_dVg;
            T11 = T8 * dAlphaz_dVb - T9 * dVdsat_dVb;
            T12 = T4 * (2.0 * T2 + T5 - 1.0); 
            double T0 = -(T10 + T11 + T12);
            d->cbgb = -(d->cggb + d->cdgb + T10);
            d->cbdb = -(d->cgdb + d->cddb + T12);
            d->cbsb = -(d->cgsb + d->cdsb + T0);
          }}
        }else if (m->xpart < 0.5) {
          /* 40/60 Charge petition model */
          {if (d->vds >= Vdsat) {
            /* saturation region */
            double T1 = Vdsat / 3.0;
            d->qgate = CoxWL * (Vgs_eff - Vfb - t->phi - T1);
            double T2 = -Two_Third_CoxWL * Vgst;
            d->qbulk = -(d->qgate + T2);
            d->qdrn = 0.4 * T2;

            d->cggb = One_Third_CoxWL * (3.0 - dVdsat_dVg) * dVgs_eff_dVg;
            T2 = -One_Third_CoxWL * dVdsat_dVb;
            d->cgsb = -(d->cggb + T2);
            d->cgdb = 0.0;
       
            double T3 = 0.4 * Two_Third_CoxWL;
            d->cdgb = -T3 * dVgs_eff_dVg;
            d->cddb = 0.0;
            double T4 = T3 * dVth_dVb;
            d->cdsb = -(T4 + d->cdgb);
            
            d->cbgb = -(d->cggb - Two_Third_CoxWL * dVgs_eff_dVg);
            T3 = -(T2 + Two_Third_CoxWL * dVth_dVb);
            d->cbsb = -(d->cbgb + T3);
            d->cbdb = 0.0;
          }else{
            /* linear region  */
            double T1 = 2.0 * Vdsat - d->vds;
            double T2 = d->vds / (3.0 * T1);
            double T3 = T2 * d->vds;
            d->qgate = CoxWL * (Vgs_eff - Vfb - t->phi - 0.5 * (d->vds - T3));
            double T5 = T3 / T1;
            d->cggb = CoxWL * (1.0 - T5 * dVdsat_dVg) * dVgs_eff_dVg;
            double tmp = -CoxWL * T5 * dVdsat_dVb;
            d->cgdb = CoxWL * (T2 - 0.5 + 0.5 * T5);
            d->cgsb = -(d->cggb + d->cgdb + tmp);

            double T6 = 1.0 / Vdsat;
            double Alphaz      =  T6 * Vgst;
            double dAlphaz_dVg =  T6 * (1.0 - Alphaz * dVdsat_dVg);
            double dAlphaz_dVb = -T6 * (dVth_dVb + Alphaz * dVdsat_dVb);
            T6 = 8.0 * Vdsat * Vdsat - 6.0 * Vdsat * d->vds 
            + 1.2 * d->vds * d->vds;
            double T8 = T2 / T1;
            double T7 = d->vds - T1 - T8 * T6;
            double T9 = 0.25 * CoxWL;
            double T4 = T9 * Alphaz;
            d->qdrn = T4 * T7;
            T7 *= T9;
            tmp = T8 / T1;
            double tmp1 = T4 * (2.0 - 4.0 * tmp * T6
                     + T8 * (16.0 * Vdsat - 6.0 * d->vds));
            d->cdgb = (T7 * dAlphaz_dVg - tmp1 * dVdsat_dVg) * dVgs_eff_dVg;
            double T10 = T7 * dAlphaz_dVb - tmp1 * dVdsat_dVb;
            d->cddb = T4 * (2.0 - (1.0 / (3.0 * T1 * T1) + 2.0 * tmp) * T6 
                        + T8 * (6.0 * Vdsat - 2.4 * d->vds));
            d->cdsb = -(d->cdgb + T10 + d->cddb);

            T7 = 2.0 * (T1 + T3);
            d->qbulk = -(d->qgate - T4 * T7);
            T7 *= T9;
            double T0 = 4.0 * T4 * (1.0 - T5);
            double T12 = (-T7 * dAlphaz_dVg - d->cdgb - T0 * dVdsat_dVg)
            * dVgs_eff_dVg;
            double T11 = -T7 * dAlphaz_dVb - T10 - T0 * dVdsat_dVb;
            T10 = -4.0 * T4 * (T2 - 0.5 + 0.5 * T5) - d->cddb;
            tmp = -(T10 + T11 + T12);
            d->cbgb = -(d->cggb + d->cdgb + T12);
            d->cbdb = -(d->cgdb + d->cddb + T11);
            d->cbsb = -(d->cgsb + d->cdsb + tmp);
          }}
        }else{
          /* 50/50 partitioning */
          {if (d->vds >= Vdsat) {
            /* saturation region */
            double T1 = Vdsat / 3.0;
            d->qgate = CoxWL * (Vgs_eff - Vfb - t->phi - T1);
            double T2 = -Two_Third_CoxWL * Vgst;
            d->qbulk = -(d->qgate + T2);
            d->qdrn = 0.5 * T2;
            
            T2 = -One_Third_CoxWL * dVdsat_dVb;
            d->cggb = One_Third_CoxWL * (3.0 - dVdsat_dVg) * dVgs_eff_dVg;
            d->cgsb = -(d->cggb + T2);
            d->cgdb = 0.0;
            
            double T4 = One_Third_CoxWL * dVth_dVb;
            d->cdgb = -One_Third_CoxWL * dVgs_eff_dVg;
            d->cddb = 0.0;
            d->cdsb = -(T4 + d->cdgb);
            
            double T3 = -(T2 + Two_Third_CoxWL * dVth_dVb);
            d->cbgb = -(d->cggb - Two_Third_CoxWL * dVgs_eff_dVg);
            d->cbsb = -(d->cbgb + T3);
            d->cbdb = 0.0;
          }else{
            /* linear region */
            double T1 = 2.0 * Vdsat - d->vds;
            double T2 = d->vds / (3.0 * T1);
            double T3 = T2 * d->vds;
            double T5 = T3 / T1;
            double tmp = -CoxWL * T5 * dVdsat_dVb;
            d->qgate = CoxWL * (Vgs_eff - Vfb - t->phi - 0.5 * (d->vds-T3));
            d->cggb = CoxWL * (1.0 - T5 * dVdsat_dVg) * dVgs_eff_dVg;
            d->cgdb = CoxWL * (T2 - 0.5 + 0.5 * T5);
            d->cgsb = -(d->cggb + d->cgdb + tmp);
            
            double T6 = 1.0 / Vdsat;
            double Alphaz =       T6 * Vgst;
            double dAlphaz_dVg =  T6 * (1.0 - Alphaz * dVdsat_dVg);
            double dAlphaz_dVb = -T6 * (dVth_dVb + Alphaz * dVdsat_dVb);
            
            double T9 = 0.25 * CoxWL;
            double T4 = T9 * Alphaz;
            double T7 = T1 + T3;
            d->qdrn = -T4 * T7;
            T7 *= T9;
            double T0 = T4 * (2.0 * T5 - 2.0);        
            double T12 = T0 * dVdsat_dVb - T7 * dAlphaz_dVb;
            d->cdgb = (T0 * dVdsat_dVg - T7 * dAlphaz_dVg) * dVgs_eff_dVg;
            d->cddb = T4 * (1.0 - 2.0 * T2 - T5);
            d->cdsb = -(d->cdgb + T12 + d->cddb);
            
            d->qbulk = - (d->qgate + d->qdrn + d->qdrn);
            d->cbgb = -(d->cggb + 2.0 * d->cdgb);
            d->cbdb = -(d->cgdb + 2.0 * d->cddb);
            d->cbsb = -(d->cgsb + 2.0 * d->cdsb);
          }}
        }}
      }} // begins 1328 this 1618
      // end of else if (m->capMod == 0) line 1309 this 1621
      }else{
      double qsrc;
      assert(m->capMod != 0);
      double VbseffCV, dVbseffCV_dVb;
      {if (Vbseff < 0.0) {
        VbseffCV = Vbseff;
        dVbseffCV_dVb = 1.0;
      }else{
        VbseffCV = t->phi - Phis;
        dVbseffCV_dVb = -dPhis_dVb;
      }}
      double CoxWL = m->cox * s->weffCV * s->leffCV;
      double Vth = d->von; // possibly wrong value -- scope problem
      double Vfb = Vth - t->phi - t->k1 * sqrtPhis;
      double dVfb_dVb = dVth_dVb - t->k1 * dsqrtPhis_dVb;
      double dVfb_dVd = dVth_dVd;

      double Vgsteff;
      {if ((VgstNVt > -EXP_THRESHOLD) && (VgstNVt < EXP_THRESHOLD)) {
        assert(ExpVgst != NOT_VALID);
        ExpVgst *= ExpVgst;
        Vgsteff = n * t->vtm * log(1.0 + ExpVgst);
        dVgsteff_dVg = ExpVgst / (1.0 + ExpVgst);
        dVgsteff_dVd = -dVgsteff_dVg * (dVth_dVd + (Vgs_eff-Vth)/n * dn_dVd)
          + Vgsteff/n * dn_dVd;
        dVgsteff_dVb = -dVgsteff_dVg * (dVth_dVb + (Vgs_eff-Vth)/n * dn_dVb)
          + Vgsteff/n * dn_dVb;
        dVgsteff_dVg *= dVgs_eff_dVg;
      }else{
        Vgsteff = d->vgst;
      }}
      {if (m->capMod == 1) {
        double Cgg, Cgd, Cgb;
        {
          double Arg1 = Vgs_eff - VbseffCV - Vfb - Vgsteff;
          {if (Arg1 <= 0.0) {
            d->qgate = CoxWL * Arg1;
            Cgg = CoxWL * (dVgs_eff_dVg - dVgsteff_dVg);
            Cgd = -CoxWL * (dVfb_dVd + dVgsteff_dVd);
            Cgb = -CoxWL * (dVfb_dVb + dVbseffCV_dVb + dVgsteff_dVb);
          }else{
            double T0 = 0.5 * t->k1;
            double T1 = sqrt(T0 * T0 + Arg1);
            double T2 = CoxWL * T0 / T1;
            d->qgate = CoxWL * t->k1 * (T1 - T0);
            Cgg = T2 * (dVgs_eff_dVg - dVgsteff_dVg);
            Cgd = -T2 * (dVfb_dVd + dVgsteff_dVd);
            Cgb = -T2 * (dVfb_dVb + dVbseffCV_dVb + dVgsteff_dVb);
          }}
        }
        d->qbulk = -d->qgate;
        double Cbg = -Cgg;
        double Cbd = -Cgd;
        double Cbb = -Cgb;
        
        double AbulkCV = Abulk0 * s->abulkCVfactor;
        double dAbulkCV_dVb = s->abulkCVfactor * dAbulk0_dVb;

        double Csg, Csb, Csd;
        {
          double VdsatCV = Vgsteff / AbulkCV;
          {if (VdsatCV < d->vds) {
            double One_Third_CoxWL = CoxWL / 3.0;
            double Two_Third_CoxWL = 2.0 * One_Third_CoxWL;
            double dVdsatCV_dVg = 1.0 / AbulkCV;
            double dVdsatCV_dVb = -VdsatCV * dAbulkCV_dVb / AbulkCV;
            {
            double T0 = Vgsteff - VdsatCV / 3.0;
            double dT0_dVg = 1.0 - dVdsatCV_dVg / 3.0;
            double dT0_dVb = -dVdsatCV_dVb / 3.0;
            d->qgate += CoxWL * T0;
            double Cgg1 = CoxWL * dT0_dVg; 
            double Cgb1 = CoxWL * dT0_dVb + Cgg1 * dVgsteff_dVb;
            double Cgd1 = Cgg1 * dVgsteff_dVd;
            Cgg1 *= dVgsteff_dVg;
            Cgg += Cgg1;
            Cgb += Cgb1;
            Cgd += Cgd1;
            }
            {
            double T0 = VdsatCV - Vgsteff;
            double dT0_dVg = dVdsatCV_dVg - 1.0;
            double dT0_dVb = dVdsatCV_dVb;
            d->qbulk += One_Third_CoxWL * T0;
            double Cbg1 = One_Third_CoxWL * dT0_dVg;
            double Cbb1 = One_Third_CoxWL * dT0_dVb + Cbg1 * dVgsteff_dVb;
            double Cbd1 = Cbg1 * dVgsteff_dVd;
            Cbg1 *= dVgsteff_dVg;
            Cbg += Cbg1;
            Cbb += Cbb1;
            Cbd += Cbd1;
            }
            double T0;
            {if (m->xpart > 0.5) {
            T0 = -Two_Third_CoxWL;
            }else if (m->xpart < 0.5) {
            T0 = -0.4 * CoxWL;
            }else{
            T0 = -One_Third_CoxWL;
            }}
            qsrc = T0 * Vgsteff;
            Csg = T0 * dVgsteff_dVg;
            Csb = T0 * dVgsteff_dVb;
            Csd = T0 * dVgsteff_dVd;
            Cgb *= dVbseff_dVb;
            Cbb *= dVbseff_dVb;
            Csb *= dVbseff_dVb;
          }else{
            double T0 = AbulkCV * d->vds;
            double T1 = 12.0 * (Vgsteff - 0.5 * T0 + 1.e-20);
            double Cgg1, Cgb1, Cgd1, Cbg1, Cbb1, Cbd1;
            {
            double T2 = d->vds / T1;
            double T3 = T0 * T2;
            double dT3_dVg = -12.0 * T2 * T2 * AbulkCV;
            double dT3_dVd = 6.0 * T0 * (4.0*Vgsteff - T0) / T1 / T1 - 0.5;
            double dT3_dVb = 12.0 * T2 * T2 * dAbulkCV_dVb * Vgsteff;
            
            d->qgate += CoxWL * (Vgsteff - 0.5 * d->vds + T3);
            Cgg1 = CoxWL * (1.0 + dT3_dVg);
            Cgb1 = CoxWL * dT3_dVb + Cgg1 * dVgsteff_dVb;
            Cgd1 = CoxWL * dT3_dVd + Cgg1 * dVgsteff_dVd;
            Cgg1 *= dVgsteff_dVg;
            Cgg += Cgg1;
            Cgb += Cgb1;
            Cgd += Cgd1;
            
            d->qbulk += CoxWL * (1.0 - AbulkCV) * (0.5 * d->vds - T3);
            Cbg1 = -CoxWL * ((1.0 - AbulkCV) * dT3_dVg);
            Cbb1 = -CoxWL * ((1.0 - AbulkCV) * dT3_dVb
                         + (0.5 * d->vds - T3) * dAbulkCV_dVb)
              + Cbg1 * dVgsteff_dVb;
            Cbd1 = -CoxWL * (1.0 - AbulkCV) * dT3_dVd
              + Cbg1 * dVgsteff_dVd;
            Cbg1 *= dVgsteff_dVg;
            Cbg += Cbg1;
            Cbb += Cbb1;
            Cbd += Cbd1;
            }
            {if (m->xpart > 0.5) {
            /* 0/100 Charge petition model */
            T1 = T1 + T1;
            qsrc = -CoxWL * (0.5 * Vgsteff + 0.25 * T0 - T0 * T0 / T1);
            Csg = -CoxWL * (0.5 + 24.0 * T0 * d->vds / T1 / T1 * AbulkCV);
            Csb = -CoxWL * (0.25 * d->vds * dAbulkCV_dVb
                    - 12.0 * T0 * d->vds / T1 / T1 * (4.0 * Vgsteff - T0)
                        * dAbulkCV_dVb) + Csg * dVgsteff_dVb;
            Csd = -CoxWL * (0.25 * AbulkCV - 12.0 * AbulkCV * T0
                        / T1 / T1 * (4.0 * Vgsteff - T0))
              + Csg * dVgsteff_dVd;
            Csg *= dVgsteff_dVg;
            }else if (m->xpart < 0.5) {
            /* 40/60 Charge petition model */
            T1 = T1 / 12.0;
            double T2 = 0.5 * CoxWL / (T1 * T1);
            double T3 = Vgsteff * (2.0 * T0 * T0 / 3.0 
                               + Vgsteff * (Vgsteff - 4.0 * T0 / 3.0))
              - 2.0 * T0 * T0 * T0 / 15.0;
            qsrc = -T2 * T3;
            double T4 = 4.0 / 3.0 * Vgsteff * (Vgsteff-T0) + 0.4 * T0 * T0;
            Csg = -2.0 * qsrc / T1 
              - T2 * (Vgsteff * (3.0 * Vgsteff - 8.0 * T0 / 3.0)
                    + 2.0 * T0 * T0 / 3.0);
            Csb = (qsrc / T1 * d->vds + T2 * T4 * d->vds) * dAbulkCV_dVb
              + Csg * dVgsteff_dVb;
            Csd = (qsrc / T1 + T2 * T4) * AbulkCV + Csg * dVgsteff_dVd;
            Csg *= dVgsteff_dVg;
            }else{
            /* 50/50 Charge petition model */
            qsrc = -0.5 * (d->qgate + d->qbulk);
            Csg = -0.5 * (Cgg1 + Cbg1);
            Csb = -0.5 * (Cgb1 + Cbb1); 
            Csd = -0.5 * (Cgd1 + Cbd1); 
            }}
            Cgb *= dVbseff_dVb;
            Cbb *= dVbseff_dVb;
            Csb *= dVbseff_dVb;
          }}
        }
        d->qdrn = -(d->qgate + d->qbulk + qsrc);
        d->cggb = Cgg;
        d->cgsb = -(Cgg + Cgd + Cgb);
        d->cgdb = Cgd;
        d->cdgb = -(Cgg + Cbg + Csg);
        d->cdsb = (Cgg + Cgd + Cgb + Cbg + Cbd + Cbb + Csg + Csd + Csb);
        d->cddb = -(Cgd + Cbd + Csd);
        d->cbgb = Cbg;
        d->cbsb = -(Cbg + Cbd + Cbb);
        d->cbdb = Cbd;
      }else if (m->capMod == 2) {
        double Qac0, dQac0_dVg, dQac0_dVd, dQac0_dVb;
        double Qsub0, dQsub0_dVg, dQsub0_dVd, dQsub0_dVb;
        {
          double Vfbeff, dVfbeff_dVd, dVfbeff_dVg, dVfbeff_dVb;
          {
            const double DELTA_3 = 0.02;
            double V3 = Vfb - Vgs_eff + VbseffCV - DELTA_3;
            double T0, T2;
            {if (Vfb <= 0.0) {
            T0 = sqrt(V3 * V3 - 4.0 * DELTA_3 * Vfb);
            T2 = -DELTA_3 / T0;
            }else{
            T0 = sqrt(V3 * V3 + 4.0 * DELTA_3 * Vfb);
            T2 = DELTA_3 / T0;
            }}
            double T1 = 0.5 * (1.0 + V3 / T0);
            Vfbeff = Vfb - 0.5 * (V3 + T0);
            dVfbeff_dVd = (1.0 - T1 - T2) * dVfb_dVd;
            dVfbeff_dVg = T1 * dVgs_eff_dVg;
            dVfbeff_dVb = (1.0 - T1 - T2) * dVfb_dVb - T1 * dVbseffCV_dVb;
          }
          //double Qac0, dQac0_dVg, dQac0_dVd, dQac0_dVb;
          {
            Qac0 = CoxWL * (Vfbeff - Vfb);
            dQac0_dVg = CoxWL * dVfbeff_dVg;
            dQac0_dVd = CoxWL * (dVfbeff_dVd - dVfb_dVd);
            dQac0_dVb = CoxWL * (dVfbeff_dVb - dVfb_dVb);
          }
          //double Qsub0, dQsub0_dVg, dQsub0_dVd, dQsub0_dVb;
          {
            double T0 = 0.5 * t->k1;
            double T3 = Vgs_eff - Vfbeff - VbseffCV - Vgsteff;
            double T1, T2;
            {if (t->k1 == 0.0) {
            T1 = 0.0;
            T2 = 0.0;
            }else if (T3 < 0.0) {
            T1 = T0 + T3 / t->k1;
            T2 = CoxWL;
            }else{
            T1 = sqrt(T0 * T0 + T3);
            T2 = CoxWL * T0 / T1;
            }}
            Qsub0 = CoxWL * t->k1 * (T1 - T0);
            dQsub0_dVg = T2 * (dVgs_eff_dVg - dVfbeff_dVg - dVgsteff_dVg);
            dQsub0_dVd = -T2 * (dVfbeff_dVd + dVgsteff_dVd);
            dQsub0_dVb = -T2 * (dVfbeff_dVb +dVbseffCV_dVb +dVgsteff_dVb);
          }
        }
        double AbulkCV = Abulk0 * s->abulkCVfactor;
        double dAbulkCV_dVb = s->abulkCVfactor * dAbulk0_dVb;

        double VdseffCV, dVdseffCV_dVg, dVdseffCV_dVd, dVdseffCV_dVb;
        {
          const double DELTA_4 = 0.02;
          double VdsatCV = Vgsteff / AbulkCV;
          double V4 = VdsatCV - d->vds - DELTA_4;
          double T0 = sqrt(V4 * V4 + 4.0 * DELTA_4 * VdsatCV);
          VdseffCV = VdsatCV - 0.5 * (V4 + T0);
          double T1 = 0.5 * (1.0 + V4 / T0);
          double T2 = DELTA_4 / T0;
          double T3 = (1.0 - T1 - T2) / AbulkCV;
          dVdseffCV_dVg = T3;
          dVdseffCV_dVd = T1;
          dVdseffCV_dVb = -T3 * VdsatCV * dAbulkCV_dVb;
        }
        double T0 = AbulkCV * VdseffCV;
        double T1 = 12.0 * (Vgsteff - 0.5 * T0 + 1e-20);
        
        double Cgg1, Cgd1, Cgb1, Cbg1, Cbd1, Cbb1;
        {
          double T2 = VdseffCV / T1;
          double T3 = T0 * T2;
          double T4 = (1.0 - 12.0 * T2 * T2 * AbulkCV);
          double T5 = (6.0 * T0 * (4.0 * Vgsteff - T0) / (T1 * T1) - 0.5);
          double T6 = 12.0 * T2 * T2 * Vgsteff;
          d->qgate = CoxWL * (Vgsteff - 0.5 * VdseffCV + T3);
          Cgg1 = CoxWL * (T4 + T5 * dVdseffCV_dVg);
          Cgd1 = CoxWL * T5 * dVdseffCV_dVd + Cgg1 * dVgsteff_dVd;
          Cgb1 = CoxWL * (T5 * dVdseffCV_dVb + T6 * dAbulkCV_dVb) 
            + Cgg1 * dVgsteff_dVb;
          Cgg1 *= dVgsteff_dVg;
          
          double T7 = 1.0 - AbulkCV;
          d->qbulk = CoxWL * T7 * (0.5 * VdseffCV - T3);
          T4 = -T7 * (T4 - 1.0);
          T5 = -T7 * T5;
          T6 = -(T7 * T6 + (0.5 * VdseffCV - T3));
          Cbg1 = CoxWL * (T4 + T5 * dVdseffCV_dVg);
          Cbd1 = CoxWL * T5 * dVdseffCV_dVd + Cbg1 * dVgsteff_dVd;
          Cbb1 = CoxWL * (T5 * dVdseffCV_dVb + T6 * dAbulkCV_dVb)
            + Cbg1 * dVgsteff_dVb;
          Cbg1 *= dVgsteff_dVg;
        }
        
        double Csg, Csd, Csb;
        {if (m->xpart > 0.5) {
          /* 0/100 Charge petition model */
          T1 = T1 + T1;
          qsrc = -CoxWL * (0.5 * Vgsteff + 0.25 * T0 - T0 * T0 / T1);
          double T7 = (4.0 * Vgsteff - T0) / (T1 * T1);
          double T4 = -(0.5 + 24.0 * T0 * T0 / (T1 * T1));
          double T5 = -(0.25 * AbulkCV - 12.0 * AbulkCV * T0 * T7);
          double T6 = -(0.25 * VdseffCV - 12.0 * T0 * VdseffCV * T7);
          Csg = CoxWL * (T4 + T5 * dVdseffCV_dVg);
          Csd = CoxWL * T5 * dVdseffCV_dVd + Csg * dVgsteff_dVd;
          Csb = CoxWL * (T5 * dVdseffCV_dVb + T6 * dAbulkCV_dVb) 
            + Csg * dVgsteff_dVb;
          Csg *= dVgsteff_dVg;
        }else if (m->xpart < 0.5) {
          /* 40/60 Charge petition model */
          T1 = T1 / 12.0;
          double T2 = 0.5 * CoxWL / (T1 * T1);
          double T3 = Vgsteff * (2.0 * T0 * T0 / 3.0 + Vgsteff
                           * (Vgsteff - 4.0 * T0 / 3.0))
            - 2.0 * T0 * T0 * T0 / 15.0;
          qsrc = -T2 * T3;
          double T7 = 4.0 / 3.0 * Vgsteff * (Vgsteff - T0) + 0.4 * T0 * T0;
          double T4 = -2.0 * qsrc / T1 
            - T2 * (Vgsteff * (3.0 * Vgsteff - 8.0 * T0 / 3.0)
                  + 2.0 * T0 * T0 / 3.0);
          double T5 = (qsrc / T1 + T2 * T7) * AbulkCV;
          double T6 = (qsrc / T1 * VdseffCV + T2 * T7 * VdseffCV);
          Csg = (T4 + T5 * dVdseffCV_dVg);
          Csd = T5 * dVdseffCV_dVd + Csg * dVgsteff_dVd;
          Csb = (T5 * dVdseffCV_dVb + T6 * dAbulkCV_dVb) 
            + Csg * dVgsteff_dVb;
          Csg *= dVgsteff_dVg;
        }else{
          /* 50/50 Charge petition model */
          qsrc = -0.5 * (d->qgate + d->qbulk);
          Csg = -0.5 * (Cgg1 + Cbg1);
          Csb = -0.5 * (Cgb1 + Cbb1); 
          Csd = -0.5 * (Cgd1 + Cbd1); 
        }}

        d->qgate += Qac0 + Qsub0;
        d->qbulk -= (Qac0 + Qsub0);
        d->qdrn = -(d->qgate + d->qbulk + qsrc);
        
        double Cgg = dQac0_dVg + dQsub0_dVg + Cgg1;
        double Cgd = dQac0_dVd + dQsub0_dVd + Cgd1;
        double Cgb = dQac0_dVb + dQsub0_dVb + Cgb1;
        
        double Cbg = Cbg1 - dQac0_dVg - dQsub0_dVg;
        double Cbd = Cbd1 - dQac0_dVd - dQsub0_dVd;
        double Cbb = Cbb1 - dQac0_dVb - dQsub0_dVb;
        
        Cgb *= dVbseff_dVb;
        Cbb *= dVbseff_dVb;
        Csb *= dVbseff_dVb;
        
        d->cggb = Cgg;
        d->cgsb = -(Cgg + Cgd + Cgb);
        d->cgdb = Cgd;
        d->cdgb = -(Cgg + Cbg + Csg);
        d->cdsb = (Cgg + Cgd + Cgb + Cbg + Cbd + Cbb + Csg + Csd + Csb);
        d->cddb = -(Cgd + Cbd + Csd);
        d->cbgb = Cbg;
        d->cbsb = -(Cbg + Cbd + Cbb);
        d->cbdb = Cbd;
      }else{
        unreachable();
        d->qbulk = d->qgate = NOT_VALID;
      }}
      /* Non-quasi-static Model */
      double tconst;
      {if (m->nqsMod) {
        //  d->gtau
        double qcheq = -d->qbulk - d->qgate;
        double T0 = s->leffCV * s->leffCV;
        tconst = t->u0temp * s->elm / CoxWL / T0;
        {if (qcheq == 0.0) {
          tconst = 0.0;
        }else if (qcheq < 0.0) {
          tconst = -tconst;
        }else{
        }}
        double gtau_drift = std::abs(tconst * qcheq);
        double gtau_diff = 16.0 * t->u0temp * t->vtm / T0;
        d->gtau =  gtau_drift + gtau_diff;
        d->cqgb = -(d->cggb + d->cbgb);
        d->cqdb = -(d->cgdb + d->cbdb);
        d->cqsb = -(d->cgsb + d->cbsb);
        d->cqbb = d->cggb +d->cgdb +d->cgsb +d->cbgb +d->cbdb +d->cbsb;
        
        d->qbulk = d->qgate = d->qdrn = qsrc = 0.0;
        d->cggb = d->cgsb = d->cgdb = 0.0;
        d->cdgb = d->cdsb = d->cddb = 0.0;
        d->cbgb = d->cbsb = d->cbdb = 0.0;
#if 0
        *(ckt->CKTstate0 + d->qcheq) = qcheq;
        if (ckt->CKTmode & MODEINITTRAN)
          *(ckt->CKTstate1 + d->qcheq) = *(ckt->CKTstate0 + d->qcheq);
        error = NIintegrate(ckt, &geq, &ceq, 0.0, d->qcheq);
        if (error) return (error);
#endif
      }else{
        d->gtau = 0.0;
        d->cqgb = d->cqdb = d->cqsb = d->cqbb = 0.0;
      }}
      }}
    }
}
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/

Generated by  Doxygen 1.6.0   Back to index