ltramisc.c

linux平台下类似著名的电路板作图软件 Spice的源代码

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,rclsqr);*/
			(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;

    case LTRA_MOD_RLC:

	    if (curtime > model->LTRAtd) {
		    tdover = 1;

		    exact = 0;
		    for (i = ckt->CKTtimeIndex ; i>= 0; i--) {
			    if (curtime - *(ckt->CKTtimePoints + i)
				    ==  model->LTRAtd)
			    {
				    exact =1;
				    break;
			    } 
			    if (curtime - *(ckt->CKTtimePoints + i)
				    > model->LTRAtd)
				break;
		    }

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

		    if (exact == 1) {
			auxindex = i-1;
		    } else {
			auxindex = i;
		    }

	    } else {
		    tdover = 0;
	    }

	    hilimit1 = curtime - *(ckt->CKTtimePoints + ckt->CKTtimeIndex);
	    lolimit1 = 0.0;
	    hivalue1 = LTRArlcH1dashTwiceIntFunc(hilimit1,model->LTRAbeta);
	    lovalue1 = 0.0;

	    f1i = hivalue1;
	    g1i = intlinfunc(lolimit1,hilimit1,lovalue1,hivalue1,
		    lolimit1,hilimit1);
	    h1dashTfirstCoeff = 0.5 * f1i *
		    (curtime - *(ckt->CKTtimePoints+ckt->CKTtimeIndex)) - g1i;

	    if (tdover) {
		hilimit1 = curtime - *(ckt->CKTtimePoints + auxindex);
		lolimit1 = *(ckt->CKTtimePoints + ckt->CKTtimeIndex) - *(ckt->CKTtimePoints + auxindex);
		lolimit1 = MAX(model->LTRAtd,lolimit1);

		/* are the following really doing the operations in the
		write-up? */
		hivalue1 = LTRArlcH2Func(hilimit1,model->LTRAtd,model->LTRAalpha,model->LTRAbeta);
		lovalue1 = LTRArlcH2Func(lolimit1,model->LTRAtd,model->LTRAalpha,model->LTRAbeta);
		f1i = twiceintlinfunc(lolimit1,hilimit1,lolimit1,lovalue1,hivalue1,lolimit1,
			hilimit1);
		g1i = thriceintlinfunc(lolimit1,hilimit1,lolimit1,lolimit1,lovalue1,
			hivalue1,lolimit1,hilimit1);
		h2TfirstCoeff = 0.5*f1i*(curtime-model->LTRAtd- *(ckt->CKTtimePoints+auxindex)) - g1i;

		hivalue1 = LTRArlcH3dashIntFunc(hilimit1,model->LTRAtd,model->LTRAbeta);
		lovalue1 = LTRArlcH3dashIntFunc(lolimit1,model->LTRAtd,model->LTRAbeta);
		f1i = intlinfunc(lolimit1,hilimit1,lovalue1,hivalue1,lolimit1,
			hilimit1);
		g1i = twiceintlinfunc(lolimit1,hilimit1,lolimit1,lovalue1,
			hivalue1,lolimit1,hilimit1);
		h3dashTfirstCoeff = 0.5*f1i*(curtime-model->LTRAtd- *(ckt->CKTtimePoints+auxindex)) - g1i;
	    }


	    /* LTEs for convolution with v1 */
	    /* get divided differences for v1 (2nd derivative estimates) */

	    /* 
	     * no need to subtract operating point values because
	     * taking differences anyway
	     */

	    dashdash = SECONDDERIV(ckt->CKTtimeIndex+1,
			    *(instance->LTRAv1+ckt->CKTtimeIndex-1),
			    *(instance->LTRAv1+ckt->CKTtimeIndex),
			    *(ckt->CKTrhsOld + instance->LTRAposNode1) -
			    *(ckt->CKTrhsOld + instance->LTRAnegNode1));
	    eq1LTE += model->LTRAadmit*FABS(dashdash * 
		    h1dashTfirstCoeff);


	    /* not bothering to interpolate since everything is approximate
	    anyway */

	    if (tdover) {

		    dashdash = SECONDDERIV(auxindex+1,
			    *(instance->LTRAv1 + auxindex - 1),
			    *(instance->LTRAv1 + auxindex),
			    *(instance->LTRAv1 + auxindex + 1)) ;

		    eq2LTE += model->LTRAadmit*FABS(dashdash * 
			    h3dashTfirstCoeff);
		    
		    }
	    /* end LTEs for convolution with v1 */

	    /* LTEs for convolution with v2 */
	    /* get divided differences for v2 (2nd derivative estimates) */

	    dashdash = SECONDDERIV(ckt->CKTtimeIndex+1,
					    *(instance->LTRAv2+ckt->CKTtimeIndex-1),
					    *(instance->LTRAv2+ckt->CKTtimeIndex),
					    *(ckt->CKTrhsOld + instance->LTRAposNode2) -
					    *(ckt->CKTrhsOld + instance->LTRAnegNode2));

	    eq2LTE += model->LTRAadmit*FABS(dashdash * 
		    h1dashTfirstCoeff);


	    if (tdover) {


		    dashdash = SECONDDERIV(auxindex+1,
			    *(instance->LTRAv2 + auxindex - 1),
			    *(instance->LTRAv2 + auxindex),
			    *(instance->LTRAv2 + auxindex + 1)) ;

		    eq1LTE += model->LTRAadmit*FABS(dashdash * 
			    h3dashTfirstCoeff);
		    
	    }

	    /* end LTEs for convolution with v2 */

	    /* LTE for convolution with i1 */
	    /* get divided differences for i1 (2nd derivative estimates) */

	    if (tdover) {


		    dashdash = SECONDDERIV(auxindex+1,
			    *(instance->LTRAi1 + auxindex - 1),
			    *(instance->LTRAi1 + auxindex),
			    *(instance->LTRAi1 + auxindex + 1)) ;

		    eq2LTE += FABS(dashdash * h2TfirstCoeff);
		    

	    }
	    /* end LTE for convolution with i1 */

	    /* LTE for convolution with i2 */
	    /* get divided differences for i2 (2nd derivative estimates) */

	    if (tdover) {

		    dashdash = SECONDDERIV(auxindex+1,