twomobdv.c

spice中支持多层次元件模型仿真的可单独运行的插件源码

/**********
Copyright 1991 Regents of the University of California.  All rights reserved.
Author:	1991 David A. Gates, U. C. Berkeley CAD Group
**********/

#include <math.h>
#include "nummacs.h"
#include "numglobs.h"
#include "twomesh.h"

/*
 * Load the derivatives of the current with respect to changes in
 * the mobility, for all edges of a silicon element.
 * It is known a priori that the element is made of semiconductor.
 * These routines work for both channel and bulk elements.
 */

void
  TWO_mobDeriv( pElem, chanType, ds )
TWOelem *pElem;			/* channel or bulk element */
int chanType;			/* flag for direction of channel */
double ds;			/* normalized hgt (len) of interface element */
{
  TWOnode *pNode;
  TWOedge *pTEdge, *pBEdge, *pLEdge, *pREdge;
  BOOLEAN channel = pElem->channel;
  double dx, dy, rDx, rDy;
  double coeffHx, coeffHy, coeffHs;
  double coeffVx, coeffVy, coeffVs;
  double dFnxDMun, dFnyDMun, dFnDEs;
  double dFpxDMup, dFpyDMup, dFpDEs;
  double dMnDEs, dMnDEx, dMnDEy, dMnDWx, dMnDWy;
  double dMpDEs, dMpDEx, dMpDEy, dMpDWx, dMpDWy;
  
  /* Initialize various quantities */
  dx = pElem->dx;
  dy = pElem->dy;
  rDx = 1.0 / dx;
  rDy = 1.0 / dy;
  
  /* compute length-dependent parameters */
  coeffHx = 0.25 * dy * rDx;	/* For horizontal edges */
  coeffHy = 0.25;
  coeffVx = 0.25;		/* For vertical edges */
  coeffVy = 0.25 * dx * rDy;
  switch ( chanType ) {
  case 0:
  case 3:
    coeffHs = 0.25 * dy / ds;
    coeffVs = 0.25 * dx / ds;
    break;
  case 1:
  case 2:
    coeffHs = - 0.25 * dy / ds;
    coeffVs = - 0.25 * dx / ds;
    break;
  }
  
  /* Get pointers to element's edges */
  pTEdge = pElem->pTopEdge;
  pBEdge = pElem->pBotEdge;
  pLEdge = pElem->pLeftEdge;
  pREdge = pElem->pRightEdge;
  
  /* Get element mobility derivatives for fast access later */
  dMnDEs = pElem->dMunDEs;
  dMnDEx = pElem->dMunDEx;
  dMnDEy = pElem->dMunDEy;
  dMnDWx = pElem->dMunDWx;
  dMnDWy = pElem->dMunDWy;
  dMpDEs = pElem->dMupDEs;
  dMpDEx = pElem->dMupDEx;
  dMpDEy = pElem->dMupDEy;
  dMpDWx = pElem->dMupDWx;
  dMpDWy = pElem->dMupDWy;
  
  /* Add mobility derivatives due to Top Edge */
  /* First compute derivatives of cont. eqn's */ 
  dFnxDMun = coeffHx * (pTEdge->wdfn * rDx);
  dFnyDMun = coeffHy * (pTEdge->wdfn * rDx);
  dFpxDMup = coeffHx * (pTEdge->wdfp * rDx);
  dFpyDMup = coeffHy * (pTEdge->wdfp * rDx);
  
  /* Do Top-Left (TL) Node of Element */
  pNode = pElem->pTLNode;
  /* n continuity wrto potential derivatives */
  *(pNode->fNPsi) +=
    dFnxDMun * ( dMnDEx - dMnDWx * pTEdge->dWnDpsiP1 )
      + dFnyDMun * ( dMnDEy - dMnDWy * pLEdge->dWnDpsiP1 );
  *(pNode->fNPsiiP1) +=
    - dFnxDMun * ( dMnDEx - dMnDWx * pTEdge->dWnDpsiP1 )
      + dFnyDMun * ( dMnDEy - dMnDWy * pREdge->dWnDpsiP1 );
  *(pNode->fNPsiiP1jP1) +=
    - dFnxDMun * ( dMnDEx - dMnDWx * pBEdge->dWnDpsiP1 )
      - dFnyDMun * ( dMnDEy - dMnDWy * pREdge->dWnDpsiP1 );
  *(pNode->fNPsijP1) +=
    dFnxDMun * ( dMnDEx - dMnDWx * pBEdge->dWnDpsiP1 )
      - dFnyDMun * ( dMnDEy - dMnDWy * pLEdge->dWnDpsiP1 );
  /* p continuity wrto potential derivatives */
  *(pNode->fPPsi) +=
    dFpxDMup * ( dMpDEx - dMpDWx * pTEdge->dWpDpsiP1 )
      + dFpyDMup * ( dMpDEy - dMpDWy * pLEdge->dWpDpsiP1 );
  *(pNode->fPPsiiP1) +=
    - dFpxDMup * ( dMpDEx - dMpDWx * pTEdge->dWpDpsiP1 )
      + dFpyDMup * ( dMpDEy - dMpDWy * pREdge->dWpDpsiP1 );
  *(pNode->fPPsiiP1jP1) +=
    - dFpxDMup * ( dMpDEx - dMpDWx * pBEdge->dWpDpsiP1 )
      - dFpyDMup * ( dMpDEy - dMpDWy * pREdge->dWpDpsiP1 );
  *(pNode->fPPsijP1) +=
    dFpxDMup * ( dMpDEx - dMpDWx * pBEdge->dWpDpsiP1 )
      - dFpyDMup * ( dMpDEy - dMpDWy * pLEdge->dWpDpsiP1 );
  /* n continuity wrto n derivatives */
  *(pNode->fNN) +=
    dFnxDMun * ( dMnDWx * pTEdge->dWnDn )
      + dFnyDMun * ( dMnDWy * pLEdge->dWnDn );
  *(pNode->fNNiP1) +=
    dFnxDMun * ( dMnDWx * pTEdge->dWnDnP1 )
      + dFnyDMun * ( dMnDWy * pREdge->dWnDn );
  *(pNode->fNNiP1jP1) +=
    dFnxDMun * ( dMnDWx * pBEdge->dWnDnP1 )
      + dFnyDMun * ( dMnDWy * pREdge->dWnDnP1 );
  *(pNode->fNNjP1) +=
    dFnxDMun * ( dMnDWx * pBEdge->dWnDn )
      + dFnyDMun * ( dMnDWy * pLEdge->dWnDnP1 );
  /* p continuity wrto p derivatives */
  *(pNode->fPP) +=
    dFpxDMup * ( dMpDWx * pTEdge->dWpDp )
      + dFpyDMup * ( dMpDWy * pLEdge->dWpDp );
  *(pNode->fPPiP1) +=
    dFpxDMup * ( dMpDWx * pTEdge->dWpDpP1 )
      + dFpyDMup * ( dMpDWy * pREdge->dWpDp );
  *(pNode->fPPiP1jP1) +=
    dFpxDMup * ( dMpDWx * pBEdge->dWpDpP1 )
      + dFpyDMup * ( dMpDWy * pREdge->dWpDpP1 );
  *(pNode->fPPjP1) +=
    dFpxDMup * ( dMpDWx * pBEdge->dWpDp )
      + dFpyDMup * ( dMpDWy * pLEdge->dWpDpP1 );
  
  /* both continuity wrto surface potential derivatives */
  if ( channel ) {
    dFnDEs   = coeffHs * (pTEdge->wdfn * rDx) * dMnDEs;
    dFpDEs   = coeffHs * (pTEdge->wdfp * rDx) * dMpDEs;
    *(pNode->fNPsiIn)   -= dFnDEs;
    *(pNode->fNPsiInP1) -= dFnDEs;
    *(pNode->fNPsiOx)   += dFnDEs;
    *(pNode->fNPsiOxP1) += dFnDEs;
    *(pNode->fPPsiIn)   -= dFpDEs;
    *(pNode->fPPsiInP1) -= dFpDEs;
    *(pNode->fPPsiOx)   += dFpDEs;
    *(pNode->fPPsiOxP1) += dFpDEs;
  }
  
  /* Do Top-Right (TR) Node of Element */
  pNode = pElem->pTRNode;
  /* n continuity wrto potential derivatives */
  *(pNode->fNPsiiM1) -=
    dFnxDMun * ( dMnDEx - dMnDWx * pTEdge->dWnDpsiP1 )
      + dFnyDMun * ( dMnDEy - dMnDWy * pLEdge->dWnDpsiP1 );
  *(pNode->fNPsi) -=
    - dFnxDMun * ( dMnDEx - dMnDWx * pTEdge->dWnDpsiP1 )
      + dFnyDMun * ( dMnDEy - dMnDWy * pREdge->dWnDpsiP1 );
  *(pNode->fNPsijP1) -=
    - dFnxDMun * ( dMnDEx - dMnDWx * pBEdge->dWnDpsiP1 )
      - dFnyDMun * ( dMnDEy - dMnDWy * pREdge->dWnDpsiP1 );
  *(pNode->fNPsiiM1jP1) -=
    dFnxDMun * ( dMnDEx - dMnDWx * pBEdge->dWnDpsiP1 )
      - dFnyDMun * ( dMnDEy - dMnDWy * pLEdge->dWnDpsiP1 );
  /* p continuity wrto potential derivatives */
  *(pNode->fPPsiiM1) -=
    dFpxDMup * ( dMpDEx - dMpDWx * pTEdge->dWpDpsiP1 )
      + dFpyDMup * ( dMpDEy - dMpDWy * pLEdge->dWpDpsiP1 );
  *(pNode->fPPsi) -=
    - dFpxDMup * ( dMpDEx - dMpDWx * pTEdge->dWpDpsiP1 )
      + dFpyDMup * ( dMpDEy - dMpDWy * pREdge->dWpDpsiP1 );
  *(pNode->fPPsijP1) -=
    - dFpxDMup * ( dMpDEx - dMpDWx * pBEdge->dWpDpsiP1 )
      - dFpyDMup * ( dMpDEy - dMpDWy * pREdge->dWpDpsiP1 );
  *(pNode->fPPsiiM1jP1) -=
    dFpxDMup * ( dMpDEx - dMpDWx * pBEdge->dWpDpsiP1 )
      - dFpyDMup * ( dMpDEy - dMpDWy * pLEdge->dWpDpsiP1 );
  /* n continuity wrto n derivatives */
  *(pNode->fNNiM1) -=
    dFnxDMun * ( dMnDWx * pTEdge->dWnDn )
      + dFnyDMun * ( dMnDWy * pLEdge->dWnDn );
  *(pNode->fNN) -=
    dFnxDMun * ( dMnDWx * pTEdge->dWnDnP1 )
      + dFnyDMun * ( dMnDWy * pREdge->dWnDn );
  *(pNode->fNNjP1) -=
    dFnxDMun * ( dMnDWx * pBEdge->dWnDnP1 )
      + dFnyDMun * ( dMnDWy * pREdge->dWnDnP1 );
  *(pNode->fNNiM1jP1) -=
    dFnxDMun * ( dMnDWx * pBEdge->dWnDn )
      + dFnyDMun * ( dMnDWy * pLEdge->dWnDnP1 );
  /* p continuity wrto p derivatives */
  *(pNode->fPPiM1) -=
    dFpxDMup * ( dMpDWx * pTEdge->dWpDp )
      + dFpyDMup * ( dMpDWy * pLEdge->dWpDp );
  *(pNode->fPP) -=
    dFpxDMup * ( dMpDWx * pTEdge->dWpDpP1 )
      + dFpyDMup * ( dMpDWy * pREdge->dWpDp );
  *(pNode->fPPjP1) -=
    dFpxDMup * ( dMpDWx * pBEdge->dWpDpP1 )
      + dFpyDMup * ( dMpDWy * pREdge->dWpDpP1 );
  *(pNode->fPPiM1jP1) -=
    dFpxDMup * ( dMpDWx * pBEdge->dWpDp )
      + dFpyDMup * ( dMpDWy * pLEdge->dWpDpP1 );
  
  /* both continuity wrto surface potential derivatives */
  if ( channel ) {
    *(pNode->fNPsiInM1) += dFnDEs;
    *(pNode->fNPsiIn)   += dFnDEs;
    *(pNode->fNPsiOxM1) -= dFnDEs;
    *(pNode->fNPsiOx)   -= dFnDEs;
    *(pNode->fPPsiInM1) += dFpDEs;
    *(pNode->fPPsiIn)   += dFpDEs;
    *(pNode->fPPsiOxM1) -= dFpDEs;
    *(pNode->fPPsiOx)   -= dFpDEs;
  }
  
  
  /* Add mobility derivatives due to Bottom Edge */
  /* First compute derivatives of cont. eqn's */ 
  dFnxDMun = coeffHx * (pBEdge->wdfn * rDx);
  dFnyDMun = coeffHy * (pBEdge->wdfn * rDx);
  dFpxDMup = coeffHx * (pBEdge->wdfp * rDx);
  dFpyDMup = coeffHy * (pBEdge->wdfp * rDx);
  
  /* Do Bottom-Left (BL) Node of Element */
  pNode = pElem->pBLNode;
  /* n continuity wrto potential derivatives */
  *(pNode->fNPsijM1) +=
    dFnxDMun * ( dMnDEx - dMnDWx * pTEdge->dWnDpsiP1 )
      + dFnyDMun * ( dMnDEy - dMnDWy * pLEdge->dWnDpsiP1 );
  *(pNode->fNPsiiP1jM1) +=
    - dFnxDMun * ( dMnDEx - dMnDWx * pTEdge->dWnDpsiP1 )
      + dFnyDMun * ( dMnDEy - dMnDWy * pREdge->dWnDpsiP1 );
  *(pNode->fNPsiiP1) +=
    - dFnxDMun * ( dMnDEx - dMnDWx * pBEdge->dWnDpsiP1 )
      - dFnyDMun * ( dMnDEy - dMnDWy * pREdge->dWnDpsiP1 );
  *(pNode->fNPsi) +=
    dFnxDMun * ( dMnDEx - dMnDWx * pBEdge->dWnDpsiP1 )
      - dFnyDMun * ( dMnDEy - dMnDWy * pLEdge->dWnDpsiP1 );
  /* p continuity wrto potential derivatives */
  *(pNode->fPPsijM1) +=
    dFpxDMup * ( dMpDEx - dMpDWx * pTEdge->dWpDpsiP1 )
      + dFpyDMup * ( dMpDEy - dMpDWy * pLEdge->dWpDpsiP1 );
  *(pNode->fPPsiiP1jM1) +=
    - dFpxDMup * ( dMpDEx - dMpDWx * pTEdge->dWpDpsiP1 )
      + dFpyDMup * ( dMpDEy - dMpDWy * pREdge->dWpDpsiP1 );
  *(pNode->fPPsiiP1) +=
    - dFpxDMup * ( dMpDEx - dMpDWx * pBEdge->dWpDpsiP1 )
      - dFpyDMup * ( dMpDEy - dMpDWy * pREdge->dWpDpsiP1 );
  *(pNode->fPPsi) +=
    dFpxDMup * ( dMpDEx - dMpDWx * pBEdge->dWpDpsiP1 )
      - dFpyDMup * ( dMpDEy - dMpDWy * pLEdge->dWpDpsiP1 );
  /* n continuity wrto n derivatives */
  *(pNode->fNNjM1) +=
    dFnxDMun * ( dMnDWx * pTEdge->dWnDn )
      + dFnyDMun * ( dMnDWy * pLEdge->dWnDn );
  *(pNode->fNNiP1jM1) +=
    dFnxDMun * ( dMnDWx * pTEdge->dWnDnP1 )
      + dFnyDMun * ( dMnDWy * pREdge->dWnDn );
  *(pNode->fNNiP1) +=
    dFnxDMun * ( dMnDWx * pBEdge->dWnDnP1 )
      + dFnyDMun * ( dMnDWy * pREdge->dWnDnP1 );
  *(pNode->fNN) +=
    dFnxDMun * ( dMnDWx * pBEdge->dWnDn )
      + dFnyDMun * ( dMnDWy * pLEdge->dWnDnP1 );
  /* p continuity wrto p derivatives */
  *(pNode->fPPjM1) +=
    dFpxDMup * ( dMpDWx * pTEdge->dWpDp )
      + dFpyDMup * ( dMpDWy * pLEdge->dWpDp );
  *(pNode->fPPiP1jM1) +=
    dFpxDMup * ( dMpDWx * pTEdge->dWpDpP1 )
      + dFpyDMup * ( dMpDWy * pREdge->dWpDp );
  *(pNode->fPPiP1) +=
    dFpxDMup * ( dMpDWx * pBEdge->dWpDpP1 )
      + dFpyDMup * ( dMpDWy * pREdge->dWpDpP1 );
  *(pNode->fPP) +=
    dFpxDMup * ( dMpDWx * pBEdge->dWpDp )
      + dFpyDMup * ( dMpDWy * pLEdge->dWpDpP1 );
  
  /* both continuity wrto surface potential derivatives */
  if (channel ) {
    dFnDEs   = coeffHs * (pBEdge->wdfn * rDx) * dMnDEs;
    dFpDEs   = coeffHs * (pBEdge->wdfp * rDx) * dMpDEs;
    *(pNode->fNPsiIn)   -= dFnDEs;
    *(pNode->fNPsiInP1) -= dFnDEs;
    *(pNode->fNPsiOx)   += dFnDEs;
    *(pNode->fNPsiOxP1) += dFnDEs;
    *(pNode->fPPsiIn)   -= dFpDEs;
    *(pNode->fPPsiInP1) -= dFpDEs;
    *(pNode->fPPsiOx)   += dFpDEs;
    *(pNode->fPPsiOxP1) += dFpDEs;
  }
  
  /* Do Bottom-Right (BR) Node of Element */
  pNode = pElem->pBRNode;
  /* n continuity wrto potential derivatives */
  *(pNode->fNPsiiM1jM1) -=
    dFnxDMun * ( dMnDEx - dMnDWx * pTEdge->dWnDpsiP1 )
      + dFnyDMun * ( dMnDEy - dMnDWy * pLEdge->dWnDpsiP1 );
  *(pNode->fNPsijM1) -=
    - dFnxDMun * ( dMnDEx - dMnDWx * pTEdge->dWnDpsiP1 )
      + dFnyDMun * ( dMnDEy - dMnDWy * pREdge->dWnDpsiP1 );
  *(pNode->fNPsi) -=
    - dFnxDMun * ( dMnDEx - dMnDWx * pBEdge->dWnDpsiP1 )
      - dFnyDMun * ( dMnDEy - dMnDWy * pREdge->dWnDpsiP1 );
  *(pNode->fNPsiiM1) -=
    dFnxDMun * ( dMnDEx - dMnDWx * pBEdge->dWnDpsiP1 )
      - dFnyDMun * ( dMnDEy - dMnDWy * pLEdge->dWnDpsiP1 );
  /* p continuity wrto potential derivatives */
  *(pNode->fPPsiiM1jM1) -=
    dFpxDMup * ( dMpDEx - dMpDWx * pTEdge->dWpDpsiP1 )
      + dFpyDMup * ( dMpDEy - dMpDWy * pLEdge->dWpDpsiP1 );
  *(pNode->fPPsijM1) -=
    - dFpxDMup * ( dMpDEx - dMpDWx * pTEdge->dWpDpsiP1 )
      + dFpyDMup * ( dMpDEy - dMpDWy * pREdge->dWpDpsiP1 );
  *(pNode->fPPsi) -=
    - dFpxDMup * ( dMpDEx - dMpDWx * pBEdge->dWpDpsiP1 )
      - dFpyDMup * ( dMpDEy - dMpDWy * pREdge->dWpDpsiP1 );
  *(pNode->fPPsiiM1) -=
    dFpxDMup * ( dMpDEx - dMpDWx * pBEdge->dWpDpsiP1 )
      - dFpyDMup * ( dMpDEy - dMpDWy * pLEdge->dWpDpsiP1 );
  /* n continuity wrto n derivatives */
  *(pNode->fNNiM1jM1) -=
    dFnxDMun * ( dMnDWx * pTEdge->dWnDn )
      + dFnyDMun * ( dMnDWy * pLEdge->dWnDn );
  *(pNode->fNNjM1) -=
    dFnxDMun * ( dMnDWx * pTEdge->dWnDnP1 )
      + dFnyDMun * ( dMnDWy * pREdge->dWnDn );
  *(pNode->fNN) -=
    dFnxDMun * ( dMnDWx * pBEdge->dWnDnP1 )
      + dFnyDMun * ( dMnDWy * pREdge->dWnDnP1 );
  *(pNode->fNNiM1) -=
    dFnxDMun * ( dMnDWx * pBEdge->dWnDn )
      + dFnyDMun * ( dMnDWy * pLEdge->dWnDnP1 );
  /* p continuity wrto p derivatives */
  *(pNode->fPPiM1jM1) -=
    dFpxDMup * ( dMpDWx * pTEdge->dWpDp )
      + dFpyDMup * ( dMpDWy * pLEdge->dWpDp );
  *(pNode->fPPjM1) -=
    dFpxDMup * ( dMpDWx * pTEdge->dWpDpP1 )
      + dFpyDMup * ( dMpDWy * pREdge->dWpDp );
  *(pNode->fPP) -=
    dFpxDMup * ( dMpDWx * pBEdge->dWpDpP1 )
      + dFpyDMup * ( dMpDWy * pREdge->dWpDpP1 );
  *(pNode->fPPiM1) -=
    dFpxDMup * ( dMpDWx * pBEdge->dWpDp )
      + dFpyDMup * ( dMpDWy * pLEdge->dWpDpP1 );