From b216db8ffd90f6d0a8d28821ac948f35ecbca260 Mon Sep 17 00:00:00 2001 From: h_vogt Date: Sat, 1 Dec 2012 13:32:27 +0100 Subject: [PATCH] cplsetup.c: indentations etc. --- src/spicelib/devices/cpl/cplsetup.c | 2525 +++++++++++++-------------- 1 file changed, 1262 insertions(+), 1263 deletions(-) diff --git a/src/spicelib/devices/cpl/cplsetup.c b/src/spicelib/devices/cpl/cplsetup.c index 86b9fbd33..f936bd374 100644 --- a/src/spicelib/devices/cpl/cplsetup.c +++ b/src/spicelib/devices/cpl/cplsetup.c @@ -29,7 +29,7 @@ Modified: 2004 Paolo Nenzi - (ng)spice integration } \ } -#define VECTOR_FREE(vec) free(vec) +#define VECTOR_FREE(vec) free(vec) #define MATRIX_FREE(mat, m, j) { \ int k; \ @@ -43,9 +43,9 @@ Modified: 2004 Paolo Nenzi - (ng)spice integration #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 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]; @@ -151,44 +151,44 @@ CPLsetup(SMPmatrix *matrix, GENmodel *inModel, CKTcircuit *ckt, int *state) for( ; model != NULL; model = model->CPLnextModel ) { if (!model->Rmgiven) { - SPfrontEnd->IFerror (ERR_FATAL, - "model %s: lossy line series resistance not given", &(model->CPLmodName)); - return(E_BADPARM); + SPfrontEnd->IFerror (ERR_FATAL, + "model %s: lossy line series resistance not given", &(model->CPLmodName)); + return(E_BADPARM); } if (!model->Gmgiven) { - SPfrontEnd->IFerror (ERR_FATAL, - "model %s: lossy line parallel conductance not given", &(model->CPLmodName)); - return(E_BADPARM); + SPfrontEnd->IFerror (ERR_FATAL, + "model %s: lossy line parallel conductance not given", &(model->CPLmodName)); + return(E_BADPARM); } if (!model->Lmgiven) { - SPfrontEnd->IFerror (ERR_FATAL, - "model %s: lossy line series inductance not given", &(model->CPLmodName)); - return (E_BADPARM); + SPfrontEnd->IFerror (ERR_FATAL, + "model %s: lossy line series inductance not given", &(model->CPLmodName)); + return (E_BADPARM); } if (!model->Cmgiven) { - SPfrontEnd->IFerror (ERR_FATAL, - "model %s: lossy line parallel capacitance not given", &(model->CPLmodName)); - return (E_BADPARM); + SPfrontEnd->IFerror (ERR_FATAL, + "model %s: lossy line parallel capacitance not given", &(model->CPLmodName)); + return (E_BADPARM); } if (!model->lengthgiven) { - SPfrontEnd->IFerror (ERR_FATAL, - "model %s: lossy line length must be given", &(model->CPLmodName)); - return (E_BADPARM); + SPfrontEnd->IFerror (ERR_FATAL, + "model %s: lossy line length must be given", &(model->CPLmodName)); + return (E_BADPARM); } /* loop through all the instances of the model */ for (here = model->CPLinstances; here != NULL ; here=here->CPLnextInstance) { - if (!here->CPLlengthGiven) - here->CPLlength=0.0; + if (!here->CPLlengthGiven) + here->CPLlength=0.0; -/* macro to make elements with built in test for out of memory */ + /* 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)) == NULL){\ return(E_NOMEM);\ } - + noL = here->dimension; here->CPLposNodes = TMALLOC(int, noL); @@ -199,11 +199,11 @@ if((here->ptr = SMPmakeElt(matrix, here->first, here->second)) == NULL){\ 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->CPLnegIbr2, noL); + VECTOR_ALLOC(here->CPLposPos, noL); VECTOR_ALLOC(here->CPLnegNeg, noL); - VECTOR_ALLOC(here->CPLnegPos, noL); - VECTOR_ALLOC(here->CPLposNeg, noL); + VECTOR_ALLOC(here->CPLnegPos, noL); + VECTOR_ALLOC(here->CPLposNeg, noL); MATRIX_ALLOC(here->CPLibr1Pos, noL, noL); MATRIX_ALLOC(here->CPLibr2Neg, noL, noL); @@ -215,70 +215,70 @@ if((here->ptr = SMPmakeElt(matrix, here->first, here->second)) == NULL){\ branchname = TMALLOC(char *, here->dimension); if (! here->CPLibr1Given) { - for (m = 0; m < here->dimension; m++) { - branchname[m] = TMALLOC(char, 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; - tfree(branchname[m]); - } - here->CPLibr1Given = 1; + for (m = 0; m < here->dimension; m++) { + branchname[m] = TMALLOC(char, 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; + tfree(branchname[m]); + } + here->CPLibr1Given = 1; } free(branchname); branchname = TMALLOC(char *, here->dimension); if (! here->CPLibr2Given) { - for (m = 0; m < here->dimension; m++) { - branchname[m] = TMALLOC(char, 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; - tfree(branchname[m]); - } - here->CPLibr2Given = 1; + for (m = 0; m < here->dimension; m++) { + branchname[m] = TMALLOC(char, 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; + tfree(branchname[m]); + } + 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 (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 (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]); + 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++) { + 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]); + 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); @@ -286,7 +286,7 @@ if((here->ptr = SMPmakeElt(matrix, here->first, here->second)) == NULL){\ } } - return(OK); + return(OK); } @@ -294,235 +294,234 @@ if((here->ptr = SMPmakeElt(matrix, here->first, here->second)) == NULL){\ 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; + 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 +static int ReadCpL(CPLinstance *here, CKTcircuit *ckt) { - int i, j, noL, counter; - double f; - char *name; - CPLine *c, *c2; - ECPLine *ec; - NODE *nd; - RLINE *lines[MAX_CP_TX_LINES]; - ERLINE *er; + int i, j, noL, counter; + double f; + char *name; + CPLine *c, *c2; + ECPLine *ec; + NODE *nd; + RLINE *lines[MAX_CP_TX_LINES]; + ERLINE *er; - c = TMALLOC(CPLine, 1); - c2 = TMALLOC(CPLine, 1); - c->vi_head = c->vi_tail = NULL; - noL = c->noL = here->dimension; - here->cplines = c; - here->cplines2 = c2; + c = TMALLOC(CPLine, 1); + c2 = TMALLOC(CPLine, 1); + 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 = TMALLOC(ECPLine, 1); - 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; + for (i = 0; i < noL; i++) { + ec = TMALLOC(ECPLine, 1); + 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 = TMALLOC(ERLINE, 1); - er->link = nd->rlptr; - nd->rlptr = er; - er->rl = lines[i] = TMALLOC(RLINE, 1); - er->rl->in_node = nd; + er = TMALLOC(ERLINE, 1); + er->link = nd->rlptr; + nd->rlptr = er; + er->rl = lines[i] = TMALLOC(RLINE, 1); + er->rl->in_node = nd; - c->dc1[i] = c->dc2[i] = 0.0; - } + c->dc1[i] = c->dc2[i] = 0.0; + } - for (i = 0; i < noL; i++) { - ec = TMALLOC(ECPLine, 1); - 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; + for (i = 0; i < noL; i++) { + ec = TMALLOC(ECPLine, 1); + 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 = TMALLOC(ERLINE, 1); - er->link = nd->rlptr; - nd->rlptr = er; - er->rl = lines[i]; - er->rl->out_node = nd; - } + er = TMALLOC(ERLINE, 1); + 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++; - } - } + + 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; + } + 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); + for (i = 0; i < noL; i++) + lines[i]->g = 1.0 / (R_m[i][i] * length); - coupled(noL); + coupled(noL); - for (i = 0; i < noL; i++) { - double d, t; - int k; + 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] = TMALLOC(TMS, 1); - 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; + 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 { - 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] = TMALLOC(TMS, 1); - 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); + c->h1t[i][j] = TMALLOC(TMS, 1); + 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->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] = TMALLOC(TMS, 1); - 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(char, strlen(message) + 1); - strcpy(errMsg,message); - return(-1); + 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; } - } + } - return(1); + 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] = TMALLOC(TMS, 1); + 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] = TMALLOC(TMS, 1); + 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(char, strlen(message) + 1); + strcpy(errMsg,message); + return(-1); + } + } + + return(1); } @@ -532,109 +531,109 @@ ReadCpL(CPLinstance *here, CKTcircuit *ckt) ****************************************************************/ -static void +static void new_memory(int dim, int deg, int deg_o) { - int i, j; + int i, j; - NG_IGNORE(deg); + NG_IGNORE(deg); - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - SiSv_1[i][j] = (double *) calloc((size_t) (deg_o+1), sizeof(double)); + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + SiSv_1[i][j] = (double *) calloc((size_t) (deg_o+1), sizeof(double)); - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - Sip[i][j] = (double *) calloc((size_t) (deg_o+1), sizeof(double)); + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + Sip[i][j] = (double *) calloc((size_t) (deg_o+1), sizeof(double)); - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - Si_1p[i][j] = (double *) calloc((size_t) (deg_o+1), sizeof(double)); + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + Si_1p[i][j] = (double *) calloc((size_t) (deg_o+1), sizeof(double)); - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - Sv_1p[i][j] = (double *) calloc((size_t) (deg_o+1), sizeof(double)); + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + Sv_1p[i][j] = (double *) calloc((size_t) (deg_o+1), sizeof(double)); - for (i = 0; i < dim; i++) - W[i] = (double *) calloc(MAX_DEG, sizeof(double)); + for (i = 0; i < dim; i++) + W[i] = (double *) calloc(MAX_DEG, sizeof(double)); } /*** ***/ /**************************************************************** - match Create a polynomial matching given data points + match Create a polynomial matching given data points ****************************************************************/ -static double +static double *vector(int nl, int nh) { - double *v; + double *v; - v = TMALLOC(double, (unsigned) (nh - nl + 1)); - if (!v) { - fprintf(stderr, "Memory Allocation Error by tmalloc in vector().\n"); - fprintf(stderr, "...now exiting to system ...\n"); - controlled_exit(EXIT_FAILURE); - } - return v-nl; + v = TMALLOC(double, (unsigned) (nh - nl + 1)); + if (!v) { + fprintf(stderr, "Memory Allocation Error by tmalloc in vector().\n"); + fprintf(stderr, "...now exiting to system ...\n"); + controlled_exit(EXIT_FAILURE); + } + return v-nl; } -static void +static void free_vector(double *v, int nl, int nh) { - NG_IGNORE(nh); + NG_IGNORE(nh); - free((void*) (v +nl)); + free((void*) (v +nl)); } -static void +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 + 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; + 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"); - controlled_exit(EXIT_FAILURE); - } - 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); + 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"); + controlled_exit(EXIT_FAILURE); + } + 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 +static int match(int n, double *cof, double *xa, double *ya) /* Given arrays xa[0..n] and ya[0..n] containing a tabulated function @@ -642,51 +641,51 @@ match(int n, double *cof, double *xa, double *ya) such that ya[i] = sum_j {cof[j]*xa[i]**j}. */ { - int k, j, i; - double xmin, dy, *x, *y, *xx; + 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, 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); - free_vector(xx, 0, n); + 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, 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); + free_vector(xx, 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; + /**** 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; } /*** @@ -718,11 +717,11 @@ match_x(int dim, double *Alfa, double *X, double *Y) } Gaussian_Elimination2(dim, 1); Alfa[0] = Y[0]; - for (i = 1; i <= dim; i++) + 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]; @@ -742,57 +741,57 @@ match_x(int dim, double *Alfa, double *X, double *Y) /*** ***/ -static int +static int Gaussian_Elimination2(int dims, int type) - /* type = 1 : to solve a linear system - -1 : to inverse a matrix */ +/* type = 1 : to solve a linear system + -1 : to inverse a matrix */ { - int i, j, k, dim; - double f; - double max; - int imax; + int i, j, k, dim; + double f; + double max; + int imax; - if (type == -1) - dim = 2 * dims; - else - dim = dims; + 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"); - controlled_exit(EXIT_FAILURE); - } - if (imax != i) - for (k = i; k <= dim; k++) { - f = A[i][k]; - A[i][k] = A[imax][k]; - A[imax][k] = f; - } + 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"); + controlled_exit(EXIT_FAILURE); + } + 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; + f = 1.0 / A[i][i]; + A[i][i] = 1.0; - for (j = i+1; j <= dim; j++) - A[i][j] *= f; + 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]; - } - } + 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); + return(1); } /*** @@ -809,7 +808,7 @@ eval_Si_Si_1(int dim, double y) if (k == j) Si_1[i][j] += Sv_1[i][k] * (y * R_m[k][j] + Scaling_F * L_m[k][j]); - else + else Si_1[i][j] += Sv_1[i][k] * L_m[k][j] * Scaling_F; / Si_1[i][j] *= Scaling_F; @@ -839,36 +838,36 @@ eval_Si_Si_1(int dim, double y) static void eval_Si_Si_1(int dim, double y) { - int i, j, k; + 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] = 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(D[i]); + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + Si_1[i][j] /= sqrt(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++) + 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]; + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + Si[i][j] = A[i][j+dim]; } /*** @@ -880,8 +879,8 @@ loop_ZY(int dim, double y) double fmin, fmin1; for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - if (i == j) + 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]; @@ -907,61 +906,61 @@ loop_ZY(int dim, double y) for (i = 0; i < dim; i++) for (j = 0; j < dim; j++) { Sv_1[i][j] = 0.0; - for (k = 0; k < dim; k++) + 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 (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++) + 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 (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) * + if (k == i) + ZY[i][j] += (Scaling_F * L_m[i][i] + R_m[i] * y) * Y5[k][j]; - else + 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++) + 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]; + 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] = 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] = 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]; + for (j = 0; j < dim; j++) + Sv[i][j] = ZY[i][j]; } ***/ @@ -969,93 +968,93 @@ loop_ZY(int dim, double y) static void loop_ZY(int dim, double y) { - int i, j, k; - double fmin, fmin1=0.0; + int i, j, k; + double fmin, fmin1=0.0; - 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; + 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"); + controlled_exit(EXIT_FAILURE); + } else { + fmin = sqrt(fmin); + fmin1 = 1 / fmin; + } + for (i = 0; i < dim; i++) + D[i] = sqrt(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] = Scaling_F * C_m[i][j]; + else + ZY[i][j] += L_m[i][k] * Y5[k][j] * Scaling_F; */ - 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"); - controlled_exit(EXIT_FAILURE); - } else { - fmin = sqrt(fmin); - fmin1 = 1 / fmin; - } - for (i = 0; i < dim; i++) - D[i] = sqrt(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++) { + 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]; - 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); - 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]; + 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]; } @@ -1063,22 +1062,22 @@ loop_ZY(int dim, double y) /*** ***/ -static void +static void poly_matrix( - double *A[MAX_DIM][MAX_DIM], - int dim, int deg) + double *A[MAX_DIM][MAX_DIM], + int dim, int deg) { - int i, j; + int i, j; - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - match(deg, A[i][j], frequency, A[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 +static int checkW(double *W, double d) { double f, y; @@ -1103,56 +1102,56 @@ checkW(double *W, double d) /*** ***/ -static void +static void poly_W(int dim, int deg) { - int i; + 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]); - */ - } + 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; - double min; + int i; + double min; - min = D[0]; + min = D[0]; - for (i = 1; i < dim; i++) - if (D[i] < min) { - min = D[i]; - } + for (i = 1; i < dim; i++) + if (D[i] < min) { + min = D[i]; + } - if (min <= 0) { - fprintf(stderr, "A mode frequency is not positive. Abort!\n"); - controlled_exit(EXIT_FAILURE); - } + if (min <= 0) { + fprintf(stderr, "A mode frequency is not positive. Abort!\n"); + controlled_exit(EXIT_FAILURE); + } - 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; - */ + 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; + frequency[0] = 0.0; - for (i = 1; i <= deg_o; i++) - frequency[i] = frequency[i-1] + min; + for (i = 1; i <= deg_o; i++) + frequency[i] = frequency[i-1] + min; - for (i = 0; i < dim; i++) - D[i] *= Scaling_F2; + for (i = 0; i < dim; i++) + D[i] *= Scaling_F2; } /*** @@ -1161,20 +1160,20 @@ eval_frequency(int dim, int deg_o) static void store(int dim, int ind) { - int i, j; + 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]; - } + 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]; + } } /*** @@ -1183,22 +1182,22 @@ store(int dim, int ind) static void store_SiSv_1(int dim, int ind) { - int i, j, k; - double temp; + 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; - } + 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 +static int check(Sip, Si_1p, Sv_1p, SiSv_1p) double *Sip[MAX_DIM][MAX_DIM]; double *Si_1p[MAX_DIM][MAX_DIM]; @@ -1278,121 +1277,88 @@ check(Sip, Si_1p, Sv_1p, SiSv_1p) /*** ***/ -static int +static int coupled(int dim) { - int deg, deg_o; - int i; + int deg, deg_o; + int i; - deg = Right_deg; - deg_o = Left_deg; - new_memory(dim, deg, deg_o); + deg = Right_deg; + deg_o = Left_deg; + new_memory(dim, deg, deg_o); - Scaling_F = Scaling_F2 = 1.0; + Scaling_F = Scaling_F2 = 1.0; - /*** y = 0 : ZY = LC ***/ - loop_ZY(dim, 0.0); - eval_frequency(dim, deg_o); - eval_Si_Si_1(dim, 0.0); - store_SiSv_1(dim, 0); - store(dim, 0); + /*** y = 0 : ZY = LC ***/ + loop_ZY(dim, 0.0); + eval_frequency(dim, deg_o); + eval_Si_Si_1(dim, 0.0); + store_SiSv_1(dim, 0); + store(dim, 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); + /*** 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); + poly_matrix(SiSv_1, dim, deg_o); - /*** - check(Sip, Si_1p, Sv_1p, SiSv_1); - ***/ + /*** + check(Sip, Si_1p, Sv_1p, SiSv_1); + ***/ - generate_out(dim, deg_o); + generate_out(dim, deg_o); - return(1); + return(1); } /*** ***/ -static int +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; + 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(sqrt(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(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] = 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]; + 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; - if (i == j && k == i) { - rtv = Pade_apx( - exp(- sqrt(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; - } + continue; + for (k = 0; k <= deg_o; k++) + p[k] /= C; + if (i == j) { + rtv = Pade_apx(sqrt(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(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; - } + rtv = Pade_apx(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; @@ -1400,223 +1366,256 @@ generate_out(int dim, int deg_o) p[4] = x2; p[5] = x3; p[6] = 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(sqrt(G_m[i][i] / R_m[i][i]) * - exp(- sqrt(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(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; - } + } + 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( + exp(- sqrt(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(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] = rtv; } - p[0] = c1; - p[1] = c2; - p[2] = c3; - p[3] = x1; - p[4] = x2; - p[5] = x3; - p[6] = rtv; - } - return(1); + + 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(sqrt(G_m[i][i] / R_m[i][i]) * + exp(- sqrt(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(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] = rtv; + } + return(1); } /**************************************************************** - mult.c Multiplication for Matrix of Polynomial + mult.c Multiplication for Matrix of Polynomial X(y) = A(y) D(y) B(y), - where D(y) is a diagonal matrix with each + 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. + 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; + 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]; - } + 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( - double *A[MAX_DIM][MAX_DIM], - double *D[MAX_DIM], - double *B[MAX_DIM][MAX_DIM], - int dim, int deg, int deg_o, - Mult_Out X[MAX_DIM][MAX_DIM]) + double *A[MAX_DIM][MAX_DIM], + double *D[MAX_DIM], + double *B[MAX_DIM][MAX_DIM], + int dim, int deg, int deg_o, + Mult_Out X[MAX_DIM][MAX_DIM]) { - int i, j, k, l; - double *p; - double *T[MAX_DIM][MAX_DIM]; - double t1; + 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((size_t) (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((size_t) (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++) - tfree(T[i][j]); - - /********** - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) { + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) { + p = T[i][j] = (double *) calloc((size_t) (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++) { - 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"); - } - ***********/ + p = X[i][j].Poly[k] = + (double *) calloc((size_t) (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++) + tfree(T[i][j]); + + /********** + 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 + mode approximation ****************************************************************/ /*** ***/ -static double +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; + 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 */ + 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; + 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(- delay * y1); + delay = sqrt(w0) * length / Scaling_F; + atten = exp(- 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] = 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; - } + a[1] *= -delay; + a[2] *= -delay; + a[3] *= -delay; + a[4] *= -delay; + a[5] *= -delay; - for (i = 0; i <= 5; i++) - b[i] *= atten; + 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; + } - return(delay); + for (i = 0; i <= 5; i++) + b[i] *= atten; + + return(delay); } /*** ***/ -static double +static double eval2(double a, double b, double c, double x) { - return(a*x*x + b*x + c); + 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, +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; + 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; + 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); + return(1); } -static int -Pade_apx(double a_b, double *b, double *c1, double *c2, double *c3, +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. @@ -1625,271 +1624,271 @@ Pade_apx(double a_b, double *b, double *c1, double *c2, double *c3, = c1 / (s - x1) + c2 / (s - x2) + c3 / (s - x3) + 1.0 */ { - double p1, p2, p3, q1, q2, q3; + 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[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[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]; + At[2][0] = b[2]; + At[2][1] = b[3]; + At[2][2] = b[4]; + At[2][3] = -b[5]; - Gaussian_Elimination(3); + 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; + 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(3.0, 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(3.0, 2.0 * p1, p2, *x1); - *c2 = eval2(q1 - p1, q2 - p2, q3 - p3, *x2) / - eval2(3.0, 2.0 * p1, p2, *x2); - *c3 = eval2(q1 - p1, q2 - p2, q3 - p3, *x3) / - eval2(3.0, 2.0 * p1, p2, *x3); - return(1); - } + 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(3.0, 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(3.0, 2.0 * p1, p2, *x1); + *c2 = eval2(q1 - p1, q2 - p2, q3 - p3, *x2) / + eval2(3.0, 2.0 * p1, p2, *x2); + *c3 = eval2(q1 - p1, q2 - p2, q3 - p3, *x3) / + eval2(3.0, 2.0 * p1, p2, *x3); + return(1); + } } -static int +static int Gaussian_Elimination(int dims) { - int i, j, k, dim; - double f; - double max; - int imax; + int i, j, k, dim; + double f; + double max; + int imax; - dim = dims; + 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"); - controlled_exit(EXIT_FAILURE); - } - 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 (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"); + controlled_exit(EXIT_FAILURE); + } + if (imax != i) + for (k = i; k <= dim; k++) { + f = At[i][k]; + At[i][k] = At[imax][k]; + At[imax][k] = f; + } - for (j = i+1; j <= dim; j++) - At[i][j] *= f; + f = 1.0 / At[i][i]; + At[i][i] = 1.0; - 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); + 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 +static double root3(double a1, double a2, double a3, double x) { - double t1, t2; + double t1, t2; - t1 = x * (x * (x + a1) + a2) + a3; - t2 = x * (2.0*a1 + 3.0*x) + a2; + t1 = x * (x * (x + a1) + a2) + a3; + t2 = x * (2.0*a1 + 3.0*x) + a2; - return(x - t1 / t2); + return(x - t1 / t2); } -static int +static int div3(double a1, double a2, double a3, double x, double *p1, double *p2) { - NG_IGNORE(a2); + NG_IGNORE(a2); - *p1 = a1 + x; + *p1 = a1 + x; - /* *p2 = a2 + (a1 + x) * x; */ + /* *p2 = a2 + (a1 + x) * x; */ - *p2 = - a3 / x; + *p2 = - a3 / x; - return(1); + return(1); } -static int +static int find_roots(double a1, double a2, double a3, double *x1, double *x2, double *x3) { - double x, t; - double p, q; + 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(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, 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, 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); - */ + /*********************************************** + 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(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, 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, 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); + *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); - } + 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; + int cmp; + NDnamePt p; - if (*ndn == NULL) { - p = *ndn = TMALLOC(NDname, 1); - 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))); - } + if (*ndn == NULL) { + p = *ndn = TMALLOC(NDname, 1); + 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; + 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); + 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) @@ -1907,30 +1906,30 @@ static int divC(double ar, double ai, double br, double bi, double *cr, double * static NODE *NEW_node(void) { - NODE *n; + NODE *n; - n = TMALLOC(NODE, 1); - 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; + n = TMALLOC(NODE, 1); + 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); + return(n); } @@ -1944,143 +1943,143 @@ static NODE static int dim; static MAXE_PTR row; -static MAXE_PTR +static MAXE_PTR sort(MAXE_PTR list, double 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); - } + 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; - double mv; + MAXE_PTR e; + int i, j, m; + double 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)) { + 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 = TMALLOC(MAXE, 1); - row = sort(row, mv, i, m, e); - } + mv = ABS(ZY[i][j]); + m = j; + } + e = TMALLOC(MAXE, 1); + row = sort(row, mv, i, m, e); + } } - -static MAXE_PTR + +static MAXE_PTR delete_1(MAXE_PTR *list, int rc) { - MAXE_PTR list1, e; + 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); + 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; - double mv; + MAXE_PTR e; + int j, m; + double mv; - m = p+1; - mv = ABS(ZY[p][m]); - for (j = m+1; j < dim; j++) + 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]); + mv = ABS(ZY[p][j]); m = j; - } - e = delete_1(&row, q); - row = sort(row, mv, q, m, e); - } + } + 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 +static void diag(int dims) { - int i, j, c; - double fmin, fmax; + int i, j, c; + double fmin, fmax; - dim = dims; - row = NULL; + 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; + 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; - } + 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(); + ordering(); - if (row) - for (c = 0; row->value > epsi2; c++) { - int p, q; + if (row) + for (c = 0; row->value > epsi2; c++) { + int p, q; - p = row->row; - q = row->col; + p = row->row; + q = row->col; - rotate(dim, p, q); - reordering(p, q); - } + rotate(dim, p, q); + reordering(p, q); + } - for (i = 0; i < dim; i++) - D[i] = ZY[i][i] / fmin; + for (i = 0; i < dim; i++) + D[i] = ZY[i][i] / fmin; - while (row) { - MAXE_PTR tmp_row = row->next; - tfree(row); - row = tmp_row; - } + while (row) { + MAXE_PTR tmp_row = row->next; + tfree(row); + row = tmp_row; + } } /**************************************************************** @@ -2090,68 +2089,68 @@ diag(int dims) static int rotate(int dim, int p, int q) { - int j; - double co, si; - double ve, mu, ld; - double T[MAX_DIM]; - double t; + 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(ld*ld + mu*mu); - co = sqrt((ve + ABS(mu)) / (2.0 * ve)); - si = SGN(mu) * ld / (2.0 * ve * co); + ld = - ZY[p][q]; + mu = 0.5 * (ZY[p][p] - ZY[q][q]); + ve = sqrt(ld*ld + mu*mu); + co = sqrt((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++) + 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; - } + 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; + 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; + ZY[p][q] = 0.0; - { - double R[MAX_DIM]; + { + 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++) { + 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; - } - } + 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); + return(1); }