/* ******************************************************************************
   **************** COPYRIGHT NOTICE(Originator: Michael Schroter)***************
   ******************************************************************************

The terms under which the HICUM/L0 software is provided are as follows:

Software is distributed as is, completely without warranty or service
support. Michael Schroter and his team members are not liable for the
condition or performance of the software.

Michael Schroter owns the copyright and grants users a perpetual,
irrevocable, worldwide, non-exclusive, royalty-free license with respect
to the software as set forth below.

Michael Schroter hereby disclaims all implied warranties.

Michael Schroter grants the users the right to modify, copy, and
redistribute the software and documentation, both within the user's
organization and externally, subject to the following restrictions.

1. The users agree not to charge for the model owner's code itself but may
   charge for additions, extensions, or support.

2. In any product based on the software, the users agree to acknowledge
   Michael Schroter who developed the model and software. This
   acknowledgment shall appear in the product documentation.

3. Redistributions to others of source code and documentation must retain
   the copyright notice, disclaimer, and list of conditions.

4. Redistributions to others in binary form must reproduce the copyright
   notice, disclaimer, and list of conditions in the documentation
   and/or other materials provided with the distribution
*/

//  HICUM Level_0 Version_1.32: A Verilog-A description
//  (A simplified version of HICUM Level2 model for BJT)
//  ## It is modified after the first version of HICUM/L0 code ##

// Minor code related changes
// 01/11: Corrected SPICE name for AHQ and ZETAIQF. Implemented a gmin between nodes CI and EI.
//        qj is now limited to positive values to improve convergence a negative bias.
//        Resolved a convergence issue for cc at calculation of voltage dependence of t0.
// 01/09: Introduction of temporary dc capacitance variable CJE_DC to call the procedure with
//        the AC and DC parameter set and assign an AC and DC result to its output variables
// 01/09: Ranges of ZE & ZEDC have been modified to a new range (0:1) from the old range (0:1]
// 12/08: gmin declaration by L. Lemaitre.
// 12/08: rth has been used instead of rth_t (dynamic variable) in the corresponding if statement
// 12/08: Macro `QJMODF has been used to compute AC as well as DC charge with corresponding AC and DC variables respectively
// 11/08: Conditional statement for calculating normalized minority charge to avoid overflow at TFH=0
// 11/08: Range of AHQ has been modified to a new range [-0.9:10] from the old range [0:10]
// 03/08: Quick Fix: Default value of TFH has been changed from infinity to zero and modification has been done to
// the default value limits to [0, inf) to include zero
// 12/06: Upper limit of FGEO is changed to infinity
//  06/06:  Thermal node "tnode" set as external
//  Flag FLSH introduced for controlling Self-heating calculation
//  all if-else blocks marked with begin-end
//  all series resistors and RTH are allowed to have a minimum value MIN_R
//  07/06: QCJMOD deleted, QJMODF introduced along with with HICJQ
//  ddx() operator used with QJMOD and QJMODF wherever needed
//  aj is kept at 2.4 except BE depletion charge
//  Substrate transistor transfer current added.
//  Gmin added to (bi,ei) and (bi,ci) branches.
//  hyperbolic smoothing used in rbi computation to avoid devide-by-zero.

// *********************************************************************************
// 06/06: Comment on NODE COLLAPSING:
// Presently this verilog code permits a minimum of 1 milli-Ohm resistance for any
// series resistance as well as for thermal resistance RTH. If any of the resistance
// values drops below this minimum value, the corresponding nodes are shorted with
// zero voltage contribution. We want the model compilers/simulators deal this
// situation in such a manner that the corresponding node is COLLAPSED.
// We expect that the simulators should permit current contribution statement
// for any branch with resistance value more than (or equal to) 1 milli-Ohm without
// any convergence problem. In fact, we wish NOT to have to use a voltage contribution
// statement in our Verilog code, except as an indication for the model compiler/simulator
// to interprete a zero branch voltage as NODE-COLLAPSING action.
// **********************************************************************************

//*****************************************************
//*****************************************************
// 08/04:(Modification by Cornelia Thiele)
// New expression for the normalized hole charge qpt and the model parameter AHQ is inserted
// The reverse Early-Effect VER is reintroduced
// A temperature dependent modeling of IQF using the model parameter ZETAIQF is included

//*****************************************************
//*****************************************************
// 11/08: Modification done at TUD
// 3 more parameters VDEDC, ZEDC, AJEDC have been introduced for DC depletion charge
// Flag FIQF has been introduced to introduce voltage dependence in the base related critical current
// ZETARTH has been introduced for temperature dependent thermal resistance
//*****************************************************

//*****************************************************
//*****************************************************
// 01/2011: Modification done at TUD
// Third order polynomial is solved for transfer current. Can be turned on by IT_MOD=1.
// Added voltage dependent Reverse Early voltage. Parameter aver describes voltage dependence.
// Parameters ZETAVER, ZETAVGBE, VGBE describes temperature dependence of VER and IQF.
// Added temperature dependence for IQFH and TFH, ALIQFH and KIQFH is used to model a second order temperature model.
// Parameter TEF_TEMP=0 turns temperature dependence for TEF0 off.
//*****************************************************


//*****************************************************
//*****************************************************
// 11/2012: Modification done at TUD
// This code contains a Verilog-A implementation of Vertical Non-Quasi-Static(NQS)
// Effects using adjunct gyrator networks. To turn on this effect please set FLNQS=1.
//
// This code contains Operating Point Information to turn-on this section please use the compiler flag CALC_OP.
//*****************************************************


/* *****************************************************
********************************************************
12/2015: Modification done at TUD and released as version 1.32
* Removal of some unrequired variables.
* Changed the type selection (npn or pnp) accordingly to HICUM/L2
   - supply a single type parameter (+1 -> npn, -1 -> pnp)
   - remove the npn and pnp parameter
* Parasitic PNP transistor now also takes the device type into account
* Usage of rth instead of rth_t for power calculation
   - compatible with node collapsing statement
******************************************************** */

//Default simulator: Spectre

`define PGIVEN(p) $param_given(p)
`ifdef insideADMS
    `define P(p) (*p*)
    `define INITIAL_MODEL @(initial_model)
`else
    `define P(p)
    `define INITIAL_MODEL
`endif

//Spectre
`include "constants.h"
`include "discipline.h"

// Comment this line, if calculation of operating point values should be omitted
`define CALC_OP

// Comment this line, if calculation of noise analysis should be omitted
//`define CALC_NOISE

`define NPN +1
`define PNP -1

`define VPT_thresh  1.0e2
`define EXPLIM      80.0
`define INF         1.0e6
`define TMAX        326.85
`define TMIN        -100.00
`define MIN_R       0.001
//`define GMIN        1.0e-12
//`define GMIN        $simparam("gmin")        //suggested by L.L
`define GMIN        $simparam("gmin",1e-12)  //suggested by L.L

`define QCMODF(vj,cj0,vd,z,aj,cjf)\
    if (cj0 > 0.0) begin\
        vf = vd*(1.0-exp(-ln(aj)/z));\
        xvf = (vf-vj)/VT;\
        xvf2 = sqrt(xvf*xvf+1.921812);\
        v_j = vf-VT*(xvf+xvf2)*0.5;\
        dvj = 0.5*(xvf+xvf2)/xvf2;\
        cjf = cj0*exp(-z*ln(1-v_j/vd))*dvj+aj*cj0*(1-dvj);\
    end else begin\
        cjf = 0.0;\
    end

// DEPLETION CHARGE CALCULATION
// Hyperbolic smoothing used; no punch-through
`define QJMODF(vj,cj0,vd,z,aj,qjf)\
    if (cj0 > 0.0) begin\
        vf = vd*(1.0-exp(-ln(aj)/z));\
        xvf = (vf-vj)/VT;\
        xvf2 = sqrt(xvf*xvf+1.921812);\
        v_j = vf-VT*(xvf+xvf2)*0.5;\
        x = 1.0-z;\
        y = 1.0-exp(x*ln(1.0-v_j/vd));\
        qjf = cj0*vd*y/x+aj*cj0*(vj-v_j);\
    end else begin\
        qjf = 0.00;\
    end

// Depletion Charge : with punch through
`define QJMOD(vj,cj0,vd,z,vpt,aj,qjf)\
    if (cj0 > 0.0) begin\
        zr = z/4.0;\
        vp = vpt-vd;\
        vf = vd*(1.0-exp(-ln(aj)/z));\
        cmax = aj*cj0;\
        cr = cj0*exp((z-zr)*ln(vd/vpt));\
        a = VT;\
        ve = (vf-vj)/a;\
        if (ve <= `EXPLIM) begin\
            ex1 = exp(ve);\
            ee1 = 1.0+ex1;\
            vj1 = vf-a*ln(ee1);\
        end else begin\
            vj1 = vj;\
        end\
        a = 0.1*vp+4.0*VT;\
        vr = (vp+vj1)/a;\
        if (vr <= `EXPLIM) begin\
            ex1 = exp(vr);\
            ee1 = 1.0+ex1;\
            vj2 = -vp+a*ln(ee1);\
        end else begin\
            vj2 = vj1;\
        end\
        vj4 = vj-vj1;\
        ez = 1.0-z;\
        ezr = 1.0-zr;\
        vdj1 = ln(1.0-vj1/vd);\
        vdj2 = ln(1.0-vj2/vd);\
        qj1 = cj0*(1.0-exp(vdj2*ez))/ez;\
        qj2 = cr*(1.0-exp(vdj1*ezr))/ezr;\
        qj3 = cr*(1.0-exp(vdj2*ezr))/ezr;\
        qjf = (qj1+qj2-qj3)*vd+cmax*vj4;\
    end else begin\
        qjf = 0.0;\
    end

// DEPLETION CHARGE CALCULATION SELECTOR:
// Dependent on junction punch-through voltage
// Important for collector related junctions
`define HICJQ(vj,cj0,vd,z,vpt,qjf)\
    if (vpt < `VPT_thresh) begin\
        `QJMOD(vj,cj0,vd,z,vpt,2.4,qjf)\
    end else begin\
        `QJMODF(vj,cj0,vd,z,2.4,qjf)\
    end

//Temperature dependence of depletion capacitance parameters
`define TMPHICJ(cj0,vd,z,vg,cj0_t,vd_t)\
    arg = 0.5*vd/vt0;\
    vdj0 = 2*vt0*ln(exp(arg)-exp(-arg));\
    vdjt = vdj0*qtt0+vg*(1-qtt0)-mg*VT*ln_qtt0;\
    vd_t = vdjt+2*VT*ln(0.5*(1+sqrt(1+4*exp(-vdjt/VT))));\
    cj0_t = cj0*exp(z*ln(vd/vd_t));

//Limiting exponential
`define LIN_EXP(le, arg)\
    if (arg > 80) begin\
        le = (1 + ((arg) - 80));\
        arg = 80;\
    end else begin\
        le=1;\
    end\
    le = le*limexp(arg);

// conductance not calculated
// INPUT:
//  IS, IST : saturation currents (model parameter related)
//  UM1     : ideality factor
//  U       : branch voltage
// IMPLICIT INPUT:
//  VT      : thermal voltage
// OUTPUT:
//  Iz      : diode current
`define HICDIO(IS,IST,UM1,U,Iz)\
    DIOY = U/(UM1*VT);\
    if (IS > 0.0) begin\
        if (DIOY > 80) begin\
            le = (1 + ((DIOY) - 80));\
            DIOY = 80;\
        end else begin\
            le = 1;\
        end\
        le = le*limexp(DIOY);\
        Iz = IST*(le-1.0);\
        if (DIOY <= -14.0) begin\
            Iz = -IST;\
        end\
    end else begin\
        Iz = 0.0;\
    end

`define qpt_mod(qpt_mod,qlow)\
    o3 = 1.0/3;\
    p2_a = -2*qj_2;\
    if (iqf == `INF && iqfh == `INF) begin\
        p2_b = 0;\
    end else begin\
        p2_b = -(qlow);\
    end\
    p2_c = -itfi*itfi/ick*tfh_t/iqfh_t;\
    tmp = p2_a*p2_a;\
    p2_p = p2_b-tmp*o3;\
    p2_q = 2*p2_a*tmp/27-p2_a*p2_b*o3+p2_c;\
    p2_D = p2_q*p2_q*0.25+p2_p*p2_p*p2_p/27;\
    if (abs(p2_D) < 1e-10) begin\
        q_p3 = 3*p2_q/p2_p-p2_a*o3;\
    end else if (p2_D > 0) begin\
        tmp2 = -p2_q*0.5;\
        tmp3 = sqrt(p2_D);\
        tmp = tmp2+tmp3;\
        if (tmp > 0) begin\
            p2_u = exp(o3*ln(tmp));\
        end else begin\
            p2_u = -exp(o3*ln(-tmp));\
        end\
        tmp = tmp2-tmp3;\
        if (tmp > 0) begin\
            p2_v = exp(o3*ln(tmp));\
        end else begin\
            p2_v = -exp(o3*ln(-tmp));\
        end\
        q_p3 = (p2_u+p2_v)-p2_a*o3;\
    end else begin\
        tmp = -p2_q*0.5*sqrt(-27.0/(p2_p*p2_p*p2_p));\
        tmp2 = tmp*tmp;\
        if (tmp >= 0) begin\
            tmp = `M_PI/2-atan(sqrt(tmp2/(1-tmp2)));\
        end else begin\
            tmp = `M_PI/2+atan(sqrt(tmp2/(1-tmp2)));\
        end\
        tmp = sqrt(-4*p2_p*o3)*cos(o3*tmp)-p2_a*o3;\
        q_p3 = tmp;\
    end\
    qpt_mod = q_p3;

module hic0_full (c,b,e,s,tnode);


//Node definitions

    inout       c,b,e,s,tnode;
    electrical  c     `P(info="external collector node");
    electrical  b     `P(info="external base node");
    electrical  e     `P(info="external emitter node");
    electrical  s     `P(info="external substrate node");
    electrical  ci    `P(info="internal collector node");
    electrical  bi    `P(info="internal base node");
    electrical  ei    `P(info="internal emitter node");
    electrical  tnode `P(info="local temperature rise node");

    electrical      xf1,xf2;
    electrical      xf;  //RC nw


    //Branch definitions
    branch      (ci,c)      br_cic_i;
    branch      (ci,c)      br_cic_v;
    branch      (ei,e)      br_eie_i;
    branch      (ei,e)      br_eie_v;
    branch      (bi,ei)     br_biei;
    branch      (bi,ci)     br_bici;
    branch      (ci,ei)     br_ciei;
    branch      (b,bi)      br_bbi_i;
    branch      (b,bi)      br_bbi_v;
    branch      (b,e)       br_be;
    branch      (b,ci)      br_bci;
    branch      (b,s)       br_bs;
    branch      (s,ci)      br_sci;
    branch      (tnode )    br_sht;

    //Phase network for ITF
    branch      (xf1 )      br_bxf1;
    branch      (xf1 )      br_cxf1;
    branch      (xf2 )      br_bxf2;
    branch      (xf2 )      br_cxf2;

    //Phase network for QF

    branch    (xf )     br_bxf;  //for RC nw
    branch    (xf )     br_cxf;  //for RC nw

//
// Parameter initialization with default values

// Collector current
    parameter real is           = 1.0e-16   from [0:1] `P(spice:name="is" info="(Modified) saturation current" m:factor="yes" unit="A");
    parameter integer it_mod    = 0         from [0:1] `P(spice:name="it_mod" info="Flag for using third order solution for transfer current");
    parameter real mcf          = 1.00      from (0:10] `P(spice:name="mcf" info="Non-ideality coefficient of forward collector current");
    parameter real mcr          = 1.00      from (0:10] `P(spice:name="mcr" info="Non-ideality coefficient of reverse collector current");
    parameter real vef          = `INF      from (0:inf) `P(spice:name="vef" info="forward Early voltage (normalization volt.)"  unit="V" default:value="infinity");
    parameter real ver          = `INF      from (0:inf) `P(spice:name="ver" info="reverse Early voltage (normalization volt.)"  unit="V" default:value="infinity");
    parameter real aver         = 0.0       from [0:100] `P(spice:name="aver" info="bias dependence for reverse Early voltage");
    parameter real iqf          = `INF      from (0:inf) `P(spice:name="iqf" info="forward d.c. high-injection toll-off current" unit="A" m:factor="yes" default:value="infinity");

    parameter real fiqf         = 0.0       from  [0:1] `P(spice:name="fiqf" info="flag for turning on base related critical current" default:value="zero");

    parameter real iqr          = `INF      from (0:inf) `P(spice:name="iqr" info="inverse d.c. high-injection roll-off current" unit="A" m:factor="yes" default:value="infinity");
    parameter real iqfh         = `INF      from (0:inf) `P(spice:name="iqfh" info="high-injection correction current" unit="A" m:factor="yes");
    parameter real tfh          = 0.0       from [0:inf) `P(spice:name="tfh" info="high-injection correction factor" test:value="2e-9" m:factor="yes");
    parameter real ahq          = 0.0       from [-0.9:inf) `P(spice:name="ahq" info="Smoothing factor for the d.c. injection width");

    // Base current
    parameter real ibes         = 1e-18     from [0:1] `P(spice:name="ibes" info="BE saturation current" unit="A" m:factor="yes");
    parameter real mbe          = 1.0       from (0:10] `P(spice:name="mbe" info="BE non-ideality factor");
    parameter real ires         = 0.0       from [0:1] `P(spice:name="ires" info="BE recombination saturation current" test:value="1e-16" unit="A" m:factor="yes");
    parameter real mre          = 2.0       from (0:10] `P(spice:name="mre" info="BE recombination non-ideality factor");
    parameter real ibcs         = 0.0       from [0:1] `P(spice:name="ibcs" info="BC saturation current" test:value="1e-16" unit="A" m:factor="yes");
    parameter real mbc          = 1.0       from (0:10] `P(spice:name="mbc" info="BC non-ideality factor");

    // BE depletion cap
    parameter real cje0         = 1.0e-20   from (0:inf) `P(spice:name="cje0" info="Zero-bias BE depletion capacitance" unit="F" test:value="2e-14" m:factor="yes");
    parameter real vde          = 0.9       from (0:10] `P(spice:name="vde" info="BE built-in voltage" unit="V");
    parameter real ze           = 0.5       from (0:1) `P(spice:name="ze" info="BE exponent factor");
    parameter real aje          = 2.5       from [1:inf) `P(spice:name="aje" info="Ratio of maximum to zero-bias value");
    parameter real vdedc        = 0.9       from (0:10] `P(spice:name="vdedc" info="BE charge built-in voltage for d.c. transfer current" unit="V");
    parameter real zedc         = 0.5       from (0:2) `P(spice:name="zedc" info="charge BE exponent factor for d.c. transfer current");
    parameter real ajedc        = 2.5       from [1:inf) `P(spice:name="ajedc" info="BE capacitance ratio Ratio maximum to zero-bias value for d.c. transfer current");

    // Transit time
    parameter real t0           = 0.0       from [0:inf) `P(spice:name="t0" info="low current transit time at Vbici=0" test:value="5e-12" unit="s");
    parameter real dt0h         = 0.0;      // from [0:inf) `P(spice:name="dt0h" info="Base width modulation contribution" test:value="2e-12" unit="s");
    parameter real tbvl         = 0.0       from [0:inf) `P(spice:name="tbvl" info="SCR width modulation contribution" test:value="4e-12" unit="s");
    parameter real tef0         = 0.0       from [0:inf) `P(spice:name="tef0" info="Storage time in neutral emitter" test:value="1e-12" unit="s");
    parameter real gte          = 1.0       from (0:20] `P(spice:name="gte" info="Exponent factor for emmiter transit time");
    parameter real thcs         = 0.0       from [0:inf) `P(spice:name="thcs" info="Saturation time at high current densities" test:value="3e-11" unit="s");
    parameter real ahc          = 0.1       from (0:10] `P(spice:name="ahc" info="Smoothing facor for current dependence");
    parameter real tr           = 0.0       from [0:inf) `P(spice:name="tr" info="Storage time at inverse operation" unit="s");

    // Critical current
    parameter real rci0         = 150       from (0:inf) `P(spice:name="rci0" info="Low-field collector resistance under emitter" test:value="50" unit="Ohm" m:inverse_factor="yes");
    parameter real vlim         = 0.5       from (0:10] `P(spice:name="vlim" info="Voltage dividing ohmic and satur.region" unit="V");
    parameter real vpt          = 100       from (0:100] `P(spice:name="vpt" info="Punch-through voltage" test:value="10" unit="V" default="infinity");
    parameter real vces         = 0.1       from [0:1] `P(spice:name="vces" info="Saturation voltage" unit="V");

    // BC depletion cap intern
    parameter real cjci0        = 1.0e-20   from (0:inf) `P(spice:name="cjci0" info="Total zero-bias BC depletion capacitance" test:value="1e-15" unit="F" m:factor="yes");
    parameter real vdci         = 0.7       from (0:10] `P(spice:name="vdci" info="BC built-in voltage" test:value="0.7" unit="V");
    parameter real zci          = 0.333     from (0:1] `P(spice:name="zci" info="BC exponent factor" test:value="0.4");
    parameter real vptci        = 100       from (0:100] `P(spice:name="vptci" info="Punch-through voltage of BC junction" test:value="50" unit="V");

    // BC depletion cap extern
    parameter real cjcx0        = 1.0e-20   from [0:inf) `P(spice:name="cjcx0" info="Zero-bias external BC depletion capacitance" unit="F" test:value="1e-15" m:factor="yes");
    parameter real vdcx         = 0.7       from (0:10] `P(spice:name="vdcx" info="External BC built-in voltage" unit="V");
    parameter real zcx          = 0.333     from (0:1] `P(spice:name="zcx" info="External BC exponent factor");
    parameter real vptcx        = 100       from (0:100] `P(spice:name="vptcx" info="Punch-through voltage" unit="V" test:value="5.0" default="infinity");
    parameter real fbc          = 1.0       from [0:1] `P(spice:name="fbc" info="Split factor = Cjci0/Cjc0" test:value="0.5");

    // Base resistance
    parameter real rbi0         = 0.0       from [0:inf) `P(spice:name="rbi0" info="Internal base resistance at zero-bias" test:value="100" unit="Ohm" m:inverse_factor="yes");
    parameter real vr0e         = 2.5       from (0:inf) `P(spice:name="vr0e" info="forward Early voltage (normalization volt.)" unit="V");
    parameter real vr0c         = `INF      from (0:inf) `P(spice:name="vr0c" info="forward Early voltage (normalization volt.)" unit="V" default="infinity" test:value="25.0");
    parameter real fgeo         = 0.656     from [0:inf) `P(spice:name="fgeo" info="Geometry factor" test:value="0.73");

    // Series resistances
    parameter real rbx          = 0.0       from [0:inf) `P(spice:name="rbx" info="External base series resistance" test:value="8.8" unit="Ohm" m:inverse_factor="yes");
    parameter real rcx          = 0.0       from [0:inf) `P(spice:name="rcx" info="Emitter series resistance" test:value="12.5" unit="Ohm" m:inverse_factor="yes");
    parameter real re           = 0.0       from [0:inf) `P(spice:name="re" info="External collector series resistance" test:value="9.16" unit="Ohm" m:inverse_factor="yes");

    // Substrate transfer current, diode current and cap
    parameter real itss         = 0.0       from [0:1.0] `P(spice:name="itss" info="Substrate transistor transfer saturation current" unit="A" test:value="1e-17" m:factor="yes");
    parameter real msf          = 1.0       from (0:10] `P(spice:name="msf" info="Substrate transistor transfer current non-ideality factor");
    parameter real iscs         = 0.0       from [0:1.0] `P(spice:name="iscs" info="SC saturation current" unit="A" test:value="1e-17" m:factor="yes");
    parameter real msc          = 1.0       from (0:10] `P(spice:name="msc" info="SC non-ideality factor");
    parameter real cjs0         = 1.0e-20   from [0:inf) `P(spice:name="cjs0" info="Zero-bias SC depletion capacitance" unit="F" test:value="1e-15" m:factor="yes");
    parameter real vds          = 0.3       from (0:10] `P(spice:name="vds" info="SC built-in voltage" unit="V");
    parameter real zs           = 0.3       from (0:1] `P(spice:name="zs" info="External SC exponent factor");
    parameter real vpts         = 100       from (0:100] `P(spice:name="vpts" info="SC punch-through voltage" unit="V" test:value="5.0" default="infinity");

    // Parasitic caps
    parameter real cbcpar       = 0.0       from [0:inf) `P(spice:name="cbcpar" info="Collector-base isolation (overlap) capacitance" unit="F" m:factor="yes" test:value="1e-15");
    parameter real cbepar       = 0.0       from [0:inf) `P(spice:name="cbepar" info="Emitter-base oxide capacitance" unit="F" m:factor="yes" test:value="2e-15");

    // BC avalanche current
    parameter real eavl         = 0.0       from [0:inf) `P(spice:name="eavl" info="Exponent factor" test:value="1e-14");
    parameter real kavl         = 0.0       from [0:inf) `P(spice:name="kavl" info="Prefactor" test:value="1.19");

    // Flicker noise
    parameter real kf           = 0.0       from [0:inf) `P(spice:name="kf" info="flicker noise coefficient" unit="M^(1-AF)");
    parameter real af           = 2.0       from (0:10] `P(spice:name="af" info="flicker noise exponent factor");

    //Non-quasi-static Effect
    parameter real alqf     = 0.167         from (0:1] `P(spice:name="alqf" info="Factor for additional delay time of minority charge");
    parameter real alit     = 0.333         from (0:1] `P(spice:name="alit" info="Factor for additional delay time of transfer current");
    parameter integer flnqs = 0             from [0:1] `P(spice:name="flnqs" info="Flag for turning on and off of vertical NQS effect");


    // Temperature dependance
    parameter real vgb          = 1.2       from (0:10] `P(spice:name="vgb" info="Bandgap-voltage" unit="V" test:value="1.17");
    parameter real vge          = 1.17      from (0:10] `P(spice:name="vge" info="Effective emitter bandgap-voltage" unit="V" test:value="1.07");
    parameter real vgc          = 1.17      from (0:10] `P(spice:name="vgc" info="Effective collector bandgap-voltage" unit="V" test:value="1.14");
    parameter real vgs          = 1.17      from (0:10] `P(spice:name="vgs" info="Effective substrate bandgap-voltage" unit="V" test:value="1.17");
    parameter real f1vg         =-1.02377e-4            `P(spice:name="f1vg" info="Coefficient K1 in T-dependent bandgap equation" unit="V/K");
    parameter real f2vg         = 4.3215e-4             `P(spice:name="f2vg" info="Coefficient K2 in T-dependent bandgap equation" unit="V/K");
    parameter real alt0         = 0.0                   `P(spice:name="alt0" info="Frist-order TC of tf0" unit="1/K");
    parameter real kt0          = 0.0                   `P(spice:name="kt0" info="Second-order TC of tf0" unit="1/K^2");
    parameter real zetact       = 3.0                   `P(spice:name="zetact" info="Exponent coefficient in transfer current temperature dependence" test:value="3.5");
    parameter real zetabet      = 3.5                   `P(spice:name="zetabet" info="Exponent coefficient in BE junction current temperature dependence" test:value="4.0");
    parameter real zetaci       = 0.0                   `P(spice:name="zetaci" info="TC of epi-collector diffusivity" test:value="1.6");
    parameter real alvs         = 0.0                   `P(spice:name="alvs" info="Relative TC of satur.drift velocity" unit="1/K" test:value="1e-3");
    parameter real alces        = 0.0                   `P(spice:name="alces" info="Relative TC of vces" unit="1/K" test:value="4e-4");
    parameter real zetarbi      = 0.0                   `P(spice:name="zetarbi" info="TC of internal base resistance" test:value="0.6");
    parameter real zetarbx      = 0.0                   `P(spice:name="zetarbx" info="TC of external base resistance" test:value="0.2");
    parameter real zetarcx      = 0.0                   `P(spice:name="zetarcx" info="TC of external collector resistance" test:value="0.2");
    parameter real zetare       = 0.0                   `P(spice:name="zetare" info="TC of emitter resistances");
    parameter real zetaiqf      = 0.0                   `P(spice:name="zetaiqf" info="TC of iqf");
    parameter real alkav        = 0.0                   `P(spice:name="alkav" info="TC of avalanche prefactor" unit="1/K");
    parameter real aleav        = 0.0                   `P(spice:name="aleav" info="TC of avalanche exponential factor" unit="1/K");

    parameter real zetarth      = 0.0                   `P(spice:name="zetarth" info="Exponent factor for temperature dependent thermal resistance" test:value="0.0");

    parameter integer tef_temp  = 1                     `P(spice:name="tef_temp" info="Flag for turning temperature dependence of tef0 on and off");
    parameter real zetaver      = -1.0                  `P(spice:name="zetaver" info="TC of Reverse Early voltage");
    parameter real zetavgbe     = 1.0                   `P(spice:name="zetavgbe" info="TC of AVER");
    parameter real dvgbe        = 0.0                   `P(spice:name="dvgbe" info="Bandgap difference between base and BE-junction");
    parameter real aliqfh       = 0.0                   `P(spice:name="aliqfh" info="Frist-order TC of iqfh" unit="1/K");
    parameter real kiqfh        = 0.0                   `P(spice:name="kiqfh" info="Second-order TC of iqfh" unit="1/K^2");

    // Self-heating
    parameter integer flsh      = 0        from [0:2]   `P(spice:name="flsh" info="Flag for self-heating calculation" test:value="2");
    parameter real rth          = 0.0      from [0:inf) `P(spice:name="rth" info="Thermal resistance" test:value="200.0" unit="K/W" m:inverse_factor="yes");
    parameter real cth          = 0.0      from [0:inf) `P(spice:name="cth" info="Thermal capacitance" test:value="0.1" unit="Ws/K" m:factor="yes");

    // Transistor type

    parameter integer npn       = 1         from [0:1]  `P(spice:isflag="yes" info="model type flag  for npn" );
    parameter integer pnp       = 1         from [0:1]  `P(info="model type flag  for pnp" );
    parameter integer type=(npn==0 ? (pnp==0 ? 0 : 1) : (pnp==0 ? -1 : 0));

    //Circuit simulator specific parameters
    parameter real tnom         = 27.0           `P(spice:name="tnom" info="Temperature for which parameters are valid" unit="C");
    parameter real dtemp        = 0.0            `P(spice:name="dtemp" type="instance" info="Temperature change for particular transistor" unit="K");

// Declaration of the variables: begin

    integer HICUMtype;

    // QCJMOD
    real zr,vp;
    real cmax,cr,ve;
    real ee1,ez,ezr,vdj1,vdj2,ex1,vr,vj1,vj2,vj4;
    real qj1,qj2,qj3;

    //cjtfun    *** tnom,VT,mg,vt0, removed: BA
    real    vdj0,vdjt;

    // temperature and drift
    real VT,Tamb,Tdev,TnomK,dT,qtt0,ln_qtt0;
    real vde_t,vdci_t,vdcx_t,vds_t,vdedc_t;
    real is_t,ires_t,ibes_t,ibcs_t,iqf_t;
    real itss_t,iscs_t,cje0_t,cjci0_t,cjcx0_t, cje0_dc_t, cje0_dc;
    real cjs0_t,rci0_t,vlim_t;
    real vces_t,thcs_t,tef0_t,rbi0_t;
    real rbx_t,rcx_t,re_t,t0_t,eavl_t,kavl_t;
    real aje_t,ajedc_t;

    // bc charge and cap
    real qjci `P(ask="yes" info="B-C internal junction charge" unit="C");
    real qjcx,qjcii,cjcii,qjcxi,qjci_int; //cjcx
    real cjci0_t_ii,cjcx0_t_ii,cjcx0_t_i,v_j;

    // be junction
    real qjei `P(ask="yes" info="B-E internal junction charge" unit="C");
    real vf,x,y;

    // transfer and internal base current
    real cc,qj_2,qj;
    real tf0,ickf,ickr,itfi,itri,qm, qml, qmh;
    real qpt,itf,itr, qpt_l, qpt_h, denom_iqf;
    real b_q;

    real it `P(ask="yes" info="Transfer Current" unit="A");
    real it_wop;
    real ibe,ire,ibi;


    // be diffusion charge
    real qf,qf0,dqfh,dqef;
    real dtef,dtfh,tf,ick;
    real vc,vceff,s3,w,wdc,a,tww, aa, a1, a2;

    // bc diffusion charge
    real qr;

    // avalanche current source
    real v_bord,a_iavl,lncc;

    // base resistance
    real rb,eta,rbi,qje,Qz_nom,fQz;

    // substrate transistor, diode and cap
    real qjs,HSa,HSb,HSI_Tsu,HSUM;

    // self heating
    real pterm;
    real rth_t;

    // new for temperature dependence
    real mg,zetabci,zetasct,zetatef,avs;
    real vgbe,vgbc,vgsc,dvg;
    real xvf,xvf2,dvj,uvc,vt0;

    // noise
    real flicker_Pwr,fourkt,twoq;

    // LIN_EXP
    real    le,arg,le1,arg1,le2,arg2;

    //HICDIO
    real    DIOY;

    // branch voltages
    real    Vbci,Vbici,Vbiei,Vciei,Vsci,Veie,Vbbi,Vcic,Vbe,Vrth;

    //Output to be seen
    real    ijbc `P(ask="yes" info="Base-collector diode current" unit="A");
    real    iavl `P(ask="yes" info="Avalanche current" unit="A");
    real    ijsc `P(ask="yes" info="Substrate-collector diode current" unit="A");
    real    Ibici `P(ask="yes" info="Base-collector diode current minus the avalanche current" unit="A");
    real    ijbe `P(ask="yes" info="Base-emitter diode current" unit="A");

    real    Qbci,Qbe,Qbici,Qbiei;
    real    aver_t,vjh,vj_z,h_vbe,ver_t,iqfh_t,tfh_t,ahq_t;
    real    p2_a,p2_b,p2_c,p2_p,p2_q,p2_D,p2_u,p2_v,q_p3;
    real    tmp,tmp2,tmp3,o3,diff_q;

    //NQS
    real Ixf1,Ixf2,Qxf1,Qxf2,Vxf1,Vxf2,Itxf,Qdeix;
    real Vxf, Ixf, Qxf,fact;

    real gmin;

// For .OP pt calculation
`ifdef CALC_OP
        (* desc="Base terminal current", units="A" multiplicity="multiply" *) real IB;
        (* desc="Collector terminal current", units="A" multiplicity="multiply" *) real IC;
        (* desc="Substrate current", units="A" multiplicity="multiply" *) real ISUB;
        (* desc="Avalanche current", units="A" multiplicity="multiply" *) real IAVL;
        (* desc="External BE voltage", units="V" multiplicity="divide" *) real VBE;
        (* desc="External BC voltage", units="V" multiplicity="divide" *) real VBC;
        (* desc="External CE voltage", units="V" multiplicity="divide" *) real VCE;
        (* desc="External SC voltage", units="V" multiplicity="divide" *) real VSC;
        (* desc="Common emitter forward current gain" multiplicity="none" *) real BETADC;
        (* desc="Internal transconductance", units="S" multiplicity="multiply" *) real GMi;
        (* desc="Internal input resistance", units="Ohm" multiplicity="divide" *) real RPIi;
        (* desc="Internal feedback resistance", units="Ohm" multiplicity="divide" *) real RMUi;
        (* desc="Internal Output resistance", units="Ohm" multiplicity="divide" *) real ROi;
        (* desc="Total BE capacitance", units="F" multiplicity="multiply" *) real CPIi;
        (* desc="Total internal BC capacitance", units="F" multiplicity="multiply" *) real CMUi;
        (* desc="Total external BC capacitance", units="F" multiplicity="multiply" *) real CBCX;
        (* desc="CS junction capacitance", units="F" multiplicity="multiply" *) real CCS;
        (* desc="Internal base resistance", units="Ohm" multiplicity="divide" *) real RBi;
        (* desc="Total base resistance", units="Ohm" multiplicity="divide" *) real RB;
        (* desc="External (saturated) collector series resistance", units="Ohm" multiplicity="divide" *) real RCXop;
        (* desc="Emitter series resistance", units="Ohm" multiplicity="divide" *) real REop;
        (* desc="Small signal current gain" multiplicity="none" *) real BETAAC;
        (* desc="Total forward transit time", units="s" multiplicity="none" *) real TF;
        (* desc="Transit frequency", units="Hz" multiplicity="none" *) real FT;
`endif


//Declaration of the variables: end


//
//======================== calculation of the transistor ===================
//

    analog begin

        gmin = `GMIN;

        `INITIAL_MODEL     // Model Initialization
        begin

            if (`PGIVEN(npn)) begin
                HICUMtype =  `NPN;
            end else if (`PGIVEN(pnp)) begin
                HICUMtype =  `PNP;
            end else begin
                HICUMtype = (`PGIVEN(type)) ? type : `NPN;
            end

        end

// assign voltages with regard to transistor type

        Vbci  = HICUMtype*V(br_bci);
        Vbici = HICUMtype*V(br_bici);
        Vbiei = HICUMtype*V(br_biei);
        Vciei = HICUMtype*V(br_ciei);
        Vsci  = HICUMtype*V(br_sci);
        Veie  = V(br_eie_v);
        Vcic  = V(br_cic_v);
        Vbbi  = V(br_bbi_v);
        Vbe   = HICUMtype*V(br_be);
        Vrth  = V(br_sht);

//
// temperature and resulting parameter drift
//

        TnomK   = tnom+273.15;
        Tamb    = $temperature;
        Tdev    = Tamb+dtemp+Vrth;

// Limit temperature to avoid FPE's in equations
        if (Tdev < `TMIN + 273.15) begin
            Tdev = `TMIN + 273.15;
        end else begin
            if (Tdev > `TMAX + 273.15) begin
                Tdev = `TMAX + 273.15;
            end
        end

        vt0     = `P_K*TnomK /`P_Q;
        VT      = `P_K*Tdev /`P_Q;
        dT      = Tdev-TnomK;
        qtt0    = Tdev/TnomK;
        ln_qtt0 = ln(qtt0);
        avs     = alvs*TnomK;
        vgbe    = (vgb+vge)/2;
        vgbc    = (vgb+vgc)/2;
        vgsc    = (vgs+vgc)/2;
        mg      = 3-`P_Q*f1vg/`P_K;
        zetabci = mg+1-zetaci;
        zetasct = mg-1.5;  //+1-m_upS with m_upS=2.5
        is_t    = is*exp(zetact*ln_qtt0+vgb/VT*(qtt0-1));
        ibes_t  = ibes*exp(zetabet*ln_qtt0+vge/VT*(qtt0-1));
        ires_t  = ires*exp(0.5*mg*ln_qtt0+0.5*vgbe/VT*(qtt0-1));
        ibcs_t  = ibcs*exp(zetabci*ln_qtt0+vgc/VT*(qtt0-1));
        itss_t  = itss*exp(zetasct*ln_qtt0+vgc/VT*(qtt0-1));
        iscs_t  = iscs*exp(zetasct*ln_qtt0+vgs/VT*(qtt0-1));

        `TMPHICJ(cje0,vde,ze,vgbe,cje0_t,vde_t)

   //  `TMPHICJ(cje0,vde,zedc,vgbe,cje0_t,vdedc_t)

        cje0_dc = cje0;

        `TMPHICJ(cje0_dc,vdedc,zedc,vgbe,cje0_dc_t,vdedc_t)   //introducing DC capacitance


        aje_t   = aje*vde_t/vde;

        ajedc_t = ajedc*vdedc_t/vdedc;

        `TMPHICJ(cjci0,vdci,zci,vgbc,cjci0_t,vdci_t)
        `TMPHICJ(cjcx0,vdcx,zcx,vgbc,cjcx0_t,vdcx_t)
        `TMPHICJ(cjs0,vds,zs,vgsc,cjs0_t,vds_t)
        iqf_t   = iqf*exp(zetaiqf*ln_qtt0-dvgbe/VT*(qtt0-1));
        rci0_t  = rci0*exp(zetaci*ln_qtt0);
        vlim_t  = vlim*exp((zetaci-avs)*ln_qtt0);
        vces_t  = vces*(1+alces*dT);
        t0_t    = t0*(1+alt0*dT+kt0*dT*dT);
        thcs_t  = thcs*exp((zetaci-1)*ln_qtt0);
        zetatef = zetabet-zetact-0.5;
        dvg     = vgb-vge;
        if (tef_temp == 1) begin
            tef0_t  = tef0*exp(zetatef*ln_qtt0-dvg/VT*(qtt0-1));
        end else begin
            tef0_t  = tef0;
        end
        rbx_t   = rbx*exp(zetarbx*ln_qtt0);
        rcx_t   = rcx*exp(zetarcx*ln_qtt0);
        rbi0_t  = rbi0*exp(zetarbi*ln_qtt0);
        re_t    = re*exp(zetare*ln_qtt0);
        eavl_t  = eavl*exp(aleav*dT);
        kavl_t  = kavl*exp(alkav*dT);


        //Temperature dependence of Thermal resistance
        if (zetarth!=0) begin
            rth_t  = rth*exp(zetarth*ln(Tdev/TnomK));
        end else begin
            rth_t=rth;
        end

        aver_t = aver*exp(zetaver*ln_qtt0);
        ver_t = ver/exp(dvgbe/VT*(exp(zetavgbe*ln_qtt0)-1));
        iqfh_t = iqfh*(1+aliqfh*dT+kiqfh*dT*dT);
        tfh_t = tfh*(1+aliqfh*dT+kiqfh*dT*dT)*exp((vgb-vge)/VT*(qtt0-1));
        ahq_t = ahq;

//
// Calculation of intrinsic transistor elements
//

// BC charge and cap (internal and external)

// The cjcx0 value is used to switch between one (cjcx0=0) and two bc parameter sets
// 1. For one parameter set only the internal bc set is partitioned by fbc
// 2. For two independent sets only the external set is partitioned by fbc

        if (cjcx0_t > 0.0) begin
            cjci0_t_ii  = cjci0_t;                 // zero bias internal portion
            cjcx0_t_ii  = cjcx0_t*fbc;
            `HICJQ(Vbici,cjcx0_t_ii,vdcx_t,zcx,vptcx,qjcxi)
            cjcx0_t_i   = cjcx0_t*(1-fbc);         // zero bias external portion
            `HICJQ(Vbci,cjcx0_t_i,vdcx_t,zcx,vptcx,qjcx)
        end else begin
            cjci0_t_ii  = cjci0_t*fbc;             // zero bias internal portion
            qjcxi       = 0;
            cjcx0_t_i   = cjci0_t*(1-fbc);         // zero bias external portion
            `HICJQ(Vbci,cjcx0_t_i,vdci_t,zci,vptci,qjcx)
        end
        `HICJQ(Vbici,cjci0_t_ii,vdci_t,zci,vptci,qjci)
        qjci_int = qjci; // int BC charge without normalization
        qjcii   = qjci+qjcxi;

//Internal bc cap without punch through for cc

        //`HICJQ(Vbici,cjci0_t_ii,vdci_t,zci,100,qjciii)
        `QCMODF(Vbici,cjci0_t_ii,vdci_t,zci,2.4,cjcii)
        //cjcii = ddx(qjciii,V(bi));

//Internal be cap and charge

//  `QJMODF(Vbiei,cje0_dc_t,vde_t,ze,aje_t,qjei)
//  cjei    = ddx(qjei,V(bi));

// Critical current: ick
        vc  = Vciei-vces_t;
        uvc = vc/VT-1;
        vceff   = VT*(1+0.5*(uvc+sqrt(uvc*uvc+1.921812)));
        x   = (vceff-vlim_t)/vpt;
        ick = vceff*(1+0.5*(x+sqrt(x*x+1e-3)))/rci0_t/sqrt(1+vceff*vceff/vlim_t/vlim_t);


// Normalized BC cap and charge

        if (cjcii > 0.0 && cjci0_t_ii > 0.0) begin

            cc  = cjci0_t_ii/cjcii;
            qjci   = qjci/cjci0_t_ii;

        end else begin

            cc      = 1.0;
            qjci = 0;

        end

    //cc    = cjci0_t_ii/cjcii;
    //qjci  = qjci/cjci0_t_ii;

        `QJMODF(Vbiei,cje0_dc_t,vdedc_t,zedc,ajedc_t,qjei)

        if (aver == 0.0) begin
            h_vbe = 1;
        end else begin
            vjh = (vdedc_t-Vbiei)/(2.0*VT);
            vjh = vdedc_t-2.0*VT*(vjh+sqrt(vjh*vjh+1.921812))*0.5;
            vjh = (vjh-VT)/VT;
            vjh = VT*(1.0+(vjh+sqrt(vjh*vjh+1.921812))*0.5);
            vj_z = (1.0-exp(zedc*ln(1.0-vjh/vdedc_t)))*aver_t;
            h_vbe = (exp(vj_z)-1.0)/vj_z;
        end

        qje   = h_vbe*qjei/cje0_dc_t;
        qj  = (1+qjci/vef+qje/ver_t);

        b_q  = 20*qj-1;
        qj_2=0.025*(1+(b_q +sqrt(b_q*b_q+1.921812))/2);

// Minority charge transit time
        tf0 = t0_t+dt0h*(cc-1)+tbvl*(1/cc-1);

    //Determination of base realted critical current

        if (fiqf==1)begin
            denom_iqf = fiqf*((tf0/t0_t)-1);
            ickf   =  iqf_t/(1+denom_iqf);
        end else begin
            ickf    = iqf_t;
        end

        ickr    = iqr;

// Ideal transfer currents
        arg1 = Vbiei/(mcf*VT);
        `LIN_EXP(le1,arg1)
        itfi=is_t*le1;

        arg2 = Vbici/(mcr*VT);
        `LIN_EXP(le2,arg2)
        itri=is_t*le2;
// Normalized minority charge at low currents (w=0) and high currents (w=1)

        if (tfh!=0)begin
            qml = itfi/ickf+itri/ickr+exp((0.6666)*ln(itfi*(itfi/ick)*((tfh_t)/iqfh_t)));
            qmh = itfi/ickf+itri/ickr+itfi/iqfh_t+exp((0.6666)*ln(itfi*(itfi/ick)*((tfh_t)/iqfh_t)));
        end else begin
            qml = itfi/ickf+itri/ickr;
            qmh = itfi/ickf+itri/ickr+itfi/iqfh_t;
        end
        qpt_l= qj_2+sqrt((qj_2)*(qj_2)+qml);
        qpt_h= qj_2+sqrt((qj_2)*(qj_2)+qmh);

// Calculation of the injection width
        diff_q = qmh-qml;
        if (abs(diff_q)>1e-8) begin
            a1= 1-ick/(1+ahq_t)/itfi*qpt_l;
            a2= 1+ick/(1+ahq_t)/itfi*(qpt_h-qpt_l);
            aa= a1/a2;

            wdc= (sqrt(aa*aa+0.01)+aa)/(1+sqrt(1+0.01));
        end else begin
            wdc = 0;
        end

// Normalized minority charge

        if (it_mod == 0) begin
            if (tfh!=0) begin
                qm  = itfi/ickf+itri/ickr+itfi/iqfh_t*wdc*wdc+exp((0.6666)*ln(itfi*(itfi/ick)*((tfh_t)/iqfh_t)));
            end else begin
                qm  = itfi/ickf+itri/ickr+itfi/iqfh_t*wdc*wdc;
            end
            // Normalized total hole charge
            qpt = qj_2+sqrt((qj_2)*(qj_2)+qm);
        end else begin
            `qpt_mod(qpt,itfi/ickf+itri/ickr+itfi/iqfh_t*wdc*wdc)
        end
        if (qpt<=1e-20) begin
            qpt=1e-20;
        end

        itf = itfi/qpt;
        itr = itri/qpt;

        // Transfer current

        if (itf<=1e-20) begin
            itf = 1e-20;
        end
        it  = itf-itr;

// BE diffusion charge

// Calculation of low-current portion
        qf0 = tf0*itf;

// Current dependent component
        a   = 1-ick/itf;
        s3  = sqrt(a*a+ahc);
        w   = (a+s3)/(1+sqrt(1+ahc));
        tww = thcs_t*w*w;
        dqfh    = tww*itf;
        dtfh    = tww*(1+2*ick/itf/s3);

// Emitter component
        dtef    = tef0_t*exp(gte*ln(itf/ick));
        dqef    = dtef*itf/(gte+1.0);

// Total minority charge and transit time
        qf  = qf0+dqef+dqfh;
        tf  = tf0+dtfh+dtef;

// BC diffusion charge
        qr  = tr*itr;

// Internal base current

// BE diode
        `HICDIO(ibes,ibes_t,mbe,Vbiei,ibe)
        `HICDIO(ires,ires_t,mre,Vbiei,ire)
        ijbe    = ibe+ire;

// BC diode
        `HICDIO(ibcs,ibcs_t,mbc,Vbici,ijbc)

// Total base current
        ibi = ijbe+ijbc;

// Avalanche current

        if (Vbici < 0) begin : HICAVL
            v_bord  = eavl_t*vdci_t;
            if (vdci_t-Vbici>v_bord) begin
                a_iavl  = kavl_t/vdci_t*exp(-cc);
                iavl    = itf*a_iavl*(v_bord+(1+cc)*(vdci_t-Vbici-v_bord));
            end else begin
                lncc    = ln(1/cc);
                iavl    = kavl_t*itf*exp(-1/zci*lncc-eavl_t*exp((1/zci-1)*lncc));
            end
        end else begin
            iavl    = 0;
        end

//
// Additional elements for external transistor
//
        `QJMODF(Vbiei,cje0_t,vde_t,ze,aje_t,qjei)  // Computation of AC charge for base resistance calculation
        qje   = qjei/cje0_t;

// Base resistance
        if (rbi0_t > 0.0) begin : HICRBI
            // Conductivity modulation with hyperbolic smoothing

            Qz_nom  = 1+qje/vr0e+qjci/vr0c+itf/ickf+itr/ickr;
            fQz = 0.5*(Qz_nom+sqrt(Qz_nom*Qz_nom+0.01));
            rbi = rbi0_t/fQz;

                //rbi= rbi0_t*(1+0.2)/(0.2+qpt);
            // Emitter current crowding
            if (ibi > 0.0) begin
                eta = fgeo*rbi*ibi/VT;
                if (eta < 1e-6) begin
                    rbi = rbi*(1-0.5*eta);
                end else begin
                    rbi = rbi*ln(eta+1)/eta;
                end
            end
        end else begin
            rbi = 0.0;
        end
    // Total base resistance
      //rbi= rbi0_t;
        rb  = rbi+rbx_t;

// Parasitic substrate transistor transfer current
        if (itss > 0.0) begin : Sub_Transfer
            HSUM    = msf*VT;
            HSa = limexp(Vbci/HSUM);
            HSb = limexp(Vsci/HSUM);
            HSI_Tsu = itss_t*(HSa-HSb);
        end else begin
            HSI_Tsu = 0.0;
        end

// Substrate diode and cap and charge

        `HICDIO(iscs,iscs_t,msc,Vsci,ijsc)

        `HICJQ(Vsci,cjs0_t,vds_t,zs,vpts,qjs)

// Self heating

        pterm = 0;

        if (flsh == 1 && rth >= `MIN_R) begin
            pterm = it*Vciei+iavl*(vdci_t-Vbici);
        end else if (flsh == 2 && rth >= `MIN_R) begin
            pterm =  Vciei*it + (vdci_t-Vbici)*iavl + ijbe*Vbiei + ijbc*Vbici + ijsc*Vsci;
            if (rb >= `MIN_R) begin
                pterm = pterm + Vbbi*Vbbi/rb;
            end
            if (re >= `MIN_R) begin
                pterm = pterm + Veie*Veie/re_t;
            end
            if (rcx >= `MIN_R) begin
                pterm = pterm + Vcic*Vcic/rcx_t;
            end
        end

        Itxf  = itf;
        Qdeix = qf;
    // Excess Phase calculation

        if (flnqs != 0 && tf != 0) begin
            Vxf1    = V(br_bxf1);
            Vxf2    = V(br_bxf2);

            Ixf1    =  (Vxf2-itf)/tf*t0;
            Ixf2    =  (Vxf2-Vxf1)/tf*t0;
            Qxf1    =  alit*Vxf1*t0;
            Qxf2    =  alit*Vxf2/3*t0;
            Itxf    =  Vxf2;

            Vxf = V(br_bxf);
            fact  = t0/tf;
            Ixf   = (Vxf - qf)*fact;
            Qxf   = alqf*Vxf*t0;
            Qdeix = Vxf;
        end else begin
            Ixf1    =  V(br_bxf1);
            Ixf2    =  V(br_bxf2);
            Qxf1    =  0;
            Qxf2    =  0;

            Ixf   = V(br_bxf);
            Qxf   = 0;
        end


//
// Compute branch sources
//

        Ibici = ijbc - iavl;

        Qbci  = cbcpar*Vbci;
        Qbe   = cbepar*Vbe;
        Qbici = qjcii+qr;
        //Qbiei = qjei+qf;
        Qbiei = qjei+Qdeix;

        ijsc  = HICUMtype*ijsc;
        qjs   = HICUMtype*qjs;
        qjcx  = HICUMtype*qjcx;
        Qbci  = HICUMtype*Qbci;
        Qbe   = HICUMtype*Qbe;

        Ibici = HICUMtype*Ibici;
        Qbici = HICUMtype*Qbici;
        ijbe  = HICUMtype*ijbe;
        Qbiei = HICUMtype*Qbiei;
        it_wop    = HICUMtype*it; // 'it' without excess phase
        it    = HICUMtype*(Itxf-itr);   // 'it' with excess phase
        iavl =  HICUMtype*iavl;
//
// Define branch sources
//
        I(br_biei)  <+ gmin*V(br_biei);
        I(br_bici)  <+ gmin*V(br_bici);
        I(br_ciei)  <+ gmin*V(br_ciei);

        I(br_bs)        <+ HICUMtype*HSI_Tsu;
        I(br_sci)       <+ ijsc + gmin*V(br_sci);//`P(spectre:gmin="add" spectre:pwl_passive="1e10");
        I(br_sci)       <+ ddt(qjs);
        I(br_bci)       <+ ddt(qjcx);
        I(br_bci)       <+ ddt(Qbci);
        I(br_be)        <+ ddt(Qbe);
        if (re >= `MIN_R) begin
            I(br_eie_i)     <+ Veie/re_t;      //`P(spectre:gmin="add");
        end else begin
            I(br_eie_i)     <+ Veie/`MIN_R;    //`P(spectre:gmin="add");
        end
        if (rcx >= `MIN_R) begin
            I(br_cic_i)     <+ Vcic/rcx_t;     //`P(spectre:gmin="add");
        end else begin
            I(br_cic_i)     <+ Vcic/`MIN_R;    //`P(spectre:gmin="add");
        end
        if (rbi0 >= `MIN_R || rbx >= `MIN_R) begin
            I(br_bbi_i)     <+ Vbbi/rb + gmin*Vbbi; //`P(spectre:gmin="add");
        end else begin
            I(br_bbi_i)     <+ Vbbi/`MIN_R;    //`P(spectre:gmin="add");
        end
        I(br_bici)  <+ Ibici;                  //`P(spectre:gmin="add" spectre:pwl_sat_current="IMAX" spectre:pwl_sat_cond="imax/0.025" spectre:pwl_rev_current="imax" spectre:pwl_rev_cond="IMAX/0.025");
        I(br_bici)  <+ ddt(Qbici);
        I(br_biei)  <+ ijbe;                   //`P(spectre:gmin="add" spectre:pwl_fwd_current="IBEIS*exp(25.0)" spectre:pwl_fwd_node="bi" spectre:pwl_fwd_cond="IBEIS*exp(25.0)/0.025" spectre:pwl_sat_current="IMAX" spectre:pwl_sat_cond="IMAX/0.025" spectre:pwl_passive="1e10");
        I(br_biei)  <+ ddt(Qbiei);
        I(br_ciei)  <+ it;                     //`P(spectre:pwl_fwd_current="IS*exp(25.0)" spectre:pwl_fwd_node="bi" spectre:pwl_fwd_cond="IS*exp(25.0)/0.025"   spectre:pwl_rev_current="IMAX" spectre:pwl_rev_cond="IMAX/0.025" spectre:pwl_passive="1e10");
        //I(br_ciei)    <+ Itxf;

        // Following code is an intermediate solution:
        // ******************************************
        if (flsh == 0 || rth < `MIN_R) begin
            I(br_sht)   <+ Vrth/`MIN_R;
        end else begin
            I(br_sht)   <+ Vrth/rth_t-pterm;     //`P(spectre:gmin="add");
            I(br_sht)   <+ ddt(cth*Vrth);
        end
        // ******************************************
        // For simulators having no problem with V(br_sht) <+ 0.0
        // with external thermal node, follwing code may be used.
        // This external thermal node should remain accessible.
        // ********************************************
        //if (flsh == 0 || rth < `MIN_R) begin
        //  V(br_sht)   <+ 0.0;
        //end else begin
        //  I(br_sht)   <+ Vrth/rth_t-pterm ;  //`P(spectre:gmin="add");
        //  I(br_sht)   <+ ddt(cth*Vrth);
        //
        //end
        // ********************************************

        // NQS effect
        I(br_bxf1) <+ Ixf1;
        I(br_cxf1) <+ ddt(Qxf1);
        I(br_bxf2) <+ Ixf2;
        I(br_cxf2) <+ ddt(Qxf2);

        I(br_bxf) <+  Ixf;
        I(br_cxf) <+  ddt(Qxf);


`ifdef CALC_NOISE
// Noise sources
// Thermal noise
            fourkt  = 4.0 * `P_K * Tdev;
            if (rbx >= `MIN_R || rbi0 >= `MIN_R) begin
                I(br_bbi_i)     <+ white_noise(fourkt/rb, "rb thermal");
            end
            if (rcx >= `MIN_R) begin
                I(br_cic_i)     <+ white_noise(fourkt/rcx_t, "rcx thermal");
            end
            if (re >= `MIN_R) begin
                I(br_eie_i)     <+ white_noise(fourkt/re_t, "re thermal");
            end

// Shot noise
            twoq    = 2.0 * `P_Q;
            I(br_biei)  <+ white_noise(twoq*ijbe, "ijbe shot");
            I(br_ciei)  <+ white_noise(twoq*it, "it shot");

// Flicker noise
            flicker_Pwr = kf*pow(ijbe,af);
            I(br_biei)  <+ flicker_noise(flicker_Pwr,1.0, "ib flicker");
`endif


`ifdef CALC_OP
//            if (analysis("static")) begin : OP_calculation
            begin : OP_calculation
                real oRPIi, oRMUi, oROi, gAVL;
                real Cdei, Cdci, Cjei, Cjci, Cjcx, CBC;
                real R_tot;

//                IB = I(<b>);
//                IC = I(<c>);
//                ISUB = I(<s>);
                IB = ibe;
                IC = it;
                ISUB = ijsc;
                IAVL = iavl;

                VBE = V(b,e);
                VBC = V(b,c);
                VCE = V(c,e);
                VSC = V(s,c);

//                GMi = ddx(it_wop,V(bi));
                GMi = ddx(it,V(bi));

                oRPIi = ddx(ijbe,V(bi));
                RPIi = 1.0/(oRPIi+1e-12);

                oRMUi = -1*ddx(Ibici,V(ci));
                RMUi = 1.0/(oRMUi+1e-12);

                gAVL = HICUMtype*ddx(iavl,V(ci));
                oROi = ddx(it_wop,V(ci));
                ROi = 1.0/(oROi+gAVL+1e-12);

                Cdei = -1*HICUMtype*ddx(qf,V(ei));
                Cjei = ddx(qjei,V(bi));
                CPIi = Cjei + Cdei + cbepar;

                Cdci = -1*HICUMtype*ddx(qr,V(ci));
                Cjci = ddx(qjci_int,V(bi));
                CMUi = Cjci + Cdci;

                Cjcx = -1*ddx(qjcxi,V(ci));
                CBCX = Cjcx + cbcpar ;

                CCS = ddx(qjs,V(s));

                RBi = rbi;
                RB = rb;
                RCXop = rcx_t;
                REop = re_t;

                BETADC = IC/IB;
                BETAAC = GMi*RPIi;

                R_tot = RCXop + REop + ((RB+REop)/BETAAC);

                TF = tf;

                CBC = CMUi+CBCX;
                FT = GMi/(2*`M_PI*(CPIi+CBC+(R_tot*CBC*GMi)));
            end
`endif


    end  // analog
endmodule
