cplsetup.c

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

/****************************************************************
     match     Create a polynomial matching given data points 
 ****************************************************************/


static double 
*vector(int nl, int nh)
{
   double *v;

   v = (double *) tmalloc((unsigned) (nh - nl +1) * sizeof(double));
   if (!v) {
      fprintf(stderr, "Memory Allocation Error by tmalloc in vector().\n");
      fprintf(stderr, "...now exiting to system ...\n");
      exit(0);
   }
   return v-nl;
}

static void 
free_vector(double *v, int nl, int nh)
{
   free((void*) (v +nl));
}

static void 
polint(double *xa, double *ya, int n, double x, double *y, double *dy)
/*
   Given arrays xa[1..n] and ya[1..n], and given a value x, this routine
   returns a value y, and an error estimate dy.  If P(x) is the 
   polynomial of degree n-1 such that P(xa) = ya, then the returned 
   value y = P(x)
 */
{
   int i, m, ns = 1;
   double den, dif, dift, ho, hp, w;
   double *c, *d;

   dif = ABS(x - xa[1]);
   c = vector(1, n);
   d = vector(1, n);
   for (i = 1; i <= n; i++) {
      if ((dift = ABS(x - xa[i])) < dif) {
	 ns = i;
	 dif = dift;
      }
      c[i] = ya[i];
      d[i] = ya[i];
   }
   *y = ya[ns--];
   for (m = 1; m < n; m++) {
      for (i = 1; i <= n-m; i++) {
	 ho = xa[i]-x;
	 hp = xa[i+m]-x;
	 w = c[i+1]-d[i];
	 if ((den=ho-hp) == 0.0) {
	    fprintf(stderr, "(Error) in routine POLINT\n");
            fprintf(stderr, "...now exiting to system ...\n");
            exit(0);
	 }
	 den = w/den;
	 d[i] = hp * den;
	 c[i] = ho * den;
      }
      *y += (*dy = (2*ns < (n-m) ? c[ns+1] : d[ns--]));
   }
   free_vector(d, 1, n);
   free_vector(c, 1, n);
}

static int 
match(int n, double *cof, double *xa, double *ya)
/*
   Given arrays xa[0..n] and ya[0..n] containing a tabulated function
   ya = f(xa), this routine returns an array of coefficients cof[0..n],
   such that ya[i] = sum_j {cof[j]*xa[i]**j}.
 */
{
   int k, j, i;
   double xmin, dy, *x, *y, *xx;

   x = vector(0, n);
   y = vector(0, n);
   xx = vector(0, n);
   for (j = 0; j <= n; j++) {
      x[j] = xa[j];
      xx[j] = y[j] = ya[j];
   }
   for (j = 0; j <= n; j++) {
      polint(x-1, y-1, n+1-j, (double) 0.0, &cof[j], &dy);
      xmin = 1.0e38;
      k = -1;
      for (i = 0; i <= n-j; i++) {
	 if (ABS(x[i]) < xmin) {
	    xmin = ABS(x[i]);
	    k = i;
	 }
	 if (x[i]) y[i] = (y[i] - cof[j]) / x[i];
      }
      for (i = k+1; i <= n-j; i++) {
	 y[i-1] = y[i];
	 x[i-1] = x[i];
      }
   }
   free_vector(y, 0, n);
   free_vector(x, 0, n);

   /****   check   ****/
   /*
   for (i = 0; i <= n; i++) {
      xmin = xa[i];
      dy = cof[0];
      for (j = 1; j <= n; j++) {
	 dy += xmin * cof[j];
         xmin *= xa[i];
      }
      printf("*** real x = %e y = %e\n", xa[i], xx[i]);
      printf("*** calculated  y = %e\n", dy);
      printf("*** error = %e \% \n", (dy-xx[i])/xx[i]);
   }
   */
   return 0;
}

/***
 ***/
/***
static int
match_x(int dim, double *Alfa, double *X, double *Y)
{
   int i, j;
   double f;
   double scale;

   ****   check   ****
   double xx[16];
   for (i = 0; i <= dim; i++)
      xx[i] = Y[i];

   if (Y[1] == Y[0])
      scale = 1.0;
   else
      scale = X[1] / (Y[1] - Y[0]);
   for (i = 0; i < dim; i++) {
      f = X[i+1];
      for (j = dim-1; j >= 0; j--) {
         A[i][j] = f;
	 f *= X[i+1];
      }
      A[i][dim] = (Y[i+1] - Y[0])*scale;
   }
   Gaussian_Elimination2(dim, 1);
   Alfa[0] = Y[0];
   for (i = 1; i <= dim; i++) 
      Alfa[i] = A[dim-i][dim] / scale;

   ****   check   ****
   * 
   for (i = 0; i <= dim; i++) {
      f = X[i];
      scale = Alfa[0];
      for (j = 1; j <= dim; j++) {
	 scale += f * Alfa[j];
         f *= X[i];
      }
      printf("*** real x = %e y = %e\n", X[i], xx[i]);
      printf("*** calculated  y = %e\n", scale);
      printf("*** error = %e \% \n", (scale-xx[i])/xx[i]);
   }
   *

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

static int 
Gaussian_Elimination2(int dims, int type)
   /*  type = 1 : to solve a linear system  
	     -1 : to inverse a matrix  */
{
   int i, j, k, dim;
   double f;
   double max;
   int imax;

   if (type == -1)
      dim = 2 * dims;
   else
      dim = dims;

   for (i = 0; i < dims; 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 < epsilon) {
         fprintf(stderr, " can not choose a pivot (misc)\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 < dims ; 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 void
eval_Si_Si_1(int dim, double y)
{
   int i, j, k;

   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) {
	 Si_1[i][j] = 0.0;
	 for (k = 0; k < dim; k++)
	    if (k == j)
			Si_1[i][j] += Sv_1[i][k] *
				(y * R_m[k][j] + Scaling_F * L_m[k][j]);
	    else 
			Si_1[i][j] += Sv_1[i][k] * L_m[k][j] * Scaling_F;
			/
			Si_1[i][j] *= Scaling_F;
			/
      }

   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++)
	 Si_1[i][j] /= sqrt((double) D[i]);

   for (i = 0; i < dim; i++) {
      for (j = 0; j < dim; j++)
	 A[i][j] = Si_1[i][j];
      for (j = dim; j < 2* dim; j++)
	 A[i][j] = 0.0;
      A[i][i+dim] = 1.0;
   }
   Gaussian_Elimination2(dim, -1);

   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++)
	 Si[i][j] = A[i][j+dim];
}

***/

static void
eval_Si_Si_1(int dim, double y)
{
   int i, j, k;

   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) {
	 Si_1[i][j] = 0.0;
	 for (k = 0; k < dim; k++)
	    Si_1[i][j] += Sv_1[i][k] * (y * R_m[k][j] + Scaling_F * L_m[k][j]);
	       /*
	    else 
	       Si_1[i][j] += Sv_1[i][k] * L_m[k][j] * Scaling_F;
	 Si_1[i][j] *= Scaling_F;
		*/
      }

   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++)
	 Si_1[i][j] /= sqrt((double) D[i]);

   for (i = 0; i < dim; i++) {
      for (j = 0; j < dim; j++)
	 A[i][j] = Si_1[i][j];
      for (j = dim; j < 2* dim; j++)
	 A[i][j] = 0.0;
      A[i][i+dim] = 1.0;
   }
   Gaussian_Elimination2(dim, -1);

   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++)
	 Si[i][j] = A[i][j+dim];
}

/***

static void
loop_ZY(int dim, double y)
{
   int i, j, k;
   double fmin, fmin1;

   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) 
	 if (i == j) 
	    ZY[i][j] = Scaling_F * C_m[i][i] + G_m[i] * y;
	 else
	    ZY[i][j] = Scaling_F * C_m[i][j];
   diag(dim);
   fmin = D[0];
   for (i = 1; i < dim; i++)
      if (D[i] < fmin)
	 fmin = D[i];
   if (fmin < 0) {
      fprintf(stderr, "(Error) The capacitance matrix of the multiconductor system is not positive definite.\n");
      exit(0);
   } else {
      fmin = sqrt(fmin);
      fmin1 = 1 / fmin;
   }
   for (i = 0; i < dim; i++)
      D[i] = sqrt((double) D[i]);
   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) {
	 Y5[i][j] = D[i] * Sv[j][i];
	 Y5_1[i][j] = Sv[j][i] / D[i];
      }
   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) {
	 Sv_1[i][j] = 0.0;
	 for (k = 0; k < dim; k++) 
	    Sv_1[i][j] += Sv[i][k] * Y5[k][j];
      }
   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) 
	 Y5[i][j] = Sv_1[i][j]; 
   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) {
	 Sv_1[i][j] = 0.0;
	 for (k = 0; k < dim; k++) 
	    Sv_1[i][j] += Sv[i][k] * Y5_1[k][j];
      }
   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) 
	 Y5_1[i][j] = Sv_1[i][j]; 

   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) {
	 ZY[i][j] = 0.0;
	 for (k = 0; k < dim; k++)
            if (k == i) 
	       ZY[i][j] += (Scaling_F *  L_m[i][i] + R_m[i] * y) * 
			       Y5[k][j];
	    else 
	       ZY[i][j] += L_m[i][k] * Y5[k][j] * Scaling_F;
      }
   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) {
	 Sv_1[i][j] = 0.0;
	 for (k = 0; k < dim; k++) 
	    Sv_1[i][j] += Y5[i][k] * ZY[k][j];
      }
   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) 
	 ZY[i][j] = Sv_1[i][j]; 

   diag(dim);

   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) {
	 Sv_1[i][j] = 0.0; 
	 for (k = 0; k < dim; k++) 
	    Sv_1[i][j] += Sv[k][i] * Y5[k][j];
	 Sv_1[i][j] *= fmin1;
      }
   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) {
	 ZY[i][j] = 0.0; 
	 for (k = 0; k < dim; k++) 
	    ZY[i][j] += Y5_1[i][k] * Sv[k][j];
	 ZY[i][j] *= fmin;
      }
   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) 
	 Sv[i][j] = ZY[i][j]; 

}
***/

static void
loop_ZY(int dim, double y)
{
   int i, j, k;
   double fmin, fmin1;

   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) 
	 ZY[i][j] = Scaling_F * C_m[i][j] + G_m[i][j] * y;
	    /*
	 else
	    ZY[i][j] = Scaling_F * C_m[i][j];
	     */
   diag(dim);
   fmin = D[0];
   for (i = 1; i < dim; i++)
      if (D[i] < fmin)
	 fmin = D[i];
   if (fmin < 0) {
      fprintf(stderr, "(Error) The capacitance matrix of the multiconductor system is not positive definite.\n");
      exit(0);
   } else {
      fmin = sqrt(fmin);
      fmin1 = 1 / fmin;
   }
   for (i = 0; i < dim; i++)
      D[i] = sqrt((double) D[i]);
   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) {
	 Y5[i][j] = D[i] * Sv[j][i];
	 Y5_1[i][j] = Sv[j][i] / D[i];
      }
   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) {
	 Sv_1[i][j] = 0.0;
	 for (k = 0; k < dim; k++) 
	    Sv_1[i][j] += Sv[i][k] * Y5[k][j];
      }
   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) 
	 Y5[i][j] = Sv_1[i][j]; 
   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) {
	 Sv_1[i][j] = 0.0;
	 for (k = 0; k < dim; k++) 
	    Sv_1[i][j] += Sv[i][k] * Y5_1[k][j];
      }
   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) 
	 Y5_1[i][j] = Sv_1[i][j]; 

   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) {
	 ZY[i][j] = 0.0;
	 for (k = 0; k < dim; k++)
	    ZY[i][j] += (Scaling_F *  L_m[i][k] + R_m[i][k] * y) * Y5[k][j];
	 /*
	    else 
	       ZY[i][j] += L_m[i][k] * Y5[k][j] * Scaling_F;
	  */
      }
   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) {
	 Sv_1[i][j] = 0.0;
	 for (k = 0; k < dim; k++) 
	    Sv_1[i][j] += Y5[i][k] * ZY[k][j];
      }
   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) 
	 ZY[i][j] = Sv_1[i][j]; 

   diag(dim);

   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) {
	 Sv_1[i][j] = 0.0; 
	 for (k = 0; k < dim; k++) 
	    Sv_1[i][j] += Sv[k][i] * Y5[k][j];
	 Sv_1[i][j] *= fmin1;
      }
   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) {
	 ZY[i][j] = 0.0; 
	 for (k = 0; k < dim; k++) 
	    ZY[i][j] += Y5_1[i][k] * Sv[k][j];
	 ZY[i][j] *= fmin;
      }
   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) 
	 Sv[i][j] = ZY[i][j]; 

}


/***
 ***/

static void 
poly_matrix(A,  dim,  deg)
   double *A[MAX_DIM][MAX_DIM];
   int dim, deg;
{
   int i, j;

   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++)
	 match(deg, A[i][j], frequency, A[i][j]);
}

/***
 ***/
/***
static int 
checkW(double *W, double d)
{
  double f, y;
  float  y1;
  int k;

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

  f = W[0];
  y = y1;
  f += y * W[1];
  for (k = 2; k < 6; k++) {
     y *= y1;
     f += y * W[k];