```
//
//	Verilog-A definition of VBIC
//

`include "disciplines.h"

`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 psibi(P,Ea,Vtv,rT) \
psiio = 2.0*(Vtv/rT)*log(exp(0.5*P*rT/Vtv)-exp(-0.5*P*rT/Vtv)); \
psiin = psiio*rT-3.0*Vtv*log(rT)-Ea*(rT-1.0); \
psibi = psiin+2.0*Vtv*log(0.5*(1.0+sqrt(1.0+4.0*exp(-psiin/Vtv))));
`define 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(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 \
qjrt  = 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; \
qjrt =  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); \
qjrt =  qlo+pow((1.0-FC),(-M))*(V-vl+vl0)-q0; \
end \
end
`define avalm(V,P,M,AV1,AV2) \
vl    = 0.5*(sqrt((P-V)*(P-V)+0.01)+(P-V)); \
avalm = AV1*vl*\$limexp(-AV2*pow(vl,(M-1.0)));

//
//	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.
//
//	The excess phase version does not need to be implemented strictly
//	using 3 extra system variables, as at first seems to be the
//	case for MNA, the inductor current and the two node voltages.
//	The resistor in the excess phase network is 1 ohm, therefore
//	the voltage at node xf2 is the same as the current through
//	the inductor. The extra system equations to solve are then:
//	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.
//

//`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
//

`ifdef ElectroThermal
`ifdef ThreeTerminal
module vbic(c,b,e,dt);
`else
module vbic(c,b,e,s,dt);
`endif
`else
`ifdef ThreeTerminal
module vbic(c,b,e);
`else
module vbic(c,b,e,s);
`endif
`endif

//
//	Node definitions
//

`ifdef ElectroThermal
`ifdef ThreeTerminal
inout		c,b,e,dt;			// external nodes
electrical	c,b,e,dt;			// external nodes
`ifdef ExcessPhase
electrical	cx,ci,bx,bi,ei,bp,xf1,xf2;	// internal nodes
`else
electrical	cx,ci,bx,bi,ei,bp;		// internal nodes
`endif
`else
inout		c,b,e,s,dt;			// external nodes
electrical	c,b,e,s,dt;			// external nodes
`ifdef ExcessPhase
electrical	cx,ci,bx,bi,ei,bp,si,xf1,xf2;	// internal nodes
`else
electrical	cx,ci,bx,bi,ei,bp,si;		// internal nodes
`endif
`endif
`else
`ifdef ThreeTerminal
inout		c,b,e;				// external nodes
electrical	c,b,e;				// external nodes
`ifdef ExcessPhase
electrical	cx,ci,bx,bi,ei,bp,xf1,xf2;	// internal nodes
`else
electrical	cx,ci,bx,bi,ei,bp;		// internal nodes
`endif
`else
inout		c,b,e,s;			// external nodes
electrical	c,b,e,s;			// external nodes
`ifdef ExcessPhase
electrical	cx,ci,bx,bi,ei,bp,si,xf1,xf2;	// internal nodes
`else
electrical	cx,ci,bx,bi,ei,bp,si;		// internal nodes
`endif
`endif
`endif

//
//	Branch definitions
//

branch (b ,e )		be;			//           base-emit
branch (b ,c )		bc;			//           base-coll
branch (bi,ei)		bei;			// intrinsic base-emit
branch (bx,ei)		bex;			// extrinsic base-emit
branch (bi,ci)		bci;			// intrinsic base-coll
branch (bi,cx)		bcx;			// extrinsic base-coll
branch (ci,ei)		cei;			// intrinsic coll-emit
branch (ei,ci)		eci;			// intrinsic emit-coll
branch (bx,bp)		bep;			// parasitic base-emit
branch (e ,ei)		re;			// emit resistance
branch (c ,cx)		rcx;			// coll resistance, extrinsic
branch (cx,ci)		rci;			// coll resistance, intrinsic
branch (b ,bx)		rbx;			// base resistance, extrinsic
branch (bx,bi)		rbi;			// base resistance, intrinsic
branch (bp,cx)		rbp;			// base resistance, parasitic
`ifdef ThreeTerminal
`else
branch (si,bp)		bcp;			// parasitic base-coll
branch (bx,si)		cep;			// parasitic coll-emit
branch (s ,si)		rs;			// subs resistance
`endif
`ifdef ElectroThermal
branch (dt)		rth;			// local thermal branch
branch (dt)		ith;			// local thermal branch
`endif
`ifdef ExcessPhase
branch (xf1)		izf;			// zero-phase   current
branch (xf1)		cxf;			// excess-phase capacitor
branch (xf1,xf2)	lxf;			// excess-phase inductor
branch (xf2)		rxf;			// excess-phase resistor
`endif

//
//	Parameter definitions
//

parameter	real	TNOM	=  27.0;
parameter	real	RCX	=   0.0		from[0.0:inf];
parameter	real	RCI	=   0.0		from[0.0:inf];
parameter	real	VO	=   0.0		from[0.0:inf];
parameter	real	GAMM	=   0.0		from[0.0:inf];
parameter	real	HRCF	=   0.0		from[0.0:inf];
parameter	real	RBX	=   0.0		from[0.0:inf];
parameter	real	RBI	=   0.0		from[0.0:inf];
parameter	real	RE	=   0.0		from[0.0:inf];
parameter	real	RS	=   0.0		from[0.0:inf];
parameter	real	RBP	=   0.0		from[0.0:inf];
parameter	real	IS	=   1.0e-16	from(0.0:inf];
parameter	real	NF	=   1.0		from(0.0:inf];
parameter	real	NR	=   1.0		from(0.0:inf];
parameter	real	FC	=   0.9		from[0.0:1.0);
parameter	real	CBEO	=   0.0		from[0.0:inf];
parameter	real	CJE	=   0.0		from[0.0:inf];
parameter	real	PE	=   0.75	from(0.0:inf];
parameter	real	ME	=   0.33	from(0.0:1.0];
parameter	real	AJE	=  -0.5;
parameter	real	CBCO	=   0.0		from[0.0:inf];
parameter	real	CJC	=   0.0		from[0.0:inf];
parameter	real	QCO	=   0.0		from[0.0:inf];
parameter	real	CJEP	=   0.0		from[0.0:inf];
parameter	real	PC	=   0.75	from(0.0:inf];
parameter	real	MC	=   0.33	from(0.0:1.0];
parameter	real	AJC	=  -0.5;
parameter	real	CJCP	=   0.0		from[0.0:inf];
parameter	real	PS	=   0.75	from(0.0:inf];
parameter	real	MS	=   0.33	from(0.0:1.0];
parameter	real	AJS	=  -0.5;
parameter	real	IBEI	=   1.0e-18	from(0.0:inf];
parameter	real	WBE	=   1.0		from[0.0:1.0];
parameter	real	NEI	=   1.0		from(0.0:inf];
parameter	real	IBEN	=   0.0		from[0.0:inf];
parameter	real	NEN	=   2.0		from(0.0:inf];
parameter	real	IBCI	=   1.0e-16	from(0.0:inf];
parameter	real	NCI	=   1.0		from(0.0:inf];
parameter	real	IBCN	=   0.0		from[0.0:inf];
parameter	real	NCN	=   2.0		from(0.0:inf];
parameter	real	AVC1	=   0.0		from[0.0:inf];
parameter	real	AVC2	=   0.0		from[0.0:inf];
parameter	real	ISP	=   0.0		from[0.0:inf];
parameter	real	WSP	=   1.0		from[0.0:1.0];
parameter	real	NFP	=   1.0		from(0.0:inf];
parameter	real	IBEIP	=   0.0		from[0.0:inf];
parameter	real	IBENP	=   0.0		from[0.0:inf];
parameter	real	IBCIP	=   0.0		from[0.0:inf];
parameter	real	NCIP	=   1.0		from(0.0:inf];
parameter	real	IBCNP	=   0.0		from[0.0:inf];
parameter	real	NCNP	=   2.0		from(0.0:inf];
parameter	real	VEF	=   0.0		from[0.0:inf];
parameter	real	VER	=   0.0		from[0.0:inf];
parameter	real	IKF	=   0.0		from[0.0:inf];
parameter	real	IKR	=   0.0		from[0.0:inf];
parameter	real	IKP	=   0.0		from[0.0:inf];
parameter	real	TF	=   0.0		from[0.0:inf];
parameter	real	QTF	=   0.0		from[0.0:inf];
parameter	real	XTF	=   0.0		from[0.0:inf];
parameter	real	VTF	=   0.0		from[0.0:inf];
parameter	real	ITF	=   0.0		from[0.0:inf];
parameter	real	TR	=   0.0		from[0.0:inf];
parameter	real	TD	=   0.0		from[0.0:inf];
parameter	real	KFN	=   0.0		from[0.0:inf];
parameter	real	AFN	=   1.0		from(0.0:inf];
parameter	real	BFN	=   1.0		from(0.0:inf];
parameter	real	XRE	=   0;
parameter	real	XRBI	=   0;
parameter	real	XRCI	=   0;
parameter	real	XRS	=   0;
parameter	real	XVO	=   0;
parameter	real	EA	=   1.12;
parameter	real	EAIE	=   1.12;
parameter	real	EAIC	=   1.12;
parameter	real	EAIS	=   1.12;
parameter	real	EANE	=   1.12;
parameter	real	EANC	=   1.12;
parameter	real	EANS	=   1.12;
parameter	real	XIS	=   3.0;
parameter	real	XII	=   3.0;
parameter	real	XIN	=   3.0;
parameter	real	TNF	=   0.0;
parameter	real	TAVC	=   0.0;
parameter	real	RTH	=   0.0		from[0.0:inf];
parameter	real	CTH	=   0.0		from[0.0:inf];
parameter	real	VRT	=   0.0		from[0.0:inf];
parameter	real	ART	=   0.1		from(0.0:inf];
parameter	real	CCSO	=   0.0		from[0.0:inf];
parameter	real	QBM	=   0.0;
parameter	real	NKF	=   0.5		from(0.0:inf];
parameter	real	XIKF	=   0;
parameter	real	XRCX	=   0;
parameter	real	XRBX	=   0;
parameter	real	XRBP	=   0;
parameter	real	ISRR	=   1.0		from(0.0:inf];
parameter	real	XISR	=   0.0;
parameter	real	DEAR	=   0.0;
parameter	real	EAP	=   1.12;
parameter	real	VBBE	=   0.0;
parameter	real	NBBE	=   1.0		from(0.0:inf];
parameter	real	IBBE	=   1.0e-6;
parameter	real	TVBBE1	=   0.0;
parameter	real	TVBBE2	=   0.0;
parameter	real	TNBBE	=   0.0;
parameter	real	EBBE	=   0.0;
parameter	real	DTEMP	=   0.0;
parameter	real	VERS	=   1.2;
parameter	real	VREV	=   0.0;

//
//	parameter	alias	TNOM	TN0M, TREF
//	parameter	alias	VO	V0
//	parameter	alias	GAMM	GAMMA
//	parameter	alias	CBEO	CBE0
//	parameter	alias	CBCO	CBC0
//	parameter	alias	CCSO	CCS0
//	parameter	alias	QCO	QC0
//	parameter	alias	XVO	XV0
//	parameter	alias	DTEMP	DTMP
//	parameter	alias	VERS	VERSION
//	parameter	alias	VREV	REV REVISION
//

real	ISatT,ISRRatT,IKFatT,IBEIatT,IBCIatT,ISPatT,IBENatT,IBCNatT;
real	IBEIPatT,IBENPatT,IBCIPatT,IBCNPatT;
real	RCXatT,RCIatT,RBXatT,RBIatT,REatT,RSatT,RBPatT;
real	PEatT,PCatT,PSatT;
real	CJEatT,CJCatT,CJEPatT,CJCPatT;
real	NFatT,NRatT,AVC2atT,VBBEatT,NBBEatT,GAMMatT,VOatT,EBBEatT;
real	Tdev,Tini,rT,dT;
real	IVEF,IVER,IIKF,IIKR,IIKP,IVO,IHRCF,IVTF,IITF,slTF,LEP,CEP;

real	Vtv,Vcbj,Ifi,Iri,Itzf,Ixzf,Itxf,Ixxf,Itzr,q1z,q1,q2,qb,Ifp,Irp,Iccp,q2p,qbp;
real	Ibe,Ibex,Ibcj,Ibc,Ibep,Ibcp,Igc,avalf;
real	Ircx,Irci,Irbx,Irbi,Ire,Irbp,Irs;
real	Kbci,Kbcx,rKp1,Iohm,derf;
real	argi,expi,argn,expn,argx,expx;
real	qdbe,qdbex,qdbc,qdbep,qdbcp;
real	sgIf,rIf,mIf,tff;
real	Qbe,Qbex,Qbc,Qbcx,Qbep,Qbcp,Qbeo,Qbco;
real	Ith,Irth,Qcth;
real	Qcxf,Flxf;
`ifdef WithHomotopy
real	LAMBDA,hf,hr,hfp,hrp;
`endif

analog begin
`ifdef WithHomotopy
LAMBDA	=  1.0;
`endif

//
//		Temperature mappings
//

Tini	=  `TABS+TNOM;
`ifdef ElectroThermal
Tdev	=  \$temperature+DTEMP+V(rth);
`else
Tdev	=  \$temperature+DTEMP;
`endif
Vtv	=  `KB*Tdev/`QQ;
rT	=  Tdev/Tini;
dT	=  Tdev-Tini;
IKFatT	=  IKF*pow(rT,XIKF);
RCXatT	=  RCX*pow(rT,XRCX);
RCIatT	=  RCI*pow(rT,XRCI);
RBXatT	=  RBX*pow(rT,XRBX);
RBIatT	=  RBI*pow(rT,XRBI);
REatT	=  RE*pow(rT,XRE);
RSatT	=  RS*pow(rT,XRS);
RBPatT	=  RBP*pow(rT,XRBP);
ISatT	=  IS*pow((pow(rT,XIS)*exp(-EA*(1.0-rT)/Vtv)),(1.0/NF));
ISRRatT	=  ISRR*pow((pow(rT,XISR)*exp(-DEAR*(1.0-rT)/Vtv)),(1.0/NR));
ISPatT	=  ISP*pow((pow(rT,XIS)*exp(-EAP*(1.0-rT)/Vtv)),(1.0/NFP));
IBEIatT	=  IBEI*pow((pow(rT,XII)*exp(-EAIE*(1.0-rT)/Vtv)),(1.0/NEI));
IBENatT	=  IBEN*pow((pow(rT,XIN)*exp(-EANE*(1.0-rT)/Vtv)),(1.0/NEN));
IBCIatT	=  IBCI*pow((pow(rT,XII)*exp(-EAIC*(1.0-rT)/Vtv)),(1.0/NCI));
IBCNatT	=  IBCN*pow((pow(rT,XIN)*exp(-EANC*(1.0-rT)/Vtv)),(1.0/NCN));
IBEIPatT=  IBEIP*pow((pow(rT,XII)*exp(-EAIC*(1.0-rT)/Vtv)),(1.0/NCI));
IBENPatT=  IBENP*pow((pow(rT,XIN)*exp(-EANC*(1.0-rT)/Vtv)),(1.0/NCN));
IBCIPatT=  IBCIP*pow((pow(rT,XII)*exp(-EAIS*(1.0-rT)/Vtv)),(1.0/NCIP));
IBCNPatT=  IBCNP*pow((pow(rT,XIN)*exp(-EANS*(1.0-rT)/Vtv)),(1.0/NCNP));
NFatT	=  NF*(1.0+dT*TNF);
NRatT	=  NR*(1.0+dT*TNF);
AVC2atT	=  AVC2*(1.0+dT*TAVC);
VBBEatT	=  VBBE*(1.0+dT*(TVBBE1+dT*TVBBE2));
NBBEatT	=  NBBE*(1.0+dT*TNBBE);
PEatT	=  `psibi(PE,EAIE,Vtv,rT);
PCatT	=  `psibi(PC,EAIC,Vtv,rT);
PSatT	=  `psibi(PS,EAIS,Vtv,rT);
CJEatT	=  CJE*pow(PE/PEatT,ME);
CJCatT	=  CJC*pow(PC/PCatT,MC);
CJEPatT	=  CJEP*pow(PC/PCatT,MC);
CJCPatT	=  CJCP*pow(PS/PSatT,MS);
GAMMatT	=  GAMM*pow(rT,XIS)*exp(-EA*(1.0-rT)/Vtv);
VOatT	=  VO*pow(rT,XVO);
EBBEatT	=  exp(-VBBEatT/(NBBEatT*Vtv));

//
//		Parameter mappings
//

IVEF	=  VEF>0.0   ? 1.0/VEF    : 0.0;
IVER	=  VER>0.0   ? 1.0/VER    : 0.0;
IIKF	=  IKF>0.0   ? 1.0/IKFatT : 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/VOatT  : 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;
`ifdef ExcessPhase
LEP	=  TD>0.0    ? TD/3.0     : 0.0;
CEP	=  TD>0.0    ? TD         : 0.0;
`endif

//
//		Electrical branch constituent relations
//

qdbe	=  `qj(V(bei),PEatT,ME,FC,AJE);
qdbex	=  `qj(V(bex),PEatT,ME,FC,AJE);
qdbc	=  `qjrt(V(bci),PCatT,MC,FC,AJC,VRT,ART);
qdbep	=  `qjrt(V(bep),PCatT,MC,FC,AJC,VRT,ART);
`ifdef ThreeTerminal
`else
if (CJCP>0.0) begin
qdbcp	=  `qj(V(bcp),PSatT,MS,FC,AJS);
end else begin
qdbcp	=  0;
end
`endif
argi	=  V(bei)/(NFatT*Vtv);
expi	=  \$limexp(argi);
Ifi	=  ISatT*(expi-1.0);
argi	=  V(bci)/(NRatT*Vtv);
expi	=  \$limexp(argi);
Iri	=  ISatT*ISRRatT*(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
qb	=  0.5*(q1+pow((pow(q1,1.0/NKF)+4.0*q2),NKF));
end else begin
qb	=  0.5*q1*(1.0+pow((1.0+4.0*q2),NKF));
end
Itzr	=  Iri/qb;
Itzf	=  Ifi/qb;
`ifdef ExcessPhase
Ixzf	= -Itzf;
Itxf	=  V(rxf);
Ixxf	=  V(rxf);
`endif
if (ISP>0.0) begin
argi	=  V(bep)/(NFP*Vtv);
expi	=  \$limexp(argi);
argx	=  V(bci)/(NFP*Vtv);
expx	=  \$limexp(argx);
Ifp	=  ISPatT*(WSP*expi+(1.0-WSP)*expx-1.0);
q2p	=  Ifp*IIKP;
qbp	=  0.5*(1.0+sqrt(1.0+4.0*q2p));		// assumes IKP>4*ISP if IKP>0
`ifdef ThreeTerminal
`else
argi	=  V(bcp)/(NFP*Vtv);
expi	=  \$limexp(argi);
Irp	=  ISPatT*(expi-1.0);
Iccp	=  (Ifp-Irp)/qbp;
`endif
end else begin
Ifp	=  0.0;
qbp	=  1.0;
`ifdef ThreeTerminal
`else
Iccp	=  0.0;
`endif
end

`ifdef WithHomotopy
hf	=  (LAMBDA*IBEIatT/ISatT)+1.0-LAMBDA;
hr	=  (LAMBDA*IBCIatT/ISatT)+1.0-LAMBDA;
if (WBE==1.0) begin
argi	=  V(bei)/(NEI*Vtv);
expi	=  \$limexp(argi);
argn	=  V(bei)/(NEN*Vtv);
expn	=  \$limexp(argn);
if (VBBE>0.0) begin
argx	=  (-VBBEatT-V(bei))/(NBBEatT*Vtv);
expx	=  \$limexp(argx);
Ibe	=  ISatT*hf*(expi-1.0)+IBENatT*(expn-1.0)-IBBE*(expx-EBBEatT);
end else begin
Ibe	=  ISatT*hf*(expi-1.0)+IBENatT*(expn-1.0);
end
Ibex	=  0.0;
end else if (WBE==0.0) begin
Ibe	=  0.0;
argi	=  V(bex)/(NEI*Vtv);
expi	=  \$limexp(argi);
argn	=  V(bex)/(NEN*Vtv);
expn	=  \$limexp(argn);
if (VBBE>0.0) begin
argx	=  (-VBBEatT-V(bex))/(NBBEatT*Vtv);
expx	=  \$limexp(argx);
Ibex	=  ISatT*hf*(expi-1.0)+IBENatT*(expn-1.0)-IBBE*(expx-EBBEatT);
end else begin
Ibex	=  ISatT*hf*(expi-1.0)+IBENatT*(expn-1.0);
end
end else begin
argi	=  V(bei)/(NEI*Vtv);
expi	=  \$limexp(argi);
argn	=  V(bei)/(NEN*Vtv);
expn	=  \$limexp(argn);
if (VBBE>0.0) begin
argx	=  (-VBBEatT-V(bei))/(NBBEatT*Vtv);
expx	=  \$limexp(argx);
Ibe	=  WBE*(ISatT*hf*(expi-1.0)+IBENatT*(expn-1.0)-IBBE*(expx-EBBEatT));
end else begin
Ibe	=  WBE*(ISatT*hf*(expi-1.0)+IBENatT*(expn-1.0));
end
argi	=  V(bex)/(NEI*Vtv);
expi	=  \$limexp(argi);
argn	=  V(bex)/(NEN*Vtv);
expn	=  \$limexp(argn);
if (VBBE>0.0) begin
argx	=  (-VBBEatT-V(bex))/(NBBEatT*Vtv);
expx	=  \$limexp(argx);
Ibex	=  (1.0-WBE)*(ISatT*hf*(expi-1.0)+IBENatT*(expn-1.0)-IBBE*(expx-EBBEatT));
end else begin
Ibex	=  (1.0-WBE)*(ISatT*hf*(expi-1.0)+IBENatT*(expn-1.0));
end
end
argi	=  V(bci)/(NCI*Vtv);
expi	=  \$limexp(argi);
argn	=  V(bci)/(NCN*Vtv);
expn	=  \$limexp(argn);
Ibcj	=  ISatT*hr*(expi-1.0)+IBCNatT*(expn-1.0);
hfp	=  (LAMBDA*IBEIPatT/ISPatT)+1.0-LAMBDA;
if ((IBEIP>0.0)||(IBENP>0.0)) begin
argi	=  V(bep)/(NCI*Vtv);
expi	=  \$limexp(argi);
argn	=  V(bep)/(NCN*Vtv);
expn	=  \$limexp(argn);
Ibep	=  ISPatT*hfp*(expi-1.0)+IBENPatT*(expn-1.0);
end else begin
Ibep	=  0.0;
end
`else
if (WBE==1.0) begin
argi	=  V(bei)/(NEI*Vtv);
expi	=  \$limexp(argi);
argn	=  V(bei)/(NEN*Vtv);
expn	=  \$limexp(argn);
if (VBBE>0.0) begin
argx	=  (-VBBEatT-V(bei))/(NBBEatT*Vtv);
expx	=  \$limexp(argx);
Ibe	=  IBEIatT*(expi-1.0)+IBENatT*(expn-1.0)-IBBE*(expx-EBBEatT);
end else begin
Ibe	=  IBEIatT*(expi-1.0)+IBENatT*(expn-1.0);
end
Ibex	=  0.0;
end else if (WBE==0.0) begin
Ibe	=  0.0;
argi	=  V(bex)/(NEI*Vtv);
expi	=  \$limexp(argi);
argn	=  V(bex)/(NEN*Vtv);
expn	=  \$limexp(argn);
if (VBBE>0.0) begin
argx	=  (-VBBEatT-V(bex))/(NBBEatT*Vtv);
expx	=  \$limexp(argx);
Ibex	=  IBEIatT*(expi-1.0)+IBENatT*(expn-1.0)-IBBE*(expx-EBBEatT);
end else begin
Ibex	=  IBEIatT*(expi-1.0)+IBENatT*(expn-1.0);
end
end else begin
argi	=  V(bei)/(NEI*Vtv);
expi	=  \$limexp(argi);
argn	=  V(bei)/(NEN*Vtv);
expn	=  \$limexp(argn);
if (VBBE>0.0) begin
argx	=  (-VBBEatT-V(bei))/(NBBEatT*Vtv);
expx	=  \$limexp(argx);
Ibe	=  WBE*(IBEIatT*(expi-1.0)+IBENatT*(expn-1.0)-IBBE*(expx-EBBEatT));
end else begin
Ibe	=  WBE*(IBEIatT*(expi-1.0)+IBENatT*(expn-1.0));
end
argi	=  V(bex)/(NEI*Vtv);
expi	=  \$limexp(argi);
argn	=  V(bex)/(NEN*Vtv);
expn	=  \$limexp(argn);
if (VBBE>0.0) begin
argx	=  (-VBBEatT-V(bex))/(NBBEatT*Vtv);
expx	=  \$limexp(argx);
Ibex	=  (1.0-WBE)*(IBEIatT*(expi-1.0)+IBENatT*(expn-1.0)-IBBE*(expx-EBBEatT));
end else begin
Ibex	=  (1.0-WBE)*(IBEIatT*(expi-1.0)+IBENatT*(expn-1.0));
end
end
argi	=  V(bci)/(NCI*Vtv);
expi	=  \$limexp(argi);
argn	=  V(bci)/(NCN*Vtv);
expn	=  \$limexp(argn);
Ibcj	=  IBCIatT*(expi-1.0)+IBCNatT*(expn-1.0);
if ((IBEIP>0.0)||(IBENP>0.0)) begin
argi	=  V(bep)/(NCI*Vtv);
expi	=  \$limexp(argi);
argn	=  V(bep)/(NCN*Vtv);
expn	=  \$limexp(argn);
Ibep	=  IBEIPatT*(expi-1.0)+IBENPatT*(expn-1.0);
end else begin
Ibep	=  0.0;
end
`endif
if (AVC1>0.0) begin
avalf	=  `avalm(V(bci),PCatT,MC,AVC1,AVC2atT);
`ifdef ExcessPhase
Igc	=  (Itxf-Itzr-Ibcj)*avalf;
`else
Igc	=  (Itzf-Itzr-Ibcj)*avalf;
`endif
end else begin
Igc	=  0.0;
end
Ibc	=  Ibcj-Igc;

if (RCX>0.0) begin
Ircx	=  V(rcx)/RCXatT;
end else begin
Ircx	=  0.0;
end
argi	=  V(bci)/Vtv;
expi	=  \$limexp(argi);
argx	=  V(bcx)/Vtv;
expx	=  \$limexp(argx);
`ifdef WithHomotopy
Kbci	=  sqrt(1.0+LAMBDA*GAMMatT*pow(expi,LAMBDA));
Kbcx	=  sqrt(1.0+LAMBDA*GAMMatT*pow(expx,LAMBDA));
`else
Kbci	=  sqrt(1.0+GAMMatT*expi);
Kbcx	=  sqrt(1.0+GAMMatT*expx);
`endif
if (RCI>0.0) begin
rKp1	=  (Kbci+1.0)/(Kbcx+1.0);
Iohm	=  (V(rci)+Vtv*(Kbci-Kbcx-log(rKp1)))/RCIatT;
derf	=  IVO*RCIatT*Iohm/(1.0+0.5*IVO*IHRCF*sqrt(V(rci)*V(rci)+0.01));
Irci	=  Iohm/sqrt(1+derf*derf);
end else begin
Irci	=  0.0;
end
if (RBX>0.0) begin
Irbx	=  V(rbx)/RBXatT;
end else begin
Irbx	=  0.0;
end
if (RBI>0.0) begin
Irbi	=  V(rbi)*qb/RBIatT;
end else begin
Irbi	=  0.0;
end
if (RE>0.0) begin
Ire	=  V(re)/REatT;
end else begin
Ire	=  0.0;
end
if (RBP>0.0) begin
Irbp	=  V(rbp)*qbp/RBPatT;
end else begin
Irbp	=  0.0;
end
`ifdef ThreeTerminal
`else
`ifdef WithHomotopy
hrp	=  (LAMBDA*IBCIPatT/ISPatT)+1.0-LAMBDA;
if ((IBCIP>0.0)||(IBCNP>0.0)) begin
argi	=  V(bcp)/(NCIP*Vtv);
expi	=  \$limexp(argi);
argn	=  V(bcp)/(NCNP*Vtv);
expn	=  \$limexp(argn);
Ibcp	=  ISPatT*hrp*(expi-1.0)+IBCNPatT*(expn-1.0);
end else begin
Ibcp	=  0.0;
end
`else
if ((IBCIP>0.0)||(IBCNP>0.0)) begin
argi	=  V(bcp)/(NCIP*Vtv);
expi	=  \$limexp(argi);
argn	=  V(bcp)/(NCNP*Vtv);
expn	=  \$limexp(argn);
Ibcp	=  IBCIPatT*(expi-1.0)+IBCNPatT*(expn-1.0);
end else begin
Ibcp	=  0.0;
end
`endif
if (RS>0.0) begin
Irs	=  V(rs)/RSatT;
end else begin
Irs	=  0.0;
end
`endif

sgIf	=  Ifi>0.0?1.0:0.0;
rIf	=  Ifi*sgIf*IITF;
mIf	=  rIf/(rIf+1.0);
tff	=  TF*(1.0+QTF*q1)*(1.0+XTF*\$limexp(V(bci)*IVTF/1.44)*(slTF+mIf*mIf)*sgIf);
Qbe	=  CJEatT*qdbe*WBE+tff*Ifi/qb;
Qbex	=  CJEatT*qdbex*(1.0-WBE);
Qbc	=  CJCatT*qdbc+TR*Iri+QCO*Kbci;
Qbcx	=  QCO*Kbcx;
Qbep	=  CJEPatT*qdbep+TR*Ifp;
`ifdef ThreeTerminal
`else
Qbcp	=  CJCPatT*qdbcp+CCSO*V(bcp);
`endif
Qbeo	=  V(be)*CBEO;
Qbco	=  V(bc)*CBCO;

`ifdef ElectroThermal
`ifdef ThreeTerminal
`ifdef ExcessPhase
Ith	= -(Ibe*V(bei)+Ibc*V(bci)+(Itxf-Itzr)*V(cei)+Ibex*V(bex)+Ibep*V(bep)+Ircx*V(rcx)+Irci*V(rci)+Irbx*V(rbx)+Irbi*V(rbi)+Ire*V(re)+Irbp*V(rbp));
`else
Ith	= -(Ibe*V(bei)+Ibc*V(bci)+(Itzf-Itzr)*V(cei)+Ibex*V(bex)+Ibep*V(bep)+Ircx*V(rcx)+Irci*V(rci)+Irbx*V(rbx)+Irbi*V(rbi)+Ire*V(re)+Irbp*V(rbp));
`endif
`else
`ifdef ExcessPhase
Ith	= -(Ibe*V(bei)+Ibc*V(bci)+(Itxf-Itzr)*V(cei)+Ibex*V(bex)+Ibep*V(bep)+Irs*V(rs)+Ibcp*V(bcp)+Iccp*V(cep)+Ircx*V(rcx)+Irci*V(rci)+Irbx*V(rbx)+Irbi*V(rbi)+Ire*V(re)+Irbp*V(rbp));
`else
Ith	= -(Ibe*V(bei)+Ibc*V(bci)+(Itzf-Itzr)*V(cei)+Ibex*V(bex)+Ibep*V(bep)+Irs*V(rs)+Ibcp*V(bcp)+Iccp*V(cep)+Ircx*V(rcx)+Irci*V(rci)+Irbx*V(rbx)+Irbi*V(rbi)+Ire*V(re)+Irbp*V(rbp));
`endif
`endif
if (RTH>0.0) begin
Irth	=  V(rth)/RTH;
end else begin
Irth	=  0.0;
end
Qcth	=  V(rth)*CTH;
`endif
`ifdef ExcessPhase
Flxf	=  LEP*Ixxf;
Qcxf	=  CEP*V(cxf);
`endif

//
//		Branch contributions to VBIC model
//

I(bei)	<+  Ibe;
I(bex)	<+  Ibex;
`ifdef ExcessPhase
I(cei)	<+  Itxf;
`else
I(cei)	<+  Itzf;
`endif
I(eci)	<+  Itzr;
I(bci)	<+  Ibc;
I(bep)	<+  Ibep;
I(rcx)	<+  Ircx;
I(rci)	<+  Irci;
I(rbx)	<+  Irbx;
I(rbi)	<+  Irbi;
I(re)	<+  Ire;
I(rbp)	<+  Irbp;
I(bei)	<+  ddt(Qbe);
I(bex)	<+  ddt(Qbex);
I(bci)	<+  ddt(Qbc);
I(bcx)	<+  ddt(Qbcx);
I(bep)	<+  ddt(Qbep);
I(be)	<+  ddt(Qbeo);
I(bc)	<+  ddt(Qbco);
`ifdef ThreeTerminal
`else
I(bcp)	<+  Ibcp;
I(cep)	<+  Iccp;
I(rs)	<+  Irs;
I(bcp)	<+  ddt(Qbcp);
`endif
`ifdef ElectroThermal
I(rth)	<+  Irth;
I(ith)	<+  Ith;
I(rth)	<+  ddt(Qcth);
`endif
`ifdef ExcessPhase
I(izf)	<+  Ixzf;
I(rxf)	<+  Ixxf;
I(cxf)	<+  ddt(Qcxf);
V(lxf)	<+  ddt(Flxf);
`endif
end
endmodule

```