ltramisc.c

linux平台下类似著名的电路板作图软件 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)
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,rclsqr); */
(hilimit1 != 0.0? temp = 2*sqrt(hilimit1/M_PI)*temp3 - temp4*temp2, (temp5*temp): 0.0);