/* mosq3.f -- translated by f2c (version of 3 February 1990  3:36:42).
   You must link the resulting object file with the libraries:
	-lF77 -lI77 -lm -lc   (in that order)
*/

#include "f2c.h"

/*<       subroutine mosq3(vds,vbs,vgs,vdsat,vth,vbin,gamasd,cox,phi, >*/
/*<      1   qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb) >*/
/* Subroutine */ int mosq3_(vds, vbs, vgs, vdsat, vth, vbin, gamasd, cox, phi,
	 qg, qc, qb, cggb, cgdb, cgsb, cbgb, cbdb, cbsb)
doublereal *vds, *vbs, *vgs, *vdsat, *vth;
doublereal *vbin, *gamasd;
doublereal *cox, *phi, *qg, *qc, *qb, *cggb, *cgdb, *cgsb, *cbgb, *cbdb, *
	cbsb;
{
    /* System generated locals */
    static doublereal equiv_0[1], equiv_1[1], equiv_2[1], equiv_3[1];

    /* Builtin functions */
    double sqrt();

    /* Local variables */
    static doublereal arga, vbix;
#define utra (equiv_3)
    static doublereal onxl;
#define uexp (equiv_2)
    static doublereal vdsx, vgsx, dadfb, dadco, gamma, dadvd;
#define alpha (equiv_0)
    static doublereal fbody;
#define theta (equiv_1)
    static doublereal phibs, gamma2, dfbdvb, dcodvb, dqbdvb, dcodvd, xl, 
	    gammas;
#define xlamda (equiv_0)
    static doublereal dcodvg, dfsdvb, onfbdy, qbonco, dsqdvb;
#define xkappa (equiv_3)
    static doublereal fbodys, dvtdvb, dvtdvd, fnarrw, sqphbs, fshort, cdo;
#define eta (equiv_2)
    static doublereal vfb, vgb, vbi;
#define vbp (equiv_1)
    static doublereal von;

/*<       implicit double precision (a-h,o-z) >*/

/*<       equivalence (xlamda,alpha),(vbp,theta),(uexp,eta),(utra,xkappa) >*/

/*     charge equations are referenced to bulk */

/*<       vgb=vgs-vbs >*/
    vgb = *vgs - *vbs;
/*<       vfb=vbi-phi >*/
    vfb = vbi - *phi;
/*<       onxl=1.0d0/xl >*/
    onxl = 1. / xl;
/*<       phibs=sqphbs*sqphbs >*/
    phibs = sqphbs * sqphbs;

/*     body effect */

/*<       gammas=gamma*fshort >*/
    gammas = gamma * fshort;
/*<       fbodys=gammas/(sqphbs+sqphbs)*0.5d0 >*/
    fbodys = gammas / (sqphbs + sqphbs) * .5;
/*<       fbody=fbodys+fnarrw >*/
    fbody = fbodys + fnarrw;
/*<       onfbdy=1.0d0/(1.0d0+fbody) >*/
    onfbdy = 1. / (fbody + 1.);
/*<       dfbdvb=-fbodys*dsqdvb/sqphbs+fbodys*dfsdvb/fshort >*/
    dfbdvb = -fbodys * dsqdvb / sqphbs + fbodys * dfsdvb / fshort;
/*<       qbonco=gammas*sqphbs+fnarrw*phibs >*/
    qbonco = gammas * sqphbs + fnarrw * phibs;
/*<       dqbdvb=gammas*dsqdvb+gamma*dfsdvb*sqphbs-fnarrw >*/
    dqbdvb = gammas * dsqdvb + gamma * dfsdvb * sqphbs - fnarrw;

/* .....static feedback effect */

/*<       vbix=vbi-eta*vds >*/
    vbix = vbi - *eta * *vds;

/* .....threshold voltage */

/*<       vth=vbix+qbonco >*/
    *vth = vbix + qbonco;
/*<       dvtdvd=-eta >*/
    dvtdvd = -(*eta);
/*<       dvtdvb=dqbdvb >*/
    dvtdvb = dqbdvb;

/*     branch according to region of operation */

/*<       if (vgs.le.vth) go to 800 >*/
    if (*vgs <= *vth) {
	goto L800;
    }
/*<       vgsx=dmax1(vgs,von) >*/
    vgsx = max(*vgs,von);

/*     branch on vds=0.0d0 */

/*<       vdsx=dmin1(vds,vdsat) >*/
    vdsx = min(*vds,*vdsat);
/*<       if ( vdsx.eq.0.0d0 ) go to 900 >*/
    if (vdsx == 0.) {
	goto L900;
    }
/*<       cdo=vgsx-vth-0.5d0*(1.0d0+fbody)*vdsx >*/
    cdo = vgsx - *vth - (fbody + 1.) * .5 * vdsx;
/*<       dcodvg=1.0d0 >*/
    dcodvg = 1.;
/*<       if (vds.lt.vdsat) dcodvd=-dvtdvd-0.5d0*(1.0d0+fbody) >*/
    if (*vds < *vdsat) {
	dcodvd = -dvtdvd - (fbody + 1.) * .5;
    }
/*<       dcodvb=-dvtdvb-0.5d0*dfbdvb*vdsx >*/
    dcodvb = -dvtdvb - dfbdvb * .5 * vdsx;

/* .....charge terms */

/*< 420   continue >*/
/* L420: */
/*<       arga=(1.0d0+fbody)*vdsx*vdsx/(12.0d0*cdo) >*/
    arga = (fbody + 1.) * vdsx * vdsx / (cdo * 12.);
/*<       dadco=-arga/cdo >*/
    dadco = -arga / cdo;
/*<       if (vds.lt.vdsat) dadvd=arga/vdsx >*/
    if (*vds < *vdsat) {
	dadvd = arga / vdsx;
    }
/*<       dadfb=arga*onfbdy >*/
    dadfb = arga * onfbdy;

/* .....gate charge */

/*<       qg=cox*(vgs-vbix-0.5d0*vdsx+arga) >*/
    *qg = *cox * (*vgs - vbix - vdsx * .5 + arga);
/*<       cggb=cox*(1.0d0+dadco*dcodvg) >*/
    *cggb = *cox * (dadco * dcodvg + 1.);
/*<       if (vds.lt.vdsat) cgdb=cox*(-dvtdvd-0.5d0+dadvd+dadco*dcodvd) >*/
    if (*vds < *vdsat) {
	*cgdb = *cox * (-dvtdvd - .5 + dadvd + dadco * dcodvd);
    }
/*<       cgsb=-cggb-cgdb-cox*(dadco*dcodvb+dadfb*dfbdvb) >*/
    *cgsb = -(*cggb) - *cgdb - *cox * (dadco * dcodvb + dadfb * dfbdvb);

/* .....bulk charge */

/*<       arga=arga*fbody >*/
    arga *= fbody;
/*<       dadco=dadco*fbody >*/
    dadco *= fbody;
/*<       if (vds.lt.vdsat) dadvd=dadvd*fbody >*/
    if (*vds < *vdsat) {
	dadvd *= fbody;
    }
/*<       dadfb=dadfb*(1.0d0+fbody+fbody) >*/
    dadfb *= fbody + 1. + fbody;

/*<       qb=-cox*(qbonco+0.5d0*fbody*vdsx-arga) >*/
    *qb = -(*cox) * (qbonco + fbody * .5 * vdsx - arga);
/*<       cbgb=cox*dadco*dcodvg >*/
    *cbgb = *cox * dadco * dcodvg;
/*<       if (vds.lt.vdsat) cbdb=-cox*(0.5d0*fbody-dadvd-dadco*dcodvd) >*/
    if (*vds < *vdsat) {
	*cbdb = -(*cox) * (fbody * .5 - dadvd - dadco * dcodvd);
    }
/*<       cbsb=-cbgb-cbdb >*/
/*<      1          +cox*(dqbdvb+(0.5d0*vdsx-dadfb)*dfbdvb-dadco*dcodvb) >*/
    *cbsb = -(*cbgb) - *cbdb + *cox * (dqbdvb + (vdsx * .5 - dadfb) * dfbdvb 
	    - dadco * dcodvb);
/*<       go to 1000 >*/
    goto L1000;

/* .....charge terms of vgs<vth */

/*< 800   continue >*/
L800:
/*<       if ( vgb.gt.vfb ) go to 810 >*/
    if (vgb > vfb) {
	goto L810;
    }
/*<       qg=cox*(vgb-vfb) >*/
    *qg = *cox * (vgb - vfb);
/*<       cggb=cox >*/
    *cggb = *cox;
/*<       go to  820 >*/
    goto L820;
/*< 810   continue >*/
L810:
/*<       gamma2=gammas*0.5d0 >*/
    gamma2 = gammas * .5;
/*<       arga=dsqrt(gamma2*gamma2+(vgb-vfb)) >*/
    arga = sqrt(gamma2 * gamma2 + (vgb - vfb));
/*<       qg=gammas*cox*(arga-gamma2) >*/
    *qg = gammas * *cox * (arga - gamma2);
/*<       cggb=0.5d0*cox*gammas/arga >*/
    *cggb = *cox * .5 * gammas / arga;
/*< 820   continue >*/
L820:
/*<       qb=-qg >*/
    *qb = -(*qg);
/*<       cbgb=-cggb >*/
    *cbgb = -(*cggb);
/*<       cgdb=0.0d0 >*/
    *cgdb = 0.;
/*<       cgsb=0.0d0 >*/
    *cgsb = 0.;
/*<       cbdb=0.0d0 >*/
    *cbdb = 0.;
/*<       cbsb=0.0d0 >*/
    *cbsb = 0.;
/*<       go to 1000 >*/
    goto L1000;

/*     special case vds=0.0d0 */

/*<   900 qg=cox*(vgs-vbi) >*/
L900:
    *qg = *cox * (*vgs - vbi);
/*<       qb=-cox*qbonco >*/
    *qb = -(*cox) * qbonco;
/*<       cggb=cox >*/
    *cggb = *cox;
/*<       cgdb=-cox*(0.5d0+dvtdvd) >*/
    *cgdb = -(*cox) * (dvtdvd + .5);
/*<       cgsb=-cox*(0.5d0-dvtdvb) >*/
    *cgsb = -(*cox) * (.5 - dvtdvb);
/*<       cbgb=0.0d0 >*/
    *cbgb = 0.;
/*<       cbdb=-0.5d0*cox*fbody >*/
    *cbdb = *cox * -.5 * fbody;
/*<       cbsb=cox*(dqbdvb+0.5d0*fbody) >*/
    *cbsb = *cox * (dqbdvb + fbody * .5);

/*     done */

/*<  1000 qc=-(qg+qb) >*/
L1000:
    *qc = -(*qg + *qb);
/*<       return >*/
    return 0;
/*<       end >*/
} /* mosq3_ */

#undef vbp
#undef eta
#undef xkappa
#undef xlamda
#undef theta
#undef alpha
#undef uexp
#undef utra


