// Model generated by Verilog-A Distiller.
// (c)2026 Arpad Buermen. All rights reserved.
// SPDX-License-Identifier: AGPL-3.0-or-later
// 
// SPICE flavour: ngspice
// Version      : ngspice-45-116-gf11f4594e
// Origin       : https://sourceforge.net/u/arpadbuermen/ngspice/ci/vadng-pre-46/tree/
// Date         : 2026-02-02 07:38:49 UTC
// Module name  : sp_mos2
// Model variant: sn
// 
// Notes:
//   Parameter m was replaced with builtin Verilog-A parameter $mfactor.
//
//   All parameters and output variables refer to a single parallel instance,
//   even if $mfactor>1.
//
//   Parameter off was removed. Bypassing is not available in Verilog-A. 
//   
//   Device type is set by the type parameter to NMOS (1) or PMOS (-1). 
//   
//   Output variables rs and rd were removed due to a conflict with model 
//   parameters. Their inverses are available as sourceconductance and 
//   drainconductance output variables. 
//   
//   Output variables cbd, cbs, ig, id, is, and ib were renamed to c_bd, 
//   c_bs, i_g, i_d, i_s, i_b due to conflicts with model parameters. 
//   
//   Meyer charge is not computed and therefore not available. 
//   Output variables qgs, qgd, and qgb were removed. 
//   
//   Parameters ic, icvds, icvgs, and icvbs were removed. Initial conditions 
//   should be set via the simulator's initial conditions mechanism. 
//   
//   Noise model excludes Meyer capacitance currents from the drain 
//   current that is used for computing noise. This is a SPICE3 feature. 
//   Simplified noise model also excludes the capacitive current of the 
//   B-D diode. Due to the way OpenVAF handles c(v)*ddt(v) the implicit 
//   equations do not go away even if simplified noise model is used.
//
//   Simplified noise model used (ignores capacitive current).
//   Appropriate only for small-signal noise analysis.
//
// SPICE variable mapping
//   scale -> $simparam("scale", 1)
//
// CKTcircuit member mapping
//   CKTdefaultMosAD -> ('defad', '$simparam("defad", 0)', 'const')
//   CKTdefaultMosAS -> ('defas', '$simparam("defas", 0)', 'const')
//    CKTdefaultMosL -> ('defl', '$simparam("defl", 1e-4)', 'const')
//    CKTdefaultMosW -> ('defw', '$simparam("defw", 1e-4)', 'const')
//       CKTfixLimit -> ('oldlimit', '$simparam("oldlimit", 0)', 'const')
//           CKTgmin -> ('gmin', '$simparam("gmin", 1e-12)')
//        CKTnomTemp -> ('tnom', '($simparam("tnom", 27)+273.15)')
//           CKTtemp -> (None, '$temperature')
//
// Source files and copyright information:
//   devsup.c
//     Copyright 1990 Regents of the University of California.  All rights reserved.
//     Author: 1985 Thomas L. Quarles
//     Modified: 2000 AlansFixes
//   mos2/mos2defs.h
//   mos2/mos2.c
//     Copyright 1990 Regents of the University of California.  All rights reserved.
//     Author: 1985 Thomas L. Quarles
//     Modified: 2000 AlansFIxes
//   mos2/mos2par.c
//     Copyright 1990 Regents of the University of California.  All rights reserved.
//     Author: 1985 Thomas L. Quarles
//     Modified: 2000 AlansFixes
//   mos2/mos2mpar.c
//     Copyright 1990 Regents of the University of California.  All rights reserved.
//     Author: 1985 Thomas L. Quarles
//   mos2/mos2ask.c
//     Copyright 1990 Regents of the University of California.  All rights reserved.
//     Author: 1987 Mathew Lew and Thomas L. Quarles
//     Modified: 2000 AlansFixes
//   mos2/mos2mask.c
//     Copyright 1990 Regents of the University of California.  All rights reserved.
//     Author: 1985 Thomas L. Quarles
//   mos2/mos2set.c
//   mos2/mos2temp.c
//     Copyright 1990 Regents of the University of California.  All rights reserved.
//     Author: 1985 Thomas L. Quarles
//     Modified: 2000 AlansFixes
//   mos2/mos2load.c
//     Copyright 1990 Regents of the University of California.  All rights reserved.
//     Author: 1985 Thomas L. Quarles
//     Modified: 2000 alansFixes
//   mos2/mos2noi.c
//     Copyright 1990 Regents of the University of California.  All rights reserved.
//     Author: 1987 Gary W. Ng
//     Modified: 2000 AlansFixes 

`include "constants.vams"
`include "disciplines.vams"

module sp_mos2(d, g, s, b);
  inout d, g, s, b;
  electrical d, g, s, b, d_int, s_int;

  // Instance parameters
  (* desc = "Length", type = "instance" *) parameter real l = 0;
  (* desc = "Width", type = "instance" *) parameter real w = 0;
  (* desc = "Drain area", type = "instance" *) parameter real ad = 0;
  (* desc = "Source area", type = "instance" *) parameter real as = 0;
  (* desc = "Drain perimeter", type = "instance" *) parameter real pd = 0;
  (* desc = "Source perimeter", type = "instance" *) parameter real ps = 0;
  (* desc = "Drain squares", type = "instance" *) parameter real nrd = 1;
  (* desc = "Source squares", type = "instance" *) parameter real nrs = 1;
  (* desc = "Instance operating temperature", type = "instance" *) parameter real temp = 0;
  (* desc = "Instance temperature difference", type = "instance" *) parameter real dtemp = 0.0;

  // Model parameters
  (* desc = "NMOS (1) or PMOS (-1)" *) parameter integer type = 1;
  (* desc = "Threshold voltage" *) parameter real vto = 0;
  aliasparam vt0 = vto;
  (* desc = "Transconductance parameter" *) parameter real kp = 0;
  (* desc = "Bulk threshold parameter" *) parameter real gamma = 0;
  (* desc = "Surface potential" *) parameter real phi = 0.6;
  (* desc = "Channel length modulation" *) parameter real lambda = 0;
  (* desc = "Drain ohmic resistance" *) parameter real rd = 0;
  (* desc = "Source ohmic resistance" *) parameter real rs = 0;
  (* desc = "B-D junction capacitance" *) parameter real cbd = 0;
  (* desc = "B-S junction capacitance" *) parameter real cbs = 0;
  (* desc = "Bulk junction sat. current" *) parameter real is = 1e-14;
  (* desc = "Bulk junction potential" *) parameter real pb = 0.8;
  (* desc = "Gate-source overlap cap." *) parameter real cgso = 0;
  (* desc = "Gate-drain overlap cap." *) parameter real cgdo = 0;
  (* desc = "Gate-bulk overlap cap." *) parameter real cgbo = 0;
  (* desc = "Sheet resistance" *) parameter real rsh = 0;
  (* desc = "Bottom junction cap per area" *) parameter real cj = 0;
  (* desc = "Bottom grading coefficient" *) parameter real mj = 0.5;
  (* desc = "Side junction cap per area" *) parameter real cjsw = 0;
  (* desc = "Side grading coefficient" *) parameter real mjsw = 0.33;
  (* desc = "Bulk jct. sat. current density" *) parameter real js = 0;
  (* desc = "Oxide thickness" *) parameter real tox = 1e-07;
  (* desc = "Lateral diffusion" *) parameter real ld = 0;
  (* desc = "Surface mobility" *) parameter real u0 = 600;
  aliasparam uo = u0;
  (* desc = "Forward bias jct. fit parm." *) parameter real fc = 0.5;
  (* desc = "Substrate doping" *) parameter real nsub = 0;
  (* desc = "Gate type" *) parameter integer tpg = 1;
  (* desc = "Surface state density" *) parameter real nss = 0;
  (* desc = "Width effect on threshold" *) parameter real delta = 0;
  (* desc = "Crit. field exp for mob. deg." *) parameter real uexp = 0;
  (* desc = "Crit. field for mob. degradation" *) parameter real ucrit = 10000.0;
  (* desc = "Maximum carrier drift velocity" *) parameter real vmax = 0;
  (* desc = "Junction depth" *) parameter real xj = 0;
  (* desc = "Total channel charge coeff." *) parameter real neff = 1;
  (* desc = "Fast surface state density" *) parameter real nfs = 0;
  (* desc = "Parameter measurement temperature" *) parameter real tnom = 0;
  (* desc = "Flicker noise coefficient" *) parameter real kf = 0;
  (* desc = "Flicker noise exponent" *) parameter real af = 1;
  (* desc = "Noise model selection" *) parameter integer nlev = 2;
  (* desc = "Channel shot noise coefficient" *) parameter real gdsnoi = 1;

  // Instance output variables
  (* desc = "Drain current" *) real i_d;
  (* desc = "B-D junction current" *) real ibd;
  (* desc = "B-S junction current" *) real ibs;
  (* desc = "Source current " *) real i_s;
  (* desc = "Gate current " *) real i_g;
  (* desc = "Bulk current " *) real i_b;
  (* desc = "Gate-Source voltage" *) real vgs;
  (* desc = "Drain-Source voltage" *) real vds;
  (* desc = "Bulk-Source voltage" *) real vbs;
  (* desc = "Bulk-Drain voltage" *) real vbd;
  (* desc = " " *) real von;
  (* desc = "Saturation drain voltage" *) real vdsat;
  (* desc = "Critical source voltage" *) real sourcevcrit;
  (* desc = "Critical drain voltage" *) real drainvcrit;
  (* desc = "Source conductance" *) real sourceconductance;
  (* desc = "Drain conductance" *) real drainconductance;
  (* desc = "Transconductance" *) real gm;
  (* desc = "Drain-Source conductance" *) real gds;
  (* desc = "Bulk-Source transconductance" *) real gmb;
  (* desc = "Bulk-Drain conductance" *) real gbd;
  (* desc = "Bulk-Source conductance" *) real gbs;
  (* desc = "Bulk-Drain capacitance" *) real c_bd;
  (* desc = "Bulk-Source capacitance" *) real c_bs;
  (* desc = "Gate-Source capacitance" *) real cgs;
  (* desc = "Gate-Drain capacitance" *) real cgd;
  (* desc = "Gate-Bulk capacitance" *) real cgb;
  (* desc = "Zero-Bias B-D junction capacitance" *) real cbd0;
  (* desc = " " *) real cbdsw0;
  (* desc = "Zero-Bias B-S junction capacitance" *) real cbs0;
  (* desc = " " *) real cbssw0;
  (* desc = "Capacitance due to gate-source charge storage" *) real cqgs;
  (* desc = "Capacitance due to gate-drain charge storage" *) real cqgd;
  (* desc = "Capacitance due to gate-bulk charge storage" *) real cqgb;
  (* desc = "Capacitance due to bulk-drain charge storage" *) real cqbd;
  (* desc = "Capacitance due to bulk-source charge storage" *) real cqbs;
  (* desc = "Bulk-Drain charge storage" *) real qbd;
  (* desc = "Bulk-Source charge storage" *) real qbs;

  // Model output variables

  // Instance fields
  integer MOS2mode = 0;
  real MOS2cgs = 0;
  real MOS2cgd = 0;
  real MOS2cgb = 0;
  real MOS2temp = 0;
  real MOS2tTransconductance = 0;
  real MOS2tSurfMob = 0;
  real MOS2tPhi = 0;
  real MOS2tVto = 0;
  real MOS2tSatCur = 0;
  real MOS2tSatCurDens = 0;
  real MOS2tCbd = 0;
  real MOS2tCbs = 0;
  real MOS2tCj = 0;
  real MOS2tCjsw = 0;
  real MOS2tBulkPot = 0;
  real MOS2tDepCap = 0;
  real MOS2tVbi = 0;
  real MOS2l = 0;
  real MOS2w = 0;
  real MOS2drainArea = 0;
  real MOS2sourceArea = 0;
  real MOS2drainPerimiter = 0;
  real MOS2sourcePerimiter = 0;
  real MOS2sourceConductance = 0;
  real MOS2drainConductance = 0;
  real MOS2von = 0;
  real MOS2vdsat = 0;
  real MOS2sourceVcrit = 0;
  real MOS2drainVcrit = 0;
  real MOS2cd = 0;
  real MOS2cbs = 0;
  real MOS2cbd = 0;
  real MOS2gmbs = 0;
  real MOS2gm = 0;
  real MOS2gds = 0;
  real MOS2gbd = 0;
  real MOS2gbs = 0;
  real MOS2capbd = 0;
  real MOS2capbs = 0;
  real MOS2Cbd = 0;
  real MOS2Cbdsw = 0;
  real MOS2Cbs = 0;
  real MOS2Cbssw = 0;
  real MOS2f2d = 0;
  real MOS2f3d = 0;
  real MOS2f4d = 0;
  real MOS2f2s = 0;
  real MOS2f3s = 0;
  real MOS2f4s = 0;
  real MOS2dcCurrent = 0;
  real meyer_scale = 0;
  real meyer_unscale = 0;
  integer MOS2lGiven = 0;
  integer MOS2wGiven = 0;
  integer MOS2drainAreaGiven = 0;
  integer MOS2sourceAreaGiven = 0;
  integer MOS2drainPerimiterGiven = 0;
  integer MOS2sourcePerimiterGiven = 0;
  integer MOS2tempGiven = 0;

  // Model fields
  real MOS2tnom = 0;
  real MOS2transconductance = 0;
  real MOS2oxideCapFactor = 0;
  real MOS2vt0 = 0;
  real MOS2bulkCapFactor = 0;
  real MOS2phi = 0;
  real MOS2gamma = 0;
  real MOS2substrateDoping = 0;
  real MOS2xd = 0;
  integer MOS2type = 0;
  integer MOS2tnomGiven = 0;

  // States
  real MOS2vbd = 0;
  real MOS2vbs = 0;
  real MOS2vgs = 0;
  real MOS2vds = 0;
  real MOS2capgs = 0;
  real MOS2qgs = 0;
  real MOS2cqgs = 0;
  real MOS2capgd = 0;
  real MOS2qgd = 0;
  real MOS2cqgd = 0;
  real MOS2capgb = 0;
  real MOS2qgb = 0;
  real MOS2cqgb = 0;
  real MOS2qbd = 0;
  real MOS2cqbd = 0;
  real MOS2qbs = 0;
  real MOS2cqbs = 0;

  // Model parameter access locals

  // Instance parameter access locals

  // Device setup locals

  // Device temperature dependence locals
  real egfet = 0;
  real wkfngs = 0;
  real wkfng = 0;
  real fermig = 0;
  real fermis = 0;
  real vfb = 0;
  real fact1 = 0;
  real fact2 = 0;
  real vt = 0;
  real vtnom = 0;
  real kt = 0;
  real kt1 = 0;
  real egfet1 = 0;
  real arg1 = 0;
  real pbfact = 0;
  real pbfact1 = 0;
  real ratio = 0;
  real ratio4 = 0;
  real phio = 0;
  real pbo = 0;
  real gmaold = 0;
  real gmanew = 0;
  real capfact = 0;
  real czbd = 0;
  real czbdsw = 0;
  real czbs = 0;
  real czbssw = 0;
  real arg = 0;
  real sarg = 0;
  real sargsw = 0;

  // Device load locals
  integer limited = 0;
  real Beta = 0;
  real DrainSatCur = 0;
  real EffectiveLength = 0;
  real GateBulkOverlapCap = 0;
  real GateDrainOverlapCap = 0;
  real GateSourceOverlapCap = 0;
  real OxideCap = 0;
  real unscaledOxideCap = 0;
  real SourceSatCur = 0;
  real load_arg = 0;
  real cdrain = 0;
  real cdreq = 0;
  real ceq = 0;
  real ceqbd = 0;
  real ceqbs = 0;
  real ceqgb = 0;
  real ceqgd = 0;
  real ceqgs = 0;
  real evbd = 0;
  real evbs = 0;
  real gcgb = 0;
  real gcgd = 0;
  real gcgs = 0;
  real geq = 0;
  real load_sarg = 0;
  real load_sargsw = 0;
  real load_vbd = 0;
  real load_vbs = 0;
  real load_vds = 0;
  real load_vdsat = 0;
  real vgb = 0;
  real vgd = 0;
  real vgdo = 0;
  real load_vgs = 0;
  real load_von = 0;
  real load_vt = 0;
  real xfact = 0;
  real capgs = 0;
  real capgd = 0;
  real capgb = 0;
  integer xnrm = 0;
  integer xrev = 0;
  real load_arg1 = 0;
  real sarg1 = 0;
  real a4_0 = 0;
  real a4_1 = 0;
  real a4_2 = 0;
  real a4_3 = 0;
  real b4_0 = 0;
  real b4_1 = 0;
  real b4_2 = 0;
  real b4_3 = 0;
  real x4_0 = 0;
  real x4_1 = 0;
  real x4_2 = 0;
  real x4_3 = 0;
  real x4_4 = 0;
  real x4_5 = 0;
  real x4_6 = 0;
  real x4_7 = 0;
  real poly4_0 = 0;
  real poly4_1 = 0;
  real poly4_2 = 0;
  real poly4_3 = 0;
  real poly4_4 = 0;
  real poly4_5 = 0;
  real poly4_6 = 0;
  real poly4_7 = 0;
  real beta1 = 0;
  real dsrgdb = 0;
  real d2sdb2 = 0;
  real sphi = 0;
  real sphi3 = 0;
  real barg = 0;
  real d2bdb2 = 0;
  real factor = 0;
  real dbrgdb = 0;
  real eta = 0;
  real vbin = 0;
  real argd = 0;
  real args = 0;
  real argss = 0;
  real argsd = 0;
  real argxs = 0;
  real argxd = 0;
  real daddb2 = 0;
  real dasdb2 = 0;
  real dbargd = 0;
  real dbargs = 0;
  real dbxwd = 0;
  real dbxws = 0;
  real dgddb2 = 0;
  real dgddvb = 0;
  real dgdvds = 0;
  real gamasd = 0;
  real xwd = 0;
  real xws = 0;
  real ddxwd = 0;
  real gammad = 0;
  real vth = 0;
  real cfs = 0;
  real cdonco = 0;
  real xn = 0;
  real argg = 0;
  real vgst = 0;
  real sarg3 = 0;
  real sbiarg = 0;
  real dgdvbs = 0;
  real body = 0;
  real gdbdv = 0;
  real dodvbs = 0;
  real dodvds = 0;
  real dxndvd = 0;
  real dxndvb = 0;
  real udenom = 0;
  real dudvgs = 0;
  real dudvds = 0;
  real dudvbs = 0;
  real gammd2 = 0;
  real argv = 0;
  real vgsx = 0;
  real ufact = 0;
  real ueff = 0;
  real dsdvgs = 0;
  real dsdvbs = 0;
  real a1 = 0;
  real a3 = 0;
  real a = 0;
  real b1 = 0;
  real b3 = 0;
  real load_b = 0;
  real c1 = 0;
  real c = 0;
  real d1 = 0;
  real fi = 0;
  real p0 = 0;
  real p2 = 0;
  real p3 = 0;
  real p4 = 0;
  real p = 0;
  real r3 = 0;
  real r = 0;
  real ro = 0;
  real s2 = 0;
  real load_s = 0;
  real v1 = 0;
  real v2 = 0;
  real xv = 0;
  real y3 = 0;
  real delta4 = 0;
  real xvalid = 0;
  real bsarg = 0;
  real dbsrdb = 0;
  real bodys = 0;
  real gdbdvs = 0;
  real sargv = 0;
  real xlfact = 0;
  real dldsat = 0;
  real xdv = 0;
  real xlv = 0;
  real vqchan = 0;
  real dqdsat = 0;
  real vl = 0;
  real dfundg = 0;
  real dfunds = 0;
  real dfundb = 0;
  real xls = 0;
  real dldvgs = 0;
  real dldvds = 0;
  real dldvbs = 0;
  real dfact = 0;
  real clfact = 0;
  real xleff = 0;
  real deltal = 0;
  real xwb = 0;
  real vdson = 0;
  real cdson = 0;
  real didvds = 0;
  real gdson = 0;
  real gmw = 0;
  real gbson = 0;
  real expg = 0;
  real xld = 0;
  real xlamda = 0;
  real lvbs = 0;
  real lvds = 0;
  real lvgs = 0;
  real phiMinVbs = 0;
  real tmp = 0;
  integer iknt = 0;
  integer jknt = 0;
  integer i = 0;
  integer j = 0;
  real sig1_0 = 1;
  real sig1_1 = -1;
  real sig1_2 = 1;
  real sig1_3 = -1;
  real sig2_0 = 1;
  real sig2_1 = 1;
  real sig2_2 = -1;
  real sig2_3 = -1;
  integer line1050 = 0;
  integer line400 = 1;
  integer line410 = 2;
  integer line500 = 3;
  integer line900 = 4;
  integer doneval = 5;
  integer goto_destination = -1;

  // Device noise locals
  real noise_vgs = 0;
  real noise_vds = 0;
  real noise_vgd = 0;
  real noise_vgst = 0;
  real alpha = 0;
  real beta = 0;
  real Sid = 0;
  real noise_dtemp = 0;
  real flicker_psd = 0;
  real flicker_exp = 1;

  // Global constants
  real CONSTroot2 = sqrt(2.0);
  real CONSTvt0 = 1.38064852e-23*(27+273.15)/1.6021766208e-19;
  real CONSTKoverQ = 1.38064852e-23/1.6021766208e-19;
  real CONSTe = exp(1.0);
  real VACONSTcharge = 1.6021766208e-19;
  real VACONSTboltz = 1.38064852e-23;
  real defad = $simparam("defad", 0);
  real defas = $simparam("defas", 0);
  real defl = $simparam("defl", 1e-4);
  real defw = $simparam("defw", 1e-4);
  real oldlimit = $simparam("oldlimit", 0);
  real VACONST_tnom = ($simparam("tnom", 27)+273.15);
  real cpscale = $simparam("scale", 1);

  // Evaluation constants
  real lc_gmin = 0;

  // Translated function definitions
  analog function real DEVqmeyer;
    input vgs, vgd, vgb, von, _vdsat;
    inout capgs, capgd, capgb;
    input phi, cox;
    real vgs;
    real vgd;
    real vgb;
    real von;
    real _vdsat;
    real capgs;
    real capgd;
    real capgb;
    real phi;
    real cox;
    real vdsat;
    real vds;
    real vddif;
    real vddif1;
    real vddif2;
    real vgst;
    begin
      DEVqmeyer = 0;
      vdsat = _vdsat;
      vgst = vgs-von;
      vdsat = ((vdsat>0.025) ? (vdsat) : (0.025));
      if (vgst<=-phi) begin
        capgb = cox/2;
        capgs = 0;
        capgd = 0;
      end else if (vgst<=-phi/2) begin
        capgb = -vgst*cox/(2*phi);
        capgs = 0;
        capgd = 0;
      end else if (vgst<=0) begin
        capgb = -vgst*cox/(2*phi);
        capgs = vgst*cox/(1.5*phi)+cox/3;
        vds = vgs-vgd;
        if (vds>=vdsat) begin
          capgd = 0;
        end else begin
          vddif = 2.0*vdsat-vds;
          vddif1 = vdsat-vds;
          vddif2 = vddif*vddif;
          capgd = capgs*(1.0-vdsat*vdsat/vddif2);
          capgs = capgs*(1.0-vddif1*vddif1/vddif2);
        end
      end else begin
        vds = vgs-vgd;
        vdsat = ((vdsat>0.025) ? (vdsat) : (0.025));
        if (vdsat<=vds) begin
          capgs = cox/3;
          capgd = 0;
          capgb = 0;
        end else begin
          vddif = 2.0*vdsat-vds;
          vddif1 = vdsat-vds;
          vddif2 = vddif*vddif;
          capgd = cox*(1.0-vdsat*vdsat/vddif2)/3;
          capgs = cox*(1.0-vddif1*vddif1/vddif2)/3;
          capgb = 0;
        end
      end
    end
  endfunction

  // Verbatim analog functions
  analog function integer initialize_limiting;
    integer il;
    begin
      il = $simparam("iniLim", -1);
      if (il<0) begin
        // iniLim not available, use heuristics
        initialize_limiting = 
          $simparam("iteration", 10)==1 && 
          (analysis("dc") || analysis("static")) && 
          analysis("nodeset");
      end else begin
        initialize_limiting = il;
      end
    end
  endfunction

  analog function real DEVlimitOldGet;
    input vnew, vold;
    real vnew, vold, result;
    begin
      DEVlimitOldGet = vold;
    end
  endfunction

  analog function real DEVlimitNewSet;
    input vnew, vold, new_value, limited;
    real vnew, vold, new_value;
    integer limited;
    begin
      DEVlimitNewSet = new_value;
      if (limited) begin
        $discontinuity(-1);
      end
    end
  endfunction

  analog function real DEVpnjlim;
    input vnew, vold, vt, vcrit;
    inout limiting_applied;
    real vnew, vold, vt, vcrit, limited, arg;
    integer limiting_applied;
    
    begin
      limited = vnew;
      if ((vnew > vcrit) && (abs(vnew - vold) > (vt + vt))) begin
        // Positive vnew above vcrit, change greater than 2*vt
        if(vold > 0) begin
          arg = (vnew - vold) / vt;
          if(arg > 0) begin
            limited = vold + vt * ln(1+arg);
            limiting_applied = 1;
          end else begin
            limited = vold - vt * ln(1-arg);
            limiting_applied = 1;
          end
        end else begin
          limited = vt *ln(vnew/vt);
          limiting_applied = 1;
        end
      end else if (vnew<0) begin
        // Negative vnew, use ngspice limiting
        if (vold > 0) begin
          arg = -1*vold-1; 
        end else begin
          arg = 2*vold-1;
        end
        if (vnew<arg) begin
          limited = arg;
          limiting_applied = 1;
        end
      end 
      DEVpnjlim = limited;
    end
  endfunction

  analog function real DEVfetlim;
    input vnew, vold, vto;
    inout limiting_applied;
    real vnew, vold, vto, vlimited;
    real vtsthi, vtstlo, vtox, delv, vtemp;
    integer limiting_applied;  
    begin
      vlimited = vnew;
      vtsthi = abs(2*(vold-vto))+2;
      // vtstlo = vtsthi/2 +2;
      vtstlo = abs(vold-vto)+1;
      vtox = vto + 3.5;
      delv = vnew-vold;
      
      if (vold >= vto) begin
        if(vold >= vtox) begin
          if(delv <= 0) begin
            /* going off */
            if(vlimited >= vtox) begin
                if(-delv > vtstlo) begin
                vlimited =  vold - vtstlo;
                end
            end else begin
              vlimited = max(vnew,vto+2);
            end
          end else begin
            /* staying on */
            if(delv >= vtsthi) begin
              vlimited = vold + vtsthi;
            end
          end
        end else begin
          /* middle region */
          if(delv <= 0) begin
            /* decreasing */
            vlimited = max(vnew,vto-0.5);
          end else begin
            /* increasing */
            vlimited = min(vnew,vto+4);
          end
        end
      end else begin
        /* off */
        if(delv <= 0) begin
          if(-delv >vtsthi) begin
            vlimited = vold - vtsthi;
          end
        end else begin
          vtemp = vto + 0.5;
          if(vnew <= vtemp) begin
            if(delv > vtstlo) begin
              vlimited = vold + vtstlo;
            end
          end else begin
            vlimited = vtemp;
          end
        end
      end
      if (vlimited!=vnew) begin
        limiting_applied = 1;
      end
      DEVfetlim = vlimited;
    end
  endfunction

  analog function real DEVlimvds;
    input vnew, vold;
    inout limiting_applied;
    real vnew, vold, vlimited;
    integer limiting_applied;
    begin
      if (vold >= 3.5) begin
        if(vnew > vold) begin
          vlimited = min(vnew, (3 * vold) + 2);
        end else begin
          if (vnew < 3.5) begin
            vlimited = max(vnew, 2);
          end else begin
            vlimited = vnew;
          end
        end
      end else begin
        if(vnew > vold) begin
          vlimited = min(vnew, 4);
        end else begin
          vlimited = max(vnew, -0.5);
        end
      end
      DEVlimvds = vlimited;
      if (vlimited!=vnew) begin
        limiting_applied = 1;
      end
    end
  endfunction

  // Analog block
  analog begin

    // Instance field assignments
    if ($param_given(l)) begin
      MOS2l = l*cpscale;
      MOS2lGiven = 1;
    end
    if ($param_given(w)) begin
      MOS2w = w*cpscale;
      MOS2wGiven = 1;
    end
    if ($param_given(ad)) begin
      MOS2drainArea = ad*cpscale*cpscale;
      MOS2drainAreaGiven = 1;
    end
    if ($param_given(as)) begin
      MOS2sourceArea = as*cpscale*cpscale;
      MOS2sourceAreaGiven = 1;
    end
    if ($param_given(pd)) begin
      MOS2drainPerimiter = pd*cpscale;
      MOS2drainPerimiterGiven = 1;
    end
    if ($param_given(ps)) begin
      MOS2sourcePerimiter = ps*cpscale;
      MOS2sourcePerimiterGiven = 1;
    end
    if ($param_given(temp)) begin
      MOS2temp = temp+273.15;
      MOS2tempGiven = 1;
    end

    // Model field assignments
    MOS2type = type;
    if ($param_given(vto)) begin
      MOS2vt0 = vto;
    end
    if ($param_given(kp)) begin
      MOS2transconductance = kp;
    end
    if ($param_given(gamma)) begin
      MOS2gamma = gamma;
    end
    if ($param_given(phi)) begin
      MOS2phi = phi;
    end
    if ($param_given(cj)) begin
      MOS2bulkCapFactor = cj;
    end
    if ($param_given(nsub)) begin
      MOS2substrateDoping = nsub;
    end
    if ($param_given(tnom)) begin
      MOS2tnom = tnom+273.15;
      MOS2tnomGiven = 1;
    end

    // Evaluation constants
    lc_gmin = $simparam("gmin", 1e-12);

    // Tag: setup, function: MOS2setup
    if (!$param_given(vto)) begin
      MOS2vt0 = 0;
    end
    if (!$param_given(phi)) begin
      MOS2phi = 0.6;
    end
    if (!$param_given(gamma)) begin
      MOS2gamma = 0;
    end
    if (!MOS2drainPerimiterGiven) begin
      MOS2drainPerimiter = 0;
    end
    if (!MOS2sourcePerimiterGiven) begin
      MOS2sourcePerimiter = 0;
    end
    MOS2vdsat = 0;
    if (!(rd!=0||nrd!=0&&rsh!=0)) begin
      V(d_int, d) <+ 0;
    end
    if (!(rs!=0||nrs!=0&&rsh!=0)) begin
      V(s_int, s) <+ 0;
    end

    // Tag: temp, function: MOS2temp
    if (!MOS2tnomGiven) begin
      MOS2tnom = VACONST_tnom;
    end
    fact1 = MOS2tnom/(27.0+273.15);
    vtnom = MOS2tnom*CONSTKoverQ;
    kt1 = 1.38064852e-23*MOS2tnom;
    egfet1 = 1.16-7.02e-4*MOS2tnom*MOS2tnom/(MOS2tnom+1108);
    arg1 = -egfet1/(kt1+kt1)+1.1150877/(1.38064852e-23*(27.0+273.15+(27.0+273.15)));
    pbfact1 = -2*vtnom*(1.5*ln(fact1)+1.6021766208e-19*arg1);
    if (MOS2phi<=0.0) begin
      $fatal(0, "Phi is not positive.");
    end
    MOS2oxideCapFactor = 3.9*8.854214871e-12/tox;
    if (!$param_given(kp)) begin
      MOS2transconductance = u0*1e-4*MOS2oxideCapFactor;
    end
    if ($param_given(nsub)) begin
      if (MOS2substrateDoping*1e6>1.45e16) begin
        if (!$param_given(phi)) begin
          MOS2phi = 2*vtnom*ln(MOS2substrateDoping*1e6/1.45e16);
          MOS2phi = ((0.1>MOS2phi) ? (0.1) : (MOS2phi));
        end
        fermis = MOS2type*0.5*MOS2phi;
        wkfng = 3.2;
        if (tpg!=0) begin
          fermig = MOS2type*tpg*0.5*egfet1;
          wkfng = 3.25+0.5*egfet1-fermig;
        end
        wkfngs = wkfng-(3.25+0.5*egfet1+fermis);
        if (!$param_given(gamma)) begin
          MOS2gamma = sqrt(2*11.70*8.854214871e-12*1.6021766208e-19*MOS2substrateDoping*1e6)/MOS2oxideCapFactor;
        end
        if (!$param_given(vto)) begin
          vfb = wkfngs-nss*1e4*1.6021766208e-19/MOS2oxideCapFactor;
          MOS2vt0 = vfb+MOS2type*(MOS2gamma*sqrt(MOS2phi)+MOS2phi);
        end else begin
          vfb = MOS2vt0-MOS2type*(MOS2gamma*sqrt(MOS2phi)+MOS2phi);
        end
        MOS2xd = sqrt((11.7*8.854214871e-12+11.7*8.854214871e-12)/(1.6021766208e-19*MOS2substrateDoping*1e6));
      end else begin
        MOS2substrateDoping = 0;
        $fatal(0, "Nsub < Ni");
      end
    end
    if (!$param_given(cj)) begin
      MOS2bulkCapFactor = sqrt(11.7*8.854214871e-12*1.6021766208e-19*MOS2substrateDoping*1e6/(2*pb));
    end
    if (!MOS2tempGiven) begin
      MOS2temp = $temperature+dtemp;
    end
    MOS2mode = 1;
    MOS2von = 0;
    vt = MOS2temp*CONSTKoverQ;
    ratio = MOS2temp/MOS2tnom;
    fact2 = MOS2temp/(27.0+273.15);
    kt = MOS2temp*1.38064852e-23;
    egfet = 1.16-7.02e-4*MOS2temp*MOS2temp/(MOS2temp+1108);
    arg = -egfet/(kt+kt)+1.1150877/(1.38064852e-23*(27.0+273.15+(27.0+273.15)));
    pbfact = -2*vt*(1.5*ln(fact2)+1.6021766208e-19*arg);
    if (!MOS2drainAreaGiven) begin
      MOS2drainArea = defad;
    end
    if (!MOS2lGiven) begin
      MOS2l = defl;
    end
    if (!MOS2sourceAreaGiven) begin
      MOS2sourceArea = defas;
    end
    if (!MOS2wGiven) begin
      MOS2w = defw;
    end
    if ($param_given(rd)) begin
      if (rd!=0) begin
        MOS2drainConductance = 1.0/rd;
      end else begin
        MOS2drainConductance = 0;
      end
    end else if ($param_given(rsh)) begin
      if (rsh!=0&&nrd!=0) begin
        MOS2drainConductance = 1.0/(rsh*nrd);
      end else begin
        MOS2drainConductance = 0;
      end
    end else begin
      MOS2drainConductance = 0;
    end
    if ($param_given(rs)) begin
      if (rs!=0) begin
        MOS2sourceConductance = 1.0/rs;
      end else begin
        MOS2sourceConductance = 0;
      end
    end else if ($param_given(rsh)) begin
      if (rsh!=0&&nrs!=0) begin
        MOS2sourceConductance = 1.0/(rsh*nrs);
      end else begin
        MOS2sourceConductance = 0;
      end
    end else begin
      MOS2sourceConductance = 0;
    end
    if (MOS2l-2*ld<=0) begin
      $warning("effective channel length less than zero");
    end
    ratio4 = ratio*sqrt(ratio);
    MOS2tTransconductance = MOS2transconductance/ratio4;
    MOS2tSurfMob = u0/ratio4;
    phio = (MOS2phi-pbfact1)/fact1;
    MOS2tPhi = fact2*phio+pbfact;
    MOS2tVbi = MOS2vt0-MOS2type*(MOS2gamma*sqrt(MOS2phi))+0.5*(egfet1-egfet)+MOS2type*0.5*(MOS2tPhi-MOS2phi);
    MOS2tVto = MOS2tVbi+MOS2type*MOS2gamma*sqrt(MOS2tPhi);
    MOS2tSatCur = is*exp(-egfet/vt+egfet1/vtnom);
    MOS2tSatCurDens = js*exp(-egfet/vt+egfet1/vtnom);
    pbo = (pb-pbfact1)/fact1;
    gmaold = (pb-pbo)/pbo;
    capfact = 1/(1+mj*(4e-4*(MOS2tnom-(27.0+273.15))-gmaold));
    MOS2tCbd = cbd*capfact;
    MOS2tCbs = cbs*capfact;
    MOS2tCj = MOS2bulkCapFactor*capfact;
    capfact = 1/(1+mjsw*(4e-4*(MOS2tnom-(27.0+273.15))-gmaold));
    MOS2tCjsw = cjsw*capfact;
    MOS2tBulkPot = fact2*pbo+pbfact;
    gmanew = (MOS2tBulkPot-pbo)/pbo;
    capfact = 1+mj*(4e-4*(MOS2temp-(27.0+273.15))-gmanew);
    MOS2tCbd = MOS2tCbd * (capfact);
    MOS2tCbs = MOS2tCbs * (capfact);
    MOS2tCj = MOS2tCj * (capfact);
    capfact = 1+mjsw*(4e-4*(MOS2temp-(27.0+273.15))-gmanew);
    MOS2tCjsw = MOS2tCjsw * (capfact);
    MOS2tDepCap = fc*MOS2tBulkPot;
    if (MOS2tSatCurDens==0||MOS2drainArea==0||MOS2sourceArea==0) begin
      MOS2drainVcrit = vt*ln(vt/(CONSTroot2*MOS2tSatCur));
      MOS2sourceVcrit = MOS2drainVcrit;
    end else begin
      MOS2drainVcrit = vt*ln(vt/(CONSTroot2*MOS2tSatCurDens*MOS2drainArea));
      MOS2sourceVcrit = vt*ln(vt/(CONSTroot2*MOS2tSatCurDens*MOS2sourceArea));
    end
    if ($param_given(cbd)) begin
      czbd = MOS2tCbd;
    end else if ($param_given(cj)) begin
      czbd = MOS2tCj*MOS2drainArea;
    end else begin
      czbd = 0;
    end
    if ($param_given(cjsw)) begin
      czbdsw = MOS2tCjsw*MOS2drainPerimiter;
    end else begin
      czbdsw = 0;
    end
    arg = 1-fc;
    sarg = exp(-mj*ln(arg));
    sargsw = exp(-mjsw*ln(arg));
    MOS2Cbd = czbd;
    MOS2Cbdsw = czbdsw;
    MOS2f2d = czbd*(1-fc*(1+mj))*sarg/arg+czbdsw*(1-fc*(1+mjsw))*sargsw/arg;
    MOS2f3d = czbd*mj*sarg/arg/MOS2tBulkPot+czbdsw*mjsw*sargsw/arg/MOS2tBulkPot;
    MOS2f4d = czbd*MOS2tBulkPot*(1-arg*sarg)/(1-mj)+czbdsw*MOS2tBulkPot*(1-arg*sargsw)/(1-mjsw)-MOS2f3d/2*(MOS2tDepCap*MOS2tDepCap)-MOS2tDepCap*MOS2f2d;
    if ($param_given(cbs)) begin
      czbs = MOS2tCbs;
    end else if ($param_given(cj)) begin
      czbs = MOS2tCj*MOS2sourceArea;
    end else begin
      czbs = 0;
    end
    if ($param_given(cjsw)) begin
      czbssw = MOS2tCjsw*MOS2sourcePerimiter;
    end else begin
      czbssw = 0;
    end
    arg = 1-fc;
    sarg = exp(-mj*ln(arg));
    sargsw = exp(-mjsw*ln(arg));
    MOS2Cbs = czbs;
    MOS2Cbssw = czbssw;
    MOS2f2s = czbs*(1-fc*(1+mj))*sarg/arg+czbssw*(1-fc*(1+mjsw))*sargsw/arg;
    MOS2f3s = czbs*mj*sarg/arg/MOS2tBulkPot+czbssw*mjsw*sargsw/arg/MOS2tBulkPot;
    MOS2f4s = czbs*MOS2tBulkPot*(1-arg*sarg)/(1-mj)+czbssw*MOS2tBulkPot*(1-arg*sargsw)/(1-mjsw)-MOS2f3s/2*(MOS2tDepCap*MOS2tDepCap)-MOS2tDepCap*MOS2f2s;

    // Tag: load, function: MOS2load
    cdrain = 0.0;
    xfact = 0.0;
    capgs = 0.0;
    capgd = 0.0;
    capgb = 0.0;
    load_vt = CONSTKoverQ*MOS2temp;
    EffectiveLength = MOS2l-2*ld;
    if (MOS2tSatCurDens==0||MOS2drainArea==0||MOS2sourceArea==0) begin
      DrainSatCur = MOS2tSatCur;
      SourceSatCur = MOS2tSatCur;
    end else begin
      DrainSatCur = MOS2tSatCurDens*MOS2drainArea;
      SourceSatCur = MOS2tSatCurDens*MOS2sourceArea;
    end
    GateSourceOverlapCap = cgso*MOS2w;
    GateDrainOverlapCap = cgdo*MOS2w;
    GateBulkOverlapCap = cgbo*EffectiveLength;
    Beta = MOS2tTransconductance*MOS2w/EffectiveLength;
    unscaledOxideCap = MOS2oxideCapFactor*EffectiveLength*MOS2w;
    OxideCap = unscaledOxideCap;
    limited = 0;
    /* Previous vgs, vds, vbs, and vbd */
    MOS2vgs = MOS2type * $limit(V(g, s_int), DEVlimitOldGet);
    MOS2vds = MOS2type * $limit(V(d_int, s_int), DEVlimitOldGet);
    MOS2vbs = MOS2type * $limit(V(b, s_int), DEVlimitOldGet);
    MOS2vbd = MOS2type * $limit(V(b, d_int), DEVlimitOldGet);
    
    /* Simplified von value used in limiting, does not require introduction of a hidden state */
    factor = 0.125*delta*2.0*3.14159265358979323846*(11.7 * 8.854214871e-12)/ 
                        OxideCap*EffectiveLength;
    vbin = MOS2tVbi*MOS2type+factor*phiMinVbs;
    load_von = vbin;
    MOS2von = MOS2type*von;
    load_vbs = MOS2type*V(b, s_int);
    load_vgs = MOS2type*V(g, s_int);
    load_vds = MOS2type*V(d_int, s_int);
    load_vbd = load_vbs-load_vds;
    vgd = load_vgs-load_vds;
    vgdo = MOS2vgs-MOS2vds;
    load_von = MOS2type*MOS2von;
    if (MOS2vds>=0) begin
      load_vgs = DEVfetlim(load_vgs, MOS2vgs, load_von, limited);
      load_vds = load_vgs-vgd;
      load_vds = DEVlimvds(load_vds, MOS2vds, limited);
      vgd = load_vgs-load_vds;
    end else begin
      vgd = DEVfetlim(vgd, vgdo, load_von, limited);
      load_vds = load_vgs-vgd;
      if (!oldlimit) begin
        load_vds = -DEVlimvds(-load_vds, -MOS2vds, limited);
      end
      load_vgs = vgd+load_vds;
    end
    if (load_vds>=0) begin
      load_vbs = DEVpnjlim(load_vbs, MOS2vbs, load_vt, MOS2sourceVcrit, limited);
      load_vbd = load_vbs-load_vds;
    end else begin
      load_vbd = DEVpnjlim(load_vbd, MOS2vbd, load_vt, MOS2drainVcrit, limited);
      load_vbs = load_vbd+load_vds;
    end
    if (initialize_limiting()) begin
        load_vbs = -1;
        load_vgs = MOS2type * MOS2tVto;
        load_vds = 0;
        load_vbd = load_vbs - load_vds;
    end
    
    load_vgs = MOS2type * $limit(V(g, s_int), DEVlimitNewSet, MOS2type * load_vgs, limited);
    load_vds = MOS2type * $limit(V(d_int, s_int), DEVlimitNewSet, MOS2type * load_vds, limited);
    load_vbs = MOS2type * $limit(V(b, s_int), DEVlimitNewSet, MOS2type * load_vbs, limited);
    load_vbd = MOS2type * $limit(V(b, d_int), DEVlimitNewSet, MOS2type * load_vbd, limited);
    load_vbd = load_vbs-load_vds;
    vgd = load_vgs-load_vds;
    vgb = load_vgs-load_vbs;
    if (load_vbs<=-3*load_vt) begin
      MOS2gbs = lc_gmin/$mfactor;
      MOS2cbs = lc_gmin/$mfactor*load_vbs-SourceSatCur;
    end else begin
      evbs = exp(((709.0<load_vbs/load_vt) ? (709.0) : (load_vbs/load_vt)));
      MOS2gbs = SourceSatCur*evbs/load_vt+lc_gmin/$mfactor;
      MOS2cbs = SourceSatCur*(evbs-1)+lc_gmin/$mfactor*load_vbs;
    end
    if (load_vbd<=-3*load_vt) begin
      MOS2gbd = lc_gmin/$mfactor;
      MOS2cbd = lc_gmin/$mfactor*load_vbd-DrainSatCur;
    end else begin
      evbd = exp(((709.0<load_vbd/load_vt) ? (709.0) : (load_vbd/load_vt)));
      MOS2gbd = DrainSatCur*evbd/load_vt+lc_gmin/$mfactor;
      MOS2cbd = DrainSatCur*(evbd-1)+lc_gmin/$mfactor*load_vbd;
    end
    if (load_vds>=0) begin
      MOS2mode = 1;
    end else begin
      MOS2mode = -1;
    end
    sphi = 0.0;
    sphi3 = 0.0;
    argd = 0.0;
    args = 0.0;
    argxs = 0.0;
    argxd = 0.0;
    xn = 0.0;
    argg = 0.0;
    dodvds = 0.0;
    dxndvd = 0.0;
    dxndvb = 0.0;
    xvalid = 0.0;
    bsarg = 0.0;
    dbsrdb = 0.0;
    bodys = 0.0;
    gdbdvs = 0.0;
    xlamda = lambda;
    lvbs = ((MOS2mode==1) ? (load_vbs) : (load_vbd));
    lvds = MOS2mode*load_vds;
    lvgs = ((MOS2mode==1) ? (load_vgs) : (vgd));
    phiMinVbs = MOS2tPhi-lvbs;
    if (lvbs<=0.0) begin
      sarg1 = sqrt(phiMinVbs);
      dsrgdb = -0.5/sarg1;
      d2sdb2 = 0.5*dsrgdb/phiMinVbs;
    end else begin
      sphi = sqrt(MOS2tPhi);
      sphi3 = MOS2tPhi*sphi;
      sarg1 = sphi/(1.0+0.5*lvbs/MOS2tPhi);
      tmp = sarg1/sphi3;
      dsrgdb = -0.5*sarg1*tmp;
      d2sdb2 = -dsrgdb*tmp;
    end
    if (lvbs-lvds<=0) begin
      barg = sqrt(phiMinVbs+lvds);
      dbrgdb = -0.5/barg;
      d2bdb2 = 0.5*dbrgdb/(phiMinVbs+lvds);
    end else begin
      sphi = sqrt(MOS2tPhi);
      sphi3 = MOS2tPhi*sphi;
      barg = sphi/(1.0+0.5*(lvbs-lvds)/MOS2tPhi);
      tmp = barg/sphi3;
      dbrgdb = -0.5*barg*tmp;
      d2bdb2 = -dbrgdb*tmp;
    end
    factor = 0.125*delta*2.0*3.1415926535897932384626433832795*(11.7*8.854214871e-12)/unscaledOxideCap*EffectiveLength;
    eta = 1.0+factor;
    vbin = MOS2tVbi*MOS2type+factor*phiMinVbs;
    if (MOS2gamma>0.0||MOS2substrateDoping>0.0) begin
      xwd = MOS2xd*barg;
      xws = MOS2xd*sarg1;
      argss = 0.0;
      argsd = 0.0;
      dbargs = 0.0;
      dbargd = 0.0;
      dgdvds = 0.0;
      dgddb2 = 0.0;
      if (xj>0) begin
        tmp = 2.0/xj;
        argxs = 1.0+xws*tmp;
        argxd = 1.0+xwd*tmp;
        args = sqrt(argxs);
        argd = sqrt(argxd);
        tmp = 0.5*xj/EffectiveLength;
        argss = tmp*(args-1.0);
        argsd = tmp*(argd-1.0);
      end
      gamasd = MOS2gamma*(1.0-argss-argsd);
      dbxwd = MOS2xd*dbrgdb;
      dbxws = MOS2xd*dsrgdb;
      if (xj>0) begin
        tmp = 0.5/EffectiveLength;
        dbargs = tmp*dbxws/args;
        dbargd = tmp*dbxwd/argd;
        dasdb2 = -MOS2xd*(d2sdb2+dsrgdb*dsrgdb*MOS2xd/(xj*argxs))/(EffectiveLength*args);
        daddb2 = -MOS2xd*(d2bdb2+dbrgdb*dbrgdb*MOS2xd/(xj*argxd))/(EffectiveLength*argd);
        dgddb2 = -0.5*MOS2gamma*(dasdb2+daddb2);
      end
      dgddvb = -MOS2gamma*(dbargs+dbargd);
      if (xj>0) begin
        ddxwd = -dbxwd;
        dgdvds = -MOS2gamma*0.5*ddxwd/(EffectiveLength*argd);
      end
    end else begin
      gamasd = MOS2gamma;
      gammad = MOS2gamma;
      dgddvb = 0.0;
      dgdvds = 0.0;
      dgddb2 = 0.0;
    end
    load_von = vbin+gamasd*sarg1;
    vth = load_von;
    load_vdsat = 0.0;
    if (nfs!=0.0&&OxideCap!=0.0) begin
      cfs = 1.6021766208e-19*nfs*1e4;
      cdonco = -(gamasd*dsrgdb+dgddvb*sarg1)+factor;
      xn = 1.0+cfs/unscaledOxideCap*MOS2w*EffectiveLength+cdonco;
      tmp = load_vt*xn;
      load_von = load_von+tmp;
      argg = 1.0/tmp;
      vgst = lvgs-load_von;
    end else begin
      vgst = lvgs-load_von;
      if (lvgs<=vbin) begin
        MOS2gds = 0.0;
        goto_destination = line1050;
      end
    end
    if (goto_destination<0) begin
      sarg3 = sarg1*sarg1*sarg1;
      sbiarg = sqrt(MOS2tBulkPot);
      gammad = gamasd;
      dgdvbs = dgddvb;
      body = barg*barg*barg-sarg3;
      gdbdv = 2.0*gammad*(barg*barg*dbrgdb-sarg1*sarg1*dsrgdb);
      dodvbs = -factor+dgdvbs*sarg1+gammad*dsrgdb;
      if (nfs==0.0) begin
        goto_destination = line400;
      end
      if (goto_destination<0) begin
        if (OxideCap==0.0) begin
          goto_destination = line410;
        end
        if (goto_destination<0) begin
          dxndvb = 2.0*dgdvbs*dsrgdb+gammad*d2sdb2+dgddb2*sarg1;
          dodvbs = dodvbs+load_vt*dxndvb;
          dxndvd = dgdvds*dsrgdb;
          dodvds = dgdvds*sarg1+load_vt*dxndvd;
        end
      end
    end
    if (goto_destination<0||goto_destination==line400) begin
      goto_destination = -1;
      if (OxideCap<=0.0) begin
        goto_destination = line410;
      end
      if (goto_destination<0) begin
        udenom = vgst;
        tmp = ucrit*100*(11.7*8.854214871e-12)/MOS2oxideCapFactor;
        if (udenom<=tmp) begin
          goto_destination = line410;
        end
        if (goto_destination<0) begin
          ufact = exp(uexp*ln(tmp/udenom));
          ueff = u0*1e-4*ufact;
          dudvgs = -ufact*uexp/udenom;
          dudvds = 0.0;
          dudvbs = uexp*ufact*dodvbs/vgst;
          goto_destination = line500;
        end
      end
    end
    if (goto_destination<0||goto_destination==line410) begin
      goto_destination = -1;
      ufact = 1.0;
      ueff = u0*1e-4;
      dudvgs = 0.0;
      dudvds = 0.0;
      dudvbs = 0.0;
    end
    if (goto_destination<0||goto_destination==line500) begin
      goto_destination = -1;
      vgsx = lvgs;
      gammad = gamasd/eta;
      dgdvbs = dgddvb;
      if (nfs!=0&&OxideCap!=0) begin
        vgsx = ((lvgs>load_von) ? (lvgs) : (load_von));
      end
      if (gammad>0) begin
        gammd2 = gammad*gammad;
        argv = (vgsx-vbin)/eta+phiMinVbs;
        if (argv<=0.0) begin
          load_vdsat = 0.0;
          dsdvgs = 0.0;
          dsdvbs = 0.0;
        end else begin
          load_arg1 = sqrt(1.0+4.0*argv/gammd2);
          load_vdsat = (vgsx-vbin)/eta+gammd2*(1.0-load_arg1)/2.0;
          load_vdsat = ((load_vdsat>0.0) ? (load_vdsat) : (0.0));
          dsdvgs = (1.0-1.0/load_arg1)/eta;
          dsdvbs = (gammad*(1.0-load_arg1)+2.0*argv/(gammad*load_arg1))/eta*dgdvbs+1.0/load_arg1+factor*dsdvgs;
        end
      end else begin
        load_vdsat = (vgsx-vbin)/eta;
        load_vdsat = ((load_vdsat>0.0) ? (load_vdsat) : (0.0));
        dsdvgs = 1.0;
        dsdvbs = 0.0;
      end
      if (vmax>0) begin
        gammd2 = gammad*gammad;
        v1 = (vgsx-vbin)/eta+phiMinVbs;
        v2 = phiMinVbs;
        xv = vmax*EffectiveLength/ueff;
        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;
        load_b = a1*c1-4.0*d1;
        c = -d1*(a1*a1-4.0*b1)-c1*c1;
        r = -a*a/3.0+load_b;
        load_s = 2.0*a*a*a/27.0-a*load_b/3.0+c;
        r3 = r*r*r;
        s2 = load_s*load_s;
        p = s2/4.0+r3/27.0;
        p0 = abs(p);
        p2 = sqrt(p0);
        if (p<0) begin
          ro = sqrt(s2/4.0+p0);
          ro = ln(ro)/3.0;
          ro = exp(ro);
          fi = atan(-2.0*p2/load_s);
          y3 = 2.0*ro*cos(fi/3.0)-a/3.0;
        end else begin
          p3 = -load_s/2.0+p2;
          p3 = exp(ln(abs(p3))/3.0);
          p4 = -load_s/2.0-p2;
          p4 = exp(ln(abs(p4))/3.0);
          y3 = p3+p4-a/3.0;
        end
        iknt = 0;
        a3 = sqrt(a1*a1/4.0-b1+y3);
        b3 = sqrt(y3*y3/4.0-d1);
        i = 1;
        if (i<=4) begin
          a4_0 = a1/2.0+sig1_0*a3;
          b4_0 = y3/2.0+sig2_0*b3;
          delta4 = a4_0*a4_0/4.0-b4_0;
          if (!(delta4<0)) begin
            iknt = iknt+1;
            tmp = sqrt(delta4);
            if (iknt-1==0) begin
              x4_0 = -a4_0/2.0+tmp;
            end else if (iknt-1==1) begin
              x4_1 = -a4_0/2.0+tmp;
            end else if (iknt-1==2) begin
              x4_2 = -a4_0/2.0+tmp;
            end else if (iknt-1==3) begin
              x4_3 = -a4_0/2.0+tmp;
            end else if (iknt-1==4) begin
              x4_4 = -a4_0/2.0+tmp;
            end else if (iknt-1==5) begin
              x4_5 = -a4_0/2.0+tmp;
            end else if (iknt-1==6) begin
              x4_6 = -a4_0/2.0+tmp;
            end else if (iknt-1==7) begin
              x4_7 = -a4_0/2.0+tmp;
            end
            iknt = iknt+1;
            if (iknt-1==0) begin
              x4_0 = -a4_0/2.0-tmp;
            end else if (iknt-1==1) begin
              x4_1 = -a4_0/2.0-tmp;
            end else if (iknt-1==2) begin
              x4_2 = -a4_0/2.0-tmp;
            end else if (iknt-1==3) begin
              x4_3 = -a4_0/2.0-tmp;
            end else if (iknt-1==4) begin
              x4_4 = -a4_0/2.0-tmp;
            end else if (iknt-1==5) begin
              x4_5 = -a4_0/2.0-tmp;
            end else if (iknt-1==6) begin
              x4_6 = -a4_0/2.0-tmp;
            end else if (iknt-1==7) begin
              x4_7 = -a4_0/2.0-tmp;
            end
          end
          i = i+1;
        end
        if (i<=4) begin
          a4_1 = a1/2.0+sig1_1*a3;
          b4_1 = y3/2.0+sig2_1*b3;
          delta4 = a4_1*a4_1/4.0-b4_1;
          if (!(delta4<0)) begin
            iknt = iknt+1;
            tmp = sqrt(delta4);
            if (iknt-1==0) begin
              x4_0 = -a4_1/2.0+tmp;
            end else if (iknt-1==1) begin
              x4_1 = -a4_1/2.0+tmp;
            end else if (iknt-1==2) begin
              x4_2 = -a4_1/2.0+tmp;
            end else if (iknt-1==3) begin
              x4_3 = -a4_1/2.0+tmp;
            end else if (iknt-1==4) begin
              x4_4 = -a4_1/2.0+tmp;
            end else if (iknt-1==5) begin
              x4_5 = -a4_1/2.0+tmp;
            end else if (iknt-1==6) begin
              x4_6 = -a4_1/2.0+tmp;
            end else if (iknt-1==7) begin
              x4_7 = -a4_1/2.0+tmp;
            end
            iknt = iknt+1;
            if (iknt-1==0) begin
              x4_0 = -a4_1/2.0-tmp;
            end else if (iknt-1==1) begin
              x4_1 = -a4_1/2.0-tmp;
            end else if (iknt-1==2) begin
              x4_2 = -a4_1/2.0-tmp;
            end else if (iknt-1==3) begin
              x4_3 = -a4_1/2.0-tmp;
            end else if (iknt-1==4) begin
              x4_4 = -a4_1/2.0-tmp;
            end else if (iknt-1==5) begin
              x4_5 = -a4_1/2.0-tmp;
            end else if (iknt-1==6) begin
              x4_6 = -a4_1/2.0-tmp;
            end else if (iknt-1==7) begin
              x4_7 = -a4_1/2.0-tmp;
            end
          end
          i = i+1;
        end
        if (i<=4) begin
          a4_2 = a1/2.0+sig1_2*a3;
          b4_2 = y3/2.0+sig2_2*b3;
          delta4 = a4_2*a4_2/4.0-b4_2;
          if (!(delta4<0)) begin
            iknt = iknt+1;
            tmp = sqrt(delta4);
            if (iknt-1==0) begin
              x4_0 = -a4_2/2.0+tmp;
            end else if (iknt-1==1) begin
              x4_1 = -a4_2/2.0+tmp;
            end else if (iknt-1==2) begin
              x4_2 = -a4_2/2.0+tmp;
            end else if (iknt-1==3) begin
              x4_3 = -a4_2/2.0+tmp;
            end else if (iknt-1==4) begin
              x4_4 = -a4_2/2.0+tmp;
            end else if (iknt-1==5) begin
              x4_5 = -a4_2/2.0+tmp;
            end else if (iknt-1==6) begin
              x4_6 = -a4_2/2.0+tmp;
            end else if (iknt-1==7) begin
              x4_7 = -a4_2/2.0+tmp;
            end
            iknt = iknt+1;
            if (iknt-1==0) begin
              x4_0 = -a4_2/2.0-tmp;
            end else if (iknt-1==1) begin
              x4_1 = -a4_2/2.0-tmp;
            end else if (iknt-1==2) begin
              x4_2 = -a4_2/2.0-tmp;
            end else if (iknt-1==3) begin
              x4_3 = -a4_2/2.0-tmp;
            end else if (iknt-1==4) begin
              x4_4 = -a4_2/2.0-tmp;
            end else if (iknt-1==5) begin
              x4_5 = -a4_2/2.0-tmp;
            end else if (iknt-1==6) begin
              x4_6 = -a4_2/2.0-tmp;
            end else if (iknt-1==7) begin
              x4_7 = -a4_2/2.0-tmp;
            end
          end
          i = i+1;
        end
        if (i<=4) begin
          a4_3 = a1/2.0+sig1_3*a3;
          b4_3 = y3/2.0+sig2_3*b3;
          delta4 = a4_3*a4_3/4.0-b4_3;
          if (!(delta4<0)) begin
            iknt = iknt+1;
            tmp = sqrt(delta4);
            if (iknt-1==0) begin
              x4_0 = -a4_3/2.0+tmp;
            end else if (iknt-1==1) begin
              x4_1 = -a4_3/2.0+tmp;
            end else if (iknt-1==2) begin
              x4_2 = -a4_3/2.0+tmp;
            end else if (iknt-1==3) begin
              x4_3 = -a4_3/2.0+tmp;
            end else if (iknt-1==4) begin
              x4_4 = -a4_3/2.0+tmp;
            end else if (iknt-1==5) begin
              x4_5 = -a4_3/2.0+tmp;
            end else if (iknt-1==6) begin
              x4_6 = -a4_3/2.0+tmp;
            end else if (iknt-1==7) begin
              x4_7 = -a4_3/2.0+tmp;
            end
            iknt = iknt+1;
            if (iknt-1==0) begin
              x4_0 = -a4_3/2.0-tmp;
            end else if (iknt-1==1) begin
              x4_1 = -a4_3/2.0-tmp;
            end else if (iknt-1==2) begin
              x4_2 = -a4_3/2.0-tmp;
            end else if (iknt-1==3) begin
              x4_3 = -a4_3/2.0-tmp;
            end else if (iknt-1==4) begin
              x4_4 = -a4_3/2.0-tmp;
            end else if (iknt-1==5) begin
              x4_5 = -a4_3/2.0-tmp;
            end else if (iknt-1==6) begin
              x4_6 = -a4_3/2.0-tmp;
            end else if (iknt-1==7) begin
              x4_7 = -a4_3/2.0-tmp;
            end
          end
          i = i+1;
        end
        jknt = 0;
        j = 1;
        if (j<=iknt) begin
          if (!(x4_0<=0)) begin
            poly4_0 = x4_0*x4_0*x4_0*x4_0+a1*x4_0*x4_0*x4_0;
            poly4_0 = poly4_0+b1*x4_0*x4_0+c1*x4_0+d1;
            if (!(abs(poly4_0)>1.0e-6)) begin
              jknt = jknt+1;
              if (jknt<=1) begin
                xvalid = x4_0;
              end
              if (!(x4_0>xvalid)) begin
                xvalid = x4_0;
              end
            end
          end
          j = j+1;
        end
        if (j<=iknt) begin
          if (!(x4_1<=0)) begin
            poly4_1 = x4_1*x4_1*x4_1*x4_1+a1*x4_1*x4_1*x4_1;
            poly4_1 = poly4_1+b1*x4_1*x4_1+c1*x4_1+d1;
            if (!(abs(poly4_1)>1.0e-6)) begin
              jknt = jknt+1;
              if (jknt<=1) begin
                xvalid = x4_1;
              end
              if (!(x4_1>xvalid)) begin
                xvalid = x4_1;
              end
            end
          end
          j = j+1;
        end
        if (j<=iknt) begin
          if (!(x4_2<=0)) begin
            poly4_2 = x4_2*x4_2*x4_2*x4_2+a1*x4_2*x4_2*x4_2;
            poly4_2 = poly4_2+b1*x4_2*x4_2+c1*x4_2+d1;
            if (!(abs(poly4_2)>1.0e-6)) begin
              jknt = jknt+1;
              if (jknt<=1) begin
                xvalid = x4_2;
              end
              if (!(x4_2>xvalid)) begin
                xvalid = x4_2;
              end
            end
          end
          j = j+1;
        end
        if (j<=iknt) begin
          if (!(x4_3<=0)) begin
            poly4_3 = x4_3*x4_3*x4_3*x4_3+a1*x4_3*x4_3*x4_3;
            poly4_3 = poly4_3+b1*x4_3*x4_3+c1*x4_3+d1;
            if (!(abs(poly4_3)>1.0e-6)) begin
              jknt = jknt+1;
              if (jknt<=1) begin
                xvalid = x4_3;
              end
              if (!(x4_3>xvalid)) begin
                xvalid = x4_3;
              end
            end
          end
          j = j+1;
        end
        if (j<=iknt) begin
          if (!(x4_4<=0)) begin
            poly4_4 = x4_4*x4_4*x4_4*x4_4+a1*x4_4*x4_4*x4_4;
            poly4_4 = poly4_4+b1*x4_4*x4_4+c1*x4_4+d1;
            if (!(abs(poly4_4)>1.0e-6)) begin
              jknt = jknt+1;
              if (jknt<=1) begin
                xvalid = x4_4;
              end
              if (!(x4_4>xvalid)) begin
                xvalid = x4_4;
              end
            end
          end
          j = j+1;
        end
        if (j<=iknt) begin
          if (!(x4_5<=0)) begin
            poly4_5 = x4_5*x4_5*x4_5*x4_5+a1*x4_5*x4_5*x4_5;
            poly4_5 = poly4_5+b1*x4_5*x4_5+c1*x4_5+d1;
            if (!(abs(poly4_5)>1.0e-6)) begin
              jknt = jknt+1;
              if (jknt<=1) begin
                xvalid = x4_5;
              end
              if (!(x4_5>xvalid)) begin
                xvalid = x4_5;
              end
            end
          end
          j = j+1;
        end
        if (j<=iknt) begin
          if (!(x4_6<=0)) begin
            poly4_6 = x4_6*x4_6*x4_6*x4_6+a1*x4_6*x4_6*x4_6;
            poly4_6 = poly4_6+b1*x4_6*x4_6+c1*x4_6+d1;
            if (!(abs(poly4_6)>1.0e-6)) begin
              jknt = jknt+1;
              if (jknt<=1) begin
                xvalid = x4_6;
              end
              if (!(x4_6>xvalid)) begin
                xvalid = x4_6;
              end
            end
          end
          j = j+1;
        end
        if (j<=iknt) begin
          if (!(x4_7<=0)) begin
            poly4_7 = x4_7*x4_7*x4_7*x4_7+a1*x4_7*x4_7*x4_7;
            poly4_7 = poly4_7+b1*x4_7*x4_7+c1*x4_7+d1;
            if (!(abs(poly4_7)>1.0e-6)) begin
              jknt = jknt+1;
              if (jknt<=1) begin
                xvalid = x4_7;
              end
              if (!(x4_7>xvalid)) begin
                xvalid = x4_7;
              end
            end
          end
          j = j+1;
        end
        if (jknt>0) begin
          load_vdsat = xvalid*xvalid-phiMinVbs;
        end
      end
      if (lvds!=0.0) begin
        gammad = gamasd;
        if (lvbs-load_vdsat<=0) begin
          bsarg = sqrt(load_vdsat+phiMinVbs);
          dbsrdb = -0.5/bsarg;
        end else begin
          sphi = sqrt(MOS2tPhi);
          sphi3 = MOS2tPhi*sphi;
          bsarg = sphi/(1.0+0.5*(lvbs-load_vdsat)/MOS2tPhi);
          dbsrdb = -0.5*bsarg*bsarg/sphi3;
        end
        bodys = bsarg*bsarg*bsarg-sarg3;
        gdbdvs = 2.0*gammad*(bsarg*bsarg*dbsrdb-sarg1*sarg1*dsrgdb);
        if (vmax<=0) begin
          if (!(MOS2substrateDoping==0.0||xlamda>0.0)) begin
            argv = (lvds-load_vdsat)/4.0;
            sargv = sqrt(1.0+argv*argv);
            load_arg1 = sqrt(argv+sargv);
            xlfact = MOS2xd/(EffectiveLength*lvds);
            xlamda = xlfact*load_arg1;
            dldsat = lvds*xlamda/(8.0*sargv);
            dldvgs = dldsat*dsdvgs;
            dldvds = -xlamda+dldsat;
            dldvbs = dldsat*dsdvbs;
          end else begin
            dldvgs = 0.0;
            dldvds = 0.0;
            dldvbs = 0.0;
          end
        end else begin
          argv = (vgsx-vbin)/eta-load_vdsat;
          xdv = MOS2xd/sqrt(neff);
          xlv = vmax*xdv/(2.0*ueff);
          vqchan = argv-gammad*bsarg;
          dqdsat = -1.0+gammad*dbsrdb;
          vl = vmax*EffectiveLength;
          dfunds = vl*dqdsat-ueff*vqchan;
          dfundg = (vl-ueff*load_vdsat)/eta;
          dfundb = -vl*(1.0+dqdsat-factor/eta)+ueff*(gdbdvs-dgdvbs*bodys/1.5)/eta;
          dsdvgs = -dfundg/dfunds;
          dsdvbs = -dfundb/dfunds;
          if (!(MOS2substrateDoping==0.0||xlamda>0.0)) begin
            argv = lvds-load_vdsat;
            argv = ((argv>0.0) ? (argv) : (0.0));
            xls = sqrt(xlv*xlv+argv);
            dldsat = xdv/(2.0*xls);
            xlfact = xdv/(EffectiveLength*lvds);
            xlamda = xlfact*(xls-xlv);
            dldsat = dldsat/EffectiveLength;
            dldvgs = dldsat*dsdvgs;
            dldvds = -xlamda+dldsat;
            dldvbs = dldsat*dsdvbs;
          end else begin
            dldvgs = 0.0;
            dldvds = 0.0;
            dldvbs = 0.0;
          end
        end
      end else begin
        dldvgs = 0.0;
        dldvds = 0.0;
        dldvbs = 0.0;
      end
      xwb = MOS2xd*sbiarg;
      xld = EffectiveLength-xwb;
      clfact = 1.0-xlamda*lvds;
      dldvds = -xlamda-dldvds;
      xleff = EffectiveLength*clfact;
      deltal = xlamda*lvds*EffectiveLength;
      if (MOS2substrateDoping==0.0) begin
        xwb = 0.25e-6;
      end
      if (xleff<xwb) begin
        xleff = xwb/(1.0+(deltal-xld)/xwb);
        clfact = xleff/EffectiveLength;
        dfact = xleff*xleff/(xwb*xwb);
        dldvgs = dfact*dldvgs;
        dldvds = dfact*dldvds;
        dldvbs = dfact*dldvbs;
      end
      beta1 = Beta*ufact/clfact;
      gammad = gamasd;
      dgdvbs = dgddvb;
      if (lvds<=1.0e-10) begin
        if (lvgs<=load_von) begin
          if (nfs==0.0||OxideCap==0.0) begin
            MOS2gds = 0.0;
            goto_destination = line1050;
          end
          if (goto_destination<0) begin
            MOS2gds = beta1*(load_von-vbin-gammad*sarg1)*exp(argg*(lvgs-load_von));
            goto_destination = line1050;
          end
        end
        if (goto_destination<0) begin
          MOS2gds = beta1*(lvgs-vbin-gammad*sarg1);
          goto_destination = line1050;
        end
      end
      if (goto_destination<0) begin
        if (nfs!=0&&OxideCap!=0) begin
          if (lvgs>load_von) begin
            goto_destination = line900;
          end
        end else begin
          if (lvgs>vbin) begin
            goto_destination = line900;
          end
          if (goto_destination<0) begin
            goto_destination = doneval;
          end
        end
        if (goto_destination<0) begin
          if (lvgs>load_von) begin
            goto_destination = line900;
          end
          if (goto_destination<0) begin
            if (load_vdsat<=0) begin
              MOS2gds = 0.0;
              if (lvgs>vth) begin
                goto_destination = doneval;
              end
              if (goto_destination<0) begin
                goto_destination = line1050;
              end
            end
            if (goto_destination<0) begin
              vdson = ((load_vdsat<lvds) ? (load_vdsat) : (lvds));
              if (lvds>load_vdsat) begin
                barg = bsarg;
                dbrgdb = dbsrdb;
                body = bodys;
                gdbdv = gdbdvs;
              end
              cdson = beta1*((load_von-vbin-eta*vdson*0.5)*vdson-gammad*body/1.5);
              didvds = beta1*(load_von-vbin-eta*vdson-gammad*barg);
              gdson = -cdson*dldvds/clfact-beta1*dgdvds*body/1.5;
              if (lvds<load_vdsat) begin
                gdson = gdson+didvds;
              end
              gbson = -cdson*dldvbs/clfact+beta1*(dodvbs*vdson+factor*vdson-dgdvbs*body/1.5-gdbdv);
              if (lvds>load_vdsat) begin
                gbson = gbson+didvds*dsdvbs;
              end
              expg = exp(argg*(lvgs-load_von));
              cdrain = cdson*expg;
              gmw = cdrain*argg;
              MOS2gm = gmw;
              if (lvds>load_vdsat) begin
                MOS2gm = gmw+didvds*dsdvgs*expg;
              end
              tmp = gmw*(lvgs-load_von)/xn;
              MOS2gds = gdson*expg-MOS2gm*dodvds-tmp*dxndvd;
              MOS2gmbs = gbson*expg-MOS2gm*dodvbs-tmp*dxndvb;
              goto_destination = doneval;
            end
          end
        end
      end
    end
    if (goto_destination<0||goto_destination==line900) begin
      goto_destination = -1;
      if (lvds<=load_vdsat) begin
        cdrain = beta1*((lvgs-vbin-eta*lvds/2.0)*lvds-gammad*body/1.5);
        load_arg1 = cdrain*(dudvgs/ufact-dldvgs/clfact);
        MOS2gm = load_arg1+beta1*lvds;
        load_arg1 = cdrain*(dudvds/ufact-dldvds/clfact);
        MOS2gds = load_arg1+beta1*(lvgs-vbin-eta*lvds-gammad*barg-dgdvds*body/1.5);
        load_arg1 = cdrain*(dudvbs/ufact-dldvbs/clfact);
        MOS2gmbs = load_arg1-beta1*(gdbdv+dgdvbs*body/1.5-factor*lvds);
      end else begin
        cdrain = beta1*((lvgs-vbin-eta*load_vdsat/2.0)*load_vdsat-gammad*bodys/1.5);
        load_arg1 = cdrain*(dudvgs/ufact-dldvgs/clfact);
        MOS2gm = load_arg1+beta1*load_vdsat+beta1*(lvgs-vbin-eta*load_vdsat-gammad*bsarg)*dsdvgs;
        MOS2gds = -cdrain*dldvds/clfact-beta1*dgdvds*bodys/1.5;
        load_arg1 = cdrain*(dudvbs/ufact-dldvbs/clfact);
        MOS2gmbs = load_arg1-beta1*(gdbdvs+dgdvbs*bodys/1.5-factor*load_vdsat)+beta1*(lvgs-vbin-eta*load_vdsat-gammad*bsarg)*dsdvbs;
      end
      goto_destination = doneval;
    end
    if (goto_destination<0||goto_destination==line1050) begin
      goto_destination = -1;
      cdrain = 0.0;
      MOS2gm = 0.0;
      MOS2gmbs = 0.0;
    end
    goto_destination = -1;
    MOS2von = MOS2type*load_von;
    MOS2vdsat = MOS2type*load_vdsat;
    MOS2dcCurrent = MOS2mode*cdrain-MOS2cbd;
    MOS2cd = MOS2dcCurrent;
    if (OxideCap==0) begin
        meyer_scale = 0;
        meyer_unscale = 1;
    end else begin
        meyer_scale = OxideCap;
        meyer_unscale = OxideCap;
    end
    if (MOS2Cbs!=0||MOS2Cbssw!=0) begin
      if (load_vbs<MOS2tDepCap) begin
        load_arg = 1-load_vbs/MOS2tBulkPot;
        if (mj==mjsw) begin
          if (mj==0.5) begin
            load_sargsw = 1/sqrt(load_arg);
            load_sarg = load_sargsw;
          end else begin
            load_sargsw = exp(-mj*ln(load_arg));
            load_sarg = load_sargsw;
          end
        end else begin
          if (mj==0.5) begin
            load_sarg = 1/sqrt(load_arg);
          end else begin
            load_sarg = exp(-mj*ln(load_arg));
          end
          if (mjsw==0.5) begin
            load_sargsw = 1/sqrt(load_arg);
          end else begin
            load_sargsw = exp(-mjsw*ln(load_arg));
          end
        end
        MOS2qbs = MOS2tBulkPot*(MOS2Cbs*(1-load_arg*load_sarg)/(1-mj)+MOS2Cbssw*(1-load_arg*load_sargsw)/(1-mjsw));
        MOS2capbs = MOS2Cbs*load_sarg+MOS2Cbssw*load_sargsw;
      end else begin
        MOS2qbs = MOS2f4s+load_vbs*(MOS2f2s+load_vbs*(MOS2f3s/2));
        MOS2capbs = MOS2f2s+MOS2f3s*load_vbs;
      end
    end else begin
      MOS2qbs = 0;
      MOS2capbs = 0;
    end
    if (MOS2Cbd!=0||MOS2Cbdsw!=0) begin
      if (load_vbd<MOS2tDepCap) begin
        load_arg = 1-load_vbd/MOS2tBulkPot;
        if (mj==0.5&&mjsw==0.5) begin
          load_sargsw = 1/sqrt(load_arg);
          load_sarg = load_sargsw;
        end else begin
          if (mj==0.5) begin
            load_sarg = 1/sqrt(load_arg);
          end else begin
            load_sarg = exp(-mj*ln(load_arg));
          end
          if (mjsw==0.5) begin
            load_sargsw = 1/sqrt(load_arg);
          end else begin
            load_sargsw = exp(-mjsw*ln(load_arg));
          end
        end
        MOS2qbd = MOS2tBulkPot*(MOS2Cbd*(1-load_arg*load_sarg)/(1-mj)+MOS2Cbdsw*(1-load_arg*load_sargsw)/(1-mjsw));
        MOS2capbd = MOS2Cbd*load_sarg+MOS2Cbdsw*load_sargsw;
      end else begin
        MOS2qbd = MOS2f4d+load_vbd*(MOS2f2d+load_vbd*MOS2f3d/2);
        MOS2capbd = MOS2f2d+load_vbd*MOS2f3d;
      end
    end else begin
      MOS2qbd = 0;
      MOS2capbd = 0;
    end
    geq = 0;
    MOS2cqbd = ddt(MOS2qbd);
    ceq = MOS2cqbd;
    MOS2gbd = MOS2gbd + (geq);
    MOS2cbd = MOS2cbd + (MOS2cqbd);
    MOS2cd = MOS2cd - (MOS2cqbd);
    geq = 0;
    MOS2cqbs = ddt(MOS2qbs);
    ceq = MOS2cqbs;
    MOS2gbs = MOS2gbs + (geq);
    MOS2cbs = MOS2cbs + (MOS2cqbs);
    MOS2vbs = load_vbs;
    MOS2vbd = load_vbd;
    MOS2vgs = load_vgs;
    MOS2vds = load_vds;
    if (MOS2mode>0) begin
      DEVqmeyer(load_vgs, vgd, vgb, load_von, load_vdsat, MOS2capgs, MOS2capgd, MOS2capgb, MOS2tPhi, OxideCap);
    end else begin
      DEVqmeyer(vgd, load_vgs, vgb, load_von, load_vdsat, MOS2capgd, MOS2capgs, MOS2capgb, MOS2tPhi, OxideCap);
    end
    capgs = MOS2capgs+MOS2capgs+GateSourceOverlapCap;
    capgd = MOS2capgd+MOS2capgd+GateDrainOverlapCap;
    capgb = MOS2capgb+MOS2capgb+GateBulkOverlapCap;
    if (capgs==0) begin
      MOS2cqgs = 0;
    end
    if (capgd==0) begin
      MOS2cqgd = 0;
    end
    if (capgb==0) begin
      MOS2cqgb = 0;
    end
    gcgs = 0;
    ceqgs = capgs*(ddt(meyer_scale*load_vgs)/meyer_unscale);
    MOS2cqgs = ceqgs;
    gcgd = 0;
    ceqgd = capgd*(ddt(meyer_scale*vgd)/meyer_unscale);
    MOS2cqgd = ceqgd;
    gcgb = 0;
    ceqgb = capgb*(ddt(meyer_scale*vgb)/meyer_unscale);
    MOS2cqgb = ceqgb;
    ceqbs = MOS2type*MOS2cbs;
    ceqbd = MOS2type*MOS2cbd;
    if (MOS2mode>=0) begin
      xnrm = 1;
      xrev = 0;
      cdreq = MOS2type*cdrain;
    end else begin
      xnrm = 0;
      xrev = 1;
      cdreq = -MOS2type*cdrain;
    end
    
    // Residuals
    I(d) <+ MOS2drainConductance*V(d, d_int);
    I(g) <+ gcgb*V(g, b)+gcgd*V(g, d_int)+gcgs*V(g, s_int)+MOS2type*(ceqgs+ceqgb+ceqgd);
    I(s) <+ MOS2sourceConductance*V(s, s_int);
    I(b) <+ gcgb*V(b, g)+(ceqbs+ceqbd-MOS2type*ceqgb);
    I(d_int) <+ MOS2drainConductance*V(d_int, d)+gcgd*V(d_int, g)+-(ceqbd-cdreq+MOS2type*ceqgd);
    I(s_int) <+ MOS2sourceConductance*V(s_int, s)+gcgs*V(s_int, g)+-(cdreq+ceqbs+MOS2type*ceqgs);

    // Tag: noise, function: MOS2noise
    if (MOS2tempGiven) begin
      noise_dtemp = MOS2temp-$temperature+(MOS2tnom-273.15);
    end else begin
      noise_dtemp = dtemp;
    end
    I(d_int, d) <+ white_noise(4*(VACONSTboltz*(($temperature+noise_dtemp)*MOS2drainConductance)), "rd");
    I(s_int, s) <+ white_noise(4*(VACONSTboltz*(($temperature+noise_dtemp)*MOS2sourceConductance)), "rs");
    if (nlev<3) begin
      Sid = 2.0/3.0*abs(MOS2gm);
    end else begin
      noise_vds = MOS2vds;
      noise_vgs = MOS2vgs;
      noise_vgd = noise_vgs-noise_vds;
      noise_vgst = ((MOS2mode==1) ? (noise_vgs) : (noise_vgd))-MOS2type*MOS2von;
      if (noise_vgst>0) begin
        if (noise_vgst<=noise_vds*MOS2mode) begin
          alpha = 0.0;
        end else begin
          alpha = 1.0-noise_vds*MOS2mode/(MOS2type*MOS2vdsat);
        end
        beta = MOS2tTransconductance*MOS2w/(MOS2l-2*ld);
        Sid = 2.0/3.0*beta*noise_vgst*(1.0+alpha+alpha*alpha)/(1.0+alpha)*gdsnoi;
      end else begin
        Sid = 0.0;
      end
    end
    I(d_int, s_int) <+ white_noise(4*(VACONSTboltz*(($temperature+noise_dtemp)*Sid)), "id");
    case (nlev)
      0: begin
        flicker_psd = ((MOS2dcCurrent>=0) ? (1) : (-1))*(kf*exp(af*ln(((abs(MOS2dcCurrent)>1E-38) ? (abs(MOS2dcCurrent)) : (1E-38))))/((MOS2l-2*ld)*(MOS2l-2*ld)*MOS2oxideCapFactor));
        flicker_exp = 1;
      end
      1: begin
        flicker_psd = ((MOS2dcCurrent>=0) ? (1) : (-1))*(kf*exp(af*ln(((abs(MOS2dcCurrent)>1E-38) ? (abs(MOS2dcCurrent)) : (1E-38))))/(MOS2w*(MOS2l-2*ld)*MOS2oxideCapFactor));
        flicker_exp = 1;
      end
      2, 
      3: begin
        flicker_psd = ((MOS2dcCurrent>=0) ? (1) : (-1))*(kf*MOS2gm*MOS2gm/(MOS2w*(MOS2l-2*ld)*MOS2oxideCapFactor));
        flicker_exp = af;
      end
    endcase

    I(d_int, s_int) <+ flicker_noise(flicker_psd, flicker_exp, "flicker");

    // Instance output variable computation
    i_d = MOS2cd;
    ibd = MOS2cbd;
    ibs = MOS2cbs;
    i_s = -MOS2cd;
    i_s = i_s - (MOS2cbd+MOS2cbs-MOS2cqgb);
    i_s = i_s - (MOS2cqgb+MOS2cqgd+MOS2cqgs);
    i_g = MOS2cqgb+MOS2cqgd+MOS2cqgs;
    i_b = MOS2cbd+MOS2cbs-MOS2cqgb;
    vgs = MOS2vgs;
    vds = MOS2vds;
    vbs = MOS2vbs;
    vbd = MOS2vbd;
    von = MOS2von;
    vdsat = MOS2vdsat;
    sourcevcrit = MOS2sourceVcrit;
    drainvcrit = MOS2drainVcrit;
    sourceconductance = MOS2sourceConductance;
    drainconductance = MOS2drainConductance;
    gm = MOS2gm;
    gds = MOS2gds;
    gmb = MOS2gmbs;
    gbd = MOS2gbd;
    gbs = MOS2gbs;
    c_bd = MOS2capbd;
    c_bs = MOS2capbs;
    cgs = 2*MOS2capgs;
    cgs = cgs + (cgso*MOS2w);
    cgd = 2*MOS2capgd;
    cgd = cgd + (cgdo*MOS2w);
    cgb = 2*MOS2capgb;
    cgb = cgb + (cgbo*(MOS2l-2*ld));
    cbd0 = MOS2Cbd;
    cbdsw0 = MOS2Cbdsw;
    cbs0 = MOS2Cbs;
    cbssw0 = MOS2Cbssw;
    cqgs = MOS2cqgs;
    cqgd = MOS2cqgd;
    cqgb = MOS2cqgb;
    cqbd = MOS2cqbd;
    cqbs = MOS2cqbs;
    qbd = MOS2qbd;
    qbs = MOS2qbs;

    // Model output variable computation

  end
endmodule
