//
// Copyright @ 2015 Chair for Electron Devices and Integrated Circuits (CEDIC) at TU Dresden.
//
// The terms under which the software and associated documentation (the Software) is provided are as the following:
//
// The Software is provided "as is", without warranty of any kind, express or implied, including but not limited to the warranties of merchantability, fitness for a particular purpose and noninfringement. In no event shall the authors or copyright holders be liable for any claim, damages or other liability, whether in an action of contract, tort or otherwise, arising from, out of or in connection with the Software or the use or other dealings in the Software.
//
// CEDIC grants, free of charge, to any users the right to modify, copy, and redistribute the Software, both within the user's organization and externally, subject to the following restrictions:
//
// 1. The users agree not to charge for the code itself but may charge for additions, extensions, or support.
//
// 2. In any product based on the Software, the users agree to acknowledge the Research Group that developed the software. This acknowledgment shall appear in the product documentation.
//
// 3. The users agree to obey all U.S. Government restrictions governing redistribution or export of the software.
//
// 4. The users agree to reproduce any copyright notice which appears on the software on any copy or modification of such made available to others.
//
// Michael Schroter (mschroter@ieee.org)
// May 27, 2015
//

//CNTFET compact model

// Changes for CCAM_v211
// *******************************************************************************************************************************************************************
// 22/08/2019: 
// COMMENT ON $simparam:
// $simparam fucntion is not supported in QUCS so a condition check for the pre-defined macro __VAMS_COMPACT_MODELING__ is implemented.
// If a simulator defines __VAMS_COMPACT_MODELING__ then $simpapram funciton is supported by the simulator
// COMMENT ON NODE COLLAPSING:
// Presently this verilog code permits a minimum of 1 micro-Ohm resistance (Rmin) for any series resistance as well as for thermal resistance rths and rthm. 
// If any of the resistances value drops below this minimum value, then the current in the node is calculated via the nodeVoltage/Rmin instead of NODE COLLAPSING.
// When testing the model with QUCS, node collapsing where 0 resistance is detected is throwing an error and no convergence can be achieved. 
// COMMENT ON THERMAL NETWORK:
// The thermal network and self heating is implemented similar to the one available in - hicumL2V2p4p0.
// COMMENT ON $vt FUNCTION:
// Without an arguement $vt function always returns thermal voltage at circuits ambient temperature. But the contribution of termapeartuer increase due to 
// self heating or from external variation in temperature shpuld be accounted so the variable Tdev is passed as an arguement to $vt function to calculate
// appropriate thermal voltage.
// COMMENT ON LINEAR EMPIRICAL CURRENT MODEL:
// Bugfix for removing the discontinuity in linear empirical current model for VGiSi = 0 V
// COMENT on INCLUDE FILES:
// 'constants.h' and 'disciplines.h' are only defined by Cadence. A standard definition files for Verlig-AMS HDL are 'constants.vams' and 'disciplines.vams'
// ********************************************************************************************************************************************************************

// Changes for CCAM_v212
// *******************************************************************************************************************************************************************
// 22/08/2019: 
// COMMENT ON width and tube density scalable model:
// Expressions of transfer current, contact resistance and tube charge of semiconducting and metallic CNTs are updated for wifth and tube density scaling
// Expression for finger resistance and parasitic capacitance are updated for width scalability.
// ********************************************************************************************************************************************************************

// Changes for CCAM_v220
// *******************************************************************************************************************************************************************
// 19/12/2019: 
// COMMENT ON TYPE SPECIFIER:
// A type specifier parameter is added to switch the model between p-type and n-type.
// ********************************************************************************************************************************************************************

`include "constants.vams"
`include "disciplines.vams"

`define a_cc    1.421e-10
`define a_0     2.4612e-10
`define t_cc    3.033
`define D_0     8/(3*`M_PI*`t_cc*`a_cc)
`define k_Bnorm `P_K/`P_Q
`ifdef __VAMS_COMPACT_MODELING__
    `define Gmin $simparam("gmin",1e-12)
`else
    `define Gmin 1e-12
`endif
`define Rmin    1e-6

module cntfet(d,g,s);

inout d,g,s;
electrical d,g,s,dx,sx,di,gi,si,dim,sim,node_trap;
electrical tnode_s,tnode_m;

// branch definition for thermal nodes
branch (tnode_s) br_sht_s;
branch (tnode_m) br_sht_m;

// -------------------------------------
// Type specifier
parameter integer devtype   = 1         from [-1:1];    // -1.. PMOS, 1 NMOS

// DC parameters s-tubes
parameter real w            = 1e-6      from [0:inf);   // total transistor width - [m]
parameter real dens         = 1         from [1:inf);	// density of CNT per um - [1/um]
parameter real pst          = 1         from (0:1];	    // semiconducting CNTs to total CNTs ratio - [-]
parameter integer currmod   = 1         from [1:2];     // 1 old model, 2 linear model - [-]
parameter real idstn        = 0.3       from [0:10];    // current factor per semiconducting tube for n-type conduction - [A]
parameter real vthin        = 0.15      from [-10:10];  // threshold voltage n-type - [V]
parameter real vth0in       = 5         from (0:100];   // fit parameter n-type current - [-]
parameter real athin        = 2         from (0:100];   // fit parameter n-type current - [-]
parameter real sthin        = 0.2       from [0:1];     // fit parameter for n-type subthreshold slope - [-]
parameter real dvthin       = 0.3       from [0:10];    // new linear model second threshold at vth0in + dvthin - [-]
parameter real m1in         = 1         from [0:100];   // new linear model fit parameter for n-type current at first threshold - [-]
parameter real m2in         = 1         from [0:100];   // new linear model fit parameter for n-type current at second threshold - [-]
parameter real m3in         = 5         from [0:20];    // new linear model fit parameter for n-type current for VGiSi <= 0 - [-]
parameter real pin          = 1         from [0:1];     // new linear model gm factor for n-type current after second threshold  - [-]
parameter real idstp        = 0.15      from [0:10];    // current factor per semiconducting tube for p-type conduction - [A]
parameter real vthip        = -0.03     from [-10:10];  // threshold voltage p-type - [V]
parameter real vth0ip       = 5         from (0:100];   // fit parameter p-type current - [-]
parameter real athip        = 2         from (0:100];   // fit parameter p-type current - [-]
parameter real sthip        = 0.2       from [0:1];     // fit parameter for subthreshold slope - [-]
parameter real dvthip       = -0.3      from [-10:10];  // new linear model second threshold at vth0ip + dvthip - [-]
parameter real m1ip         = 1         from [0:100];   // new linear model fit parameter for p-type current at first threshold - [-]
parameter real m2ip         = 1         from [0:100];   // new linear model fit parameter for p-type current at second threshold - [-]
parameter real m3ip         = 5         from [0:20];    // new linear model fit parameter for p-type current for VGiSi >= 0 - [-]
parameter real pip          = 1         from [0:1];     // new linear model gm factor for for p-type current after second threshold  - [-]
parameter real vdcr         = 2         from (0:50];    // fit parameter for current calculation - [-]
parameter real betdcr       = 1         from (0:50];    // fit parameter for current calculation - [-]
parameter real betds        = 1         from [0:10];    // new linear model fit parameter for VDS dependency - [-]
parameter real smss         = 10        from [0:100];   // fit parameter transition to subthreshold region - [-]
parameter real facss        = 1         from [1:10];    // fit parameter DIBL - [-]
parameter real rscs         = 20        from [0:10M];   // source contact resistance per s-tube - [Ohm]
parameter real rdcs         = 20        from [0:10M];   // drain contact resistance per s-tube - [Ohm]

// -------------------------------------
// AC parameters s-tubes
parameter integer chargemod = 1         from [1:2];     // 0 old model, 1 linear model  - [-]
parameter real ctn0         = 2e-13     from [0:1u];    // charge factor per semiconducting tube for electrons - [C]
parameter real ctp0         = 2e-13     from [0:1u];    // charge factor per semiconducting tube for holes - [C]
parameter real vthqs        = -0.1      from [-10:10];  // threshhold voltage for Qs - [V]
parameter real vth0qs       = 0.7       from (0:100];   // fit parameter for Qs calculation - [-]
parameter real athqs        = 3         from (0:100];   // fit parameter for Qs calculation - [-]
parameter real mqs          = 5         from [0:1000];  // new linear model charge factor - [-]
parameter real vthqd        = -1        from [-10:10];  // threshhold voltage for Qd - [V]
parameter real vth0qd       = 3         from (0:100];   // fit parameter for Qd calculation - [-]
parameter real athqd        = 3         from (0:100];   // fit parameter for Qd calculation - [-]
parameter real mqd          = 5         from [0:1000];  // new linear model charge factor - [-]
parameter real pqsd         = 0.15      from [0:1];     // source/drain partitioning factor - [-]

// -------------------------------------
// parameters m-tubes
parameter real rmta         = 0         from [0:1G];    // low voltage resistance per metallic CNT  - [Ohm]
parameter real amto         = 25        from [0:1M];    // saturation parameter of metallic CNTs - [-]
parameter real cmt          = 0         from [0:1u];    // capacitance factor per metallic tube - [F]
parameter real pqmt         = 0.5       from [0:1];     // source/drain partitioning factor - [-]
parameter real rscm         = 8         from [0:10M];   // source contact resistance per metallic tube  - [Ohm]
parameter real rdcm         = 8         from [0:10M];   // drain contact resistance per metallic tube - [Ohm]

// -------------------------------------
// parameters for noise model
parameter real fanof        = 1.0       from [0:10];    // Fano factor for channel shot noise calculation - [-]
parameter real fanofmt      = 1.0       from [0:10];    // Fano factor for metallic tube shot noise calculation - [-]
parameter real hoogef       = 2e-3      from [0:1];     // Hooge parameter for 1/f-noise calcultation - [-]
parameter real beta_fn      = 1         from [0:10];    // slope of flicker frequency distribution - [-]

// -------------------------------------
// parameters for trap model
parameter real atrap        = 0.5;                      // trap model parameter - [-]
parameter real btrap        = 0.5;                      // trap model parameter - [-]
parameter real ctrap        = 2;                        // trap model parameter - [-]
parameter real w0trap       = 0;                        // lower characteristic frequency of traps - [-]
parameter real strap        = 0         from [0:1];     // trap model parameter - [-]

// -------------------------------------
// external parameters
parameter real rsf          = 0         from [0:1k];    // width normalized source finger resistance - [Ohm/um]
parameter real rdf          = 0         from [0:1k];    // width normalized drain finger resistance - [Ohm/um]
parameter real rg           = 10        from [0:1k];    // width normalized gate finger resistance - [Ohm/um]
parameter real cgspar1      = 0         from [0:1u];    // parasitic gate-source capacitance 1 - [F]
parameter real cgspar2      = 3e-14     from [0:1u];    // width normalized parasitic gate-source capacitance 2 - [F/um]
parameter real cgdpar1      = 0         from [0:1u];    // parasitic gate-drain capacitance 1 - [F]
parameter real cgdpar2      = 3e-14     from [0:1u];    // width normalized parasitic gate-drain capacitance 2 - [F/um]
parameter real cdspar       = 2e-14     from [0:1u];    // width normalized parasitic drain-source capacitance - [F/um]


// -------------------------------------
// parameters for temperature dependence
parameter real tnom         = 27.0;                     // Temperature at which parameters are specified - [C]
parameter real dt           = 0.0;                      // Temperature change w.r.t. chip temperature for particular transistor - [K]
parameter real zetrmta      = 0         from [0:10];    // Temperature dependence parameter for rmta - [-]
parameter real zetamto      = 0         from [0:10];    // Temperature dependence parameter for amto - [-]
parameter real alvthin      = 0         from [0:10];    // Temperature dependence parameter for vthin - [1/K]
parameter real alvth0in     = 0         from [0:10];    // Temperature dependence parameter for vth0in - [1/K]
parameter real alvthip      = 0         from [0:10];    // Temperature dependence parameter for vthip - [1/K]
parameter real alvth0ip     = 0         from [0:10];    // Temperature dependence parameter for vth0ip  - [1/K]
parameter real alidst       = 0         from [0:10];    // Temperature dependence parameter for idstn and idstp - [1/K]
parameter real alct0        = 0         from [0:10];    // Temperature dependence parameter for ctn0 and ctp0 - [1/K]
parameter real alvthqs      = 0         from [0:10];    // Temperature dependence parameter for vthqs - [1/K]
parameter real alvth0qs     = 0         from [0:10];    // Temperature dependence parameter for vth0qs  - [1/K]
parameter real alvthqd      = 0         from [0:10];    // Temperature dependence parameter for vthqd  - [1/K]
parameter real alvth0qd     = 0         from [0:10];    // Temperature dependence parameter for vth0qd  - [1/K]
parameter real alrscon      = 0         from [0:10];    // Temperature dependence parameter for rscs and rdcs - [1/K]
parameter real alrmcon      = 0         from [0:10];    // Temperature dependence parameter for rscm and rdcm - [1/K]

// -------------------------------------
// parameters for self heating
parameter integer flsh      = 0         from [0:1];     // Flag for turning on and off self-heating effect - [-]
parameter real rths         = 1         from [0:1M];    // Thermal resistance of semiconducting CNTs - [K/W]
parameter real rthm         = 1         from [0:1M];    // Thermal resistance of metallic CNTs - [K/W]

// -------------------------------------
real scalingFactor;
(*desc="charge"*) real Qs_cnt,Qd_cnt;
real fds, fgs, fgsp, fgsn, u;
real vg_t, vg_eff_n, vg_eff_p, athn_con, athp_con;
real curr, rmt;
real fnamp;
real rmta_t, amto_t;
real VT, T0, Tamb, Pdiss_s, Pdiss_m;
(*desc="tempSH"*) real delta_T_s, Tdev_s, delta_T_m, Tdev_m;
real idstn_t, idstp_t, ctn0_t, ctp0_t;
real vth0in_t, vthin_t, vth0ip_t, vthip_t;
real vthqs_t, vth0qs_t, vthqd_t, vth0qd_t;
real rscs_t, rdcs_t, rscm_t, rdcm_t;
real w_s, w_u;
real i_trap, R_trap, Cap_trap;
real a_ss, b_ss, c_ss;
real arg1, arg2, smooth;
real vthn2, vthp2;
real fqsn, fqsp, fqdn, fqdp;
real nom,denom;
real V_gisi,V_disi,V_gidi;

//////////////////////////////////////////////////////////////////
analog begin


scalingFactor = 1e6; // <--- Scaling factor for CNT density


////////////////////////////////////////////////////////////////
/////    trap model parameter calculation  //////////////////////
////////////////////////////////////////////////////////////////

// if strap==0 the trap model is turned off
if (strap>0) begin
    R_trap=1;
    Cap_trap=1/pow(w0trap,10);
end
else begin
    R_trap=0;
    Cap_trap=0;
end


////////////////////////////////////////////////////////////////
/////   DIBL modeling  /////////////////////////////////////////
////////////////////////////////////////////////////////////////

if (facss<1) begin
    a_ss=1;
    b_ss=2/(1-facss)-1;
    c_ss=2/(1-facss);
end
else if (facss==1) begin
    a_ss=0;
    b_ss=1;
    c_ss=1;
end
else begin
    a_ss=-1;
    b_ss=(1+facss)/(facss-1);
    c_ss=2/(facss-1);
end


////////////////////////////////////////////////////////////////
/////   Temperature dependant prefactors  //////////////////////
////////////////////////////////////////////////////////////////

T0      = tnom + `P_CELSIUS0; // <--- transformation to Kelvin
Tamb    = $temperature;

if (flsh == 0) begin
    Tdev_s  = Tamb+dt;
    Tdev_m  = Tamb+dt;
end 
else begin
    Tdev_s  = Tamb+dt+V(br_sht_s);
    Tdev_m  = Tamb+dt+V(br_sht_m);
end

delta_T_s   = Tdev_s-T0;
delta_T_m   = Tdev_m-T0;

vth0in_t    = vth0in*(1 + alvth0in*delta_T_s);
vthin_t     = vthin*(1 + alvthin*delta_T_s);
vth0ip_t    = vth0ip*(1 + alvth0ip*delta_T_s);
vthip_t     = vthip*(1 + alvthip*delta_T_s);

idstn_t     = idstn*(1 + alidst*delta_T_s);
idstp_t     = idstp*(1 + alidst*delta_T_s);

vthqs_t     = vthqs*(1 + alvthqs*delta_T_s);
vth0qs_t    = vth0qs*(1 + alvth0qs*delta_T_s);
vthqd_t     = vthqd*(1 + alvthqd*delta_T_s);
vth0qd_t    = vth0qd*(1 + alvth0qd*delta_T_s);

ctn0_t  = ctn0/(1 + alct0*delta_T_s);
ctp0_t  = ctp0/(1 + alct0*delta_T_s);

rmta_t  = rmta*pow(Tdev_m/T0,zetrmta);
amto_t  = amto*pow(Tdev_m/T0,zetamto);

rscs_t  = rscs*(1 + alrscon*delta_T_s);
rdcs_t  = rdcs*(1 + alrscon*delta_T_s);
rscm_t  = rscm*(1 + alrmcon*delta_T_m);
rdcm_t  = rdcm*(1 + alrmcon*delta_T_m);

// Thermal voltage is calculated with the updated device temperature
VT      = $vt(Tdev_s);


////////////////////////////////////////////////////////////////
/////   Network for trap modelling    //////////////////////////
////////////////////////////////////////////////////////////////

i_trap=atrap*V(g,s)+btrap*V(g,s)*V(d,s)+ctrap;

I(node_trap) <+ Cap_trap*ddt(V(node_trap));
I(node_trap) <+ V(node_trap)*`Gmin;

if (R_trap>0) begin
    I(node_trap) <+ -i_trap;
    I(node_trap) <+ V(node_trap)/R_trap;
end
else begin
    V(node_trap) <+ 0;
end

vthin_t     = vthin_t+V(node_trap);
vthip_t     = vthip_t+V(node_trap);

////////////////////////////////////////////////////////////////
/////   type modifier changes   ////////////////////////////////
////////////////////////////////////////////////////////////////

V_gisi = devtype*V(gi,si);
V_disi = devtype*V(di,si);
V_gidi = devtype*V(gi,di);


////////////////////////////////////////////////////////////////
/////   Current model    ///////////////////////////////////////
////////////////////////////////////////////////////////////////

//electron current   fgsn

vg_eff_n=V_gisi-0.5*(V_disi - hypot(V_disi,VT));

if (currmod==1) begin // empirical model
    vg_t = ((vg_eff_n-vthin_t)+hypot(vg_eff_n-vthin_t,sthin))/2;

    u   = 1-vth0in_t/vg_t;

    athn_con = (ln(1+exp(1/athin)));

    if (u/athin>-10) begin
        fgsn    = ln(1+exp(u/athin))/athn_con;
    end
    else begin
        fgsn    = exp(u/athin)/athn_con;
    end
end
else begin // linear empirical model
    vthn2 = vthin_t+dvthin;
    
    nom   = 1+limexp((V_gisi-vthin_t)/(m1in*VT));
    denom = 1+limexp(         -vthin_t /(m1in*VT));
    arg1  = nom*(1/(1+limexp((V_gisi-m3in*VT)/VT))) + (nom/denom)*(1-1/(1+limexp((V_gisi-m3in*VT)/VT)));

    nom   = 1+limexp((V_gisi-pow(abs(V_disi),betds)-vthn2)/(m2in*VT));
    denom = 1+limexp((        -pow(abs(V_disi),betds)-vthn2)/(m2in*VT));
    arg2  = nom/denom;
    
    smooth  = (V_gisi/vthin_t)/pow((1+pow(V_gisi/vthin_t,2)),0.5);
    
    fgsn    = m1in*ln(arg1)-pin*m2in*ln(arg2)*smooth;
end


//hole current   fgsp

vg_eff_p = V_gidi+0.5*(V_disi-hypot(V_disi,VT));

if (currmod==1) begin // empirical model
	vg_t = ((-vg_eff_p+vthip_t)+hypot(-vg_eff_p+vthip_t,sthip))/2;

	u   = 1-vth0ip_t/vg_t;

	athp_con    = (ln(1+exp(1/athip)));

	if (u/athip>-10) begin
		fgsp    = ln(1+exp(u/athip))/athp_con;
	end
	else begin
		fgsp    = exp(u/athip)/athp_con;
	end
end 
else begin // linear empirical model
    vthp2 = vthip_t+dvthip;
    
    nom   = 1+limexp((-V_gisi+vthip_t)/(m1ip*VT));
    denom = 1+limexp(           vthip_t /(m1ip*VT));
    arg1  = nom*(1/(1+limexp((-V_gisi-m3in*VT)/VT))) + (nom/denom)*(1-1/(1+limexp((-V_gisi-m3in*VT)/VT)));

    nom   = 1+limexp((-V_gisi+pow(abs(V_disi),betds)+vthp2)/(m2ip*VT));
    denom = 1+limexp((          pow(abs(V_disi),betds)+vthp2)/(m2ip*VT));
    arg2  = nom/denom;
    
    smooth  = (V_gisi/vthip_t)/pow((1+pow(V_gisi/vthip_t,2)),0.5);
   
    fgsp    = m1ip*ln(arg1)-pip*m2ip*ln(arg2)*smooth;
end


// --> current from holes and electrons

fgs = idstp_t*fgsp+idstn_t*fgsn;


//current  fds

if (idstp==0) begin
    u = V_disi/((tanh(a_ss*smss*(vg_eff_n-vthin_t))+b_ss)/c_ss*vdcr);
end
else if (idstn==0) begin
    u = V_disi/((tanh(a_ss*smss*(-vg_eff_p+vthip_t))+b_ss)/c_ss*vdcr);
end
else begin
    u = V_disi/vdcr;
end

fds     = u/pow(1+pow(pow((u*u + (VT/vdcr)),0.5),betdcr),(1/betdcr));

// --> final current 

curr    = devtype*w*scalingFactor*pst*dens*fgs*fds;


////////////////////////////////////////////////////////////////
/////   Charge  model    ///////////////////////////////////////
////////////////////////////////////////////////////////////////


//Qs

if (chargemod==1) begin // empirical model
	if (V_gisi>vthqs_t) begin
		vg_t    = V_gisi-vthqs_t;
		if (vg_t>1e-6*VT) begin
			u   = 1-vth0qs_t/vg_t;
		end
		else begin
			u   = 1-vth0qs_t/(1e-6*VT);
		end
		w_s     = sqrt(u*u+athqs);
		w_u     = (u+w_s)/(1+sqrt(1+athqs));
		Qs_cnt  = devtype*pqsd*w*scalingFactor*pst*dens*ctn0_t*vg_t*w_u*w_u;
	end
	else begin
		vg_t    = -V_gisi+vthqs_t;
		if (vg_t > 1e-6*VT) begin
			u   = 1-vth0qs_t/vg_t;
		end
		else begin
			u   = 1-vth0qs_t/(1e-6*VT);
		end
		w_s     = sqrt(u*u+athqs);
		w_u     = (u+w_s)/(1+sqrt(1+athqs));
		Qs_cnt  = -devtype*pqsd*w*scalingFactor*pst*dens*ctp0_t*vg_t*w_u*w_u;
	end
end
else begin // linear empirical model
    vg_t   = V_gisi-vthqs_t;

    nom    = 1+exp(    vg_t/mqs/VT);
    denom  = 1+exp(-vthqs_t/mqs/VT);
    fqsn   = mqs*VT*ln(nom/denom );

    nom    = 1+exp(-   vg_t/mqs/VT);
    denom  = 1+exp(-vthqs_t/mqs/VT);
    fqsp   = mqs*VT*ln(nom/denom );

    Qs_cnt = devtype*pqsd*w*scalingFactor*pst*dens*(ctn0_t*fqsn-ctp0_t*fqsp);
end


//Qd

if (chargemod==1) begin // empirical model
	if (V_gidi>vthqd_t) begin
		vg_t=V_gidi-vthqd_t;
		if (vg_t>1e-6*VT) begin
			u=1-vth0qd_t/vg_t;
		end
		else begin
			u=1-vth0qd_t/(1e-6*VT);
		end
		w_s=sqrt(u*u+athqd);
		w_u=(u+w_s)/(1+sqrt(1+athqd));
		Qd_cnt = devtype*(1-pqsd)*w*scalingFactor*pst*dens*ctn0_t*vg_t*w_u*w_u;
	end
	else begin
		vg_t=-V_gidi+vthqd_t;
		if (vg_t>1e-6*VT) begin
			u=1-vth0qd_t/vg_t;
		end
		else begin
			u=1-vth0qd_t/(1e-6*VT);
		end
		w_s=sqrt(u*u+athqd);
		w_u=(u+w_s)/(1+sqrt(1+athqd));
		Qd_cnt = -devtype*(1-pqsd)*w*scalingFactor*pst*dens*ctp0_t*vg_t*w_u*w_u;
	end
end	
else begin // linear empirical model
    vg_t   = V_gidi-vthqd_t;
    
    nom    = 1+exp(    vg_t /mqd/VT);
    denom  = 1+exp(-vthqd_t /mqd/VT);
    fqdn   = mqd*VT*ln( nom / denom );

    nom    = 1+exp(-   vg_t /mqd/VT);
    denom  = 1+exp(-vthqd_t /mqd/VT);
    fqdp   = mqd*VT*ln( nom / denom );

    Qd_cnt = devtype*(1-pqsd)*w*scalingFactor*pst*dens*(ctn0_t*fqdn-ctp0_t*fqdp);
end


////////////////////////////////////////////////////////////////
/////   Circuit elements   /////////////////////////////////////
////////////////////////////////////////////////////////////////


// Transfer current

I(di,si) <+ curr;

// metallic tube current
// is disabled if rmta is set to zero

if(rmta != 0) begin
	rmt = rmta_t+amto_t*hypot(V(dim,sim),sqrt(VT));
	I(dim,sim) <+ (V(dim,sim)/rmt)*w*scalingFactor*(1-pst)*dens;
end
else begin 
	rmt = 0;
	I(dim,sim) <+ 0;
end


// Dynamic currents

I(gi,di) <+ ddt(Qd_cnt);
I(gi,si) <+ ddt(Qs_cnt);

I(gi,dim) <+ pqmt*w*scalingFactor*(1-pst)*dens*cmt*ddt(V(gi,dim));
I(gi,sim) <+ (1-pqmt)*w*scalingFactor*(1-pst)*dens*cmt*ddt(V(gi,sim));


// contact resistance

if (rscs > `Rmin) begin
    I(si,sx) <+ (V(si,sx)/rscs_t)*w*scalingFactor*pst*dens;
end
else begin
    I(si,sx) <+ V(si,sx)/`Rmin;
end

if (rdcs > `Rmin) begin
    I(di,dx) <+ (V(di,dx)/rdcs_t)*w*scalingFactor*pst*dens;
end
else begin
    I(di,dx) <+ V(di,dx)/`Rmin;
end


if (rscm>`Rmin & pst<1) begin
    I(sim,sx) <+ (V(sim,sx)/rscm_t)*w*scalingFactor*(1-pst)*dens;
end
else begin
    I(sim,sx) <+ V(sim,sx)/`Rmin;
end

if (rdcm>`Rmin & pst<1) begin
    I(dim,dx) <+ (V(dim,dx)/rdcm_t)*w*scalingFactor*(1-pst)*dens;
end
else begin
    I(dim,dx) <+ V(dim,dx)/`Rmin;
end


// parasitic finger resistance

if (rsf > `Rmin) begin
    I(sx,s) <+ V(sx,s)/(rsf*w*scalingFactor);
end
else begin
    I(sx,s) <+ V(sx,s)/`Rmin;
end

if (rdf > `Rmin) begin
    I(dx,d) <+ V(dx,d)/(rdf*w*scalingFactor);
end
else begin
    I(dx,d) <+ V(dx,d)/`Rmin;
end

if (rg > `Rmin) begin
    I(gi,g) <+ V(gi,g)/(rg*w*scalingFactor);
end
else begin
    I(gi,g) <+ V(gi,g)/`Rmin;
end


// parasitic finger capacitance

I(g,dx)     <+ cgdpar1*ddt(V(g,dx));
I(g,sx)     <+ cgspar1*ddt(V(g,sx));
I(gi,dx)    <+ cgdpar2*w*scalingFactor*ddt(V(gi,dx));
I(gi,sx)    <+ cgspar2*w*scalingFactor*ddt(V(gi,sx));
I(d,s)      <+ cdspar*w*scalingFactor*ddt(V(d,s));


////////////////////////////////////////////////////////////////
/////   Transistor noise contribution   ////////////////////////
////////////////////////////////////////////////////////////////


// thermal noise of res. elements

if (rscs>0) begin
    I(si,sx)    <+ white_noise(4*`P_K*Tdev_s/rscs_t*w*scalingFactor*pst*dens,"th_r_scs");
end

if (rdcs>0) begin
    I(di,dx)    <+ white_noise(4*`P_K*Tdev_s/rdcs_t*w*scalingFactor*pst*dens,"th_r_dcs");
end

if (rscm>0) begin
    I(sim,sx)   <+ white_noise(4*`P_K*Tdev_m/rscm_t*w*scalingFactor*(1-pst)*dens,"th_r_scm");
end

if (rdcm>0) begin
    I(dim,dx)   <+ white_noise(4*`P_K*Tdev_m/rdcm_t*w*scalingFactor*(1-pst)*dens,"th_r_dcm");
end

if (rsf>0) begin
    I(sx,s)     <+ white_noise(4*`P_K*$temperature/(rsf*w*scalingFactor),"th_rsf");
end

if (rdf>0) begin
    I(dx,d)     <+ white_noise(4*`P_K*$temperature/(rdf*w*scalingFactor),"th_rdf");
end

if (rg>0) begin
    I(gi,g)     <+ white_noise(4*`P_K*$temperature/(rg*w*scalingFactor),"th_rg");
end


// thermal noise semiconducting tubes
// description still missing (taking advantage of transmission probability)

// thermal noise metallic tubes
//V(dim,sim) <+ white_noise(4*`P_K*Tdev_m*rmt,"th_rmt");

// shot noise semiconducting tubes
I(di,si) <+ white_noise(2*`P_Q*curr*fanof,"sh_ch");

// shot noise of metallic tubes
//I(dim,sim) <+ white_noise(2*`P_Q*I(dim,sim)*fanofmt,"sh_mt");

// 1/f noise at channel semiconducting tubes

if ((Qs_cnt+Qd_cnt)>0) begin
    fnamp = hoogef/(abs((Qs_cnt+Qd_cnt)/`P_Q));
    I(di,si) <+ flicker_noise(fnamp*pow(abs(curr),2),beta_fn,"1of_ch");
end

// 1/f noise of metallic tubes
// description as for semiconducting tubes, but number of carriers required


////////////////////////////////////////////////////////////////
/////   Thermal network and self-heating   ////////////////////////
////////////////////////////////////////////////////////////////

// For simulators having problem with V(branch) <+ 0.0
// with external thermal node, 
if(flsh==0 || rths<`Rmin) begin
    I(br_sht_s)  <+ V(br_sht_s)/`Rmin;
end
else begin
    Pdiss_s     = abs(curr*V_disi);
    I(br_sht_s)  <+ V(br_sht_s)/rths - Pdiss_s;
end


if(flsh==0 || rthm<`Rmin) begin
    I(br_sht_m)  <+ V(br_sht_m)/`Rmin;
end
else begin
    if(rmt!= 0) begin
        Pdiss_m     = abs(V(dim,sim)*V(dim,sim)/rmt);
        I(br_sht_m)  <+ V(br_sht_m)/rthm - Pdiss_m;
    end
    else begin
        I(br_sht_m)  <+ V(br_sht_m)/`Rmin;
    end    
end


end

endmodule