cplsetup.c

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

{
   double p1, p2, p3, q1, q2, q3;

   At[0][0] = 1.0 - a_b;
   At[0][1] = b[1];
   At[0][2] = b[2];
   At[0][3] = -b[3];

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

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

   Gaussian_Elimination(3);

   p3 = At[0][3];
   p2 = At[1][3];
   p1 = At[2][3];
   /*
   if (p3 < 0.0 || p2 < 0.0 || p1 < 0.0 || p1*p2 <= p3)
      return(0);
    */
   q1 = p1 + b[1];
   q2 = b[1] * p1 + p2 + b[2];
   q3 = p3 * a_b;

   if (find_roots(p1, p2, p3, x1, x2, x3)) {
      /*
      printf("complex roots : %e %e %e \n", *x1, *x2, *x3);
       */
      *c1 = eval2(q1 - p1, q2 - p2, q3 - p3, *x1) /
                   eval2((double) 3.0, (double) 2.0 * p1, p2, *x1);
      get_c(q1 - p1, q2 - p2, q3 - p3, p1, p2, *x2, *x3, c2, c3);
      return(2);
   } else {
      /* new
	  printf("roots are %e %e %e \n", *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);
      return(1);
   }
}

static int 
Gaussian_Elimination(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(At[i][i]);
      for (j = i+1; j < dim; j++)
	 if (ABS(At[j][i]) > max) {
	    imax = j;
	    max = ABS(At[j][i]);
         }
      if (max < epsi_mult) {
         fprintf(stderr, " can not choose a pivot (mult)\n");
         exit(0);
      }
      if (imax != i)
	 for (k = i; k <= dim; k++) {
	    f = At[i][k];
	    At[i][k] = At[imax][k];
	    At[imax][k] = f;
	 }
      
      f = 1.0 / At[i][i];
      At[i][i] = 1.0;

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

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

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

   t1 = x * (x * (x + a1) + a2) + a3;
   t2 = x * (2.0*a1 + 3.0*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 = a2 + (a1 + x) * x; */

	*p2 = - a3 / x;

   return(1);
}


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

   /***********************************************
   double m,n;
   p = a1*a1/3.0 - a2;
   q = a1*a2/3.0 - a3 - 2.0*a1*a1*a1/27.0;
   p = p*p*p/27.0;
   t = q*q - 4.0*p;
   if (t < 0.0) {
      if (q != 0.0) {
	 t = atan(sqrt((double)-t)/q);
	 if (t < 0.0)
	    t += 3.141592654;
	 t /= 3.0;
	 x = 2.0 * pow(p, 0.16666667) * cos(t) - a1 / 3.0;
      } else {
         t /= -4.0;
         x = pow(t, 0.16666667) * 1.732 - a1 / 3.0;
      }
   } else {
      t = sqrt(t);
      m = 0.5*(q - t);
      n = 0.5*(q + t);
      if (m < 0.0)
         m = -pow((double) -m, (double) 0.3333333);
      else
         m = pow((double) m, (double) 0.3333333);
      if (n < 0.0)
         n = -pow((double) -n, (double) 0.3333333);
      else
         n = pow((double) n, (double) 0.3333333);
      x = m + n - a1 / 3.0;
   }
    ************************************************/
   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;
      }
   }
   /*
   fprintf(stderr, "..1.. %e\n", x*x*x+a1*x*x+a2*x+a3);
    */
   { 
      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;
   }
   /*
   fprintf(stderr, "..2.. %e\n", x*x*x+a1*x*x+a2*x+a3);
    */


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

   t = a1 * a1 - 4.0 * a2;
   if (t < 0) {
      /*
      fprintf(stderr, "***** Two Imaginary Roots.\n Update.\n");
      *x2 = -0.5 * a1;
      *x3 = a2 / *x2;
       */
      *x3 = 0.5 * sqrt(-t);
      *x2 = -0.5 * a1;
      return(1);
   } else {
      t = sqrt(t);
      if (a1 >= 0.0)
	 *x2 = t = -0.5 * (a1 + t);
      else
	 *x2 = t = -0.5 * (a1 - t);
      *x3 = a2 / t;
      return(0);
   }
}


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 int divC(double ar, double ai, double br, double bi, double *cr, double *ci)
{
	double t;

	t = br*br + bi*bi;
	*cr = (ar*br + ai*bi) / t;
	*ci = (ai*br - ar*bi) / t;

	return(1);
}
***/

static NODE
*NEW_node(void)
{
   /*
   char *tmalloc();
    */
   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.001;
   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);
}



/****************************************************************
     diag.c      This file contains the main().
 ****************************************************************/

#define  epsi2    1.0e-8

static  int         dim;
static  MAXE_PTR    row;

static MAXE_PTR 
sort(MAXE_PTR list, float val, int r, int c, MAXE_PTR e)
{
   if (list == NULL || list->value < val) {
      e->row = r;
      e->col = c;
      e->value = val;
      e->next = list;
      return(e);
   } else {
      list->next = sort(list->next, val, r, c, e);
      return(list);
   }
}


static void
ordering(void)
{
   MAXE_PTR e;
   int i, j, m;
   float mv;

   for (i = 0; i < dim-1; i++) {
      m = i+1;
      mv = ABS(ZY[i][m]);
      for (j = m+1; j < dim; j++) 
	if ((int)(ABS(ZY[i][j]) * 1e7) > (int) (1e7 *mv)) {

	    mv = ABS(ZY[i][j]);
	    m = j;
	 }
      e = (MAXE_PTR) tmalloc(sizeof (MAXE));
      row = sort(row, mv, i, m, e);
   }
}
      

static MAXE_PTR 
delete_1(MAXE_PTR *list, int rc)
{
   MAXE_PTR list1, e;

   list1 = *list;
   if ((*list)->row == rc) {
      *list = (*list)->next;
      return(list1);
   }
   for (e = list1->next; e->row != rc; e = e->next) 
      list1 = e;
   list1->next = e->next;
   return(e);
}


static void
reordering(int p, int q)
{
   MAXE_PTR e;
   int j, m;
   float mv;

   m = p+1;
   mv = ABS(ZY[p][m]);
   for (j = m+1; j < dim; j++) 
	if ((int)(ABS(ZY[p][j]) * 1e7) > (int) (1e7 *mv)) {
         mv = ABS(ZY[p][j]);
         m = j;
      }
   e = delete_1(&row, p);
   row = sort(row, mv, p, m, e);
   
   m = q+1;
   if (m != dim) {
      mv = ABS(ZY[q][m]);
      for (j = m+1; j < dim; j++)
	if ((int)(ABS(ZY[q][j]) * 1e7) > (int) (1e7 *mv)) {
	
            mv = ABS(ZY[q][j]);
            m = j;
         }
      e = delete_1(&row, q);
      row = sort(row, mv, q, m, e);
   }

}

static void 
diag(int dims)
{
   int i, j, c;
   double fmin, fmax;

   dim = dims;
   row = NULL;

   fmin = fmax = ABS(ZY[0][0]);
   for (i = 0; i < dim; i++) 
      for (j = i; j < dim; j++) 
	 if (ABS(ZY[i][j]) > fmax)
	    fmax = ABS(ZY[i][j]);
	 else if (ABS(ZY[i][j]) < fmin)
	    fmin = ABS(ZY[i][j]);
   fmin = 2.0 / (fmin + fmax);
   for (i = 0; i < dim; i++)
      for (j = i; j < dim; j++)
         ZY[i][j] *= fmin;

   for (i = 0; i < dim; i++) {
      for (j = 0; j < dim; j++) 
		if (i == j)
			Sv[i][i] = 1.0;
		else 
			Sv[i][j] = 0.0;
   }

   ordering();

   if (row)
	   for (c = 0; row->value > epsi2; c++) {
   		   int p, q;

   		   p = row->row;
   		   q = row->col;

   		   rotate(dim, p, q);
   		   reordering(p, q);
   		}

   for (i = 0; i < dim; i++)
      D[i] = ZY[i][i] / fmin;
}

/****************************************************************
     rotate()      rotation of the Jacobi's method
 ****************************************************************/

static int
rotate(int dim, int p, int q)
{
   int j;
   double co, si;
   double ve, mu, ld;
   double T[MAX_DIM];
   double t;

   ld = - ZY[p][q];
   mu = 0.5 * (ZY[p][p] - ZY[q][q]);
   ve = sqrt((double) ld*ld + mu*mu);
   co = sqrt((double) (ve + ABS(mu)) / (2.0 * ve));
   si = SGN(mu) * ld / (2.0 * ve * co);

   for (j = p+1; j < dim; j++)
      T[j] = ZY[p][j];
   for (j = 0; j < p; j++)
      T[j] = ZY[j][p];

   for (j = p+1; j < dim; j++) {
      if (j == q)
	 continue;
      if (j > q)
         ZY[p][j] = T[j] * co - ZY[q][j] * si;
      else
         ZY[p][j] = T[j] * co - ZY[j][q] * si;
   }
   for (j = q+1; j < dim; j++) {
      if (j == p)
	 continue;
      ZY[q][j] = T[j] * si + ZY[q][j] * co;
   }
   for (j = 0; j < p; j++) {
      if (j == q)
	 continue;
      ZY[j][p] = T[j] * co - ZY[j][q] * si;
   }
   for (j = 0; j < q; j++) {
      if (j == p)
	 continue;
      ZY[j][q] = T[j] * si + ZY[j][q] * co;
   }

   t = ZY[p][p];
   ZY[p][p] = t * co * co + ZY[q][q] * si * si - 2.0 * ZY[p][q] * si * co;
   ZY[q][q] = t * si * si + ZY[q][q] * co * co + 2.0 * ZY[p][q] * si * co;

   ZY[p][q] = 0.0;

   {
      double R[MAX_DIM];

      for (j = 0; j < dim; j++) {
         T[j] = Sv[j][p];
         R[j] = Sv[j][q];
      }

      for (j = 0; j < dim; j++) {
         Sv[j][p] = T[j] * co - R[j] * si;
	 Sv[j][q] = T[j] * si + R[j] * co;
      }
   }

   return(1);

}