twopcont.c

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

    pLEdge = pElem->pLeftEdge;
    pREdge = pElem->pRightEdge;
    
    /* load for all i,j */
    for ( index = 0; index <= 3; index++ ) {
      pNode = pElem->pNodes[ index ];
      if ( pNode->nodeType ISNOT CONTACT ) {
	*(pNode->fPsiPsi) += dyOverDx + dxOverDy;
	if ( pElem->elemType IS SEMICON ) {
	  if ( index <= 1 ) {
	    pHEdge = pTEdge;
	  } else {
	    pHEdge = pBEdge;
	  }
	  if ( index IS 0 OR index IS 3 ) {
	    pVEdge = pLEdge;
	  } else {
	    pVEdge = pREdge;
	  }
	  nConc = *(pDevice->devState0 + pNode->nodeN);
	  *(pNode->fPsiPsi) += dxdy * nConc;
	  *(pNode->fPsiP) -= dxdy;
	  *(pNode->fPPsi) -= dy * pHEdge->dJpDpsiP1 + dx * pVEdge->dJpDpsiP1;
	  
	  /* Handle generation terms */
	  *(pNode->fPP) += dxdy * pNode->dUdP;
	  *(pNode->fPPsi) += dxdy * pNode->dUdN * nConc;
	}
      }
    }
    
    /* Handle neighbor and edge dependent terms */
    pNode = pElem->pTLNode;
    if ( pNode->nodeType ISNOT CONTACT ) {
      *(pNode->fPsiPsiiP1) -= dyOverDx;
      *(pNode->fPsiPsijP1) -= dxOverDy;
      if ( pElem->elemType IS SEMICON ) {
	*(pNode->fPP) += dy * pTEdge->dJpDp + dx * pLEdge->dJpDp;
	*(pNode->fPPsiiP1) += dy * pTEdge->dJpDpsiP1;
	*(pNode->fPPiP1) += dy * pTEdge->dJpDpP1;
	*(pNode->fPPsijP1) += dx * pLEdge->dJpDpsiP1;
	*(pNode->fPPjP1) += dx * pLEdge->dJpDpP1;
      }
    }
    pNode = pElem->pTRNode;
    if ( pNode->nodeType ISNOT CONTACT ) {
      *(pNode->fPsiPsiiM1) -= dyOverDx;
      *(pNode->fPsiPsijP1) -= dxOverDy;
      if ( pElem->elemType IS SEMICON ) {
	*(pNode->fPP) += -dy * pTEdge->dJpDpP1 + dx * pREdge->dJpDp;
	*(pNode->fPPsiiM1) += dy * pTEdge->dJpDpsiP1;
	*(pNode->fPPiM1) -= dy * pTEdge->dJpDp;
	*(pNode->fPPsijP1) += dx * pREdge->dJpDpsiP1;
	*(pNode->fPPjP1) += dx * pREdge->dJpDpP1;
      }
    }
    pNode = pElem->pBRNode;
    if ( pNode->nodeType ISNOT CONTACT ) {
      *(pNode->fPsiPsiiM1) -= dyOverDx;
      *(pNode->fPsiPsijM1) -= dxOverDy;
      if ( pElem->elemType IS SEMICON ) {
	*(pNode->fPP) += -dy * pBEdge->dJpDpP1 - dx * pREdge->dJpDpP1;
	*(pNode->fPPsiiM1) += dy * pBEdge->dJpDpsiP1;
	*(pNode->fPPiM1) -= dy * pBEdge->dJpDp;
	*(pNode->fPPsijM1) += dx * pREdge->dJpDpsiP1;
	*(pNode->fPPjM1) -= dx * pREdge->dJpDp;
      }
    }
    pNode = pElem->pBLNode;
    if ( pNode->nodeType ISNOT CONTACT ) {
      *(pNode->fPsiPsiiP1) -= dyOverDx;
      *(pNode->fPsiPsijM1) -= dxOverDy;
      if ( pElem->elemType IS SEMICON ) {
	*(pNode->fPP) += dy * pBEdge->dJpDp - dx * pLEdge->dJpDpP1;
	*(pNode->fPPsiiP1) += dy * pBEdge->dJpDpsiP1;
	*(pNode->fPPiP1) += dy * pBEdge->dJpDpP1;
	*(pNode->fPPsijM1) += dx * pLEdge->dJpDpsiP1;
	*(pNode->fPPjM1) -= dx * pLEdge->dJpDp;
      }
    }
  }

  /* Calculate the Inversion-Layer Mobility Dependent Terms in Jac. */
  if ( MobDeriv AND SurfaceMobility ) {
    for ( pCh = pDevice->pChannel; pCh ISNOT NIL(TWOchannel);
	 pCh = pCh->next ) {
      /* Find effective height of oxide element at interface. */
      if ( pCh->type%2 == 0 ) { /* Vertical slice */
	ds = pCh->pNElem->dy / pCh->pNElem->epsRel;
      } else {			/* Horizontal slice */
	ds = pCh->pNElem->dx / pCh->pNElem->epsRel;
      }
      pElem = pCh->pSeed;
      nextIndex = (pCh->type + 2)%4;
      while (pElem && pElem->channel == pCh->id) {  
	TWOPmobDeriv( pElem, pCh->type, ds );
	pElem = pElem->pElems[ nextIndex ];
      }
    } /* endfor pCh ISNOT NIL */
  } /* endif MobDeriv and SurfaceMobility */
}

/* load only the Rhs vector */
void 
  TWOPrhsLoad( pDevice, tranAnalysis, info ) 
TWOdevice *pDevice;
BOOLEAN tranAnalysis;
TWOtranInfo *info;
{
  TWOelem *pElem;
  TWOnode *pNode;
  TWOedge *pHEdge, *pVEdge;
  TWOedge *pTEdge, *pBEdge, *pLEdge, *pREdge;
  TWOchannel *pCh;
  int index, eIndex;
  double *pRhs = pDevice->rhs;
  double dx, dy, dxdy, dyOverDx, dxOverDy;
  double dPsiT, dPsiB, dPsiL, dPsiR;
  double rhsP;
  double nConc, pConc;
  double perTime;
  void TWOPcommonTerms();
  
  /* first compute the currents */
  TWOPcommonTerms( pDevice, TRUE, tranAnalysis, info );
  
  /* find reciprocal timestep */
  if ( tranAnalysis ) {
    perTime = info->intCoeff[0];
  }
  
  /* zero the rhs vector */
  for ( index = 1 ; index <= pDevice->numEqns ; index++ ) {
    pRhs[ index ] = 0.0;
  }
  
  for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) {
    pElem = pDevice->elements[ eIndex ];

    dx = 0.5 * pElem->dx;
    dy = 0.5 * pElem->dy;
    dxdy = dx * dy;
    dxOverDy = 0.5 * pElem->epsRel * pElem->dxOverDy;
    dyOverDx = 0.5 * pElem->epsRel * pElem->dyOverDx;
    
    pTEdge = pElem->pTopEdge;
    pBEdge = pElem->pBotEdge;
    pLEdge = pElem->pLeftEdge;
    pREdge = pElem->pRightEdge;
    dPsiT = pTEdge->dPsi;
    dPsiB = pBEdge->dPsi;
    dPsiL = pLEdge->dPsi;
    dPsiR = pREdge->dPsi;
    
    /* load for all i,j */
    for ( index = 0; index <= 3; index++ ) {
      pNode = pElem->pNodes[ index ];
      if ( pNode->nodeType ISNOT CONTACT ) {
	if ( index <= 1 ) {
	  pHEdge = pTEdge;
	} else {
	  pHEdge = pBEdge;
	}
	if ( index IS 0 OR index IS 3 ) {
	  pVEdge = pLEdge;
	} else {
	  pVEdge = pREdge;
	}
	/* Add surface state charges. */
	pRhs[ pNode->psiEqn ] += dx * pHEdge->qf;
	pRhs[ pNode->psiEqn ] += dy * pVEdge->qf;
	if ( pElem->elemType IS SEMICON ) {
	  nConc = *(pDevice->devState0 + pNode->nodeN);
	  pConc = *(pDevice->devState0 + pNode->nodeP);
	  pRhs[ pNode->psiEqn ] += dxdy * (pNode->netConc + pConc - nConc);
	  
	  /* Handle generation terms */
	  rhsP =   dxdy * pNode->uNet;
	  pRhs[ pNode->pEqn ] -= rhsP;
	  
	  /* Handle dXdT continuity terms */
	  if ( tranAnalysis ) {
	    pRhs[ pNode->pEqn ] -= dxdy * pNode->dPdT;
	  }
	}
      }
    }
    
    /* Handle neighbor and edge dependent terms */
    pNode = pElem->pTLNode;
    if ( pNode->nodeType ISNOT CONTACT ) {
      pRhs[ pNode->psiEqn ] -= -dyOverDx * dPsiT - dxOverDy * dPsiL;
      if ( pElem->elemType IS SEMICON ) {
	pRhs[ pNode->pEqn ] -= dy * pTEdge->jp + dx * pLEdge->jp;
      }
    }
    pNode = pElem->pTRNode;
    if ( pNode->nodeType ISNOT CONTACT ) {
      pRhs[ pNode->psiEqn ] -= dyOverDx * dPsiT - dxOverDy * dPsiR;
      if ( pElem->elemType IS SEMICON ) {
	pRhs[ pNode->pEqn ] -= -dy * pTEdge->jp + dx * pREdge->jp;
      }
    }
    pNode = pElem->pBRNode;
    if ( pNode->nodeType ISNOT CONTACT ) {
      pRhs[ pNode->psiEqn ] -= dyOverDx * dPsiB + dxOverDy * dPsiR;
      if ( pElem->elemType IS SEMICON ) {
	pRhs[ pNode->pEqn ] -= -dy * pBEdge->jp - dx * pREdge->jp;
      }
    }
    pNode = pElem->pBLNode;
    if ( pNode->nodeType ISNOT CONTACT ) {
      pRhs[ pNode->psiEqn ] -= -dyOverDx * dPsiB + dxOverDy * dPsiL;
      if ( pElem->elemType IS SEMICON ) {
	pRhs[ pNode->pEqn ] -= dy * pBEdge->jp - dx * pLEdge->jp;
      }
    }
  }
}

/*
 * computation of current densities, recombination rates,
 *  mobilities and their derivatives
 */
void 
  TWOPcommonTerms( pDevice, currentOnly, tranAnalysis, info )
TWOdevice *pDevice;
BOOLEAN currentOnly, tranAnalysis;
TWOtranInfo *info;
{
  TWOelem *pElem, *pElem1;
  TWOedge *pEdge;
  TWOnode *pNode;
  int index, eIndex;
  int nextIndex;			/* index of node to find next element */
  double psi1, psi2, refPsi, pC, pP1;
  double dPsiP;
  double bPsiP, dbPsiP, bMPsiP, dbMPsiP;
  double muP, dMuP, rDx, rDy;
  double psi, nConc, pConc;
  double cnAug, cpAug;
  double eSurf;			/* For channel mobilities */
  double qInt;
  TWOchannel *pCh;
  
  /* evaluate all node (including recombination) and edge quantities */
  for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) {
    pElem = pDevice->elements[ eIndex ];
    refPsi = pElem->matlInfo->refPsi;
    cnAug = pElem->matlInfo->cAug[ELEC];
    cpAug = pElem->matlInfo->cAug[HOLE];
    for ( index = 0; index <= 3; index++ ) {
      if ( pElem->evalNodes[ index ] ) { 
	pNode = pElem->pNodes[ index ]; 
	if ( pNode->nodeType ISNOT CONTACT ) { 
	  psi = pDevice->dcSolution[ pNode->psiEqn ];
	  if ( pElem->elemType IS SEMICON ) { 
	    nConc = pNode->nie * exp( psi - refPsi );
	    pConc = pDevice->dcSolution[ pNode->pEqn ];
	    if ( Srh ) {
	      recomb(nConc, pConc, 
		     pNode->tn, pNode->tp, cnAug, cpAug, pNode->nie, 
		     &pNode->uNet, &pNode->dUdN, &pNode->dUdP);
	    } else {
	      pNode->uNet = 0.0;
	      pNode->dUdN = 0.0;
	      pNode->dUdP = 0.0;
	    }
	  }
	} else {
	  /* a contact node */
	  psi = pNode->psi;
	  if ( pElem->elemType IS SEMICON ) { 
	    nConc = pNode->nConc;
	    pConc = pNode->pConc;
	  }
	}
	
	/* store info in the state tables */
	*(pDevice->devState0 + pNode->nodePsi) = psi;
	if ( pElem->elemType IS SEMICON ) {
	  *(pDevice->devState0 + pNode->nodeN) = nConc;
	  *(pDevice->devState0 + pNode->nodeP) = pConc;
	  if ( tranAnalysis AND pNode->nodeType ISNOT CONTACT ) {
	    pNode->dNdT = integrate( pDevice->devStates, info, pNode->nodeN );
	    pNode->dPdT = integrate( pDevice->devStates, info, pNode->nodeP );
	  }
	}
      }
    }
    for ( index = 0; index <= 3; index++ ) {
      if ( pElem->evalEdges[ index ] ) { 
	pEdge = pElem->pEdges[ index ];
	pNode = pElem->pNodes[ index ];
	if ( pNode->nodeType ISNOT CONTACT ) {
	  psi1 = pDevice->dcSolution[pNode->psiEqn]; 
	} else {
	  psi1 = pNode->psi;
	}
	pNode = pElem->pNodes[ (index + 1) % 4 ];
	if ( pNode->nodeType ISNOT CONTACT ) {
	  psi2 = pDevice->dcSolution[pNode->psiEqn]; 
	} else {
	  psi2 = pNode->psi;
	}
	if ( index <= 1 ) {
	  pEdge->dPsi = psi2 - psi1;
	} else {
	  pEdge->dPsi = psi1 - psi2;
	}
	*(pDevice->devState0 + pEdge->edgeDpsi) = pEdge->dPsi;
	
	if ( pElem->elemType IS SEMICON ) {
	  /* Calculate weighted driving forces - wdfn & wdfp for the edge */
	  dPsiP = pEdge->dPsi - pEdge->dVBand;
	  bernoulli( dPsiP, &bPsiP, &dbPsiP,
		    &bMPsiP, &dbMPsiP, NOT currentOnly);
	  if ( index <= 1 ) {
	    pC = *(pDevice->devState0 + pElem->pNodes[ index ]->nodeP);
	    pP1 = *(pDevice->devState0 + pElem->pNodes[ index+1 ]->nodeP);
	  } else {
	    pC = *(pDevice->devState0 + pElem->pNodes[(index+1)%4]->nodeP);
	    pP1 = *(pDevice->devState0 + pElem->pNodes[ index ]->nodeP);
	  }
	  pEdge->wdfp = bPsiP * pC - bMPsiP * pP1;
	  pEdge->jp = 0.0;
	  if ( NOT currentOnly ) {
	    pEdge->dWpDpsiP1 = dbPsiP * pC - dbMPsiP * pP1;
	    pEdge->dWpDp     = bPsiP;
	    pEdge->dWpDpP1   = - bMPsiP;
	    pEdge->dJpDpsiP1 = 0.0;
	    pEdge->dJpDp     = 0.0;
	    pEdge->dJpDpP1   = 0.0;
	  }
	}
      }
    }
  }
  
  /* DAG: calculate mobilities for channel elems */
  if ( SurfaceMobility ) {
    for ( pCh = pDevice->pChannel;
	 pCh ISNOT NIL(TWOchannel); pCh = pCh->next ) {
      pElem = pCh->pNElem;
      switch (pCh->type) {
      case 0:
	eSurf = - 0.5 * (pElem->pLeftEdge->dPsi + pElem->pRightEdge->dPsi )
	  * pElem->epsRel / pElem->dy;
	qInt = 0.5 * pElem->pBotEdge->qf;
	break;
      case 1:
	eSurf = - 0.5 * (pElem->pTopEdge->dPsi + pElem->pBotEdge->dPsi )
	  * pElem->epsRel / pElem->dx;
	qInt = 0.5 * pElem->pLeftEdge->qf;
	break;
      case 2:
	eSurf = - 0.5 * (pElem->pLeftEdge->dPsi + pElem->pRightEdge->dPsi )
	  * pElem->epsRel / pElem->dy;
	qInt = 0.5 * pElem->pTopEdge->qf;
	break;
      case 3:
	eSurf = - 0.5 * (pElem->pTopEdge->dPsi + pElem->pBotEdge->dPsi )
	  * pElem->epsRel / pElem->dx;
	qInt = 0.5 * pElem->pRightEdge->qf;
	break;
      }
      eSurf += qInt;
      pElem = pCh->pSeed;
      nextIndex = (pCh->type + 2)%4;
      while (pElem && pElem->channel == pCh->id) {
	TWOPmobility( pElem, eSurf );
	pElem = pElem->pElems[ nextIndex ];
      }
    } /* endfor pCH ISNOT NIL */
  } /* endif SurfaceMobility */
  
  /* calculate the current densities assuming mobility value depend on Ex,Ey*/
  for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) {
    pElem = pDevice->elements[ eIndex ];
    
    rDx = 1.0 / pElem->dx;
    rDy = 1.0 / pElem->dy;
    for ( index = 0; index <= 3; index++ ) {
      pEdge = pElem->pEdges[ index ];
      /* calculate conductive currents */
      if ( pElem->elemType IS SEMICON ) {
	/* get mobility for this edge */
	if ( NOT pElem->channel ) {
	  /* Calculate mobility for non-channel elements */
	  muP = pElem->mup0;
	  dMuP = 0.0;
	  dPsiP = pEdge->dPsi - pEdge->dVBand;
	  if ( index%2 IS 0 ) {
	    MOBfieldDep( pElem->matlInfo, HOLE, - dPsiP * rDx, &muP, &dMuP );
	  } else {
	    MOBfieldDep( pElem->matlInfo, HOLE, - dPsiP * rDy, &muP, &dMuP );
	  }
	} else {
	  /* Retrieve previously calculated value. */
	  muP = pElem->mup;
	  dMuP = 0.0;
	}
	switch ( index ) {
	case 0:
	  muP *= pEdge->kPos * rDx;
	  dMuP *= pEdge->kPos * rDx * rDx;
	  break;
	case 1:
	  muP *= pEdge->kNeg * rDy;
	  dMuP *= pEdge->kNeg * rDy * rDy;
	  break;
	case 2:
	  muP *= pEdge->kNeg * rDx;
	  dMuP *= pEdge->kNeg * rDx * rDx;
	  break;
	case 3:
	  muP *= pEdge->kPos * rDy;
	  dMuP *= pEdge->kPos * rDy * rDy;
	  break;
	}
	/* Now that the mobility for this edge is known, do current */
	pEdge->jp += muP * pEdge->wdfp;
	if ( NOT currentOnly ) {
	  pEdge->dJpDpsiP1 += muP * pEdge->dWpDpsiP1;
	  pEdge->dJpDp += muP * pEdge->dWpDp;
	  pEdge->dJpDpP1 += muP * pEdge->dWpDpP1;
	  if ( MobDeriv AND (NOT pElem->channel) ) {
	    pEdge->dJpDpsiP1 -= dMuP * pEdge->wdfp;
	  }
	}
      }
      /* calculate displacement current only once */
      if ( pElem->evalEdges[ index ] ) { 
	if ( tranAnalysis ) {
	  if ( index IS 0 OR index IS 2 ) {
	    /* horizontal edges */
	    pEdge->jd = -integrate(pDevice->devStates, info,
				   pEdge->edgeDpsi) * rDx;
	  } else {
	    /* vertical edges */
	    pEdge->jd = -integrate(pDevice->devStates, info,
				   pEdge->edgeDpsi) * rDy;
	  }
	}
      }
    }
  }
}