cplsetup.c

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

  }
  printf("W[i]= %e\n ", f*exp((double)-d/y1));

  return(1);
}
***/
/***
 ***/

static void 
poly_W(int dim, int deg)
{
   int i;

   for (i = 0; i < dim; i++) {
      match(deg, W[i], frequency, W[i]);
      TAU[i] = approx_mode(W[i], W[i], length);
	  /*
      checkW(W[i], TAU[i]);
	  */
   }
}

/***
 ***/
   
static void
eval_frequency(int dim, int deg_o)
{
   int i, im;
   double min;

   min = D[0];
   im = 0;

   for (i = 1; i < dim; i++) 
      if (D[i] < min) {
	 min = D[i];
	 im = i;
      }

   if (min <= 0) {
      fprintf(stderr, "A mode frequency is not positive.  Abort!\n");
      exit(0);
   }

   Scaling_F2 = 1.0 / min;
   Scaling_F = sqrt(Scaling_F2);
   min = length * 8.0;
   /*
   min *= 1.0e18;
   min = sqrt(min)*1.0e-9*length/8.0;
    */

   frequency[0] = 0.0;

   for (i = 1; i <= deg_o; i++)
      frequency[i] = frequency[i-1] + min;

   for (i = 0; i < dim; i++)
      D[i] *= Scaling_F2;
}

/***
 ***/

static void
store(int dim, int ind)
{
   int i, j;

   for (i = 0; i < dim; i++) {
      for (j = 0; j < dim; j++) {
         /*  store_Sip  */
	 Sip[i][j][ind] = Si[i][j];
         /*  store_Si_1p  */
	 Si_1p[i][j][ind] = Si_1[i][j];
         /*  store_Sv_1p  */
	 Sv_1p[i][j][ind] = Sv_1[i][j];
      }
      /*  store_W  */
      W[i][ind] = D[i];
   }
}

/***
 ***/

static void
store_SiSv_1(int dim, int ind)
{
   int i, j, k;
   double temp;

   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) {
	 temp = 0.0;
	 for (k = 0; k < dim; k++)
	   temp += Si[i][k] * Sv_1[k][j];
      SiSv_1[i][j][ind] = temp;
   }
}

/***
 ***/
/***
static int 
check(Sip, Si_1p, Sv_1p, SiSv_1p)
   double *Sip[MAX_DIM][MAX_DIM];
   double *Si_1p[MAX_DIM][MAX_DIM];
   double *Sv_1p[MAX_DIM][MAX_DIM];
   double *SiSv_1p[MAX_DIM][MAX_DIM];
{
   double f, y;
   float  y1;
   int i, j, k;

   printf("y =");
   scanf("%f", &y1);

   printf("\n");
   printf("Si =\n");
   for (i = 0; i < 4; i++) {
      for (j = 0; j < 4; j++) {
	 f = Sip[i][j][0];
         y = y1;
	 f += y * Sip[i][j][1];
	 for (k = 2; k < 8; k++) {
	    y *= y1;
	    f += y * Sip[i][j][k];
	 }
	 printf("%e ", f);
      }
      printf("\n");
   }
   printf("\n");
   printf("Si_1 =\n");
   for (i = 0; i < 4; i++) {
      for (j = 0; j < 4; j++) {
	 f = Si_1p[i][j][0];
         y = y1;
	 f += y * Si_1p[i][j][1];
	 for (k = 2; k < 8; k++) {
	    y *= y1;
	    f += y * Si_1p[i][j][k];
	 }
	 printf("%e ", f);
      }
      printf("\n");
   }
   printf("\n");
   printf("Sv_1 =\n");
   for (i = 0; i < 4; i++) {
      for (j = 0; j < 4; j++) {
	 f = Sv_1p[i][j][0];
         y = y1;
	 f += y * Sv_1p[i][j][1];
	 for (k = 2; k < 8; k++) {
	    y *= y1;
	    f += y * Sv_1p[i][j][k];
	 }
	 printf("%e ", f);
      }
      printf("\n");
   }
   printf("\n");
   printf("SiSv_1 =\n");
   for (i = 0; i < 4; i++) {
      for (j = 0; j < 4; j++) {
	 f = SiSv_1p[i][j][0];
         y = y1;
	 f += y * SiSv_1p[i][j][1];
	 for (k = 2; k < 8; k++) {
	    y *= y1;
	    f += y * SiSv_1p[i][j][k];
	 }
	 printf("%e ", f);
      }
      printf("\n");
   }
   return(1);
}
***/
/***
 ***/

static int 
coupled(int dim)
{
   int deg, deg_o;
   int i;

   deg = Right_deg;
   deg_o =  Left_deg;
   new_memory(dim, deg, deg_o);

   Scaling_F = Scaling_F2 = 1.0;

   /***     y = 0 : ZY = LC    ***/
   loop_ZY(dim, (double) 0.0);
   eval_frequency(dim, deg_o);
   eval_Si_Si_1(dim, (double) 0.0);
   store_SiSv_1(dim, (int) 0);
   store(dim, (int) 0);

   /***     Step  1     ***/
   /***     Step  2     ***/
   for (i = 1; i <= deg_o; i++) { 
       loop_ZY(dim, frequency[i]);
       eval_Si_Si_1(dim, frequency[i]);
       store_SiSv_1(dim, i);
       store(dim, i);
   }
   poly_matrix(Sip, dim, deg_o);
   poly_matrix(Si_1p, dim, deg_o);
   poly_matrix(Sv_1p, dim, deg_o);
   poly_W(dim, deg_o);
   matrix_p_mult(Sip, W, Si_1p, dim, deg_o, deg_o, IWI);
   matrix_p_mult(Sip, W, Sv_1p, dim, deg_o, deg_o, IWV);

   poly_matrix(SiSv_1, dim, deg_o);

   /***
   check(Sip, Si_1p, Sv_1p, SiSv_1);
   ***/

   generate_out(dim, deg_o);

   return(1);
}

/***
 ***/

static int 
generate_out(int dim, int deg_o)
{
   int i, j, k, rtv;
   double C;
   double *p;
   double c1, c2, c3, x1, x2, x3;

   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) {
	 p = SiSv_1[i][j];
	 SIV[i][j].C_0 = C = p[0];
	 if (C == 0.0)
	    continue;
	 for (k = 0; k <= deg_o; k++) 
	    p[k] /= C;
	 if (i == j) {
	    rtv = Pade_apx((double) sqrt((double) G_m[i][i] / R_m[i][i]) / C, 
		  p, &c1, &c2, &c3, &x1, &x2, &x3);
	    if (rtv == 0) {
	       SIV[i][j].C_0 = 0.0;
	       printf("SIV\n");
	       continue;
	    }
	 } else {
	    rtv = Pade_apx((double) 0.0, 
		  p, &c1, &c2, &c3, &x1, &x2, &x3);
	    if (rtv == 0) {
	       SIV[i][j].C_0 = 0.0;
	       printf("SIV\n");
	       continue;
	    }
	 }
         p = SIV[i][j].Poly = (double *) calloc(7, sizeof(double));
	 p[0] = c1;
	 p[1] = c2;
	 p[2] = c3;
	 p[3] = x1;
	 p[4] = x2;
	 p[5] = x3;
	 p[6] = (double) rtv;
      }
   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) 
	 for (k = 0; k < dim; k++) {
	    p = IWI[i][j].Poly[k];
	    C = IWI[i][j].C_0[k];
	    if (C == 0.0)
	       continue;
	    if (i == j && k == i) {
	       rtv = Pade_apx((double) 
	             exp(- sqrt((double) G_m[i][i] * R_m[i][i]) * length) / C,
		     p, &c1, &c2, &c3, &x1, &x2, &x3);
	       if (rtv == 0) {
		  IWI[i][j].C_0[k] = 0.0;
		  printf("IWI %d %d %d\n", i, j, k);
		  continue;
	       }
	    } else {
	       rtv = Pade_apx((double) 0.0, 
		  p, &c1, &c2, &c3, &x1, &x2, &x3);
	       if (rtv == 0) {
		  IWI[i][j].C_0[k] = 0.0;
		  printf("IWI %d %d %d\n", i, j, k);
		  continue;
	       }
	    }
	    p[0] = c1;
	    p[1] = c2;
	    p[2] = c3;
	    p[3] = x1;
	    p[4] = x2;
	    p[5] = x3;
	    p[6] = (double) rtv;
         }

   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) 
	 for (k = 0; k < dim; k++) {
	    p = IWV[i][j].Poly[k];
	    C = IWV[i][j].C_0[k];
	    if (C == 0.0)
	       continue;
	    if (i == j && k == i) {
	       rtv = Pade_apx((double) sqrt((double) G_m[i][i] / R_m[i][i]) * 
		     exp(- sqrt((double) G_m[i][i] * R_m[i][i]) * length) / C,
		     p, &c1, &c2, &c3, &x1, &x2, &x3);
	       if (rtv == 0) {
		  IWV[i][j].C_0[k] = 0.0;
		  printf("IWV %d %d %d\n", i, j, k);
		  continue;
	       }
	    } else {
	       rtv = Pade_apx((double) 0.0, 
		  p, &c1, &c2, &c3, &x1, &x2, &x3);
	       if (rtv == 0) {
		  IWV[i][j].C_0[k] = 0.0;
		  printf("IWV %d %d %d\n", i, j, k);
		  continue;
	       }
            }
	    p[0] = c1;
	    p[1] = c2;
	    p[2] = c3;
	    p[3] = x1;
	    p[4] = x2;
	    p[5] = x3;
	    p[6] = (double) rtv;
         }
   return(1);
}

/****************************************************************
     mult.c     Multiplication for Matrix of Polynomial 
		   X(y) = A(y) D(y) B(y),
		   where D(y) is a diagonal matrix with each 
		   diagonal entry of the form
			   e^{-a_i s}d(y), for which s = 1/y
				                 and i = 1..N.
		   Each entry of X(y) will be of the form
		      \sum_{i=1}^N c_i e^{-a_i s} b_i(y), where
		      b_i(0) = 1; therefore, those
		      b_i(y)'s will be each entry's output. 
 ****************************************************************/

static void
mult_p(double *p1, double *p2, double *p3, int d1, int d2, int d3)
/*   p3 = p1 * p2   */
{
   int i, j, k;

   for (i = 0; i <= d3; i++)
      p3[i] = 0.0;
   for (i = 0; i <= d1; i++)
      for (j = i, k = 0; k <= d2; j++, k++) {
	 if (j > d3)
	    break;
	 p3[j] += p1[i] * p2[k];
      }
}


static void
matrix_p_mult(A, D, B, dim, deg, deg_o, X)
   double  *A[MAX_DIM][MAX_DIM]; 
   double  *D[MAX_DIM]; 
   double  *B[MAX_DIM][MAX_DIM]; 
   int     dim, deg, deg_o;
   Mult_Out  X[MAX_DIM][MAX_DIM]; 
{
   int i, j, k, l;
   double *p;
   double *T[MAX_DIM][MAX_DIM];
   double t1;

   for (i  = 0; i < dim; i++) 
      for (j = 0; j < dim; j++) {
	 p = T[i][j] = (double *) calloc(deg_o+1, sizeof(double));
	 mult_p(B[i][j], D[i], p, deg, deg_o, deg_o);
      }
   for (i  = 0; i < dim; i++)
      for (j = 0; j < dim; j++) 
	 for (k = 0; k < dim; k++) {
	    p = X[i][j].Poly[k] = 
			    (double *) calloc(deg_o+1, sizeof(double));
	    mult_p(A[i][k], T[k][j], p, deg, deg_o, deg_o);
	    t1 = X[i][j].C_0[k] = p[0];
	    if (t1 != 0.0) {
	       p[0] = 1.0;
	       for (l = 1; l <= deg_o; l++)
	          p[l] /= t1;
	    }
	 }
      
      /**********
   for (i  = 0; i < dim; i++)
      for (j = 0; j < dim; j++) {
	    for (k = 0; k < dim; k++) {
	    fprintf(outFile, "(%5.3f)", X[i][j].C_0[k]);
	    p = X[i][j].Poly[k];
	    for (l = 0; l <= deg_o; l++)
	       fprintf(outFile, "%5.3f ", p[l]); 
	    fprintf(outFile, "\n");
	 }
         fprintf(outFile, "\n");
      }
       ***********/
}

/****************************************************************
	          mode  approximation  

 ****************************************************************/

/***
 ***/

static double 
approx_mode(double *X, double *b, double length)
{
   double w0, w1, w2, w3, w4, w5;
   double a[8];
   double delay, atten;
   double y1, y2, y3, y4, y5, y6;
   int    i, j;

   w0 = X[0];
   w1 = X[1] / w0;  /* a */
   w2 = X[2] / w0;  /* b */
   w3 = X[3] / w0;  /* c */
   w4 = X[4] / w0;  /* d */
   w5 = X[5] / w0;  /* e */

   y1 = 0.5 * w1;
   y2 = w2 - y1 * y1;
   y3 = 3 * w3 - 3.0 * y1 * y2;
   y4 = 12.0 * w4 - 3.0 * y2 * y2 - 4.0 * y1 * y3;
   y5 = 60.0 * w5 - 5.0 * y1 * y4 -10.0 * y2 * y3;
   y6 = -10.0 * y3 * y3 - 15.0 * y2 * y4 - 6.0 * y1 * y5;

   delay = sqrt(w0) * length / Scaling_F;
   atten = exp((double) - delay * y1);

   a[1] = y2 / 2.0;
   a[2] = y3 / 6.0;
   a[3] = y4 / 24.0;
   a[4] = y5 / 120.0;
   a[5] = y6 / 720.0;

   a[1] *= -delay;
   a[2] *= -delay;
   a[3] *= -delay;
   a[4] *= -delay;
   a[5] *= -delay;
   
   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;
   }

   for (i = 0; i <= 5; i++)
      b[i] *= atten;

   return(delay);
}

/***
 ***/

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

/***
 ***/

static int 
get_c(double q1, double q2, double q3, double p1, double p2, double a, double b, 
      double *cr, double *ci)
{
   double d, n;

   d = (3.0*(a*a-b*b)+2.0*p1*a+p2)*(3.0*(a*a-b*b)+2.0*p1*a+p2);
   d += (6.0*a*b+2.0*p1*b)*(6.0*a*b+2.0*p1*b);
   n = -(q1*(a*a-b*b)+q2*a+q3)*(6.0*a*b+2.0*p1*b);
   n += (2.0*q1*a*b+q2*b)*(3.0*(a*a-b*b)+2.0*p1*a+p2);
   *ci = n/d;
   n = (3.0*(a*a-b*b)+2.0*p1*a+p2)*(q1*(a*a-b*b)+q2*a+q3);
   n += (6.0*a*b+2.0*p1*b)*(2.0*q1*a*b+q2*b);
   *cr = n/d;

   return(1);
}


static int 
Pade_apx(double a_b, double *b, double *c1, double *c2, double *c3, 
         double *x1, double *x2, double *x3)
/*     
        b[0] + b[1]*y + b[2]*y^2 + ... + b[5]*y^5 + ...
      = (q3*y^3 + q2*y^2 + q1*y + 1) / (p3*y^3 + p2*y^2 + p1*y + 1)
       
	where b[0] is always equal to 1.0 and neglected,
	  and y = 1/s.

	(q3*y^3 + q2*y^2 + q1*y + 1) / (p3*y^3 + p2*y^2 + p1*y + 1)
      = (s^3 + q1*s^2 + q2*s + q3) / (s^3 + p1*s^2 + p2*s + p3)
      = c1 / (s - x1) + c2 / (s - x2) + c3 / (s - x3) + 1.0
 */