1 /**********
2 Copyright 1990 Regents of the University of California. All rights reserved.
3 Author: 1985 Thomas L. Quarles
4 **********/
5
6 #include "ngspice/ngspice.h"
7 #include "ngspice/cktdefs.h"
8
9 #define ccap (qcap+1)
10
11
12 void
CKTterr(int qcap,CKTcircuit * ckt,double * timeStep)13 CKTterr(int qcap, CKTcircuit *ckt, double *timeStep)
14 {
15 double volttol;
16 double chargetol;
17 double tol;
18 double del;
19 double diff[8];
20 double deltmp[8];
21 double factor=0;
22 int i;
23 int j;
24 static double gearCoeff[] = {
25 .5,
26 .2222222222,
27 .1363636364,
28 .096,
29 .07299270073,
30 .05830903790
31 };
32 static double trapCoeff[] = {
33 .5,
34 .08333333333
35 };
36
37 volttol = ckt->CKTabstol + ckt->CKTreltol *
38 MAX( fabs(ckt->CKTstate0[ccap]), fabs(ckt->CKTstate1[ccap]));
39
40 chargetol = MAX(fabs(ckt->CKTstate0[qcap]),fabs(ckt->CKTstate1[qcap]));
41 chargetol = ckt->CKTreltol * MAX(chargetol,ckt->CKTchgtol)/ckt->CKTdelta;
42 tol = MAX(volttol,chargetol);
43 /* now divided differences */
44 for(i=ckt->CKTorder+1;i>=0;i--) {
45 diff[i] = ckt->CKTstates[i][qcap];
46 }
47 for(i=0 ; i <= ckt->CKTorder ; i++) {
48 deltmp[i] = ckt->CKTdeltaOld[i];
49 }
50 j = ckt->CKTorder;
51 for (;;) {
52 for(i=0;i <= j;i++) {
53 diff[i] = (diff[i] - diff[i+1])/deltmp[i];
54 }
55 if (--j < 0) break;
56 for(i=0;i <= j;i++) {
57 deltmp[i] = deltmp[i+1] + ckt->CKTdeltaOld[i];
58 }
59 }
60 switch(ckt->CKTintegrateMethod) {
61 case GEAR:
62 factor = gearCoeff[ckt->CKTorder-1];
63 break;
64
65 case TRAPEZOIDAL:
66 factor = trapCoeff[ckt->CKTorder - 1] ;
67 break;
68 }
69 del = ckt->CKTtrtol * tol/MAX(ckt->CKTabstol,factor * fabs(diff[0]));
70 if(ckt->CKTorder == 2) {
71 del = sqrt(del);
72 } else if (ckt->CKTorder > 2) {
73 del = exp(log(del)/ckt->CKTorder);
74 }
75 *timeStep = MIN(*timeStep,del);
76 return;
77 }
78