/* dev_mos3  07/30/91  (XZ)  12/19/93, 94.01.13, 94.02.11
 * mos model equations: spice level 3 equivalent
 */
#include "ecah.h"
#include "branch.h"
#include "d_diode.h"
#include "d_mos.h"
#include "error.h"
#include "io.h"
#include "declare.h"
/*--------------------------------------------------------------------------*/
extern const struct ioctrl io;
/*--------------------------------------------------------------------------*/
void eval_mos3(branch_t *brh)
{
 static struct mos *x;
 static struct mmod *m;
 double phi_0, phi_vbs, phi_vbd, vbd;
 double beta;
 double sigma;
 double ids;
 double us, ueff;
 double fb, fs, fn;
 double wc, wp;
 double va, vb;
 double ep, ep_xd;
 double theta_vgs;
 double cox_w_le;
 double fb_1;
 double xd_2;
 double gamma_fs;
 double us2_vmax_le;
 double us_vds_vmax, us_vds_vmax2;
 double us_vdsat_vmax, us_vdsat_vmax2;
 double xjwp_sq, xjwp_sqrt;
 double vgs_vth_vdsat, ueff_sat, dueff_dvdsat, gdsat;
 double del_le, lfac;
 double vgs_vth_vds;
 double dus_dvgs;
 double dfs_dvbs;
 double dueff_dvds, dueff_dvgs;
 double dwp_dvbs, dwc_dvbs;
 double dvth_dvbs;
 double dbeta_dvds, dbeta_dvgs;
 double dids_dvds, dids_dvgs, dids_dvbs;
 double a = 0.0631354;
 double b = 0.8013292;
 double c = 0.01110777;
 
 double vdsx, vgsx;

 int use_vmax;
 int use_xj;

 x = (struct mos*)brh->x;
 m = x->m;

 use_vmax = (m->vmax != NOT_INPUT  &&  m->vmax != 0.);
 use_xj = (m->xj != NOT_INPUT  &&  m->xj > 0.);

 phi_0 = sqrt(m->phi);
 if (x->vbs <= 0.){
    phi_vbs = sqrt(m->phi - x->vbs);
    x->sbfwd = NO;
 }else{
    phi_vbs = sqrt(m->phi) / (1. + .5 * x->vbs / m->phi);
    x->sbfwd = YES;
    if (!io.suppresserrors){
       error((x->vbs < -.01) ? bPICKY : bTRACE,
          "%s: source fwd biased.  vbs=%g\n", printlabel(brh,NO), x->vbs);
    }
 }

 vbd = x->vbs - x->vds;
 if (vbd <= 0.){
    phi_vbd = sqrt(m->phi - vbd);
    x->dbfwd = NO;
 }else{
    phi_vbd = sqrt(m->phi) / (1. + .5 * vbd / m->phi);
    x->dbfwd = YES;
    if (!io.suppresserrors){
       error((vbd < -.01) ? bPICKY : bTRACE,
          "%s: drain fwd biased.  vbd=%g\n", printlabel(brh,NO), vbd);
    }
 }

						       /* threshold voltage */
 sigma = 8.15e-22 * m->eta / (m->cox * x->le * x->le * x->le);
 fn = (kPI * E_SI * m->delta) / (4. * m->cox * x->we);
 wp = m->xd * sqrt(m->d->pb - x->vbs);
 wc = use_xj ? a*m->xj + b*wp - c*wp*wp/m->xj : 0.;
 xjwp_sq = wp / (m->xj + wp);
 xjwp_sq *= xjwp_sq;
 xjwp_sqrt = sqrt(1. - xjwp_sq);
 fs = use_xj ? 1. - m->xj/x->le * ((m->ld + wc)/m->xj * xjwp_sqrt - m->ld/m->xj)
             : 1.;
 fb = m->gamma * fs / (2. * phi_vbs) + fn;
 fb_1 = .5 * (1. + fb);
 gamma_fs = m->gamma*fs*phi_vbs + fn*(m->phi - x->vbs);
 x->von = m->vfb + m->phi - sigma*x->vds + gamma_fs;
 x->vgst = x->vgs - x->von;

						     /* mobility modulation */
 theta_vgs = 1. + m->theta * x->vgst;
 us = m->uo / theta_vgs;
 us_vds_vmax = use_vmax ? 1. + us * x->vds / (m->vmax * x->le)
			: 1.;
 us_vds_vmax2 = us_vds_vmax * us_vds_vmax;
 us2_vmax_le = use_vmax ? -us*us / (m->vmax*x->le)
			: 1.;
 ueff = us / us_vds_vmax;
 cox_w_le = m->cox * x->we / x->le;
 beta = ueff * cox_w_le;

						      /* saturation voltage */
 va = x->vgst / (1. + fb);
 vb = use_vmax ? m->vmax * x->le / us : 0.;
 x->vdsat = use_vmax ? va + vb - sqrt(va*va + vb*vb) : va;
 vdsx = MIN(x->vds, x->vdsat);
 vgsx = MAX(x->vgs, x->von);
 ids = beta * (vgsx - x->von - fb_1*vdsx) * vdsx;

 x->saturated = (x->vds > x->vdsat);
 x->cutoff = (x->vgst < 0.);

 if (x->cutoff){
    lfac = 0.;
 }else if (!x->saturated){
    lfac = 1.;
 }else{					       /* channel-length modulation */
    us_vdsat_vmax = use_vmax ? 1. + us * x->vdsat/(m->vmax * x->le)
			     : 1.;
    us_vdsat_vmax2 = us_vdsat_vmax * us_vdsat_vmax;
    ueff_sat = us / us_vdsat_vmax;
    dueff_dvdsat = use_vmax ? us2_vmax_le / us_vdsat_vmax2
			    : 0.;
    vgs_vth_vdsat = vgsx - x->von - fb_1*x->vdsat;
    gdsat = (cox_w_le*dueff_dvdsat*x->vdsat + beta) * vgs_vth_vdsat  -
            beta * x->vdsat * (fb_1 - sigma);
    xd_2 = m->xd * m->xd;
    ep = m->kappa * ids / (gdsat * x->le);
    ep_xd = .5 * ep * xd_2;
    del_le = sqrt(m->kappa*xd_2*(x->vds - x->vdsat) + ep_xd*ep_xd) - ep_xd;
    lfac = 1. / (1. - (del_le/x->le));
 }

						     /* (trans)conductances */
 vgs_vth_vds = vgsx - x->von - fb_1*x->vds;
 dwp_dvbs = -.5 * m->xd / phi_vbs;
 dwc_dvbs = use_xj ? (b - 2.*c*wp/m->xj) * dwp_dvbs : 0.;
 dfs_dvbs = use_xj ? (m->ld + wc) * xjwp_sq / (x->le * m->xj * xjwp_sqrt) *
		     dwp_dvbs - xjwp_sqrt / x->le * dwc_dvbs
		   : 0.;
 dus_dvgs = -m->uo * m->theta / (theta_vgs * theta_vgs);
 dueff_dvgs = dus_dvgs / us_vds_vmax2;
 dueff_dvds = use_vmax ? us2_vmax_le / us_vds_vmax2 : 0.;
 dbeta_dvgs = dueff_dvgs * cox_w_le;
 dbeta_dvds = cox_w_le * dueff_dvds;
 dvth_dvbs = m->gamma * (phi_vbs*dfs_dvbs - .5*fs/phi_vbs);

 dids_dvds = (dbeta_dvds + beta)*vgs_vth_vds - beta*x->vds*(fb_1 - sigma);
 dids_dvgs = (dbeta_dvgs*vgs_vth_vds + beta) * x->vds;
 dids_dvbs = -beta * x->vds * (dvth_dvbs + x->vds*m->gamma/(4.*phi_vbs)*
             dfs_dvbs);

 x->ids = ids * lfac;
 x->gds = dids_dvds * lfac;
 x->gm  = dids_dvgs * lfac;
 x->gmb = dids_dvbs * lfac;
}  
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
