// 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:39:31 UTC
// Module name  : sp_bjt
// Model variant: full
// 
// 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 NPN (1) or PNP (-1). 
//   
//   Model parameter subs chooses whether the device is vertical 
//   (1, default) or lateral (-1). 
//   
//   Output variable is and csub were renamed to isubs and csubs due to 
//   conflicts with model parameters is and csub. 
//   
//   Input instance parameter ic was removed. Initial conditions should be 
//   set via the simulator's initial conditions mechanism. 
//   
//   Output variables geqcb, geqbx, and gcsub were removed as equivalent 
//   conductances of capacitors are not available in Verilog-A. 
//   
//   Output variables for accessing node numbers (*node) were removed. 
//   
//   The Weil-McNamee excess phase model has been implemented using the 
//   McAndrew-Huszka-Coram approach. The model adds two internal nodes (xf1 and 
//   xf2) when the excess phase parameter (ptf) is nonzero.
//
//   This model exposes output variables that can lead to creation of
//   implicit equations by OpenVAF.
//
//   Full noise model used (appropriate for all noise analyses).
//   OpenVAF may create implicit equations.
//
// CKTcircuit member mapping
//    CKTabstol -> ('curtol', '$simparam("abstol", 1e-12)', 'const')
//      CKTgmin -> ('gmin', '$simparam("gmin", 1e-12)')
//   CKTnomTemp -> ('tnom', '($simparam("tnom", 27)+273.15)')
//      CKTtemp -> (None, '$temperature')
//   CKTvoltTol -> ('voltol', '$simparam("vntol", 1e-6)', 'const')
//
// Source files and copyright information:
//   bjt/bjtdefs.h
//   bjt/bjt.c
//   bjt/bjtparam.c
//   bjt/bjtmpar.c
//     Copyright 1990 Regents of the University of California.  All rights reserved.
//     Author: 1985 Thomas L. Quarles
//   bjt/bjtask.c
//     Copyright 1990 Regents of the University of California.  All rights reserved.
//     Author: 1985 Mathew Lew and Thomas L. Quarles
//   bjt/bjtmask.c
//     Copyright 1990 Regents of the University of California.  All rights reserved.
//     Author: 1987 Mathew Lew and Thomas L. Quarles
//   bjt/bjtsetup.c
//   bjt/bjttemp.c
//   bjt/bjtload.c
//     Copyright 1990 Regents of the University of California.  All rights reserved.
//     Author: 1985 Thomas L. Quarles
//     Modified: 2000 AlansFixes
//   bjt/bjtnoise.c
//     Copyright 1990 Regents of the University of California.  All rights reserved.
//     Author: 1987 Gary W. Ng

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

module sp_bjt(c, b, e, sub);
  inout c, b, e, sub;
  electrical c, b, e, sub, cx_int, c_int, b_int, e_int, sub_con, xf1, xf2;

  // Instance parameters
  (* desc = "(Emitter) Area factor", type = "instance" *) parameter real area = 1.0;
  (* desc = "Base area factor", type = "instance" *) parameter real areab = 0;
  (* desc = "Collector area factor", type = "instance" *) parameter real areac = 0;
  (* desc = "instance temperature", type = "instance" *) parameter real temp = 0;
  (* desc = "instance temperature delta from circuit", type = "instance" *) parameter real dtemp = 0.0;

  // Model parameters
  (* desc = "N-type (1) or P-type (-1)" *) parameter integer type = 1;
  (* desc = "vertical (1) or lateral (-1)" *) parameter integer subs = 1;
  (* desc = "Parameter measurement temperature" *) parameter real tnom = 0;
  aliasparam tref = tnom;
  (* desc = "Saturation Current" *) parameter real is = 1e-16;
  (* desc = "Base-Emitter saturation Current" *) parameter real ibe = 0.0;
  (* desc = "Base-Collector saturation Current" *) parameter real ibc = 0.0;
  (* desc = "Ideal forward beta" *) parameter real bf = 100;
  (* desc = "Forward emission coefficient" *) parameter real nf = 1;
  (* desc = "Forward Early voltage" *) parameter real vaf = 0;
  aliasparam va = vaf;
  (* desc = "Forward beta roll-off corner current" *) parameter real ikf = 0;
  aliasparam ik = ikf;
  (* desc = "B-E leakage saturation current" *) parameter real ise = 0;
  aliasparam c2 = ise;
  (* desc = "B-E leakage emission coefficient" *) parameter real ne = 1.5;
  (* desc = "Ideal reverse beta" *) parameter real br = 1;
  (* desc = "Reverse emission coefficient" *) parameter real nr = 1;
  (* desc = "Reverse Early voltage" *) parameter real var = 0;
  aliasparam vb = var;
  (* desc = "reverse beta roll-off corner current" *) parameter real ikr = 0;
  (* desc = "B-C leakage saturation current" *) parameter real isc = 0;
  aliasparam c4 = isc;
  (* desc = "B-C leakage emission coefficient" *) parameter real nc = 2;
  (* desc = "Zero bias base resistance" *) parameter real rb = 0;
  (* desc = "Current for base resistance=(rb+rbm)/2" *) parameter real irb = 0;
  (* desc = "Minimum base resistance" *) parameter real rbm = 0;
  (* desc = "Emitter resistance" *) parameter real re = 0;
  (* desc = "Collector resistance" *) parameter real rc = 0;
  (* desc = "Zero bias B-E depletion capacitance" *) parameter real cje = 0;
  (* desc = "B-E built in potential" *) parameter real vje = 0.75;
  aliasparam pe = vje;
  (* desc = "B-E junction grading coefficient" *) parameter real mje = 0.33;
  aliasparam me = mje;
  (* desc = "Ideal forward transit time" *) parameter real tf = 0;
  (* desc = "Coefficient for bias dependence of TF" *) parameter real xtf = 0;
  (* desc = "Voltage giving VBC dependence of TF" *) parameter real vtf = 0;
  (* desc = "High current dependence of TF" *) parameter real itf = 0;
  (* desc = "Excess phase" *) parameter real ptf = 0;
  (* desc = "Zero bias B-C depletion capacitance" *) parameter real cjc = 0;
  (* desc = "B-C built in potential" *) parameter real vjc = 0.75;
  aliasparam pc = vjc;
  (* desc = "B-C junction grading coefficient" *) parameter real mjc = 0.33;
  aliasparam mc = mjc;
  (* desc = "Fraction of B-C cap to internal base" *) parameter real xcjc = 1;
  (* desc = "Ideal reverse transit time" *) parameter real tr = 0;
  (* desc = "Zero bias Substrate capacitance" *) parameter real cjs = 0;
  aliasparam csub = cjs;
  aliasparam ccs = cjs;
  (* desc = "Substrate junction built in potential" *) parameter real vjs = 0.75;
  aliasparam ps = vjs;
  (* desc = "Substrate junction grading coefficient" *) parameter real mjs = 0;
  aliasparam ms = mjs;
  (* desc = "Forward and reverse beta temp. exp." *) parameter real xtb = 0;
  (* desc = "Energy gap for IS temp. dependency" *) parameter real eg = 1.11;
  (* desc = "Temp. exponent for IS" *) parameter real xti = 3;
  (* desc = "Forward bias junction fit parameter" *) parameter real fc = 0;
  (* desc = "Flicker Noise Coefficient" *) parameter real kf = 0;
  (* desc = "Flicker Noise Exponent" *) parameter real af = 1;
  (* desc = "Substrate Jct. Saturation Current" *) parameter real iss = 0.0;
  (* desc = "Substrate current emission coefficient" *) parameter real ns = 1.0;
  (* desc = "Intrinsic coll. resistance" *) parameter real rco = 0;
  (* desc = "Epi drift saturation voltage" *) parameter real vo = 10.0;
  (* desc = "Epi doping parameter" *) parameter real gamma = 1e-11;
  (* desc = "Epi Charge parameter" *) parameter real qco = 0.0;
  (* desc = "Temperature equation selector" *) parameter integer tlev = 0;
  (* desc = "Temperature equation selector" *) parameter integer tlevc = 0;
  (* desc = "BF 1. temperature coefficient" *) parameter real tbf1 = 0.0;
  (* desc = "BF 2. temperature coefficient" *) parameter real tbf2 = 0.0;
  (* desc = "BR 1. temperature coefficient" *) parameter real tbr1 = 0.0;
  (* desc = "BR 2. temperature coefficient" *) parameter real tbr2 = 0.0;
  (* desc = "IKF 1. temperature coefficient" *) parameter real tikf1 = 0.0;
  (* desc = "IKF 2. temperature coefficient" *) parameter real tikf2 = 0.0;
  (* desc = "IKR 1. temperature coefficient" *) parameter real tikr1 = 0.0;
  (* desc = "IKR 2. temperature coefficient" *) parameter real tikr2 = 0.0;
  (* desc = "IRB 1. temperature coefficient" *) parameter real tirb1 = 0.0;
  (* desc = "IRB 2. temperature coefficient" *) parameter real tirb2 = 0.0;
  (* desc = "NC 1. temperature coefficient" *) parameter real tnc1 = 0.0;
  (* desc = "NC 2. temperature coefficient" *) parameter real tnc2 = 0.0;
  (* desc = "NE 1. temperature coefficient" *) parameter real tne1 = 0.0;
  (* desc = "NE 2. temperature coefficient" *) parameter real tne2 = 0.0;
  (* desc = "NF 1. temperature coefficient" *) parameter real tnf1 = 0.0;
  (* desc = "NF 2. temperature coefficient" *) parameter real tnf2 = 0.0;
  (* desc = "NR 1. temperature coefficient" *) parameter real tnr1 = 0.0;
  (* desc = "NR 2. temperature coefficient" *) parameter real tnr2 = 0.0;
  (* desc = "RB 1. temperature coefficient" *) parameter real trb1 = 0.0;
  aliasparam trb = trb1;
  (* desc = "RB 2. temperature coefficient" *) parameter real trb2 = 0.0;
  (* desc = "RC 1. temperature coefficient" *) parameter real trc1 = 0.0;
  aliasparam trc = trc1;
  (* desc = "RC 2. temperature coefficient" *) parameter real trc2 = 0.0;
  (* desc = "RE 1. temperature coefficient" *) parameter real tre1 = 0.0;
  aliasparam tre = tre1;
  (* desc = "RE 2. temperature coefficient" *) parameter real tre2 = 0.0;
  (* desc = "RBM 1. temperature coefficient" *) parameter real trm1 = 0.0;
  (* desc = "RBM 2. temperature coefficient" *) parameter real trm2 = 0.0;
  (* desc = "VAF 1. temperature coefficient" *) parameter real tvaf1 = 0.0;
  (* desc = "VAF 2. temperature coefficient" *) parameter real tvaf2 = 0.0;
  (* desc = "VAR 1. temperature coefficient" *) parameter real tvar1 = 0.0;
  (* desc = "VAR 2. temperature coefficient" *) parameter real tvar2 = 0.0;
  (* desc = "CJC temperature coefficient" *) parameter real ctc = 0.0;
  (* desc = "CJE temperature coefficient" *) parameter real cte = 0.0;
  (* desc = "CJS temperature coefficient" *) parameter real cts = 0.0;
  (* desc = "VJC temperature coefficient" *) parameter real tvjc = 0.0;
  (* desc = "VJE temperature coefficient" *) parameter real tvje = 0.0;
  (* desc = "VJS temperature coefficient" *) parameter real tvjs = 0.0;
  (* desc = "ITF 1. temperature coefficient" *) parameter real titf1 = 0.0;
  (* desc = "ITF 2. temperature coefficient" *) parameter real titf2 = 0.0;
  (* desc = "TF 1. temperature coefficient" *) parameter real ttf1 = 0.0;
  (* desc = "TF 2. temperature coefficient" *) parameter real ttf2 = 0.0;
  (* desc = "TR 1. temperature coefficient" *) parameter real ttr1 = 0.0;
  (* desc = "TR 2. temperature coefficient" *) parameter real ttr2 = 0.0;
  (* desc = "MJE 1. temperature coefficient" *) parameter real tmje1 = 0.0;
  (* desc = "MJE 2. temperature coefficient" *) parameter real tmje2 = 0.0;
  (* desc = "MJC 1. temperature coefficient" *) parameter real tmjc1 = 0.0;
  (* desc = "MJC 2. temperature coefficient" *) parameter real tmjc2 = 0.0;
  (* desc = "MJS 1. temperature coefficient" *) parameter real tmjs1 = 0.0;
  (* desc = "MJS 2. temperature coefficient" *) parameter real tmjs2 = 0.0;
  (* desc = "NS 1. temperature coefficient" *) parameter real tns1 = 0.0;
  (* desc = "NS 2. temperature coefficient" *) parameter real tns2 = 0.0;
  (* desc = "NKF High current beta rolloff exponent" *) parameter real nkf = 0.5;
  aliasparam nk = nkf;
  (* desc = "IS 1. temperature coefficient" *) parameter real tis1 = 0.0;
  (* desc = "IS 2. temperature coefficient" *) parameter real tis2 = 0.0;
  (* desc = "ISE 1. temperature coefficient" *) parameter real tise1 = 0.0;
  (* desc = "ISE 2. temperature coefficient" *) parameter real tise2 = 0.0;
  (* desc = "ISC 1. temperature coefficient" *) parameter real tisc1 = 0.0;
  (* desc = "ISC 2. temperature coefficient" *) parameter real tisc2 = 0.0;
  (* desc = "ISS 1. temperature coefficient" *) parameter real tiss1 = 0.0;
  (* desc = "ISS 2. temperature coefficient" *) parameter real tiss2 = 0.0;
  (* desc = "Temperature equation selector" *) parameter integer quasimod = 0;
  (* desc = "Energy gap for QS temp. dependency" *) parameter real vg = 1.206;
  (* desc = "Temperature exponent of RCI" *) parameter real cn = 0;
  (* desc = "Temperature exponent of VO" *) parameter real d = 0;
  (* desc = "maximum voltage B-E junction" *) parameter real vbe_max = 1e+99;
  (* desc = "maximum voltage B-C junction" *) parameter real vbc_max = 1e+99;
  (* desc = "maximum voltage C-E branch" *) parameter real vce_max = 1e+99;
  (* desc = "maximum device power dissipation" *) parameter real pd_max = 1e+99;
  (* desc = "maximum collector current" *) parameter real ic_max = 1e+99;
  (* desc = "maximum base current" *) parameter real ib_max = 1e+99;
  (* desc = "maximum temperature" *) parameter real te_max = 1e+99;
  (* desc = "thermal resistance juntion to ambient" *) parameter real rth0 = 0;

  // Instance output variables
  (* desc = "Current at collector node" *) real ic;
  (* desc = "Current at base node" *) real ib;
  (* desc = "Emitter current" *) real ie;
  (* desc = "Substrate current" *) real isubs;
  (* desc = "B-E voltage" *) real vbe;
  (* desc = "B-C voltage" *) real vbc;
  (* desc = "Small signal transconductance" *) real gm;
  (* desc = "Small signal input conductance - pi" *) real gpi;
  (* desc = "Small signal conductance - mu" *) real gmu;
  (* desc = "Conductance from base to internal base" *) real gx;
  (* desc = "Small signal output conductance" *) real go;
  (* desc = "Internal Subs. Diode equiv. cond." *) real gdsub;
  (* desc = "Internal base to emitter capacitance" *) real cpi;
  (* desc = "Internal base to collector capacitance" *) real cmu;
  (* desc = "Base to collector capacitance" *) real cbx;
  (* desc = "Substrate capacitance" *) real csubs;
  (* desc = "Cap. due to charge storage in B-E jct." *) real cqbe;
  (* desc = "Cap. due to charge storage in B-C jct." *) real cqbc;
  (* desc = "Cap. due to charge storage in Subs. jct." *) real cqsub;
  (* desc = "Cap. due to charge storage in B-X jct." *) real cqbx;
  (* desc = "Total Capacitance in B-X junction" *) real cexbc;
  (* desc = "Charge storage B-E junction" *) real qbe;
  (* desc = "Charge storage B-C junction" *) real qbc;
  (* desc = "Charge storage Subs. junction" *) real qsub;
  (* desc = "Charge storage B-X junction" *) real qbx;
  (* desc = "Power dissipation" *) real p;

  // Model output variables
  (* desc = "Inverse early voltage:forward" *) real invearlyvoltf;
  (* desc = "Inverse early voltage:reverse" *) real invearlyvoltr;
  (* desc = "Inverse roll off - forward" *) real invrollofff;
  (* desc = "Inverse roll off - reverse" *) real invrolloffr;
  (* desc = "Collector conductance" *) real collectorconduct;
  (* desc = "Emitter conductance" *) real emitterconduct;
  (* desc = "Transit time VBC factor" *) real transtimevbcfact;
  (* desc = "Excess phase fact." *) real excessphasefactor;

  // Instance fields
  real BJTareab = 0;
  real BJTareac = 0;
  real BJTtemp = 0;
  real BJTtSatCur = 0;
  real BJTBEtSatCur = 0;
  real BJTBCtSatCur = 0;
  real BJTtBetaF = 0;
  real BJTtBetaR = 0;
  real BJTtBEleakCur = 0;
  real BJTtBCleakCur = 0;
  real BJTtBEcap = 0;
  real BJTtBEpot = 0;
  real BJTtBCcap = 0;
  real BJTtBCpot = 0;
  real BJTtSubcap = 0;
  real BJTtSubpot = 0;
  real BJTtDepCap = 0;
  real BJTtf1 = 0;
  real BJTtf4 = 0;
  real BJTtf5 = 0;
  real BJTtf2 = 0;
  real BJTtf3 = 0;
  real BJTtf6 = 0;
  real BJTtf7 = 0;
  real BJTtVcrit = 0;
  real BJTtSubVcrit = 0;
  real BJTtSubSatCur = 0;
  real BJTtcollectorConduct = 0;
  real BJTtemitterConduct = 0;
  real BJTtbaseResist = 0;
  real BJTtbaseCurrentHalfResist = 0;
  real BJTtminBaseResist = 0;
  real BJTtinvEarlyVoltF = 0;
  real BJTtinvEarlyVoltR = 0;
  real BJTtinvRollOffF = 0;
  real BJTtinvRollOffR = 0;
  real BJTtemissionCoeffF = 0;
  real BJTtemissionCoeffR = 0;
  real BJTtleakBEemissionCoeff = 0;
  real BJTtleakBCemissionCoeff = 0;
  real BJTttransitTimeHighCurrentF = 0;
  real BJTttransitTimeF = 0;
  real BJTttransitTimeR = 0;
  real BJTtjunctionExpBE = 0;
  real BJTtjunctionExpBC = 0;
  real BJTtjunctionExpSub = 0;
  real BJTtemissionCoeffS = 0;
  real BJTtintCollResist = 0;
  real BJTtepiSatVoltage = 0;
  real BJTtepiDoping = 0;
  real BJTcapbe = 0;
  real BJTcapbc = 0;
  real BJTcapsub = 0;
  real BJTcapbx = 0;
  real BJTcapbcx = 0;
  real BJTdcCollectorCurrent = 0;
  real BJTdcBaseCurrent = 0;
  real xfinput = 0;
  real xfoutput = 0;
  real tdrel = 0;
  real tdscale = 0;
  integer BJTtempGiven = 0;

  // Model fields
  integer BJTsubs = 0;
  real BJTtnom = 0;
  real BJTleakBEcurrent = 0;
  real BJTleakBCcurrent = 0;
  real BJTminBaseResist = 0;
  real BJTdepletionCapCoeff = 0;
  real BJTintCollResist = 0;
  real BJTnkf = 0;
  real BJTinvEarlyVoltF = 0;
  real BJTinvEarlyVoltR = 0;
  real BJTinvRollOffF = 0;
  real BJTinvRollOffR = 0;
  real BJTcollectorConduct = 0;
  real BJTemitterConduct = 0;
  real BJTtransitTimeVBCFactor = 0;
  real BJTexcessPhaseFactor = 0;
  real BJTf2 = 0;
  real BJTf3 = 0;
  real BJTf6 = 0;
  real BJTf7 = 0;
  real BJTtempExpRCI = 0;
  real BJTtempExpVO = 0;
  integer BJTtype = 0;
  integer BJTtnomGiven = 0;

  // States
  real BJTvbe = 0;
  real BJTvbc = 0;
  real BJTvbcx = 0;
  real BJTvrci = 0;
  real BJTcc = 0;
  real BJTcb = 0;
  real BJTgpi = 0;
  real BJTgmu = 0;
  real BJTgm = 0;
  real BJTgo = 0;
  real BJTqbe = 0;
  real BJTcqbe = 0;
  real BJTqbc = 0;
  real BJTcqbc = 0;
  real BJTqsub = 0;
  real BJTcqsub = 0;
  real BJTqbx = 0;
  real BJTcqbx = 0;
  real BJTgx = 0;
  real BJTcexbc = 0;
  real BJTgeqcb = 0;
  real BJTgcsub = 0;
  real BJTgeqbx = 0;
  real BJTvsub = 0;
  real BJTcdsub = 0;
  real BJTgdsub = 0;
  real BJTirci = 0;
  real BJTirci_Vrci = 0;
  real BJTirci_Vbci = 0;
  real BJTirci_Vbcx = 0;
  real BJTqbcx = 0;
  real BJTcqbcx = 0;
  real BJTgbcx = 0;

  // Model parameter access locals

  // Instance parameter access locals

  // Device setup locals

  // Device temperature dependence locals
  real xfc = 0;
  real vt = 0;
  real vtnom = 0;
  real ratlog = 0;
  real ratio1 = 0;
  real factlog = 0;
  real bfactor = 0;
  real factor = 0;
  real fact1 = 0;
  real fact2 = 0;
  real pbo = 0;
  real pbfact = 0;
  real gmaold = 0;
  real gmanew = 0;
  real egfet = 0;
  real arg = 0;
  real dt = 0;
  real arg1 = 0;
  real pbfact1 = 0;
  real egfet1 = 0;
  real rT = 0;
  real xvar1 = 0;
  real xvar2 = 0;
  real xvar3 = 0;

  // Device load locals
  integer limited = 0;
  real load_arg1 = 0;
  real arg2 = 0;
  real arg3 = 0;
  real load_arg = 0;
  real argtf = 0;
  real capbc = 0;
  real capbe = 0;
  real capbx = 0;
  real capsub = 0;
  real cb = 0;
  real cbc = 0;
  real cbcn = 0;
  real cbe = 0;
  real cben = 0;
  real cc = 0;
  real cdis = 0;
  real ceq = 0;
  real ceqbc = 0;
  real ceqbe = 0;
  real ceqbx = 0;
  real geqsub = 0;
  real ceqsub = 0;
  real cex = 0;
  real ctot = 0;
  real czbc = 0;
  real czbcf2 = 0;
  real czbe = 0;
  real czbef2 = 0;
  real czbx = 0;
  real czbxf2 = 0;
  real czsub = 0;
  real dqbdvc = 0;
  real dqbdve = 0;
  real evbc = 0;
  real evbcn = 0;
  real evbe = 0;
  real evben = 0;
  real f1 = 0;
  real f2 = 0;
  real f3 = 0;
  real fcpc = 0;
  real fcpe = 0;
  real gbc = 0;
  real gbcn = 0;
  real gbe = 0;
  real gben = 0;
  real gcsub = 0;
  real geq = 0;
  real geqbx = 0;
  real geqcb = 0;
  real gex = 0;
  real load_gm = 0;
  real load_gmu = 0;
  real load_go = 0;
  real load_gpi = 0;
  real load_gx = 0;
  real ovtf = 0;
  real load_pc = 0;
  real load_pe = 0;
  real load_ps = 0;
  real q1 = 0;
  real q2 = 0;
  real qb = 0;
  real rbpi = 0;
  real sarg = 0;
  real sqarg = 0;
  real load_temp = 0;
  real load_tf = 0;
  real load_tr = 0;
  real load_vbc = 0;
  real vbcx = 0;
  real load_vbe = 0;
  real vbx = 0;
  real vsub = 0;
  real load_vt = 0;
  real vtc = 0;
  real vte = 0;
  real vtn = 0;
  real vts = 0;
  real xjtf = 0;
  real xmc = 0;
  real xme = 0;
  real xms = 0;
  real load_xtf = 0;
  real evsub = 0;
  real load_gdsub = 0;
  real cdsub = 0;
  real vrci = 0;
  real Irci = 0;
  real Irci_Vrci = 0;
  real Irci_Vbci = 0;
  real Irci_Vbcx = 0;
  real Qbci = 0;
  real Qbci_Vbci = 0;
  real Qbcx = 0;
  real Qbcx_Vbcx = 0;
  real gbcx = 0;
  real cbcx = 0;
  integer ttype = 0;
  real Kbci = 0;
  real Kbci_Vbci = 0;
  real Kbcx = 0;
  real Kbcx_Vbcx = 0;
  real rKp1 = 0;
  real rKp1_Vbci = 0;
  real rKp1_Vbcx = 0;
  real load_xvar1 = 0;
  real xvar1_Vbci = 0;
  real xvar1_Vbcx = 0;
  real Vcorr = 0;
  real Vcorr_Vbci = 0;
  real Vcorr_Vbcx = 0;
  real Iohm = 0;
  real Iohm_Vrci = 0;
  real Iohm_Vbci = 0;
  real Iohm_Vbcx = 0;
  real quot = 0;
  real quot_Vrci = 0;
  real rhs_current = 0;

  // Device noise locals
  real noise_dtemp = 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 curtol = $simparam("abstol", 1e-12);
  real VACONST_tnom = ($simparam("tnom", 27)+273.15);
  real voltol = $simparam("vntol", 1e-6);

  // Evaluation constants
  real lc_gmin = 0;

  // Translated function definitions

  // 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 block
  analog begin

    // Instance field assignments
    if ($param_given(areab)) begin
      BJTareab = areab;
    end
    if ($param_given(areac)) begin
      BJTareac = areac;
    end
    if ($param_given(temp)) begin
      BJTtemp = temp+273.15;
      BJTtempGiven = 1;
    end

    // Model field assignments
    BJTtype = type;
    BJTsubs = subs;
    if ($param_given(tnom)) begin
      BJTtnom = tnom+273.15;
      BJTtnomGiven = 1;
    end
    if ($param_given(ise)) begin
      BJTleakBEcurrent = ise;
    end
    if ($param_given(isc)) begin
      BJTleakBCcurrent = isc;
    end
    if ($param_given(rbm)) begin
      BJTminBaseResist = rbm;
    end
    if ($param_given(fc)) begin
      BJTdepletionCapCoeff = fc;
    end
    if ($param_given(rco)) begin
      BJTintCollResist = rco;
    end
    if ($param_given(nkf)) begin
      BJTnkf = nkf;
    end
    if ($param_given(cn)) begin
      BJTtempExpRCI = cn;
    end
    if ($param_given(d)) begin
      BJTtempExpVO = d;
    end

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

    // Tag: setup, function: BJTsetup
    if (BJTtype!=1&&BJTtype!=-1) begin
      BJTtype = 1;
    end
    if (!$param_given(subs)||BJTsubs!=1&&BJTsubs!=-1) begin
      if (BJTtype==1) begin
        BJTsubs = 1;
      end else begin
        BJTsubs = -1;
      end
    end
    if (!$param_given(ise)) begin
      BJTleakBEcurrent = 0;
    end else if (BJTleakBEcurrent>1e-04) begin
      BJTleakBEcurrent = is*BJTleakBEcurrent;
    end
    if (!$param_given(isc)) begin
      BJTleakBCcurrent = 0;
    end else if (BJTleakBCcurrent>1e-04) begin
      BJTleakBCcurrent = is*BJTleakBCcurrent;
    end
    if (!$param_given(rco)||BJTintCollResist<0.01) begin
      BJTintCollResist = 0.01;
    end
    if (!$param_given(nkf)) begin
      BJTnkf = 0.5;
    end else if (BJTnkf>1.0) begin
      $warning("NKF has been set to its maximum value: 1.0");
      BJTnkf = 1.0;
    end
    if (!$param_given(cn)) begin
      if (BJTtype==1) begin
        BJTtempExpRCI = 2.42;
      end else begin
        BJTtempExpRCI = 2.2;
      end
    end
    if (!$param_given(d)) begin
      if (BJTtype==1) begin
        BJTtempExpVO = 0.87;
      end else begin
        BJTtempExpVO = 0.52;
      end
    end
    if (!$param_given(areab)) begin
      BJTareab = area;
    end
    if (!$param_given(areac)) begin
      BJTareac = area;
    end
    if (rc==0) begin
      V(cx_int, c) <+ 0;
    end
    if (!$param_given(rco)) begin
      V(c_int, cx_int) <+ 0;
    end
    if (rb==0) begin
      V(b_int, b) <+ 0;
    end
    if (re==0) begin
      V(e_int, e) <+ 0;
    end
    if (BJTsubs==-1) begin
      V(sub_con, b_int) <+ 0;
    end else begin
      V(sub_con, c_int) <+ 0;
    end

    // Tag: temp, function: BJTtemp
    bfactor = 1.0;
    if (!BJTtnomGiven) begin
      BJTtnom = VACONST_tnom;
    end
    vtnom = CONSTKoverQ*BJTtnom;
    fact1 = BJTtnom/(27.0+273.15);
    if (!$param_given(rbm)) begin
      BJTminBaseResist = rb;
    end
    if ($param_given(vtf)&&vtf!=0) begin
      BJTtransitTimeVBCFactor = 1/(vtf*1.44);
    end else begin
      BJTtransitTimeVBCFactor = 0;
    end
    BJTexcessPhaseFactor = ptf/(180.0/3.1415926535897932384626433832795)*tf;
    if ($param_given(fc)) begin
      if (BJTdepletionCapCoeff>0.9999) begin
        BJTdepletionCapCoeff = 0.9999;
        $warning("parameter fc limited to 0.9999");
      end
    end else begin
      BJTdepletionCapCoeff = 0.5;
    end
    xfc = ln(1-BJTdepletionCapCoeff);
    if (!BJTtempGiven) begin
      BJTtemp = $temperature+dtemp;
    end
    dt = BJTtemp-BJTtnom;
    if ($param_given(vaf)&&vaf!=0) begin
      BJTtinvEarlyVoltF = 1/(vaf*(1+tvaf1*dt+tvaf2*dt*dt));
    end else begin
      BJTtinvEarlyVoltF = 0;
    end
    if ($param_given(ikf)&&ikf!=0) begin
      BJTtinvRollOffF = 1/(ikf*(1+tikf1*dt+tikf2*dt*dt));
      BJTtinvRollOffF = BJTtinvRollOffF / (area);
    end else begin
      BJTtinvRollOffF = 0;
    end
    if ($param_given(var)&&var!=0) begin
      BJTtinvEarlyVoltR = 1/(var*(1+tvar1*dt+tvar2*dt*dt));
    end else begin
      BJTtinvEarlyVoltR = 0;
    end
    if ($param_given(ikr)&&ikr!=0) begin
      BJTtinvRollOffR = 1/(ikr*(1+tikr1*dt+tikr2*dt*dt));
      BJTtinvRollOffR = BJTtinvRollOffR / (area);
    end else begin
      BJTtinvRollOffR = 0;
    end
    if ($param_given(rc)&&rc!=0) begin
      BJTtcollectorConduct = 1/(rc*(1+trc1*dt+trc2*dt*dt));
      BJTtcollectorConduct = BJTtcollectorConduct * (area);
    end else begin
      BJTtcollectorConduct = 0;
    end
    if ($param_given(re)&&re!=0) begin
      BJTtemitterConduct = 1/(re*(1+tre1*dt+tre2*dt*dt));
      BJTtemitterConduct = BJTtemitterConduct * (area);
    end else begin
      BJTtemitterConduct = 0;
    end
    BJTtbaseResist = rb*(1+trb1*dt+trb2*dt*dt);
    BJTtbaseResist = BJTtbaseResist / (area);
    BJTtminBaseResist = BJTminBaseResist*(1+trm1*dt+trm2*dt*dt);
    BJTtminBaseResist = BJTtminBaseResist / (area);
    BJTtbaseCurrentHalfResist = irb*(1+tirb1*dt+tirb2*dt*dt);
    BJTtbaseCurrentHalfResist = BJTtbaseCurrentHalfResist * (area);
    BJTtemissionCoeffF = nf*(1+tnf1*dt+tnf2*dt*dt);
    BJTtemissionCoeffR = nr*(1+tnr1*dt+tnr2*dt*dt);
    BJTtleakBEemissionCoeff = ne*(1+tne1*dt+tne2*dt*dt);
    BJTtleakBCemissionCoeff = nc*(1+tnc1*dt+tnc2*dt*dt);
    BJTttransitTimeHighCurrentF = itf*(1+titf1*dt+titf2*dt*dt);
    BJTttransitTimeHighCurrentF = BJTttransitTimeHighCurrentF * (area);
    BJTttransitTimeF = tf*(1+ttf1*dt+ttf2*dt*dt);
    BJTttransitTimeR = tr*(1+ttr1*dt+ttr2*dt*dt);
    BJTtjunctionExpBE = mje*(1+tmje1*dt+tmje2*dt*dt);
    if (BJTtjunctionExpBE>0.999) begin
      BJTtjunctionExpBE = 0.999;
      $warning("parameter mje (including tempco) is limited to 0.999");
    end
    BJTtjunctionExpBC = mjc*(1+tmjc1*dt+tmjc2*dt*dt);
    if (BJTtjunctionExpBC>0.999) begin
      BJTtjunctionExpBC = 0.999;
      $warning("parameter mjc (including tempco) is limited to 0.999");
    end
    BJTtjunctionExpSub = mjs*(1+tmjs1*dt+tmjs2*dt*dt);
    if (BJTtjunctionExpSub>0.999) begin
      BJTtjunctionExpSub = 0.999;
      $warning("parameter mjs (including tempco) is limited to 0.999");
    end
    BJTtemissionCoeffS = ns*(1+tns1*dt+tns2*dt*dt);
    vt = BJTtemp*CONSTKoverQ;
    fact2 = BJTtemp/(27.0+273.15);
    egfet = 1.16-7.02e-4*BJTtemp*BJTtemp/(BJTtemp+1108);
    arg = -egfet/(2*1.38064852e-23*BJTtemp)+1.1150877/(1.38064852e-23*(27.0+273.15+(27.0+273.15)));
    pbfact = -2*vt*(1.5*ln(fact2)+1.6021766208e-19*arg);
    egfet1 = 1.16-7.02e-4*BJTtnom*BJTtnom/(BJTtnom+1108);
    arg1 = -egfet1/(2*1.38064852e-23*BJTtnom)+1.1150877/(1.38064852e-23*(27.0+273.15+(27.0+273.15)));
    pbfact1 = -2*vtnom*(1.5*ln(fact1)+1.6021766208e-19*arg1);
    ratlog = ln(BJTtemp/BJTtnom);
    ratio1 = BJTtemp/BJTtnom-1;
    factlog = ratio1*eg/vt+xti*ratlog;
    if (tlev==0||tlev==1) begin
      factor = exp(factlog);
      BJTtSatCur = area*is*factor;
      if ($param_given(ibe)&&$param_given(ibc)) begin
        factor = exp(factlog/nf);
        BJTBEtSatCur = area*ibe*factor;
      end else begin
        BJTBEtSatCur = BJTtSatCur;
      end
      if ($param_given(ibe)&&$param_given(ibc)) begin
        factor = exp(factlog/nr);
        BJTBCtSatCur = ibc*factor;
      end else begin
        BJTBCtSatCur = BJTtSatCur;
      end
      if ($param_given(iss)) begin
        BJTtSubSatCur = iss*factor;
      end
    end else if (tlev==3) begin
      BJTtSatCur = area*pow(is, 1+tis1*dt+tis2*dt*dt);
      if ($param_given(ibe)&&$param_given(ibc)) begin
        BJTBEtSatCur = area*pow(ibe, 1+tis1*dt+tis2*dt*dt);
      end else begin
        BJTBEtSatCur = BJTtSatCur;
      end
      if ($param_given(ibe)&&$param_given(ibc)) begin
        BJTBCtSatCur = pow(ibc, 1+tis1*dt+tis2*dt*dt);
      end else begin
        BJTBCtSatCur = BJTtSatCur;
      end
      if ($param_given(iss)) begin
        BJTtSubSatCur = pow(iss, 1+tiss1*dt+tiss2*dt*dt);
      end
    end
    if (BJTsubs==1) begin
      BJTBCtSatCur = BJTBCtSatCur * (BJTareab);
    end else begin
      BJTBCtSatCur = BJTBCtSatCur * (BJTareac);
    end
    if ($param_given(iss)) begin
      if ($param_given(ibe)&&$param_given(ibc)) begin
        if (BJTsubs==1) begin
          BJTtSubSatCur = BJTtSubSatCur * (BJTareac);
        end else begin
          BJTtSubSatCur = BJTtSubSatCur * (BJTareab);
        end
      end else begin
        BJTtSubSatCur = BJTtSubSatCur * (area);
      end
    end
    if ($param_given(rco)) begin
      if (quasimod==1) begin
        rT = BJTtemp/BJTtnom;
        BJTtintCollResist = BJTintCollResist*pow(rT, BJTtempExpRCI);
        BJTtepiSatVoltage = vo*pow(rT, BJTtempExpVO);
        xvar1 = pow(rT, xti);
        xvar2 = -vg*(1.0-rT)/vt;
        xvar3 = exp(xvar2);
        BJTtepiDoping = gamma*xvar1*xvar3;
      end else begin
        BJTtintCollResist = BJTintCollResist;
        BJTtepiSatVoltage = vo;
        BJTtepiDoping = gamma;
      end
    end
    if (tlev==0) begin
      bfactor = exp(ratlog*xtb);
    end else if (tlev==1) begin
      bfactor = 1+xtb*dt;
    end
    if ($param_given(tbf1)||$param_given(tbf2)) begin
      BJTtBetaF = bf*(1+tbf1*dt+tbf2*dt*dt);
    end else begin
      BJTtBetaF = bf*bfactor;
    end
    if ($param_given(tbr1)||$param_given(tbr2)) begin
      BJTtBetaR = br*(1+tbr1*dt+tbr2*dt*dt);
    end else begin
      BJTtBetaR = br*bfactor;
    end
    if (tlev==0||tlev==1) begin
      BJTtBEleakCur = area*BJTleakBEcurrent*exp(factlog/ne)/bfactor;
      BJTtBCleakCur = BJTleakBCcurrent*exp(factlog/nc)/bfactor;
    end else if (tlev==3) begin
      BJTtBEleakCur = area*pow(BJTleakBEcurrent, 1+tise1*dt+tise2*dt*dt);
      BJTtBCleakCur = pow(BJTleakBCcurrent, 1+tisc1*dt+tisc2*dt*dt);
    end
    if (BJTsubs==1) begin
      BJTtBCleakCur = BJTtBCleakCur * (BJTareab);
    end else begin
      BJTtBCleakCur = BJTtBCleakCur * (BJTareac);
    end
    if (tlevc==0) begin
      pbo = (vje-pbfact1)/fact1;
      gmaold = (vje-pbo)/pbo;
      BJTtBEcap = cje/(1+BJTtjunctionExpBE*(4e-4*(BJTtnom-(27.0+273.15))-gmaold));
      BJTtBEpot = fact2*pbo+pbfact;
      gmanew = (BJTtBEpot-pbo)/pbo;
      BJTtBEcap = BJTtBEcap * (1+BJTtjunctionExpBE*(4e-4*(BJTtemp-(27.0+273.15))-gmanew));
    end else if (tlevc==1) begin
      BJTtBEcap = cje*(1+cte*dt);
      BJTtBEpot = vje-tvje*dt;
    end
    BJTtBEcap = BJTtBEcap * (area);
    if (tlevc==0) begin
      pbo = (vjc-pbfact1)/fact1;
      gmaold = (vjc-pbo)/pbo;
      BJTtBCcap = cjc/(1+BJTtjunctionExpBC*(4e-4*(BJTtnom-(27.0+273.15))-gmaold));
      BJTtBCpot = fact2*pbo+pbfact;
      gmanew = (BJTtBCpot-pbo)/pbo;
      BJTtBCcap = BJTtBCcap * (1+BJTtjunctionExpBC*(4e-4*(BJTtemp-(27.0+273.15))-gmanew));
    end else if (tlevc==1) begin
      BJTtBCcap = cjc*(1+ctc*dt);
      BJTtBCpot = vjc-tvjc*dt;
    end
    if (BJTsubs==1) begin
      BJTtBCcap = BJTtBCcap * (BJTareab);
    end else begin
      BJTtBCcap = BJTtBCcap * (BJTareac);
    end
    if (tlevc==0) begin
      pbo = (vjs-pbfact1)/fact1;
      gmaold = (vjs-pbo)/pbo;
      BJTtSubcap = cjs/(1+BJTtjunctionExpSub*(4e-4*(BJTtnom-(27.0+273.15))-gmaold));
      BJTtSubpot = fact2*pbo+pbfact;
      gmanew = (BJTtSubpot-pbo)/pbo;
      BJTtSubcap = BJTtSubcap * (1+BJTtjunctionExpSub*(4e-4*(BJTtemp-(27.0+273.15))-gmanew));
    end else if (tlevc==1) begin
      BJTtSubcap = cjs*(1+cts*dt);
      BJTtSubpot = vjs-tvjs*dt;
    end
    if (BJTsubs==1) begin
      BJTtSubcap = BJTtSubcap * (BJTareac);
    end else begin
      BJTtSubcap = BJTtSubcap * (BJTareab);
    end
    BJTtDepCap = BJTdepletionCapCoeff*BJTtBEpot;
    BJTtf1 = BJTtBEpot*(1-exp((1-BJTtjunctionExpBE)*xfc))/(1-BJTtjunctionExpBE);
    BJTtf4 = BJTdepletionCapCoeff*BJTtBCpot;
    BJTtf5 = BJTtBCpot*(1-exp((1-BJTtjunctionExpBC)*xfc))/(1-BJTtjunctionExpBC);
    BJTtVcrit = vt*ln(vt/(CONSTroot2*BJTtSatCur));
    if ($param_given(iss)) begin
      BJTtSubVcrit = vt*ln(vt/(CONSTroot2*BJTtSubSatCur));
    end
    BJTtf2 = exp((1+BJTtjunctionExpBE)*xfc);
    BJTtf3 = 1-BJTdepletionCapCoeff*(1+BJTtjunctionExpBE);
    BJTtf6 = exp((1+BJTtjunctionExpBC)*xfc);
    BJTtf7 = 1-BJTdepletionCapCoeff*(1+BJTtjunctionExpBC);
    if (!$port_connected(sub)) begin
      V(sub) <+ 0;
    end
    if (BJTexcessPhaseFactor==0) begin
      V(xf1) <+ 0;
      V(xf2) <+ 0;
    end

    // Tag: load, function: BJTload
    capbx = 0;
    capsub = 0;
    vbx = 0.0;
    vsub = 0.0;
    vrci = 0.0;
    Irci = 0.0;
    Irci_Vrci = 0.0;
    Irci_Vbci = 0.0;
    Irci_Vbcx = 0.0;
    Qbci = 0.0;
    Qbci_Vbci = 0.0;
    Qbcx_Vbcx = 0.0;
    ttype = BJTtype*BJTsubs;
    load_vt = BJTtemp*CONSTKoverQ;
    gcsub = 0;
    geqbx = 0;
    geqcb = 0;
    gbcx = 0;
    cbcx = 0;
    rbpi = BJTtbaseResist-BJTtminBaseResist;
    vte = BJTtleakBEemissionCoeff*load_vt;
    vtc = BJTtleakBCemissionCoeff*load_vt;
    limited = 0;
    BJTvbe = BJTtype * $limit(V(b_int, e_int), DEVlimitOldGet);
    BJTvbc = BJTtype * $limit(V(b_int, c_int), DEVlimitOldGet);
    BJTvsub = ttype * $limit(V(sub, sub_con), DEVlimitOldGet);
    load_vbe = BJTtype*V(b_int, e_int);
    load_vbc = BJTtype*V(b_int, c_int);
    vbcx = BJTtype*V(b_int, cx_int);
    vrci = BJTtype*V(cx_int, c_int);
    vbx = BJTtype*V(b, c_int);
    vsub = ttype*V(sub, sub_con);
    load_vbe = DEVpnjlim(load_vbe, BJTvbe, load_vt, BJTtVcrit, limited);
    load_vbc = DEVpnjlim(load_vbc, BJTvbc, load_vt, BJTtVcrit, limited);
    if ($param_given(iss)) begin
      vsub = DEVpnjlim(vsub, BJTvsub, load_vt, BJTtSubVcrit, limited);
    end else begin
      vsub = DEVpnjlim(vsub, BJTvsub, load_vt, 50, limited);
    end
    vrci = load_vbc-vbcx;
    if (initialize_limiting()) begin
        load_vbe = BJTtVcrit;
        vbcx = 0;
        load_vbc = vbcx;
        vbx = 0;
        vsub = vbx;
        vrci = 0;
    end
    load_vbe = BJTtype * $limit(V(b_int, e_int), DEVlimitNewSet, BJTtype * load_vbe, limited);
    load_vbc = BJTtype * $limit(V(b_int, c_int), DEVlimitNewSet, BJTtype * load_vbc, limited);
    vsub = ttype * $limit(V(sub, sub_con), DEVlimitNewSet, ttype * vsub, limited);
    vtn = load_vt*BJTtemissionCoeffF;
    if (load_vbe>=-3*vtn) begin
      evbe = exp(load_vbe/vtn);
      cbe = BJTBEtSatCur*(evbe-1);
      gbe = BJTBEtSatCur*evbe/vtn;
    end else begin
      load_arg = 3*vtn/(load_vbe*CONSTe);
      load_arg = load_arg*load_arg*load_arg;
      cbe = -BJTBEtSatCur*(1+load_arg);
      gbe = BJTBEtSatCur*3*load_arg/load_vbe;
    end
    if (BJTtBEleakCur==0) begin
      cben = 0;
      gben = 0;
    end else if (load_vbe>=-3*vte) begin
      evben = exp(load_vbe/vte);
      cben = BJTtBEleakCur*(evben-1);
      gben = BJTtBEleakCur*evben/vte;
    end else begin
      load_arg = 3*vte/(load_vbe*CONSTe);
      load_arg = load_arg*load_arg*load_arg;
      cben = -BJTtBEleakCur*(1+load_arg);
      gben = BJTtBEleakCur*3*load_arg/load_vbe;
    end
    gben = gben + (lc_gmin);
    cben = cben + (lc_gmin*load_vbe);
    vtn = load_vt*BJTtemissionCoeffR;
    if (load_vbc>=-3*vtn) begin
      evbc = exp(load_vbc/vtn);
      cbc = BJTBCtSatCur*(evbc-1);
      gbc = BJTBCtSatCur*evbc/vtn;
    end else begin
      load_arg = 3*vtn/(load_vbc*CONSTe);
      load_arg = load_arg*load_arg*load_arg;
      cbc = -BJTBCtSatCur*(1+load_arg);
      gbc = BJTBCtSatCur*3*load_arg/load_vbc;
    end
    if (BJTtBCleakCur==0) begin
      cbcn = 0;
      gbcn = 0;
    end else if (load_vbc>=-3*vtc) begin
      evbcn = exp(load_vbc/vtc);
      cbcn = BJTtBCleakCur*(evbcn-1);
      gbcn = BJTtBCleakCur*evbcn/vtc;
    end else begin
      load_arg = 3*vtc/(load_vbc*CONSTe);
      load_arg = load_arg*load_arg*load_arg;
      cbcn = -BJTtBCleakCur*(1+load_arg);
      gbcn = BJTtBCleakCur*3*load_arg/load_vbc;
    end
    gbcn = gbcn + (lc_gmin);
    cbcn = cbcn + (lc_gmin*load_vbc);
    if ($param_given(iss)) begin
      vts = load_vt*BJTtemissionCoeffS;
      if (vsub<=-3*vts) begin
        load_arg = 3*vts/(vsub*CONSTe);
        load_arg = load_arg*load_arg*load_arg;
        load_gdsub = BJTtSubSatCur*3*load_arg/vsub+lc_gmin;
        cdsub = -BJTtSubSatCur*(1+load_arg)+lc_gmin*vsub;
      end else begin
        evsub = exp(((709.0<vsub/vts) ? (709.0) : (vsub/vts)));
        load_gdsub = BJTtSubSatCur*evsub/vts+lc_gmin;
        cdsub = BJTtSubSatCur*(evsub-1)+lc_gmin*vsub;
      end
    end else begin
      load_gdsub = lc_gmin;
      cdsub = lc_gmin*vsub;
    end
    if ($param_given(rco)) begin
      if (vrci>0.) begin
        Kbci = sqrt(1+BJTtepiDoping*exp(load_vbc/load_vt));
        Kbci_Vbci = BJTtepiDoping*exp(load_vbc/load_vt)/(2*load_vt*Kbci);
        Kbcx = sqrt(1+BJTtepiDoping*exp(vbcx/load_vt));
        Kbcx_Vbcx = BJTtepiDoping*exp(vbcx/load_vt)/(2*load_vt*Kbcx);
        rKp1 = (1+Kbci)/(1+Kbcx);
        rKp1_Vbci = Kbci_Vbci/(1+Kbci);
        rKp1_Vbcx = -(1+Kbci)*Kbcx_Vbcx/((Kbcx+1)*(Kbcx+1));
        load_xvar1 = ln(rKp1);
        xvar1_Vbci = rKp1_Vbci/rKp1;
        xvar1_Vbcx = rKp1_Vbcx/rKp1;
        Vcorr = load_vt*(Kbci-Kbcx-load_xvar1);
        Vcorr_Vbci = load_vt*(Kbci_Vbci-xvar1_Vbci);
        Vcorr_Vbcx = load_vt*(-Kbcx_Vbcx-xvar1_Vbcx);
        Iohm = (vrci+Vcorr)/BJTtintCollResist;
        Iohm_Vrci = 1/BJTtintCollResist;
        Iohm_Vbci = Vcorr_Vbci/BJTtintCollResist;
        Iohm_Vbcx = Vcorr_Vbcx/BJTtintCollResist;
        quot = 1+abs(vrci)/BJTtepiSatVoltage;
        quot_Vrci = vrci/(BJTtepiSatVoltage*abs(vrci));
        Irci = Iohm/quot+lc_gmin*vrci;
        Irci_Vrci = Iohm_Vrci/quot-Iohm*quot_Vrci/(quot*quot)+lc_gmin;
        Irci_Vbci = Iohm_Vbci/quot;
        Irci_Vbcx = Iohm_Vbcx/quot;
        Qbci = qco*Kbci;
        Qbci_Vbci = qco*Kbci_Vbci;
        Qbcx = qco*Kbcx;
        Qbcx_Vbcx = qco*Kbcx_Vbcx;
        BJTqbcx = Qbcx;
        BJTcapbcx = Qbcx_Vbcx;
      end else begin
        Irci = vrci/BJTtintCollResist+lc_gmin*vrci;
        Irci_Vrci = 1/BJTtintCollResist+lc_gmin;
        Irci_Vbci = 0.0;
        Irci_Vbcx = 0.0;
        Qbci = 0.0;
        Qbci_Vbci = 0.0;
        Qbcx = 0.0;
        Qbcx_Vbcx = 0.0;
        BJTqbcx = Qbcx;
        BJTcapbcx = Qbcx_Vbcx;
      end
    end
    q1 = 1/(1-BJTtinvEarlyVoltF*load_vbc-BJTtinvEarlyVoltR*load_vbe);
    if (BJTtinvRollOffF==0&&BJTtinvRollOffR==0) begin
      qb = q1;
      dqbdve = q1*qb*BJTtinvEarlyVoltR;
      dqbdvc = q1*qb*BJTtinvEarlyVoltF;
    end else begin
      q2 = BJTtinvRollOffF*cbe+BJTtinvRollOffR*cbc;
      load_arg = ((0>1+4*q2) ? (0) : (1+4*q2));
      sqarg = 1;
      if (!$param_given(nkf)) begin
        if (load_arg!=0) begin
          sqarg = sqrt(load_arg);
        end
      end else if (load_arg!=0) begin
        sqarg = pow(load_arg, BJTnkf);
      end
      qb = q1*(1+sqarg)/2;
      if (!$param_given(nkf)) begin
        dqbdve = q1*(qb*BJTtinvEarlyVoltR+BJTtinvRollOffF*gbe/sqarg);
        dqbdvc = q1*(qb*BJTtinvEarlyVoltF+BJTtinvRollOffR*gbc/sqarg);
      end else begin
        dqbdve = q1*(qb*BJTtinvEarlyVoltR+BJTtinvRollOffF*gbe*2*sqarg*BJTnkf/load_arg);
        dqbdvc = q1*(qb*BJTtinvEarlyVoltF+BJTtinvRollOffR*gbc*2*sqarg*BJTnkf/load_arg);
      end
    end
    //
    // Weil-McNamee model of excess phase
    // Implemented as suggested by McAndrew, Huszka, and Coram
    // Excess group delay [s]: BJTexcessPhaseFactor = TD0 = ptf*pi/180*tf
    //
    xfinput = cbe/qb;
    tdscale = voltol/curtol;
    tdrel = 1.0;
    I(xf1) <+ (tdscale*xfinput-V(xf2))/tdrel - ddt(BJTexcessPhaseFactor*V(xf1));
    I(xf2) <+ (V(xf1)-V(xf2)/tdrel-ddt(BJTexcessPhaseFactor*V(xf2)/3));
    xfoutput = V(xf2)/tdscale;
    if (BJTexcessPhaseFactor == 0) begin
      // Zero excess phase
      cex = cbe;
      gex = gbe;
      cc = (cbe-cbc)/qb - cbc/BJTtBetaR - cbcn;
    end else begin
      // Nonzero excess phase
      cex = cbe;
      gex = gbe;
      cc = xfoutput-cbc/qb - cbc/BJTtBetaR - cbcn;
    end
    //
    // End of Weil-McNamee model
    //
    cb = cbe/BJTtBetaF+cben+cbc/BJTtBetaR+cbcn;
    load_gx = BJTtminBaseResist+rbpi/qb;
    if (BJTtbaseCurrentHalfResist!=0) begin
      load_arg1 = ((cb/BJTtbaseCurrentHalfResist>1e-9) ? (cb/BJTtbaseCurrentHalfResist) : (1e-9));
      arg2 = (-1+sqrt(1+14.59025*load_arg1))/2.4317/sqrt(load_arg1);
      load_arg1 = tan(arg2);
      load_gx = BJTtminBaseResist+3*rbpi*(load_arg1-arg2)/arg2/load_arg1/load_arg1;
    end
    if (load_gx!=0) begin
      load_gx = 1/load_gx;
    end
    load_gpi = gbe/BJTtBetaF+gben;
    load_gmu = gbc/BJTtBetaR+gbcn;
    load_go = (gbc+(cex-cbc)*dqbdvc/qb)/qb;
    load_gm = (gex-(cex-cbc)*dqbdve/qb)/qb-load_go;
    load_tf = BJTttransitTimeF;
    load_tr = BJTttransitTimeR;
    czbe = BJTtBEcap;
    load_pe = BJTtBEpot;
    xme = BJTtjunctionExpBE;
    cdis = xcjc;
    ctot = BJTtBCcap;
    czbc = ctot*cdis;
    czbx = ctot-czbc;
    load_pc = BJTtBCpot;
    xmc = BJTtjunctionExpBC;
    fcpe = BJTtDepCap;
    czsub = BJTtSubcap;
    load_ps = BJTtSubpot;
    xms = BJTtjunctionExpSub;
    load_xtf = xtf;
    ovtf = BJTtransitTimeVBCFactor;
    xjtf = BJTttransitTimeHighCurrentF;
    if (load_tf!=0&&load_vbe>0) begin
      argtf = 0;
      arg2 = 0;
      arg3 = 0;
      if (load_xtf!=0) begin
        argtf = load_xtf;
        if (ovtf!=0) begin
          argtf = argtf*exp(load_vbc*ovtf);
        end
        arg2 = argtf;
        if (xjtf!=0) begin
          load_temp = cbe/(cbe+xjtf);
          argtf = argtf*load_temp*load_temp;
          arg2 = argtf*(3-load_temp-load_temp);
        end
        arg3 = cbe*argtf*ovtf;
      end
      cbe = cbe*(1+argtf)/qb;
      gbe = (gbe*(1+arg2)-cbe*dqbdve)/qb;
      geqcb = load_tf*(arg3-cbe*dqbdvc)/qb;
    end
    if (load_vbe<fcpe) begin
      load_arg = 1-load_vbe/load_pe;
      sarg = exp(-xme*ln(load_arg));
      BJTqbe = load_tf*cbe+load_pe*czbe*(1-load_arg*sarg)/(1-xme);
      capbe = load_tf*gbe+czbe*sarg;
    end else begin
      f1 = BJTtf1;
      f2 = BJTtf2;
      f3 = BJTtf3;
      czbef2 = czbe/f2;
      BJTqbe = load_tf*cbe+czbe*f1+czbef2*(f3*(load_vbe-fcpe)+xme/(load_pe+load_pe)*(load_vbe*load_vbe-fcpe*fcpe));
      capbe = load_tf*gbe+czbef2*(f3+xme*load_vbe/load_pe);
    end
    fcpc = BJTtf4;
    f1 = BJTtf5;
    f2 = BJTtf6;
    f3 = BJTtf7;
    if (load_vbc<fcpc) begin
      load_arg = 1-load_vbc/load_pc;
      sarg = exp(-xmc*ln(load_arg));
      BJTqbc = load_tr*cbc+load_pc*czbc*(1-load_arg*sarg)/(1-xmc);
      capbc = load_tr*gbc+czbc*sarg;
    end else begin
      czbcf2 = czbc/f2;
      BJTqbc = load_tr*cbc+czbc*f1+czbcf2*(f3*(load_vbc-fcpc)+xmc/(load_pc+load_pc)*(load_vbc*load_vbc-fcpc*fcpc));
      capbc = load_tr*gbc+czbcf2*(f3+xmc*load_vbc/load_pc);
    end
    if (vbx<fcpc) begin
      load_arg = 1-vbx/load_pc;
      sarg = exp(-xmc*ln(load_arg));
      BJTqbx = load_pc*czbx*(1-load_arg*sarg)/(1-xmc);
      capbx = czbx*sarg;
    end else begin
      czbxf2 = czbx/f2;
      BJTqbx = czbx*f1+czbxf2*(f3*(vbx-fcpc)+xmc/(load_pc+load_pc)*(vbx*vbx-fcpc*fcpc));
      capbx = czbxf2*(f3+xmc*vbx/load_pc);
    end
    if (vsub<0) begin
      load_arg = 1-vsub/load_ps;
      sarg = exp(-xms*ln(load_arg));
      BJTqsub = load_ps*czsub*(1-load_arg*sarg)/(1-xms);
      capsub = czsub*sarg;
    end else begin
      BJTqsub = vsub*czsub*(1+xms*vsub/(2*load_ps));
      capsub = czsub*(1+xms*vsub/load_ps);
    end
    BJTcapbe = capbe;
    BJTcapbc = capbc;
    if ($param_given(rco)) begin
      BJTqbc = BJTqbc + (Qbci);
      BJTcapbc = BJTcapbc + (Qbci_Vbci);
      capbc = capbc + (Qbci_Vbci);
    end
    BJTcapsub = capsub;
    BJTcapbx = capbx;
    geq = 0;
    BJTcqbe = ddt(BJTqbe);
    ceq = BJTcqbe;
    load_gpi = load_gpi+geq;
    BJTdcBaseCurrent = cb;
    cb = cb+BJTcqbe;
    geq = 0;
    BJTcqbc = ddt(BJTqbc);
    ceq = BJTcqbc;
    load_gmu = load_gmu+geq;
    cb = cb+BJTcqbc;
    BJTdcCollectorCurrent = cc;
    cc = cc-BJTcqbc;
    geq = 0;
    BJTcqbcx = ddt(BJTqbcx);
    ceq = BJTcqbcx;
    gbcx = geq;
    cbcx = BJTcqbcx;
    gcsub = 0;
    BJTcqsub = ddt(BJTqsub);
    ceq = BJTcqsub;
    geqbx = 0;
    BJTcqbx = ddt(BJTqbx);
    ceq = BJTcqbx;
    BJTvbe = load_vbe;
    BJTvbc = load_vbc;
    BJTvbcx = vbcx;
    BJTvrci = vrci;
    BJTcc = cc;
    BJTcb = cb;
    BJTgpi = load_gpi;
    BJTgmu = load_gmu;
    BJTgm = load_gm;
    BJTgo = load_go;
    BJTgx = load_gx;
    BJTgeqcb = geqcb;
    BJTgcsub = gcsub;
    BJTgeqbx = geqbx;
    BJTvsub = vsub;
    BJTgdsub = load_gdsub;
    BJTcdsub = cdsub;
    BJTirci = Irci;
    BJTirci_Vrci = Irci_Vrci;
    BJTirci_Vbci = Irci_Vbci;
    BJTirci_Vbcx = Irci_Vbcx;
    geqsub = gcsub+load_gdsub;
    ceqsub = ttype*(BJTcqsub+cdsub);
    ceqbx = BJTtype*BJTcqbx;
    ceqbe = BJTtype*(cc+cb);
    ceqbc = BJTtype*-cc;
    if ($param_given(rco)) begin
      rhs_current = BJTtype*Irci;
    end
    
    // Residuals
    I(c) <+ BJTtcollectorConduct*V(c, cx_int);
    I(b) <+ load_gx*V(b, b_int)+ceqbx;
    I(e) <+ BJTtemitterConduct*V(e, e_int);
    I(sub) <+ ceqsub;
    I(cx_int) <+ BJTtcollectorConduct*V(cx_int, c)+gbcx*V(cx_int, b_int)+(rhs_current+-cbcx);
    I(c_int) <+ -(ceqbx+ceqbc)+-rhs_current;
    I(b_int) <+ gbcx*V(b_int, cx_int)+load_gx*V(b_int, b)+(-(-ceqbe-ceqbc)+cbcx);
    I(e_int) <+ BJTtemitterConduct*V(e_int, e)+-ceqbe;
    I(sub_con) <+ -ceqsub;

    // Tag: noise, function: BJTnoise
    if (BJTtempGiven) begin
      noise_dtemp = BJTtemp-$temperature+(BJTtnom-273.15);
    end else begin
      noise_dtemp = dtemp;
    end
    I(cx_int, c) <+ white_noise(4*(VACONSTboltz*(($temperature+noise_dtemp)*BJTtcollectorConduct)), "rc");
    I(b_int, b) <+ white_noise(4*(VACONSTboltz*(($temperature+noise_dtemp)*BJTgx)), "rb");
    I(e_int, e) <+ white_noise(4*(VACONSTboltz*(($temperature+noise_dtemp)*BJTtemitterConduct)), "re");
    I(c_int, e_int) <+ white_noise(2*(VACONSTcharge*BJTcc), "ic");
    I(b_int, e_int) <+ white_noise(2*(VACONSTcharge*BJTcb), "ib");
    I(b_int, e_int) <+ flicker_noise(((BJTcb>=0) ? (1) : (-1))*(kf*exp(af*ln(((abs(BJTcb)>1E-38) ? (abs(BJTcb)) : (1E-38))))), 1, "flicker");

    // Instance output variable computation
    ic = BJTcc * BJTtype;
    ib = BJTcb;
        if (BJTsubs==-1) begin // lateral
        ib = ib - BJTcdsub - BJTcqsub;
        end
        ib = ib * BJTtype;
    ie = -BJTcc -BJTcb;
        if (BJTsubs==1) begin // vertical
        ie = ie + BJTcdsub + BJTcqsub;
        end
        ie = ie * BJTtype;
    isubs = -BJTsubs * (BJTcqsub + BJTcdsub) * BJTtype;
    vbe = BJTvbe;
    vbc = BJTvbc;
    gm = BJTgm;
    gpi = BJTgpi;
    gmu = BJTgmu;
    gx = BJTgx;
    go = BJTgo;
    gdsub = BJTgdsub;
    cpi = BJTcapbe;
    cmu = BJTcapbc;
    cbx = BJTcapbx;
    csubs = BJTcapsub;
    cqbe = BJTcapbe;
    cqbc = BJTcapbc;
    cqsub = BJTcapsub;
    cqbx = BJTcapbx;
    cexbc = geqcb;
    qbe = BJTqbe;
    qbc = BJTqbc;
    qsub = BJTqsub;
    qbx = BJTqbx;
    p = abs(BJTcc*V(c, e));
    p = p + (abs(BJTcb*V(b, e)));
    p = p + (abs(BJTcdsub*V(sub_con, sub)));
    p = p + (BJTcqsub*abs(V(sub_con, sub)));

    // Model output variable computation

  end
endmodule
