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