txlsetup.c

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

static double f3(double a, double b, double z)
{
   double t4, t3, t2, t1;
   double t14, t13, t12, t11;
   double sqt11;

   t1 = 1 / (1.0 + b * z);
   t2 = t1 * t1;
   t3 = t2 * t1;
   t4 = t3 * t1;

   t11 = (1.0 + a * z) * t1;
   t12 = (1.0 + a * z) * t2;
   t13 = (1.0 + a * z) * t3;
   t14 = (1.0 + a * z) * t4;

   sqt11 = sqrt(t11);


   return(
     -0.5 * (-2.0*a*b*t2 + 2.0*b*b*t13) * (a*t1 - b*t12) / (t11*sqt11)
     +3.0/8.0 * (a*t1-b*t12)*(a*t1-b*t12)*(a*t1-b*t12) / (t11*t11*sqt11)
     +0.5 * (4.0*a*b*b*t3 + 2.0*a*b*b*t3 - 6.0*b*b*b*t14) / sqt11
     -0.25 * (-2.0*a*b*t2 + 2.0*b*b*t13) * (a*t1-b*(1.0+a*z)) / 
     (t11*sqt11)
     );
}
***/

/***
static double f2(a, b, z)
   double a, b, z;
{
   double t3, t2, t1;
   double t13, t12, t11;
   double sqt11;

   t1 = 1 / (1.0 + b * z);
   t2 = t1 * t1;
   t3 = t2 * t1;

   t11 = (1.0 + a * z) * t1;
   t12 = (1.0 + a * z) * t2;
   t13 = (1.0 + a * z) * t3;

   sqt11 = sqrt(t11);

   return(
     -0.25 * (a*t1-b*t12) * (a*t1-b*t12) / (t11*sqt11)
     +0.5 * (-2.0*a*b*t2 + 2.0*b*b*t13) / sqt11
   );
}   
***/

static int 
mac(double at, double bt, double *b1, double *b2, double *b3, double *b4, double *b5)
   /* float at, bt; */
{
   double a, b;
   double y1, y2, y3, y4, y5;

   a = at;
   b = bt;

   y1 = *b1 = 0.5 * (a - b);
   y2 = 0.5 * (3.0 * b * b - 2.0 * a * b - a * a) * y1 / (a - b);
   y3 = ((3.0 * b * b + a * a) * y1 * y1 + 0.5 * (3.0 * b * b 
	- 2.0 * a * b - a * a) * y2) / (a - b);
   y4 = ((3.0 * b * b - 3.0 * a * a) * y1 * y1 * y1 + (9.0 * b * b
	+ 3.0 * a * a) * y1 * y2 + 0.5 * (3.0 * b * b - 2.0 * a * b
	- a * a) * y3) / (a - b);
   y5 = (12.0 * a * a * y1 * y1 * y1 * y1 + y1 * y1 * y2 * (
	 18.0 * b * b - 18.0 * a * a) + (9.0 * b * b + 3.0 * a * a) *
	 (y2 * y2 + y1 * y3) + (3.0 * b * b + a * a) * y1 * y3 +
	 0.5 * (3.0 * b * b - 2.0 * a * b - a * a) * y4) / (a - b);

   *b2 = y2 / 2.0;
   *b3 = y3 / 6.0;
   *b4 = y4 / 24.0;
   *b5 = y5 / 120.0;

   return(1);
}


/****************************************************
 exp.c
 ****************************************************/

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

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

   return(exp((double) - st * tau - a0) * 
	(s3 + eq1*s2 + eq2*s1 + eq3) / (s3 + ep1*s2 + ep2*s1 + ep3));
}
***/

static int 
get_c(double eq1, double eq2, double eq3, double ep1, double ep2, double a, double b, 
      double *cr, double *ci)
{
   double d, n;

   d = (3.0*(a*a-b*b)+2.0*ep1*a+ep2)*(3.0*(a*a-b*b)+2.0*ep1*a+ep2);
   d += (6.0*a*b+2.0*ep1*b)*(6.0*a*b+2.0*ep1*b);
   n = -(eq1*(a*a-b*b)+eq2*a+eq3)*(6.0*a*b+2.0*ep1*b);
   n += (2.0*eq1*a*b+eq2*b)*(3.0*(a*a-b*b)+2.0*ep1*a+ep2);
   *ci = n/d;
   n = (3.0*(a*a-b*b)+2.0*ep1*a+ep2)*(eq1*(a*a-b*b)+eq2*a+eq3);
   n += (6.0*a*b+2.0*ep1*b)*(2.0*eq1*a*b+eq2*b);
   *cr = n/d;

   return(1);
}

/***
static double exp_approx2(double st)
{
   if (ifImg) 
      return(1.0 + ec1/(st - ex1) + 2.0*(ec2*(st-ex2)-ec3*ex3) /
	 ((st-ex2)*(st-ex2) + ex3*ex3));
   else 
      return(1.0 + ec1/(st - ex1) + ec2/(st - ex2) + ec3/(st - ex3));
}
***/

static int 
exp_pade(float R, float L, float G, float C, float l, TXLine *h)
{
  
   tau = sqrt((double) L*C);
   RdL = R / L;
   GdC = G / C;
   RG  = R * G;
   RC = R * C;
   GL = G * L;

   {
      double a, b, t;
      double y1, y2, y3, y4, y5, y6;

      a = RdL;
      b = GdC;
      t = tau;

      /*
      y1 = 0.5 * (a + b);
      y2 = a * b - y1 * y1;
      y3 = - a * b * y1 - 2.0 * y1 * y2 + y1 * y1 * y1;
      y4 = 2.0 * a * b * y1 * y1 - a * b * y2 - 2.0 * y2 * y2
          - 2.0 * y1 * y3 + 5.0 * y1 * y1 * y2
          - 2.0 * y1 * y1 * y1 * y1;
      y5 = 6.0 * a * b * (y1 * y2 - y1 * y1 * y1) - a * b * y3
           - 2.0 * y1 * y4
           - 6.0 * y2 * y3 + 12.0 * y2 * y2 * y1 + 7.0 * y1 * y1 * y3
           -10.0 * y1 * y1 * y1 * y2 - 8.0 * y1 * y1 * y1 * y2
           + 6.0 * y1 * y1 * y1 * y1 * y1;
   y6 = 24.0 * a * b * y1 * y1 * y1 * y1 - 36.0 * a * b * y1 * y1 * y2
        + 6.0 * a * b * y2 * y2 + 8.0 * a * b * y1 * y3 - 2.0 * y2 * y4
        - 2.0 * y1 * y5 + 2.0 * y1 * y1 * y4 - a * b * y4 -6.0 * y3 * y3
        + 44.0 * y1 * y2 * y3 + 60.0 * y1 * y1 * y1 * y1 * y2
	-24.0 * y1 * y1 * y1 * y1 * y1 * y1 + 12.0 * y2 * y2 * y2
	-54.0 * y1 * y1 * y2 * y2 + 7.0 * y1 * y1 * y4
	-24.0 * y1 * y1 * y1 * y3 - 24.0 * y1 * y1 * y2 * y2
	-8.0 * y1 * y1 * y1 * y3 + 24.0 * y1 * y1 * y1 * y1 * y2
	- 6.0 * y2 * y4;
	*/

       y1 = 0.5 * (a + b);
       y2 = a * b - y1 * y1;
       y3 = -3.0 * y1 * y2;
       y4 = -3.0 * y2 * y2 - 4.0 * y1 * y3;
       y5 = - 5.0 * y1 * y4 -10.0 * y2 * y3;
       y6 = -10.0 * y3 * y3 - 15.0 * y2 * y4 - 6.0 * y1 * y5;

       a0 = y1 * t;
       a1 = y2 * t * t / 2.0;
       a2 = y3 * t * t * t / 6.0;
       a3 = y4 * t * t * t * t / 24.0;
       a4 = y5 * t * t * t * t * t / 120.0;
       a5 = y6 * t * t * t * t * t * t / 720.0;

   }

   a0 *= l;
   a1 *= l;
   a2 *= l;
   a3 *= l;
   a4 *= l;
   a5 *= l;

   pade(l);

      h->taul = tau * l;
      h->h2_aten = exp(- a0);
      h->h2_term[0].c = ec1;
      h->h2_term[1].c = ec2;
      h->h2_term[2].c = ec3;
      h->h2_term[0].x = ex1;
      h->h2_term[1].x = ex2;
      h->h2_term[2].x = ex3;

   return(ifImg);
}

static int pade(float l)
{
   int i, j;
   double a[6];
   double b[6];

   a[1] = -a1;
   a[2] = -a2;
   a[3] = -a3;
   a[4] = -a4;
   a[5] = -a5;
   
   b[0] = 1.0;
   b[1] = a[1];
   for (i = 2; i <= 5; i++) {
      b[i] = 0.0;
      for (j = 1; j <= i; j++)
	 b[i] += j * a[j] * b[i-j];
      b[i] = b[i] / (double) i;
   }

   AA[0][0] = 1.0 - exp((double) a0 - l * sqrt(RG));
   AA[0][1] = b[1];
   AA[0][2] = b[2];
   AA[0][3] = -b[3];

   AA[1][0] = b[1];
   AA[1][1] = b[2];
   AA[1][2] = b[3];
   AA[1][3] = -b[4];

   AA[2][0] = b[2];
   AA[2][1] = b[3];
   AA[2][2] = b[4];
   AA[2][3] = -b[5];

   Gaussian_Elimination2(3);

   ep3 = AA[0][3];
   ep2 = AA[1][3];
   ep1 = AA[2][3];

   eq1 = ep1 + b[1];
   eq2 = b[1] * ep1 + ep2 + b[2];
   eq3 = ep3 * exp((double) a0 - l * sqrt(RG));

   ep3 = ep3 / (tau*tau*tau);
   ep2 = ep2 / (tau*tau);
   ep1 = ep1 / tau;
   eq3 = eq3 / (tau*tau*tau);
   eq2 = eq2 / (tau*tau);
   eq1 = eq1 / tau;
    /*
   printf("factor = %e\n", exp(-a0));
   printf("ep1 = %e ep2 = %e ep3 = %e\n", ep1, ep2, ep3);
     */
   exp_find_roots(ep1, ep2, ep3, &ex1, &ex2, &ex3);
    /*
   printf("roots are %e %e %e \n", ex1, ex2, ex3);
     */
   ec1 = eval2(eq1 - ep1, eq2 - ep2, eq3 - ep3, ex1) / 
		   eval2((double) 3.0, (double) 2.0 * ep1, ep2, ex1);
   if (ifImg) 
      get_c(eq1 - ep1, eq2 - ep2, eq3 - ep3, ep1, ep2, ex2, ex3, &ec2, &ec3);
   else {
      ec2 = eval2(eq1 - ep1, eq2 - ep2, eq3 - ep3, ex2) / 
		   eval2((double) 3.0, (double) 2.0 * ep1, ep2, ex2);
      ec3 = eval2(eq1 - ep1, eq2 - ep2, eq3 - ep3, ex3) / 
		   eval2((double) 3.0, (double) 2.0 * ep1, ep2, ex3);
   }
   return (1);
}

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

   dim = dims;

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

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

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

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

   return(1);
}
***/

/***
 ***/

static int 
exp_find_roots(double a1, double a2, double a3, double *ex1, double *ex2, double *ex3)
{
   double x, t;
   double p, q;

   q = (a1*a1-3.0*a2) / 9.0;
   p = (2.0*a1*a1*a1-9.0*a1*a2+27.0*a3) / 54.0;
   t = q*q*q - p*p;
   if (t >= 0.0) {
      t = acos((double) p /(q * sqrt(q)));
      x = -2.0*sqrt(q)*cos(t / 3.0) - a1/3.0;
   } else {
      if (p > 0.0) {
         t = pow(sqrt(-t)+p, (double) 1.0 / 3.0);
         x = -(t + q / t) - a1/3.0;
      } else if (p == 0.0) {
         x = -a1/3.0;
      } else {
         t = pow(sqrt(-t)-p, (double) 1.0 / 3.0);
         x = (t + q / t) - a1/3.0;
      }
   }
   {
      double ex1;
      int i = 0;
      ex1 = x;
      for (t = root3(a1, a2, a3, x); ABS(t-x) > 5.0e-4;
           t = root3(a1, a2, a3, x))
	 if (++i == 32) {
            x = ex1;
            break;
         } else
            x = t;
   }
   /***
   x = a1;
   for (t = root3(a1, a2, a3, x); ABS(t-x) > epsi2; 
			    t = root3(a1, a2, a3, x)) {
      x = t;
      i++;
      if (i > 1000) {
         x = 0.5 * (x + root3(a1, a2, a3, x));
         j++;
         if (j == 3)
            break;
         i = 0;
      }
   }
    ***/
   *ex1 = x;
   div3(a1, a2, a3, x, &a1, &a2);

   t = a1 * a1 - 4.0 * a2;
   if (t < 0) {
      ifImg = 1;
      printf("***** Two Imaginary Roots.\n");
      *ex3 = 0.5 * sqrt(-t);
      *ex2 = -0.5 * a1;
   } else {
      ifImg = 0;
      t *= 1.0e-16;
      t = sqrt(t)*1.0e8;
      if (a1 >= 0.0)
         *ex2 = t = -0.5 * (a1 + t);
      else
         *ex2 = t = -0.5 * (a1 - t);
      *ex3 = a2 / t;
      /*
      *ex2 = 0.5 * (-a1 + t);
      *ex3 = 0.5 * (-a1 - t);
       */
   }

   return(1);
}
static NDnamePt
insert_ND(char *name, NDnamePt *ndn)
{
   int       cmp;
   NDnamePt  p;

   if (*ndn == NULL) {
      p = *ndn = (NDnamePt) tmalloc(sizeof (NDname));
      p->nd = NULL;
      p->right = p->left = NULL;
      strcpy(p->id, name);
      return(p);
   }
   cmp = strcmp((*ndn)->id, name);
   if (cmp == 0) 
      return(*ndn);
   else {
      if (cmp < 0)
         return(insert_ND(name, &((*ndn)->left)));
      else
         return(insert_ND(name, &((*ndn)->right)));
   }
}


static NODE 
*insert_node(char *name)
{
   NDnamePt n;
   NODE    *p;

   n = insert_ND(name, &ndn);
   if (n->nd == NULL) {
      p = NEW_node();
      p->name = n;
      n->nd = p;
      p->next = node_tab;
      node_tab = p;
      return(p);
   } else
      return(n->nd);
}

static NODE
*NEW_node(void)
{
   NODE *n;

   n = (NODE *) tmalloc (sizeof (NODE));
   n->mptr = NULL;
   n->gptr = NULL;
   n->cptr = NULL;
   n->rptr = NULL;
   n->tptr = NULL;
   n->cplptr = NULL;
   n->rlptr = NULL;
   n->ddptr = NULL;
   n->cvccsptr = NULL;
   n->vccsptr = NULL;
   n->CL = 0.0;
   n->V = n->dv = 0.0;
   n->gsum = n->cgsum = 0;
   n->is   = 0;
   n->tag  = 0;
   n->flag = 0;
   n->region = NULL;
   n->ofile = NULL;
   n->dvtag = 0;

   return(n);
}

static int 
find_roots(double a1, double a2, double a3, double *x1, double *x2, double *x3)
{
   double x, t;
   double p, q;

   q = (a1*a1-3.0*a2) / 9.0;
   p = (2.0*a1*a1*a1-9.0*a1*a2+27.0*a3) / 54.0;
   t = q*q*q - p*p;
   if (t >= 0.0) {
      t = acos((double) p /(q * sqrt(q)));
      x = -2.0*sqrt(q)*cos(t / 3.0) - a1/3.0;
   } else {
      if (p > 0.0) {
	 t = pow(sqrt(-t)+p, (double) 1.0 / 3.0);
         x = -(t + q / t) - a1/3.0;
      } else if (p == 0.0) {
	 x = -a1/3.0;
      } else {
 	 t = pow(sqrt(-t)-p, (double) 1.0 / 3.0);
	 x = (t + q / t) - a1/3.0;
      }
   }
   {
      double x1;
      int i = 0;
      x1 = x;
      for (t = root3(a1, a2, a3, x); ABS(t-x) > 5.0e-4;
			t = root3(a1, a2, a3, x))
         if (++i == 32) {
	    x = x1;
	    break;
	 } else
	    x = t;
   }
   /*
   x = a1;
   i = 0;
   j = 0;
   for (t = root3(a1, a2, a3, x); ABS(t-x) > epsi; 
			    t = root3(a1, a2, a3, x)) {
      x = t;
      i++;
      if (i > 1000) {
	 x = 0.5 * (x + root3(a1, a2, a3, x));
	 j++;
	 if (j == 3)
	    break;
	 i = 0;
      }
   }
    */

   *x1 = x;
   div3(a1, a2, a3, x, &a1, &a2);

   t = a1 * a1 - 4.0 * a2;
   if (t < 0) {
      printf("***** Two Imaginary Roots in Characteristic Admittance.\n");
      exit(0);
   }

   t *= 1.0e-18;
   t = sqrt(t) * 1.0e9;
   if (a1 >= 0.0)
      *x2 = t = -0.5 * (a1 + t);
   else
      *x2 = t = -0.5 * (a1 - t);
   *x3 = a2 / t;
   /*
   *x2 = 0.5 * (-a1 + t);
   *x3 = 0.5 * (-a1 - t);
    */
   return(1);
}