ltramisc.c

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

/**********
Copyright 1990 Regents of the University of California.  All rights reserved.
Author: 1990 Jaijeet S. Roychowdhury
**********/

#include "spice.h"
#include "ltradefs.h"
#include "util.h"
#include "suffix.h"

/*
 * Miscellaneous functions to do with lossy lines
 */

/*
 * LTRAquadInterp - quadratic interpolation function t = timepoint where
 * value wanted t1, t2, t3 are three timepoints where the value is known c1,
 * c2, c3 are set to the proper coefficients by the function the interpolated
 * value is c1*v1 + c2*v2 + c3*v3; this should be done in the calling
 * program; (v1,v2,v3 are the known values at t1,t2,t3)
 */


int
LTRAquadInterp(t, t1, t2, t3, c1, c2, c3)
  double t, t1, t2, t3;
  double *c1, *c2, *c3;
{
  double f1, f2, f3;

  if (t == t1) {
    *c1 = 1.0;
    *c2 = 0.0;
    *c3 = 0.0;
    return (0);
  }
  if (t == t2) {
    *c1 = 0.0;
    *c2 = 1.0;
    *c3 = 0.0;
    return (0);
  }
  if (t == t3) {
    *c1 = 0.0;
    *c2 = 0.0;
    *c3 = 1.0;
    return (0);
  }
  if ((t2 - t1) == 0 || (t3 - t2) == 0 || (t1 - t3) == 0)
    return (1);

  f1 = (t - t2) * (t - t3);
  f2 = (t - t1) * (t - t3);
  f3 = (t - t1) * (t - t2);
  if ((t2 - t1) == 0) {		/* should never happen, but don't want to
				 * divide by zero, EVER... */
    f1 = 0;
    f2 = 0;
  } else {
    f1 /= (t1 - t2);
    f2 /= (t2 - t1);
  }
  if ((t3 - t2) == 0) {		/* should never happen, but don't want to
				 * divide by zero, EVER... */
    f2 = 0;
    f3 = 0;
  } else {
    f2 /= (t2 - t3);
    f3 /= (t2 - t3);
  }
  if ((t3 - t1) == 0) {		/* should never happen, but don't want to
				 * divide by zero, EVER... */
    f1 = 0;
    f2 = 0;
  } else {
    f1 /= (t1 - t3);
    f3 /= (t1 - t3);
  }
  *c1 = f1;
  *c2 = f2;
  *c3 = f3;
  return (0);
}

/* linear interpolation */

int
LTRAlinInterp(t, t1, t2, c1, c2)
  double t, t1, t2;
  double *c1, *c2;

{
  double temp;

  if (t1 == t2)
    return (1);

  if (t == t1) {
    *c1 = 1.0;
    *c2 = 0.0;
    return (0);
  }
  if (t == t2) {
    *c1 = 0.0;
    *c2 = 1.0;
    return (0);
  }
  temp = (t - t1) / (t2 - t1);
  *c2 = temp;
  *c1 = 1 - temp;
  return (0);
}

/*
 * intlinfunc returns \int_lolimit^hilimit h(\tau) d \tau, where h(\tau) is
 * assumed to be linear, with values lovalue and hivalue \tau = t1 and t2
 * respectively this is used only locally
 */

double
intlinfunc(lolimit, hilimit, lovalue, hivalue, t1, t2)
  double lolimit, hilimit, lovalue, hivalue, t1, t2;
{
  double width, m;


  width = t2 - t1;
  if (width == 0.0)
    return (0.0);
  m = (hivalue - lovalue) / width;

  return ((hilimit - lolimit) * lovalue + 0.5 * m * ((hilimit - t1) * (hilimit - t1)
	  - (lolimit - t1) * (lolimit - t1)));
}


/*
 * twiceintlinfunc returns \int_lolimit^hilimit \int_otherlolimit^\tau
 * h(\tau') d \tau' d \tau , where h(\tau') is assumed to be linear, with
 * values lovalue and hivalue \tau = t1 and t2 respectively this is used only
 * locally
 */

double
twiceintlinfunc(lolimit, hilimit, otherlolimit, lovalue, hivalue, t1, t2)
  double lolimit, hilimit, lovalue, hivalue, t1, t2, otherlolimit;
{
  double width, m, dummy;
  double temp1, temp2, temp3;


  width = t2 - t1;
  if (width == 0.0)
    return (0.0);
  m = (hivalue - lovalue) / width;

  temp1 = hilimit - t1;
  temp2 = lolimit - t1;
  temp3 = otherlolimit - t1;
  dummy = lovalue * ((hilimit - otherlolimit) * (hilimit - otherlolimit) -
      (lolimit - otherlolimit) * (lolimit - otherlolimit));
  dummy += m * ((temp1 * temp1 * temp1 - temp2 * temp2 * temp2) / 3.0 -
      temp3 * temp3 * (hilimit - lolimit));
  return (dummy * 0.5);
}

/*
 * thriceintlinfunc returns \int_lolimit^hilimit \int_secondlolimit^\tau
 * \int_thirdlolimit^\tau' h(\tau'') d \tau'' d \tau' d \tau , where
 * h(\tau'') is assumed to be linear, with values lovalue and hivalue \tau =
 * t1 and t2 respectively this is used only locally
 */

double
thriceintlinfunc(lolimit, hilimit, secondlolimit, thirdlolimit, lovalue,
    hivalue, t1, t2)
  double lolimit, hilimit, lovalue, hivalue, t1, t2, secondlolimit, thirdlolimit;
{
  double width, m, dummy;
  double temp1, temp2, temp3, temp4;
  double temp5, temp6, temp7, temp8, temp9, temp10;


  width = t2 - t1;
  if (width == 0.0)
    return (0.0);
  m = (hivalue - lovalue) / width;

  temp1 = hilimit - t1;
  temp2 = lolimit - t1;
  temp3 = secondlolimit - t1;
  temp4 = thirdlolimit - t1;
  temp5 = hilimit - thirdlolimit;
  temp6 = lolimit - thirdlolimit;
  temp7 = secondlolimit - thirdlolimit;
  temp8 = hilimit - lolimit;
  temp9 = hilimit - secondlolimit;
  temp10 = lolimit - secondlolimit;
  dummy = lovalue * ((temp5 * temp5 * temp5 - temp6 * temp6 * temp6) / 3 -
      temp7 * temp5 * temp8);
  dummy += m * (((temp1 * temp1 * temp1 * temp1 - temp2 * temp2 * temp2 * temp2) * 0.25 -
	  temp3 * temp3 * temp3 * temp8) / 3 - temp4 * temp4 * 0.5 * (temp9 * temp9 -
	  temp10 * temp10));
  return (dummy * 0.5);
}


/*
 * These are from the book Numerical Recipes in C
 * 
 */

double
bessI0(x)
  double x;
{
  double ax, ans;
  double y;

  if ((ax = fabs(x)) < 3.75) {
    y = x / 3.75;
    y *= y;
    ans = 1.0 + y * (3.5156229 + y * (3.0899424 + y * (1.2067492
		+ y * (0.2659732 + y * (0.360768e-1 + y * 0.45813e-2)))));
  } else {
    y = 3.75 / ax;
    ans = (exp(ax) / sqrt(ax)) * (0.39894228 + y * (0.1328592e-1
	    + y * (0.225319e-2 + y * (-0.157565e-2 + y * (0.916281e-2
			+ y * (-0.2057706e-1 + y * (0.2635537e-1 + y * (-0.1647633e-1
				    + y * 0.392377e-2))))))));
  }
  return (ans);
}

double
bessI1(x)
  double x;
{
  double ax, ans;
  double y;

  if ((ax = fabs(x)) < 3.75) {
    y = x / 3.75;
    y *= y;
    ans = ax * (0.5 + y * (0.87890594 + y * (0.51498869 + y * (0.15084934
		    + y * (0.2658733e-1 + y * (0.301532e-2 + y * 0.32411e-3))))));
  } else {
    y = 3.75 / ax;
    ans = 0.2282967e-1 + y * (-0.2895312e-1 + y * (0.1787654e-1
	    - y * 0.420059e-2));
    ans = 0.39894228 + y * (-0.3988024e-1 + y * (-0.362018e-2
	    + y * (0.163801e-2 + y * (-0.1031555e-1 + y * ans))));
    ans *= (exp(ax) / sqrt(ax));
  }
  return (x < 0.0 ? -ans : ans);
}

double
bessI1xOverX(x)
  double x;
{
  double ax, ans;
  double y;

  if ((ax = fabs(x)) < 3.75) {
    y = x / 3.75;
    y *= y;
    ans = 0.5 + y * (0.87890594 + y * (0.51498869 + y * (0.15084934
		+ y * (0.2658733e-1 + y * (0.301532e-2 + y * 0.32411e-3)))));
  } else {
    y = 3.75 / ax;
    ans = 0.2282967e-1 + y * (-0.2895312e-1 + y * (0.1787654e-1
	    - y * 0.420059e-2));
    ans = 0.39894228 + y * (-0.3988024e-1 + y * (-0.362018e-2
	    + y * (0.163801e-2 + y * (-0.1031555e-1 + y * ans))));
    ans *= (exp(ax) / (ax * sqrt(ax)));
  }
  return (ans);
}

/* LTRArlcH1dashFunc - the first impulse response function */

double
LTRArlcH1dashFunc(time, T, alpha, beta)
  double time, T, alpha, beta;

{
  double besselarg, exparg, returnval;
  /* T is not used in this function */

  /*
   * result = alpha * e^{- beta*time} * {I_1(alpha*time) - I_0(alpha*time)}
   */

  if (alpha == 0.0)
    return (0.0);

  exparg = -beta * time;
  besselarg = alpha * time;

  returnval = (bessI1(besselarg) - bessI0(besselarg)) * alpha * exp(exparg);
  return (returnval);
}

double
LTRArlcH2Func(time, T, alpha, beta)
  double time, T, alpha, beta;

{
  double besselarg, exparg, returnval;

  /*
   * result = 0, time < T = (alpha*T*e^{-beta*time})/sqrt(t^2 - T^2) *
   * I_1(alpha*sqrt(t^2 - T^2)), time >= T
   */

  if (alpha == 0.0)
    return (0.0);
  if (time < T)
    return (0.0);

  if (time != T) {
    besselarg = alpha * sqrt(time * time - T * T);
  } else {
    besselarg = 0.0;
  }
  exparg = -beta * time;

  returnval = alpha * alpha * T * exp(exparg) * bessI1xOverX(besselarg);
  return (returnval);
}

double
LTRArlcH3dashFunc(time, T, alpha, beta)
  double time, T, alpha, beta;

{
  double exparg, besselarg, returnval;

  /*
   * result = 0, time < T = alpha*e^{-beta*time}*(t/sqrt(t^2-T^2)*
   * I_1(alpha*sqrt(t^2-T^2)) - I_0(alpha*sqrt(t^2-T^2)))
   */

  if (alpha == 0.0)
    return (0.0);
  if (time < T)
    return (0.0);

  exparg = -beta * time;
  if (time != T) {
    besselarg = alpha * sqrt(time * time - T * T);
  } else {
    besselarg = 0.0;
  }

  returnval = alpha * time * bessI1xOverX(besselarg) - bessI0(besselarg);
  returnval *= alpha * exp(exparg);
  return (returnval);
}

/*
 * LTRArlcH1dashTwiceIntFunc - twice repeated integral of h1dash for the
 * special case of G = 0
 */

double
LTRArlcH1dashTwiceIntFunc(time, beta)
  double time, beta;
{
  double arg, returnval;

  /*
   * result = time * e^{- beta*time} * {I_0(beta*time) + I_1(beta*time)} -
   * time
   */

  if (beta == 0.0)
    return (time);
  arg = beta * time;
  if (arg == 0.0)
    return (0.0);

  returnval = (bessI1(arg) + bessI0(arg)) * time * exp(-arg) - time;
  return (returnval);
}

/*
 * LTRArlcH3dashIntFunc - twice repeated integral of h1dash for the special
 * case of G = 0
 */

double
LTRArlcH3dashIntFunc(time, T, beta)
  double time, T, beta;

{
  double exparg, besselarg;
  double returnval;

  if (time <= T)
    return (0.0);
  if (beta == 0.0)
    return (0.0);
  exparg = -beta * time;
  besselarg = beta * sqrt(time * time - T * T);
  returnval = exp(exparg) * bessI0(besselarg) - exp(-beta * T);
  return (returnval);
}


double
LTRArcH1dashTwiceIntFunc(time, cbyr)
  double time, cbyr;
{

  return (sqrt(4 * cbyr * time / M_PI));
}




double
LTRArcH2TwiceIntFunc(time, rclsqr)
  double time, rclsqr;
{
  double temp;

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



double
LTRArcH3dashTwiceIntFunc(time, cbyr, rclsqr)