// EPFL-EKV version 2.6: A Verilog-A description.
// The intrinsic device is coded according to the official manual
// (revision II) available at http://legwww.epfl.ch/ekv.
// contribution of Ivan Riis Nielsen 11/2006, modified by Dietmar Warning 01/2009

//Default simulator: Spectre

`ifdef insideADMS
  `define P(txt) (*txt*)
  `define PGIVEN(p)	$given(p)
  `define INITIAL_MODEL	@(initial_model)
  `define INSTANCE	@(initial_instance)
  `define NOISE 	@(noise)
`else
  `define P(txt) (txt)
  `define PGIVEN(p) 	p
  `define INITIAL_MODEL
  `define INSTANCE
  `define NOISE
`endif

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

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

`define NMOS 1
`define PMOS -1

`define EPSSI `P_EPS0*11.7
`define EPSOX `P_EPS0*3.9
`define TREF 300.15

`define SQR(x) ((x)*(x))

`define VT(temp) (`P_K*temp/`P_Q)
`define EG(temp) (1.16-0.000702*`SQR(temp)/(temp+1108))
`define NI(temp) (1.45e16*(temp/`TREF)*exp(`EG(`TREF)/(2*`VT(`TREF))-`EG(temp)/(2*`VT(temp))))


`define oneThird              3.3333333333333333e-01

// Constants needed in safe exponential function (called "expl")
`define se05                  2.3025850929940458e+02
`define ke05                  1.0e-100
`define ke05inv               1.0e100

//  P3       3rd order polynomial expansion of exp()
`define P3(u) (1.0 + (u) * (1.0 + 0.5 * ((u) * (1.0 + (u) * `oneThird))))

//  expl     exp() with 3rd order polynomial extrapolation
//           to avoid overflows and underflows and retain C-3 continuity
`define expl(x, res) \
if (abs(x) < `se05) begin\
    res       = exp(x); \
end else begin \
    if ((x) < -`se05) begin\
        res       = `ke05 / `P3(-`se05 - (x)); \
    end else begin\
        res       =  `ke05inv * `P3((x) - `se05); \
    end \
end


module ekv (d,g,s,b);

   // Node definitions

   inout d,g,s,b;
   electrical d,g,s,b,di,si;

   // Model parameters
   
   parameter integer nmos=1 from [0:1] `P(info="MOS type : nmos:0");
   parameter integer pmos=1 from [0:1] `P(info="MOS type : pmos:0");
   parameter integer MTYPE=(nmos==0 ? (pmos==0 ? 0 : 1) : (pmos==0 ? -1 : 1));
   parameter real TNOM=27 from (-273.15:inf) 
                  `P(info="Nominal temperature [degC]");
   parameter real IMAX=1 from (0:inf) 
                  `P(info="Maximum forward junction current before linearization [A]");
   
   // - intrinsic model (optional, section 4.2.1)
   parameter real TOX=0 from [0:inf) 
                  `P(info="Oxide thickness [m]");
   parameter real NSUB=0 from [0:inf)
                  `P(info="Channel doping [cm^-3]");
   parameter real VFB=1001.0 from (-inf:inf)        // use 1001V as "not specified"
                  `P(info="Flat-band voltage [V]");
   parameter real UO=0 from [0:inf) 
                  `P(info="Low-field mobility [cm^2/Vs]");
   parameter real VMAX=0 from [0:inf) 
                  `P(info="Saturation velocity [m/s]");
   parameter real THETA=0 from [0:inf) 
                  `P(info="Mobility reduction coefficient [V^-1]");
   
   // - intrinsic model (process related, section 4.1)
   parameter real COX=((TOX>0) ? (`EPSOX/TOX) : 0.7m) from [0:inf) 
                  `P(info="Oxide capacitance [F/m^2]");
   parameter real XJ=0.1u from [1n:inf) 
                  `P(info="Junction depth [m]");
   parameter real DL=0    from (-inf:inf) 
                  `P(info="Length correction [m]");
   parameter real DW=0    from (-inf:inf) 
                  `P(info="Width correction [m]");
   
   // - intrinsic model (basic, section 4.2)
   parameter real GAMMA=((NSUB>0) ? (sqrt(2*`P_Q*`EPSSI*NSUB*1e6)/COX) : 1) from [0:inf) 
                  `P(info="Body effect parameter [V^0.5]");
   parameter real PHI=((NSUB>0) ? (2*`VT((TNOM+273.15))*ln(max(NSUB,1)*1e6/`NI((TNOM+273.15)))) : 0.7) from [0.1:inf) 
                  `P(info="Bulk Fermi potential (*2) [V]");
   parameter real VTO=((VFB<1000.0) ? (VFB+MTYPE*(PHI+GAMMA*sqrt(PHI))) : 0.5) from (-inf:inf) 
                  `P(info="Long-channel threshold voltage [V]");
   parameter real KP=((UO>0) ? (UO*1e-4*COX) : 50u) from (0:inf) 
                  `P(info="Transconductance parameter [A/V^2]");
   parameter real UCRIT=(((VMAX>0) && (UO>0)) ? (VMAX/(UO*1e-4)) : 2e6 ) from [100k:inf) 
                  `P(info="Longitudinal critical field [V/m]");
   parameter real E0=((THETA>0) ? 0 : 1e12) from [100k:inf) 
                  `P(info="Mobility reduction coefficient [V/m]");
   
   // - intrinsic model (channel length modulation and charge sharing, section 4.3)
   parameter real LAMBDA=0.5 from [0:inf) 
                  `P(info="Depletion length coefficient (CLM)");
   parameter real WETA=0.25 from (-inf:inf) 
                  `P(info="Narrow-channel effect coefficient");
   parameter real LETA=0.1 from (-inf:inf) 
                  `P(info="Short-channel effect coefficient");
   
   // - intrinsic model (reverse short channel effect, section 4.4)
   parameter real Q0=0 from (-inf:inf) 
                  `P(info="RSCE peak charge density [C/m^2]");
   parameter real LK=0.29u from [10n:inf) 
                  `P(info="RSCE characteristic length [m]");
   
   // - intrinsic model (impact ionization, section 4.5)
   parameter real IBA=0 from (-inf:inf) 
                  `P(info="First impact ionization coefficient [m^-1]");
   parameter real IBB=3e8 from [1e8:inf) 
                  `P(info="Second impact ionization coefficient [V/m]");
   parameter real IBN=1 from [0.1:inf) 
                  `P(info="Saturation voltage factor for impact ionization");
   
   // - intrinsic model (temperature, section 4.6)
   parameter real TCV=1m from (-inf:inf) 
                  `P(info="Threshold voltage TC [V/K]");
   parameter real BEX=-1.5 from (-inf:inf)
                  `P(info="Mobility temperature exponent");
   parameter real UCEX=0.8 from (-inf:inf) 
                  `P(info="Longitudinal critical field temperature exponent");
   parameter real IBBT=9e-4 from (-inf:inf) 
                  `P(info="Temperature coefficient for IBB [K^-1]");
   
   // - intrinsic model (matching, section 4.7)
   parameter real AVTO=0 from (-inf:inf)
                  `P(info="Area related VTO mismatch parameter [Vm]");
   parameter real AKP=0 from (-inf:inf)
                  `P(info="Area related KP mismatch parameter [m]");
   parameter real AGAMMA=0 from (-inf:inf) 
                  `P(info="Area related GAMMA mismatch parameter [V^0.5*m]");
   
   // - intrinsic model (flicker noise, section 4.8)
   parameter real KF=0 from [0:inf)
                  `P(info="Flicker noise coefficient");
   parameter real AF=1 from (-inf:inf)
                  `P(info="Flicker noise exponent");
   
   // - intrinsic model (setup, section 4.9)
   parameter real NQS=0 from [0:1]
                  `P(info="Non-quasi-static operation switch");
   parameter real SATLIM=exp(4) from (0:inf) 
                  `P(info="Saturation limit (if/ir)");
   parameter real XQC=0.4 from [0:1]
                  `P(info="Charge/capacitance model selector");
   
   // - external parasitic parameters
   parameter real HDIF=0 from [0:inf) 
                  `P(info="S/D diffusion length (/2) [m]");
   parameter real RSH=0 from [0:inf) 
                  `P(info="S/D sheet resistance [ohm]");
   parameter real JS=0 from [0:inf) 
                  `P(info="S/D junction saturation current density [A/m^2]");
   parameter real JSW=0 from [0:inf) 
                  `P(info="S/D junction sidewall saturation current density [A/m]");
   parameter real XTI=0 from [0:inf) 
                  `P(info="S/D diode saturation current temperature exponent");
   parameter real N=1 from [0.5:10]
                  `P(info="S/D diode emission coefficient");
   parameter real CJ=0 from [0:inf)
                  `P(info="S/D zero-bias junction capacitance per area [F/m^2]");
   parameter real CJSW=0 from [0:inf)
                  `P(info="S/D zero-bias junction capacitance per perimeter [F/m]");
   parameter real PB=0.8 from (0:inf)
                  `P(info="S/D bottom junction builtin potential [V]");
   parameter real PBSW=PB from (0:inf)
                  `P(info="S/D sidewall junction builtin potential [V]");
   parameter real MJ=0.5 from (0:inf)
                  `P(info="S/D bottom junction grading coefficient");
   parameter real MJSW=0.333 from (0:inf)
                  `P(info="S/D sidewall junction grading coefficient");
   parameter real FC=0.5 from (0:inf)
                  `P(info="S/D bottom junction forward-bias threshold");
   parameter real FCSW=FC from (0:inf)
                  `P(info="S/D sidewall junction forward-bias threshold");
   parameter real CGSO=0 from [0:inf)
                  `P(info="Gate-source overlap capacitance per width [F/m]");
   parameter real CGDO=0 from [0:inf)
                  `P(info="Gate-drain overlap capacitance per width [F/m]");
   parameter real CGBO=0 from [0:inf)
                  `P(info="Gate-bulk overlap capacitance per length [F/m]");


   // Instance parameters

   // - intrinsic model
   parameter real L=10u from [0:inf] 
                  `P(type="instance" info="Drawn length [m]" unit="m");
   parameter real W=10u from [0:inf] 
                  `P(type="instance" info="Drawn width [m]" unit="m");
   parameter real M=1   from [0:inf] 
                  `P(type="instance" info="Parallel multiplier" unit="m");
//   parameter real N=1   from [0:inf] 
//                  `P(type="instance" info="Series multiplier" unit="m");

   // - external parasitics
   parameter real AD=((HDIF>0) ? (2*HDIF*W) : 0) from [0:inf) 
                  `P(info="Drain area [m^2]" type="instance");
   parameter real AS=((HDIF>0) ? (2*HDIF*W) : 0) from [0:inf) 
                  `P(info="Source area [m^2]" type="instance");
   parameter real PD=((HDIF>0) ? (4*HDIF+2*W) : 0) from [0:inf) 
                  `P(info="Drain perimeter [m]" type="instance");
   parameter real PS=((HDIF>0) ? (4*HDIF+2*W) : 0) from [0:inf) 
                  `P(info="Source perimeter [m]" type="instance");
   parameter real NRD=((HDIF>0) ? (HDIF/W) : 0) from [0:inf)
                  `P(info="Drain no. squares" type="instance");
   parameter real NRS=((HDIF>0) ? (HDIF/W) : 0) from [0:inf)
                  `P(info="Source no. squares" type="instance");
   parameter real RS=((RSH>0) ? (RSH*NRS) : 0) from [0:inf)
                  `P(info="Source resistance [ohms]" type="instance");
   parameter real RD=((RSH>0) ? (RSH*NRD) : 0) from [0:inf)
                  `P(info="Drain resistance [ohms]" type="instance");


   // Declaration of variables
   integer  mode;
   real     lc,isat_s,vexp_s,gexp_s,isat_d,vexp_d,gexp_d,fact,
            weff,leff,np,ns,lmin,rd,rs,ceps,ca,xsi,dvrsce,
            tempk,vt,sqrt_A,vto_a,kp_a,gamma_a,ucrit,phi,ibb,vc,qb0,
            vg,vd,vs,tmp,vgprime,vp0,vsprime,vdprime,gamma0,gammaprime,vp,n,ifwd,
            vdss,vdssprime,dv,vds,vip,dl,lprime,leq,irprime,irev,beta0,nau,
            nq,xf,xr,qd,qs,qi,qb,qg,beta0prime,beta,vpprime,is,ids,vib,
            idb,ibdj,ibsj,coxt,qdt,qst,qgt,qbt,
            cbs0,cbs0sw,cbs,cbd0,cbd0sw,cbd,
            fv,z0,z1,y;

   real     cgso,cgdo,cgbo;
   

   analog begin

      `INITIAL_MODEL begin      // Model Initialization

         lc = sqrt(`EPSSI/COX*XJ);
         
      end // INITIAL_MODEL

      `INSTANCE begin // temperature independent device initialization

         weff = W+DW;
         leff = L+DL;

         np = M;
         ns = 1;

         // eq. 54
         lmin = 0.1*ns*leff;
         
         rs = ns/np*RS;
         rd = ns/np*RD;

         ceps = 4*22e-3*22e-3;
         ca = 0.028;
         xsi = ca*(10*leff/LK-1);
         dvrsce = 2*Q0/COX/`SQR(1+0.5*(xsi+sqrt(xsi*xsi+ceps)));

         coxt = np*ns*COX*weff*leff;
         
      end // temperature independent
      
      `INSTANCE begin // temperature dependent device initialization
         tempk = $temperature;
         vt = `VT(tempk);
         
         sqrt_A = sqrt(np*weff*ns*leff);
         
         vto_a   = MTYPE*(VTO+TCV*(tempk-(TNOM+273.15)))+AVTO/sqrt_A;
         kp_a    = KP*pow(tempk/(TNOM+273.15),BEX)*(1+AKP/sqrt_A);
         gamma_a = GAMMA+AGAMMA/sqrt_A;
         ucrit   = UCRIT*pow(tempk/(TNOM+273.15),UCEX);
         phi     = PHI*tempk/(TNOM+273.15)-3*vt*ln(tempk/(TNOM+273.15))-`EG(TNOM+273.15)*tempk/(TNOM+273.15)+`EG(tempk);
         ibb     = IBB*(1+IBBT*(tempk-(TNOM+273.15)));

         vc = ucrit*ns*leff;
         
         // eq. 60
         qb0 = gamma_a*sqrt(phi);

         fact = (`EG(TNOM+273.15)/`VT(TNOM+273.15)-`EG(tempk)/vt) * pow(tempk/(TNOM+273.15),XTI);
         `expl(fact,tmp)
         isat_s = np*ns*(JS*AS+JSW*PS)*tmp;
         isat_d = np*ns*(JS*AD+JSW*PD)*tmp;

         if (isat_s>0) begin
            vexp_s = vt*ln(IMAX/isat_s+1);
            gexp_s = (IMAX+isat_s)/vt;
         end else begin
            vexp_s = -1e9;
            gexp_s = 0;
         end

         if (isat_d>0) begin
            vexp_d = vt*ln(IMAX/isat_d+1);
            gexp_d = (IMAX+isat_d)/vt;
         end else begin
            vexp_d = -1e9;
            gexp_d = 0;
         end
         
         cbs0 = np*ns*CJ*AS;
         cbd0 = np*ns*CJ*AD;
         cbs0sw = np*ns*CJSW*PS;
         cbd0sw = np*ns*CJSW*PD;

         cgso = np*ns*CGSO*weff;
         cgdo = np*ns*CGDO*weff;
         cgbo = np*ns*CGBO*leff;
         
      end // temperature dependent

      
      begin     //Bias-dependent model evaluation

         vg = MTYPE*V(g,b);
         vd = MTYPE*V(di,b);
         vs = MTYPE*V(si,b);
         // $strobe("vg=%e vd=%e vs=%e",vg,vd,vs);
         
         if (vd>=vs)
            mode = 1;
         else begin
            mode = -1;
            tmp = vs;
            vs = vd;
            vd = tmp;
         end

         // eq. 33
         vgprime = vg-vto_a-dvrsce+phi+gamma_a*sqrt(phi);
         // eq. 35
         vsprime = 0.5*(vs+phi+sqrt(`SQR(vs+phi)+16*`SQR(vt)));
         vdprime = 0.5*(vd+phi+sqrt(`SQR(vd+phi)+16*`SQR(vt)));
         // $strobe("vgprime=%e vdprime=%e vsprime=%e",vgprime,vdprime,vsprime);
         // eq. 34
         if (vgprime>=0) begin
            vp0 = vgprime-phi-gamma_a*(sqrt(vgprime+0.25*`SQR(gamma_a))-0.5*gamma_a);
            // eq. 36
            gamma0 = gamma_a-`EPSSI/COX*(LETA/leff*(sqrt(vsprime)+sqrt(vdprime))-3*WETA/weff*sqrt(vp0+phi));
         end else begin
            vp0 = -phi;
            // eq. 36 - skipped sqrt(vp0+phi) here, it produces inf on derivative
            gamma0 = gamma_a-`EPSSI/COX*(LETA/leff*(sqrt(vsprime)+sqrt(vdprime))  );
         end
         // eq. 37
         gammaprime = 0.5*(gamma0+sqrt(`SQR(gamma0)+0.1*vt));
         // eq. 38
         if (vgprime>=0)
           vp = vgprime-phi-gammaprime*(sqrt(vgprime+0.25*`SQR(gammaprime))-0.5*gammaprime);
         else
           vp = -phi;
         // $strobe("vp0=%e vp=%e gamma0=%e gammaprime=%e",vp0,vp,gamma0,gammaprime);
         // eq. 39
         n = 1+gamma_a*0.5/sqrt(vp+phi+4*vt);

	 // Forward  current (43-44)
   	 fv=(vp-vs)/vt;
	 
     	 if (fv > -0.35) begin
           z0 = 2.0/(1.3 + fv - ln(fv+1.6));
    	   z1 = (2.0 + z0) / (1.0 + fv + ln(z0));
	   y = (1.0 + fv + ln(z1)) / (2.0 + z1);
         end
    	 else if (fv > -15) begin
	   `expl(-fv,tmp)
           z0 = 1.55 + tmp;
    	   z1 = (2.0 + z0) / (1.0 + fv + ln(z0));
	   y = (1.0 + fv + ln(z1)) / (2.0 + z1);
         end 
         else if (fv > -23.0) begin
	   `expl(-fv,tmp)
	   y = 1.0 / (2.0 + tmp);
         end
   	 else begin
 	   `expl(fv,tmp)
	   y = tmp + 1.0e-64;
	 end
	 
   	 ifwd = y*(1.0 + y);
   	 z0  = 1;
   	 z1  = 1;
         
         // eq. 46
         vdss = vc*(sqrt(0.25+vt/vc*sqrt(ifwd))-0.5);
         // eq. 47
         vdssprime = vc*(sqrt(0.25+vt/vc*(sqrt(ifwd)-0.75*ln(ifwd)))-0.5)+vt*(ln(0.5*vc/vt)-0.6);
         // $strobe("ifwd=%e vdss=%e vdssprime=%e",ifwd,vdss,vdssprime);
         // eq. 48
         dv = 4*vt*sqrt(LAMBDA*(sqrt(ifwd)-vdss/vt)+1.0/64);
         // eq. 49
         vds = 0.5*(vd-vs);
         // eq. 50
         vip = sqrt(`SQR(vdss)+`SQR(dv))-sqrt(`SQR(vds-vdss)+`SQR(dv));
         // eq. 52
         dl = LAMBDA*lc*ln(1+(vds-vip)/(lc*ucrit));

         // eq. 53
         lprime = ns*leff-dl+(vds+vip)/ucrit;
         // eq. 55
         leq = 0.5*(lprime+sqrt(`SQR(lprime)+`SQR(lmin)));

         // eq. 56
   	 fv=(vp-vds-vs-sqrt(`SQR(vdssprime)+`SQR(dv))+sqrt(`SQR(vds-vdssprime)+`SQR(dv)))/vt;
	 
     	 if (fv > -0.35) begin
           z0 = 2.0/(1.3 + fv - ln(fv+1.6));
    	   z1 = (2.0 + z0) / (1.0 + fv + ln(z0));
	   y = (1.0 + fv + ln(z1)) / (2.0 + z1);
         end
    	 else if (fv > -15) begin
	   `expl(-fv,tmp)
           z0 = 1.55 + tmp;
    	   z1 = (2.0 + z0) / (1.0 + fv + ln(z0));
	   y = (1.0 + fv + ln(z1)) / (2.0 + z1);
         end 
         else if (fv > -23.0) begin
	   `expl(-fv,tmp)
	   y = 1.0 / (2.0 + tmp);
         end
   	 else begin
 	   `expl(fv,tmp)
	   y = tmp + 1.0e-64;
	 end
	 
   	 irprime = y*(1.0 + y);
   	 z0  = 1;
   	 z1  = 1;

         // eq. 57
   	 fv=(vp-vd)/vt;
	 
     	 if (fv > -0.35) begin
           z0 = 2.0/(1.3 + fv - ln(fv+1.6));
    	   z1 = (2.0 + z0) / (1.0 + fv + ln(z0));
	   y = (1.0 + fv + ln(z1)) / (2.0 + z1);
         end
    	 else if (fv > -15) begin
	   `expl(-fv,tmp)
           z0 = 1.55 + tmp;
    	   z1 = (2.0 + z0) / (1.0 + fv + ln(z0));
	   y = (1.0 + fv + ln(z1)) / (2.0 + z1);
         end 
         else if (fv > -23.0) begin
	   `expl(-fv,tmp)
	   y = 1.0 / (2.0 + tmp);
         end
   	 else begin
 	   `expl(fv,tmp)
	   y = tmp + 1.0e-64;
	 end
	 
   	 irev = y*(1.0 + y);

         // eq. 58
         beta0 = kp_a*np*weff/leq;
         // eq. 59
         nau = (5+MTYPE)/12.0;

         // eq. 69
         nq = 1+0.5*gamma_a/sqrt(vp+phi+1e-6);
         
         // eq. 70
         xf = sqrt(0.25+ifwd);
         // eq. 71
         xr = sqrt(0.25+irev);
         // eq. 72
         qd = -nq*(4.0/15*(3*`SQR(xr)*(xr+2*xf)+2*`SQR(xf)*(xf+2*xr))/`SQR(xf+xr)-0.5);
         // eq. 73
         qs = -nq*(4.0/15*(3*`SQR(xf)*(xf+2*xr)+2*`SQR(xr)*(xr+2*xf))/`SQR(xf+xr)-0.5);
         // eq. 74
         qi = qs+qd;
         // eq. 75
         if (vgprime>=0)
           qb = (-gamma_a*sqrt(vp+phi+1e-6))/vt-(nq-1)/nq*qi;
         else
           qb = -vgprime/vt;
         // eq. 76 (qox removed since it is assumed to be zero)
         qg = -qi-qb;

         if (E0!=0) begin 
           // eq. 61
           beta0prime = beta0*(1+COX/(E0*`EPSSI)*qb0);
           // eq. 62
           beta = beta0prime/(1+COX/(E0*`EPSSI)*vt*abs(qb+nau*qi));
         end else begin
            // eq. 63
            vpprime = 0.5*(vp+sqrt(`SQR(vp)+2*`SQR(vt)));
            // eq. 64
            beta = beta0/(1+THETA*vpprime);
         end // else: !if(e0!=0)
         // eq. 65
         is = 2*n*beta*`SQR(vt);
         // $strobe("beta0=%e beta0prime=%e beta=%e E0=%e qb0=%e qb=%e qi=%e",beta0,beta0prime,beta,E0,qb0,qb,qi);
         // eq. 66
         ids = is*(ifwd-irprime);
         // eq. 67
         vib = vd-vs-IBN*2*vdss;
         // eq. 68
         if (vib>0) begin
           `expl((-ibb*lc)/vib,tmp)
           idb = ids*IBA/ibb*vib*tmp;
         end else
           idb = 0;
         // $strobe("ids=%e idb=%e",ids,idb);

         if (mode>1) begin
            if (isat_s>0) begin
               if (-vs>vexp_s)
                 ibsj = IMAX+gexp_s*(-vs-vexp_s);
               else begin
                 `expl(-vs/(N*vt),tmp)
                 ibsj = isat_s*(tmp-1);
               end
            end else
              ibsj = 0;
            
            if (isat_d>0) begin
               if (-vd>vexp_d)
                 ibdj = IMAX+gexp_d*(-vd-vexp_d);
               else begin
                 `expl(-vd/(N*vt),tmp)
                 ibdj = isat_d*(tmp-1);
               end
            end else
              ibdj = 0;
            
         end else begin // if (mode>1)
            if (isat_s>0) begin
               if (-vd>vexp_s)
                 ibsj = IMAX+gexp_s*(-vd-vexp_s);
               else begin
                 `expl(-vd/(N*vt),tmp)
                 ibsj = isat_s*(tmp-1);
               end
            end else
              ibsj = 0;
            
            if (isat_d>0) begin
               if (-vs>vexp_d)
                 ibdj = IMAX+gexp_d*(-vs-vexp_d);
               else begin
                 `expl(-vs/(N*vt),tmp)
                 ibdj = isat_d*(tmp-1);
               end
            end else
              ibdj = 0;
            
         end // else: !if(mode>1)

         qdt = coxt*vt*qd;
         qst = coxt*vt*qs;
         qgt = coxt*vt*qg;
         qbt = coxt*vt*qb;

         cbs = 0;
         cbd = 0;
         
         if (cbs0>0) begin
            if (MTYPE*V(b,si)>FC*PB)
              cbs = cbs+cbs0/pow(1-FC,MJ)*(1+MJ*(MTYPE*V(b,si)-PB*FC))/(PB*(1-FC));
            else
              cbs = cbs+cbs0/pow(1-MTYPE*V(b,si),MJ);
         end
         if (cbd0>0) begin
            if (MTYPE*V(b,di)>FC*PB)
              cbd = cbd+cbd0/pow(1-FC,MJ)*(1+MJ*(MTYPE*V(b,di)-PB*FC))/(PB*(1-FC));
            else
              cbd = cbd+cbd0/pow(1-MTYPE*V(b,di),MJ);
         end
         if (cbs0sw>0) begin
            if (MTYPE*V(b,si)>FCSW*PBSW)
              cbs = cbs+cbs0sw/pow(1-FCSW,MJSW)*(1+MJSW*(MTYPE*V(b,si)-PBSW*FCSW))/(PBSW*(1-FCSW));
            else
              cbs = cbs+cbs0sw/pow(1-MTYPE*V(b,si),MJSW);
         end
         if (cbd0sw>0) begin
            if (MTYPE*V(b,di)>FCSW*PBSW)
              cbd = cbd+cbd0sw/pow(1-FCSW,MJSW)*(1+MJSW*(MTYPE*V(b,di)-PBSW*FCSW))/(PBSW*(1-FCSW));
            else
              cbd = cbd+cbd0sw/pow(1-MTYPE*V(b,di),MJSW);
         end
         
      end       //Bias-dependent model evaluation

      begin     //Define branch sources

         I(di,si) <+ MTYPE*mode*ids;
         if (mode>0) begin
            I(di,b) <+ MTYPE*idb;

            I(di,g) <+ MTYPE*ddt(qdt);
            I(si,g) <+ MTYPE*ddt(qst);
            
         end else begin
            I(si,b) <+ MTYPE*idb;

            I(si,g) <+ MTYPE*ddt(qdt);
            I(di,g) <+ MTYPE*ddt(qst);
            
         end // else: !if(mode>0)

         I(b,si) <+ MTYPE*ibsj;
         I(b,di) <+ MTYPE*ibdj;
         
         I(b,g)  <+ MTYPE*ddt(qbt);
         
         I(g,si) <+ cgso*ddt(V(g,si));
         I(g,di) <+ cgdo*ddt(V(g,di));
         I(g,b)  <+ cgbo*ddt(V(g,b));
         
         if (RD>0)
           I(d,di) <+ V(d,di)/rd;
         else   
           V(d,di) <+ 0.0;
         if (RS>0)
           I(s,si) <+ V(s,si)/rs;
         else
           V(s,si) <+ 0.0;

         I(b,si) <+ cbs*ddt(V(b,si));
         I(b,di) <+ cbd*ddt(V(b,di));
         
      end // begin

//      `NOISE begin    //Define noise sources
//       
//      end // noise
      
   end //analog

endmodule
