txlsetup.c

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

/**********
Copyright 1992 Regents of the University of California.  All rights
reserved.
Author: 1992 Charles Hough
**********/


#include "ngspice.h"
#include "smpdefs.h"
#include "txldefs.h"
#include "sperror.h"
#include "suffix.h"

static int 		ReadTxL(TXLinstance*, CKTcircuit*);
/*static int 		multC();*/
static int 		main_pade(double, double, double, double, double, TXLine*);
static int 		mac(double, double, double*, double*, double*, double*, double*);
/*static int 		divC();*/
static int 		div_C(double, double, double, double, double*, double*);
static int 		div3(double, double, double, double, double*, double*);
/*static double 	approx1();*/
/*static double 	approx2();*/
static int 		find_roots(double, double, double, double*, double*, double*);
/*static double 	f3();*/
/*static double 	f2();*/
/*static int 		expC();*/
/*static double 	exp_approx1();*/
/*static double 	exp_approx2();*/
static int 		exp_pade(float, float, float, float, float, TXLine*);
/*static int 		exp_div3();*/
static int 		exp_find_roots(double, double, double, double*, double*, double* );
static double		eval2(double, double, double, double);
static int 		get_c(double, double, double, double, double, double, double, double*, double*);
static int 		get_h3(TXLine*);
static int 		Gaussian_Elimination2(int);
static int 		Gaussian_Elimination1(int);
static int 		pade(float);
static int 		update_h1C_c(TXLine *);
static void		y_pade(double, double, double, double, TXLine*);
static double		root3(double, double, double, double);
static NDnamePt 	insert_ND(char*, NDnamePt*);
static NODE 		*insert_node(char*);
static NODE 		*NEW_node(void);
/*static VI_list_txl *new_vi_txl();*/

NODE     		*node_tab = NULL;
NDnamePt 		ndn = NULL;
VI_list_txl 	*pool_vi_txl = NULL;

/* pade.c */
/**
static double xx1, xx2, xx3, xx4, xx5, xx6;
static double cc1, cc2, cc3, cc4, cc5, cc6;
**/

/* y.c */
static double sqtCdL;
static double b1, b2, b3, b4, b5;
static double p1, p2, p3, q1, q2, q3;
static double c1, c2, c3, x1, x2, x3;
static double A[3][4];

/* exp.c */
static double RdL, GdC, RG, tau, RC, GL;
static double a0, a1, a2, a3, a4, a5;
static double ep1, ep2, ep3, eq1, eq2, eq3;
static double ec1, ec2, ec3, ex1, ex2, ex3;
static int    ifImg;
static double AA[3][4];

#define epsi 1.0e-16
#define epsi2 1.0e-28

/* ARGSUSED */
int
TXLsetup(SMPmatrix *matrix, GENmodel *inModel, CKTcircuit*ckt, int *state)
{
  TXLmodel *model = (TXLmodel *)inModel;
  TXLinstance *here;     
  CKTnode *tmp;
  int error;

    /*  loop through all the models */
    for( ; model != NULL; model = model->TXLnextModel ) {

        /* loop through all the instances of the model */
        for (here = model->TXLinstances; here != NULL ;
                here=here->TXLnextInstance) {
            
/* macro to make elements with built in test for out of memory */
#define TSTALLOC(ptr,first,second) \
if((here->ptr = SMPmakeElt(matrix,here->first,here->second))==(double *)NULL){\
    return(E_NOMEM);\
}

			if (! here->TXLibr1Given) {
				error = CKTmkCur(ckt, &tmp, here->TXLname, "branch1");
				if (error) return (error);
				here->TXLibr1 = tmp->number;
			}
			if (! here->TXLibr2Given) {
				error = CKTmkCur(ckt, &tmp, here->TXLname, "branch2");
				if (error) return (error);
				here->TXLibr2 = tmp->number;
			}

            TSTALLOC(TXLposPosptr, TXLposNode, TXLposNode);
            TSTALLOC(TXLposNegptr, TXLposNode, TXLnegNode);
            TSTALLOC(TXLnegPosptr, TXLnegNode, TXLposNode);
            TSTALLOC(TXLnegNegptr, TXLnegNode, TXLnegNode);
            TSTALLOC(TXLibr1Posptr, TXLibr1, TXLposNode);
            TSTALLOC(TXLibr2Negptr, TXLibr2, TXLnegNode);
            TSTALLOC(TXLnegIbr2ptr, TXLnegNode, TXLibr2);
            TSTALLOC(TXLposIbr1ptr, TXLposNode, TXLibr1);
            TSTALLOC(TXLibr1Ibr1ptr, TXLibr1, TXLibr1);
            TSTALLOC(TXLibr2Ibr2ptr, TXLibr2, TXLibr2);
            TSTALLOC(TXLibr1Negptr, TXLibr1, TXLnegNode);
            TSTALLOC(TXLibr2Posptr, TXLibr2, TXLposNode);
            TSTALLOC(TXLibr1Ibr2ptr, TXLibr1, TXLibr2);
            TSTALLOC(TXLibr2Ibr1ptr, TXLibr2, TXLibr1);

			here->in_node_name = CKTnodName(ckt,here->TXLposNode);
			here->out_node_name = CKTnodName(ckt,here->TXLnegNode);
			ReadTxL(here, ckt);

        }
    }

    return(OK);
}

int
TXLunsetup(GENmodel *inModel, CKTcircuit *ckt)
{
  TXLmodel *model;
  TXLinstance *here;
  
  for (model = (TXLmodel *) inModel; model != NULL;
      model = model->TXLnextModel) {
    for (here = model->TXLinstances; here != NULL;
        here = here->TXLnextInstance) {
  
          if (here->TXLibr1) {
               CKTdltNNum(ckt, here->TXLibr1);
               here->TXLibr1 = 0;
          }

          if (here->TXLibr2) {
               CKTdltNNum(ckt, here->TXLibr2);
               here->TXLibr2 = 0;
          }

    here->TXLdcGiven=0;

    }
  }
  return OK;
}

/***
static VI_list_txl
*new_vi_txl(void)
{
	VI_list_txl *q;

	if (pool_vi_txl) {
		q = pool_vi_txl;
		pool_vi_txl = pool_vi_txl->pool;
		return(q);
	} else
		return((VI_list_txl *) tmalloc (sizeof (VI_list_txl)));
}
***/

static int 
ReadTxL(TXLinstance *tx, CKTcircuit *ckt)
{
   double R, L, G, C, l;
   char *p, *n;
   NODE *nd;
   ETXLine *et;
   TXLine *t, *t2;
   RLINE *line;
   ERLINE *er;
   double LL = 1e-12;
	
   p = tx->in_node_name;
   n = tx->out_node_name;

   line = (RLINE *) tmalloc(sizeof (RLINE));
   er = (ERLINE *) tmalloc(sizeof (ERLINE));
   et = (ETXLine *) tmalloc(sizeof (ETXLine));
   t = (TXLine *) tmalloc(sizeof (TXLine));
   t2 = (TXLine *) tmalloc(sizeof (TXLine));
   tx->txline = t;
   tx->txline2 = t2;
   t->newtp = 0;
   t2->newtp = 0;
   t->vi_head = t->vi_tail = NULL;
   nd = insert_node(p);
   et->link = nd->tptr;
   nd->tptr = et;
   et->line = t;
   t->in_node = nd;
   t2->in_node = nd;
   er->link = nd->rlptr;
   nd->rlptr = er;
   er->rl = line;
   line->in_node = nd;
   et = (ETXLine *) tmalloc(sizeof (ETXLine));
   nd = insert_node(n);
   et->link = nd->tptr;
   nd->tptr = et;
   et->line = t;
   t->out_node = nd;
   t2->out_node = nd;
   er = (ERLINE *) tmalloc(sizeof (ERLINE));
   er->link = nd->rlptr;
   nd->rlptr = er;
   er->rl = line;
   line->out_node = nd;
   t->dc1 = t->dc2 = 0.0;
   t2->dc1 = t2->dc2 = 0.0;
   t->lsl = 0;
   t2->lsl = 0;
   l = 0.0;

   R = tx->TXLmodPtr->R;
   L = tx->TXLmodPtr->L;
   L = MAX(L, LL);
   C = tx->TXLmodPtr->C;
   G = tx->TXLmodPtr->G;
   if (tx->TXLlengthgiven == TRUE)
   		l = tx->TXLlength;
   else l = tx->TXLmodPtr->length;	


   if (l == 0.0) {
	   fprintf(stderr, "(Error) transmission line of zero length\n");
	   exit(0);
   }
   else {
		if (R / L  < 5.0e+5) {
			line->g = 1.0e+2;
			if (G < 1.0e-2) {
				t->lsl = 1;  /* lossless line */
				t->taul = sqrt((double) C * L) * l * 1.0e12;
				t->h3_aten = t->sqtCdL = sqrt((double) C / L);
				t->h2_aten = 1.0;
				t->h1C = 0.0;
			}
		}
		else line->g = 1.0 / (R * l);
	}

   if (! t->lsl)
   		main_pade(R, L, G, C, l, t);

   return(1);
}


/****************************************************************
     pade.c  :  Calculate the Pade Approxximation of Y(s)
 ****************************************************************/


static int 
main_pade(double R, double L, double G, double C, double l, TXLine *h)
{
   y_pade(R, L, G, C, h);
   h->ifImg = exp_pade(R, L, G, C, l, h);
   get_h3(h);
   h->taul *= 1.0e12;
   update_h1C_c(h);

   return(1);
}

static int 
div_C(double ar, double ai, double br, double bi, double *cr, double *ci)
{
   *cr = ar * br + ai * bi;
   *ci = - ar * bi + ai * br;
   *cr = *cr / (br * br + bi * bi);
   *ci = *ci / (br * br + bi * bi);
   return (1);
}

/***
static int expC(ar, ai, h, cr, ci)
   double ar, ai, *cr, *ci;
   float h;
{
   double e, cs, si;

   e = exp((double) ar * h);
   cs = cos((double) ai * h);
   si = sin((double) ai * h);
   *cr = e * cs;
   *ci = e * si;

   return(1);
}
***/

/***
static int multC(ar, ai, br, bi, cr, ci)
   double ar, ai, br, bi;
   double *cr, *ci;
{
   *cr = ar*br - ai*bi;
   *ci = ar*bi + ai*br;

   return(1);
}
***/

/***
static int divC(ar, ai, br, bi, cr, ci)
   double ar, ai, br, bi;
   double *cr, *ci;
{
   double t;
   t = br*br + bi*bi;
   *cr = (ar*br + ai*bi) / t;
   *ci = (ai*br - ar*bi) / t;

   return(1);
}
***/

static int 
get_h3(TXLine *h)
{
   double cc1, cc2, cc3, cc4, cc5, cc6;
   double xx1, xx2, xx3, xx4, xx5, xx6;

   h->h3_aten = h->h2_aten * h->sqtCdL;
   h->h3_term[0].x = xx1 = h->h1_term[0].x;
   h->h3_term[1].x = xx2 = h->h1_term[1].x;
   h->h3_term[2].x = xx3 = h->h1_term[2].x;
   h->h3_term[3].x = xx4 = h->h2_term[0].x;
   h->h3_term[4].x = xx5 = h->h2_term[1].x;
   h->h3_term[5].x = xx6 = h->h2_term[2].x;
   cc1 = h->h1_term[0].c;
   cc2 = h->h1_term[1].c;
   cc3 = h->h1_term[2].c;
   cc4 = h->h2_term[0].c;
   cc5 = h->h2_term[1].c;
   cc6 = h->h2_term[2].c;

   if (h->ifImg) {
      double r, i;

      h->h3_term[0].c = cc1 + cc1 * (cc4/(xx1-xx4) + 
         2.0*(cc5*xx1-xx6*cc6-xx5*cc5)/(xx1*xx1-2.0*xx5*xx1+xx5*xx5+xx6*xx6));
      h->h3_term[1].c = cc2 + cc2 * (cc4/(xx2-xx4) + 
         2.0*(cc5*xx2-xx6*cc6-xx5*cc5)/(xx2*xx2-2.0*xx5*xx2+xx5*xx5+xx6*xx6));
      h->h3_term[2].c = cc3 + cc3 * (cc4/(xx3-xx4) + 
         2.0*(cc5*xx3-xx6*cc6-xx5*cc5)/(xx3*xx3-2.0*xx5*xx3+xx5*xx5+xx6*xx6));

      h->h3_term[3].c = cc4 + cc4 * (cc1/(xx4-xx1) + cc2/(xx4-xx2) + cc3/(xx4-xx3));

      h->h3_term[4].c = cc5;
      h->h3_term[5].c = cc6;
      div_C(cc5, cc6, xx5-xx1, xx6, &r, &i);
      h->h3_term[4].c += r * cc1;
      h->h3_term[5].c += i * cc1;
      div_C(cc5, cc6, xx5-xx2, xx6, &r, &i);
      h->h3_term[4].c += r * cc2;
      h->h3_term[5].c += i * cc2;
      div_C(cc5, cc6, xx5-xx3, xx6, &r, &i);
      h->h3_term[4].c += r * cc3;
      h->h3_term[5].c += i * cc3;
   } else {
      h->h3_term[0].c = cc1 + cc1 * (cc4/(xx1-xx4) + cc5/(xx1-xx5) + cc6/(xx1-xx6));
      h->h3_term[1].c = cc2 + cc2 * (cc4/(xx2-xx4) + cc5/(xx2-xx5) + cc6/(xx2-xx6));
      h->h3_term[2].c = cc3 + cc3 * (cc4/(xx3-xx4) + cc5/(xx3-xx5) + cc6/(xx3-xx6));

      h->h3_term[3].c = cc4 + cc4 * (cc1/(xx4-xx1) + cc2/(xx4-xx2) + cc3/(xx4-xx3));
      h->h3_term[4].c = cc5 + cc5 * (cc1/(xx5-xx1) + cc2/(xx5-xx2) + cc3/(xx5-xx3));
      h->h3_term[5].c = cc6 + cc6 * (cc1/(xx6-xx1) + cc2/(xx6-xx2) + cc3/(xx6-xx3));
   }
      
   return(1);
}

static int 
update_h1C_c(TXLine *h)
{
   int i;
   double d = 0;

   for (i = 0; i < 3; i++) {
      h->h1_term[i].c *= h->sqtCdL;
      d += h->h1_term[i].c;
   }
   h->h1C = d;

   for (i = 0; i < 3; i++) 
      h->h2_term[i].c *= h->h2_aten;
   
   for (i = 0; i < 6; i++) 
      h->h3_term[i].c *= h->h3_aten;

   return(1);
}
/****************************************************************
     y.c  :  Calculate the Pade Approximation of Y(s)
 ****************************************************************/


static double 
eval2(double a, double b, double c, double x)
{
   return(a*x*x + b*x + c);
}

/***
static double approx1(st)
   double st;
{
   double s3, s2, s1;

   s1 = st;
   s2 = s1 * s1;
   s3 = s2 * s1;

   return((s3 + q1*s2 + q2*s1 + q3) / (s3 + p1*s2 + p2*s1 + p3));
}
***/
/***
static double approx2(st)
   double st;
{
   return(1.0 + c1/(st - x1) + c2/(st - x2) + c3/(st - x3));
}
***/

static void 
y_pade(double R, double L, double G, double C, TXLine *h)
{

   /* float RdL, GdC; */
   double RdL, GdC; 

   sqtCdL = sqrt((double) C / L);
   RdL = R / L;
   GdC = G / C;

   mac(GdC, RdL, &b1, &b2, &b3, &b4, &b5);

   A[0][0] = 1.0 - sqrt((double) (GdC / RdL));
   A[0][1] = b1;
   A[0][2] = b2;
   A[0][3] = -b3;

   A[1][0] = b1;
   A[1][1] = b2;
   A[1][2] = b3;
   A[1][3] = -b4;

   A[2][0] = b2;
   A[2][1] = b3;
   A[2][2] = b4;
   A[2][3] = -b5;

   Gaussian_Elimination1(3);

   p3 = A[0][3];
   p2 = A[1][3];
   p1 = A[2][3];

   q1 = p1 + b1;
   q2 = b1 * p1 + p2 + b2;
   q3 = p3 * sqrt((double) (GdC / RdL));

   find_roots(p1, p2, p3, &x1, &x2, &x3);
   c1 = eval2(q1 - p1, q2 - p2, q3 - p3, x1) / 
		   eval2((double) 3.0, (double) 2.0 * p1, p2, x1);
   c2 = eval2(q1 - p1, q2 - p2, q3 - p3, x2) / 
		   eval2((double) 3.0, (double) 2.0 * p1, p2, x2);
   c3 = eval2(q1 - p1, q2 - p2, q3 - p3, x3) / 
		   eval2((double) 3.0, (double) 2.0 * p1, p2, x3);

   h->sqtCdL = sqtCdL;
   h->h1_term[0].c = c1;
   h->h1_term[1].c = c2;
   h->h1_term[2].c = c3;
   h->h1_term[0].x = x1;
   h->h1_term[1].x = x2;
   h->h1_term[2].x = x3;

}

static int 
Gaussian_Elimination1(int dims)
{
   int i, j, k, dim;
   double f;
   int imax;
   double max;

   dim = dims;

   for (i = 0; i < dim; i++) {
      imax = i;
      max = ABS(A[i][i]);
      for (j = i+1; j < dim; j++)
         if (ABS(A[j][i]) > max) {
	    imax = j;
	    max = ABS(A[j][i]);
	 } 
      if (max < epsi) {
         fprintf(stderr, " can not choose a pivot \n");
         exit(0);
      }
      if (imax != i)
	 for (k = i; k <= dim; k++) {
	    f = A[i][k];
	    A[i][k] = A[imax][k];
	    A[imax][k] = f;
	 }
      
      f = 1.0 / A[i][i];
      A[i][i] = 1.0;

      for (j = i+1; j <= dim; j++)
	 A[i][j] *= f;

      for (j = 0; j < dim ; j++) {
	 if (i == j)
	    continue;
	 f = A[j][i];
	 A[j][i] = 0.0;
	 for (k = i+1; k <= dim; k++)
	    A[j][k] -= f * A[i][k];
      }
   }
   return(1);
}

static double 
root3(double a1, double a2, double a3, double x)
{
   double t1, t2;

   t1 = x*x*x + a1*x*x + a2*x + a3;
   t2 = 3.0*x*x + 2.0*a1*x + a2;

   return(x - t1 / t2);
}

static int 
div3(double a1, double a2, double a3, double x, double *p1, double *p2)
{
   *p1 = a1 + x;
   *p2 = - a3 / x;

   return(1);
}


/****************************************************************
         Calculate the Maclaurin series of F(z)

   F(z) = sqrt((1+az) / (1+bz))
	= 1 + b1 z + b2 z^2 + b3 z^3 + b4 z^4 + b5 z^5 
 ****************************************************************/

/***