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

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

`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 MIN_R 0.001
`define VEXLIM  200.0

`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 y_fv(fv,y)\
    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 \
        tmp = exp(-fv);\
        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 \
        tmp = exp(-fv);\
        y = 1.0 / (2.0 + tmp);\
    end \
    else begin \
        tmp = exp(fv);\
        y = tmp + 1.0e-64;\
    end

`define expLin(result, x)\
    if (x < `VEXLIM)\
        result = exp(x);\
    else begin\
        result = exp(`VEXLIM) * (1.0 + (x - `VEXLIM));\
    end

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

    // Node definitions

    inout           d,g,s,b;     // external nodes
    electrical      d,g,s,b;     // external nodes
    electrical      dp,sp;       // internal nodes

    branch (dp,sp) dpsp;
    branch (dp,b)  dpb;
    branch (sp,b)  spb;
    branch (g,dp)  gpdp;
    branch (g,sp)  gpsp;
    branch (g,b)   gb;
    branch (d,dp)  ddp;
    branch (s,sp)  ssp;

    // Instance parameters

    // - intrinsic model
    parameter real l=10e-6 from [0:inf) `P(type="instance" desc="Drawn length [m]" units="m");
    parameter real w=10e-6 from [0:inf) `P(type="instance" desc="Drawn width [m]" units="m");
    parameter real m=1.0   from [1:inf) `P(type="instance" desc="Parallel multiplier" units="m");
    parameter real ns=1.0  from [1:inf) `P(type="instance" desc="Series multiplier" units="m");
    parameter real dtemp = 0.0 from (-inf:inf) `P(type="instance" desc="Difference sim. temp and device temp" units="C");

    // - external parasitics
    parameter real ad=0.0 from [0:inf) `P(type="instance" desc="Drain area" units="m^2");
    parameter real as=0.0 from [0:inf) `P(type="instance" desc="Source area" units="m^2");
    parameter real pd=0.0 from [0:inf) `P(type="instance" desc="Drain perimeter" units="m");
    parameter real ps=0.0 from [0:inf) `P(type="instance" desc="Source perimeter" units="m");
    parameter real nrd=0.0 from [0:inf) `P(type="instance" desc="Drain no. squares");
    parameter real nrs=0.0 from [0:inf) `P(type="instance" desc="Source no. squares");

    // Model parameters

    parameter integer nmos=1 from [0:1] `P(desc="MOS channel type");
    parameter integer pmos=1 from [0:1] `P(desc="MOS channel type");
    parameter integer type=1 from [-1:1] exclude 0;
    parameter real tnom=27.0 from [-273.15:inf) `P(desc="Nominal temperature" units="C");
    parameter real imax=1.0 from [0:inf) `P(desc="Maximum forward junction current before linearization" units="A");

    // - intrinsic model (optional, section 4.2.1)
    parameter real tox=0.0 from [0:inf) `P(desc="Oxide thickness" units="m");
    parameter real nsub=0.0 from [0:inf) `P(desc="Channel doping" units="cm^-3");
    parameter real vfb=-1.0 from (-inf:inf) `P(desc="Flat-band voltage" units="V");

    parameter real uo=0.0 from [0:inf) `P(desc="Low-field mobility" units="cm^2/Vs");
    parameter real vmax=0.0 from [0:inf) `P(desc="Saturation velocity" units="m/s");
    parameter real theta=0.0 from [0:inf) `P(desc="Mobility reduction coefficient" units="V^-1");

    // - intrinsic model (process related, section 4.1)
    parameter real cox=0.7e-3 from [1e-6:inf) `P(desc="Oxide capacitance" units="F/m^2");
    parameter real xj=0.1e-6 from [1n:inf) `P(desc="Junction depth" units="m");
    parameter real dl=0.0    from (-inf:inf) `P(desc="Length correction" units="m");
    parameter real dw=0.0    from (-inf:inf) `P(desc="Width correction" units="m");

    // - intrinsic model (basic, section 4.2)
    parameter real gamma=0.7 from [0:inf) `P(desc="Body effect parameter" units="V^0.5");
    parameter real phi=0.5 from [0.1:inf) `P(desc="Bulk Fermi potential (*2)" units="V");
    parameter real vto=0.5 from (-inf:inf) `P(desc="Long-channel threshold voltage" units="V");
    parameter real kp=20e-6 from [0.0:inf) `P(desc="Transconductance parameter" units="A/V^2");
    parameter real ucrit=1.0e+6 from [100e+3:inf) `P(desc="Longitudinal critical field" units="V/m");
    parameter real e0=1.0e-9 from [1e-12:inf) `P(desc="Mobility reduction coefficient" units="V/m");

    // - intrinsic model (channel length modulation and charge sharing, section 4.3)
    parameter real lambda=0.5 from [0:inf) `P(desc="Depletion length coefficient (CLM)");
    parameter real weta=0.25 from (-inf:inf) `P(desc="Narrow-channel effect coefficient");
    parameter real leta=0.1 from (-inf:inf) `P(desc="Short-channel effect coefficient");

    // - intrinsic model (reverse short channel effect, section 4.4)
    parameter real q0=0.0 from (-inf:inf) `P(desc="RSCE peak charge density" units="C/m^2");
    parameter real lk=0.29e-6 from [10n:inf) `P(desc="RSCE characteristic length" units="m");

    // - intrinsic model (impact ionization, section 4.5)
    parameter real iba=0.0 from (-inf:inf) `P(desc="First impact ionization coefficient" units="m^-1");
    parameter real ibb=3e8 from [1e8:inf) `P(desc="Second impact ionization coefficient" units="V/m");
    parameter real ibn=1.0 from [0.1:inf) `P(desc="Saturation voltage factor for impact ionization");

    // - intrinsic model (temperature, section 4.6)
    parameter real tcv=1e-3 from (-inf:inf) `P(desc="Threshold voltage TC" units="V/deg");
    parameter real bex=-1.5 from (-inf:inf) `P(desc="Mobility temperature exponent");
    parameter real ucex=0.8 from (-inf:inf) `P(desc="Longitudinal critical field temperature exponent");
    parameter real ibbt=9e-4 from (-inf:inf) `P(desc="Temperature coefficient for ibb" units="K^-1");

    // - intrinsic model (matching, section 4.7)
    parameter real avto=0.0 from (-inf:inf) `P(desc="Area related vto mismatch parameter" units="Vm");
    parameter real akp=0.0 from (-inf:inf) `P(desc="Area related kp mismatch parameter" units="m");
    parameter real agamma=0.0 from (-inf:inf) `P(desc="Area related gamma mismatch parameter" units="V^0.5*m");

    // - intrinsic model (flicker noise, section 4.8)
    parameter real kf=0.0 from [0:inf) `P(desc="Flicker noise coefficient");
    parameter real af=1.0 from (-inf:inf) `P(desc="Flicker noise exponent");

    // - external parasitic parameters
    parameter real hdif=0.0 from [0:inf) `P(desc="S/D diffusion length (/2)" units="m");
    parameter real rsh=0.0 from [0:inf) `P(desc="S/D sheet resistance" units="Ohm");
    parameter real js=0.0 from [0:inf) `P(desc="S/D junction saturation current density" units="A/m^2");
    parameter real jsw=0.0 from [0:inf) `P(desc="S/D junction sidewall saturation current density" units="A/m");
    parameter real xti=0.0 from [0:inf) `P(desc="S/D diode saturation current temperature exponent");
    parameter real n=1 from [0.5:10] `P(desc="S/D diode emission coefficient");
    parameter real cj=0.0 from [0:inf) `P(desc="S/D zero-bias junction capacitance per area" units="F/m^2");
    parameter real cjsw=0.0 from [0:inf) `P(desc="S/D zero-bias junction capacitance per perimeter" units="F/m");
    parameter real pb=0.8 from (0:inf) `P(desc="S/D bottom junction builtin potential" units="V");
    parameter real pbsw=pb from (0:inf) `P(desc="S/D sidewall junction builtin potential" units="V");
    parameter real mj=0.5 from (0:inf) `P(desc="S/D bottom junction grading coefficient");
    parameter real mjsw=0.333 from (0:inf) `P(desc="S/D sidewall junction grading coefficient");
    parameter real fc=0.5 from (0:inf) `P(desc="S/D bottom junction forward-bias threshold");
    parameter real fcsw=fc from (0:inf) `P(desc="S/D sidewall junction forward-bias threshold");
    parameter real cgso=1.5e-10 from [0:inf) `P(desc="Gate-source overlap capacitance per width" units="F/m");
    parameter real cgdo=1.5e-10 from [0:inf) `P(desc="Gate-drain overlap capacitance per width" units="F/m");
    parameter real cgbo=4.0e-10 from [0:inf) `P(desc="Gate-bulk overlap capacitance per length" units="F/m");

    // Declaration of variables
    integer mode, MOStype;
    real lc,isat_s,vexp_s,gexp_s,isat_d,vexp_d,gexp_d,fact,
    weff,leff,lmin,RDeff,RSeff,ceps,ca,xsi,dvrsce,
    TempK,Vt,sqrt_A,vto_a,kp_a,gamma_a,ucrit_a,phi_a,ibb_a,vc,qb0,
    vg,vd,vs,tmp,vgprime,vp0,vsprime,vdprime,gamma0,gammaprime,vp,nslope,ifwd,
    vdss,vdssprime,dv,vds,vip,dl_a,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,qdtx,qstx,qgt,qjs,qjd,
    cbs0,cbs0sw,cbs,cbd0,cbd0sw,cbd,v_bp_dp, v_bp_sp,
    fv,z0,z1,y;

    real ADeff, ASeff, PDeff, PSeff;
    real cgso_s,cgdo_s,cgbo_s;
    real gmin, TnomK;
    real cox_p, gamma_p, phi_p, kp_p, vto_p, ucrit_p;

    analog begin

        MOStype = 1;

        gmin = $simparam("gmin");

        `INITIAL_MODEL     // Model Initialization
        begin
            if (`PGIVEN(nmos)) begin
                MOStype = `NMOS;
            end else if (`PGIVEN(pmos)) begin
                MOStype = `PMOS;
            end else begin
                MOStype = (`PGIVEN(type)) ? type : `NMOS;
            end
            //$strobe("MOStype %d", MOStype);

            if (`PGIVEN(cox)) begin
                cox_p = cox;
            end else begin
                cox_p = (tox>0) ? (`EPSOX/tox) : 0.7e-3;
            end

            if (`PGIVEN(gamma)) begin
                gamma_p = gamma;
            end else begin
                gamma_p = (nsub>0) ? (sqrt(2*`P_Q*`EPSSI*nsub*1e+6)/cox_p) : 0.7;
            end

            if (`PGIVEN(phi)) begin
                phi_p = phi;
            end else begin
                phi_p = (nsub>0) ? (2*`VT(tnom+273.15)*ln(nsub*1e+6/`NI(tnom+273.15))) : 0.5;
            end

            if (`PGIVEN(kp)) begin
                kp_p = kp;
            end else begin
                kp_p = (uo>0) ? (uo*1e-4*cox_p) : 20e-6;
            end

            if (`PGIVEN(vto)) begin
                vto_p = vto;
            end else begin
                vto_p = (`PGIVEN(vfb)) ? (vfb+phi_p+gamma_p*sqrt(phi_p)) : 0.5;
            end

            if (`PGIVEN(ucrit)) begin
                ucrit_p = ucrit;
            end else begin
                ucrit_p = (vmax>0 && uo>0) ? (vmax/(uo*1e-4)) : 100e+3;
            end

            lc = sqrt(`EPSSI/cox_p*xj);

        end // INITIAL_MODEL

        `INSTANCE // temperature independent device initialization
        begin
            weff = w+dw;
            leff = l+dl;

            // eq. 54
            lmin = 0.1*ns*leff;

            if (hdif > 0) begin
                RSeff = ns/m*rsh*hdif/weff;
                RDeff = ns/m*rsh*hdif/weff;
                ADeff = 2*hdif*weff;
                ASeff = 2*hdif*weff;
                PDeff = 4*hdif+2*weff;
                PSeff = 4*hdif+2*weff;
            end else begin
                RSeff = ns/m*rsh*nrs;
                RDeff = ns/m*rsh*nrd;
                ADeff = ad;
                ASeff = as;
                PDeff = pd;
                PSeff = ps;
            end
            if (RDeff < `MIN_R) begin
                RDeff = `MIN_R;
            end
            if (RSeff < `MIN_R) begin
                RSeff = `MIN_R;
            end

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

            coxt = m*ns*cox_p*weff*leff;

            cbs0 = m*ns*cj*ASeff;
            cbd0 = m*ns*cj*ADeff;
            cbs0sw = m*ns*cjsw*PSeff;
            cbd0sw = m*ns*cjsw*PDeff;

            cgso_s = m*ns*cgso*weff;
            cgdo_s = m*ns*cgdo*weff;
            cgbo_s = m*ns*cgbo*leff;

        end // temperature independent

        `INSTANCE // temperature dependent device initialization
        begin
            if (dtemp > 0.0) begin
                TempK = $temperature + dtemp;
            end else begin
                TempK = $temperature;
            end

            TnomK = tnom + 273.15;

            Vt = `VT(TempK);

            sqrt_A = sqrt(m*weff*ns*leff);

            vto_a   = MOStype*(vto_p+tcv*(TempK-TnomK))+avto/sqrt_A;
            kp_a    = m*kp_p*pow(TempK/TnomK,bex)*(1+akp/sqrt_A);
            gamma_a = gamma_p+agamma/sqrt_A;
            ucrit_a = ucrit_p*pow(TempK/TnomK,ucex);
            phi_a   = phi_p*TempK/TnomK-3*Vt*ln(TempK/TnomK)-`EG(TnomK)*TempK/TnomK+`EG(TempK);
            ibb_a   = ibb*(1+ibbt*(TempK-TnomK));

            vc = ucrit_a*ns*leff;

            // eq. 60
            qb0 = gamma_a*sqrt(phi_a);

            fact = (`EG(TnomK)/`VT(TnomK)-`EG(TempK)/Vt) * pow(TempK/TnomK,xti);
            `expLin(tmp,fact)
            isat_s = m*ns*(js*ASeff+jsw*PSeff)*tmp;
            isat_d = m*ns*(js*ADeff+jsw*PDeff)*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

        end // temperature dependent


        begin     //Bias-dependent model evaluation

            vg = MOStype*V(gb);
            vd = MOStype*V(dpb);
            vs = MOStype*V(spb);
            // $strobe("MOStype %d vg=%e vd=%e vs=%e",MOStype,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_a+gamma_a*sqrt(phi_a);
            // eq. 35
            vsprime = 0.5*(vs+phi_a+sqrt(`SQR(vs+phi_a)+16*`SQR(Vt)));
            vdprime = 0.5*(vd+phi_a+sqrt(`SQR(vd+phi_a)+16*`SQR(Vt)));
            // $strobe("vgprime=%e vdprime=%e vsprime=%e",vgprime,vdprime,vsprime);
            // eq. 34
            if (vgprime>=0) begin
                vp0 = vgprime-phi_a-gamma_a*(sqrt(vgprime+0.25*`SQR(gamma_a))-0.5*gamma_a);
                // eq. 36
                gamma0 = gamma_a-`EPSSI/cox_p*(leta/leff*(sqrt(vsprime)+sqrt(vdprime))-3*weta/weff*sqrt(vp0+phi_a));
            end else begin
                vp0 = -phi_a;
                // eq. 36 - skipped sqrt(vp0+phi_a) here, it produces inf on derivative
                gamma0 = gamma_a-`EPSSI/cox_p*(leta/leff*(sqrt(vsprime)+sqrt(vdprime))  );
            end
            // eq. 37
            gammaprime = 0.5*(gamma0+sqrt(`SQR(gamma0)+0.1*Vt));
            // eq. 38
            if (vgprime>=0) begin
                vp = vgprime-phi_a-gammaprime*(sqrt(vgprime+0.25*`SQR(gammaprime))-0.5*gammaprime);
            end else begin
                vp = -phi_a;
            end
            // $strobe("vp0=%e vp=%e gamma0=%e gammaprime=%e",vp0,vp,gamma0,gammaprime);
            // eq. 39
            nslope = 1+gamma_a*0.5/sqrt(vp+phi_a+4*Vt);

            // Forward  current (43-44)
            fv=(vp-vs)/Vt;

            `y_fv(fv,y)

            ifwd = y*(1.0 + y);

            // 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_a = lambda*lc*ln(1+(vds-vip)/(lc*ucrit_a));

            // eq. 53
            lprime = ns*leff-dl_a+(vds+vip)/ucrit_a;
            // 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;

            `y_fv(fv,y)

            irprime = y*(1.0 + y);

            // eq. 57
            fv=(vp-vd)/Vt;

            `y_fv(fv,y)

            irev = y*(1.0 + y);

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

            // eq. 69
            nq = 1+0.5*gamma_a/sqrt(vp+phi_a+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) begin
                qb = (-gamma_a*sqrt(vp+phi_a+1e-6))/Vt-(nq-1)/nq*qi;
            end else begin
                qb = -vgprime/Vt;
            end
            // eq. 76 (qox removed since it is assumed to be zero)
            qg = -qi-qb;

            if (e0!=0) begin
                // eq. 61
                beta0prime = beta0*(1+cox_p/(e0*`EPSSI)*qb0);
                // eq. 62
                beta = beta0prime/(1+cox_p/(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*nslope*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
                `expLin(tmp,(-ibb_a*lc)/vib)
                idb = ids*iba/ibb_a*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
                        `expLin(tmp,-vs/(n*Vt))
                        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
                        `expLin(tmp,-vd/(n*Vt))
                        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
                        `expLin(tmp,-vd/(n*Vt))
                        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
                        `expLin(tmp,-vs/(n*Vt))
                        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;
            v_bp_dp = MOStype*V(b,dp);
            v_bp_sp = MOStype*V(b,sp);
            if (cbs0>0) begin
                if (v_bp_sp>fc*pb) begin
                    cbs = cbs+cbs0/pow(1-fc,mj)*(1+mj*(v_bp_sp-pb*fc))/(pb*(1-fc));
                end else begin
                    cbs = cbs+cbs0/pow(1-v_bp_sp,mj);
                end
            end
            if (cbd0>0) begin
                if (v_bp_dp>fc*pb) begin
                    cbd = cbd+cbd0/pow(1-fc,mj)*(1+mj*(v_bp_dp-pb*fc))/(pb*(1-fc));
                end else begin
                    cbd = cbd+cbd0/pow(1-v_bp_dp,mj);
                end
            end
            if (cbs0sw>0) begin
                if (v_bp_sp>fcsw*pbsw) begin
                    cbs = cbs+cbs0sw/pow(1-fcsw,mjsw)*(1+mjsw*(v_bp_sp-pbsw*fcsw))/(pbsw*(1-fcsw));
                end else begin
                    cbs = cbs+cbs0sw/pow(1-v_bp_sp,mjsw);
                end
            end
            if (cbd0sw>0) begin
                if (v_bp_dp>fcsw*pbsw) begin
                    cbd = cbd+cbd0sw/pow(1-fcsw,mjsw)*(1+mjsw*(v_bp_dp-pbsw*fcsw))/(pbsw*(1-fcsw));
                end else begin
                    cbd = cbd+cbd0sw/pow(1-v_bp_dp,mjsw);
                end
            end

        end       //Bias-dependent model evaluation

        begin     //Define branch sources

            I(dpsp) <+ MOStype*mode*ids;
            I(dpsp) <+ gmin*V(dpsp);

            if (mode>0) begin
                I(dpb) <+ MOStype*idb;
                I(dpb) <+ gmin*V(dpb);

                qdtx = qdt;
                qstx = qst;

            end else begin
                I(spb) <+ MOStype*idb;
                I(spb) <+ gmin*V(spb);

                qdtx = qst;
                qstx = qdt;

            end // else: !if(mode>0)

            I(dpb) <+ MOStype*ddt(qdtx);
            I(spb) <+ MOStype*ddt(qstx);
            I(gb)  <+ MOStype*ddt(qgt);
            //$strobe("V(dpb): %e qdtx=%e V(spb): %e qstx=%e V(gb): %e qgt: %e" ,V(dpb),qdtx,V(spb),qstx,V(gb),qgt);

            I(b,sp) <+ MOStype*ibsj;
            I(b,dp) <+ MOStype*ibdj;

            qjs = cbs * v_bp_sp;
            qjd = cbd * v_bp_dp;
            I(b,sp) <+ MOStype*ddt(qjs);
            I(b,dp) <+ MOStype*ddt(qjd);
            //$strobe("v_bp_sp: %e cbs=%e v_bp_dp: %e cbd=%e" ,v_bp_sp,cbs,v_bp_dp,cbd);

            I(gpsp) <+ ddt(cgso_s*V(gpsp));
            I(gpdp) <+ ddt(cgdo_s*V(gpdp));
            I(gb)   <+ ddt(cgbo_s*V(gb)) + gmin*V(dpsp);
            //$strobe("V(gpsp): %e cgso_s=%e V(gpdp): %e cgdo_s=%e V(gb): %e cgbo_s: %e" ,V(gpsp),cgso_s,V(gpdp),cgdo_s,V(gb),cgbo_s);

            I(ddp) <+ V(ddp)/RDeff;
            I(ssp) <+ V(ssp)/RSeff;

        end // begin

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

    end //analog

endmodule
