oneadmit.c

ngspice又一个电子CAD仿真软件代码.功能更全

/**********
Copyright 1992 Regents of the University of California.  All rights reserved.
Author:	1987 Kartikeya Mayaram, U. C. Berkeley CAD Group
$Id: oneadmit.c,v 1.3 2005/05/23 06:58:42 sjborley Exp $
**********/

/* Functions to compute small-signal parameters of 1D devices */

#include "ngspice.h"
#include "numglobs.h"
#include "numenum.h"
#include "numconst.h"
#include "onedev.h"
#include "onemesh.h"
#include "complex.h"
#include "spmatrix.h"
#include "ifsim.h"

#include "onedext.h"
#include "oneddefs.h"
#include "cidersupt.h"


extern IFfrontEnd *SPfrontEnd;


/* 
 * mmhhh this may cause troubles
 * Paolo Nenzi 2002
 */
 SPcomplex yAc;


/* Forward Declarations */
SPcomplex *computeAdmittance(ONEnode *, BOOLEAN, double *, 
                  double *, SPcomplex *);
BOOLEAN ONEsorSolve(ONEdevice *pDevice, double *xReal, 
                    double *xImag, double omega);

int
NUMDadmittance(ONEdevice *pDevice, double omega, SPcomplex *yd)
{
  ONEnode *pNode;
  ONEelem *pElem;
  ONEedge *pEdge;
  int index, i;
  double yReal, yImag;
  double *solutionReal, *solutionImag;
  SPcomplex yAc, cOmega;
  SPcomplex *y;
  BOOLEAN SORFailed;
  double startTime;


  /* Each time we call this counts as one AC iteration. */
  pDevice->pStats->numIters[STAT_AC] += 1;

  /*
   * Change context names of solution vectors for ac analysis.
   * dcDeltaSolution stores the real part and copiedSolution stores the
   * imaginary part of the ac solution vector.
   */
  solutionReal = pDevice->dcDeltaSolution;
  solutionImag = pDevice->copiedSolution;
  pDevice->solverType = SLV_SMSIG;

  /* use a normalized radian frequency */
  omega *= TNorm;
  CMPLX_ASSIGN_VALUE(cOmega, 0.0, omega);
  yReal = 0.0;
  yImag = 0.0;

  if ((AcAnalysisMethod == SOR) || (AcAnalysisMethod == SOR_ONLY)) {
    /* LOAD */
    startTime = SPfrontEnd->IFseconds();
    /* zero the rhs before loading in the new rhs */
    for (index = 1; index <= pDevice->numEqns; index++) {
      pDevice->rhs[index] = 0.0;
      pDevice->rhsImag[index] = 0.0;
    }
    /* store the new rhs vector */
    pElem = pDevice->elemArray[pDevice->numNodes - 1];
    pNode = pElem->pLeftNode;
    pDevice->rhs[pNode->psiEqn] = pElem->epsRel * pElem->rDx;
    if (pElem->elemType == SEMICON) {
      pEdge = pElem->pEdge;
      pDevice->rhs[pNode->nEqn] -= pEdge->dJnDpsiP1;
      pDevice->rhs[pNode->pEqn] -= pEdge->dJpDpsiP1;
    }
    pDevice->pStats->loadTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;

    /* SOLVE */
    startTime = SPfrontEnd->IFseconds();
    SORFailed = ONEsorSolve(pDevice, solutionReal, solutionImag, omega);
    pDevice->pStats->solveTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;

    if (SORFailed && AcAnalysisMethod == SOR) {
      AcAnalysisMethod = DIRECT;
      printf("SOR failed at %g Hz, switching to direct-method ac analysis.\n",
	  omega / (TWO_PI * TNorm) );
    } else if (SORFailed) {	/* Told to only do SOR, so give up. */
      printf("SOR failed at %g Hz, returning null admittance.\n",
	  omega / (TWO_PI * TNorm) );
      CMPLX_ASSIGN_VALUE(*yd, 0.0, 0.0);
      return (AcAnalysisMethod);
    }
  }
  if (AcAnalysisMethod == DIRECT) {
    /* LOAD */
    startTime = SPfrontEnd->IFseconds();
    /* solve the system of equations directly */
    for (index = 1; index <= pDevice->numEqns; index++) {
      pDevice->rhs[index] = 0.0;
      pDevice->rhsImag[index] = 0.0;
    }
    pElem = pDevice->elemArray[pDevice->numNodes - 1];
    pNode = pElem->pLeftNode;
    pDevice->rhs[pNode->psiEqn] = pElem->epsRel * pElem->rDx;
    if (pElem->elemType == SEMICON) {
      pEdge = pElem->pEdge;
      pDevice->rhs[pNode->nEqn] -= pEdge->dJnDpsiP1;
      pDevice->rhs[pNode->pEqn] -= pEdge->dJpDpsiP1;
    }
    ONE_jacLoad(pDevice);
    spSetComplex(pDevice->matrix);
    for (index = 1; index < pDevice->numNodes; index++) {
      pElem = pDevice->elemArray[index];
      if (pElem->elemType == SEMICON) {
	for (i = 0; i <= 1; i++) {
	  pNode = pElem->pNodes[i];
	  if (pNode->nodeType != CONTACT) {
	    spADD_COMPLEX_ELEMENT(pNode->fNN, 0.0, -0.5 * pElem->dx * omega);
	    spADD_COMPLEX_ELEMENT(pNode->fPP, 0.0, 0.5 * pElem->dx * omega);
	  }
	}
      }
    }
    pDevice->pStats->loadTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;

    /* FACTOR */
    startTime = SPfrontEnd->IFseconds();
    spFactor(pDevice->matrix);
    pDevice->pStats->factorTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;

    /* SOLVE */
    startTime = SPfrontEnd->IFseconds();
    spSolve(pDevice->matrix, pDevice->rhs, solutionReal,
	pDevice->rhsImag, solutionImag);
    pDevice->pStats->solveTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
  }
  /* MISC */
  startTime = SPfrontEnd->IFseconds();
  pNode = pDevice->elemArray[1]->pLeftNode;
  y = computeAdmittance(pNode, FALSE, solutionReal, solutionImag, &cOmega);
  CMPLX_ASSIGN_VALUE(yAc, -y->real, -y->imag);
  CMPLX_ASSIGN(*yd, yAc);
  CMPLX_MULT_SELF_SCALAR(*yd, GNorm * pDevice->area);
  pDevice->pStats->miscTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;

  return (AcAnalysisMethod);
}


int
NBJTadmittance(ONEdevice *pDevice, double omega, SPcomplex *yIeVce, 
               SPcomplex *yIcVce, SPcomplex *yIeVbe, SPcomplex *yIcVbe)
{
  ONEelem *pCollElem = pDevice->elemArray[pDevice->numNodes - 1];
  ONEelem *pBaseElem = pDevice->elemArray[pDevice->baseIndex - 1];
  ONEelem *pElem;
  ONEedge *pEdge;
  ONEnode *pNode;
  int index, i;
  double area = pDevice->area;
  double *solutionReal, *solutionImag;
  BOOLEAN SORFailed;
  SPcomplex *y;
  SPcomplex cOmega, pIeVce, pIcVce, pIeVbe, pIcVbe;
  double startTime;

  /* Each time we call this counts as one AC iteration. */
  pDevice->pStats->numIters[STAT_AC] += 1;

  /*
   * change context names of solution vectors for ac analysis dcDeltaSolution
   * stores the real part and copiedSolution stores the imaginary part of the
   * ac solution vector
   */
  solutionReal = pDevice->dcDeltaSolution;
  solutionImag = pDevice->copiedSolution;
  pDevice->solverType = SLV_SMSIG;

  /* use a normalized radian frequency */
  omega *= TNorm;
  CMPLX_ASSIGN_VALUE(cOmega, 0.0, omega);

  if ((AcAnalysisMethod == SOR) || (AcAnalysisMethod == SOR_ONLY)) {
    /* LOAD */
    startTime = SPfrontEnd->IFseconds();
    /* zero the rhs before loading in the new rhs */
    for (index = 1; index <= pDevice->numEqns; index++) {
      pDevice->rhs[index] = 0.0;
      pDevice->rhsImag[index] = 0.0;
    }
    /* store the new rhs vector */
    pNode = pCollElem->pLeftNode;
    pDevice->rhs[pNode->psiEqn] = pCollElem->epsRel * pCollElem->rDx;
    if (pCollElem->elemType == SEMICON) {
      pEdge = pCollElem->pEdge;
      pDevice->rhs[pNode->nEqn] -= pEdge->dJnDpsiP1;
      pDevice->rhs[pNode->pEqn] -= pEdge->dJpDpsiP1;
    }
    pDevice->pStats->loadTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;

    /* SOLVE */
    startTime = SPfrontEnd->IFseconds();
    SORFailed = ONEsorSolve(pDevice, solutionReal, solutionImag, omega);
    pDevice->pStats->solveTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
    if (SORFailed && (AcAnalysisMethod == SOR)) {
      AcAnalysisMethod = DIRECT;
      printf("SOR failed at %g Hz, switching to direct-method ac analysis.\n",
	  omega / (TWO_PI * TNorm) );
    } else if (SORFailed) {	/* Told to only do SOR, so give up. */
      printf("SOR failed at %g Hz, returning null admittance.\n",
	  omega / (TWO_PI * TNorm) );
      CMPLX_ASSIGN_VALUE(*yIeVce, 0.0, 0.0);
      CMPLX_ASSIGN_VALUE(*yIcVce, 0.0, 0.0);
      CMPLX_ASSIGN_VALUE(*yIeVbe, 0.0, 0.0);
      CMPLX_ASSIGN_VALUE(*yIcVbe, 0.0, 0.0);
      return (AcAnalysisMethod);
    } else {
      /* MISC */
      startTime = SPfrontEnd->IFseconds();
      pElem = pDevice->elemArray[1];
      pNode = pElem->pLeftNode;
      y = computeAdmittance(pNode, FALSE, solutionReal, solutionImag, &cOmega);
      CMPLX_ASSIGN_VALUE(pIeVce, -y->real, -y->imag);
      pNode = pCollElem->pRightNode;
      y = computeAdmittance(pNode, TRUE, solutionReal, solutionImag, &cOmega);
      CMPLX_ASSIGN_VALUE(pIcVce, -y->real, -y->imag);
      pDevice->pStats->miscTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;

       /* LOAD */
      startTime = SPfrontEnd->IFseconds();
      /* load in the base contribution to the rhs */
      for (index = 1; index <= pDevice->numEqns; index++) {
	pDevice->rhs[index] = 0.0;
      }
      pNode = pBaseElem->pRightNode;
      if (pNode->baseType == N_TYPE) {
	pDevice->rhs[pNode->nEqn] = pNode->nConc * pNode->eg;
      } else if (pNode->baseType == P_TYPE) {
	pDevice->rhs[pNode->pEqn] = pNode->pConc * pNode->eg;
      } else {
	printf("projectBJTsolution: unknown base type\n");
      }
      pDevice->pStats->loadTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;

      /* SOLVE */
      startTime = SPfrontEnd->IFseconds();
      SORFailed = ONEsorSolve(pDevice, solutionReal, solutionImag, omega);
      pDevice->pStats->solveTime[STAT_AC] +=
	  SPfrontEnd->IFseconds() - startTime;
      if (SORFailed && (AcAnalysisMethod == SOR)) {
	AcAnalysisMethod = DIRECT;
	printf("SOR failed at %g Hz, switching to direct-method ac analysis.\n",
	    omega / (TWO_PI * TNorm) );
      } else if (SORFailed) {	/* Told to only do SOR, so give up. */
	printf("SOR failed at %g Hz, returning null admittance.\n",
	    omega / (TWO_PI * TNorm) );
	CMPLX_ASSIGN_VALUE(*yIeVce, 0.0, 0.0);
	CMPLX_ASSIGN_VALUE(*yIcVce, 0.0, 0.0);
	CMPLX_ASSIGN_VALUE(*yIeVbe, 0.0, 0.0);
	CMPLX_ASSIGN_VALUE(*yIcVbe, 0.0, 0.0);
	return (AcAnalysisMethod);
      }
    }
  }
  if (AcAnalysisMethod == DIRECT) {
    /* LOAD */
    startTime = SPfrontEnd->IFseconds();
    for (index = 1; index <= pDevice->numEqns; index++) {
      pDevice->rhs[index] = 0.0;
      pDevice->rhsImag[index] = 0.0;
    }
    /* solve the system of equations directly */
    ONE_jacLoad(pDevice);
    pNode = pCollElem->pLeftNode;
    pDevice->rhs[pNode->psiEqn] = pCollElem->epsRel * pCollElem->rDx;
    if (pCollElem->elemType == SEMICON) {
      pEdge = pCollElem->pEdge;
      pDevice->rhs[pNode->nEqn] -= pEdge->dJnDpsiP1;
      pDevice->rhs[pNode->pEqn] -= pEdge->dJpDpsiP1;
    }
    spSetComplex(pDevice->matrix);
    for (index = 1; index < pDevice->numNodes; index++) {
      pElem = pDevice->elemArray[index];
      if (pElem->elemType == SEMICON) {
	for (i = 0; i <= 1; i++) {
	  pNode = pElem->pNodes[i];
	  if (pNode->nodeType != CONTACT) {
	    spADD_COMPLEX_ELEMENT(pNode->fNN, 0.0, -0.5 * pElem->dx * omega);
	    spADD_COMPLEX_ELEMENT(pNode->fPP, 0.0, 0.5 * pElem->dx * omega);
	  }
	}
      }
    }
    pDevice->pStats->loadTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;

    /* FACTOR */
    startTime = SPfrontEnd->IFseconds();
    spFactor(pDevice->matrix);
    pDevice->pStats->factorTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;

    /* SOLVE */
    startTime = SPfrontEnd->IFseconds();
    spSolve(pDevice->matrix, pDevice->rhs, solutionReal,
	pDevice->rhsImag, solutionImag);
    pDevice->pStats->solveTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;

    /* MISC */
    startTime = SPfrontEnd->IFseconds();
    pElem = pDevice->elemArray[1];
    pNode = pElem->pLeftNode;
    y = computeAdmittance(pNode, FALSE, solutionReal, solutionImag, &cOmega);
    CMPLX_ASSIGN_VALUE(pIeVce, -y->real, -y->imag);
    pNode = pCollElem->pRightNode;
    y = computeAdmittance(pNode, TRUE, solutionReal, solutionImag, &cOmega);
    CMPLX_ASSIGN_VALUE(pIcVce, -y->real, -y->imag);
    pDevice->pStats->miscTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;

    /* LOAD */
    startTime = SPfrontEnd->IFseconds();
    /* load in the base contribution in the rhs */
    for (index = 1; index <= pDevice->numEqns; index++) {
      pDevice->rhs[index] = 0.0;
    }
    pNode = pBaseElem->pRightNode;
    if (pNode->baseType == N_TYPE) {
      pDevice->rhs[pNode->nEqn] = pNode->nConc * pNode->eg;
    } else if (pNode->baseType == P_TYPE) {
      pDevice->rhs[pNode->pEqn] = pNode->pConc * pNode->eg;
    } else {
      printf("\n BJTadmittance: unknown base type");
    }
    pDevice->pStats->loadTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;

    /* FACTOR: already done, no need to repeat. */

    /* SOLVE */
    startTime = SPfrontEnd->IFseconds();
    spSolve(pDevice->matrix, pDevice->rhs, solutionReal,
	pDevice->rhsImag, solutionImag);
    pDevice->pStats->solveTime[STAT_AC] += SPfrontEnd->IFseconds() - startTime;
  }
  /* MISC */
  startTime = SPfrontEnd->IFseconds();
  pElem = pDevice->elemArray[1];
  pNode = pElem->pLeftNode;
  y = computeAdmittance(pNode, FALSE, solutionReal, solutionImag, &cOmega);
  CMPLX_ASSIGN_VALUE(pIeVbe, -y->real, -y->imag);
  pNode = pCollElem->pRightNode;
  y = computeAdmittance(pNode, FALSE, solutionReal, solutionImag, &cOmega);
  CMPLX_ASSIGN_VALUE(pIcVbe, -y->real, -y->imag);

  CMPLX_ASSIGN(*yIeVce, pIeVce);
  CMPLX_ASSIGN(*yIcVce, pIcVce);
  CMPLX_ASSIGN(*yIeVbe, pIeVbe);
  CMPLX_ASSIGN(*yIcVbe, pIcVbe);
  CMPLX_MULT_SELF_SCALAR(*yIeVce, GNorm * area);
  CMPLX_MULT_SELF_SCALAR(*yIeVbe, GNorm * area);
  CMPLX_MULT_SELF_SCALAR(*yIcVce, GNorm * area);