ltramisc.c

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

  double time, cbyr, rclsqr;
{
  double temp;

  if (time != 0.0) {
    temp = rclsqr / (4 * time);
    temp = 2 * sqrt(time / M_PI) * exp(-temp) - sqrt(rclsqr) * erfc(sqrt(temp));
    return (sqrt(cbyr) * temp);
  } else {
    return (0.0);
  }
}


/*
 * LTRArcCoeffsSetup sets up the all coefficient lists for the special case
 * where L=G=0
 */

void
LTRArcCoeffsSetup(h1dashfirstcoeff, h2firstcoeff, h3dashfirstcoeff, h1dashcoeffs, h2coeffs, h3dashcoeffs, listsize, cbyr, rclsqr, curtime, timelist, timeindex, reltol)
  double *h1dashcoeffs, *h2coeffs, *h3dashcoeffs, *timelist;
  int listsize, timeindex;
  double cbyr, rclsqr, curtime, *h1dashfirstcoeff, *h2firstcoeff, *h3dashfirstcoeff;
  double reltol;

{
  double delta1, delta2;
  double h1dummy1, h1dummy2;
  double h2dummy1, h2dummy2;
  double h3dummy1, h3dummy2;
  double lolimit1, lolimit2, hilimit1, hilimit2;
  double h1lovalue1, h1lovalue2, h1hivalue1, h1hivalue2;
  double h2lovalue1, h2lovalue2, h2hivalue1, h2hivalue2;
  double h3lovalue1, h3lovalue2, h3hivalue1, h3hivalue2;
  double temp, temp2, temp3, temp4, temp5;
  double h1relval, h2relval, h3relval;
  int doh1 = 1, doh2 = 1, doh3 = 1;
  int i, auxindex;

  /* coefflists should already have been allocated to the necessary size */

#ifdef LTRAdebug
  if (listsize <= timeindex) {
    printf("LTRAcoeffSetup: not enough space in coefflist\n");
  }
#endif



  auxindex = timeindex;

  /* the first coefficients */

  delta1 = curtime - *(timelist + auxindex);
  lolimit1 = 0.0;
  hilimit1 = delta1;

  h1lovalue1 = 0.0;
  h1hivalue1 =			/* LTRArcH1dashTwiceIntFunc(hilimit1,cbyr); */
      sqrt(4 * cbyr * hilimit1 / M_PI);
  h1dummy1 = h1hivalue1 / delta1;
  *h1dashfirstcoeff = h1dummy1;
  h1relval = FABS(h1dummy1 * reltol);

  temp = rclsqr / (4 * hilimit1);
  temp2 = (temp >= 100.0 ? 0.0 : erfc(sqrt(temp)));
  temp3 = exp(-temp);
  temp4 = sqrt(rclsqr);
  temp5 = sqrt(cbyr);

  h2lovalue1 = 0.0;
  h2hivalue1 =			/* LTRArcH2TwiceIntFunc(hilimit1,rclsqr); */
      (hilimit1 != 0.0 ? (hilimit1 + rclsqr * 0.5) * temp2 - sqrt(hilimit1 * rclsqr / M_PI) * temp3 : 0.0);


  h2dummy1 = h2hivalue1 / delta1;
  *h2firstcoeff = h2dummy1;
  h2relval = FABS(h2dummy1 * reltol);

  h3lovalue1 = 0.0;
  h3hivalue1 =			/* LTRArcH3dashTwiceIntFunc(hilimit1,cbyr,rcls
				 * qr); */
      (hilimit1 != 0.0 ? temp = 2 * sqrt(hilimit1 / M_PI) * temp3 - temp4 * temp2, (temp5 * temp) : 0.0);


  h3dummy1 = h3hivalue1 / delta1;
  *h3dashfirstcoeff = h3dummy1;
  h3relval = FABS(h3dummy1 * reltol);

  /* the coefficients for the rest of the timepoints */

  for (i = auxindex; i > 0; i--) {

    delta2 = delta1;		/* previous delta1 */
    lolimit2 = lolimit1;	/* previous lolimit1 */
    hilimit2 = hilimit1;	/* previous hilimit1 */

    delta1 = *(timelist + i) - *(timelist + i - 1);
    lolimit1 = hilimit2;
    hilimit1 = curtime - *(timelist + i - 1);

    if (doh1) {
      h1lovalue2 = h1lovalue1;	/* previous lovalue1 */
      h1hivalue2 = h1hivalue1;	/* previous hivalue1 */
      h1dummy2 = h1dummy1;	/* previous dummy1 */

      h1lovalue1 = h1hivalue2;
      h1hivalue1 =		/* LTRArcH1dashTwiceIntFunc(hilimit1,cbyr); */
	  sqrt(4 * cbyr * hilimit1 / M_PI);
      h1dummy1 = (h1hivalue1 - h1lovalue1) / delta1;
      *(h1dashcoeffs + i) = h1dummy1 - h1dummy2;
      if (FABS(*(h1dashcoeffs + i)) < h1relval)
	doh1 = 0;
    } else
      *(h1dashcoeffs + i) = 0.0;

    if (doh2 || doh3) {
      temp = rclsqr / (4 * hilimit1);
      temp2 = (temp >= 100.0 ? 0.0 : erfc(sqrt(temp)));
      temp3 = exp(-temp);
    }
    if (doh2) {
      h2lovalue2 = h2lovalue1;	/* previous lovalue1 */
      h2hivalue2 = h2hivalue1;	/* previous hivalue1 */
      h2dummy2 = h2dummy1;	/* previous dummy1 */

      h2lovalue1 = h2hivalue2;
      h2hivalue1 =		/* LTRArcH2TwiceIntFunc(hilimit1,rclsqr); */
	  (hilimit1 != 0.0 ? (hilimit1 + rclsqr * 0.5) * temp2 - sqrt(hilimit1 * rclsqr / M_PI) * temp3 : 0.0);
      h2dummy1 = (h2hivalue1 - h2lovalue1) / delta1;
      *(h2coeffs + i) = h2dummy1 - h2dummy2;
      if (FABS(*(h2coeffs + i)) < h2relval)
	doh2 = 0;
    } else
      *(h2coeffs + i) = 0.0;

    if (doh3) {
      h3lovalue2 = h3lovalue1;	/* previous lovalue1 */
      h3hivalue2 = h3hivalue1;	/* previous hivalue1 */
      h3dummy2 = h3dummy1;	/* previous dummy1 */

      h3lovalue1 = h3hivalue2;
      h3hivalue1 =		/* LTRArcH3dashTwiceIntFunc(hilimit1,cbyr,rcls
				 * qr); */
	  (hilimit1 != 0.0 ? temp = 2 * sqrt(hilimit1 / M_PI) * temp3 - temp4 * temp2, (temp5 * temp) : 0.0);
      h3dummy1 = (h3hivalue1 - h3lovalue1) / delta1;
      *(h3dashcoeffs + i) = h3dummy1 - h3dummy2;
      if (FABS(*(h3dashcoeffs + i)) < h3relval)
	doh3 = 0;
    } else
      *(h3dashcoeffs + i) = 0.0;
  }
}

void
LTRArlcCoeffsSetup(h1dashfirstcoeff, h2firstcoeff, h3dashfirstcoeff, h1dashcoeffs, h2coeffs, h3dashcoeffs, listsize, T, alpha, beta, curtime, timelist, timeindex, reltol, auxindexptr)
  double *h1dashfirstcoeff, *h2firstcoeff, *h3dashfirstcoeff;
  double *h1dashcoeffs, *h2coeffs, *h3dashcoeffs, *timelist;
  double T, alpha, beta, curtime, reltol;
  int listsize, timeindex, *auxindexptr;
{
  unsigned exact;
  double lolimit1, lolimit2, hilimit1, hilimit2;
  double delta1, delta2;

  double h1dummy1, h1dummy2;
  double h1lovalue1, h1lovalue2, h1hivalue1, h1hivalue2;

  double h2dummy1, h2dummy2;
  double h2lovalue1, h2lovalue2, h2hivalue1, h2hivalue2;

  double h3dummy1, h3dummy2;
  double h3lovalue1, h3lovalue2, h3hivalue1, h3hivalue2;

  double exparg, besselarg, expterm, bessi1overxterm, bessi0term;
  double expbetaTterm, alphasqTterm;
  double h1relval, h2relval, h3relval;
  int doh1 = 1, doh2 = 1, doh3 = 1;

  int i, auxindex;

  /* coefflists should already have been allocated to the necessary size */

#ifdef LTRAdebug
  if (listsize <= timeindex) {
    printf("LTRArlcCoeffsSetup: not enough space in coefflist\n");
  }
#endif


  /*
   * we assume a piecewise linear function, and we calculate the coefficients
   * using this assumption in the integration of the function
   */

  if (T == 0.0) {
    auxindex = timeindex;
  } else {

    if (curtime - T <= 0.0) {
      auxindex = 0;
    } else {
      exact = 0;
      for (i = timeindex; i >= 0; i--) {
	if (curtime - *(timelist + i) == T) {
	  exact = 1;
	  break;
	}
	if (curtime - *(timelist + i) > T)
	  break;
      }

#ifdef LTRADEBUG
      if ((i < 0) || ((i == 0) && (exact == 1)))
	printf("LTRAcoeffSetup: i <= 0: some mistake!\n");
#endif

      if (exact == 1) {
	auxindex = i - 1;
      } else {
	auxindex = i;
      }
    }
  }
  /* the first coefficient */

  if (auxindex != 0) {
    lolimit1 = T;
    hilimit1 = curtime - *(timelist + auxindex);
    delta1 = hilimit1 - lolimit1;

    h2lovalue1 = LTRArlcH2Func(T, T, alpha, beta);
    besselarg = (hilimit1 > T) ? alpha * sqrt(hilimit1 * hilimit1 - T * T) : 0.0;
    exparg = -beta * hilimit1;
    expterm = exp(exparg);
    bessi1overxterm = bessI1xOverX(besselarg);
    alphasqTterm = alpha * alpha * T;
    h2hivalue1 =		/* LTRArlcH2Func(hilimit1,T,alpha,beta); */
	((alpha == 0.0) || (hilimit1 < T)) ? 0.0 : alphasqTterm * expterm * bessi1overxterm;

    h2dummy1 = twiceintlinfunc(lolimit1, hilimit1, lolimit1, h2lovalue1,
	h2hivalue1, lolimit1, hilimit1) / delta1;
    *h2firstcoeff = h2dummy1;
    h2relval = FABS(reltol * h2dummy1);

    h3lovalue1 = 0.0;		/* E3dash should be consistent with this */
    bessi0term = bessI0(besselarg);
    expbetaTterm = exp(-beta * T);
    h3hivalue1 =		/* LTRArlcH3dashIntFunc(hilimit1,T,beta); */
	((hilimit1 <= T) || (beta == 0.0)) ? 0.0 : expterm * bessi0term - expbetaTterm;
    h3dummy1 = intlinfunc(lolimit1, hilimit1, h3lovalue1,
	h3hivalue1, lolimit1, hilimit1) / delta1;
    *h3dashfirstcoeff = h3dummy1;
    h3relval = FABS(h3dummy1 * reltol);
  } else {
    *h2firstcoeff = *h3dashfirstcoeff = 0.0;
  }

  lolimit1 = 0.0;
  hilimit1 = curtime - *(timelist + timeindex);
  delta1 = hilimit1 - lolimit1;
  exparg = -beta * hilimit1;
  expterm = exp(exparg);

  h1lovalue1 = 0.0;
  h1hivalue1 =			/* LTRArlcH1dashTwiceIntFunc(hilimit1,beta); */
      (beta == 0.0) ? hilimit1 : ((hilimit1 == 0.0) ? 0.0 : (bessI1(-exparg) + bessI0(-exparg)) * hilimit1 * expterm - hilimit1);
  h1dummy1 = h1hivalue1 / delta1;
  *h1dashfirstcoeff = h1dummy1;
  h1relval = FABS(h1dummy1 * reltol);


  /* the coefficients for the rest of the timepoints */

  for (i = timeindex; i > 0; i--) {

    if (doh1 || doh2 || doh3) {
      lolimit2 = lolimit1;	/* previous lolimit1 */
      hilimit2 = hilimit1;	/* previous hilimit1 */
      delta2 = delta1;		/* previous delta1 */

      lolimit1 = hilimit2;
      hilimit1 = curtime - *(timelist + i - 1);
      delta1 = *(timelist + i) - *(timelist + i - 1);

      exparg = -beta * hilimit1;
      expterm = exp(exparg);
    }
    if (doh1) {
      h1lovalue2 = h1lovalue1;	/* previous lovalue1 */
      h1hivalue2 = h1hivalue1;	/* previous hivalue1 */
      h1dummy2 = h1dummy1;	/* previous dummy1 */

      h1lovalue1 = h1hivalue2;
      h1hivalue1 =		/* LTRArlcH1dashTwiceIntFunc(hilimit1,beta); */
	  (beta == 0.0) ? hilimit1 : ((hilimit1 == 0.0) ? 0.0 : (bessI1(-exparg) + bessI0(-exparg)) * hilimit1 * expterm - hilimit1);
      h1dummy1 = (h1hivalue1 - h1lovalue1) / delta1;

      *(h1dashcoeffs + i) = h1dummy1 - h1dummy2;
      if (FABS(*(h1dashcoeffs + i)) <= h1relval)
	doh1 = 0;
    } else
      *(h1dashcoeffs + i) = 0.0;

    if (i <= auxindex) {

      /*
       * if (i == auxindex) { lolimit2 = T; delta2 = hilimit2 - lolimit2; }
       */

      if (doh2 || doh3)
	besselarg = (hilimit1 > T) ? alpha * sqrt(hilimit1 * hilimit1 - T * T) : 0.0;

      if (doh2) {
	h2lovalue2 = h2lovalue1;/* previous lovalue1 */
	h2hivalue2 = h2hivalue1;/* previous hivalue1 */
	h2dummy2 = h2dummy1;	/* previous dummy1 */

	h2lovalue1 = h2hivalue2;
	bessi1overxterm = bessI1xOverX(besselarg);
	h2hivalue1 =		/* LTRArlcH2Func(hilimit1,T,alpha,beta); */
	    ((alpha == 0.0) || (hilimit1 < T)) ? 0.0 : alphasqTterm * expterm * bessi1overxterm;
	h2dummy1 = twiceintlinfunc(lolimit1, hilimit1, lolimit1,
	    h2lovalue1, h2hivalue1, lolimit1, hilimit1) / delta1;

	*(h2coeffs + i) = h2dummy1 - h2dummy2 + intlinfunc(lolimit2, hilimit2,
	    h2lovalue2, h2hivalue2, lolimit2, hilimit2);
	if (FABS(*(h2coeffs + i)) <= h2relval)
	  doh2 = 0;
      } else
	*(h2coeffs + i) = 0.0;

      if (doh3) {
	h3lovalue2 = h3lovalue1;/* previous lovalue1 */
	h3hivalue2 = h3hivalue1;/* previous hivalue1 */
	h3dummy2 = h3dummy1;	/* previous dummy1 */

	h3lovalue1 = h3hivalue2;
	bessi0term = bessI0(besselarg);
	h3hivalue1 =		/* LTRArlcH3dashIntFunc(hilimit1,T,beta); */
	    ((hilimit1 <= T) || (beta == 0.0)) ? 0.0 : expterm * bessi0term - expbetaTterm;
	h3dummy1 = intlinfunc(lolimit1, hilimit1, h3lovalue1, h3hivalue1, lolimit1, hilimit1) / delta1;

	*(h3dashcoeffs + i) = h3dummy1 - h3dummy2;
	if (FABS(*(h3dashcoeffs + i)) <= h3relval)
	  doh3 = 0;
      } else
	*(h3dashcoeffs + i) = 0.0;
    }
  }
  *auxindexptr = auxindex;
}

/*
 * LTRAstraightLineCheck - takes the co-ordinates of three points, finds the
 * area of the triangle enclosed by these points and compares this area with
 * the area of the quadrilateral formed by the line between the first point
 * and the third point, the perpendiculars from the first and third points to
 * the x-axis, and the x-axis. If within reltol, then it returns 1, else 0.
 * The purpose of this function is to determine if three points lie
 * acceptably close to a straight line. This area criterion is used because
 * it is related to integrals and convolution
 */

int
LTRAstraightLineCheck(x1, y1, x2, y2, x3, y3, reltol, abstol)
  double x1, y1, x2, y2, x3, y3, reltol, abstol;

{
  /*
   * double asqr, bsqr, csqr, c, c1sqr; double htsqr;
   */
  double TRarea, QUADarea1, QUADarea2, QUADarea3, area;

  /*
   * asqr = (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1); bsqr = (x3-x2)*(x3-x2) +
   * (y3-y2)*(y3-y2); csqr = (x3-x1)*(x3-x1) + (y3-y1)*(y3-y1); c =
   * sqrt(csqr); c1sqr = (asqr - bsqr + csqr)/(2*c); c1sqr *= c1sqr; htsqr =
   * asqr - c1sqr; TRarea = c*sqrt(htsqr)*0.5;
   */
  /*
   * this should work if y1,y2,y3 all have the same sign and x1,x2,x3 are in
   * increasing order
   */

  QUADarea1 = (FABS(y2) + FABS(y1)) * 0.5 * FABS(x2 - x1);
  QUADarea2 = (FABS(y3) + FABS(y2)) * 0.5 * FABS(x3 - x2);
  QUADarea3 = (FABS(y3) + FABS(y1)) * 0.5 * FABS(x3 - x1);
  TRarea = FABS(QUADarea3 - QUADarea1 - QUADarea2);
  area = QUADarea1 + QUADarea2;
  if (area * reltol + abstol > TRarea)
    return (1);
  else
    return (0);
}


/*
 * i is the index of the latest value, a,b,c values correspond to values at
 * t_{i-2}, t{i-1} and t_i
 */

#define SECONDDERIV(i,a,b,c) (oof = (i==ckt->CKTtimeIndex+1?curtime:\
*(ckt->CKTtimePoints+i)),\
(( c - b )/(oof-*(ckt->CKTtimePoints+i-1)) -\
( b - a )/(*(ckt->CKTtimePoints+i-1)-\
*(ckt->CKTtimePoints+i-2)))/(oof - \
*(ckt->CKTtimePoints+i-2)))

/*
 * LTRAlteCalculate - returns sum of the absolute values of the total local
 * truncation error of the 2 equations for the LTRAline
 */


double
LTRAlteCalculate(ckt, genmodel, geninstance, curtime)
  register CKTcircuit *ckt;
  GENmodel *genmodel;
  GENinstance *geninstance;
  double curtime;
{
  LTRAmodel *model = (LTRAmodel *) genmodel;
  LTRAinstance *instance = (LTRAinstance *) geninstance;
  double h1dashTfirstCoeff;
  double h2TfirstCoeff;
  double h3dashTfirstCoeff;
  double dashdash;
  double oof;
  double hilimit1, lolimit1, hivalue1, lovalue1, f1i, g1i;
  double eq1LTE = 0.0, eq2LTE = 0.0;
  int auxindex, tdover, i, exact;

  switch (model->LTRAspecialCase) {

  case LTRA_MOD_LC:
  case LTRA_MOD_RG:
    return (0.0);
    break;