//HICUM Level_2 Version_2.24: A Verilog-A Description

//dw 09/10: Modifications for ngspice and adms:
// backup to ddx for capacitance calculation
// all V(...) <+ 0.0; are replaced by I(...) < V(...)/`MIN_R;
// using GMIN from ngspice: _circuit_gmin
// switch of section for correlated noise (see below)
// removed obsolete variables S_avl, f_p
// don't like internal variable declaration

//**************New Implementations*****************

//******************************************************************************
//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. 
//Although Vertical NQS effects have been taken into account in HICUM from the very 
//beginning (see original FTN code and built-in v2.1 HICUM model inside most of the 
//existing circuit simulators) their implementation has been based on Weil's approach. 
//However, using Verilog, it is presently not possible to implement Weil's approach, 
//since there does not exist access to previous time-steps of the simulatior. 
//The nearly available Verilog-A solution reproduces the results of previous 
//HICUM versions (cf. documentation).
//******************************************************************************

//******************************************************************************
// Implementation of noise correlation
// Please turn-off (by front slash “//”) the noise correlation code section for simulation with Spectre
//******************************************************************************

// ***************Bug fix and optimization*************
// 12/09: Elimination of 'ddx' operator for internal capacitance determination
//        Capacitances will be calculated by Analytical Equations 
// 10/09: Temperature coefficient (ZETAxyz) range modified to [-10:10]
// 09/09: VPT ranges modified from [0:inf] to (0:inf]
// 07/09: DT0H rabge modified from no range to (-inf:inf) 
// 03/09: Simplified input NQS adjacent circuit (RC network)
// 04/08: New range has been defined for FDQR0.
// 11/07: Bugs have been fixed in macro HICFCI and HICQFC
// 10/06: in @(initial_model), external if-block for HICTUN_T removed
// 11/06: within HICQFC, minor changes made for LATB<=0.01;
// also HICFCI and HICFCT are changed accordingly
// to ensure correct derivatives
// Upper limit of FGEO parameter was changed to infinity.
// 12/06: expressions for Cdei and Cdci are corrected not to include
// Ccdei and Cbdci respectively (used in Crbi expression).

// 01/06: FCdf1_dw assigned expression (missing in v2.21)
// FCa and FCa1 are found to have same expression: FCa is omitted in those cases
// FCa1 written instead of FCa in the expression for FCf_ci
// Thermal node "tnode" set as external
// zetasct = mg+1-2.5 changed to zetasct = mg-1.5;
// Code optimization: Temperature dependent parts are moduled in two separate blocks:
// within @(initial_model) when self-heating is OFF
// outside @(initial_model) when self-heating is ON
// 03/06 : Further fix
// vlim_t,ibcis_t,ibcxs_t,itss_t,iscs_t considered in compatibility block
// ddt() operators are separated in contribution expressions.
// FLCOMP parameter is given different values
// 05/06:
// all if-else blocks marked with begin-end
// unused variables deleted
// all series resistors and RTH are allowed to have a minimum value MIN_R
// only tunelling current source contribution within if-then-else
// 06/06: HICRBI deleted and instead the code changed (hyperbolic smoothing in
// conductivity modulation part) and put in relevant portion of the code.
// 07/06: ddx() operator used to find out capacitances from charges:
// QJMODF,QJMOD,HICJQ changed accordingly
// Lateral NQS effect modified with ddx() operator.
// HICFCT included for downward compatibility reason.
// Few macros are taken inside the code: HICICK, HICAVL, HICTUN (more optimized),
// internal base resistance (Qjci included under conductivity modulation, hyperbolic smoothing used)
// Gmin added at (bi,ei) and (bi,ci) branches.
// 08/06: Units added in the parameter descriptions.

// *********************************************************************************
// 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.
// **********************************************************************************


//Default simulator: Spectre

`ifdef insideADMS
 `define MODEL @(initial_model)
 `define NOISE @(noise)
 `define ATTR(txt) (*txt*)
`else
 `define MODEL
 `define NOISE
 `define ATTR(txt)
`endif


`define VPT_thresh      1.0e2
`define Dexp_lim        80.0
`define Cexp_lim        80.0
`define DFa_fj          1.921812
`define RTOLC           1.0e-5
`define l_itmax         100
`define TMAX            326.85
`define TMIN            -100.0
`define LN_EXP_LIMIT    11.0
`define MIN_R           0.001
//`define Gmin            1.0e-12
//`define Gmin      $simparam("gmin",1e-12)  //suggested by L.L

//ADS
//`include "constants.vams"
//`include "disciplines.vams"
//`include "compact.vams"

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

//////////////Explicit Capacitance and Charge Expression///////////////

// DEPLETION CHARGE CALCULATION
// Hyperbolic smoothing used; no punch-through
// INPUT:
//  c_0		: zero-bias capacitance
//  u_d		: built-in voltage
//  z		: exponent coefficient
//  a_j		: control parameter for C peak value at high forward bias
//  U_cap	: voltage across junction
// IMPLICIT INPUT:
//  VT		: thermal voltage
// OUTPUT:
//  Qz		: depletion Charge 
//  C		: depletion capacitance
`define QJMODF(c_0,u_d,z,a_j,U_cap,C,Qz)\
    if(c_0 > 0.0) begin\
        DFV_f	 = u_d*(1.0-exp(-ln(a_j)/z));\
        DFv_e	 = (DFV_f-U_cap)/VT;\
        DFs_q	 = sqrt(DFv_e*DFv_e+`DFa_fj);\
        DFs_q2	 = (DFv_e+DFs_q)*0.5;\
        DFv_j	 = DFV_f-VT*DFs_q2;\
        DFdvj_dv = DFs_q2/DFs_q;\
        DFb	     = ln(1.0-DFv_j/u_d);\
        DFC_j1   = c_0*exp(-z*DFb)*DFdvj_dv;\
        C		 = DFC_j1+a_j*c_0*(1.0-DFdvj_dv);\
        DFQ_j	 = c_0*u_d*(1.0-exp(DFb*(1.0-z)))/(1.0-z);\
        Qz	     = DFQ_j+a_j*c_0*(U_cap-DFv_j);\
    end else begin\
        C        = 0.0;\
        Qz	     = 0.0;\
    end
   
////////////////////////////////////////////////////////////////


//////////////Explicit Capacitance and Charge Expression///////////////


// DEPLETION CHARGE CALCULATION CONSIDERING PUNCH THROUGH
// smoothing of reverse bias region (punch-through) 
// and limiting to a_j=Cj,max/Cj0 for forward bias.
// Important for base-collector and collector-substrate junction
// INPUT:
//  c_0		: zero-bias capacitance
//  u_d		: built-in voltage
//  z 		: exponent coefficient
//  a_j		: control parameter for C peak value at high forward bias
//  v_pt	: punch-through voltage (defined as qNw^2/2e)
//  U_cap	: voltage across junction
// IMPLICIT INPUT:
//  VT		: thermal voltage
// OUTPUT:
//  Qz		: depletion charge
//  C		: depletion capacitance
`define QJMOD(c_0,u_d,z,a_j,v_pt,U_cap,C,Qz)\
    if(c_0 > 0.0) begin\
        Dz_r	= z/4.0;\
        Dv_p	= v_pt-u_d;\
        DV_f	= u_d*(1.0-exp(-ln(a_j)/z));\
        DC_max	= a_j*c_0;\
        DC_c	= c_0*exp((Dz_r-z)*ln(v_pt/u_d));\
        Dv_e	= (DV_f-U_cap)/VT;\
        if(Dv_e < `Cexp_lim) begin\
            De	    = exp(Dv_e);\
            De_1	= De/(1.0+De);\
            Dv_j1	= DV_f-VT*ln(1.0+De);\
        end else begin\
            De_1	= 1.0;\
            Dv_j1	= U_cap;\
        end\
        Da	    = 0.1*Dv_p+4.0*VT;\
        Dv_r	= (Dv_p+Dv_j1)/Da;\
        if(Dv_r < `Cexp_lim) begin\
            De	    = exp(Dv_r);\
            De_2	= De/(1.0+De);\
            Dv_j2	= -Dv_p+Da*ln(1.0+De);\
        end else begin\
            De_2	= 1.0;\
            Dv_j2	= Dv_j1;\
        end\
        Dv_j4	= U_cap-Dv_j1;\
        DCln1	= ln(1.0-Dv_j1/u_d);\
        DCln2	= ln(1.0-Dv_j2/u_d);\
        Dz1	    = 1.0-z;\
        Dzr1	= 1.0-Dz_r;\
        DC_j1	= c_0*exp(DCln2*(-z))*De_1*De_2;\
        DC_j2	= DC_c*exp(DCln1*(-Dz_r))*(1.0-De_2);\
        DC_j3	= DC_max*(1.0-De_1);\
        C		= DC_j1+DC_j2+DC_j3;\
        DQ_j1	= c_0*(1.0-exp(DCln2*Dz1))/Dz1;\
        DQ_j2	= DC_c*(1.0-exp(DCln1*Dzr1))/Dzr1;\
        DQ_j3	= DC_c*(1.0-exp(DCln2*Dzr1))/Dzr1;\
        Qz	    = (DQ_j1+DQ_j2-DQ_j3)*u_d+DC_max*Dv_j4;\
    end else begin\
        C	    = 0.0;\
        Qz	    = 0.0;\
    end


// DEPLETION CHARGE & CAPACITANCE CALCULATION SELECTOR
// Dependent on junction punch-through voltage
// Important for collector related junctions
`define HICJQ(c_0,u_d,z,v_pt,U_cap,C,Qz)\
    if(v_pt < `VPT_thresh) begin\
        `QJMOD(c_0,u_d,z,2.4,v_pt,U_cap,C,Qz)\
    end else begin\
        `QJMODF(c_0,u_d,z,2.4,U_cap,C,Qz)\
    end


// A CALCULATION NEEDED FOR COLLECTOR MINORITY CHARGE FORMULATION
// INPUT:
//  zb,zl       : zeta_b and zeta_l (model parameters, TED 10/96)
//  w           : normalized injection width
// OUTPUT:
// hicfcio      : function of equation (2.1.17-10)
`define HICFCI(zb,zl,w,hicfcio,dhicfcio_dw)\
    z       = zb*w;\
    lnzb    = ln(1+zb*w);\
    if(z > 1.0e-6) begin\
        x               = 1.0+z;\
        a               = x*x;\
        a2              = 0.250*(a*(2.0*lnzb-1.0)+1.0);\
        a3              = (a*x*(3.0*lnzb-1.0)+1.0)/9.0;\
        r               = zl/zb;\
        hicfcio         = ((1.0-r)*a2+r*a3)/zb;\
        dhicfcio_dw     = ((1.0-r)*x+r*a)*lnzb;\
    end else begin\
        a               = z*z;\
        a2              = 3.0+z-0.25*a+0.10*z*a;\
        a3              = 2.0*z+0.75*a-0.20*a*z;\
        hicfcio         = (zb*a2+zl*a3)*w*w/6.0;\
        dhicfcio_dw     = (1+zl*w)*(1+z)*lnzb;\
    end


// NEEDED TO CALCULATE WEIGHTED ICCR COLLECTOR MINORITY CHARGE
// INPUT:
//  z : zeta_b or zeta_l
//  w : normalized injection width
// OUTPUT:
//  hicfcto     : output
//  dhicfcto_dw : derivative of output wrt w
`define HICFCT(z,w,hicfcto,dhicfcto_dw)\
    a = z*w;\
    lnz = ln(1+z*w);\
    if (a > 1.0e-6) begin\
        hicfcto     = (a - lnz)/z;\
        dhicfcto_dw = a / (1.0 + a);\
    end else begin\
        hicfcto     = 0.5 * a * w;\
        dhicfcto_dw = a;\
    end


// COLLECTOR CURRENT SPREADING CALCULATION
// collector minority charge incl. 2D/3D current spreading (TED 10/96)
// INPUT:
//  Ix                          : forward transport current component (itf)
//  I_CK                        : critical current
//  FFT_pcS                     : dependent on fthc and thcs (parameters)
// IMPLICIT INPUT:
//  ahc, latl, latb             : model parameters
//  VT                          : thermal voltage
// OUTPUT:
//  Q_fC, Q_CT: actual and ICCR (weighted) hole charge
//  T_fC, T_cT: actual and ICCR (weighted) transit time
//  Derivative dfCT_ditf not properly implemented yet
`define HICQFC(Ix,I_CK,FFT_pcS,Q_fC,Q_CT,T_fC,T_cT)\
    Q_fC            = FFT_pcS*Ix;\
    FCa             = 1.0-I_CK/Ix;\
    FCrt            = sqrt(FCa*FCa+ahc);\
    FCa_ck          = 1.0-(FCa+FCrt)/(1.0+sqrt(1.0+ahc));\
    FCdaick_ditf    = (FCa_ck-1.0)*(1-FCa)/(FCrt*Ix);\
    if(latb > latl) begin\
        FCz             = latb-latl;\
        FCxl            = 1.0+latl;\
        FCxb            = 1.0+latb;\
        if(latb > 0.01) begin\
            FCln            = ln(FCxb/FCxl);\
            FCa1            = exp((FCa_ck-1.0)*FCln);\
            FCd_a           = 1.0/(latl-FCa1*latb);\
            FCw             = (FCa1-1.0)*FCd_a;\
            FCdw_daick      = -FCz*FCa1*FCln*FCd_a*FCd_a;\
            FCa1            = ln((1.0+latb*FCw)/(1.0+latl*FCw));\
            FCda1_dw        = latb/(1.0+latb*FCw) - latl/(1.0+latl*FCw);\
        end else begin\
            FCf1            = 1.0-FCa_ck;\
            FCd_a           = 1.0/(1.0+FCa_ck*latb);\
            FCw             = FCf1*FCd_a;\
            FCdw_daick      = -1.0*FCd_a*FCd_a*FCxb*FCd_a;\
            FCa1            = FCz*FCw;\
            FCda1_dw        = FCz;\
        end\
        FCf_CT          = 2.0/FCz;\
        FCw2            = FCw*FCw;\
        FCf1            = latb*latl*FCw*FCw2/3.0+(latb+latl)*FCw2/2.0+FCw;\
        FCdf1_dw        = latb*latl*FCw2 + (latb+latl)*FCw + 1.0;\
        `HICFCI(latb,latl,FCw,FCf2,FCdf2_dw)\
        `HICFCI(latl,latb,FCw,FCf3,FCdf3_dw)\
        FCf_ci          = FCf_CT*(FCa1*FCf1-FCf2+FCf3);\
        FCdfc_dw        = FCf_CT*(FCa1*FCdf1_dw+FCda1_dw*FCf1-FCdf2_dw+FCdf3_dw);\
        FCdw_ditf       = FCdw_daick*FCdaick_ditf;\
        FCdfc_ditf      = FCdfc_dw*FCdw_ditf;\
        if(flcomp == 0.0 || flcomp == 2.1) begin\
            `HICFCT(latb,FCw,FCf2,FCdf2_dw)\
            `HICFCT(latl,FCw,FCf3,FCdf3_dw)\
            FCf_CT          = FCf_CT*(FCf2-FCf3);\
            FCdfCT_dw       = FCf_CT*(FCdf2_dw-FCdf3_dw);\
            FCdfCT_ditf     = FCdfCT_dw*FCdw_ditf;\
        end else begin\
            FCf_CT          = FCf_ci;\
            FCdfCT_ditf     = FCdfc_ditf;\
        end\
    end else begin\
        if(latb > 0.01) begin\
            FCd_a           = 1.0/(1.0+FCa_ck*latb);\
            FCw             = (1.0-FCa_ck)*FCd_a;\
            FCdw_daick      = -(1.0+latb)*FCd_a*FCd_a;\
        end else begin\
            FCw             = 1.0-FCa_ck-FCa_ck*latb;\
            FCdw_daick      = -(1.0+latb);\
        end\
        FCw2            = FCw*FCw;\
        FCz             = latb*FCw;\
        FCz_1           = 1.0+FCz;\
        FCd_f           = 1.0/(FCz_1);\
        FCf_ci          = FCw2*(1.0+FCz/3.0)*FCd_f;\
        FCdfc_dw        = 2.0*FCw*(FCz_1+FCz*FCz/3.0)*FCd_f*FCd_f;\
        FCdw_ditf       = FCdw_daick*FCdaick_ditf;\
        FCdfc_ditf      = FCdfc_dw*FCdw_ditf;\
        if(flcomp == 0.0 || flcomp == 2.1) begin\
            if (FCz > 0.001) begin\
                FCf_CT          = 2.0*(FCz_1*ln(FCz_1)-FCz)/(latb*latb*FCz_1);\
                FCdfCT_dw       = 2.0*FCw*FCd_f*FCd_f;\
            end else begin\
                FCf_CT          = FCw2*(1.0-FCz/3.0)*FCd_f;\
                FCdfCT_dw       = 2.0*FCw*(1.0-FCz*FCz/3.0)*FCd_f*FCd_f;\
            end\
            FCdfCT_ditf     = FCdfCT_dw*FCdw_ditf;\
        end else begin\
            FCf_CT          = FCf_ci;\
            FCdfCT_ditf     = FCdfc_ditf;\
        end\
    end\
    Q_CT    = Q_fC*FCf_CT;\
    Q_fC    = Q_fC*FCf_ci;\
    T_fC    = FFT_pcS*(FCf_ci+Ix*FCdfc_ditf);\
    T_cT    = FFT_pcS*(FCf_CT+Ix*FCdfCT_ditf);

// TRANSIT-TIME AND STORED MINORITY CHARGE
// INPUT:
//  itf         : forward transport current
//  I_CK        : critical current
//  T_f         : transit time    \
//  Q_f         : minority charge / for low current
// IMPLICIT INPUT:
//  tef0, gtfe, fthc, thcs, ahc, latl, latb     : model parameters
// OUTPUT:
//  T_f         : transit time    \
//  Q_f         : minority charge / transient analysis
//  T_fT        : transit time    \
//  Q_fT        : minority charge / ICCR (transfer current)
//  Q_bf        : excess base charge
`define HICQFF(itf,I_CK,T_f,Q_f,T_fT,Q_fT,Q_bf)\
    if(itf < 1.0e-6*I_CK) begin\
        Q_fT            = Q_f;\
        T_fT            = T_f;\
    end else begin\
        FFa             = I_CK/itf;\
        FFd_TfE         = tef0_t*exp(-gtfe*ln(FFa));\
        FFd_QfE         = FFd_TfE*itf/(gtfe+1.0);\
        FFT_fbS         = (1.0-fthc)*thcs_t;\
        FFx             = 1.0-FFa;\
        FFs             = sqrt(FFx*FFx+ahc);\
        FFw             = (FFx+FFs)/(1.0+sqrt(1.0+ahc));\
        FFw_2           = FFw*FFw;\
        FFd_QfB         = FFT_fbS*itf*FFw_2;\
        Q_bf            = FFd_QfB;\
        FFa_w           = FFw_2*(1.0+2.0*FFa/FFs);\
        FFd_TfB         = FFT_fbS*FFa_w;\
        FFT_pcS         = fthc*thcs_t;\
        if(latb <= 0.0 && latl <= 0.0) begin\
            FFQ_fC        = FFT_pcS*itf*FFw_2;\
            FFT_fC        = FFT_pcS*FFa_w;\
            FFQ_cT        = FFQ_fC;\
            FFT_cT        = FFT_fC;\
        end else begin\
            `HICQFC(itf,I_CK,FFT_pcS,FFQ_fC,FFQ_cT,FFT_fC,FFT_cT)\
        end\
        Q_f             = Q_f+FFd_QfB;\
        T_f             = T_f+FFd_TfB;\
        Q_fT            = Q_f+hfe*FFd_QfE+hfc*FFQ_cT;\
        T_fT            = T_f+hfe*FFd_TfE+hfc*FFT_cT;\
        Q_f             = Q_f+FFd_QfE+FFQ_fC;\
        T_f             = T_f+FFd_TfE+FFT_fC;\
    end




// IDEAL DIODE (WITHOUT CAPACITANCE):
// conductance calculation not required
// 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 > `Dexp_lim) begin\
            le      = (1 + (DIOY - `Dexp_lim));\
            DIOY    = `Dexp_lim;\
        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



// TEMPERATURE UPDATE OF JUNCTION CAPACITANCE RELATED PARAMETERS
// INPUT:
//  mostly model parameters
//  x           : zero bias junction capacitance
//  y           : junction built-in potencial
//  z           : grading co-efficient
//  w           : ratio of maximum to zero-bias value of capacitance or punch-through voltage
//  is_al       : condition factor to check what "w" stands for
//  vgeff       : band-gap voltage
// IMPLICIT INPUT:
//  VT          : thermal voltage
//  vt0,qtt0,ln_qtt0,mg : other model variables
// OUTPUT:
//  c_j_t               : temperature update of "c_j"
//  vd_t                : temperature update of "vd0"
//  w_t                 : temperature update of "w"
`define TMPHICJ(c_j,vd0,z,w,is_al,vgeff,c_j_t,vd_t,w_t)\
    if (c_j > 0.0) begin\
        vdj0    = 2*vt0*ln(exp(vd0*0.5/vt0)-exp(-0.5*vd0/vt0));\
        vdjt    = vdj0*qtt0+vgeff*(1-qtt0)-mg*VT*ln_qtt0;\
        vdt     = vdjt+2*VT*ln(0.5*(1+sqrt(1+4*exp(-vdjt/VT))));\
        vd_t    = vdt;\
        c_j_t   = c_j*exp(z*ln(vd0/vd_t));\
        if (is_al == 1) begin\
            w_t = w*vd_t/vd0;\
        end else begin\
            w_t = w;\
        end\
    end else begin\
        c_j_t   = c_j;\
        vd_t    = vd0;\
        w_t     = w;\
    end

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

//Node definitions

inout           c,b,e,s,tnode;
electrical      c,b,e,s,ci,ei,bp,bi,si;
electrical	xf1,xf2;
electrical      xf;  //RC nw

electrical      tnode;
electrical 	    n1,n2;

//Branch definitions
branch          (b,bp)          br_bbp_i;
branch          (b,bp)          br_bbp_v;
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          (bp,bi)         br_bpbi_i;
branch          (bp,bi)         br_bpbi_v;
branch          (si,s)          br_sis_i;
branch          (si,s)          br_sis_v;
branch          (bi,ei)         br_biei;
branch          (bi,ci)         br_bici;
branch          (ci,bi)         br_cibi;
branch          (ci,ei)         br_ciei;
branch          (ei,ci)         br_eici;
branch          (bp,e)          br_bpe;
branch          (b,e)           br_be;
branch          (bp,ei)         br_bpei;
branch          (bp,ci)         br_bpci;
branch          (b,ci)          br_bci;
branch          (si,ci)         br_sici;
branch          (bp,si)         br_bpsi;
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

//Noise

branch 		(n1 ) 	b_n1;
branch 		(n2 ) 	b_n2;	
	


// -- ###########################################################
// -- ###########     Parameters initialization  ################
// -- ###########################################################


//Transfer current
parameter real c10      = 2.0E-30       from [0:1]      `ATTR(info="GICCR constant" unit="A^2s");
parameter real qp0      = 2.0E-14       from (0:1]      `ATTR(info="Zero-bias hole charge" unit="Coul");
parameter real ich      = 0.0           from [0:inf)    `ATTR(info="High-current correction for 2D and 3D effects" unit="A"); //`0' signifies infinity
parameter real hfe      = 1.0           from [0:inf]    `ATTR(info="Emitter minority charge weighting factor in HBTs");
parameter real hfc      = 1.0           from [0:inf]    `ATTR(info="Collector minority charge weighting factor in HBTs");
parameter real hjei     = 1.0           from [0:100]    `ATTR(info="B-E depletion charge weighting factor in HBTs");
parameter real hjci     = 1.0           from [0:100]    `ATTR(info="B-C depletion charge weighting factor in HBTs");

//Base-Emitter diode currents
parameter real ibeis    = 1.0E-18       from [0:1]      `ATTR(info="Internal B-E saturation current" unit="A");
parameter real mbei     = 1.0           from (0:10]     `ATTR(info="Internal B-E current ideality factor");
parameter real ireis    = 0.0           from [0:1]      `ATTR(info="Internal B-E recombination saturation current" unit="A");
parameter real mrei     = 2.0           from (0:10]     `ATTR(info="Internal B-E recombination current ideality factor");
parameter real ibeps    = 0.0           from [0:1]      `ATTR(info="Peripheral B-E saturation current" unit="A");
parameter real mbep     = 1.0           from (0:10]     `ATTR(info="Peripheral B-E current ideality factor");
parameter real ireps    = 0.0           from [0:1]      `ATTR(info="Peripheral B-E recombination saturation current" unit="A");
parameter real mrep     = 2.0           from (0:10]     `ATTR(info="Peripheral B-E recombination current ideality factor");
parameter real mcf      = 1.0           from (0:10]     `ATTR(info="Non-ideality factor for III-V HBTs");

//Transit time for excess recombination current at b-c barrier
parameter real tbhrec   = 0.0           from [0:inf)    `ATTR(info="Base current recombination time constant at B-C barrier for high forward injection" unit="s");

//Base-Collector diode currents
parameter real ibcis    = 1.0E-16       from [0:1.0]    `ATTR(info="Internal B-C saturation current" unit="A");
parameter real mbci     = 1.0           from (0:10]     `ATTR(info="Internal B-C current ideality factor");
parameter real ibcxs    = 0.0           from [0:1.0]    `ATTR(info="External B-C saturation current" unit="A");
parameter real mbcx     = 1.0           from (0:10]     `ATTR(info="External B-C current ideality factor");

//Base-Emitter tunneling current
parameter real ibets    = 0.0           from [0:1]      `ATTR(info="B-E tunneling saturation current" unit="A");
parameter real abet     = 40            from [0:inf)    `ATTR(info="Exponent factor for tunneling current");
parameter integer tunode= 1             from [0:1]      `ATTR(info="Specifies the base node connection for the tunneling current"); // =1 signifies perimeter node

//Base-Collector avalanche current
parameter real favl     = 0.0           from [0:inf)    `ATTR(info="Avalanche current factor" unit="1/V");
parameter real qavl     = 0.0           from [0:inf)    `ATTR(info="Exponent factor for avalanche current" unit="Coul");
parameter real alfav    = 0.0                           `ATTR(info="Relative TC for FAVL" unit="1/K");
parameter real alqav    = 0.0                           `ATTR(info="Relative TC for QAVL" unit="1/K");

//Series resistances
parameter real rbi0     = 0.0           from [0:inf)    `ATTR(info="Zero bias internal base resistance" unit="Ohm");
parameter real rbx      = 0.0           from [0:inf)    `ATTR(info="External base series resistance" unit="Ohm");
parameter real fgeo     = 0.6557        from [0:inf]    `ATTR(info="Factor for geometry dependence of emitter current crowding");
parameter real fdqr0    = 0.0           from [-0.5:100]    `ATTR(info="Correction factor for modulation by B-E and B-C space charge layer");
parameter real fcrbi    = 0.0           from [0:1]      `ATTR(info="Ratio of HF shunt to total internal capacitance (lateral NQS effect)");
parameter real fqi      = 1.0           from [0:1]      `ATTR(info="Ration of internal to total minority charge");
parameter real re       = 0.0           from [0:inf)    `ATTR(info="Emitter series resistance" unit="Ohm");
parameter real rcx      = 0.0           from [0:inf)    `ATTR(info="External collector series resistance" unit="Ohm");

//Substrate transistor
parameter real itss     = 0.0           from [0:1.0]    `ATTR(info="Substrate transistor transfer saturation current" unit="A");
parameter real msf      = 1.0           from (0:10]     `ATTR(info="Forward ideality factor of substrate transfer current");
parameter real iscs     = 0.0           from [0:1.0]    `ATTR(info="C-S diode saturation current" unit="A");
parameter real msc      = 1.0           from (0:10]     `ATTR(info="Ideality factor of C-S diode current");
parameter real tsf      = 0.0           from [0:inf)    `ATTR(info="Transit time for forward operation of substrate transistor" unit="s");

//Intra-device substrate coupling
parameter real rsu      = 0.0           from [0:inf)    `ATTR(info="Substrate series resistance" unit="Ohm");
parameter real csu      = 0.0           from [0:inf)    `ATTR(info="Substrate shunt capacitance" unit="F");

//Depletion Capacitances
parameter real cjei0    = 1.0E-20       from [0:inf)    `ATTR(info="Internal B-E zero-bias depletion capacitance" unit="F");
parameter real vdei     = 0.9           from (0:10]     `ATTR(info="Internal B-E built-in potential" unit="V");
parameter real zei      = 0.5	          from (0:1]      `ATTR(info="Internal B-E grading coefficient");
parameter real ajei     = 2.5           from [1:inf)    `ATTR(info="Ratio of maximum to zero-bias value of internal B-E capacitance");
parameter real cjep0    = 1.0E-20       from [0:inf)    `ATTR(info="Peripheral B-E zero-bias depletion capacitance" unit="F");
parameter real vdep     = 0.9           from (0:10]     `ATTR(info="Peripheral B-E built-in potential" unit="V");
parameter real zep      = 0.5           from (0:1]      `ATTR(info="Peripheral B-E grading coefficient");
parameter real ajep     = 2.5           from [1:inf)    `ATTR(info="Ratio of maximum to zero-bias value of peripheral B-E capacitance");
parameter real cjci0    = 1.0E-20       from [0:inf)    `ATTR(info="Internal B-C zero-bias depletion capacitance" unit="F");
parameter real vdci     = 0.7           from (0:10]     `ATTR(info="Internal B-C built-in potential" unit="V");
parameter real zci      = 0.4           from (0:1]      `ATTR(info="Internal B-C grading coefficient");
parameter real vptci    = 100           from (0:100]    `ATTR(info="Internal B-C punch-through voltage" unit="V");
parameter real cjcx0    = 1.0E-20       from [0:inf)    `ATTR(info="External B-C zero-bias depletion capacitance" unit="F");
parameter real vdcx     = 0.7           from (0:10]     `ATTR(info="External B-C built-in potential" unit="V");
parameter real zcx      = 0.4           from (0:1]      `ATTR(info="External B-C grading coefficient");
parameter real vptcx    = 100           from (0:100]    `ATTR(info="External B-C punch-through voltage" unit="V");
parameter real fbcpar   = 0.0           from [0:1]      `ATTR(info="Partitioning factor of parasitic B-C cap");
parameter real fbepar   = 1.0           from [0:1]      `ATTR(info="Partitioning factor of parasitic B-E cap");
parameter real cjs0     = 0.0           from [0:inf)    `ATTR(info="C-S zero-bias depletion capacitance" unit="F");
parameter real vds      = 0.6           from (0:10]     `ATTR(info="C-S built-in potential" unit="V");
parameter real zs       = 0.5           from (0:1]      `ATTR(info="C-S grading coefficient");
parameter real vpts     = 100           from (0:100]    `ATTR(info="C-S punch-through voltage" unit="V"); 

//Diffusion Capacitances
parameter real t0       = 0.0           from [0:inf)    `ATTR(info="Low current forward transit time at VBC=0V" unit="s");
parameter real dt0h     = 0.0           from (-inf:inf) `ATTR(info="Time constant for base and B-C space charge layer width modulation" unit="s");
parameter real tbvl     = 0.0           from [0:inf)    `ATTR(info="Time constant for modelling carrier jam at low VCE" unit="s");
parameter real tef0     = 0.0           from [0:inf)    `ATTR(info="Neutral emitter storage time" unit="s");
parameter real gtfe     = 1.0           from (0:10]     `ATTR(info="Exponent factor for current dependence of neutral emitter storage time");
parameter real thcs     = 0.0           from [0:inf)    `ATTR(info="Saturation time constant at high current densities" unit="s");
parameter real ahc      = 0.1           from (0:10]     `ATTR(info="Smoothing factor for current dependence of base and collector transit time");
parameter real fthc     = 0.0           from [0:1]      `ATTR(info="Partitioning factor for base and collector portion");
parameter real rci0     = 150           from (0:inf)    `ATTR(info="Internal collector resistance at low electric field" unit="Ohm");
parameter real vlim     = 0.5           from (0:10]     `ATTR(info="Voltage separating ohmic and saturation velocity regime" unit="V");
parameter real vces     = 0.1           from [0:1]      `ATTR(info="Internal C-E saturation voltage" unit="V");
parameter real vpt      = 100.0         from (0:inf]    `ATTR(info="Collector punch-through voltage" unit="V"); // `0' signifies infinity
parameter real tr       = 0.0           from [0:inf)    `ATTR(info="Storage time for inverse operation" unit="s");
  
//Isolation Capacitances
parameter real cbepar   = 0.0           from [0:inf)    `ATTR(info="Total parasitic B-E capacitance" unit="F");
parameter real cbcpar   = 0.0           from [0:inf)    `ATTR(info="Total parasitic B-C capacitance" unit="F");

//Non-quasi-static Effect
parameter real alqf     = 0.0           from [0:1]      `ATTR(info="Factor for additional delay time of minority charge");
parameter real alit     = 0.0           from [0:1]      `ATTR(info="Factor for additional delay time of transfer current");
parameter integer flnqs = 0             from [0:1]      `ATTR(info="Flag for turning on and off of vertical NQS effect");
  
//Noise
parameter real kf       = 0.0           from [0:inf)    `ATTR(info="Flicker noise coefficient");
parameter real af       = 2.0           from (0:10]     `ATTR(info="Flicker noise exponent factor");
parameter integer cfbe  = -1            from [-2:-1]    `ATTR(info="Flag for determining where to tag the flicker noise source");


//Lateral Geometry Scaling (at high current densities)
parameter real latb     = 0.0           from [0:inf)    `ATTR(info="Scaling factor for collector minority charge in direction of emitter width");
parameter real latl     = 0.0           from [0:inf)    `ATTR(info="Scaling factor for collector minority charge in direction of emitter length");

//Temperature dependence
parameter real vgb      = 1.17          from (0:10]     `ATTR(info="Bandgap voltage extrapolated to 0 K" unit="V");
parameter real alt0     = 0.0                           `ATTR(info="First order relative TC of parameter T0" unit="1/K");
parameter real kt0      = 0.0                           `ATTR(info="Second order relative TC of parameter T0");
parameter real zetaci   = 0.0           from [-10:10]   `ATTR(info="Temperature exponent for RCI0");
parameter real alvs     = 0.0                           `ATTR(info="Relative TC of saturation drift velocity" unit="1/K");
parameter real alces    = 0.0                           `ATTR(info="Relative TC of VCES" unit="1/K");
parameter real zetarbi  = 0.0           from [-10:10]   `ATTR(info="Temperature exponent of internal base resistance");
parameter real zetarbx  = 0.0           from [-10:10]   `ATTR(info="Temperature exponent of external base resistance");
parameter real zetarcx  = 0.0           from [-10:10]   `ATTR(info="Temperature exponent of external collector resistance");
parameter real zetare   = 0.0           from [-10:10]   `ATTR(info="Temperature exponent of emitter resistance");
parameter real zetacx   = 1.0           from [-10:10]   `ATTR(info="Temperature exponent of mobility in substrate transistor transit time");
parameter real vge      = 1.17          from (0:10]     `ATTR(info="Effective emitter bandgap voltage" unit="V");
parameter real vgc      = 1.17          from (0:10]     `ATTR(info="Effective collector bandgap voltage" unit="V");
parameter real vgs      = 1.17          from (0:10]     `ATTR(info="Effective substrate bandgap voltage" unit="V");
parameter real f1vg     =-1.02377e-4                    `ATTR(info="Coefficient K1 in T-dependent band-gap equation");
parameter real f2vg     = 4.3215e-4                     `ATTR(info="Coefficient K2 in T-dependent band-gap equation");
parameter real zetact   = 3.0           from [-10:10]   `ATTR(info="Exponent coefficient in transfer current temperature dependence");
parameter real zetabet  = 3.5           from [-10:10]   `ATTR(info="Exponent coefficient in B-E junction current temperature dependence");
parameter real alb      = 0.0                           `ATTR(info="Relative TC of forward current gain for V2.1 model" unit="1/K");

//Self-Heating
parameter integer flsh  = 0             from [0:2]      `ATTR(info="Flag for turning on and off self-heating effect");
parameter real rth      = 0.0           from [0:inf)    `ATTR(info="Thermal resistance" unit="K/W");
parameter real cth      = 0.0           from [0:inf)    `ATTR(info="Thermal capacitance" unit="J/W");

//Compatibility with V2.1
parameter real flcomp = 0.0             from [0:inf)    `ATTR(info="Flag for compatibility with v2.1 model (0=v2.1)");

//Circuit simulator specific parameters
parameter real tnom     = 27.0                          `ATTR(info="Temperature at which parameters are specified" unit="C");
parameter real dt       = 0.0                           `ATTR(info="Temperature change w.r.t. chip temperature for particular transistor" unit="K");


//
//======================== Transistor model formulation ===================
//

    //Declaration of variables

    real _circuit_gmin;

    //Temperature and drift
    real VT,Tdev,qtt0,ln_qtt0,r_VgVT,V_gT,dT,k;
    real ireis_t,ibeis_t,ibcxs_t,ibcis_t,iscs_t,cjci0_t;
    real cjs0_t,rci0_t,vlim_t,vces_t,thcs_t,tef0_t,rbi0_t;
    real rbx_t,rcx_t,re_t,t0_t,vdei_t,vdci_t,vpts_t,itss_t,tsf_t;
    real c10_t,cjei0_t,qp0_t,vdcx_t,vptcx_t,cjcx01_t,cjcx02_t;
    real qjcx0_t_i,qjcx0_t_ii,cratio_t,c_dummy;
    real ibeps_t,ireps_t,cjep0_t;
    real ajei_t,qavl_t,favl_t,ibets_t,abet_t,vptci_t,vdep_t,ajep_t,zetatef;
    real k1,k2,dvg0,vge_t,vgb_t,vgbe_t,vds_t,vt0,Tnom,Tamb,a,avs;
    real zetabci,zetabcxt,zetasct,vgbe0,mg,vgb_t0,vge_t0,vgbe_t0,vgbc0,vgsc0;
    real cbcpar1,cbcpar2,cbepar2,cbepar1,Oich,Ovpt,Otbhrec;

    //Charges, capacitances and currents
    real Qjci,Qjei,Qjep;
    real it,ibei,irei,ibci,ibep,irep,ibh_rec;
    real Qdei,Qdci,qrbi;
    real ibet,iavl;
    real ijbcx,ijsc,Qjs,HSUM,HSI_Tsu,Qdsu;

    //Base resistance and self-heating power
    real rbi,pterm;

    //Variables for macro TMPHICJ
    real  vdj0,vdjt,vdt;

    //Model initialization
    real k10,k20,C_1;

    //Model evaluation
    real Cjci,Cjcit,cc,Cjei,Cjep;
    real itf,itr,Tf,Tr,VT_f,i_0f,i_0r,a_bpt,Q_0,Q_p,Q_bpt;
    real Orci0_t,b_q,Q_fC,T_fC,T_cT,I_Tf1,T_f0,Q_fT,T_fT,Q_bf;
    real ICKa,d1;
    real A,a_h,Q_pT,d_Q,d_Q0;
    real Qf,Cdei,Qr,Cdci,Crbi;
    real ick,vc,vceff,cjcx01,cjcx02,HSa,HSb;
    integer l_it;

    //Variables for macros
    real  DIOY,le;//HICDIO
    real  FFT_fbS,FFa,FFx,FFs,FFw,FFw_2,FFd_QfB,FFd_TfB,FFT_pcS,FFQ_fC,FFT_fC,FFQ_cT,FFT_cT,FFd_TfE,FFd_QfE,FFa_w;//HICQFF
    real  FCz,FCw2,FCf1,FCf2,FCf3,FCf_ci,FCz_1,FCa1,FCa_ck,FCxl,FCxb;//HICQFC
    real  FCd_a,FCdaick_ditf,FCa,FCw,FCdw_daick,FCdfc_dw,FCdw_ditf,FCdfc_ditf,FCf_CT,FCdfCT_ditf,FCrt,FCln,lnz,FCda1_dw,FCdf1_dw,FCdf2_dw,FCdf3_dw,FCd_f,FCdfCT_dw;//HICQFC
    real  Dz_r,Dv_p,DV_f,DC_max,DC_c,Da,Dv_e,De,De_1,Dv_j1,Dv_r,De_2,Dv_j2,Dv_j4,DQ_j1,DQ_j2,DQ_j3,DCln1,DCln2,Dz1,Dzr1,DC_j1,DC_j2,DC_j3;//QJMOD
    real  DFV_f,DFv_e,DFv_j,DFb,DFQ_j,DFs_q,DFs_q2,DFdvj_dv,DFC_j1;//QJMODF
    real  z,a2,a3,r,x;//HICFCI
    real zb,zl,lnzb,w,hicfcio,dhicfcio_dw; //HICFCT

    //Noise
    real fourkt,twoq,flicker_Pwr;
    real thermal_Rbx,thermal_Rbi,thermal_Rcx,thermal_Re,betad,betan,betadin,betadc,icn,icn1,icn2;

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

    real pocce,czz;
    real Qz_nom,f_QR,ETA,Qz0,fQz;
    real v_bord,v_q,U0,av,avl;
    real cV_f,cv_e,cs_q,cs_q2,cv_j,cdvj_dv;
    real a_eg,ab,aa;

    //end of variables

analog begin

`MODEL begin : Model_initialization

    Tnom    = tnom+`P_CELSIUS0;
    Tamb    = $temperature;
    vt0     = `P_K*Tnom /`P_Q;
    k10     = f1vg*Tnom*ln(Tnom);
    k20     = f2vg*Tnom;
    avs     = alvs*Tnom;
    vgb_t0  = vgb+k10+k20;
    vge_t0  = vge+k10+k20;
    vgbe_t0 = (vgb_t0+vge_t0)/2;
    vgbe0   = (vgb+vge)/2;
    vgbc0   = (vgb+vgc)/2;
    vgsc0   = (vgs+vgc)/2;
    mg      = 3-`P_Q*f1vg/`P_K;
    zetabci = mg+1-zetaci;
    zetabcxt= mg+1-zetacx;
    zetasct = mg-1.5;

    //Depletion capacitance splitting at b-c junction
    //Capacitances at peripheral and external base node
    C_1    = (1.0-fbcpar)*(cjcx0+cbcpar);
    if (C_1 >= cbcpar) begin
        cbcpar1 = cbcpar;
        cbcpar2 = 0.0;
        cjcx01  = C_1-cbcpar;
        cjcx02  = cjcx0-cjcx01;
    end else begin
        cbcpar1 = C_1;
        cbcpar2 = cbcpar-cbcpar1;
        cjcx01  = 0.0;
        cjcx02  = cjcx0;
    end

    //Parasitic b-e capacitance partitioning: No temperature dependence
    cbepar2 = fbepar*cbepar;
    cbepar1 = cbepar-cbepar2;

    //Avoid devide-by-zero and define infinity other way
    //High current correction for 2D and 3D effects
    if (ich != 0.0) begin
        Oich    = 1.0/ich;
    end else begin
        Oich    = 0.0;
    end

    //Base current recombination time constant at b-c barrier
    if (tbhrec != 0.0) begin
        Otbhrec = 1.0/tbhrec;
    end else begin
        Otbhrec = 0.0;
    end

    // Temperature and resulting parameter drift
    if (flsh==0 || rth < `MIN_R) begin : Thermal_updat_without_self_heating
        Tdev    = Tamb+dt;
        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
        VT      = `P_K*Tdev /`P_Q;
        dT      = Tdev-Tnom;
        qtt0    = Tdev/Tnom;
        ln_qtt0 = ln(qtt0);
        k1      = f1vg*Tdev*ln(Tdev);
        k2      = f2vg*Tdev;
        vgb_t   = vgb+k1+k2;
        vge_t   = vge+k1+k2;
        vgbe_t  = (vgb_t+vge_t)/2;

        //Internal b-e junction capacitance
        `TMPHICJ(cjei0,vdei,zei,ajei,1,vgbe0,cjei0_t,vdei_t,ajei_t)

        if (flcomp == 0.0 || flcomp == 2.1) begin
            V_gT     = 3.0*VT*ln_qtt0 + vgb*(qtt0-1.0);
            r_VgVT   = V_gT/VT;
            //Internal b-e diode saturation currents
            a        = mcf*r_VgVT/mbei - alb*dT;
            ibeis_t  = ibeis*exp(a);
            a        = mcf*r_VgVT/mrei - alb*dT;
            ireis_t  = ireis*exp(a);
            a        = mcf*r_VgVT/mbep - alb*dT;
            //Peripheral b-e diode saturation currents
            ibeps_t  = ibeps*exp(a);
            a        = mcf*r_VgVT/mrep - alb*dT;
            ireps_t  = ireps*exp(a);
            //Internal b-c diode saturation current
            a       = r_VgVT/mbci;
            ibcis_t = ibcis*exp(a);
            //External b-c diode saturation currents
            a       = r_VgVT/mbcx;
            ibcxs_t = ibcxs*exp(a);
            //Saturation transfer current for substrate transistor
            a       = r_VgVT/msf;
            itss_t  = itss*exp(a);
            //Saturation current for c-s diode
            a       = r_VgVT/msc;
            iscs_t  = iscs*exp(a);
            //Zero bias hole charge
            a        = vdei_t/vdei;
            qp0_t    = qp0*(1.0+0.5*zei*(1.0-a));
            //Voltage separating ohmic and saturation velocity regime
            a = vlim*(1.0-alvs*dT)*exp(zetaci*ln_qtt0);
            k = (a-VT)/VT;
            if (k < `LN_EXP_LIMIT) begin
                vlim_t = VT + VT*ln(1.0+exp(k));
            end else begin
                vlim_t = a;
            end
            //Neutral emitter storage time
            a        = 1.0+alb*dT;
            k        = 0.5*(a+sqrt(a*a+0.01));
            tef0_t   = tef0*qtt0/k;
        end else begin
            //Internal b-e diode saturation currents
            ibeis_t  = ibeis*exp(zetabet*ln_qtt0+vge/VT*(qtt0-1));
            ireis_t  = ireis*exp(0.5*mg*ln_qtt0+0.5*vgbe0/VT*(qtt0-1));
            //Peripheral b-e diode saturation currents
            ibeps_t  = ibeps*exp(zetabet*ln_qtt0+vge/VT*(qtt0-1));
            ireps_t  = ireps*exp(0.5*mg*ln_qtt0+0.5*vgbe0/VT*(qtt0-1));
            //Internal b-c diode saturation currents
            ibcis_t = ibcis*exp(zetabci*ln_qtt0+vgc/VT*(qtt0-1));
            //External b-c diode saturation currents
            ibcxs_t = ibcxs*exp(zetabcxt*ln_qtt0+vgc/VT*(qtt0-1));
            //Saturation transfer current for substrate transistor
            itss_t  = itss*exp(zetasct*ln_qtt0+vgc/VT*(qtt0-1));
            //Saturation current for c-s diode
            iscs_t  = iscs*exp(zetasct*ln_qtt0+vgs/VT*(qtt0-1));
            //Zero bias hole charge
            a       = exp(zei*ln(vdei_t/vdei));
            qp0_t   = qp0*(2.0-a);
            //Voltage separating ohmic and saturation velocity regime
            vlim_t  = vlim*exp((zetaci-avs)*ln_qtt0);
            //Neutral emitter storage time
            zetatef = zetabet-zetact-0.5;
            dvg0    = vgb-vge;
            tef0_t  = tef0*exp(zetatef*ln_qtt0-dvg0/VT*(qtt0-1));
        end

        //GICCR prefactor
        c10_t   = c10*exp(zetact*ln_qtt0+vgb/VT*(qtt0-1));

        // Low-field internal collector resistance
        rci0_t  = rci0*exp(zetaci*ln_qtt0);

        //Voltage separating ohmic and saturation velocity regime
        //vlim_t  = vlim*exp((zetaci-avs)*ln_qtt0);

        //Internal c-e saturation voltage
        vces_t  = vces*(1+alces*dT);


        //Internal b-c diode saturation current
        //ibcis_t = ibcis*exp(zetabci*ln_qtt0+vgc/VT*(qtt0-1));

        //Internal b-c junction capacitance
        `TMPHICJ(cjci0,vdci,zci,vptci,0,vgbc0,cjci0_t,vdci_t,vptci_t)

        //Low-current forward transit time
        t0_t    = t0*(1+alt0*dT+kt0*dT*dT);

        //Saturation time constant at high current densities
        thcs_t  = thcs*exp((zetaci-1)*ln_qtt0);


        //Avalanche caurrent factors
        favl_t  = favl*exp(alfav*dT);
        qavl_t  = qavl*exp(alqav*dT);

        //Zero bias internal base resistance
        rbi0_t  = rbi0*exp(zetarbi*ln_qtt0);


        //Peripheral b-e junction capacitance
        `TMPHICJ(cjep0,vdep,zep,ajep,1,vgbe0,cjep0_t,vdep_t,ajep_t)

        //Tunneling current factors
        begin : HICTUN_T
//            real a_eg,ab,aa;
            ab      = 1.0;
            aa      = 1.0;
            a_eg=vgbe_t0/vgbe_t;
            if(tunode==1 && cjep0 > 0.0 && vdep >0.0) begin
                ab      = (cjep0_t/cjep0)*sqrt(a_eg)*vdep_t*vdep_t/(vdep*vdep);
                aa      = (vdep/vdep_t)*(cjep0/cjep0_t)*pow(a_eg,-1.5);
            end else if (tunode==0 && cjei0 > 0.0 && vdei >0.0) begin
                ab      = (cjei0_t/cjei0)*sqrt(a_eg)*vdei_t*vdei_t/(vdei*vdei);
                aa      = (vdei/vdei_t)*(cjei0/cjei0_t)*pow(a_eg,-1.5);
            end
            ibets_t = ibets*ab;
            abet_t  = abet*aa;
        end


        //Temperature mapping for tunneling current is done inside HICTUN

        `TMPHICJ(1.0,vdcx,zcx,vptcx,0,vgbc0,cratio_t,vdcx_t,vptcx_t)
        cjcx01_t=cratio_t*cjcx01;
        cjcx02_t=cratio_t*cjcx02;


        //External b-c diode saturation currents
        //ibcxs_t       = ibcxs*exp(zetabcxt*ln_qtt0+vgc/VT*(qtt0-1));


        //Constant external series resistances
        rcx_t   = rcx*exp(zetarcx*ln_qtt0);
        rbx_t   = rbx*exp(zetarbx*ln_qtt0);
        re_t    = re*exp(zetare*ln_qtt0);

        //Forward transit time in substrate transistor
        tsf_t   = tsf*exp((zetacx-1.0)*ln_qtt0);

        //Capacitance for c-s junction
        `TMPHICJ(cjs0,vds,zs,vpts,0,vgsc0,cjs0_t,vds_t,vpts_t)

    end // of Thermal_update_without_self_heating

end     //of Model_initialization

if (flsh!=0 && rth >= `MIN_R) begin : Thermal_update_with_self_heating
    Tdev    = Tamb+dt+V(br_sht);
    // Limit temperature to avoid FPEs 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
    VT      = `P_K*Tdev /`P_Q;
    dT      = Tdev-Tnom;
    qtt0    = Tdev/Tnom;
    ln_qtt0 = ln(qtt0);
    k1      = f1vg*Tdev*ln(Tdev);
    k2      = f2vg*Tdev;
    vgb_t   = vgb+k1+k2;
    vge_t   = vge+k1+k2;
    vgbe_t  = (vgb_t+vge_t)/2;

    //Internal b-e junction capacitance
    `TMPHICJ(cjei0,vdei,zei,ajei,1,vgbe0,cjei0_t,vdei_t,ajei_t)

    if (flcomp == 0.0 || flcomp == 2.1) begin
        V_gT     = 3.0*VT*ln_qtt0 + vgb*(qtt0-1.0);
        r_VgVT   = V_gT/VT;
        //Internal b-e diode saturation currents
        a        = mcf*r_VgVT/mbei - alb*dT;
        ibeis_t  = ibeis*exp(a);
        a        = mcf*r_VgVT/mrei - alb*dT;
        ireis_t  = ireis*exp(a);
        a        = mcf*r_VgVT/mbep - alb*dT;
        //Peripheral b-e diode saturation currents
        ibeps_t  = ibeps*exp(a);
        a        = mcf*r_VgVT/mrep - alb*dT;
        ireps_t  = ireps*exp(a);
        //Internal b-c diode saturation current
        a       = r_VgVT/mbci;
        ibcis_t = ibcis*exp(a);
        //External b-c diode saturation currents
        a       = r_VgVT/mbcx;
        ibcxs_t = ibcxs*exp(a);
        //Saturation transfer current for substrate transistor
        a       = r_VgVT/msf;
        itss_t  = itss*exp(a);
        //Saturation current for c-s diode
        a       = r_VgVT/msc;
        iscs_t  = iscs*exp(a);
        //Zero bias hole charge
        a        = vdei_t/vdei;
        qp0_t    = qp0*(1.0+0.5*zei*(1.0-a));
        //Voltage separating ohmic and saturation velocity regime
        a = vlim*(1.0-alvs*dT)*exp(zetaci*ln_qtt0);
        k = (a-VT)/VT;
        if (k < `LN_EXP_LIMIT) begin
            vlim_t = VT + VT*ln(1.0+exp(k));
        end else begin
            vlim_t = a;
        end
        //Neutral emitter storage time
        a        = 1.0+alb*dT;
        k        = 0.5*(a+sqrt(a*a+0.01));
        tef0_t   = tef0*qtt0/k;
    end else begin
        //Internal b-e diode saturation currents
        ibeis_t  = ibeis*exp(zetabet*ln_qtt0+vge/VT*(qtt0-1));
        ireis_t  = ireis*exp(0.5*mg*ln_qtt0+0.5*vgbe0/VT*(qtt0-1));
        //Peripheral b-e diode saturation currents
        ibeps_t  = ibeps*exp(zetabet*ln_qtt0+vge/VT*(qtt0-1));
        ireps_t  = ireps*exp(0.5*mg*ln_qtt0+0.5*vgbe0/VT*(qtt0-1));
        //Internal b-c diode saturation currents
        ibcis_t = ibcis*exp(zetabci*ln_qtt0+vgc/VT*(qtt0-1));
        //External b-c diode saturation currents
        ibcxs_t = ibcxs*exp(zetabcxt*ln_qtt0+vgc/VT*(qtt0-1));
        //Saturation transfer current for substrate transistor
        itss_t  = itss*exp(zetasct*ln_qtt0+vgc/VT*(qtt0-1));
        //Saturation current for c-s diode
        iscs_t  = iscs*exp(zetasct*ln_qtt0+vgs/VT*(qtt0-1));
        //Zero bias hole charge
        a       = exp(zei*ln(vdei_t/vdei));
        qp0_t   = qp0*(2.0-a);
        //Voltage separating ohmic and saturation velocity regime
        vlim_t  = vlim*exp((zetaci-avs)*ln_qtt0);
        //Neutral emitter storage time
        zetatef = zetabet-zetact-0.5;
        dvg0    = vgb-vge;
        tef0_t  = tef0*exp(zetatef*ln_qtt0-dvg0/VT*(qtt0-1));
    end

    //GICCR prefactor
    c10_t   = c10*exp(zetact*ln_qtt0+vgb/VT*(qtt0-1));

    // Low-field internal collector resistance
    rci0_t  = rci0*exp(zetaci*ln_qtt0);

    //Voltage separating ohmic and saturation velocity regime
    //vlim_t  = vlim*exp((zetaci-avs)*ln_qtt0);

    //Internal c-e saturation voltage
    vces_t  = vces*(1+alces*dT);


    //Internal b-c diode saturation current
    //ibcis_t = ibcis*exp(zetabci*ln_qtt0+vgc/VT*(qtt0-1));

    //Internal b-c junction capacitance
    `TMPHICJ(cjci0,vdci,zci,vptci,0,vgbc0,cjci0_t,vdci_t,vptci_t)

    //Low-current forward transit time
    t0_t    = t0*(1+alt0*dT+kt0*dT*dT);

    //Saturation time constant at high current densities
    thcs_t  = thcs*exp((zetaci-1)*ln_qtt0);


    //Avalanche caurrent factors
    favl_t  = favl*exp(alfav*dT);
    qavl_t  = qavl*exp(alqav*dT);

    //Zero bias internal base resistance
    rbi0_t  = rbi0*exp(zetarbi*ln_qtt0);


    //Peripheral b-e junction capacitance
    `TMPHICJ(cjep0,vdep,zep,ajep,1,vgbe0,cjep0_t,vdep_t,ajep_t)

    //Tunneling current factors
    if (V(br_bpei) < 0.0 || V(br_biei) < 0.0) begin : HICTUN_T
//        real a_eg,ab,aa;
        ab      = 1.0;
        aa      = 1.0;
        a_eg=vgbe_t0/vgbe_t;
        if(tunode==1 && cjep0 > 0.0 && vdep >0.0) begin
            ab      = (cjep0_t/cjep0)*sqrt(a_eg)*vdep_t*vdep_t/(vdep*vdep);
            aa      = (vdep/vdep_t)*(cjep0/cjep0_t)*pow(a_eg,-1.5);
        end else if (tunode==0 && cjei0 > 0.0 && vdei >0.0) begin
            ab      = (cjei0_t/cjei0)*sqrt(a_eg)*vdei_t*vdei_t/(vdei*vdei);
            aa      = (vdei/vdei_t)*(cjei0/cjei0_t)*pow(a_eg,-1.5);
        end
        ibets_t = ibets*ab;
        abet_t  = abet*aa;
    end


    //Temperature mapping for tunneling current is done inside HICTUN

    `TMPHICJ(1.0,vdcx,zcx,vptcx,0,vgbc0,cratio_t,vdcx_t,vptcx_t)
    cjcx01_t=cratio_t*cjcx01;
    cjcx02_t=cratio_t*cjcx02;


    //External b-c diode saturation currents
    //ibcxs_t       = ibcxs*exp(zetabcxt*ln_qtt0+vgc/VT*(qtt0-1));


    //Constant external series resistances
    rcx_t   = rcx*exp(zetarcx*ln_qtt0);
    rbx_t   = rbx*exp(zetarbx*ln_qtt0);
    re_t    = re*exp(zetare*ln_qtt0);

    //Forward transit time in substrate transistor
    tsf_t   = tsf*exp((zetacx-1.0)*ln_qtt0);

    //Capacitance for c-s junction
    `TMPHICJ(cjs0,vds,zs,vpts,0,vgsc0,cjs0_t,vds_t,vpts_t)

end     //of Thermal_update_with_self_heating


begin : Model_evaluation

    //Intrinsic transistor
    //Internal base currents across b-e junction
    `HICDIO(ibeis,ibeis_t,mbei,V(br_biei),ibei)
    `HICDIO(ireis,ireis_t,mrei,V(br_biei),irei)

    //HICCR: begin

    //Inverse of low-field internal collector resistance: needed in HICICK
    Orci0_t = 1.0/rci0_t;

    //Initialization
    //Transfer current, minority charges and transit times

    Tr      = tr;
    VT_f    = mcf*VT;
    i_0f    = c10_t * limexp(V(br_biei)/VT_f);
    i_0r    = c10_t * limexp(V(br_bici)/VT);

    //Internal b-e and b-c junction capacitances and charges
    //`QJMODF(cjei0_t,vdei_t,zei,ajei_t,V(br_biei),Qjei)
    //Cjei    = ddx(Qjei,V(bi));
    `QJMODF(cjei0_t,vdei_t,zei,ajei_t,V(br_biei),Cjei,Qjei)

    //`HICJQ(cjci0_t,vdci_t,zci,vptci_t,V(br_bici),Qjci)
    //Cjci    = ddx(Qjci,V(bi));
    `HICJQ(cjci0_t,vdci_t,zci,vptci_t,V(br_bici),Cjci,Qjci)

    //Hole charge at low bias
    a_bpt   = 0.05;
    Q_0     = qp0_t + hjei*Qjei + hjci*Qjci;
    Q_bpt   = a_bpt*qp0_t;
    b_q     = Q_0/Q_bpt-1;
    Q_0     = Q_bpt*(1+(b_q +sqrt(b_q*b_q+1.921812))/2);

    //Transit time calculation at low current density
    if(cjci0_t > 0.0) begin : CJMODF
//        real cV_f,cv_e,cs_q,cs_q2,cv_j,cdvj_dv;
        cV_f    = vdci_t*(1.0-exp(-ln(2.4)/zci));
        cv_e    = (cV_f-V(br_bici))/VT;
        cs_q    = sqrt(cv_e*cv_e+1.921812);
        cs_q2   = (cv_e+cs_q)*0.5;
        cv_j    = cV_f-VT*cs_q2;
        cdvj_dv = cs_q2/cs_q;
        Cjcit   = cjci0_t*exp(-zci*ln(1.0-cv_j/vdci_t))*cdvj_dv+2.4*cjci0_t*(1.0-cdvj_dv);
    end else begin
        Cjcit   = 0.0;
    end
    if(Cjcit > 0.0) begin
        cc      = cjci0_t/Cjcit;
    end else begin
        cc      = 1.0;
    end
    T_f0    = t0_t+dt0h*(cc-1.0)+tbvl*(1/cc-1.0);

    //Effective collector voltage
    vc      = V(br_ciei)-vces_t;

    //Critical current for onset of high-current effects
    begin : HICICK
        Ovpt    = 1.0/vpt;
        a       = vc/VT;
        d1      = a-1;
        vceff   = (1.0+((d1+sqrt(d1*d1+1.921812))/2))*VT;
        a       = vceff/vlim_t;
        ick     = vceff*Orci0_t/sqrt(1.0+a*a);
        ICKa    = (vceff-vlim_t)*Ovpt;
        ick     = ick*(1.0+0.5*(ICKa+sqrt(ICKa*ICKa+1.0e-3)));
    end

    //Initial formulation of forward and reverse component of transfer current
    Q_p     = Q_0;
    if (T_f0 > 0.0 || Tr > 0.0) begin
        A       = 0.5*Q_0;
        Q_p     = A+sqrt(A*A+T_f0*i_0f+Tr*i_0r);
    end
    I_Tf1   =i_0f/Q_p;
    a_h     = Oich*I_Tf1;
    itf     = I_Tf1*(1.0+a_h);
    itr     = i_0r/Q_p;

    //Initial formulation of forward transit time, diffusion, GICCR and excess b-c charge
    Q_bf    = 0.0;
    Tf      = T_f0;
    Qf      = T_f0*itf;
    `HICQFF(itf,ick,Tf,Qf,T_fT,Q_fT,Q_bf)

    //Initial formulation of reverse diffusion charge
    Qr      = Tr*itr;

    //Preparation for iteration to get total hole charge and related variables
    l_it    = 0;
    if(Qf > `RTOLC*Q_p || a_h > `RTOLC) begin
        //Iteration for Q_pT is required for improved initial solution
        Qf      = sqrt(T_f0*itf*Q_fT);
        Q_pT    = Q_0+Qf+Qr;
        d_Q     = Q_pT;
        while (abs(d_Q) >= `RTOLC*abs(Q_pT) && l_it <= `l_itmax) begin
            d_Q0    = d_Q;
            I_Tf1   = i_0f/Q_pT;
            a_h     = Oich*I_Tf1;
            itf     = I_Tf1*(1.0+a_h);
            itr     = i_0r/Q_pT;
            Tf      = T_f0;
            Qf      = T_f0*itf;
            `HICQFF(itf,ick,Tf,Qf,T_fT,Q_fT,Q_bf)
            Qr      = Tr*itr;
            if(Oich == 0.0) begin
                a       = 1.0+(T_fT*itf+Qr)/Q_pT;
            end else begin
                a       = 1.0+(T_fT*I_Tf1*(1.0+2.0*a_h)+Qr)/Q_pT;
            end
            d_Q     = -(Q_pT-(Q_0+Q_fT+Qr))/a;
            //Limit maximum change of Q_pT
            a       = abs(0.3*Q_pT);
            if(abs(d_Q) > a) begin
                if (d_Q>=0) begin
                    d_Q     = a;
                end else begin
                    d_Q     = -a;
                end
            end
            Q_pT    = Q_pT+d_Q;
            l_it    = l_it+1;
        end //while

        I_Tf1   = i_0f/Q_pT;
        a_h     = Oich*I_Tf1;
        itf     = I_Tf1*(1.0+a_h);
        itr     = i_0r/Q_pT;

        //Final transit times, charges and transport current components
        Tf      = T_f0;
        Qf      = T_f0*itf;
        `HICQFF(itf,ick,Tf,Qf,T_fT,Q_fT,Q_bf)
        Qr      = Tr*itr;

    end //if

    //NQS effect implemented with LCR networks
    //Once the delay in ITF is considered, IT_NQS is calculated afterwards

    it      = itf-itr;

    //Diffusion charges for further use
    Qdei    = Qf;
    Qdci    = Qr;


    //High-frequency emitter current crowding (lateral NQS)
    Cdei    = -1*ddx(Qdei,V(ei));
    Cdci    = -1*ddx(Qdci,V(ci));
    Crbi    = fcrbi*(Cjei+Cjci+Cdei+Cdci);
    qrbi    = Crbi*V(br_bpbi_v);

//    qrbi	= fcrbi*(Qjei+Qjci+Qdei+Qdci);

    //HICCR: end

    //Internal base current across b-c junction
    `HICDIO(ibcis,ibcis_t,mbci,V(br_bici),ibci)

    //Avalanche current
    if((V(br_bici) < 0.0) && (favl_t > 0.0) && (cjci0_t > 0.0)) begin : HICAVL
//        real v_bord,v_q,U0,av,avl;
        v_bord  = vdci_t-V(br_bici);
        v_q     = qavl_t/Cjci;
        U0      = qavl_t/cjci0_t;
        if(v_bord > U0) begin
                av      = favl_t*exp(-v_q/U0);
                avl     = av*(U0+(1.0+v_q/U0)*(v_bord-U0));
        end else begin
            avl     = favl_t*v_bord*exp(-v_q/v_bord);
        end
        iavl    = itf*avl;
    end else begin
        iavl    = 0.0;
    end

    //Excess base current from recombination at the b-c barrier
    ibh_rec = Q_bf*Otbhrec;

    //Internal base resistance
    if(rbi0_t > 0.0) begin : HICRBI
//        real Qz_nom,f_QR,ETA,Qz0,fQz;
        // Consideration of conductivity modulation
        // To avoid convergence problem hyperbolic smoothing used
        f_QR    = (1+fdqr0)*qp0_t;
        Qz0     = Qjei+Qjci+Qf;
        Qz_nom  = 1+Qz0/f_QR;
        fQz     = 0.5*(Qz_nom+sqrt(Qz_nom*Qz_nom+0.01));
        rbi     = rbi0_t/fQz;
        // Consideration of emitter current crowding
        if( ibei > 0.0) begin
            ETA     = rbi*ibei*fgeo/VT;
            if(ETA < 1.0e-6) begin
                rbi     = rbi*(1.0-0.5*ETA);
            end else begin
                rbi     = rbi*ln(1.0+ETA)/ETA;
            end
        end
        // Consideration of peripheral charge
        if(Qf > 0.0) begin
            rbi     = rbi*(Qjei+Qf*fqi)/(Qjei+Qf);
        end
    end else begin
        rbi     = 0.0;
    end

    //Base currents across peripheral b-e junction
    `HICDIO(ibeps,ibeps_t,mbep,V(br_bpei),ibep)
    `HICDIO(ireps,ireps_t,mrep,V(br_bpei),irep)

    //Peripheral b-e junction capacitance and charge
    `QJMODF(cjep0_t,vdep_t,zep,ajep_t,V(br_bpei),Cjep,Qjep)

    //Tunelling current
    if (V(br_bpei) <0.0 || V(br_biei) < 0.0) begin : HICTUN
//        real pocce,czz;
        if(tunode==1 && cjep0_t > 0.0 && vdep_t >0.0) begin
            pocce   = exp((1-1/zep)*ln(Cjep/cjep0_t));
            czz     = -(V(br_bpei)/vdep_t)*ibets_t*pocce;
            ibet    = czz*exp(-abet_t/pocce);
            end else if (tunode==0 && cjei0_t > 0.0 && vdei_t >0.0) begin
                pocce   = exp((1-1/zei)*ln(Cjei/cjei0_t));
                czz     = -(V(br_biei)/vdei_t)*ibets_t*pocce;
                ibet    = czz*exp(-abet_t/pocce);
            end else begin
            ibet    = 0.0;
        end
    end else begin
        ibet    = 0.0;
    end


    //Depletion capacitance and charge at peripheral b-c junction (bp,ci)
    `HICJQ(cjcx02_t,vdcx_t,zcx,vptcx_t,V(br_bpci),c_dummy,qjcx0_t_ii)

    //Base currents across peripheral b-c junction (bp,ci)
    `HICDIO(ibcxs,ibcxs_t,mbcx,V(br_bpci),ijbcx)

    //Depletion capacitance and charge at external b-c junction (b,ci)
    `HICJQ(cjcx01_t,vdcx_t,zcx,vptcx_t,V(br_bci),c_dummy,qjcx0_t_i)

    //Depletion substrate capacitance and charge at s-c junction (si,ci)
    `HICJQ(cjs0_t,vds_t,zs,vpts_t,V(br_sici),c_dummy,Qjs)

    //Parasitic substrate transistor transfer current and diffusion charge
    if(itss > 0.0) begin : Sub_Transfer
        HSUM    = msf*VT;
        HSa     = limexp(V(br_bpci)/HSUM);
        HSb     = limexp(V(br_sici)/HSUM);
        HSI_Tsu = itss_t*(HSa-HSb);
        if(tsf > 0.0) begin
            Qdsu    = tsf_t*itss_t*HSa;
        end else begin
            Qdsu    = 0.0;
        end
    end else begin
        HSI_Tsu = 0.0;
        Qdsu    = 0.0;
    end

    // Current gain computation for correlated noise implementation
    betad=ibei;
    if (betad > 0.0) begin
        betadin=betad;
        betan=it;
        betadc=betan/betad;
    end else begin
        betadc=0.0;
    end

    //Diode current for s-c junction (si,ci)
    `HICDIO(iscs,iscs_t,msc,V(br_sici),ijsc)

    //Self-heating calculation
    if (flsh == 1 && rth >= `MIN_R) begin
        pterm   =  V(br_ciei)*it + (vdci_t-V(br_bici))*iavl;
    end else if (flsh == 2 && rth >= `MIN_R) begin
        pterm   =  V(br_ciei)*it + (vdci_t-V(br_bici))*iavl + ibei*V(br_biei) + ibci*V(br_bici) + ibep*V(br_bpei) + ijbcx*V(br_bpci) + ijsc*V(br_sici);
        if (rbi >= `MIN_R) begin
            pterm   = pterm + V(br_bpbi_i)*V(br_bpbi_i)/rbi;
        end
        if (re_t >= `MIN_R) begin
            pterm   = pterm + V(br_eie_i)*V(br_eie_i)/re_t;
        end
        if (rcx_t >= `MIN_R) begin
            pterm   = pterm + V(br_cic_i)*V(br_cic_i)/rcx_t;
        end
        if (rbx_t >= `MIN_R) begin
            pterm   = pterm + V(br_bbp_i)*V(br_bbp_i)/rbx_t;
        end
    end

    Itxf	= itf;
    Qdeix	= Qdei;
    // 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);               //for RC nw
        fact  = t0/Tf;                 //for RC nw  
        Ixf   = (Vxf - Qdei)*fact;     //for RC nw
        Qxf   = alqf*Vxf*t0;           //for RC nw
        Qdeix = Vxf;                   //for RC nw
    end else begin
        Ixf1	=  V(br_bxf1);
        Ixf2	=  V(br_bxf2);
        Qxf1	=  0;
        Qxf2	=  0;

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

end     //of Model_evaluation

begin : Load_sources

    I(br_biei)      <+ _circuit_gmin*V(br_biei);
    I(br_bici)      <+ _circuit_gmin*V(br_bici);

    I(br_bci)       <+ ddt(qjcx0_t_i);
    I(br_bci)       <+ ddt(cbcpar1*V(br_bci));
    I(br_bpci)      <+ ddt(cbcpar2*V(br_bpci));
    if (rbx >= `MIN_R) begin
        I(br_bbp_i)     <+ V(br_bbp_i)/rbx_t;
    end else begin
//        V(br_bbp_v)     <+ 0.0;
        I(br_bbp_i)     <+ V(br_bbp_i)/`MIN_R;
    end
    if(rbi0 >= `MIN_R) begin
        I(br_bpbi_i)    <+ V(br_bpbi_i)/rbi;
        I(br_bpbi_i)    <+ ddt(qrbi);
    end else begin
//        V(br_bpbi_v)    <+ 0.0;
        I(br_bpbi_i)    <+ V(br_bpbi_i)/`MIN_R;
    end
    if (tunode==1.0) begin
        I(br_bpei)      <+ -ibet;
    end else begin
        I(br_biei)      <+ -ibet;
    end
    I(br_bpei)      <+ ibep;
    I(br_bpei)      <+ irep;
    I(br_bpei)      <+ ddt(Qjep);
    I(br_biei)      <+ ibei;
    I(br_biei)      <+ irei;
    I(br_biei)      <+ ibh_rec;
    I(br_biei)      <+ ddt(Qdeix+Qjei);
    I(br_bpsi)      <+ HSI_Tsu;
    I(br_bpci)      <+ ijbcx;
    I(br_bpci)      <+ ddt(qjcx0_t_ii+Qdsu);
    I(br_be)        <+ ddt(cbepar1*V(br_be));
    I(br_bpe)       <+ ddt(cbepar2*V(br_bpe));
    I(br_bici)      <+ ibci-iavl;
    I(br_bici)      <+ ddt(Qdci+Qjci);
    I(br_sici)      <+ ijsc;
    I(br_sici)      <+ ddt(Qjs);
    I(br_ciei)      <+ Itxf;
    I(br_eici)      <+ itr;
    if (rcx >= `MIN_R) begin
        I(br_cic_i)     <+ V(br_cic_i)/rcx_t;
    end else begin
//        V(br_cic_v)     <+ 0.0;
        I(br_cic_i)     <+ V(br_cic_i)/`MIN_R;
    end
    if (re >= `MIN_R) begin
        I(br_eie_i)     <+ V(br_eie_i)/re_t;
    end else begin
//        V(br_eie_v)     <+ 0.0;
        I(br_eie_i)     <+ V(br_eie_i)/`MIN_R;
    end
    if(rsu >= `MIN_R) begin
        I(br_sis_i)     <+ V(br_sis_i)/rsu;
        I(br_sis_i)     <+ ddt(csu*V(br_sis_i));
    end else begin
//        V(br_sis_v)     <+ 0.0;
        I(br_sis_i)     <+ V(br_sis_i)/`MIN_R;
    end

    // Following code is an intermediate solution (if branch contribution is not supported):
    // ******************************************
    if(flsh == 0 || rth < `MIN_R) begin
          I(br_sht)       <+ V(br_sht)/`MIN_R;
    end else begin
          I(br_sht)       <+ V(br_sht)/rth-pterm;
          I(br_sht)       <+ ddt(cth*V(br_sht));
    end

    // ******************************************

    // For simulators having no problem with V(br_sht) <+ 0.0
    // with external thermal node, follwing code may be used.
    // Note that external thermal node should remain accessible
    // even without self-heating.
    // ********************************************
    //if(flsh == 0 || rth < `MIN_R) begin
    //    V(br_sht)       <+ 0.0;
    //end else begin
    //    I(br_sht)       <+ V(br_sht)/rth-pterm;
    //    I(br_sht)       <+ ddt(cth*V(br_sht));
    //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;         //for RC nw
    I(br_cxf) <+  ddt(Qxf);    //for RC nw

end //of Load_sources


`NOISE begin : Noise_sources

    //Thermal noise
    fourkt  = 4.0 * `P_K * Tdev;
    if(rbx >= `MIN_R) begin
        I(br_bbp_i)     <+ white_noise(fourkt/rbx_t,    "thermal");
    end
    if(rbi0 >= `MIN_R) begin
        I(br_bpbi_i)    <+ white_noise(fourkt/rbi,      "thermal");
    end
    if(rcx >= `MIN_R) begin
        I(br_cic_i)     <+ white_noise(fourkt/rcx_t,    "thermal");
    end
    if(re >= `MIN_R) begin
        I(br_eie_i)     <+ white_noise(fourkt/re_t,     "thermal");
    end
    if(rsu >= `MIN_R) begin
        I(br_sis_i)     <+ white_noise(fourkt/rsu,      "thermal");
    end

    //Flicker noise : Fully correlated between the perimeter and internal base-node
    flicker_Pwr     = kf*pow((ibei+ibep),af);
    if (cfbe == -1) begin
        I(br_biei)      <+ flicker_noise(flicker_Pwr,1.0);
    end else begin
        I(br_bpei)      <+ flicker_noise(flicker_Pwr,1.0);
    end

    //Shot noise

    twoq    = 2.0 * `P_Q;
    //     I(br_ciei)      <+ white_noise(twoq*it,         "shot");
    I(br_cibi)      <+ white_noise(twoq*iavl,       "shot");

    //     I(br_biei)      <+ white_noise(twoq*ibei,       "shot");

    I(br_bici)      <+ white_noise(twoq*abs(ibci),  "shot");

    I(br_bpei)      <+ white_noise(twoq*ibep,       "shot");

    I(br_bpci)      <+ white_noise(twoq*abs(ijbcx), "shot");

    I(br_sici)      <+ white_noise(twoq*abs(ijsc),  "shot");

    // Code section for correlated noise
    // Please turn-off this code section by "//" in order to run the code with Spectre

//    I(b_n1)		<+ white_noise(2 * `P_Q * ibei,		"shot");
//    I(b_n1)		<+ V(b_n1);
//    I(b_n2) 		<+ white_noise(2 * `P_Q * it,		"shot");
//    I(b_n2)		<+ V(b_n2);
//
//    I(bi,ei) 		<+ V(b_n1); 
//    I(ci,ei) 		<+ V(b_n2)+ddt((betadc/2)*alit*Tf*alit*Tf*ddt(V(b_n2))); 
//    I(ci,ei)		<+ betadc*ddt(-(Tf*alit)*V(b_n1));                       
end     //of Noise_sources

end //analog
endmodule
