/* d_mos2.c  94.04.28
 * Copyright 1983-1992   Albert Davis
 * mos model equations: spice level 2 equivalent
 */
#include "ecah.h"
#include "branch.h"
#include "d_mos.h"
#include "error.h"
#include "io.h"
#include "options.h"
#include "declare.h"
/*--------------------------------------------------------------------------*/
	void	eval_mos2(branch_t*);
/*--------------------------------------------------------------------------*/
#define short_channel	(m->xj != NOT_INPUT  &&  m->xj > 0.)
#define do_subthreshold	(m->nfs != 0.)

extern const struct ioctrl io;
extern const struct options opt;
extern const double temp;
/*--------------------------------------------------------------------------*/
void eval_mos2(branch_t *brh)
{
static struct mos *x;
static struct mmod *m;
 double sarg, dsarg_dvbs/*ds*/;
 double bargx;
 double gamma_s/*mb*/, dgamma_s_dvds/*ds*/, dgamma_s_dvbs/*mb*/;
 double beta;
 double body;
 double dvdsat_dvgs = 0./*m*/;
 double dvdsat_dvbs = 0./*mb*/;
 double vgsx;
 double vdsx;
 double ufact;
 double dudvgs/*m*/, dudvds/*ds*/, dudvbs/*mb*/;
 double clfact;
 double dldvgs/*m*/, dldvds/*ds*/, dldvbs/*mb*/;
 double dodvbs/*mb*/;
 double xn = 0.;
 double vtxn = 0.;
 double expg;
 double dxndvb =0./*mb*/;	/* subthreshold only */
 double ids_on, didvds;
 double vt/*ds*/;
 double v_phi_d, v_phi_s, v_phi_ds;
 double d2sdb2;
 double barg, dbarg_dvbs, d2bdb2;
 double bsarg, dbsarg_dvbs;
 double vbdsat;
 double dgddb2;
 double vc, vbin, vc_eta;
 double bodys;
 double lambda;
 double vbd, vbs;
 double ueff = 0.;

 int use_vmax;

 x = (struct mos*)brh->x;
 m = x->m;
  
 use_vmax = m->vmax != NOT_INPUT;
 vt = (K/Q) * temp;
 
 vbs = x->vbs;
 v_phi_s = m->phi - vbs;
 if (vbs <= 0.){
    sarg = sqrt(v_phi_s);
    dsarg_dvbs = -.5 / sarg;
    d2sdb2 = .5 * dsarg_dvbs / v_phi_s;
    x->sbfwd = NO;
 }else{
    double sphi3 = pow(m->phi, 1.5);
    if (opt.mosflags & 0001  &&  vbs > m->phi){
       error(bDANGER,"%s: vbs(%g) > phi(%g)\n", printlabel(brh,0),vbs,m->phi);
       vbs = m->phi;
       v_phi_s = 0.;
    }
    sarg = sqrt(m->phi) / (1. + .5 * vbs / m->phi);
    dsarg_dvbs = -.5 * sarg * sarg / sphi3;
    d2sdb2 = -dsarg_dvbs * sarg / sphi3;
    x->sbfwd = YES;
    if (!io.suppresserrors){
       error((x->vbs > m->phi) ? bPICKY : bTRACE,
          "%s: source fwd biased. vbs=%g\n", printlabel(brh,NO), x->vbs);
    }
 }
 
 if (opt.mosflags & 0004){
    vbd = vbs - x->vds;
 }else{
    vbd = x->vbs - x->vds;
 }
 v_phi_d = m->phi - vbd;
 if (vbd <= 0.){
    barg = sqrt(v_phi_d);
    dbarg_dvbs = -.5 / barg;
    d2bdb2 = .5 * dbarg_dvbs / v_phi_d;
    x->dbfwd = NO;
 }else{
    double sphi3 = pow(m->phi, 1.5);
    if (opt.mosflags & 0002  &&  vbd > m->phi){
       error(bDANGER,"%s: vbd(%g) > phi(%g)\n", printlabel(brh,0),vbd,m->phi);
       vbd = m->phi;
       v_phi_d = 0.;
    }
    barg = sqrt(m->phi) / (1. + .5 * vbd / m->phi);
    dbarg_dvbs = -.5 * barg * barg / sphi3;
    d2bdb2 = -dbarg_dvbs * barg / sphi3;
    x->dbfwd = YES;
    if (!io.suppresserrors){
       error((vbd > m->phi) ? bPICKY : bTRACE,
          "%s: drain fwd biased. vbd=%g\n",printlabel(brh,NO), x->vbs-x->vds);
    }
 }
 
 if (short_channel){
    double wd, ws;
    double alpha_d, alpha_s;
    double dalpha_d_dvds, dalpha_d_dvbs;
    double dalpha_s_dvbs;
    double argxs, argxd;
    double args, argd;
    double dasdb2, daddb2;

    wd = m->xd * barg;
    argxd = 1. + 2.*wd/m->xj;
    argd = sqrt(argxd);
    alpha_d = x->relxj * (argd - 1.);
    dalpha_d_dvds = m->xd / (4. * x->le * argd * barg);
    dalpha_d_dvbs = -dalpha_d_dvds;
    
    ws = m->xd * sarg;
    argxs = 1. + 2.*ws/m->xj;
    args = sqrt(argxs);
    alpha_s = x->relxj * (args - 1.);
    dalpha_s_dvbs = -m->xd / (4. * x->le * args * sarg);
    
    gamma_s = m->gamma * (1. - alpha_s - alpha_d);
    dgamma_s_dvds = -m->gamma *  dalpha_d_dvds;
    dgamma_s_dvbs = -m->gamma * (dalpha_d_dvbs + dalpha_s_dvbs);

    dasdb2 = -m->xd*(d2sdb2 + dsarg_dvbs*dsarg_dvbs*m->xd / (m->xj*argxs))
    				/ (x->le*args);
    daddb2 = -m->xd*(d2bdb2 + dbarg_dvbs*dbarg_dvbs*m->xd / (m->xj*argxd))
    				/ (x->le*argd);
    dgddb2 = -.5 * m->gamma * (dasdb2 + daddb2);

    if (gamma_s <= 0. && m->gamma > 0. && !io.suppresserrors){
       error(bWARNING, "%s: gamma is negative\n", printlabel(brh,NO));
       error(bTRACE  , "+   gamma_s=%g, alpha_s=%g, alpha_d=%g\n",
       			    gamma_s,    alpha_s,    alpha_d);
    }
 }else{
    gamma_s = m->gamma;
    dgamma_s_dvds = dgamma_s_dvbs = 0.;
    dgddb2 = 0.;
 }
 
 
 vbin = m->vbi + x->eta_1 * v_phi_s;
 x->von = vbin + gamma_s * sarg;
 
 dodvbs = -x->eta_1 + dgamma_s_dvbs * sarg + gamma_s * dsarg_dvbs;
 if (do_subthreshold){
    double cdonco;
    dxndvb = 2. * dgamma_s_dvbs * dsarg_dvbs
       			+ gamma_s * d2sdb2 + dgddb2 * sarg;
    dodvbs += vt * dxndvb;
    cdonco = - (gamma_s * dsarg_dvbs + dgamma_s_dvbs * sarg) + x->eta_1;
    xn = 1. + m->cfsox + cdonco;
    vtxn = vt * xn;
    x->von += vtxn;
    x->subthreshold = (x->vgs < x->von);
    x->cutoff = NO;
 }else if (x->vgs < x->von){
    x->cutoff = YES;
    x->subthreshold = NO;
    x->ids = 0.;
    x->gm = 0.;
    if (opt.mosflags == 0100){
       x->gds = opt.gmin;
    }else{
       x->gds = 0.;
    }
    x->gmb = 0.;
    x->vgst = x->vgs - x->von;
    return;
 }else{
    x->cutoff = x->subthreshold = NO;
 }

 vgsx = (x->subthreshold) ? x->von : x->vgs;
 vc = vgsx - vbin;
 vc_eta = vc / x->eta;
 x->vgst = x->vgs - x->von;

 if (m->uexp != NOT_INPUT  &&  x->vgst > m->vbp){
    ufact = pow(m->vbp/x->vgst, m->uexp);
    dudvgs = -ufact * m->uexp / x->vgst;
    dudvds = 0.;	/* wrong, but as per spice2 */
    dudvbs = dodvbs * ufact * m->uexp / x->vgst;
 }else{
    ufact = 1.;
    dudvgs = dudvds = dudvbs = 0.;
 }

 if (use_vmax){
    double sarg3, gammad;
    double v1, v2, xv, a1, b1, c1, d1;
    double a, b, c, r, s, r3, s2, p, p0, p2, y3;
    double xvalid = 0.;
    double x4[8];
    int iknt, j;
    int root_count;
    sarg3 = sarg*sarg*sarg;
    gammad = gamma_s / x->eta;
    ueff = m->uo * ufact;
    v1 = vc_eta + v_phi_s;
    v2 = v_phi_s;
    xv = m->vmax * x->le / ueff;
    a1 = gammad * (4./3.);
    b1 = -2. * (v1+xv);
    c1 = -2. * gammad * xv;				/* end of scope */
    d1 = 2.*v1*(v2+xv) - v2*v2 - (4./3.)*gammad*sarg3; 	/* xv, v1, v2, sarg3 */
    a = -b1;
    b = a1 * c1 - 4. * d1;
    c = -d1 * (a1*a1 - 4.*b1) - c1*c1;
    r = -a*a / 3. + b;
    s = 2. * a*a*a / 27. - a*b / 3. + c;		/* b, c */
    r3 = r*r*r;						/* r */
    s2 = s*s;
    p = s2 / 4. + r3 / 27.;				/* r3 */
    p0 = fabs(p);
    p2 = sqrt(p0);
    if (p < 0.){					/* p */
       double ro, fi;
       ro = pow((s2 / 4. + p0), (1./6.));		/* s2, p0 */
       fi = atan(-2. * p2 / s);
       y3 = 2. * ro * cos(fi/3.) - a / 3.;
    }else{
       double p3, p4;
       p3 = pow((fabs(-s/2.+p2)), (1./3.));
       p4 = pow((fabs(-s/2.-p2)), (1./3.));		/* s, p2 */
       y3 = p3 + p4 - a / 3.;				/* a */
    }

    iknt = 0;
    if (a1*a1 / 4. - b1 + y3  < 0.  &&  y3*y3 / 4. - d1  < 0.){
       error(bWARNING,
          "%s: internal error: a3,b4, a1=%g, b1=%g, y3=%g, d1=%g\n",
           printlabel(brh,NO),	      a1,    b1,    y3,    d1);
    }else{
       double a3, b3;
       int i;
       a3 = sqrt(a1*a1 / 4. - b1 + y3);
       b3 = sqrt(y3*y3 / 4. - d1);
       for (i = 0;   i < 4;   i++){
          double delta4;
          double a4, b4;
          static const double sig1[4] = {1., -1., 1., -1.};
          static const double sig2[4] = {1., 1., -1., -1.};
          a4 = a1 / 2. + sig1[i] * a3;
          b4 = y3 / 2. + sig2[i] * b3;			/* y3 */
          delta4 = a4*a4 / 4. - b4;
          if (delta4 >= 0.){
             x4[iknt++] = - a4 / 2. + sqrt(delta4);
             x4[iknt++] = - a4 / 2. - sqrt(delta4);	/* i */
          }
       }
    }

    root_count = 0;
    for (j = 0;   j < iknt;   j++){			/* iknt */
       if (x4[j] > 0.){
          double poly4;
          poly4 =      x4[j]*x4[j]*x4[j]*x4[j]	/* ~= 0, used as check	*/
		+ a1 * x4[j]*x4[j]*x4[j]	/* roundoff error not	*/
		+ b1 * x4[j]*x4[j]		/* propagated, so ok	*/
		+ c1 * x4[j]
		+ d1;					/* a1, b1, c1, d1 */
          if (fabs(poly4) <= 1e-6){
             root_count++;
             if (root_count <= 1)	/* xvalid = min(x4[j]) */
                xvalid=x4[j];
             if (x4[j] <= xvalid)
                xvalid=x4[j];				/* x4[], j */
	  }else{
	  }
       }
    }
    if (root_count <= 0){				/* root_count */
       if (!io.suppresserrors)
          error(bWARNING, "%s: Baum's theory rejected\n", printlabel(brh,NO));
       use_vmax = NO;
    }else{
       x->vdsat = xvalid*xvalid - v_phi_s;
    }
 }

 if (!use_vmax){
    if (gamma_s > 0.){
       double argv;
       argv = vc_eta + v_phi_s;
       if (argv > 0.){
          double gammad, gammd2, arg;
          gammad = gamma_s / x->eta;
          gammd2 = gammad * gammad;
          arg = sqrt(1. + 4. * argv / gammd2);
          x->vdsat = vc_eta  +  gammd2 * (1.-arg) / 2.;
          dvdsat_dvgs = (1. - 1./arg) / x->eta;
          dvdsat_dvbs = (gammad * (1.-arg) + 2.*argv / (gammad*arg))
	  		/ x->eta * dgamma_s_dvbs
			+ 1./arg + x->eta_1 * dvdsat_dvgs;
       }else{
          x->vdsat = 0.;
          dvdsat_dvgs = dvdsat_dvbs = 0.;
	  if (!io.suppresserrors){
             error(bWARNING, "%s: argv is negative\n", printlabel(brh,NO));
	     error(bTRACE  , "+   vc=%g, argv=%g, vdsat=%g\n",
				  vc,    argv,    x->vdsat);
	     error(bTRACE  , "+   vds=%g, vgs=%g, vbs=%g\n",
				  x->vds, x->vgs, x->vbs);
	  }
       }
    }else{
       x->vdsat = vc_eta;
       dvdsat_dvgs = 1.;
       dvdsat_dvbs = 0.;
    }
 }

 if (x->vdsat < 0.){
    error(bWARNING, "%s: calculated vdsat (%g) < 0.  using vdsat = 0.\n",
    		printlabel(brh,NO), x->vdsat);
    x->vdsat = 0.;
 }
 
 vbdsat = vbs - x->vdsat;
 v_phi_ds = m->phi - vbdsat;
 if (vbdsat <= 0.){
    bsarg = sqrt(v_phi_ds);
    dbsarg_dvbs = -.5 / bsarg;
 }else{
    double sphi3 = pow(m->phi, 1.5);
    bsarg = sqrt(m->phi) / (1. + .5 * vbdsat / m->phi);
    dbsarg_dvbs = -.5 * bsarg * bsarg / sphi3;
 }
 bodys = bsarg*bsarg*bsarg - sarg*sarg*sarg;

 if (use_vmax){
    double argv, vqchan, dqdsat, vl, dfunds, dfundg, dfundb, gdbdvs;
    gdbdvs = 2. * gamma_s * (bsarg*bsarg*dbsarg_dvbs - sarg*sarg*dsarg_dvbs);
    argv = vc_eta - x->vdsat;
    vqchan = argv - gamma_s * bsarg;
    dqdsat = -1. + gamma_s * dbsarg_dvbs;
    vl = m->vmax * x->le;
    dfunds = vl * dqdsat - ueff * vqchan;
    dfundg = (vl - ueff * x->vdsat) / x->eta;
    dfundb = -vl * (1. + dqdsat - x->eta_1 / x->eta)
       		+ ueff * (gdbdvs - dgamma_s_dvbs * bodys / 1.5) / x->eta;
    dvdsat_dvgs = -dfundg / dfunds;
    dvdsat_dvbs = -dfundb / dfunds;
 }
 
 if (m->lambda == NOT_INPUT){
    if (x->vds != 0.){
       double dldsat;
       if (use_vmax){
	  double xdv, xlv, argv, xls;
          xdv = m->xd / sqrt(m->neff);
          xlv = m->vmax * xdv / (2. * ueff);
          argv = x->vds - x->vdsat;
          if (argv < 0.)
             argv = 0.;
          xls = sqrt(xlv*xlv + argv);
          lambda = (xls-xlv) * xdv / (x->le * x->vds);
          dldsat = xdv / (2. * xls * x->le);
       }else{
          double argv, sargv, dl;
          argv = (x->vds - x->vdsat) / 4.;
	  sargv = sqrt(1. + argv*argv);
	  if (argv + sargv >= 0.){
	     dl = m->xd * sqrt(argv + sargv);
             lambda = dl / (x->le * x->vds);
             dldsat = lambda * x->vds / (8. * sargv);
       	  }else{
	     lambda = 0.;
	     dldsat = 0.;
	     error(bWARNING, "%s: internal error: vds(%g) < vdsat(%g)\n",
	     		      printlabel(brh,NO), x->vds,   x->vdsat);
	  }
       }
       dldvgs =  dvdsat_dvgs * dldsat;
       dldvds =              - dldsat;
       dldvbs =  dvdsat_dvbs * dldsat;
    }else{
       lambda = 0.;
       dldvgs = dldvds = dldvbs = 0.;
    }
 }else{
    lambda = m->lambda;
    dldvgs = dldvbs = 0.;
    dldvds = -lambda;
 }

 clfact = (1. - lambda * x->vds);
 if (clfact < m->xwb/x->le){
    double leff, dfact;
    leff = m->xwb / (2. - (clfact * x->le / m->xwb));
    clfact = leff / x->le;
    dfact = (leff * leff) / (m->xwb * m->xwb);
    dldvgs *= dfact;
    dldvds *= dfact;
    dldvbs *= dfact;
 }

 x->saturated = (x->vds > x->vdsat);
 vdsx =  (x->saturated) ? x->vdsat : x->vds;
 bargx = (x->saturated) ? bsarg : barg;
 body = bargx*bargx*bargx - sarg*sarg*sarg;
 expg = (x->subthreshold) ? exp(x->vgst / vtxn) : 1.;
 beta = x->beta * ufact / clfact;

 ids_on = beta * ((vc - x->eta_2 * vdsx) * vdsx  - (2./3.) * gamma_s * body);
 didvds = beta * (vc  -  x->eta * vdsx  -  gamma_s * bargx);

 x->ids = ids_on * expg;

 x->gm = beta * vdsx;
 x->gm += ids_on * (dudvgs/ufact - dldvgs/clfact);
 if (x->saturated)
    x->gm += didvds * dvdsat_dvgs;
 if (x->subthreshold){
    x->gm = ids_on / vtxn;
    if (x->saturated)
       x->gm += didvds * dvdsat_dvgs;
    x->gm *= expg;
 }
 x->gds = (x->saturated) ? 0.: didvds;
 x->gds += ids_on * (dudvds/ufact - dldvds/clfact);
 if (short_channel)
    x->gds -= beta * (2./3.) * body * dgamma_s_dvds;
 if (x->subthreshold){
    double dodvds, dxndvd, gmw;
    dxndvd = dgamma_s_dvds * dsarg_dvbs;
    dodvds = dgamma_s_dvds * sarg + vt * dxndvd;
    gmw = x->ids * x->vgst / (vtxn * xn);
    x->gds *= expg;
    x->gds -= x->gm * dodvds + gmw * dxndvd;
 }
 x->gmb = beta * (x->eta_1 * vdsx - gamma_s * (sarg - bargx));
 x->gmb += ids_on * (dudvbs/ufact - dldvbs/clfact);
 if (short_channel)
    x->gmb -= beta * (2./3.) * body * dgamma_s_dvbs;
 if (x->saturated)
    x->gmb += didvds * dvdsat_dvbs;
 if (x->subthreshold){
    double gmw;
    gmw = x->ids * x->vgst / (vtxn * xn);
    x->gmb += beta * dodvbs * vdsx;
    x->gmb *= expg;
    x->gmb -= x->gm * dodvbs + gmw * dxndvb;
 }
}
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
