// 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:37:05 UTC
// Module name  : sp_vdmos
// Model variant: default
// 
// 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 cgd, cgs, ig, id, and is were renamed to c_gd, 
//   c_gs, i_g, i_d, and i_s due to conflicts with model parameters. 
//   Parameter thermal renamed to thermal_enabled due to conflict with 
//   a Verilog-A discipline. 
//   
//   Meyer charge is not computed and therefore not available. 
//   Output variables qgs and qgd were removed. 
//   
//   Parameters ic, icvds, and icvgs 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 
//   D-S 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.
//
//   Full noise model used (appropriate for all noise analyses).
//
// SPICE variable mapping
//   scale -> $simparam("scale", 1)
//
// CKTcircuit member mapping
//     CKTepsmin -> ('epsmin', '$simparam("epsmin", 1e-28)')
//   CKTfixLimit -> ('oldlimit', '$simparam("oldlimit", 0)', 'const')
//       CKTgmin -> ('gmin', '$simparam("gmin", 1e-12)')
//    CKTnomTemp -> ('tnom', '($simparam("tnom", 27)+273.15)')
//     CKTreltol -> ('reltol', '$simparam("reltol", 1e-3)', 'const')
//       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
//   vdmos/vdmosdefs.h
//     Copyright 1990 Regents of the University of California.  All rights reserved.
//     Author: 1985 Thomas L. Quarles
//     Modified: 2000 AlansFixes
//     VDMOS: 2018 Holger Vogt, 2020 Dietmar Warning
//   vdmos/vdmos.c
//     Copyright 1990 Regents of the University of California.  All rights reserved.
//     Author: 1987 Thomas L. Quarles
//     Modified: 2000 AlansFixes
//     VDMOS: 2018 Holger Vogt, 2020 Dietmar Warning
//   vdmos/vdmospar.c
//     Copyright 1990 Regents of the University of California.  All rights reserved.
//     Author: 1985 Thomas L. Quarles
//     Modified: 2000 AlansFixes
//     VDMOS: 2018 Holger Vogt, 2020 Dietmar Warning
//   vdmos/vdmosmpar.c
//     Copyright 1990 Regents of the University of California.  All rights reserved.
//     Author: 1985 Thomas L. Quarles
//     VDMOS: 2018 Holger Vogt, 2020 Dietmar Warning
//   vdmos/vdmosask.c
//     Copyright 1990 Regents of the University of California.  All rights reserved.
//     Author: 1987 Thomas L. Quarles
//     Modified: 2000 AlansFixes
//   vdmos/vdmosmask.c
//     Copyright 1990 Regents of the University of California.  All rights reserved.
//     Author: 1987 Thomas L. Quarles
//   vdmos/vdmosset.c
//   vdmos/vdmostemp.c
//   vdmos/vdmosload.c
//     Copyright 1990 Regents of the University of California.  All rights reserved.
//     Author: 1985 Thomas L. Quarles
//     Modified: 2000 AlansFixes
//     VDMOS: 2018 Holger Vogt, 2020 Dietmar Warning
//   vdmos/vdmosnoi.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_vdmos(d, g, s, t, tc);
  inout d, g, s, t, tc;
  electrical d, g, s, t, tc, d_int, g_int, s_int, bd_int, t_int;
  branch (t_int) tbr;

  // Instance parameters
  (* desc = "Instance temperature", type = "instance" *) parameter real temp = 0;
  (* desc = "Instance temperature difference", type = "instance" *) parameter real dtemp = 0.0;
  (* desc = "Thermal model switch on/off", type = "instance" *) parameter integer thermal_enabled = 0;

  // Model parameters
  (* desc = "NMOS (1) or PMOS (-1)" *) parameter integer type = 1;
  (* desc = "Threshold voltage" *) parameter real vto = 0;
  aliasparam vth0 = vto;
  (* desc = "Transconductance parameter" *) parameter real kp = 0;
  (* desc = "Surface potential" *) parameter real phi = 0.6;
  (* desc = "Channel length modulation" *) parameter real lambda = 0;
  (* desc = "Vgs dependence on mobility" *) parameter real theta = 0;
  (* desc = "Drain ohmic resistance" *) parameter real rd = 0.0;
  (* desc = "Source ohmic resistance" *) parameter real rs = 0.0;
  (* desc = "Gate ohmic resistance" *) parameter real rg = 0.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 = "DMOS transistor" *) parameter real vdmos = 1;
  (* desc = "Quasi saturation resistance fitting parameter" *) parameter real rq = 0.0;
  (* desc = "Quasi saturation voltage fitting parameter" *) parameter real vq = 0.0;
  (* desc = "Conductance multiplier in triode region" *) parameter real mtriode = 1.0;
  (* desc = "Linear Vth0 temperature coefficient" *) parameter real tcvth = 0.0;
  aliasparam vtotc = tcvth;
  (* desc = "Exponent of gain temperature dependency" *) parameter real mu = 0;
  aliasparam bex = mu;
  (* desc = "Drain resistance rd0 temperature exponent" *) parameter real texp0 = 1.5;
  (* desc = "Drain resistance rd1 temperature exponent" *) parameter real texp1 = 0.3;
  (* desc = "Drain resistance linear temperature coefficient" *) parameter real trd1 = 0.0;
  (* desc = "Drain resistance quadratic temperature coefficient" *) parameter real trd2 = 0.0;
  (* desc = "Gate resistance linear temperature coefficient" *) parameter real trg1 = 0.0;
  (* desc = "Gate resistance quadratic temperature coefficient" *) parameter real trg2 = 0.0;
  (* desc = "Source resistance linear temperature coefficient" *) parameter real trs1 = 0.0;
  (* desc = "Source resistance quadratic temperature coefficient" *) parameter real trs2 = 0.0;
  (* desc = "Body resistance linear temperature coefficient" *) parameter real trb1 = 0.0;
  (* desc = "Body resistance quadratic temperature coefficient" *) parameter real trb2 = 0.0;
  (* desc = "Shift of weak inversion plot on the vgs axis" *) parameter real subshift = 0;
  (* desc = "Slope of weak inversion log current versus vgs" *) parameter real ksubthres = 0.1;
  (* desc = "Linear temperature coefficient of ksubthres" *) parameter real tksubthres1 = 0.0;
  (* desc = "Quadratic temperature coefficient of ksubthres" *) parameter real tksubthres2 = 0.0;
  (* desc = "Vds breakdown voltage" *) parameter real bv = 1e+99;
  (* desc = "Current at Vds=bv" *) parameter real ibv = 1e-10;
  (* desc = "Vds breakdown emission coefficient" *) parameter real nbv = 1.0;
  (* desc = "Drain-source shunt resistance" *) parameter real rds = 1000000000000000.0;
  (* desc = "Body diode ohmic resistance" *) parameter real rb = 0.0;
  (* desc = "Body diode emission coefficient" *) parameter real n = 1.0;
  (* desc = "Body diode transit time" *) parameter real tt = 0.0;
  (* desc = "Body diode activation energy for temperature effect on Is" *) parameter real eg = 1.11;
  (* desc = "Body diode saturation current temperature exponent" *) parameter real xti = 3.0;
  (* desc = "Body diode saturation current" *) parameter real is = 1e-14;
  (* desc = "Body diode junction potential" *) parameter real vj = 0.8;
  (* desc = "Zero-bias body diode junction capacitance" *) parameter real cjo = 5e-10;
  (* desc = "Body diode grading coefficient" *) parameter real m = 0.5;
  (* desc = "Body diode coefficient for forward-bias depletion capacitance formula" *) parameter real fc = 0.5;
  (* desc = "Minimum non-linear G-D capacitance" *) parameter real cgdmin = 2e-11;
  (* desc = "Maximum non-linear G-D capacitance" *) parameter real cgdmax = 2e-09;
  (* desc = "Non-linear Cgd capacitance parameter" *) parameter real a = 1.0;
  (* desc = "Gate-source capacitance" *) parameter real cgs = 1.4e-09;
  (* desc = "Self-heating thermal resistance, junction-to-case" *) parameter real rthjc = 1.0;
  (* desc = "Self-heating thermal resistance, case-to-ambient" *) parameter real rthca = 23;
  (* desc = "Self-heating thermal capacitance" *) parameter real cthj = 1e-05;
  (* desc = "maximum voltage G-S branch" *) parameter real vgs_max = 1e+99;
  (* desc = "maximum voltage G-D branch" *) parameter real vgd_max = 1e+99;
  (* desc = "maximum voltage D-S branch" *) parameter real vds_max = 1e+99;
  (* desc = "maximum voltage G-S branch" *) parameter real vgsr_max = 1e+99;
  (* desc = "maximum voltage G-D branch" *) parameter real vgdr_max = 1e+99;
  (* desc = "maximum device power dissipation" *) parameter real pd_max = 1e+99;
  (* desc = "maximum drain/source current" *) parameter real id_max = 0;
  (* desc = "maximum drain/source reverse current" *) parameter real idr_max = 1e+99;
  (* desc = "maximum temperature" *) parameter real te_max = 1e+99;
  (* desc = "thermal resistance case to ambient, incl. heat sink" *) parameter real rth_ext = 0;
  (* desc = "thermal derating for power" *) parameter real derating = 0;

  // Instance output variables
  (* desc = "Gate-Source voltage" *) real vgs;
  (* desc = "Drain-Source voltage" *) real vds;
  (* desc = "Gate-Source capacitance" *) real c_gs;
  (* desc = "Gate-Drain capacitance" *) real c_gd;
  (* desc = "Drain-Source capacitance" *) real cds;
  (* desc = "Device on state voltage" *) real von;
  (* desc = "Conductance of source" *) real sourceconductance;
  (* desc = "Conductance of drain" *) real drainconductance;
  (* desc = "Transconductance" *) real gm;
  (* desc = "Drain-Source conductance" *) real gds;
  (* desc = "Capacitance due to gate-source charge storage" *) real cqgs;
  (* desc = "Capacitance due to gate-drain charge storage" *) real cqgd;

  // Model output variables

  // Instance fields
  real VDMOSsourceConductance = 0;
  real VDMOSdrainConductance = 0;
  real i_VDMOSdrainResistance = 0;
  real i_VDMOSqsResistance = 0;
  real VDMOSgateConductance = 0;
  real VDMOSdsConductance = 0;
  real VDMOStemp = 0;
  integer VDMOSthermal = 0;
  real VDMOStTransconductance = 0;
  real VDMOStPhi = 0;
  real VDMOStVth = 0;
  real VDMOStksubthres = 0;
  real VDMOSvon = 0;
  real VDMOSvdsat = 0;
  real VDMOScd = 0;
  real VDMOSgm = 0;
  real VDMOSgds = 0;
  real VDIOcap = 0;
  real VDIOtSatCur = 0;
  real VDIOtSatCur_dT = 0;
  real VDIOinitCond = 0;
  real VDIOtVcrit = 0;
  real VDIOconductance = 0;
  real VDIOtConductance = 0;
  real VDIOtConductance_dT = 0;
  real VDIOtBrkdwnV = 0;
  real VDIOtJctCap = 0;
  real VDIOtDepCap = 0;
  real VDIOtJctPot = 0;
  real VDIOtGradingCoeff = 0;
  real VDIOtTransitTime = 0;
  real VDIOtF1 = 0;
  real VDIOtF2 = 0;
  real VDIOtF3 = 0;
  real VDMOSTempSH = 0;
  real VDMOSgmT = 0;
  real VDMOSgtempg = 0;
  real VDMOSgtempd = 0;
  real VDMOSgtempT = 0;
  real VDMOScgT = 0;
  real VDMOScdT = 0;
  real VDMOScth = 0;
  integer VDMOSmode = 0;
  real VDMOSdcCurrent = 0;
  real capmodel_scale = 0;
  real capmodel_unscale = 0;
  integer VDMOStempGiven = 0;

  // Model fields
  integer VDMOStype = 0;
  real VDMOStnom = 0;
  real VDMOStransconductance = 0;
  real VDMOSoxideCapFactor = 0;
  real VDMOSvth0 = 0;
  real VDIOdepletionCapCoeff = 0;
  real VDIOjctSatCur = 0;
  real VDIOtranTimeTemp1 = 0;
  real VDIOtranTimeTemp2 = 0;
  real VDIOeg = 0;
  real VDIOgradCoeffTemp1 = 0;
  real VDIOgradCoeffTemp2 = 0;
  real VDMOSmu = 0;
  real VDMOSpd_max = 0;
  real VDMOSrth_ext = 0;
  integer VDMOSqsGiven = 0;
  integer VDMOStnomGiven = 0;
  integer VDMOStypeGiven = 0;

  // States
  real VDMOSvgs = 0;
  real VDMOSvds = 0;
  real VDMOSdelTemp = 0;
  real VDMOScapgs = 0;
  real VDMOSqgs = 0;
  real VDMOScqgs = 0;
  real VDMOScapgd = 0;
  real VDMOSqgd = 0;
  real VDMOScqgd = 0;
  real VDIOvoltage = 0;
  real VDIOcurrent = 0;
  real VDIOconduct = 0;
  real VDIOcapCharge = 0;
  real VDIOcapCurrent = 0;
  real VDMOScapth = 0;
  real VDMOSqth = 0;
  real VDMOScqth = 0;
  real VDIOdIdio_dT = 0;

  // Model parameter access locals

  // Instance parameter access locals

  // Device setup locals

  // Device temperature dependence locals
  real egfet = 0;
  real egfet1 = 0;
  real fact1 = 0;
  real fact2 = 0;
  real kt = 0;
  real kt1 = 0;
  real arg1 = 0;
  real ratio = 0;
  real phio = 0;
  real pbfact1 = 0;
  real pbfact = 0;
  real vt = 0;
  real vtnom = 0;
  real xfc = 0;
  real arg = 0;
  real dt = 0;
  real pbo = 0;
  real gmaold = 0;
  real gmanew = 0;
  real factor = 0;
  real tBreakdownVoltage = 0;
  real vte = 0;
  real cbv = 0;
  real xbv = 0;
  real xcbv = 0;
  real tol = 0;
  real iter = 0;
  real arg1_dT = 0;
  real arg2 = 0;
  real arg2_dT = 0;

  // Device load locals
  integer limited = 0;
  real Beta = 0;
  real load_arg = 0;
  real cdrain = 0;
  real cdreq = 0;
  real ceq = 0;
  real ceqgd = 0;
  real ceqgs = 0;
  real gcgd = 0;
  real gcgs = 0;
  real geq = 0;
  real sarg = 0;
  real load_vds = 0;
  real vdsat = 0;
  real vgd = 0;
  real vgdo = 0;
  real load_vgs = 0;
  real load_von = 0;
  real load_vt = 0;
  real xfact = 0;
  integer xnrm = 0;
  integer xrev = 0;
  real capgs = 0;
  real capgd = 0;
  real capth = 0;
  integer selfheat = 0;
  real rd0T = 0;
  real rd1T = 0;
  real dBeta_dT = 0;
  real dIds_dT = 0;
  real Vrd = 0;
  real dIth_dVrd = 0;
  real dIrd_dT = 0;
  real drd0T_dT = 0;
  real drd1T_dT = 0;
  real dgdrain_dT = 0;
  real rsT = 0;
  real Vrs = 0;
  real dIth_dVrs = 0;
  real dIrs_dT = 0;
  real dgsource_dT = 0;
  real delTemp = 0;
  real delTemp1 = 0;
  real Temp = 0;
  real Vds = 0;
  real Vgs = 0;
  real ceqqth = 0;
  real GmT = 0;
  real gTtg = 0;
  real gTtdp = 0;
  real gTtt = 0;
  real gTtsp = 0;
  real gcTt = 0;
  real load_cgdmin = 0;
  real load_cgdmax = 0;
  real load_a = 0;
  real load_cgs = 0;
  real vd = 0;
  real cd = 0;
  real load_vte = 0;
  real vtebrk = 0;
  real vbrknp = 0;
  real cdb = 0;
  real cdeq = 0;
  real capd = 0;
  real gd = 0;
  real gdb = 0;
  real gbpr = 0;
  real delvd = 0;
  real evrev = 0;
  real Ith = 0;
  real dIth_dT = 0;
  real dIdio_dT = 0;
  real dIth_dVdio = 0;
  real vrb = 0;
  real dIrb_dT = 0;
  real dIth_dVrb = 0;
  real load_egfet = 0;
  real load_egfet1 = 0;
  real load_fact1 = 0;
  real load_fact2 = 0;
  real load_kt = 0;
  real load_kt1 = 0;
  real load_arg1 = 0;
  real load_ratio = 0;
  real load_phio = 0;
  real load_pbfact1 = 0;
  real load_pbfact = 0;
  real load1_vt = 0;
  real load_vtnom = 0;
  real load_xfc = 0;
  real load1_arg = 0;
  real load_dt = 0;
  real load_pbo = 0;
  real load_gmaold = 0;
  real load_gmanew = 0;
  real load_factor = 0;
  real load_tBreakdownVoltage = 0;
  real load1_vte = 0;
  real load_cbv = 0;
  real load_xbv = 0;
  real load_xcbv = 0;
  real load_tol = 0;
  real load_iter = 0;
  real load_arg1_dT = 0;
  real load_arg2 = 0;
  real load_arg2_dT = 0;
  real drsT_dT = 0;
  real vgst = 0;
  real slope = 0;
  real load_lambda = 0;
  real load_theta = 0;
  real shift = 0;
  real mtr = 0;
  real vdss = 0;
  real t0 = 0;
  real t1 = 0;
  real betap = 0;
  real dbetapdvgs = 0;
  real dbetapdvds = 0;
  real t2 = 0;
  real dvgstdvgs = 0;
  real dvgst_dT = 0;
  real dvdsat_dT = 0;
  real dt1_dT = 0;
  real dbetap_dT = 0;
  real vdsn = 0;
  real load_rd = 0;
  real drd_dT = 0;
  real dIrd_dgdrain = 0;
  real dIrs_dgsource = 0;
  real vDevTemp = 0;
  real gthjc = 0;
  real gthca = 0;
  real vdtemp = 0;
  real evd = 0;
  real arg3 = 0;
  real darg3_dT = 0;
  real czero = 0;
  real diffcharge = 0;
  real deplcharge = 0;
  real diffcap = 0;
  real deplcap = 0;
  real czof2 = 0;
  real dIrb_dgbpr = 0;
  real dIth_dIrb = 0;

  // Device noise locals
  real coxSquared = 0;
  real tempRatioSH = 0;

  // 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 epsmin = $simparam("epsmin", 1e-28);
  real oldlimit = $simparam("oldlimit", 0);
  real VACONST_tnom = ($simparam("tnom", 27)+273.15);
  real reltol = $simparam("reltol", 1e-3);
  real cpscale = $simparam("scale", 1);

  // Evaluation constants
  real lc_gmin = 0;

  // Translated function definitions
  analog function real DevCapVDMOS;
    input vgd, cgdmin, cgdmax, a, cgs;
    inout capgs, capgd;
    real vgd;
    real cgdmin;
    real cgdmax;
    real a;
    real cgs;
    real capgs;
    real capgd;
    real s;
    real y;
    begin
      DevCapVDMOS = 0;
      s = (cgdmax-cgdmin)/(1+3.1415926535897932384626433832795/2);
      y = cgdmax-s;
      if (vgd>0) begin
        capgd = 0.5*(s*tanh(a*vgd)+y);
      end else begin
        capgd = 0.5*(s*atan(a*vgd)+y);
      end
      capgs = 0.5*cgs;
    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 function real DEVlimitlog;
    input deltemp, deltemp_old, LIM_TOL;
    inout limiting_applied;
    real deltemp, deltemp_old, LIM_TOL, vlimited;
    integer limiting_applied;
    begin
      /* In Verilog-A one cannot cehck for nan values. 
        The first part of this function was removed. */
  
      /* Logarithmic damping of deltemp beyond LIM_TOL */
      if (deltemp > deltemp_old + LIM_TOL) begin
        vlimited = deltemp_old + LIM_TOL + log((deltemp-deltemp_old)/LIM_TOL)/log(10);
        limiting_applied = 1;
      end else if (deltemp < deltemp_old - LIM_TOL) begin
        vlimited = deltemp_old - LIM_TOL - log((deltemp_old-deltemp)/LIM_TOL)/log(10);
        limiting_applied = 1;
      end else begin
        vlimited = deltemp;
      end
      DEVlimitlog = vlimited;
    end
  endfunction

  // Analog block
  analog begin

    // Instance field assignments
    if ($param_given(temp)) begin
      VDMOStemp = temp+273.15;
      VDMOStempGiven = 1;
    end
    if ($param_given(thermal_enabled)) begin
      VDMOSthermal = thermal_enabled!=0;
    end

    // Model field assignments
    VDMOStype = type;
    VDMOStypeGiven = 1;
    if ($param_given(vto)) begin
      VDMOSvth0 = vto;
    end
    if ($param_given(kp)) begin
      VDMOStransconductance = kp;
    end
    if ($param_given(tnom)) begin
      VDMOStnom = tnom+273.15;
      VDMOStnomGiven = 1;
    end
    if ($param_given(vdmos)) begin
      if (vdmos) begin
        VDMOStype = 1;
        VDMOStypeGiven = 1;
      end
    end
    if ($param_given(mu)) begin
      VDMOSmu = mu;
    end
    if ($param_given(eg)) begin
      VDIOeg = eg;
    end
    if ($param_given(is)) begin
      VDIOjctSatCur = is;
    end
    if ($param_given(fc)) begin
      VDIOdepletionCapCoeff = fc;
    end
    if ($param_given(pd_max)) begin
      VDMOSpd_max = pd_max;
    end
    if ($param_given(rth_ext)) begin
      VDMOSrth_ext = rth_ext;
    end

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

    // Tag: setup, function: VDMOSsetup
    if (!VDMOStypeGiven) begin
      VDMOStype = 1;
    end
    if (!$param_given(kp)) begin
      VDMOStransconductance = 15+5*VDMOStype;
    end
    if (!$param_given(vto)) begin
      VDMOSvth0 = 3*VDMOStype;
    end
    if (!$param_given(is)) begin
      VDIOjctSatCur = 1e-14;
    end
    if (!$param_given(fc)) begin
      VDIOdepletionCapCoeff = 0.5;
    end
    if (!$param_given(eg)) begin
      VDIOeg = 1.11;
    end
    if (!$param_given(mu)) begin
      VDMOSmu = -1.5;
    end
    if (!$param_given(pd_max)) begin
      VDMOSpd_max = 1e99;
    end
    if (!$param_given(id_max)) begin
      VDMOSpd_max = 1e99;
    end
    if (!$param_given(rth_ext)) begin
      VDMOSrth_ext = rthca;
    end
    if ($param_given(rq)&&$param_given(vq)) begin
      VDMOSqsGiven = 1;
    end else begin
      VDMOSqsGiven = 0;
    end
    if (!VDMOStnomGiven) begin
      VDMOStnom = VACONST_tnom;
    end
    if (phi<=0.0) begin
      $fatal(0, "Phi is not positive.");
    end
    VDMOSoxideCapFactor = 3.9*8.854214871e-12/1e-07;
    if (VDIOeg<0.1) begin
      $warning("body diode activation energy too small, limited to 0.1");
      VDIOeg = 0.1;
    end
    if (VDIOdepletionCapCoeff>0.95) begin
      $warning("coefficient Fc too large, limited to 0.95");
      VDIOdepletionCapCoeff = 0.95;
    end
    if (VDIOjctSatCur<epsmin) begin
      VDIOjctSatCur = epsmin;
    end
    if (VDMOStempGiven&&$param_given(dtemp)) begin
      $warning("temp and dtemp instance parameter given - using dtemp");
      VDMOStempGiven = 0;
    end
    VDMOSvdsat = 0;
    VDMOSvon = 0;
    if (rd>0) begin
      VDMOSdrainConductance = 1.0/rd;
    end else begin
      VDMOSdrainConductance = 0.0;
    end
    if (rs>0) begin
      VDMOSsourceConductance = 1.0/rs;
    end else begin
      VDMOSsourceConductance = 0.0;
    end
    if (rg>0) begin
      VDMOSgateConductance = 1.0/rg;
    end else begin
      VDMOSgateConductance = 0.0;
    end
    if ($param_given(rds)) begin
      if (rds>0) begin
        VDMOSdsConductance = 1.0/rds;
      end else begin
        VDMOSdsConductance = 1e-15;
      end
    end else begin
      VDMOSdsConductance = 1e-15;
    end
    if (rb>0) begin
      VDIOconductance = 1.0/rb;
    end else begin
      VDIOconductance = 0.0;
    end
    if (!(rd>0)) begin
      V(d_int, d) <+ 0;
    end
    if (!(rg>0)) begin
      V(g_int, g) <+ 0;
    end
    if (!(rs>0)) begin
      V(s_int, s) <+ 0;
    end
    if (!(rb>0)) begin
      V(bd_int, s) <+ 0;
    end
    if (!(VDMOSthermal&&$param_given(rthjc))) begin
      V(t) <+ 0;
      V(tc) <+ 0;
    end

    // Tag: temp, function: VDMOStemp
    if (!VDMOStempGiven) begin
      VDMOStemp = $temperature+dtemp;
    end
    // inlined: VDMOStempUpdate(model, here, here->VDMOStemp, ckt)
    fact1 = VDMOStnom/(27.0+273.15);
    vtnom = VDMOStnom*CONSTKoverQ;
    kt1 = 1.38064852e-23*VDMOStnom;
    egfet1 = 1.16-7.02e-4*VDMOStnom*VDMOStnom/(VDMOStnom+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);
    xfc = ln(1-VDIOdepletionCapCoeff);
    dt = VDMOStemp-VDMOStnom;
    ratio = VDMOStemp/VDMOStnom;
    VDMOStTransconductance = VDMOStransconductance*pow(ratio, VDMOSmu);
    VDMOStVth = VDMOSvth0+tcvth*dt;
    VDMOStksubthres = ksubthres*(1.0+tksubthres1*dt+tksubthres2*dt*dt);
    if ($param_given(texp0)) begin
      i_VDMOSdrainResistance = rd*pow(ratio, texp0);
    end else begin
      i_VDMOSdrainResistance = rd*(1.0+trd1*dt+trd2*dt*dt);
    end
    if (rg>0) begin
      VDMOSgateConductance = 1.0/rg/(1.0+trg1*dt+trg2*dt*dt);
    end else begin
      VDMOSgateConductance = 0.0;
    end
    if (rs>0) begin
      VDMOSsourceConductance = 1.0/rs/(1.0+trs1*dt+trs2*dt*dt);
    end else begin
      VDMOSsourceConductance = 0.0;
    end
    if (VDMOSqsGiven) begin
      i_VDMOSqsResistance = rq*pow(ratio, texp1);
    end
    vt = VDMOStemp*CONSTKoverQ;
    fact2 = VDMOStemp/(27.0+273.15);
    kt = VDMOStemp*1.38064852e-23;
    egfet = 1.16-7.02e-4*VDMOStemp*VDMOStemp/(VDMOStemp+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);
    phio = (phi-pbfact1)/fact1;
    VDMOStPhi = fact2*phio+pbfact;
    factor = 1.0+VDIOgradCoeffTemp1*dt+VDIOgradCoeffTemp2*dt*dt;
    VDIOtGradingCoeff = m*factor;
    pbo = (vj-pbfact1)/fact1;
    gmaold = (vj-pbo)/pbo;
    VDIOtJctCap = cjo/(1+VDIOtGradingCoeff*(400e-6*(VDMOStnom-(27.0+273.15))-gmaold));
    VDIOtJctPot = pbfact+fact2*pbo;
    gmanew = (VDIOtJctPot-pbo)/pbo;
    VDIOtJctCap = VDIOtJctCap * (1+VDIOtGradingCoeff*(400e-6*(VDMOStemp-(27.0+273.15))-gmanew));
    vte = n*vt;
    arg1 = (VDMOStemp/VDMOStnom-1)*VDIOeg/vte;
    arg1_dT = VDIOeg/(vte*VDMOStnom)-VDIOeg*(VDMOStemp/VDMOStnom-1)/(vte*VDMOStemp);
    arg2 = xti/n*ln(VDMOStemp/VDMOStnom);
    arg2_dT = xti/n/VDMOStemp;
    VDIOtSatCur = VDIOjctSatCur*exp(arg1+arg2);
    VDIOtSatCur_dT = VDIOjctSatCur*exp(arg1+arg2)*(arg1_dT+arg2_dT);
    VDIOtF1 = VDIOtJctPot*(1-exp((1-VDIOtGradingCoeff)*xfc))/(1-VDIOtGradingCoeff);
    VDIOtDepCap = VDIOdepletionCapCoeff*VDIOtJctPot;
    VDIOtVcrit = vte*ln(vte/(CONSTroot2*VDIOtSatCur));
    if (VDIOtDepCap>2.5) begin
      VDIOtJctPot = 2.5/n;
      VDIOtDepCap = n*VDIOtJctPot;
      $warning("junction potential VJ too large, limited to %f", VDIOtJctPot);
    end
    if ($param_given(bv)) begin
      tBreakdownVoltage = abs(bv);
      cbv = ibv;
      if (cbv<VDIOtSatCur*tBreakdownVoltage/vt) begin
        xbv = tBreakdownVoltage;
      end else begin
        tol = reltol*cbv;
        xbv = tBreakdownVoltage-nbv*vt*ln(1+cbv/VDIOtSatCur);
        for(iter = 0; iter<25; iter = iter+1) begin
          xbv = tBreakdownVoltage-nbv*vt*ln(cbv/VDIOtSatCur+1-xbv/vt);
          xcbv = VDIOtSatCur*(exp((tBreakdownVoltage-xbv)/(nbv*vt))-1+xbv/vt);
          if (abs(xcbv-cbv)<=tol) begin
            iter = 25;
          end
        end
        if (abs(xcbv-cbv)<=tol) begin
          $warning("unable to match forward and reverse diode regions: bv = %g, ibv = %g", xbv, xcbv);
        end
      end
      VDIOtBrkdwnV = xbv;
    end
    factor = 1.0+VDIOtranTimeTemp1*dt+VDIOtranTimeTemp2*dt*dt;
    VDIOtTransitTime = tt*factor;
    factor = 1.0+trb1*dt+trb2*dt*dt;
    VDIOtConductance = VDIOconductance/factor;
    VDIOtConductance_dT = -VDIOconductance*(trb1+trb2*dt)/(factor*factor);
    VDIOtF2 = exp((1+VDIOtGradingCoeff)*xfc);
    VDIOtF3 = 1-VDIOdepletionCapCoeff*(1+VDIOtGradingCoeff);
    // end of inlined: VDMOStempUpdate(model, here, here->VDMOStemp, ckt)

    // Tag: load, function: VDMOSload
    xfact = 0.0;
    capgs = 0.0;
    capgd = 0.0;
    capth = 0.0;
    dIds_dT = 0.0;
    Vrd = 0.0;
    dIth_dVrd = 0.0;
    dIrd_dT = 0.0;
    dgdrain_dT = 0.0;
    Vrs = 0.0;
    dIth_dVrs = 0.0;
    dIrs_dT = 0.0;
    dgsource_dT = 0.0;
    ceqqth = 0.0;
    gcTt = 0.0;
    load_cgdmin = cgdmin;
    load_cgdmax = cgdmax;
    load_a = a;
    load_cgs = cgs;
    Temp = VDMOStemp;
    selfheat = VDMOSthermal&&$param_given(rthjc);
    limited = 0;
    /* Previous vgs, vds */
    VDMOSvgs = VDMOStype * $limit(V(g_int, s_int), DEVlimitOldGet);
    VDMOSvds = VDMOStype * $limit(V(d_int, s_int), DEVlimitOldGet);
    VDMOSdelTemp = $limit(V(t), DEVlimitOldGet);
    
    /* vgdo value */
    vgdo = VDMOSvgs - VDMOSvds;
    
    /* Von value */
    load_von = VDMOStVth * VDMOStype;
    load_vgs = VDMOStype*V(g_int, s_int);
    load_vds = VDMOStype*V(d_int, s_int);
    if (selfheat) begin
      delTemp = V(t);
    end else begin
      delTemp = 0.0;
    end
    vgd = load_vgs-load_vds;
    vgdo = VDMOSvgs-VDMOSvds;
    load_von = VDMOStype*VDMOSvon;
    if (VDMOSvds>=0) begin
      load_vgs = DEVfetlim(load_vgs, VDMOSvgs, load_von, limited);
      load_vds = load_vgs-vgd;
      load_vds = DEVlimvds(load_vds, VDMOSvds, limited);
    end else begin
      vgd = DEVfetlim(vgd, vgdo, load_von, limited);
      load_vds = load_vgs-vgd;
      if (!oldlimit) begin
        load_vds = -DEVlimvds(-load_vds, -VDMOSvds, limited);
      end
      load_vgs = vgd+load_vds;
    end
    if (selfheat) begin
      delTemp = DEVlimitlog(delTemp, VDMOSdelTemp, 10, limited);
    end else begin
      delTemp = 0.0;
    end
    if (initialize_limiting()) begin
        load_vgs = VDMOStype * VDMOSvth0 + 0.1;
        load_vds = 0;
        delTemp = $temperature - 273.15;
    end
    
    load_vgs = VDMOStype * $limit(V(g_int, s_int), DEVlimitNewSet, VDMOStype * load_vgs, limited);
    load_vds = VDMOStype * $limit(V(d_int, s_int), DEVlimitNewSet, VDMOStype * load_vds, limited);
    delTemp = $limit(V(t), DEVlimitNewSet, delTemp, limited);
    if (selfheat) begin
      Temp = delTemp+273.15;
      // inlined: VDMOStempUpdate(model, here, Temp, ckt)
      load_fact1 = VDMOStnom/(27.0+273.15);
      load_vtnom = VDMOStnom*CONSTKoverQ;
      load_kt1 = 1.38064852e-23*VDMOStnom;
      load_egfet1 = 1.16-7.02e-4*VDMOStnom*VDMOStnom/(VDMOStnom+1108);
      load_arg1 = -load_egfet1/(load_kt1+load_kt1)+1.1150877/(1.38064852e-23*(27.0+273.15+(27.0+273.15)));
      load_pbfact1 = -2*load_vtnom*(1.5*ln(load_fact1)+1.6021766208e-19*load_arg1);
      load_xfc = ln(1-VDIOdepletionCapCoeff);
      load_dt = Temp-VDMOStnom;
      load_ratio = Temp/VDMOStnom;
      VDMOStTransconductance = VDMOStransconductance*pow(load_ratio, VDMOSmu);
      VDMOStVth = VDMOSvth0+tcvth*load_dt;
      VDMOStksubthres = ksubthres*(1.0+tksubthres1*load_dt+tksubthres2*load_dt*load_dt);
      if ($param_given(texp0)) begin
        i_VDMOSdrainResistance = rd*pow(load_ratio, texp0);
      end else begin
        i_VDMOSdrainResistance = rd*(1.0+trd1*load_dt+trd2*load_dt*load_dt);
      end
      if (rg>0) begin
        VDMOSgateConductance = 1.0/rg/(1.0+trg1*load_dt+trg2*load_dt*load_dt);
      end else begin
        VDMOSgateConductance = 0.0;
      end
      if (rs>0) begin
        VDMOSsourceConductance = 1.0/rs/(1.0+trs1*load_dt+trs2*load_dt*load_dt);
      end else begin
        VDMOSsourceConductance = 0.0;
      end
      if (VDMOSqsGiven) begin
        i_VDMOSqsResistance = rq*pow(load_ratio, texp1);
      end
      load1_vt = Temp*CONSTKoverQ;
      load_fact2 = Temp/(27.0+273.15);
      load_kt = Temp*1.38064852e-23;
      load_egfet = 1.16-7.02e-4*Temp*Temp/(Temp+1108);
      load1_arg = -load_egfet/(load_kt+load_kt)+1.1150877/(1.38064852e-23*(27.0+273.15+(27.0+273.15)));
      load_pbfact = -2*load1_vt*(1.5*ln(load_fact2)+1.6021766208e-19*load1_arg);
      load_phio = (phi-load_pbfact1)/load_fact1;
      VDMOStPhi = load_fact2*load_phio+load_pbfact;
      load_factor = 1.0+VDIOgradCoeffTemp1*load_dt+VDIOgradCoeffTemp2*load_dt*load_dt;
      VDIOtGradingCoeff = m*load_factor;
      load_pbo = (vj-load_pbfact1)/load_fact1;
      load_gmaold = (vj-load_pbo)/load_pbo;
      VDIOtJctCap = cjo/(1+VDIOtGradingCoeff*(400e-6*(VDMOStnom-(27.0+273.15))-load_gmaold));
      VDIOtJctPot = load_pbfact+load_fact2*load_pbo;
      load_gmanew = (VDIOtJctPot-load_pbo)/load_pbo;
      VDIOtJctCap = VDIOtJctCap * (1+VDIOtGradingCoeff*(400e-6*(Temp-(27.0+273.15))-load_gmanew));
      load1_vte = n*load1_vt;
      load_arg1 = (Temp/VDMOStnom-1)*VDIOeg/load1_vte;
      load_arg1_dT = VDIOeg/(load1_vte*VDMOStnom)-VDIOeg*(Temp/VDMOStnom-1)/(load1_vte*Temp);
      load_arg2 = xti/n*ln(Temp/VDMOStnom);
      load_arg2_dT = xti/n/Temp;
      VDIOtSatCur = VDIOjctSatCur*exp(load_arg1+load_arg2);
      VDIOtSatCur_dT = VDIOjctSatCur*exp(load_arg1+load_arg2)*(load_arg1_dT+load_arg2_dT);
      VDIOtF1 = VDIOtJctPot*(1-exp((1-VDIOtGradingCoeff)*load_xfc))/(1-VDIOtGradingCoeff);
      VDIOtDepCap = VDIOdepletionCapCoeff*VDIOtJctPot;
      VDIOtVcrit = load1_vte*ln(load1_vte/(CONSTroot2*VDIOtSatCur));
      if (VDIOtDepCap>2.5) begin
        VDIOtJctPot = 2.5/n;
        VDIOtDepCap = n*VDIOtJctPot;
        $warning("junction potential VJ too large, limited to %f", VDIOtJctPot);
      end
      if ($param_given(bv)) begin
        load_tBreakdownVoltage = abs(bv);
        load_cbv = ibv;
        if (load_cbv<VDIOtSatCur*load_tBreakdownVoltage/load1_vt) begin
          load_xbv = load_tBreakdownVoltage;
        end else begin
          load_tol = reltol*load_cbv;
          load_xbv = load_tBreakdownVoltage-nbv*load1_vt*ln(1+load_cbv/VDIOtSatCur);
          for(load_iter = 0; load_iter<25; load_iter = load_iter+1) begin
            load_xbv = load_tBreakdownVoltage-nbv*load1_vt*ln(load_cbv/VDIOtSatCur+1-load_xbv/load1_vt);
            load_xcbv = VDIOtSatCur*(exp((load_tBreakdownVoltage-load_xbv)/(nbv*load1_vt))-1+load_xbv/load1_vt);
            if (abs(load_xcbv-load_cbv)<=load_tol) begin
              load_iter = 25;
            end
          end
          if (abs(load_xcbv-load_cbv)<=load_tol) begin
            $warning("unable to match forward and reverse diode regions: bv = %g, ibv = %g", load_xbv, load_xcbv);
          end
        end
        VDIOtBrkdwnV = load_xbv;
      end
      load_factor = 1.0+VDIOtranTimeTemp1*load_dt+VDIOtranTimeTemp2*load_dt*load_dt;
      VDIOtTransitTime = tt*load_factor;
      load_factor = 1.0+trb1*load_dt+trb2*load_dt*load_dt;
      VDIOtConductance = VDIOconductance/load_factor;
      VDIOtConductance_dT = -VDIOconductance*(trb1+trb2*load_dt)/(load_factor*load_factor);
      VDIOtF2 = exp((1+VDIOtGradingCoeff)*load_xfc);
      VDIOtF3 = 1-VDIOdepletionCapCoeff*(1+VDIOtGradingCoeff);
      // end of inlined: VDMOStempUpdate(model, here, Temp, ckt)
    end else begin
      Temp = VDMOStemp;
    end
    VDMOSTempSH = Temp;
    if (selfheat) begin
      Beta = VDMOStTransconductance;
      dBeta_dT = Beta*VDMOSmu/Temp;
      rd0T = i_VDMOSdrainResistance;
      if ($param_given(texp0)) begin
        drd0T_dT = rd0T*texp0/Temp;
      end else begin
        drd0T_dT = rd*(trd1+2*trd2*(Temp-VDMOStnom));
      end
      rd1T = 0.0;
      drd1T_dT = 0.0;
      if (VDMOSqsGiven) begin
        rd1T = i_VDMOSqsResistance;
        drd1T_dT = rd1T*texp1/Temp;
      end
      rsT = 1/VDMOSsourceConductance;
      drsT_dT = rs*(trs1+2*trs2*(Temp-VDMOStnom));
      dgsource_dT = -drsT_dT/(rsT*rsT);
    end else begin
      Beta = VDMOStTransconductance;
      dBeta_dT = 0.0;
      rd0T = i_VDMOSdrainResistance;
      drd0T_dT = 0.0;
      rd1T = 0.0;
      drd1T_dT = 0.0;
      if (VDMOSqsGiven) begin
        rd1T = i_VDMOSqsResistance;
      end
      dgsource_dT = 0.0;
    end
    vgd = load_vgs-load_vds;
    if (load_vds>=0) begin
      VDMOSmode = 1;
      Vds = load_vds;
      Vgs = load_vgs;
    end else begin
      VDMOSmode = -1;
      Vds = -load_vds;
      Vgs = vgd;
    end
    load_von = VDMOStVth*VDMOStype;
    vgst = ((VDMOSmode==1) ? (load_vgs) : (vgd))-load_von;
    vdsat = ((vgst>0) ? (vgst) : (0));
    slope = VDMOStksubthres;
    load_lambda = lambda;
    load_theta = theta;
    shift = subshift;
    mtr = mtriode;
    vdss = load_vds*mtr*VDMOSmode;
    t0 = 1+load_lambda*load_vds;
    t1 = 1+load_theta*vdsat;
    betap = Beta*t0/t1;
    dbetapdvgs = -Beta*load_theta*t0/(t1*t1);
    dbetapdvds = Beta*load_lambda/t1;
    t2 = exp((vgst-shift)/slope);
    vgst = slope*ln(1+t2);
    dvgstdvgs = t2/(t2+1);
    if (vgst<=vdss) begin
      cdrain = betap*vgst*vgst*0.5;
      VDMOSgm = betap*vgst*dvgstdvgs+0.5*dbetapdvgs*vgst*vgst;
      VDMOSgds = 0.5*dbetapdvds*vgst*vgst;
    end else begin
      cdrain = betap*vdss*(vgst-0.5*vdss);
      VDMOSgm = betap*vdss*dvgstdvgs+vdss*dbetapdvgs*(vgst-0.5*vdss);
      VDMOSgds = vdss*dbetapdvds*(vgst-0.5*vdss)+betap*mtr*(vgst-0.5*vdss)-0.5*vdss*betap*mtr;
    end
    if (selfheat) begin
      dvgst_dT = VDMOStype*tcvth;
      dvdsat_dT = 0.0;
      if (vgst>0) begin
        dvdsat_dT = dvgst_dT;
      end
      dt1_dT = load_theta*dvdsat_dT;
      dbetap_dT = t0*(t1*dBeta_dT-Beta*dt1_dT)/(t1*t1);
      if (vgst<=vdss) begin
        dIds_dT = 0.5*dbetap_dT*vgst*vgst+betap*vgst*dvgst_dT;
      end else begin
        dIds_dT = vdss*(dbetap_dT*vgst+betap*dvgst_dT)-dbetap_dT*vdss*0.5*vdss;
      end
    end
    VDMOSvon = VDMOStype*load_von;
    VDMOSvdsat = VDMOStype*vdsat;
    VDMOScd = VDMOSmode*cdrain;
    capmodel_scale = max(load_cgdmin, load_cgdmax);
    if (capmodel_scale==0) begin
      capmodel_scale = 0;
      capmodel_unscale = 1;
    end else begin
      capmodel_unscale = capmodel_scale;
    end
    VDMOSvgs = load_vgs;
    VDMOSvds = load_vds;
    VDMOSdelTemp = delTemp;
    if (VDMOSqsGiven&&VDMOSmode==1) begin
      vdsn = VDMOStype*V(d, s);
      load_rd = rd0T+rd1T*(vdsn/(vdsn+abs(vq)));
      drd_dT = drd0T_dT+drd1T_dT*(vdsn/(vdsn+abs(vq)));
      if (load_rd>0) begin
        VDMOSdrainConductance = 1/load_rd+lc_gmin/$mfactor;
        dgdrain_dT = -drd_dT/(load_rd*load_rd);
      end else begin
        VDMOSdrainConductance = 1/rd0T;
        dgdrain_dT = -drd0T_dT/(rd0T*rd0T);
      end
    end else if (rd0T>0) begin
      VDMOSdrainConductance = 1/rd0T;
      dgdrain_dT = -drd0T_dT/(rd0T*rd0T);
    end
    if (selfheat) begin
      GmT = dIds_dT;
      VDMOSgmT = GmT;
    end else begin
      GmT = 0.0;
      VDMOSgmT = 0.0;
    end
    if (selfheat) begin
      VDMOSgtempg = VDMOStype*VDMOSgm*Vds;
      VDMOSgtempT = GmT*Vds;
      VDMOSgtempd = VDMOStype*(VDMOSgds*Vds+cdrain);
      VDMOScth = cdrain*Vds;
      Vrd = V(d, d_int);
      dIth_dVrd = VDMOSdrainConductance*2*Vrd;
      dIrd_dgdrain = Vrd;
      dIrd_dT = dIrd_dgdrain*dgdrain_dT;
      VDMOScth = VDMOScth + (VDMOSdrainConductance*Vrd*Vrd);
      Vrs = V(s, s_int);
      dIth_dVrs = VDMOSsourceConductance*2*Vrs;
      dIrs_dgsource = Vrs;
      dIrs_dT = dIrs_dgsource*dgsource_dT;
      VDMOScth = VDMOScth + (VDMOSsourceConductance*Vrs*Vrs);
    end
    DevCapVDMOS(vgd, load_cgdmin, load_cgdmax, load_a, load_cgs, VDMOScapgs, VDMOScapgd);
    VDMOScapth = cthj/2;
    delTemp1 = VDMOSdelTemp;
    capgs = VDMOScapgs+VDMOScapgs;
    capgd = VDMOScapgd+VDMOScapgd;
    capth = VDMOScapth+VDMOScapth;
    VDMOSqth = delTemp*capth;
    if (capgs==0) begin
      VDMOScqgs = 0;
    end
    if (capgd==0) begin
      VDMOScqgd = 0;
    end
    if (capth==0) begin
      VDMOScqth = 0;
    end
    gcgs = 0;
    ceqgs = capgs*(ddt(capmodel_scale*load_vgs)/capmodel_unscale);
    VDMOScqgs = ceqgs;
    gcgd = 0;
    ceqgd = capgd*(ddt(capmodel_scale*vgd)/capmodel_unscale);
    VDMOScqgd = ceqgd;
    gcTt = 0;
    VDMOScqth = ddt(VDMOSqth);
    ceqqth = VDMOScqth;
    if (selfheat) begin
      if (VDMOSmode>=0) begin
        GmT = VDMOStype*VDMOSgmT;
        gTtg = VDMOSgtempg;
        gTtdp = VDMOSgtempd;
        gTtt = VDMOSgtempT;
        gTtsp = -(gTtg+gTtdp);
      end else begin
        GmT = -VDMOStype*VDMOSgmT;
        gTtg = VDMOSgtempg;
        gTtsp = VDMOSgtempd;
        gTtt = VDMOSgtempT;
        gTtdp = -(gTtg+gTtsp);
      end
    end else begin
      GmT = 0.0;
      gTtg = 0.0;
      gTtdp = 0.0;
      gTtt = 0.0;
      gTtsp = 0.0;
    end
    if (VDMOSmode>=0) begin
      xnrm = 1;
      xrev = 0;
      cdreq = VDMOStype*cdrain;
    end else begin
      xnrm = 0;
      xrev = 1;
      cdreq = -VDMOStype*cdrain;
    end
    if (selfheat) begin
      vDevTemp = $temperature-273.15;
      if (VDMOStempGiven) begin
        vDevTemp = VDMOStemp-273.15;
      end
    end
    if (selfheat) begin
      gthjc = 1.0/rthjc;
      gthca = 1.0/rthca;
    end
    Ith = 0.0;
    dIth_dT = 0.0;
    dIdio_dT = 0.0;
    dIth_dVdio = 0.0;
    vrb = 0.0;
    dIrb_dT = 0.0;
    dIth_dVrb = 0.0;
    gbpr = VDIOtConductance;
    load_vt = CONSTKoverQ*Temp;
    load_vte = n*load_vt;
    vtebrk = nbv*load_vt;
    vbrknp = VDIOtBrkdwnV;
    vd = VDMOStype*V(bd_int, d);
    delvd = vd-VDIOvoltage;
    /* Previous vd */
    VDIOvoltage = VDMOStype * $limit(V(bd_int, d), DEVlimitOldGet);
    if ($param_given(bv)&&vd<((0<-vbrknp+10*vtebrk) ? (0) : (-vbrknp+10*vtebrk))) begin
      vdtemp = -(vd+vbrknp);
      vdtemp = DEVpnjlim(vdtemp, -(VDIOvoltage+vbrknp), vtebrk, VDIOtVcrit, limited);
      vd = -(vdtemp+vbrknp);
    end else begin
      vd = DEVpnjlim(vd, VDIOvoltage, load_vte, VDIOtVcrit, limited);
    end
    if (initialize_limiting()) begin
        vd = VDIOtVcrit;
    end
    
    vd  = VDMOStype * $limit(V(bd_int, d), DEVlimitNewSet, VDMOStype * vd, limited);
    if (vd>=-3*load_vte) begin
      evd = exp(vd/load_vte);
      cdb = VDIOtSatCur*(evd-1);
      dIdio_dT = VDIOtSatCur_dT*(evd-1)-VDIOtSatCur*vd*evd/(load_vte*Temp);
      gdb = VDIOtSatCur*evd/load_vte;
    end else if (!$param_given(bv)||vd>=-vbrknp) begin
      load_arg = 3*load_vte/(vd*CONSTe);
      arg3 = load_arg*load_arg*load_arg;
      darg3_dT = 3*arg3/Temp;
      cdb = -VDIOtSatCur*(1+arg3);
      dIdio_dT = -VDIOtSatCur_dT*(arg3+1)-VDIOtSatCur*darg3_dT;
      gdb = VDIOtSatCur*3*load_arg/vd;
    end else begin
      evrev = exp(-(vbrknp+vd)/vtebrk);
      cdb = -VDIOtSatCur*evrev;
      dIdio_dT = VDIOtSatCur*(-vbrknp-vd)*evrev/vtebrk/Temp-VDIOtSatCur_dT*evrev;
      gdb = VDIOtSatCur*evrev/vtebrk;
    end
    cd = cdb+lc_gmin/$mfactor*vd;
    gd = gdb+lc_gmin/$mfactor;
    czero = VDIOtJctCap;
    if (vd<VDIOtDepCap) begin
      load_arg = 1-vd/VDIOtJctPot;
      sarg = exp(-VDIOtGradingCoeff*ln(load_arg));
      deplcharge = VDIOtJctPot*czero*(1-load_arg*sarg)/(1-VDIOtGradingCoeff);
      deplcap = czero*sarg;
    end else begin
      czof2 = czero/VDIOtF2;
      deplcharge = czero*VDIOtF1+czof2*(VDIOtF3*(vd-VDIOtDepCap)+VDIOtGradingCoeff/(VDIOtJctPot+VDIOtJctPot)*(vd*vd-VDIOtDepCap*VDIOtDepCap));
      deplcap = czof2*(VDIOtF3+VDIOtGradingCoeff*vd/VDIOtJctPot);
    end
    diffcharge = VDIOtTransitTime*cdb;
    VDIOcapCharge = diffcharge+deplcharge;
    diffcap = VDIOtTransitTime*gdb;
    capd = diffcap+deplcap;
    VDIOcap = capd;
    geq = 0;
    VDIOcapCurrent = ddt(VDIOcapCharge);
    ceq = VDIOcapCurrent;
    gd = gd+geq;
    cd = cd+VDIOcapCurrent;
    VDIOvoltage = vd;
    VDIOcurrent = cd;
    VDIOconduct = gd;
    VDIOdIdio_dT = dIdio_dT;
    if (selfheat) begin
      vrb = V(s, bd_int);
      Ith = vd*cd;
      dIth_dVdio = cd+vd*gd;
      dIth_dVrb = vrb*gbpr;
      dIrb_dgbpr = vrb;
      dIrb_dT = dIrb_dgbpr*VDIOtConductance_dT;
      dIth_dIrb = vrb;
      dIth_dT = dIth_dIrb*dIrb_dT+dIdio_dT*vd;
    end
    cdeq = cd;
    /* Internal MOSFET */
    I(g_int) <+ VDMOStype * (ceqgs + ceqgd);
    I(d_int) <+ cdreq - VDMOStype * ceqgd; 
    I(s_int) <+ -cdreq - VDMOStype * ceqgs;
    
    /* Body diode */
    I(bd_int, s) <+ V(bd_int, s) * gbpr;
    I(bd_int, d) <+ VDMOStype * cdeq;
    
    /* Series resistors */
    I(g_int, g) <+ V(g_int, g) * VDMOSgateConductance;
    I(d_int, d) <+ V(d_int, d) * VDMOSdrainConductance;
    I(s_int, s) <+ V(s_int, s) * VDMOSsourceConductance;
    
    /* D-S conductance */
    I(d, s) <+ V(d, s) * VDMOSdsConductance;
    
    /* Thermal model */
    V(tbr) <+ vDevTemp * $simparam("sourceScaleFactor", 1.0);
    I(t, tc) <+ V(t, tc) * gthjc;
    I(tc, t_int) <+ V(tc, t_int) * gthca;
    I(t) <+ -VDMOScth + ceqqth;

    // Tag: noise, function: VDMOSnoise
    if (VDMOSoxideCapFactor==0.0) begin
      coxSquared = 3.9*8.854214871e-12/1e-7;
    end else begin
      coxSquared = VDMOSoxideCapFactor;
    end
    coxSquared = coxSquared * (coxSquared);
    if (VDMOSthermal&&$param_given(rthjc)) begin
      tempRatioSH = VDMOSTempSH/$temperature;
    end else begin
      tempRatioSH = 1.0;
    end
    I(d_int, d) <+ white_noise(4*(VACONSTboltz*($temperature*(VDMOSdrainConductance*tempRatioSH))), "rd");
    I(s_int, s) <+ white_noise(4*(VACONSTboltz*($temperature*(VDMOSsourceConductance*tempRatioSH))), "rs");
    I(d_int, s_int) <+ white_noise(4*(VACONSTboltz*($temperature*(2.0/3.0*abs(VDMOSgm)*tempRatioSH))), "id");
    I(d_int, s_int) <+ flicker_noise(((VDMOScd>=0) ? (1) : (-1))*(kf*exp(af*ln(((abs(VDMOScd)>1E-38) ? (abs(VDMOScd)) : (1E-38))))/coxSquared), 1, "flicker");

    // Instance output variable computation
    vgs = VDMOSvgs;
    vds = VDMOSvds;
    c_gs = 2*VDMOScapgs;
    c_gd = 2*VDMOScapgd;
    cds = VDIOcap;
    von = VDMOSvon;
    sourceconductance = VDMOSsourceConductance;
    drainconductance = VDMOSdrainConductance;
    gm = VDMOSgm;
    gds = VDMOSgds;
    cqgs = VDMOScqgs;
    cqgd = VDMOScqgd;

    // Model output variable computation

  end
endmodule
