cplsetup.c

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

/**********
Copyright 1992 Regents of the University of California.  All rights
reserved.
Author: 1992 Charles Hough
Modified: 2004 Paolo Nenzi - (ng)spice integration
**********/


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

#include "multi_line.h"


#define VECTOR_ALLOC(vec, n) { \
	vec = (double **) tmalloc(n * sizeof(double *)); \
}

#define MATRIX_ALLOC(mat, m, j) { \
	int k; \
	mat = (double ***) tmalloc(m * sizeof(double **)); \
	for (k = 0; k < m; k++) {  \
		VECTOR_ALLOC(mat[k], j); \
	} \
}

#define VECTOR_FREE(vec) free(vec) 

#define MATRIX_FREE(mat, m, j) { \
	int k; \
	for (k = 0; k < m; k++) {  \
		free(mat[k]); \
	} \
        free(mat); \
}

#define MAX_DEG 8
#define epsilon 1.0e-88
#define MAX_STRING 128

static  double  ZY[MAX_DIM][MAX_DIM]; 
static  double  Sv[MAX_DIM][MAX_DIM]; 
static  double  D[MAX_DIM]; 
static  double  Y5[MAX_DIM][MAX_DIM];
static  double  Y5_1[MAX_DIM][MAX_DIM];
static  double  Sv_1[MAX_DIM][MAX_DIM];

static  double R_m[MAX_DIM][MAX_DIM];
static  double G_m[MAX_DIM][MAX_DIM];
static  double L_m[MAX_DIM][MAX_DIM];
static  double C_m[MAX_DIM][MAX_DIM];
static  double length;
static  double TAU[MAX_DIM];

static  double  A[MAX_DIM][2*MAX_DIM];

static  double  frequency[MAX_DEG];

static  double Si[MAX_DIM][MAX_DIM];
static  double Si_1[MAX_DIM][MAX_DIM];

/*  MacLaurin Series  */
static  double *SiSv_1[MAX_DIM][MAX_DIM];
static  double *Sip[MAX_DIM][MAX_DIM];
static  double *Si_1p[MAX_DIM][MAX_DIM];
static  double *Sv_1p[MAX_DIM][MAX_DIM];
static  double *W[MAX_DIM];

static  Mult_Out IWI[MAX_DIM][MAX_DIM];
static  Mult_Out IWV[MAX_DIM][MAX_DIM];
static  Single_Out SIV[MAX_DIM][MAX_DIM];
static  double  At[4][4];
static  double  Scaling_F;
static  double  Scaling_F2;

/* misc.c match */
static void new_memory(int, int, int);
static double *vector(int, int);
static void free_vector(double*, int, int);
static void polint(double*, double*, int, double, double*, double*);
static int match(int, double*, double*, double*);
/* static int match_x(int, double*, double*, double*); */
static int Gaussian_Elimination2(int, int);
static void eval_Si_Si_1(int, double);
static void loop_ZY(int, double);
static void poly_matrix(); /* quale è il prototipo ? */
/* static int checkW(double*, double); */
static void poly_W(int, int);
static void eval_frequency(int, int);
static void store(int, int);
static void store_SiSv_1(int, int);
/*static int check(); quale è il prototipo ?*/
static int coupled(int);
static int generate_out(int, int);
static int ReadCpL(CPLinstance*, CKTcircuit*);
/* static int divC(double, double, double, double, double*, double*); */

/* mult */
static void mult_p(double*, double*, double*, int, int, int);
static void matrix_p_mult(); /* quale è il prototipo ?*/
static double approx_mode(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 Pade_apx(double, double*, double*, double*, double*, double*, double*, double*);
static int Gaussian_Elimination(int);
static double root3(double, double, double, double);
static int div3(double, double, double, double, double*, double*);
static int find_roots(double, double, double, double*, double*,double*);

static NODE* insert_node(char*);
static NDnamePt insert_ND(char*, NDnamePt*);
static NODE* NEW_node(void);
static NDnamePt ndn;
static NODE *node_tab;
#define epsi_mult 1e-28

/* diag */
static MAXE_PTR sort(MAXE_PTR, float, int, int, MAXE_PTR);
static void ordering(void);
static MAXE_PTR delete_1(MAXE_PTR*, int);
static void reordering(int, int);
static void diag(int);
static int rotate(int, int, int);

#define epsi 1.0e-16
static char *message = "tau of coupled lines is larger than max time step";

/* ARGSUSED */
int
CPLsetup(SMPmatrix *matrix, GENmodel *inModel, CKTcircuit *ckt, int *state)
{
    CPLmodel *model = (CPLmodel *)inModel;
    CPLinstance *here;
    CKTnode *tmp, *node;
    int error, m, p;
    char **branchname;
    int noL;


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

        /* loop through all the instances of the model */
        for (here = model->CPLinstances; here != NULL ;
                here=here->CPLnextInstance) {
            


        if (!here->CPLlengthGiven)
            here->CPLlength=0.0;

/* 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);\
}
                    
			noL = here->dimension;

			here->CPLposNodes = (int *) tmalloc(noL * sizeof(int));
			here->CPLnegNodes = (int *) tmalloc(noL * sizeof(int));
			here->CPLibr1 = (int *) tmalloc(noL * sizeof(int));
			here->CPLibr2 = (int *) tmalloc(noL * sizeof(int));

			VECTOR_ALLOC(here->CPLibr1Ibr1, noL);
			VECTOR_ALLOC(here->CPLibr2Ibr2, noL);
			VECTOR_ALLOC(here->CPLposIbr1, noL);
			VECTOR_ALLOC(here->CPLnegIbr2, noL); 
			VECTOR_ALLOC(here->CPLposPos, noL); 
			VECTOR_ALLOC(here->CPLnegNeg, noL);
			VECTOR_ALLOC(here->CPLnegPos, noL); 
			VECTOR_ALLOC(here->CPLposNeg, noL); 

			MATRIX_ALLOC(here->CPLibr1Pos, noL, noL);
			MATRIX_ALLOC(here->CPLibr2Neg, noL, noL);
			MATRIX_ALLOC(here->CPLibr1Neg, noL, noL);
			MATRIX_ALLOC(here->CPLibr2Pos, noL, noL);
			MATRIX_ALLOC(here->CPLibr1Ibr2, noL, noL);
			MATRIX_ALLOC(here->CPLibr2Ibr1, noL, noL);


			branchname = (char **) tmalloc(sizeof(char *) * here->dimension);
			if (! here->CPLibr1Given) {
				for (m = 0; m < here->dimension; m++) {
					branchname[m] = tmalloc(MAX_STRING);
					sprintf(branchname[m], "branch1_%d", m);
					error = 
					CKTmkCur(ckt, &tmp, here->CPLname, branchname[m]);
					if (error) return (error);
					here->CPLibr1[m] = tmp->number;
				}
				here->CPLibr1Given = 1;
			}
			free(branchname);
			branchname = (char **) tmalloc(sizeof(char *) * here->dimension);
			if (! here->CPLibr2Given) {
				for (m = 0; m < here->dimension; m++) {
					branchname[m] = tmalloc(MAX_STRING);
					sprintf(branchname[m], "branch2_%d", m);
					error = 
					CKTmkCur(ckt, &tmp, here->CPLname, branchname[m]);
					if (error) return (error);
					here->CPLibr2[m] = tmp->number;
				}
				here->CPLibr2Given = 1;
			}
			free(branchname);

			for (m = 0; m < here->dimension; m++) {
				for (node = ckt->CKTnodes; node; node = node->next) {
					if (strcmp(here->in_node_names[m], 
						node->name) == 0){
						here->CPLposNodes[m] = node->number;
					}
				}
			}
			for (m = 0; m < here->dimension; m++) {
				for (node = ckt->CKTnodes; node; node = node->next) {
					if (strcmp(here->out_node_names[m],
						node->name) == 0){
						here->CPLnegNodes[m] = node->number;
					}
				}
			}

			for (m = 0; m < here->dimension; m++) {
				TSTALLOC(CPLibr1Ibr1[m],CPLibr1[m],CPLibr1[m]);
				TSTALLOC(CPLibr2Ibr2[m],CPLibr2[m],CPLibr2[m]);
				TSTALLOC(CPLposIbr1[m],CPLposNodes[m],CPLibr1[m]);
	            TSTALLOC(CPLnegIbr2[m],CPLnegNodes[m],CPLibr2[m]);
	            TSTALLOC(CPLposPos[m],CPLposNodes[m],CPLposNodes[m]);
	            TSTALLOC(CPLnegNeg[m],CPLnegNodes[m],CPLnegNodes[m]);
	            TSTALLOC(CPLnegPos[m],CPLnegNodes[m],CPLposNodes[m]);
	            TSTALLOC(CPLposNeg[m],CPLposNodes[m],CPLnegNodes[m]);

				for (p = 0; p < here->dimension; p++) {

            TSTALLOC(CPLibr1Pos[m][p],CPLibr1[m],CPLposNodes[p]);
            TSTALLOC(CPLibr2Neg[m][p],CPLibr2[m],CPLnegNodes[p]);
            TSTALLOC(CPLibr1Neg[m][p],CPLibr1[m],CPLnegNodes[p]);
            TSTALLOC(CPLibr2Pos[m][p],CPLibr2[m],CPLposNodes[p]);
            TSTALLOC(CPLibr1Ibr2[m][p],CPLibr1[m],CPLibr2[p]);
            TSTALLOC(CPLibr2Ibr1[m][p],CPLibr2[m],CPLibr1[p]);

				}
			}

			ReadCpL(here, ckt);

        }
    }

	return(OK);
}



int
CPLunsetup(GENmodel *inModel, CKTcircuit *ckt)
{
  CPLmodel *model;
  CPLinstance *here;
  int m;
  int noL; 
  
  for (model = (CPLmodel *) inModel; model != NULL;
      model = model->CPLnextModel) {
    for (here = model->CPLinstances; here != NULL;
        here = here->CPLnextInstance) {

        noL = here->dimension;

        VECTOR_FREE(here->CPLibr1Ibr1); 
        VECTOR_FREE(here->CPLibr2Ibr2);
        VECTOR_FREE(here->CPLposIbr1);
        VECTOR_FREE(here->CPLnegIbr2);
        VECTOR_FREE(here->CPLposPos);
        VECTOR_FREE(here->CPLnegNeg); 
        VECTOR_FREE(here->CPLnegPos);
        VECTOR_FREE(here->CPLposNeg);

      
     MATRIX_FREE(here->CPLibr1Pos, noL, noL);
     MATRIX_FREE(here->CPLibr2Neg, noL, noL);
     MATRIX_FREE(here->CPLibr1Neg, noL, noL);
     MATRIX_FREE(here->CPLibr2Pos, noL, noL);
     MATRIX_FREE(here->CPLibr1Ibr2, noL, noL);
     MATRIX_FREE(here->CPLibr2Ibr1, noL, noL);

      
      for (m = 0; m < noL; m++) {
          if (here->CPLibr1[m]) {
               CKTdltNNum(ckt, here->CPLibr1[m]);
               here->CPLibr1[m] = 0;
          }
      }

      for (m = 0; m < noL; m++) {
          if (here->CPLibr2[m]) {
               CKTdltNNum(ckt, here->CPLibr2[m]);
               here->CPLibr2[m] = 0;
          }
      }

  free(here->CPLposNodes);
  free(here->CPLnegNodes);
  free(here->CPLibr1);
  free(here->CPLibr2);

  /* reset switches */
  here->CPLdcGiven=0;
  here->CPLibr1Given = 0;
  here->CPLibr2Given = 0;
    }
  }
  return OK;
}



static int 
ReadCpL(CPLinstance *here, CKTcircuit *ckt)
{
   int i, j, noL, counter;
   float f;
   char *name;
   CPLine *c, *c2;
   ECPLine *ec;
   NODE *nd;
   RLINE *lines[MAX_CP_TX_LINES];
   ERLINE *er;

   c = (CPLine *) tmalloc(sizeof (CPLine));
   c2 = (CPLine *) tmalloc(sizeof (CPLine));
   c->vi_head = c->vi_tail = NULL;
   noL = c->noL = here->dimension;
   here->cplines = c;
   here->cplines2 = c2;

   for (i = 0; i < noL; i++) {
      ec = (ECPLine *) tmalloc(sizeof (ECPLine));
	  name = here->in_node_names[i];
      nd = insert_node(name);
      ec->link = nd->cplptr;
      nd->cplptr = ec;
      ec->line = c;
      c->in_node[i] = nd;
      c2->in_node[i] = nd;

      er = (ERLINE *) tmalloc(sizeof (ERLINE));
      er->link = nd->rlptr;
      nd->rlptr = er;
      er->rl = lines[i] = (RLINE *) tmalloc(sizeof (RLINE));
      er->rl->in_node = nd;

      c->dc1[i] = c->dc2[i] = 0.0;
   }

   for (i = 0; i < noL; i++) {
      ec = (ECPLine *) tmalloc(sizeof (ECPLine));
	  name = here->out_node_names[i];
      nd = insert_node(name);
      ec->link = nd->cplptr;
      nd->cplptr = ec;
      ec->line = c;
      c->out_node[i] = nd;
      c2->out_node[i] = nd;

      er = (ERLINE *) tmalloc(sizeof (ERLINE));
      er->link = nd->rlptr;
      nd->rlptr = er;
      er->rl = lines[i];
      er->rl->out_node = nd;
   }

   
	counter = 0;
	for (i = 0; i < noL; i++) {
		for (j = 0; j < noL; j++) {
			if (i > j) {
				R_m[i][j] = R_m[j][i];
				G_m[i][j] = G_m[j][i];
				C_m[i][j] = C_m[j][i];
				L_m[i][j] = L_m[j][i];
			}
			else {
				f = here->CPLmodPtr->Rm[counter];
				R_m[i][j] = here->CPLmodPtr->Rm[counter] = MAX(f, 1.0e-4);
				G_m[i][j] = here->CPLmodPtr->Gm[counter];
				L_m[i][j] = here->CPLmodPtr->Lm[counter];
				C_m[i][j] = here->CPLmodPtr->Cm[counter];
				counter++;
			}
		}
	}
	if (here->CPLlengthGiven)
		length = here->CPLlength;
	else length = here->CPLmodPtr->length;

   for (i = 0; i < noL; i++) 
      lines[i]->g = 1.0 / (R_m[i][i] * length);

   coupled(noL);

   for (i = 0; i < noL; i++) {
      double d, t;
      int k;

      c->taul[i] = TAU[i] * 1.0e+12;
      for (j  = 0; j < noL; j++) {
         if (SIV[i][j].C_0 == 0.0)
            c->h1t[i][j] = NULL;
	 else {
	    c->h1t[i][j] = (TMS *) tmalloc(sizeof (TMS));
	    d = c->h1t[i][j]->aten = SIV[i][j].C_0;
	    c->h1t[i][j]->ifImg = (int) SIV[i][j].Poly[6] - 1.0;
	        /* since originally 2 for img 1 for noimg */
	    c->h1t[i][j]->tm[0].c = SIV[i][j].Poly[0] * d;
	    c->h1t[i][j]->tm[1].c = SIV[i][j].Poly[1] * d;
	    c->h1t[i][j]->tm[2].c = SIV[i][j].Poly[2] * d;
	    c->h1t[i][j]->tm[0].x = SIV[i][j].Poly[3];
	    c->h1t[i][j]->tm[1].x = SIV[i][j].Poly[4];
	    c->h1t[i][j]->tm[2].x = SIV[i][j].Poly[5];
	    if (c->h1t[i][j]->ifImg)
	       c->h1C[i][j] = c->h1t[i][j]->tm[0].c + 2.0 * c->h1t[i][j]->tm[1].c;
            else {
	       t = 0.0;
               for (k = 0; k < 3; k++)
                  t += c->h1t[i][j]->tm[k].c;
	       c->h1C[i][j] = t;
	    }
	 }

         for (k = 0; k < noL; k++) {
            if (IWI[i][j].C_0[k] == 0.0)
               c->h2t[i][j][k] = NULL;
            else {
               c->h2t[i][j][k] = (TMS *) tmalloc(sizeof (TMS));
               d = c->h2t[i][j][k]->aten = IWI[i][j].C_0[k];
               c->h2t[i][j][k]->ifImg = (int) IWI[i][j].Poly[k][6] - 1.0;
                /* since originally 2 for img 1 for noimg */
               c->h2t[i][j][k]->tm[0].c = IWI[i][j].Poly[k][0] * d;
               c->h2t[i][j][k]->tm[1].c = IWI[i][j].Poly[k][1] * d;
               c->h2t[i][j][k]->tm[2].c = IWI[i][j].Poly[k][2] * d;
               c->h2t[i][j][k]->tm[0].x = IWI[i][j].Poly[k][3];
               c->h2t[i][j][k]->tm[1].x = IWI[i][j].Poly[k][4];
               c->h2t[i][j][k]->tm[2].x = IWI[i][j].Poly[k][5];
	       if (c->h2t[i][j][k]->ifImg)
	          c->h2C[i][j][k] = c->h2t[i][j][k]->tm[0].c + 2.0 * 
			  c->h2t[i][j][k]->tm[1].c;
               else 
		  c->h2C[i][j][k] =  c->h2t[i][j][k]->tm[0].c +
				     c->h2t[i][j][k]->tm[1].c +
				     c->h2t[i][j][k]->tm[2].c;
            }
            if (IWV[i][j].C_0[k] == 0.0)
               c->h3t[i][j][k] = NULL;
            else {
               c->h3t[i][j][k] = (TMS *) tmalloc(sizeof (TMS));
               d = c->h3t[i][j][k]->aten = IWV[i][j].C_0[k];
               c->h3t[i][j][k]->ifImg = (int) IWV[i][j].Poly[k][6] - 1.0;
                /* since originally 2 for img 1 for noimg */
               c->h3t[i][j][k]->tm[0].c = IWV[i][j].Poly[k][0] * d;
               c->h3t[i][j][k]->tm[1].c = IWV[i][j].Poly[k][1] * d;
               c->h3t[i][j][k]->tm[2].c = IWV[i][j].Poly[k][2] * d;
               c->h3t[i][j][k]->tm[0].x = IWV[i][j].Poly[k][3];
               c->h3t[i][j][k]->tm[1].x = IWV[i][j].Poly[k][4];
               c->h3t[i][j][k]->tm[2].x = IWV[i][j].Poly[k][5];
	       if (c->h3t[i][j][k]->ifImg)
	          c->h3C[i][j][k] = c->h3t[i][j][k]->tm[0].c + 2.0 * 
			  c->h3t[i][j][k]->tm[1].c;
               else 
		  c->h3C[i][j][k] =  c->h3t[i][j][k]->tm[0].c +
				     c->h3t[i][j][k]->tm[1].c +
				     c->h3t[i][j][k]->tm[2].c;
            }
         }
      }
   }

   for (i = 0; i < noL; i++) {
		if (c->taul[i] < ckt->CKTmaxStep) {
			errMsg = tmalloc(strlen(message)+1);
			strcpy(errMsg,message);
			return(-1);
		}
	}

   return(1);
}


/****************************************************************
     misc.c      Miscellaneous procedures for simulation of
		 coupled transmission lines.
 ****************************************************************/


static void 
new_memory(int dim, int deg, int deg_o)
{
   int i, j;

   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) 
         SiSv_1[i][j] = (double *) calloc(deg_o+1, sizeof(double));

   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) 
         Sip[i][j] = (double *) calloc(deg_o+1, sizeof(double));

   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) 
         Si_1p[i][j] = (double *) calloc(deg_o+1, sizeof(double));

   for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++) 
         Sv_1p[i][j] = (double *) calloc(deg_o+1, sizeof(double));

   for (i = 0; i < dim; i++) 
      W[i] = (double *) calloc(MAX_DEG, sizeof(double));
}

/***
 ***/