text: /**********
  Copyright 1991 Regents of the University of California.  All rights reserved.
  Author:	1987 Kartikeya Mayaram, U. C. Berkeley CAD Group
  **********/
  
  /* Functions to compute terminal conductances & currents. */
  
  #include "ngspice/ngspice.h"
  #include "ngspice/numglobs.h"
  #include "ngspice/numenum.h"
  #include "ngspice/twomesh.h"
  #include "ngspice/twodev.h"
  #include "ngspice/bool.h"
  #include "ngspice/spmatrix.h"
  #include "twoddefs.h"
  #include "twodext.h"
  
  
  void 
     NUMD2conductance(TWOdevice *pDevice, BOOLEAN tranAnalysis, 
                      double *intCoeff, double *gd)
  {
    TWOcontact *pContact = pDevice->pFirstContact;
    double *incVpn;
    BOOLEAN deltaVContact = FALSE;
    
    /* 
     * store the new rhs for computing the incremental quantities
     * with the second to last node. solve the system of equations
     */
    incVpn = pDevice->dcDeltaSolution;
    storeNewRhs( pDevice, pDevice->pLastContact );
    spSolve( pDevice->matrix, pDevice->rhs, incVpn, NULL, NULL);
    
    incVpn = pDevice->dcDeltaSolution;
    *gd = contactConductance( pDevice, pContact, deltaVContact, incVpn,
  			  tranAnalysis, intCoeff );
    *gd *= - GNorm * pDevice->width * LNorm;
  }
  
  void 
     NBJT2conductance(TWOdevice *pDevice, BOOLEAN tranAnalysis, 
                      double *intCoeff, double *dIeDVce, double *dIcDVce, 
  	   	    double *dIeDVbe, double *dIcDVbe)
  {
    TWOcontact *pEmitContact = pDevice->pLastContact;
    TWOcontact *pColContact = pDevice->pFirstContact;
    TWOcontact *pBaseContact = pDevice->pFirstContact->next;
    double width = pDevice->width;
    double *incVce, *incVbe;
    
    /* 
     * store the new rhs for computing the incremental quantities
     * incVce (dcDeltaSolution) and incVbe (copiedSolution) are used to
     * store the incremental quantities associated with Vce and Vbe
     */
    
    incVce = pDevice->dcDeltaSolution;
    incVbe = pDevice->copiedSolution;
    storeNewRhs( pDevice, pColContact );
    spSolve( pDevice->matrix, pDevice->rhs, incVce, NULL, NULL);
    storeNewRhs( pDevice, pBaseContact );
    spSolve( pDevice->matrix, pDevice->rhs, incVbe, NULL, NULL);
    
    *dIeDVce = contactConductance( pDevice, pEmitContact, FALSE, incVce,
  				tranAnalysis, intCoeff );
    *dIeDVbe = contactConductance( pDevice, pEmitContact, FALSE, incVbe,
  				tranAnalysis, intCoeff );
    
    *dIcDVce = contactConductance( pDevice, pColContact, TRUE, incVce,
  				tranAnalysis, intCoeff );
    *dIcDVbe = contactConductance( pDevice, pColContact, FALSE, incVbe,
  				tranAnalysis, intCoeff );
    *dIeDVce *= GNorm * width * LNorm;
    *dIcDVce *= GNorm * width * LNorm;
    *dIeDVbe *= GNorm * width * LNorm;
    *dIcDVbe *= GNorm * width * LNorm;
  }
  
  void 
     NUMOSconductance(TWOdevice *pDevice, BOOLEAN tranAnalysis, double *intCoeff, 
                      struct mosConductances *dIdV)
  {
    TWOcontact *pDContact = pDevice->pFirstContact;
    TWOcontact *pGContact = pDevice->pFirstContact->next;
    TWOcontact *pSContact = pDevice->pFirstContact->next->next;
    double width = pDevice->width;
    double *incVdb, *incVsb, *incVgb;
    
    /* 
     * store the new rhs for computing the incremental quantities
     * incVdb (dcDeltaSolution)
     */
    
    incVdb = pDevice->dcDeltaSolution;
    incVsb = pDevice->copiedSolution;
    incVgb = pDevice->rhsImag;
    storeNewRhs( pDevice, pDContact );
    spSolve( pDevice->matrix, pDevice->rhs, incVdb, NULL, NULL);
    storeNewRhs( pDevice, pSContact );
    spSolve( pDevice->matrix, pDevice->rhs, incVsb, NULL, NULL);
    storeNewRhs( pDevice, pGContact );
    spSolve( pDevice->matrix, pDevice->rhs, incVgb, NULL, NULL);
    
    dIdV->dIdDVdb = contactConductance( pDevice, pDContact, TRUE,
  				     incVdb, tranAnalysis, intCoeff );
    dIdV->dIsDVdb = contactConductance( pDevice, pSContact, FALSE,
  				     incVdb, tranAnalysis, intCoeff );
    dIdV->dIgDVdb = GateTypeConductance( pDevice, pGContact, FALSE,
  				      incVdb, tranAnalysis, intCoeff );
    dIdV->dIdDVsb = contactConductance( pDevice, pDContact, FALSE,
  				     incVsb, tranAnalysis, intCoeff );
    dIdV->dIsDVsb = contactConductance( pDevice, pSContact, TRUE,
  				     incVsb, tranAnalysis, intCoeff );
    dIdV->dIgDVsb = GateTypeConductance( pDevice, pGContact, FALSE,
  				      incVsb, tranAnalysis, intCoeff );
    dIdV->dIdDVgb = contactConductance( pDevice, pDContact, FALSE,
  				     incVgb, tranAnalysis, intCoeff );
    dIdV->dIsDVgb = contactConductance( pDevice, pSContact, FALSE,
  				     incVgb, tranAnalysis, intCoeff );
    dIdV->dIgDVgb = GateTypeConductance( pDevice, pGContact, TRUE,
  				      incVgb, tranAnalysis, intCoeff );
    
    dIdV->dIdDVdb *= GNorm * width * LNorm;
    dIdV->dIdDVsb *= GNorm * width * LNorm;
    dIdV->dIdDVgb *= GNorm * width * LNorm;
    dIdV->dIsDVdb *= GNorm * width * LNorm;
    dIdV->dIsDVsb *= GNorm * width * LNorm;
    dIdV->dIsDVgb *= GNorm * width * LNorm;
    dIdV->dIgDVdb *= GNorm * width * LNorm;
    dIdV->dIgDVsb *= GNorm * width * LNorm;
    dIdV->dIgDVgb *= GNorm * width * LNorm;
    
  }
  
  double
    contactCurrent(TWOdevice *pDevice, TWOcontact *pContact)
  {
    /* computes the current through the contact given in pContact */
    int index, i, numContactNodes;
    TWOnode *pNode;
    TWOelem *pElem;
    TWOedge *pHEdge, *pVEdge;
    double dx, dy;
    double jTotal = 0.0;
    
    numContactNodes = pContact->numNodes;
    for ( index = 0; index < numContactNodes; index++ ) {
      pNode = pContact->pNodes[ index ];
      for ( i = 0; i <= 3; i++ ) {
        pElem = pNode->pElems[ i ];
        if ( pElem != NULL ) {
  	dx = 0.5 * pElem->dx;
  	dy = 0.5 * pElem->dy;
  	switch ( i ) {
  	case 0:
  	  /* Bottom Right node */
  	  pHEdge = pElem->pBotEdge;
  	  pVEdge = pElem->pRightEdge;
  	  jTotal += pElem->epsRel * ( -dy * pHEdge->jd - dx * pVEdge->jd );
  	  if ( pElem->elemType == SEMICON ) {
  	    jTotal += -dy * (pHEdge->jn + pHEdge->jp)
  	      -dx * (pVEdge->jn + pVEdge->jp);
  	  }
  	  break;
  	case 1:
  	  /* Bottom Left node */
  	  pHEdge = pElem->pBotEdge;
  	  pVEdge = pElem->pLeftEdge;
  	  jTotal += pElem->epsRel * ( dy * pHEdge->jd - dx * pVEdge->jd );
  	  if ( pElem->elemType == SEMICON ) {
  	    jTotal += dy * (pHEdge->jn + pHEdge->jp)
  	      -dx * (pVEdge->jn + pVEdge->jp);
  	  }
  	  break;
  	case 2:
  	  /* Top Left node */
  	  pHEdge = pElem->pTopEdge;
  	  pVEdge = pElem->pLeftEdge;
  	  jTotal += pElem->epsRel * ( dy * pHEdge->jd + dx * pVEdge->jd );
  	  if ( pElem->elemType == SEMICON ) {
  	    jTotal += dy * (pHEdge->jn + pHEdge->jp)
  	      + dx * (pVEdge->jn + pVEdge->jp);
  	  }
  	  break;
  	case 3:
  	  /* Top Right Node */
  	  pHEdge = pElem->pTopEdge;
  	  pVEdge = pElem->pRightEdge;
  	  jTotal += pElem->epsRel * ( -dy * pHEdge->jd + dx * pVEdge->jd );
  	  if ( pElem->elemType == SEMICON ) {
  	    jTotal += -dy * (pHEdge->jn + pHEdge->jp)
  	      + dx * (pVEdge->jn + pVEdge->jp);
  	  }
  	  break;
  	}
        }
      }
    }
    
    return( jTotal * pDevice->width * LNorm * JNorm );
  }
  
  double 
    oxideCurrent(TWOdevice *pDevice, TWOcontact *pContact, 
                 BOOLEAN tranAnalysis)
  {
    /* computes the current through the contact given in pContact */
    int index, i, numContactNodes;
    TWOnode *pNode;
    TWOelem *pElem;
    TWOedge *pHEdge, *pVEdge;
    double dx, dy;
    double jTotal = 0.0;
    
    if ( !tranAnalysis ) {
      return( jTotal );
    }
    
    numContactNodes = pContact->numNodes;
    for ( index = 0; index < numContactNodes; index++ ) {
      pNode = pContact->pNodes[ index ];
      for ( i = 0; i <= 3; i++ ) {
        pElem = pNode->pElems[ i ];
        if ( pElem != NULL ) {
  	dx = 0.5 * pElem->dx;
  	dy = 0.5 * pElem->dy;
  	switch ( i ) {
  	case 0:
  	  /* Bottom Right node */
  	  pHEdge = pElem->pBotEdge;
  	  pVEdge = pElem->pRightEdge;
  	  jTotal += pElem->epsRel * ( -dy * pHEdge->jd - dx * pVEdge->jd );
  	  break;
  	case 1:
  	  /* Bottom Left node */
  	  pHEdge = pElem->pBotEdge;
  	  pVEdge = pElem->pLeftEdge;
  	  jTotal += pElem->epsRel * ( dy * pHEdge->jd - dx * pVEdge->jd );
  	  break;
  	case 2:
  	  /* Top Left node */
  	  pHEdge = pElem->pTopEdge;
  	  pVEdge = pElem->pLeftEdge;
  	  jTotal += pElem->epsRel * ( dy * pHEdge->jd + dx * pVEdge->jd );
  	  break;
  	case 3:
  	  /* Top Right Node */
  	  pHEdge = pElem->pTopEdge;
  	  pVEdge = pElem->pRightEdge;
  	  jTotal += pElem->epsRel * ( -dy * pHEdge->jd + dx * pVEdge->jd );
  	  break;
  	}
        }
      }
    }
    
    return( jTotal * pDevice->width * LNorm * JNorm );
  }
  
  
  double 
     contactConductance(TWOdevice *pDevice, TWOcontact *pContact, 
                        BOOLEAN delVContact, double *dxDv, 
  	              BOOLEAN tranAnalysis, double *intCoeff)
  {
    /* computes the conductance of the contact given in pContact */
    int index, i, numContactNodes;
    TWOnode *pNode, *pHNode = NULL, *pVNode = NULL;
    TWOelem *pElem;
    TWOedge *pHEdge = NULL, *pVEdge = NULL;
    double dPsiDv, dnDv, dpDv;
    double gTotal = 0.0;
    int nInc, pInc;
    
    NG_IGNORE(pDevice);
  
    /* for one carrier the rest of this code relies on appropriate 
       current derivative term to be zero */
    if ( !OneCarrier ) {
      nInc = 1;
      pInc = 2;
    } else {
      nInc = 1;
      pInc = 1;
    } 
    
    numContactNodes = pContact->numNodes;
    for ( index = 0; index < numContactNodes; index++ ) {
      pNode = pContact->pNodes[ index ];
      for ( i = 0; i <= 3; i++ ) {
        pElem = pNode->pElems[ i ];
        if ( pElem != NULL ) {
  	switch ( i ) {
  	case 0:
  	  /* the TL element */
  	  pHNode = pElem->pBLNode;
  	  pVNode = pElem->pTRNode;
  	  pHEdge = pElem->pBotEdge;
  	  pVEdge = pElem->pRightEdge;
  	  if ( pElem->elemType == SEMICON ) {
  	    /* compute the derivatives with n,p */
  	    if ( pHNode->nodeType != CONTACT ) { 
  	      dnDv = dxDv[ pHNode->nEqn ]; 
  	      dpDv = dxDv[ pHNode->pEqn ];
  	      gTotal -= 0.5 * pElem->dy * (pHEdge->dJnDn * dnDv
  					   + pHEdge->dJpDp * dpDv);
  	    }
  	    if ( pVNode->nodeType != CONTACT ) { 
  	      dnDv = dxDv[ pVNode->nEqn ]; 
  	      dpDv = dxDv[ pVNode->pEqn ];
  	      gTotal -= 0.5 * pElem->dx * (pVEdge->dJnDn * dnDv
  					   + pVEdge->dJpDp * dpDv);
  	    }
  	  }
  	  break;
  	case 1:
  	  /* the TR element */
  	  pHNode = pElem->pBRNode;
  	  pVNode = pElem->pTLNode;
  	  pHEdge = pElem->pBotEdge;
  	  pVEdge = pElem->pLeftEdge;
  	  if ( pElem->elemType == SEMICON ) {
  	    /* compute the derivatives with n,p */
  	    if ( pHNode->nodeType != CONTACT ) { 
  	      dnDv = dxDv[ pHNode->nEqn ]; 
  	      dpDv = dxDv[ pHNode->pEqn ];
  	      gTotal += 0.5 * pElem->dy * (pHEdge->dJnDnP1 * dnDv
  					   + pHEdge->dJpDpP1 * dpDv);
  	    }
  	    if ( pVNode->nodeType != CONTACT ) { 
  	      dnDv = dxDv[ pVNode->nEqn ]; 
  	      dpDv = dxDv[ pVNode->pEqn ];
  	      gTotal -= 0.5 * pElem->dx * (pVEdge->dJnDn * dnDv
  					   + pVEdge->dJpDp * dpDv);
  	    }
  	  }
  	  break;
  	case 2:
  	  /* the BR element*/
  	  pHNode = pElem->pTRNode;
  	  pVNode = pElem->pBLNode;
  	  pHEdge = pElem->pTopEdge;
  	  pVEdge = pElem->pLeftEdge;
  	  if ( pElem->elemType == SEMICON ) {
  	    /* compute the derivatives with n,p */
  	    if ( pHNode->nodeType != CONTACT ) { 
  	      dnDv = dxDv[ pHNode->nEqn ]; 
  	      dpDv = dxDv[ pHNode->pEqn ];
  	      gTotal += 0.5 * pElem->dy * (pHEdge->dJnDnP1 * dnDv
  					   + pHEdge->dJpDpP1 * dpDv);
  	    }
  	    if ( pVNode->nodeType != CONTACT ) { 
  	      dnDv = dxDv[ pVNode->nEqn ]; 
  	      dpDv = dxDv[ pVNode->pEqn ];
  	      gTotal += 0.5 * pElem->dx * (pVEdge->dJnDnP1 * dnDv
  					   + pVEdge->dJpDpP1 * dpDv);
  	    }
  	  }
  	  break;
  	case 3:
  	  /* the BL element */
  	  pHNode = pElem->pTLNode;
  	  pVNode = pElem->pBRNode;
  	  pHEdge = pElem->pTopEdge;
  	  pVEdge = pElem->pRightEdge;
  	  if ( pElem->elemType == SEMICON ) {
  	    /* compute the derivatives with n,p */
  	    if ( pHNode->nodeType != CONTACT ) { 
  	      dnDv = dxDv[ pHNode->nEqn ]; 
  	      dpDv = dxDv[ pHNode->pEqn ];
  	      gTotal -= 0.5 * pElem->dy * (pHEdge->dJnDn * dnDv
  					   + pHEdge->dJpDp * dpDv);
  	    }
  	    if ( pVNode->nodeType != CONTACT ) { 
  	      dnDv = dxDv[ pVNode->nEqn ]; 
  	      dpDv = dxDv[ pVNode->pEqn ];
  	      gTotal += 0.5 * pElem->dx * (pVEdge->dJnDnP1 * dnDv
  					   + pVEdge->dJpDpP1 * dpDv);
  	    }
  	  }
  	  break;
  	}
  	if ( pElem->elemType == SEMICON ) {
  	  if ( pHNode->nodeType != CONTACT ) { 
  	    dPsiDv = dxDv[ pHNode->psiEqn ];
  	    gTotal += 0.5 * pElem->dy * dPsiDv * (pHEdge->dJnDpsiP1 + pHEdge->dJpDpsiP1 ); 
  	    if ( delVContact ) {
  	      gTotal -= 0.5 * pElem->dy * (pHEdge->dJnDpsiP1 + pHEdge->dJpDpsiP1 );
  	    }
  	  }
  	  if ( pVNode->nodeType != CONTACT ) { 
  	    dPsiDv = dxDv[ pVNode->psiEqn ];
  	    gTotal += 0.5 * pElem->dx * dPsiDv * (pVEdge->dJnDpsiP1 + pVEdge->dJpDpsiP1 ); 
  	    if ( delVContact ) {
  	      gTotal -= 0.5 * pElem->dx * (pVEdge->dJnDpsiP1 + pVEdge->dJpDpsiP1 );
  	    }
  	  }
  	}
  	if ( tranAnalysis ) {
  	  /* add the displacement current terms */
  	  if ( pHNode->nodeType != CONTACT ) { 
  	    dPsiDv = dxDv[ pHNode->psiEqn ];
  	    gTotal -= intCoeff[0] * pElem->epsRel * 0.5 * pElem->dyOverDx * dPsiDv;
  	    if ( delVContact ) {
  	      gTotal += intCoeff[0] * pElem->epsRel * 0.5 * pElem->dyOverDx;
  	    }
  	  }
  	  if ( pVNode->nodeType != CONTACT ) { 
  	    dPsiDv = dxDv[ pVNode->psiEqn ];
  	    gTotal -= intCoeff[0] * pElem->epsRel * 0.5 * pElem->dxOverDy * dPsiDv;
  	    if ( delVContact ) {
  	      gTotal += intCoeff[0] * pElem->epsRel * 0.5 * pElem->dxOverDy;
  	    }
  	  }
  	}
        }
      }
    }
    
    return( gTotal );
  }
  
  
  double 
     oxideConductance(TWOdevice *pDevice, TWOcontact *pContact, 
                      BOOLEAN delVContact, double *dxDv, 
  	            BOOLEAN tranAnalysis, double *intCoeff)
  {
    /* computes the conductance of the contact given in pContact */
    int index, i, numContactNodes;
    TWOnode *pNode, *pHNode = NULL, *pVNode = NULL;
    TWOelem *pElem;
    double dPsiDv;
    double gTotal = 0.0;
    
    NG_IGNORE(pDevice);
  
    if ( !tranAnalysis ) {
      return( gTotal );
    }
    
    numContactNodes = pContact->numNodes;
    for ( index = 0; index < numContactNodes; index++ ) {
      pNode = pContact->pNodes[ index ];
      for ( i = 0; i <= 3; i++ ) {
        pElem = pNode->pElems[ i ];
        if ( pElem != NULL ) {
  	switch ( i ) {
  	case 0:
  	  /* the TL element */
  	  pHNode = pElem->pBLNode;
  	  pVNode = pElem->pTRNode;
  	  break;
  	case 1:
  	  /* the TR element */
  	  pHNode = pElem->pBRNode;
  	  pVNode = pElem->pTLNode;
  	  break;
  	case 2:
  	  /* the BR element*/
  	  pHNode = pElem->pTRNode;
  	  pVNode = pElem->pBLNode;
  	  break;
  	case 3:
  	  /* the BL element */
  	  pHNode = pElem->pTLNode;
  	  pVNode = pElem->pBRNode;
  	  break;
  	}
  	/* add the displacement current terms */
  	if ( pHNode->nodeType != CONTACT ) { 
  	  dPsiDv = dxDv[ pHNode->psiEqn ];
  	  gTotal -= intCoeff[0] * pElem->epsRel * 0.5 * pElem->dyOverDx * dPsiDv;
  	  if ( delVContact ) {
  	    gTotal += intCoeff[0] * pElem->epsRel * 0.5 * pElem->dyOverDx;
  	  }
  	}
  	if ( pVNode->nodeType != CONTACT ) { 
  	  dPsiDv = dxDv[ pVNode->psiEqn ];
  	  gTotal -= intCoeff[0] * pElem->epsRel * 0.5 * pElem->dxOverDy * dPsiDv;
  	  if ( delVContact ) {
  	    gTotal += intCoeff[0] * pElem->epsRel * 0.5 * pElem->dxOverDy;
  	  }
  	}
        }
      }
    }
    
    return( gTotal );
  }
  
  /* these functions are used for solving the complete system of 
   * equations directly using LU decomposition   1/22/88
   */
  
  void 
     NUMD2current(TWOdevice *pDevice, BOOLEAN tranAnalysis, 
                  double *intCoeff, double *id)
  {
    TWOcontact *pPContact = pDevice->pFirstContact;
  /*  TWOcontact *pNContact = pDevice->pLastContact; */
    double ip, ipPrime, *solution;
  /*  double in;*/
    BOOLEAN deltaVContact = FALSE;
    
    solution = pDevice->dcDeltaSolution;
    ip = contactCurrent( pDevice, pPContact );
    /*
    in = contactCurrent( pDevice, pNContact );
    fprintf(stdout, "DIO current: ( %11.4e error )\n", ip+in );
    fprintf(stdout, "     Ip = %11.4e     In = %11.4e\n", ip, in );
    */
  
    /* 
     * for the additional contribution to id will make use of 
     * contactConductance. This function will be called
     * with the dcDeltaSolution vector instead of the incremental quantities
     */
    ipPrime = contactConductance( pDevice, pPContact, deltaVContact,
  			       solution, tranAnalysis, intCoeff );
    
    ipPrime *= JNorm * pDevice->width * LNorm;
    ip += ipPrime;
    *id = ip;
  }
  
  void 
     NBJT2current(TWOdevice *pDevice, BOOLEAN tranAnalysis, double *intCoeff, 
                  double *ie, double *ic)
  {
    TWOcontact *pEmitContact = pDevice->pLastContact;
    TWOcontact *pColContact = pDevice->pFirstContact;
  /*  TWOcontact *pBaseContact = pDevice->pFirstContact->next; */
    double *solution, iePrime, icPrime;
  /*  double ib; */
    
    solution = pDevice->dcDeltaSolution;
    
    *ie = contactCurrent( pDevice, pEmitContact );
    *ic = contactCurrent( pDevice, pColContact );
    /*
    ib = contactCurrent( pDevice, pBaseContact );
    fprintf(stdout, "BJT current: ( %11.4e error )\n", *ic+ib+*ie );
    fprintf(stdout, "     Ic = %11.4e     Ib = %11.4e\n", *ic, ib );
    fprintf(stdout, "     Ie = %11.4e\n", *ie );
    */
    
    iePrime = contactConductance( pDevice, pEmitContact, FALSE, solution,
  			       tranAnalysis, intCoeff );
    
    icPrime = contactConductance( pDevice, pColContact, FALSE, solution,
  			       tranAnalysis, intCoeff );
    
    iePrime *= JNorm * pDevice->width * LNorm;
    icPrime *= JNorm * pDevice->width * LNorm;
    
    *ie += iePrime;
    *ic += icPrime;
  }
  
  void 
     NUMOScurrent(TWOdevice *pDevice, BOOLEAN tranAnalysis, double *intCoeff, 
                  double *id, double *is, double *ig)
  {
    TWOcontact *pDContact = pDevice->pFirstContact;
    TWOcontact *pGContact = pDevice->pFirstContact->next;
    TWOcontact *pSContact = pDevice->pFirstContact->next->next;
  /*  TWOcontact *pBContact = pDevice->pLastContact; */
    double *solution, idPrime, isPrime, igPrime;
  /*  double ib; */
    
    solution = pDevice->dcDeltaSolution;
    
    *id = contactCurrent( pDevice, pDContact );
  
  /* 
   * This is a terrible hack
   */  
    
  #ifdef NORMAL_GATE  
    *ig = GateTypeCurrent( pDevice, pGContact, tranAnalysis );
  #else
    *ig = GateTypeCurrent( pDevice, pGContact);
  #endif  
    
    *is = contactCurrent( pDevice, pSContact );
    /*
    ib = contactCurrent( pDevice, pBContact );
    fprintf(stdout, "MOS current: ( %11.4e error )\n", *id+*ig+*is+ib );
    fprintf(stdout, "     Id = %11.4e     Is = %11.4e\n", *id, *is );
    fprintf(stdout, "     Ig = %11.4e     Ib = %11.4e\n", *ig, ib );
    */
  
    idPrime = contactConductance( pDevice, pDContact, FALSE,
  			       solution, tranAnalysis, intCoeff );
    
    isPrime = contactConductance( pDevice, pSContact, FALSE,
  			       solution, tranAnalysis, intCoeff );
    
    igPrime = GateTypeConductance( pDevice, pGContact, FALSE,
  				solution, tranAnalysis, intCoeff );
    
    idPrime *= JNorm * pDevice->width * LNorm;
    isPrime *= JNorm * pDevice->width * LNorm;
    igPrime *= JNorm * pDevice->width * LNorm;
    
    *id += idPrime;
    *is += isPrime;
    *ig += igPrime;
  }
  