//
//	Verilog-A definition of VBIC
//
//	Version:	1.2
//	Revision:	1.1
//
//	 1. Functions redefined to work properly in Verilog-A, previously
//	    they were substituted by special code to generate .va files
//	    without function calls.
//	 2. Noise declarations added, consistent with 1.1.5 version.
//	 3. qb and qbp calculations protected from having argument of the
//	    square-root go less than zero (this required very unusual
//	    conditions and parameters to occur, and needless to say
//	    this happened some times).
//	 4. Parameter names changed to lower case, this is what they are
//	    in simulator implementations.
//	 5. log mapped to ln, this was done in generated code, but in
//	    Verilog-A log is log_10, not log_e.
//	 6. The avalanche function was modified to protect against
//	    numerical problems for forward and low reverse bias.
//	 7. Hard limits on the local temperature were added. This is
//	    required to avoid numerical problems.
//	 8. The excess phase network has been defined in a simpler
//	    manner, but is exactly equivalent to the previous version.
//	    It is now done without an inductor. This is the transformation:
//	    Itzf-V(xf2)-j*w*C*V(xf1)=0
//	    j*w*L*V(xf2)+V(xf2)-V(xf1)=0
//	    where C=TD and L=TD/3 is what was done.
//	    Itzf=V(xf2)+j*w*TD*V(xf1)
//	    V(xf1)=j*w*(TD/3)*V(xf2)+V(xf2)
//	    node xf1: two VCCS (Itzf, V(xf2)) plus TD capacitor
//	    node xf2: one VCCS (V(xf1)), TD/3 capacitor, 1Ohm resistor
//	 9. gmin added explicitly
//	10. pnjmaxi added explicitly, diode type currents linearized
//	    for values greater than pnjmaxi
//	11. Fixed implementation of shrink factors
//

`include "discipline.h"

`ifdef insideADMS
  `define P(prop) (* prop *)
  `define PGIVEN(p)		$given(p)
  `define GMIN			$options("OPTgmin")
  `define PNJMAXI		$options("OPTpnjmaxi")
  `define SHRINKL		$shrinkl("q")
  `define SHRINKA		$shrinka("q")
  `define SCALE			$scale
  `define DERIVATE(x,y)		$derivate(x,y)
  `define INITIAL_MODEL		@(initial_model)
  `define INITIAL_INSTANCE	@(initial_instance)
  `define FINAL_STEP
`else
  `define P(prop)
  `define PGIVEN(p)		p
  `define GMIN			gmin
  `define PNJMAXI		pnjmaxi
  `define SHRINKL		1.0
  `define SHRINKA		1.0
  `define SCALE			1.0
  `define DERIVATE(x,y)		1.0
  `define INITIAL_MODEL		@(initial_step("op","dc","ac","noise"))
  `define INITIAL_INSTANCE	@(initial_step("op","dc","ac","noise"))
  `define FINAL_STEP		@(final_step)
`endif

`define	KB	1.380662e-23		// Boltzmann constant (J/K)
`define	QQ	1.602189e-19		// mag. of electronic charge (C)
`define	TABS	2.731500e+02		// 0C in K

`define	NPN	+1
`define	PNP	-1

`define expLinA(expv,V,Vmax,a) \
	if (V<Vmax) begin \
		mac1	=  V*a; \
		expv	=  exp(mac1); \
	end else begin \
		mac1	=  Vmax*a; \
		expl	=  exp(mac1); \
		expv	=  expl*(1.0+(V-Vmax)*a); \
	end
`define expLin(expv,V,Vmax) \
	if (V<Vmax) \
		expv	=  exp(V); \
	else begin \
		expl	=  exp(Vmax); \
		expv	=  expl*(1.0+(V-Vmax)); \
	end
`define psibi(P_T,P,EA,Vtv,rT) \
	psiio = 2.0*(Vtv/rT)*ln(exp(0.5*P*rT/Vtv)-exp(-0.5*P*rT/Vtv)); \
	psiin = psiio*rT-3.0*Vtv*ln(rT)-EA*(rT-1.0); \
	P_T   = psiin+2.0*Vtv*ln(0.5*(1.0+sqrt(1.0+4.0*exp(-psiin/Vtv))));
`define qj(qj,V,P,M,FC,A) \
	dv0   = -P*FC; \
	if (A<=0.0) begin \
		dvh =  V+dv0; \
		if (dvh>0.0) begin \
			pwq = pow((1.0-FC),(-1.0-M)); \
			qlo = P*(1.0-pwq*(1.0-FC)*(1.0-FC))/(1.0-M); \
			qhi = dvh*(1.0-FC+0.5*M*dvh/P)*pwq; \
		end else begin \
			qlo = P*(1.0-pow((1.0-V/P),(1.0-M)))/(1.0-M); \
			qhi = 0.0; \
		end \
		qj  = qlo+qhi; \
	end else begin \
		mv0 =  sqrt(dv0*dv0+4*A*A); \
		vl0 = -0.5*(dv0+mv0); \
		q0  = -P*pow((1.0-vl0/P),(1.0-M))/(1.0-M); \
		dv  =  V+dv0; \
		mv  =  sqrt(dv*dv+4*A*A); \
		vl  =  0.5*(dv-mv)-dv0; \
		qlo = -P*pow((1.0-vl/P),(1.0-M))/(1.0-M); \
		qj  =  qlo+pow((1.0-FC),(-M))*(V-vl+vl0)-q0; \
	end
`define qjrt(qj,V,P,M,FC,A,VRT,ART) \
	dv0   = -P*FC; \
	if (A<=0.0) begin \
		dvh =  V+dv0; \
		if (dvh>0.0) begin \
			pwq = pow((1.0-FC),(-1.0-M)); \
			qlo = P*(1.0-pwq*(1.0-FC)*(1.0-FC))/(1.0-M); \
			qhi = dvh*(1.0-FC+0.5*M*dvh/P)*pwq; \
		end else begin \
			if ((VRT>0.0)&&(V<-VRT)) begin \
				qlo = P*(1.0-pow((1.0+VRT/P),(1.0-M))*(1.0-((1.0-M)*(V+VRT))/(P+VRT)))/(1.0-M); \
			end else begin \
				qlo = P*(1.0-pow((1.0-V/P),(1.0-M)))/(1.0-M); \
			end \
			qhi = 0.0; \
		end \
		qj  = qlo+qhi; \
	end else begin \
		if ((VRT>0.0)&&(ART>0.0)) begin \
			vn0  =  (VRT+dv0)/(VRT-dv0); \
			vnl0 =  2.0*vn0/(sqrt((vn0-1.0)*(vn0-1)+4*A*A)+sqrt((vn0+1.0)*(vn0+1)+4*ART*ART)); \
			vl0  =  0.5*(vnl0*(VRT-dv0)-VRT-dv0); \
			qlo0 =  P*(1.0-pow((1.0-vl0/P),(1.0-M)))/(1.0-M); \
			vn   =  (2*V+VRT+dv0)/(VRT-dv0); \
			vnl  =  2.0*vn/(sqrt((vn-1.0)*(vn-1)+4*A*A)+sqrt((vn+1.0)*(vn+1)+4*ART*ART)); \
			vl   =  0.5*(vnl*(VRT-dv0)-VRT-dv0); \
			qlo  =  P*(1.0-pow((1.0-vl/P),(1.0-M)))/(1.0-M); \
			sel  =  0.5*(vnl+1.0); \
			crt  =  pow((1.0+VRT/P),(-M)); \
			cmx  =  pow((1.0+dv0/P),(-M)); \
			cl   =  (1.0-sel)*crt+sel*cmx; \
			ql   =  (V-vl+vl0)*cl; \
			qj   =  ql+qlo-qlo0; \
		end else begin \
			mv0  =  sqrt(dv0*dv0+4*A*A); \
			vl0  = -0.5*(dv0+mv0); \
			q0   = -P*pow((1.0-vl0/P),(1.0-M))/(1.0-M); \
			dv   =  V+dv0; \
			mv   =  sqrt(dv*dv+4*A*A); \
			vl   =  0.5*(dv-mv)-dv0; \
			qlo  = -P*pow((1.0-vl/P),(1.0-M))/(1.0-M); \
			qj   =  qlo+pow((1.0-FC),(-M))*(V-vl+vl0)-q0; \
		end \
	end
`define avalm(avalm,V,P,M,AV1,AV2,VmaxExp) \
	vminm = pow((0.02*(AV2+1.0)),(1.0/(1.01-M))); \
	vl    = 0.5*(sqrt((P-V-vminm)*(P-V-vminm)+0.01)+(P-V-vminm))+vminm; \
	mac1  = -AV2*pow(vl,(M-1.0)); \
	if (mac1<VmaxExp) \
		expi	=  exp(mac1); \
	else begin \
		expl	=  exp(VmaxExp); \
		expi	=  expl*(1.0+(mac1-VmaxExp)); \
	end \
	avalm = AV1*vl*expi;

//
//	There are 8 separate versions of VBIC, defined as follows:
//
//				#Elect	Electro	Excess
//		Name		Terms	Thermal	Phase
//		=============	=======	=======	======
//		vbic_3T_it_cf	3	no	no
//		vbic_3T_it_xf	3	no	yes
//		vbic_3T_et_cf	3	yes	no
//		vbic_3T_et_xf	3	yes	yes
//		vbic_4T_it_cf	4	no	no
//		vbic_4T_it_xf	4	no	yes
//		vbic_4T_et_cf	4	yes	no
//		vbic_4T_et_xf	4	yes	yes
//
//	These can be selected by appropriate specification of
//	the following `define text macros. Note that the specification
//	of a 3- or 4-terminal model relates to the number of
//	electrical terminals, and does not include the local temperature
//	node for the electrothermal version of the model.
//
//	There are two separate versions of each of the above,
//	with and without homotopy. When the homotopy is included
//	this code will not work in Verilog-A, but will be
//	handled by the VBIC code generator properly.
//

//`define ThreeTerminal		// default is FourTerminal
//`define ExcessPhase		// default is ConstantPhase
//`define ElectroThermal	// default is IsoThermal
//`define WithHomotopy		// default is NoHomotopy

//
//	Start of VBIC model code
//

	////////////////////////////////////////////////////
	// General Information for Simulator
	////////////////////////////////////////////////////

`ifdef VERSION
`else
  `define VERSION "unknown"
`endif

module VBICSELFT(c,b,e,s)
  (*
    version=`VERSION
    spice:prefix="q"
    spice:ntype="npn"
    spice:ptype="pnp"
    spice:level="3"
    spice:terminal:min="4"
    spice:terminal:max="4"
    description="Vertical Bipolar Transistor"
    info="Level 3 BJT: vbic selft"
    mica:isdefault="0"
  *);

//
//	Node definitions
//

	inout		c,b;				// external nodes
	inout		e,s;				// external nodes

	// external nodes
	electrical c (*
		spice:name="cnode"
		mica:reprint="yes"
		info="external collector node                   " *);
	electrical b (*
		spice:name="bnode"
		mica:reprint="yes"
		info="external base node                        " *);
	electrical e (*
		spice:name="enode"
		mica:reprint="yes"
		info="external emitter node                     " *);
	electrical s (*
		spice:name="snode"
		mica:reprint="yes"
		info="external substrate node                   " *);
	// electrothermal node
	electrical dt (*
		spice:name="dtnode"
		info="local temperature rise node               " *);
	// internal nodes
	electrical cx (*
		spice:name="cxnode"
		info="internal extrinsic collector node         " *);
	electrical ci (*
		spice:name="cinode"
		info="internal intrinsic collector node         " *);
	electrical bx (*
		spice:name="bxnode"
		info="internal extrinsic base node              " *);
	electrical bi (*
		spice:name="binode"
		info="internal intrinsic base node              " *);
	electrical ei (*
		spice:name="einode"
		info="internal intrinsic emitter node           " *);
	electrical bp (*
		spice:name="bpnode"
		info="internal parasitic base node              " *);
	electrical si (*
		spice:name="sinode"
		info="internal substrate node                   " *);
	electrical xf1 (*
		spice:name="xf1node"
		info="excess phase node xf1                     " *);
	electrical xf2 (*
		spice:name="xf2node"
		info="excess phase node xf2                     " *);

//
//	Branch definitions
//

	branch (b ,e )		b_be;			//           base-emit
	branch (b ,c )		b_bc;			//           base-coll
	branch (bi,ei)		b_bei;			// intrinsic base-emit
	branch (bx,ei)		b_bex;			// extrinsic base-emit
	branch (bi,ci)		b_bci;			// intrinsic base-coll
	branch (bi,cx)		b_bcx;			// extrinsic base-coll
	branch (ci,ei)		b_cei;			// intrinsic coll-emit
	branch (ei,ci)		b_eci;			// intrinsic emit-coll
	branch (bx,bp)		b_bep;			// parasitic base-emit
	branch (e ,ei)		b_re;			// emit resistance
	branch (c ,cx)		b_rcx;			// coll resistance, extrinsic
	branch (cx,ci)		b_rci;			// coll resistance, intrinsic
	branch (b ,bx)		b_rbx;			// base resistance, extrinsic
	branch (bx,bi)		b_rbi;			// base resistance, intrinsic
	branch (bp,cx)		b_rbp;			// base resistance, parasitic
	branch (si,bp)		b_bcp;			// parasitic base-coll
	branch (bx,si)		b_cep;			// parasitic coll-emit
	branch (s ,si)		b_rs;			// subs resistance
	branch (dt)		b_rth;			// local thermal branch
	branch (dt)		b_ith;			// local thermal branch
	branch (xf1)		b_xf1;
	branch (xf1)		c_xf1;
	branch (xf2)		b_xf2;
	branch (xf2)		c_xf2;

//
//	Parameter definitions
//

	parameter	real	tnom	=  27.0
	  `P(info="nominal parameter measurement temperature " unit="C");
	parameter	real	rcx	=   0.0		from[0.0:inf]
	  `P(info="extrinsic collector resistance            " unit="Ohm");
	parameter	real	rci	=   0.0		from[0.0:inf]
	  `P(info="intrinsic collector resistance            " unit="Ohm");
	parameter	real	vo	=   0.0		from[0.0:inf]
	  `P(info="epi drift saturation voltage              " unit="V");
	parameter	real	gamm	=   0.0		from[0.0:inf]
	  `P(info="epi doping parameter                      " );
	parameter	real	hrcf	=   0.0		from[0.0:inf]
	  `P(info="high current collector resistance factor  " );
	parameter	real	rbx	=   0.0		from[0.0:inf]
	  `P(info="extrinsic base resistance                 " unit="Ohm");
	parameter	real	rbi	=   0.0		from[0.0:inf]
	  `P(info="intrinsic base resistance                 " unit="Ohm");
	parameter	real	re	=   0.0		from[0.0:inf]
	  `P(info="extrinsic emitter resistance              " unit="Ohm");
	parameter	real	rs	=   0.0		from[0.0:inf]
	  `P(info="extrinsic substrate resistance            " unit="Ohm");
	parameter	real	rbp	=   0.0		from[0.0:inf]
	  `P(info="parasitic transistor base resistance      " unit="Ohm");
	parameter	real	is	=   1.0e-16	from(0.0:inf]
	  `P(info="transport saturation current              " unit="A");
	parameter	real	nf	=   1.0		from(0.0:inf]
	  `P(info="fwd emission coefficient (ideality factor)" );
	parameter	real	nr	=   1.0		from(0.0:inf]
	  `P(info="rev emission coefficient (ideality factor)" );
	parameter	real	fc	=   0.9		from[0.0:1.0)
	  `P(info="forward bias depletion capacitance limit  " );
	parameter	real	cbeo	=   0.0		from[0.0:inf]
	  `P(info="extrinsic b-e overlap capacitance         " unit="F");
	parameter	real	cje	=   0.0		from[0.0:inf]
	  `P(info="zero-bias b-e depletion capacitance       " unit="F");
	parameter	real	pe	=   0.75	from(0.0:inf]
	  `P(info="b-e built-in potential                    " unit="V");
	parameter	real	me	=   0.33	from(0.0:1.0]
	  `P(info="b-e grading coefficient                   " );
	parameter	real	aje	=  -0.5
	  `P(info="b-e capacitance smoothing factor          " );
	parameter	real	cbco	=   0.0		from[0.0:inf]
	  `P(info="extrinsic b-c overlap capacitance         " unit="F");
	parameter	real	cjc	=   0.0		from[0.0:inf]
	  `P(info="zero-bias intrinsic b-c depletion cap     " unit="F");
	parameter	real	qco	=   0.0		from[0.0:inf]
	  `P(info="epi charge parameter                      " unit="C");
	parameter	real	cjep	=   0.0		from[0.0:inf]
	  `P(info="zero-bias extrinsic b-c depletion cap     " unit="F");
	parameter	real	pc	=   0.75	from(0.0:inf]
	  `P(info="b-c built-in potential                    " unit="V");
	parameter	real	mc	=   0.33	from(0.0:1.0]
	  `P(info="b-c grading coefficient                   " );
	parameter	real	ajc	=  -0.5
	  `P(info="b-c capacitance smoothing factor          " );
	parameter	real	cjcp	=   0.0		from[0.0:inf]
	  `P(info="zero-bias c-s depletion capacitance       " unit="F");
	parameter	real	ps	=   0.75	from(0.0:inf]
	  `P(info="c-s built-in potential                    " unit="V");
	parameter	real	ms	=   0.33	from(0.0:1.0]
	  `P(info="c-s grading coefficient                   " );
	parameter	real	ajs	=  -0.5
	  `P(info="c-s capacitance smoothing factor          " );
	parameter	real	ibei	=   1.0e-18	from(0.0:inf]
	  `P(info="ideal b-e saturation current              " unit="A");
	parameter	real	wbe	=   1.0		from[0.0:1.0]
	  `P(info="partitioning of Ibe/Ibex and Qbe/Qbex     " );
	parameter	real	nei	=   1.0		from(0.0:inf]
	  `P(info="ideal b-e emission coefficient            " );
	parameter	real	iben	=   0.0		from[0.0:inf]
	  `P(info="non-ideal b-e saturation current          " unit="A");
	parameter	real	nen	=   2.0		from(0.0:inf]
	  `P(info="non-ideal b-e emission coefficient        " );
	parameter	real	ibci	=   1.0e-16	from(0.0:inf]
	  `P(info="ideal b-c saturation current              " unit="A");
	parameter	real	nci	=   1.0		from(0.0:inf]
	  `P(info="ideal b-c emission coefficient            " );
	parameter	real	ibcn	=   0.0		from[0.0:inf]
	  `P(info="non-ideal b-c saturation current          " unit="A");
	parameter	real	ncn	=   2.0		from(0.0:inf]
	  `P(info="non-ideal b-c emission coefficient        " );
	parameter	real	avc1	=   0.0		from[0.0:inf]
	  `P(info="b-c weak avalanche parameter 1            " unit="/V");
	parameter	real	avc2	=   0.0		from[0.0:inf]
	  `P(info="b-c weak avalanche parameter 2            " );
	parameter	real	isp	=   0.0		from[0.0:inf]
	  `P(info="parasitic transport saturation current    " unit="A");
	parameter	real	wsp	=   1.0		from[0.0:1.0]
	  `P(info="partitioning of Iccp between Vbep and Vbci" );
	parameter	real	nfp	=   1.0		from(0.0:inf]
	  `P(info="parasitic emission coeff (ideality fctr)  " );
	parameter	real	ibeip	=   0.0		from[0.0:inf]
	  `P(info="ideal parasitic b-e saturation current    " unit="A");
	parameter	real	ibenp	=   0.0		from[0.0:inf]
	  `P(info="non-ideal parasitic b-e saturation current" unit="A");
	parameter	real	ibcip	=   0.0		from[0.0:inf]
	  `P(info="ideal parasitic b-c saturation current    " unit="A");
	parameter	real	ncip	=   1.0		from(0.0:inf]
	  `P(info="ideal parasitic b-c emission coefficient  " );
	parameter	real	ibcnp	=   0.0		from[0.0:inf]
	  `P(info="non-ideal parasitic b-c saturation current" unit="A");
	parameter	real	ncnp	=   2.0		from(0.0:inf]
	  `P(info="non-ideal parasitic b-c emission coeff    " );
	parameter	real	vef	=   0.0		from[0.0:inf]
	  `P(info="forward Early voltage (zero=infinite)     " unit="V");
	parameter	real	ver	=   0.0		from[0.0:inf]
	  `P(info="reverse Early voltage (zero=infinite)     " unit="V");
	parameter	real	ikf	=   0.0		from[0.0:inf]
	  `P(info="forward knee current (zero=infinite)      " unit="A");
	parameter	real	ikr	=   0.0		from[0.0:inf]
	  `P(info="reverse knee current (zero=infinite)      " unit="A");
	parameter	real	ikp	=   0.0		from[0.0:inf]
	  `P(info="parasitic knee current (zero=infinite)    " unit="A");
	parameter	real	tf	=   0.0		from[0.0:inf]
	  `P(info="forward transit time                      " unit="s");
	parameter	real	qtf	=   0.0		from[0.0:inf]
	  `P(info="variation of tf with base-width modulation" );
	parameter	real	xtf	=   0.0		from[0.0:inf]
	  `P(info="tf bias dependence coefficient            " );
	parameter	real	vtf	=   0.0		from[0.0:inf]
	  `P(info="tf coefficient of Vbci dependence         " unit="V");
	parameter	real	itf	=   0.0		from[0.0:inf]
	  `P(info="tf coefficient of Ic dependence           " unit="A");
	parameter	real	tr	=   0.0		from[0.0:inf]
	  `P(info="reverse transit time                      " unit="s");
	parameter	real	td	=   0.0		from[0.0:inf]
	  `P(info="forward excess-phase delay time           " unit="s");
	parameter	real	kfn	=   0.0		from[0.0:inf]
	  `P(info="b-e flicker noise constant                " );
	parameter	real	afn	=   1.0		from(0.0:inf]
	  `P(info="b-e flicker noise current exponent        " );
	parameter	real	bfn	=   1.0		from(0.0:inf]
	  `P(info="b-e flicker noise 1/f dependence          " );
	parameter	real	xre	=   0.0
	  `P(info="temperature exponent of re                " );
	parameter	real	xrbi	=   0.0
	  `P(info="temperature exponent of rbi               " );
	parameter	real	xrci	=   0.0
	  `P(info="temperature exponent of rci               " );
	parameter	real	xrs	=   0.0
	  `P(info="temperature exponent of rs                " );
	parameter	real	xvo	=   0.0
	  `P(info="temperature exponent of vo                " );
	parameter	real	ea	=   1.12
	  `P(info="activation energy for is                  " unit="V");
	parameter	real	eaie	=   1.12
	  `P(info="activiation energy for ibei               " unit="V");
	parameter	real	eaic	=   1.12
	  `P(info="activiation energy for ibci and ibeip     " unit="V");
	parameter	real	eais	=   1.12
	  `P(info="activiation energy for ibcip              " unit="V");
	parameter	real	eane	=   1.12
	  `P(info="activiation energy for iben               " unit="V");
	parameter	real	eanc	=   1.12
	  `P(info="activiation energy for ibcn and ibenp     " unit="V");
	parameter	real	eans	=   1.12
	  `P(info="activiation energy for ibcnp              " unit="V");
	parameter	real	xis	=   3.0
	  `P(info="temperature exponent of is                " );
	parameter	real	xii	=   3.0
	  `P(info="temp exponent of ibei, ibci, ibeip, ibcip " );
	parameter	real	xin	=   3.0
	  `P(info="temp exponent of iben, ibcn, ibenp, ibcnp " );
	parameter	real	tnf	=   0.0
	  `P(info="temperature exponent of nf and nr         " unit="/C");
	parameter	real	tavc	=   0.0
	  `P(info="temperature exponent of avc2              " unit="/C");
	parameter	real	rth	=   0.0		from[0.0:inf]
	  `P(info="thermal resistance                        " unit="C/W");
	parameter	real	cth	=   0.0		from[0.0:inf]
	  `P(info="thermal capacitance                       " unit="J/C");
	parameter	real	vrt	=   0.0		from[0.0:inf]
	  `P(info="reach-through voltage for Cbc limiting    " unit="V");
	parameter	real	art	=   0.1		from(0.0:inf]
	  `P(info="smoothing parameter for reach-through     " );
	parameter	real	ccso	=   0.0		from[0.0:inf]
	  `P(info="extrinsic c-s overlap capacitance         " unit="F");
	parameter	real	qbm	=   0.0
	  `P(info="base charge model selection parameter     " );
	parameter	real	nkf	=   0.5		from(0.0:inf]
	  `P(info="high current beta roll-off parameter      " );
	parameter	real	xikf	=   0.0
	  `P(info="temperature exponent of ikf               " );
	parameter	real	xrcx	=   0.0
	  `P(info="temperature exponent of rcx               " );
	parameter	real	xrbx	=   0.0
	  `P(info="temperature exponent of rbx               " );
	parameter	real	xrbp	=   0.0
	  `P(info="temperature exponent of rbp               " );
	parameter	real	isrr	=   1.0		from(0.0:inf]
	  `P(info="ratio of is(reverse) to is(forward)       " );
	parameter	real	xisr	=   0.0
	  `P(info="temperature exponent for isrr             " );
	parameter	real	dear	=   0.0
	  `P(info="delta activation energy for isrr          " unit="V");
	parameter	real	eap	=   1.12
	  `P(info="activiation energy for isp                " unit="V");
	parameter	real	vbbe	=   0.0
	  `P(info="b-e breakdown voltage                     " unit="V");
	parameter	real	nbbe	=   1.0		from(0.0:inf]
	  `P(info="b-e breakdown emission coefficient        " );
	parameter	real	ibbe	=   1.0e-6
	  `P(info="b-e breakdown current                     " unit="A");
	parameter	real	tvbbe1	=   0.0
	  `P(info="linear temperature coefficient of vbbe    " unit="/C");
	parameter	real	tvbbe2	=   0.0
	  `P(info="quadratic temperature coefficient of vbbe " unit="/C^2");
	parameter	real	tnbbe	=   0.0
	  `P(info="temperature coefficient of nbbe           " );
	parameter	real	ebbe	=   0.0
	  `P(info="calculated exp(-vbbe/(nbbe*Vtv))          " );
	parameter	real	dtemp	=   0.0
	  `P(info="local temperature rise                    " unit="C");
	parameter	real	vers	=   1.2
	  `P(info="version number                            " );
	parameter	real	vrev	=   1.0
	  `P(info="revision number                           " );

	parameter	real	xrb	=   0.0
	  `P(info="temp exponent of rbx/i, xrbx/i not given  " );
	parameter	real	xrc	=   0.0
	  `P(info="temp exp rcx/i&rbp, xrcx/i&xrbp not given " );
	parameter	real	npn	=   0.0
	  `P(spice:isflag="yes" info="model type flag  for npn                  " );
	parameter	real	pnp	=   0.0
	  `P(info="model type flag  for pnp                  " );
	parameter	real	m	=   1.0		from(0.0:inf]
	  `P(spice:isflag="yes" type="instance" info="multiplicity scale factor                 " );
	parameter	real	mag	=   1.0		from(0.0:inf]
	  `P(type="instance" info="multiplicity scale factor                 " );
`ifdef insideADMS
`else
	parameter	real	gmin	=   1.0e-12	from[0.0:inf];
	parameter	real	pnjmaxi	=   1.0		from[0.0:inf];
`endif
	parameter	real	tmin	=  -1.0e+02	from[-250:27]
	  `P(info="min temperature                           " unit="K");
	parameter	real	tmax	=   5.0e+02	from[27:1000]
	  `P(info="max temperature                           " unit="K");
	parameter	real	shrink	=   0.0		from[0.0:100.0]
	  `P(info="shrink factor                             " );
	parameter	real	shrink2	=   0.0		from[0.0:100.0]
	  `P(info="shrink2 factor                            " );
	parameter	integer	off	=   1.0		from[0:1]
	  `P(spice:isflag="yes" type="instance" info="Device initially off?                     " );
	parameter	real	area	=   1.0		from(0.0:inf]
	  `P(type="instance" info="area factor                               " );
	parameter	real	maxexp	=   1.0e22	from[1.0e10:inf]
	  `P(info="maximum allowed value of exponential      " );

	real	is_t     `P(ask="yes" info="transport saturation current at temp      " unit="A");
	real	isrr_t   `P(ask="yes" info="ratio of is(rev) to is(fwd) at temperature" );
	real	ikf_t    `P(ask="yes" info="fwd knee current (zero=infinite) at temp  " unit="A");
	real	ibei_t   `P(ask="yes" info="ideal b-e saturation current at temp      " unit="A");
	real	ibci_t   `P(ask="yes" info="ideal b-c saturation current at temp      " unit="A");
	real	isp_t    `P(ask="yes" info="parasitic transport satn current at temp  " unit="A");
	real	iben_t   `P(ask="yes" info="non-ideal b-e saturation current at temp  " unit="A");
	real	ibcn_t   `P(ask="yes" info="non-ideal b-c saturation current at temp  " unit="A");
	real	ibeip_t  `P(ask="yes" info="ideal parasitic b-e satn current at temp  " unit="A");
	real	ibenp_t  `P(ask="yes" info="non-ideal para b-e satn current at temp   " unit="A");
	real	ibcip_t  `P(ask="yes" info="ideal parasitic b-c satn current at temp  " unit="A");
	real	ibcnp_t  `P(ask="yes" info="non-ideal para b-c satn current at temp   " unit="A");
	real	rcx_t    `P(ask="yes" info="extrinsic collector resistance at temp    " unit="Ohm");
	real	rci_t    `P(ask="yes" info="intrinsic collector resistance at temp    " unit="Ohm");
	real	rbx_t    `P(ask="yes" info="extrinsic base resistance at temperature  " unit="Ohm");
	real	rbi_t    `P(ask="yes" info="intrinsic base resistance at temperature  " unit="Ohm");
	real	re_t     `P(ask="yes" info="extrinsic emitter resistance at temp      " unit="Ohm");
	real	rs_t     `P(ask="yes" info="extrinsic substrate resistance at temp    " unit="Ohm");
	real	rbp_t    `P(ask="yes" info="parasitic base resistance at temperature  " );
	real	pe_t     `P(ask="yes" info="b-e built-in potential at temperature     " unit="V");
	real	pc_t     `P(ask="yes" info="b-c built-in potential at temperature     " unit="V");
	real	ps_t     `P(ask="yes" info="c-s built-in potential at temperature     " unit="V");
	real	cje_t    `P(ask="yes" info="zero-bias b-e depletion cap at temp       " unit="F");
	real	cjc_t    `P(ask="yes" info="zero-bias intrinsic b-c depl cap at temp  " unit="F");
	real	cjep_t   `P(ask="yes" info="zero-bias extrinsic b-c depl cap at temp  " unit="F");
	real	cjcp_t   `P(ask="yes" info="zero-bias c-s depletion cap at temperature" unit="F");
	real	nf_t     `P(ask="yes" info="forward emission coefficient at temp      " );
	real	nr_t     `P(ask="yes" info="reverse emission coefficient at temp      " );
	real	avc2_t   `P(ask="yes" info="b-c weak avalanche parameter 2 at temp    " );
	real	vbbe_t   `P(ask="yes" info="b-e breakdown voltage at temperature      " unit="V");
	real	nbbe_t   `P(ask="yes" info="b-e breakdown emission coeff at temp      " );
	real	gamm_t   `P(ask="yes" info="epi doping parameter at temperature       " );
	real	vo_t     `P(ask="yes" info="epi drift saturation voltage at temp      " unit="V");
	real	ebbe_t;
	real	Tdev     `P(ask="yes" spice:name="tdev" info="device temperature (K)                    " );
	real	Tini, rT, dT, Ivef, Iver, Iikf, Iikr, Iikp;
	real	Ivo, Ihrcf, Ivtf, Iitf, sltf, psiio, psiin, dv0;
	real	dvh, pwq, qlo, qhi, mv0, vl0, q0, dv;
	real	mv, vl, vminm, vn0, vnl0, qlo0, vn, vnl;
	real	sel, crt, cmx, cl, ql, Gcx, Gci, Gbx;
	real	Gbi, Ge, Gbp, Gs, Gth, maxvIfi, maxvIri;
	real	maxvIp, maxvIbbe, maxvIbei, maxvIben, maxvIbci, maxvIbcn;
	real	maxvIbeip, maxvIbenp, maxvIbcip, maxvIbcnp;
	real	Vtv, Ifi, Iri, Itzf;
	real	Itxf `P(ask="yes" spice:name="itzf" info="forward transport current                 " unit="A");
	real	Itzr `P(ask="yes" spice:name="itzr" info="reverse transport current                 " unit="A");
	real	q1z, q1, q2, qb, Ifp, Irp;
	real	Iccp `P(ask="yes" spice:name="iccp" info="parasitic transport current               " unit="A");
	real	q2p, qbp;
	real	Ibe  `P(ask="yes" spice:name="ibe" info="intrinsic b-e current                     " unit="A");
	real	Ibex `P(ask="yes" spice:name="ibex" info="extrinsic b-e current                     " unit="A");
	real	Ibcj;
	real	Ibc  `P(ask="yes" spice:name="ibc" info="intrinsic b-c current                     " unit="A");
	real	Ibep `P(ask="yes" spice:name="ibep" info="extrinsic b-c current                     " unit="A");
	real	Ibcp `P(ask="yes" spice:name="ibcp" info="parasitic b-c current                     " unit="A");
	real	Igc `P(ask="yes" spice:name="igc" info="c-b avalanche current                     " unit="A");
	real	avalf;
	real	Ircx `P(ask="yes" spice:name="ircx" info="current in rcx element                    " unit="A");
	real	Irci `P(ask="yes" spice:name="irci" info="current in rci element                    " unit="A");
	real	Irbx `P(ask="yes" spice:name="irbx" info="current in rbx element                    " unit="A");
	real	Irbi `P(ask="yes" spice:name="irbi" info="current in rbi element                    " unit="A");
	real	Ire  `P(ask="yes" spice:name="ire" info="current in re  element                    " unit="A");
	real	Irbp `P(ask="yes" spice:name="irbp" info="current in rbp element                    " unit="A");
	real	Irs  `P(ask="yes" spice:name="irs" info="current in rs  element                    " unit="A");
	real	Kbci, Kbcx, rKp1, Iohm, derf, arg, expl;
	real	expi, expn, expx, afac, mac1, VmaxExp, qdbe, qdbex;
	real	qdbc, qdbep, qdbcp;

	// conductances (for backward compatibility with level 2)
	real	ibc_vbci `P(ask="yes" info="Derivative of ibc with respect to vbci" unit="S");
	real	 ibcp_vbcp `P(ask="yes" info="Derivative of ibcp with respect to vbcp" unit="S");
	real	 ibe_vbei `P(ask="yes" info="Derivative of ibe with respect to vbei" unit="S");
	real	 ibep_vbep `P(ask="yes" info="Derivative of ibep with respect to vbep" unit="S");
	real	 ibex_vbex `P(ask="yes" info="Derivative of ibex with respect to Vbex" unit="S");
	real	 iccp_vbci `P(ask="yes" info="Derivative of iccp with respect to vbci" unit="S");
	real	iccp_vbcp `P(ask="yes" info="Derivative of iccp with respect to vbcp" unit="S");
	real	 iccp_vbep `P(ask="yes" info="Derivative of iccp with respect to vbep" unit="S");
	real	 igc_vbci `P(ask="yes" info="Derivative of igc with respect to vbci" unit="S");
	real	 igc_vbei `P(ask="yes" info="Derivative of igc with respect to vbei" unit="S");
	real	 irbi_vbci `P(ask="yes" info="Derivative of irbi with respect to vbci" unit="S");
	real	 irbi_vbei `P(ask="yes" info="Derivative of irbi with respect to vbei" unit="S");
	real	irbp_vbci `P(ask="yes" info="Derivative of irbp with respect to vbci" unit="S");
	real	 irci_vbci `P(ask="yes" info="Derivative of irci with respect to vbci" unit="S");
	real	 itzf_vbci `P(ask="yes" info="Derivative of itzf with respect to vbci" unit="S");
	real	 itzf_vbei `P(ask="yes" info="Derivative of itzf with respect to vbei" unit="S");
	real	 itzr_vbci `P(ask="yes" info="Derivative of itzr with respect to vbci" unit="S");
	real	 itzr_vbei `P(ask="yes" info="Derivative of itzr with respect to vbei" unit="S");
	// capacitances (for backward compatibility with level 2)
	real	qbc_vbci `P(ask="yes" info="Derivative of qbc with respect to vbci    " unit="C/V");
	real	 qbco_vbc `P(ask="yes" info="Derivative of qbco with respect to vbc    " unit="C/V");
	real	 qbcp_vbcp `P(ask="yes" info="Derivative of qbcp with respect to vbcp   " unit="C/V");
	real	 qbcx_vbcx `P(ask="yes" info="Derivative of qbcx with respect to vbcx   " unit="C/V");
	real	 qbe_vbci `P(ask="yes" info="Derivative of qbe with respect to vbci    " unit="C/V");
	real	qbe_vbei `P(ask="yes" info="Derivative of qbe with respect to vbei    " unit="C/V");
	real	 qbeo_vbe `P(ask="yes" info="Derivative of qbeo with respect to vbe    " unit="C/V");
	real	 qbep_vbci `P(ask="yes" info="Derivative of qbep with respect to vbci   " unit="C/V");
	real	 qbep_vbep `P(ask="yes" info="Derivative of qbep with respect to vbep   " unit="C/V");
	real	 qbex_vbex `P(ask="yes" info="Derivative of qbex with respect to vbex   " unit="C/V");
	// power current (for backward compatibility with level 2)
	real	p `P(ask="yes" info="thermal power generation                  " unit="W");
	// temperature (for backward compatibility with level 2)
	real	trise `P(ask="yes" info="local temperature rise                    " unit="C");
	real	ft_int;

	real	sgIf, rIf, mIf, tff;
	real	Qbe `P(ask="yes" spice:name="qbe" info="intrinsic b-e charge                      " unit="C");
	real	Qbex `P(ask="yes" spice:name="qbex" info="extrinsic (side-wall) b-e charge          " unit="C");
	real	Qbc `P(ask="yes" spice:name="qbc" info="intrinsic b-c charge                      " unit="C");
	real	Qbcx `P(ask="yes" spice:name="qbcx" info="extrinsic b-c charge                      " unit="C");
	real	Qbep `P(ask="yes" spice:name="qbep" info="parasitic b-e charge                      " unit="C");
	real	Qbcp `P(ask="yes" spice:name="qbcp" info="parasitic b-c charge                      " unit="C");
	real	Qbeo `P(ask="yes" spice:name="qbeo" info="extrinsic b-e overlap charge              " unit="C");
	real	Qbco `P(ask="yes" spice:name="qbco" info="extrinsic b-c overlap charge              " unit="C");
	real	Ixf1, Ixf2, Qxf1, Qxf2, Ith, Irth, Qcth;
	real	Vbei, Vbci, Vbex, Vbep, Vbcp, Vbcx, Vxf1, Vxf2;
	real	Bvbe, Vrth;
	real	Vbe `P(ask="yes" spice:name="vbe" info="external b-e voltage                      " unit="V");
	real	Vbc `P(ask="yes" spice:name="vbc" info="external b-c voltage                      " unit="V");
	real	Vrcx, Vrci, Vrbx, Vrbi, Vre, Vrbp, Vrs, Vcei, Vcep;
	real	VBICtype `P(spice:name="type" info="Device type from npn or pnp flags         " unit="no");
	real	meff `P(ask="yes" info="multiplicity scale factor (incl subckt)   " );
	real	areaef `P(ask="yes" info="Effective area factor                     " unit="m^2");

	real	shrinklFactor `P(ask="yes" spice:name="shrinkl" info="lineal shrink factor                      " );
	real	shrinkaFactor `P(ask="yes" spice:name="shrinka" info="areal shrink factor                       " );
	real	instanceShrinklFactor `P(ask="yes" spice:name="shrinkl" info="lineal shrink factor                      " );
	real	instanceShrinkaFactor `P(ask="yes" spice:name="shrinka" info="areal shrink factor                       " );
	real	ic `P(ask="yes" info="current into collector                    " unit="A");
	real	ib `P(ask="yes" info="current into base                         " unit="A");
	real	ie `P(ask="yes" info="current into emittor                      " unit="A");
	real	isub `P(ask="yes" spice:name="is" info="current into substrate                    " unit="A");
	real	powerT `P(ask="yes" info="thermal power generation                  " unit="W");
	real	powerD `P(ask="yes" info="thermal power dissipation in rth          " unit="W");
	real	Vce `P(ask="yes" spice:name="vce" info="external c-e voltage                      " unit="V");
	real	Ircx_Vrcx `P(ask="yes" spice:name="ircx_vrcx" info="rcx conductance                           " unit="S");
	real	Irci_Vrci `P(ask="yes" spice:name="irci_vrci" info="rci effective conductance                 " unit="S");
	real	Irbx_Vrbx `P(ask="yes" spice:name="irbx_vrbx" info="rbx conductance                           " unit="S");
	real	Irbi_Vrbi `P(ask="yes" spice:name="irbi_vrbi" info="rbi effective conductance                 " unit="S");
	real	Irbp_Vrbp `P(ask="yes" spice:name="irbp_vrbp" info="rbp effective conductance                 " unit="S");
	real	Ire_Vre `P(ask="yes" spice:name="ire_vre" info="re conductance                            " unit="S");
	real	Irs_Vrs `P(ask="yes" spice:name="irs_vrs" info="rs conductance                            " unit="S");
	real	gm `P(ask="yes" info="intrinsic transconductance                " unit="S");
	real	go `P(ask="yes" info="intrinsic output conductance              " unit="S");
	real	gpi `P(ask="yes" info="intrinsic b-e conductance                 " unit="S");
	real	gpix `P(ask="yes" info="extrinsic b-e conductance                 " unit="S");
	real	gmu `P(ask="yes" info="intrinsic b-c conductance                 " unit="S");
	real	gmux `P(ask="yes" info="extrinsic b-c conductance                 " unit="S");
	real	cpi `P(ask="yes" info="intrinsic b-e capacitance d(Qbe)/d(Vbei)  " unit="F");
	real	cpix `P(ask="yes" info="extrinsic b-e capacitance d(Qbex)/d(Vbex) " unit="F");
	real	cmu `P(ask="yes" info="intrinsic b-c capacitance d(Qbc)/d(Vbci)  " unit="F");
	real	cmux `P(ask="yes" info="extrinsic b-c capacitance d(Qbcx)/d(Vbcx) " unit="F");
	real	cpip `P(ask="yes" info="parasitic b-e capacitance d(Qbep)/d(Vbep) " unit="F");
	real	ccs `P(ask="yes" info="c-s capacitance d(Qbcp)/d(Vbcp)           " unit="F");
	real	cbeoXX `P(ask="yes" spice:name="cbeo" info="b-e overlap capacitance d(Qbeo)/d(Vbe)    " unit="F");
	real	cbcoXX `P(ask="yes" spice:name="cbco" info="b-c overlap capacitance d(Qbco)/d(Vbc)    " unit="F");
	real	rci_mod `P(ask="yes" spice:name="rcimod" info="rci modulatn fctr, rci_atBias/rci_zeroBias" );
	real	rci_eff `P(ask="yes" info="effective (small-signal) intr collctr res " unit="Ohm");
	real	rbi_eff `P(ask="yes" info="modulated intrinsic base resistance       " unit="Ohm");
	real	rbp_eff `P(ask="yes" info="modulated parasitic base resistance       " unit="Ohm");
	real	tauTh `P(ask="yes" info="thermal time constant                     " unit="s");

	parameter	integer level=3    `P(info="level number                              " );
	parameter	integer libset = 0 `P(info="model parameter set ID number             " );
	string	ft            `P(ask="yes" info="Cutoff frequency                          " unit="Hz");
	string	VBICtypeName  `P(ask="yes" spice:name="type" info="Device type from npn or pnp flags         " unit="no");
	string	INSTANCEmodel `P(ask="yes" spice:name="model" info="model name                                " );
	string	INSTANCEtype  `P(ask="yes" spice:name="type" info="bjt transistor type                       " );
	real	INSTANCElevel `P(ask="yes" spice:name="level" info="level number                              " );

	real	inoise        `P(ask="yes" spice:name="inoise" info="Total noise"
	  mica:ask="return(TOTnoiseAsk(&value->rValue, ASK_IN_NOISE, instance->VBICSELFTtotNoiseIndex));"
	);
	real	inoise_itzf   `P(ask="yes" spice:name="inoise.itzf" info="itzf shot thermal noise"
	  mica:ask="return(SHOTnoiseAsk(&value->rValue, ASK_IN_NOISE, instance->VBICSELFTshotSrcIndex, ITZF_NOISE));"
	);
	real	inoise_ibe    `P(ask="yes" spice:name="inoise.ibe" info="ibe shot thermal noise"
	  mica:ask="return(SHOTnoiseAsk(&value->rValue, ASK_IN_NOISE, instance->VBICSELFTshotSrcIndex, IBE_NOISE));"
	);
	real	inoise_ibefn  `P(ask="yes" spice:name="inoise.ibefn" info="ibe flicker thermal noise"
	  mica:ask="return(FLICKERnoiseAsk(&value->rValue, ASK_IN_NOISE, instance->VBICSELFTflickerSrcIndex, IBE_FLICKER));"
	);
	real	inoise_rcx    `P(ask="yes" spice:name="inoise.rcx" info="rcx collector thermal noise"
	  mica:ask="return(THERMALnoiseAsk(&value->rValue, ASK_IN_NOISE, instance->VBICSELFTthermalSrcIndex, RCX_NOISE));"
	);
	real	inoise_rci    `P(ask="yes" spice:name="inoise.rci" info="rci collector thermal noise"
	  mica:ask="return(THERMALnoiseAsk(&value->rValue, ASK_IN_NOISE, instance->VBICSELFTthermalSrcIndex, RCI_NOISE));"
	);
	real	inoise_rbx    `P(ask="yes" spice:name="inoise.rbx" info="rbx base thermal noise"
	  mica:ask="return(THERMALnoiseAsk(&value->rValue, ASK_IN_NOISE, instance->VBICSELFTthermalSrcIndex, RBX_NOISE));"
	);
	real	inoise_rbi    `P(ask="yes" spice:name="inoise.rbi" info="rbi base thermal noise"
	  mica:ask="return(THERMALnoiseAsk(&value->rValue, ASK_IN_NOISE, instance->VBICSELFTthermalSrcIndex, RBI_NOISE));"
	);
	real	inoise_re     `P(ask="yes" spice:name="inoise.re" info="re emitter thermal noise"
	  mica:ask="return(THERMALnoiseAsk(&value->rValue, ASK_IN_NOISE, instance->VBICSELFTthermalSrcIndex, RE_NOISE));"
	);
	real	inoise_rs     `P(ask="yes" spice:name="inoise.rs" info="rs thermal noise"
	  mica:ask="return(THERMALnoiseAsk(&value->rValue, ASK_IN_NOISE, instance->VBICSELFTthermalSrcIndex, RS_NOISE));"
	);
	real	inoise_iccp   `P(ask="yes" spice:name="inoise.iccp" info="iccp shot noise"
	  mica:ask="return(SHOTnoiseAsk(&value->rValue, ASK_IN_NOISE, instance->VBICSELFTshotSrcIndex, ICCP_NOISE));"
	);
	real	inoise_ibep   `P(ask="yes" spice:name="inoise.ibep" info="ibep shot noise"
	  mica:ask="return(SHOTnoiseAsk(&value->rValue, ASK_IN_NOISE, instance->VBICSELFTshotSrcIndex, IBEP_NOISE));   "
	);
	real	inoise_ibepfn `P(ask="yes" spice:name="inoise.ibepfn" info="ibep flicker thermal noise"
	  mica:ask="return(FLICKERnoiseAsk(&value->rValue, ASK_IN_NOISE, instance->VBICSELFTflickerSrcIndex, IBEP_FLICKER));"
	);
	real	inoise_rbp    `P(ask="yes" spice:name="inoise.rbp" info="rbp thermal noise"
	  mica:ask="return(THERMALnoiseAsk(&value->rValue, ASK_IN_NOISE, instance->VBICSELFTthermalSrcIndex, RBP_NOISE));"
	);
	real	onoise        `P(ask="yes" spice:name="onoise" info="Total noise"
	  mica:ask="return(TOTnoiseAsk(&value->rValue, ASK_OUT_NOISE, instance->VBICSELFTtotNoiseIndex));"
	);
	real	onoise_itzf   `P(ask="yes" spice:name="onoise.itzf" info="itzf shot thermal noise"
	  mica:ask="return(SHOTnoiseAsk(&value->rValue, ASK_OUT_NOISE, instance->VBICSELFTshotSrcIndex, ITZF_NOISE));"
	);
	real	onoise_ibe    `P(ask="yes" spice:name="onoise.ibe" info="ibe shot thermal noise"
	  mica:ask="return(SHOTnoiseAsk(&value->rValue, ASK_OUT_NOISE, instance->VBICSELFTshotSrcIndex, IBE_NOISE)); "
	);
	real	onoise_ibefn  `P(ask="yes" spice:name="onoise.ibefn" info="ibe flicker thermal noise"
	  mica:ask="return(FLICKERnoiseAsk(&value->rValue, ASK_OUT_NOISE, instance->VBICSELFTflickerSrcIndex, IBE_FLICKER));"
	);
	real	onoise_rcx    `P(ask="yes" spice:name="onoise.rcx" info="rcx collector thermal noise"
	  mica:ask="return(THERMALnoiseAsk(&value->rValue, ASK_OUT_NOISE, instance->VBICSELFTthermalSrcIndex, RCX_NOISE));"
	);
	real	onoise_rci    `P(ask="yes" spice:name="onoise.rci" info="rci collector thermal noise"
	  mica:ask="return(THERMALnoiseAsk(&value->rValue, ASK_OUT_NOISE, instance->VBICSELFTthermalSrcIndex, RCI_NOISE));"
	);
	real	onoise_rbx    `P(ask="yes" spice:name="onoise.rbx" info="rbx base thermal noise"
	  mica:ask="return(THERMALnoiseAsk(&value->rValue, ASK_OUT_NOISE, instance->VBICSELFTthermalSrcIndex, RBX_NOISE));"
	);
	real	onoise_rbi    `P(ask="yes" spice:name="onoise.rbi" info="rbi base thermal noise"
	  mica:ask="return(THERMALnoiseAsk(&value->rValue, ASK_OUT_NOISE, instance->VBICSELFTthermalSrcIndex, RBI_NOISE));"
	);
	real	onoise_re     `P(ask="yes" spice:name="onoise.re" info="re emitter thermal noise"
	  mica:ask="return(THERMALnoiseAsk(&value->rValue, ASK_OUT_NOISE, instance->VBICSELFTthermalSrcIndex, RE_NOISE));"
	);
	real	onoise_rs     `P(ask="yes" spice:name="onoise.rs" info="rs thermal noise"
	  mica:ask="return(THERMALnoiseAsk(&value->rValue, ASK_OUT_NOISE, instance->VBICSELFTthermalSrcIndex, RS_NOISE));"
	);
	real	onoise_iccp   `P(ask="yes" spice:name="onoise.iccp" info="iccp shot noise"
	  mica:ask="return(SHOTnoiseAsk(&value->rValue, ASK_OUT_NOISE, instance->VBICSELFTshotSrcIndex, ICCP_NOISE));"
	);
	real	onoise_ibep   `P(ask="yes" spice:name="onoise.ibep" info="ibep shot noise"
	  mica:ask="return(SHOTnoiseAsk(&value->rValue, ASK_OUT_NOISE, instance->VBICSELFTshotSrcIndex, IBEP_NOISE)); "
	);
	real	onoise_ibepfn `P(ask="yes" spice:name="onoise.ibepfn" info="ibep flicker thermal noise"
	  mica:ask="return(FLICKERnoiseAsk(&value->rValue, ASK_OUT_NOISE, instance->VBICSELFTflickerSrcIndex, IBEP_FLICKER));"
	);
	real	onoise_rbp    `P(ask="yes" spice:name="onoise.rbp" info="rbp thermal noise"
	  mica:ask="return(THERMALnoiseAsk(&value->rValue, ASK_OUT_NOISE, instance->VBICSELFTthermalSrcIndex, RBP_NOISE));"
	);

	analog begin

	`INITIAL_MODEL
	begin	// begin loadmodel block: code independent of instance parameters or bias
		if (`PGIVEN(shrink)) begin
			shrinklFactor=1.0-shrink/100.0;
			if (`PGIVEN(shrink2))
				shrinkaFactor=1.0-shrink2/100.0;
			else
				shrinkaFactor=shrinklFactor*shrinklFactor;
		end else begin
			shrinklFactor=`SHRINKL;
			if (`PGIVEN(shrink2))
				shrinkaFactor=1.0-shrink2/100.0;
			else
				shrinkaFactor=`SHRINKA;
		end

		if (`PGIVEN(npn))
			VBICtype	=  `NPN;
		else if (`PGIVEN(pnp))
			VBICtype	=  `PNP;
		else
			VBICtype	=  `NPN;
`ifdef insideADMS
		VBICtypeName = (VBICtype==`NPN)?"Vertical npn":"Vertical pnp";
`endif
	end	// end loadmodel block

//
//	Temperature mappings
//

	`INITIAL_INSTANCE
	begin	// begin loadinstance block: code independent of bias
`ifdef insideADMS
		INSTANCEmodel	= $model;
		INSTANCElevel	= level;
		INSTANCEtype=VBICtypeName;
`endif
		instanceShrinklFactor=shrinklFactor;
		instanceShrinkaFactor=shrinkaFactor;
		if (`PGIVEN(shrink2))
			meff	=  m*mag*area*(1.0-shrink2/100.0);
		else if (`PGIVEN(shrink))
			meff	=  m*mag*area*(1.0-shrink/100.0)*(1.0-shrink/100.0);
		else
			meff	=  m*mag*area*`SHRINKA;
	end	// end loadinstance block
		Tini	=  `TABS+tnom;
		Vrth	=  V(b_rth);
		Tdev	=  $temperature+dtemp+Vrth;
		Tdev	=  Tdev-`TABS;
		if (Tdev<tmin)
			Tdev	=  tmin;
		if (Tdev>tmax)
			Tdev	=  tmax;
		Tdev	=  Tdev+`TABS;
		Vtv	=  `KB*Tdev/`QQ;
		rT	=  Tdev/Tini;
		dT	=  Tdev-Tini;
		ikf_t	=  ikf*pow(rT,xikf);
		if (`PGIVEN(xrcx))
			rcx_t	=  rcx*pow(rT,xrcx);
		else
			rcx_t	=  rcx*pow(rT,xrc);
		if (`PGIVEN(xrci))
			rci_t	=  rci*pow(rT,xrci);
		else
			rci_t	=  rci*pow(rT,xrc);
		if (`PGIVEN(xrbx))
			rbx_t	=  rbx*pow(rT,xrbx);
		else
			rbx_t	=  rbx*pow(rT,xrb);
		if (`PGIVEN(xrbi))
			rbi_t	=  rbi*pow(rT,xrbi);
		else
			rbi_t	=  rbi*pow(rT,xrb);
		re_t	=  re*pow(rT,xre);
		rs_t	=  rs*pow(rT,xrs);
		if (`PGIVEN(xrbp))
			rbp_t	=  rbp*pow(rT,xrbp);
		else
			rbp_t	=  rbp*pow(rT,xrc);
		is_t	=  is*pow((pow(rT,xis)*exp(-ea*(1.0-rT)/Vtv)),(1.0/nf));
		isrr_t	=  isrr*pow((pow(rT,xisr)*exp(-dear*(1.0-rT)/Vtv)),(1.0/nr));
		isp_t	=  isp*pow((pow(rT,xis)*exp(-eap*(1.0-rT)/Vtv)),(1.0/nfp));
		ibei_t	=  ibei*pow((pow(rT,xii)*exp(-eaie*(1.0-rT)/Vtv)),(1.0/nei));
		iben_t	=  iben*pow((pow(rT,xin)*exp(-eane*(1.0-rT)/Vtv)),(1.0/nen));
		ibci_t	=  ibci*pow((pow(rT,xii)*exp(-eaic*(1.0-rT)/Vtv)),(1.0/nci));
		ibcn_t	=  ibcn*pow((pow(rT,xin)*exp(-eanc*(1.0-rT)/Vtv)),(1.0/ncn));
		ibeip_t	=  ibeip*pow((pow(rT,xii)*exp(-eaic*(1.0-rT)/Vtv)),(1.0/nci));
		ibenp_t	=  ibenp*pow((pow(rT,xin)*exp(-eanc*(1.0-rT)/Vtv)),(1.0/ncn));
		ibcip_t	=  ibcip*pow((pow(rT,xii)*exp(-eais*(1.0-rT)/Vtv)),(1.0/ncip));
		ibcnp_t	=  ibcnp*pow((pow(rT,xin)*exp(-eans*(1.0-rT)/Vtv)),(1.0/ncnp));
		nf_t	=  nf*(1.0+dT*tnf);
		nr_t	=  nr*(1.0+dT*tnf);
		avc2_t	=  avc2*(1.0+dT*tavc);
		vbbe_t	=  vbbe*(1.0+dT*(tvbbe1+dT*tvbbe2));
		nbbe_t	=  nbbe*(1.0+dT*tnbbe);
		`psibi(pe_t,pe,eaie,Vtv,rT);
		`psibi(pc_t,pc,eaic,Vtv,rT);
		`psibi(ps_t,ps,eais,Vtv,rT);
		cje_t	=  cje*pow(pe/pe_t,me);
		cjc_t	=  cjc*pow(pc/pc_t,mc);
		cjep_t	=  cjep*pow(pc/pc_t,mc);
		cjcp_t	=  cjcp*pow(ps/ps_t,ms);
		gamm_t	=  gamm*pow(rT,xis)*exp(-ea*(1.0-rT)/Vtv);
		vo_t	=  vo*pow(rT,xvo);
		ebbe_t	=  exp(-vbbe_t/(nbbe_t*Vtv));
		if (ibbe>0.0)
			maxvIbbe=  nbbe_t*Vtv*ln(ebbe_t+`PNJMAXI/ibbe);
		else
			maxvIbbe=  0.0;
		if (is_t>0.0)
			maxvIfi	=  nf_t*Vtv*ln(1.0+`PNJMAXI/is_t);
		else
			maxvIfi	=  0.0;
		if (is_t>0.0&&isrr_t>0.0)
			maxvIri	=  nr_t*Vtv*ln(1.0+`PNJMAXI/(is_t*isrr_t));
		else
			maxvIri	=  0.0;
		if (isp_t>0.0)
			maxvIp	=  nfp*Vtv*ln(1.0+`PNJMAXI/isp_t);
		else
			maxvIp	=  0.0;
		if (ibei_t>0.0)
			maxvIbei=  nei*Vtv*ln(1.0+`PNJMAXI/ibei_t);
		else
			maxvIbei=  0.0;
		if (iben_t>0.0)
			maxvIben=  nen*Vtv*ln(1.0+`PNJMAXI/iben_t);
		else
			maxvIben=  0.0;
		if (ibci_t>0.0)
			maxvIbci=  nci*Vtv*ln(1.0+`PNJMAXI/ibci_t);
		else
			maxvIbci=  0.0;
		if (ibcn_t>0.0)
			maxvIbcn=  ncn*Vtv*ln(1.0+`PNJMAXI/ibcn_t);
		else
			maxvIbcn=  0.0;
		if (ibeip_t>0.0)
			maxvIbeip=  nci*Vtv*ln(1.0+`PNJMAXI/ibeip_t);
		else
			maxvIbeip=  0.0;
		if (ibenp_t>0.0)
			maxvIbenp=  ncn*Vtv*ln(1.0+`PNJMAXI/ibenp_t);
		else
			maxvIbenp=  0.0;
		if (ibcip_t>0.0)
			maxvIbcip=  ncip*Vtv*ln(1.0+`PNJMAXI/ibcip_t);
		else
			maxvIbcip=  0.0;
		if (ibcnp_t>0.0)
			maxvIbcnp=  ncnp*Vtv*ln(1.0+`PNJMAXI/ibcnp_t);
		else
			maxvIbcnp=  0.0;
		VmaxExp	=  ln(maxexp);

//
//	Parameter mappings
//

		Gcx	=  rcx_t>1.0e-3 ? 1.0/rcx_t  : 1.0e3;
		Gci	=  rci_t>1.0e-3 ? 1.0/rci_t  : 1.0e3;
		Gbx	=  rbx_t>1.0e-3 ? 1.0/rbx_t  : 1.0e3;
		Gbi	=  rbi_t>1.0e-3 ? 1.0/rbi_t  : 1.0e3;
		Ge	=  re_t >1.0e-3 ? 1.0/re_t   : 1.0e3;
		Gbp	=  rbp_t>1.0e-3 ? 1.0/rbp_t  : 1.0e3;
		Gs	=  rs_t >1.0e-3 ? 1.0/rs_t   : 1.0e3;
		Gth	=  rth  >1.0e-3 ? 1.0/rth    : 1.0e3;
		Ivef	=  vef  >0.0    ? 1.0/vef    : 0.0;
		Iver	=  ver  >0.0    ? 1.0/ver    : 0.0;
		Iikf	=  ikf  >0.0    ? 1.0/ikf_t  : 0.0;
		Iikr	=  ikr  >0.0    ? 1.0/ikr    : 0.0;
		Iikp	=  ikp  >0.0    ? 1.0/ikp    : 0.0;
		Ivo	=  vo   >0.0    ? 1.0/vo_t   : 0.0;
		Ihrcf	=  hrcf >0.0    ? 1.0/hrcf   : 0.0;
		Ivtf	=  vtf  >0.0    ? 1.0/vtf    : 0.0;
		Iitf	=  itf  >0.0    ? 1.0/itf    : 0.0;
		sltf	=  itf  >0.0    ? 0.0        : 1.0;

//
//	Branch voltages
//

		Vbei	=  VBICtype*V(b_bei);
		Vbex	=  VBICtype*V(b_bex);
		Vbci	=  VBICtype*V(b_bci);
		Vbcx	=  VBICtype*V(b_bcx);
		Vcei	=  VBICtype*V(b_cei);
		Vbep	=  VBICtype*V(b_bep);
		Vbcp	=  VBICtype*V(b_bcp);
		Vcep	=  VBICtype*V(b_cep);
		Vxf1	=           V(b_xf1);
		Vxf2	=           V(b_xf2);
		Vbe	=           V(b_be);
		Vbc	=           V(b_bc);
		Vrcx	=           V(b_rcx);
		Vrci	=  VBICtype*V(b_rci);
		Vrbx	=           V(b_rbx);
		Vrbi	=           V(b_rbi);
		Vre	=           V(b_re);
		Vrbp	=           V(b_rbp);
		Vrs	=           V(b_rs);

//
//	Electrical branch constituent relations
//

//
//	Depletion charges
//

		`qj(qdbe,Vbei,pe_t,me,fc,aje);
		`qj(qdbex,Vbex,pe_t,me,fc,aje);
		`qjrt(qdbc,Vbci,pc_t,mc,fc,ajc,vrt,art);
		`qjrt(qdbep,Vbep,pc_t,mc,fc,ajc,vrt,art);
		if (cjcp>0.0) begin
			`qj(qdbcp,Vbcp,ps_t,ms,fc,ajs);
		end else
			qdbcp	=  0.0;

//
//	Transport current of main transistor
//

		afac	=  1.0/(nf_t*Vtv);
		`expLinA(expi,Vbei,maxvIfi,afac);
		Ifi	=  is_t*(expi-1.0);
		afac	=  1.0/(nr_t*Vtv);
		`expLinA(expi,Vbci,maxvIri,afac);
		Iri	=  is_t*isrr_t*(expi-1.0);
		q1z	=  1.0+qdbe*Iver+qdbc*Ivef;
		q1	=  0.5*(sqrt((q1z-1.0e-4)*(q1z-1.0e-4)+1.0e-8)+q1z-1.0e-4)+1.0e-4;
		q2	=  Ifi*Iikf+Iri*Iikr;
		if (qbm<0.5) begin
			arg	=  pow(q1,1.0/nkf)+4.0*q2;
			if (arg>1.0e-8)
				qb	=  0.5*(q1+pow(arg,nkf));
			else
				qb	=  0.5*(q1+pow(1.0e-8,nkf));
		end else begin
			arg	=  1.0+4.0*q2;
			if (arg>1.0e-8)
				qb	=  0.5*q1*(1.0+pow(arg,nkf));
			else
				qb	=  0.5*q1*(1.0+pow(1.0e-8,nkf));
		end
		Itzr	=  Iri/qb;
		Itzf	=  Ifi/qb;
		Itxf	=  Vxf2;

//
//	Transport current of parasitic transistor
//

		if (isp>0.0) begin
			afac	=  1.0/(nfp*Vtv);
			`expLinA(expi,Vbep,maxvIp,afac);
			`expLinA(expx,Vbci,maxvIp,afac);
			Ifp	=  isp_t*(wsp*expi+(1.0-wsp)*expx-1.0);
			q2p	=  Ifp*Iikp;
			arg	=  1.0+4.0*q2p;
			if (arg>1.0e-8)
				qbp	=  0.5*(1.0+sqrt(arg));
			else
				qbp	=  0.5*(1.0+sqrt(1.0e-8));
			`expLinA(expi,Vbcp,maxvIp,afac);
			Irp	=  isp_t*(expi-1.0);
			Iccp	=  (Ifp-Irp)/qbp;
		end else begin
			Ifp	=  0.0;
			qbp	=  1.0;
			Iccp	=  0.0;
		end

//
//	Diode-like currents for main transistor (includes
//	exponential-like breakdown for base-emitter)
//	and base-emitter of parasitic transistor
//

		if (wbe==1.0) begin
			afac	=  1.0/(nei*Vtv);
			`expLinA(expi,Vbei,maxvIbei,afac);
			afac	=  1.0/(nen*Vtv);
			`expLinA(expn,Vbei,maxvIben,afac);
			if (vbbe>0.0) begin
				Bvbe	=  -vbbe_t-Vbei;
				afac	=  1.0/(nbbe_t*Vtv);
				`expLinA(expx,Bvbe,maxvIbbe,afac);
				Ibe	=  ibei_t*(expi-1.0)+iben_t*(expn-1.0)-ibbe*(expx-ebbe_t);
			end else
				Ibe	=  ibei_t*(expi-1.0)+iben_t*(expn-1.0);
			Ibex	=  0.0;
		end else if (wbe==0.0) begin
			Ibe	=  0.0;
			afac	=  1.0/(nei*Vtv);
			`expLinA(expi,Vbex,maxvIbei,afac);
			afac	=  1.0/(nen*Vtv);
			`expLinA(expn,Vbex,maxvIben,afac);
			if (vbbe>0.0) begin
				Bvbe	=  -vbbe_t-Vbei;
				afac	=  1.0/(nbbe_t*Vtv);
				`expLinA(expx,Bvbe,maxvIbbe,afac);
				Ibex	=  ibei_t*(expi-1.0)+iben_t*(expn-1.0)-ibbe*(expx-ebbe_t);
			end else
				Ibex	=  ibei_t*(expi-1.0)+iben_t*(expn-1.0);
		end else begin
			afac	=  1.0/(nei*Vtv);
			`expLinA(expi,Vbei,maxvIbei,afac);
			afac	=  1.0/(nen*Vtv);
			`expLinA(expn,Vbei,maxvIben,afac);
			if (vbbe>0.0) begin
				Bvbe	=  -vbbe_t-Vbei;
				afac	=  1.0/(nbbe_t*Vtv);
				`expLinA(expx,Bvbe,maxvIbbe,afac);
				Ibe	=  wbe*(ibei_t*(expi-1.0)+iben_t*(expn-1.0)-ibbe*(expx-ebbe_t));
			end else
				Ibe	=  wbe*(ibei_t*(expi-1.0)+iben_t*(expn-1.0));
			afac	=  1.0/(nei*Vtv);
			`expLinA(expi,Vbex,maxvIbei,afac);
			afac	=  1.0/(nen*Vtv);
			`expLinA(expn,Vbex,maxvIben,afac);
			if (vbbe>0.0) begin
				Bvbe	=  -vbbe_t-Vbei;
				afac	=  1.0/(nbbe_t*Vtv);
				`expLinA(expx,Bvbe,maxvIbbe,afac);
				Ibex	=  (1.0-wbe)*(ibei_t*(expi-1.0)+iben_t*(expn-1.0)-ibbe*(expx-ebbe_t));
			end else
				Ibex	=  (1.0-wbe)*(ibei_t*(expi-1.0)+iben_t*(expn-1.0));
		end
		afac	=  1.0/(nci*Vtv);
		`expLinA(expi,Vbci,maxvIbci,afac);
		afac	=  1.0/(ncn*Vtv);
		`expLinA(expn,Vbci,maxvIbcn,afac);
		Ibcj	=  ibci_t*(expi-1.0)+ibcn_t*(expn-1.0);
		if ((ibeip>0.0)||(ibenp>0.0)) begin
			afac	=  1.0/(nci*Vtv);
			`expLinA(expi,Vbep,maxvIbeip,afac);
			afac	=  1.0/(ncn*Vtv);
			`expLinA(expn,Vbep,maxvIbenp,afac);
			Ibep	=  ibeip_t*(expi-1.0)+ibenp_t*(expn-1.0);
		end else
			Ibep	=  0.0;

//
//	Avalanche current
//

		if (avc1>0.0) begin
			`avalm(avalf,Vbci,pc_t,mc,avc1,avc2_t,VmaxExp);
			Igc	=  (Itxf-Itzr-Ibcj)*avalf;
		end else
			Igc	=  0.0;
		Ibc	=  Ibcj-Igc;

//
//	Base pushout charge factors
//

		arg	=  Vbci/Vtv;
		`expLin(expi,arg,VmaxExp);
		arg	=  Vbcx/Vtv;
		`expLin(expx,arg,VmaxExp);
		Kbci	=  sqrt(1.0+gamm_t*expi);
		Kbcx	=  sqrt(1.0+gamm_t*expx);

//
//	Diode-like current for base-collector of parasitic transistor
//

		if ((ibcip>0.0)||(ibcnp>0.0)) begin
			afac	=  1.0/(ncip*Vtv);
			`expLinA(expi,Vbcp,maxvIbcip,afac);
			afac	=  1.0/(ncnp*Vtv);
			`expLinA(expn,Vbcp,maxvIbcnp,afac);
			Ibcp	=  ibcip_t*(expi-1.0)+ibcnp_t*(expn-1.0);
		end else
			Ibcp	=  0.0;

//
//	Transit time
//

		sgIf	=  Ifi>0.0?1.0:0.0;
		rIf	=  Ifi*sgIf*Iitf;
		mIf	=  rIf/(rIf+1.0);
		arg	=  Vbci*Ivtf/1.44;
		`expLin(expi,arg,VmaxExp);
		tff	=  tf*(1.0+qtf*q1)*(1.0+xtf*expi*(sltf+mIf*mIf)*sgIf);

//
//	Charge elements
//

		Qbe	=  cje_t*qdbe*wbe+tff*Ifi/qb;
		Qbex	=  cje_t*qdbex*(1.0-wbe);
		Qbc	=  cjc_t*qdbc+tr*Iri+qco*Kbci;
		Qbcx	=  qco*Kbcx;
		Qbep	=  cjep_t*qdbep+tr*Ifp;
		Qbcp	=  cjcp_t*qdbcp+ccso*Vbcp;
		Qbeo	=  Vbe*cbeo;
		Qbco	=  Vbc*cbco;

//
//	Resistance branch currents
//

		Ircx	=  Vrcx*Gcx;
		rKp1	=  (Kbci+1.0)/(Kbcx+1.0);
		Iohm	=  (Vrci+Vtv*(Kbci-Kbcx-ln(rKp1)))*Gci;
		derf	=  Ivo*Iohm/(Gci*(1.0+0.5*Ivo*Ihrcf*sqrt(Vrci*Vrci+0.01)));
		Irci	=  Iohm/sqrt(1+derf*derf);
		Irbx	=  Vrbx*Gbx;
		Irbi	=  Vrbi*qb*Gbi;
		Ire	=  Vre*Ge;
		Irbp	=  Vrbp*qbp*Gbp;
		Irs	=  Vrs*Gs;

//
//	Thermal network elements
//

			Ith	= -(Ibe*Vbei+Ibc*Vbci+(Itxf-Itzr)*Vcei+Ibex*Vbex+Ibep*Vbep+Irs*Vrs+Ibcp*Vbcp+Iccp*Vcep+Ircx*Vrcx+Irci*Vrci+Irbx*Vrbx+Irbi*Vrbi+Ire*Vre+Irbp*Vrbp);
		Irth	=  Vrth*Gth;
		Qcth	=  Vrth*cth;

//
//	Excess phase elements
//

		Ixf1	=  Vxf2-Itzf;
		Ixf2	=  Vxf2-Vxf1;
		Qxf1	=  td*Vxf1;
		Qxf2	=  td*Vxf2/3;

//
//	Add gmin current to diode-like branches
//

		Ibe	=  Ibe  + `GMIN * Vbei;
		Ibex	=  Ibex + `GMIN * Vbex;
		Ibep	=  Ibep + `GMIN * Vbep;
		Ibc	=  Ibc  + `GMIN * Vbci;
		Ibcp	=  Ibcp + `GMIN * Vbcp;

//
//	Apply multiplicity factor and device type (+1 or -1, polarity)
//	(resistors apart from Irci do not depend on polarity)
//

		Ibe	=  VBICtype * meff * Ibe;
		Ibex	=  VBICtype * meff * Ibex;
		Itzf	=  VBICtype * meff * Itzf;
		Itxf	=  VBICtype * meff * Itxf;
		Itzr	=  VBICtype * meff * Itzr;
		Ibc	=  VBICtype * meff * Ibc;
		Ibep	=  VBICtype * meff * Ibep;
		Ircx	=             meff * Ircx;
		Irci	=  VBICtype * meff * Irci;
		Irbx	=             meff * Irbx;
		Irbi	=             meff * Irbi;
		Ire	=             meff * Ire;
		Irbp	=             meff * Irbp;
		Qbe	=  VBICtype * meff * Qbe;
		Qbex	=  VBICtype * meff * Qbex;
		Qbc	=  VBICtype * meff * Qbc;
		Qbcx	=  VBICtype * meff * Qbcx;
		Qbep	=  VBICtype * meff * Qbep;
		Qbeo	=             meff * Qbeo;
		Qbco	=             meff * Qbco;
		Ibcp	=  VBICtype * meff * Ibcp;
		Iccp	=  VBICtype * meff * Iccp;
		Irs	=             meff * Irs;
		Qbcp	=  VBICtype * meff * Qbcp;
		Ith	=             meff * Ith;
		Irth	=             meff * Irth;
		Qcth	=             meff * Qcth;

//
//	Branch contributions to VBIC model
//

		I(b_bei) <+ Ibe;
		I(b_bex) <+ Ibex;
		I(b_cei) <+ Itxf;
		I(b_eci) <+ Itzr;
		I(b_bci) <+ Ibc;
		I(b_bep) <+ Ibep;
		I(b_rcx) <+ Ircx;
		I(b_rci) <+ Irci;
		I(b_rbx) <+ Irbx;
		I(b_rbi) <+ Irbi;
		I(b_re)  <+ Ire;
		I(b_rbp) <+ Irbp;
		I(b_bei) <+ ddt(Qbe);
		I(b_bex) <+ ddt(Qbex);
		I(b_bci) <+ ddt(Qbc);
		I(b_bcx) <+ ddt(Qbcx);
		I(b_bep) <+ ddt(Qbep);
		I(b_be)  <+ ddt(Qbeo);
		I(b_bc)  <+ ddt(Qbco);
		I(b_bcp) <+ Ibcp;
		I(b_cep) <+ Iccp;
		I(b_rs)  <+ Irs;
		I(b_bcp) <+ ddt(Qbcp);
		I(b_xf1) <+ Ixf1;
		I(c_xf1) <+ ddt(Qxf1);
		I(b_xf2) <+ Ixf2;
		I(c_xf2) <+ ddt(Qxf2);
		I(b_rth) <+ Irth;
		I(b_rth) <+ ddt(Qcth);
		I(b_ith) <+ Ith;

`ifdef insideADMS
`else
	// <noise>
	begin	// begin noise block
		I(b_bei) <+ white_noise(2*`QQ*abs(Ibe))+flicker_noise(meff*kfn*pow(abs(Ibe/meff),afn),bfn);
		I(b_bex) <+ white_noise(2*`QQ*abs(Ibex))+flicker_noise(meff*kfn*pow(abs(Ibex/meff),afn),bfn);
		I(b_cei) <+ white_noise(2*`QQ*abs(Itzf));
		I(b_bep) <+ white_noise(2*`QQ*abs(Ibep))+flicker_noise(meff*kfn*pow(abs(Ibep/meff),afn),bfn);
		I(b_rcx) <+ white_noise(4*`KB*Tdev*Gcx*meff);
		I(b_rci) <+ white_noise(4*`KB*Tdev*((abs(Irci)+1.0e-10*Gci)/(abs(Vrci)+1.0e-10))*meff);
		I(b_rbx) <+ white_noise(4*`KB*Tdev*Gbx*meff);
		I(b_rbi) <+ white_noise(4*`KB*Tdev*qb*Gbi*meff);
		I(b_re)  <+ white_noise(4*`KB*Tdev*Ge*meff);
		I(b_rbp) <+ white_noise(4*`KB*Tdev*qbp*Gbp*meff);
		I(b_cep) <+ white_noise(2*`QQ*abs(Iccp));
		I(b_rs)  <+ white_noise(4*`KB*Tdev*Gs*meff);
	end	// end noise block
	// </noise>
`endif

	`FINAL_STEP
	begin

		// conductances of resistance elements
		Ircx_Vrcx = `DERIVATE(Ircx,V(b_rcx));
		Irci_Vrci = `DERIVATE(Irci,V(b_rci));
		Irbx_Vrbx = `DERIVATE(Irbx,V(b_rbx));
		Irbi_Vrbi = `DERIVATE(Irbi,V(b_rbi));
		Irbp_Vrbp = `DERIVATE(Irbp,V(b_rbp));
		Ire_Vre   = `DERIVATE(Ire,V(b_re));
		Irs_Vrs   = `DERIVATE(Irs,V(b_rs));

		// conductances
		gm     = `DERIVATE(Itzf,V(b_bei))+`DERIVATE(Itzf,V(b_bci))-`DERIVATE(Itzr,V(b_bei))-`DERIVATE(Itzr,V(b_bci));
		go     = `DERIVATE(Itzr,V(b_bci))-`DERIVATE(Itzf,V(b_bci));
		gpi    = `DERIVATE(Ibe,V(b_bei));
		gpix   = `DERIVATE(Ibex,V(b_bex));
		gmu    = `DERIVATE(Ibc,V(b_bci));
		gmux   = `DERIVATE(Ibep,V(b_bep));

		// conductances (for backward compatibility with level 2)
		ibc_vbci	= gmu;
		ibcp_vbcp	= `DERIVATE(Ibcp,V(b_bcp));
		ibe_vbei	= gpi;
		ibep_vbep	= gmux;
		ibex_vbex	= gpix;
		iccp_vbci	= `DERIVATE(Iccp,V(b_bci));
		iccp_vbcp	= `DERIVATE(Iccp,V(b_bcp));
		iccp_vbep	= `DERIVATE(Iccp,V(b_bep));
		irbi_vbci	= `DERIVATE(Irbi,V(b_bci));
		irbi_vbei	= `DERIVATE(Irbi,V(b_bei));
		irbp_vbci	= `DERIVATE(Irbp,V(b_bci));
		irci_vbci	= `DERIVATE(Irci,V(b_bci));
		itzf_vbci	= `DERIVATE(Itzf,V(b_bci));
		itzf_vbei	= `DERIVATE(Itzf,V(b_bei));
		itzr_vbci	= `DERIVATE(Itzr,V(b_bci));
		itzr_vbei	= `DERIVATE(Itzr,V(b_bei));

		// capacitances
		cpi    = `DERIVATE(Qbe,V(b_bei));
		cpix   = `DERIVATE(Qbex,V(b_bex));
		cmu    = `DERIVATE(Qbc,V(b_bci));
		cmux   = `DERIVATE(Qbcx,V(b_bcx));
		cpip   = `DERIVATE(Qbep,V(b_bep));
		ccs    = `DERIVATE(Qbcp,V(b_bcp));
		cbeoXX = `DERIVATE(Qbeo,V(b_be));
		cbcoXX = `DERIVATE(Qbco,V(b_bc));

		// crude estimate of ft, intrinsic component just has
		// overlap capacitances added, and the extrinsic ft
		// will be affected by extrinsic elements like Rc*Cc,
		// proper terminal simulations are needed to get the
		// actual ft, this is just an estimate: gm/(2*pi*Cbase)
		ft_int = 0.1591549*abs(gm)/(1.0e-20+abs(cpi+cpix+cmu+cmux+cpip+cbeoXX+cbcoXX));

		// capacitances (for backward compatibility with level 2)
		qbc_vbci  = cmu;
		qbco_vbc  = cbcoXX;
		qbcp_vbcp = ccs;
		qbcx_vbcx = `DERIVATE(Qbcx,V(b_bcx));	// not equivalent to level 2, Vbcx=Vbci-Vrcx
		qbe_vbci  = `DERIVATE(Qbe,V(b_bci));
		qbe_vbei  = cpi;
		qbeo_vbe  = cbeoXX;
		qbep_vbci = `DERIVATE(Qbep,V(b_bci));
		qbep_vbep = cpip;
		qbex_vbex = cpix;


		// resistances
		rci_mod = Gci/((abs(Irci)+1.0e-10*Gci)/(abs(Vrci)+1.0e-10));
		rci_eff = `DERIVATE(Irci,V(b_rci));
                rci_eff = (rci_eff==0)?rci_eff:(1/rci_eff);
		rbi_eff = (qb*Gbi);
                rbi_eff = (rbi_eff==0)?rbi_eff:((1/rbi_eff)/meff);
		rbp_eff = (qbp*Gbp);
                rbp_eff = (rbp_eff==0)?rbp_eff:((1/rbp_eff)/meff);

		// terminal currents
		ic   = Itzf - Itzr - Ibc - Ibcp - Ibep;
		ib   = Ibe + Ibex + Ibc + Ibep + Iccp;
		ie   = -Itzf + Itzr - Ibe - Ibex;
		isub = -ic - ib - ie;
		ic   = VBICtype * ic;
		ib   = VBICtype * ib;
		ie   = VBICtype * ie;
		isub = VBICtype * isub;

		// power currents
		powerT =-Ith;
		powerD = Irth;
		tauTh  = cth/Gth;

		// power currents (for backward compatibility with level 2)
		p = powerT;

		// power currents (for backward compatibility with level 2)
		trise = dtemp;

		// terminal voltages
		Vce = V(b_be)-V(b_bc);

	end

	end
endmodule