-- File: [Thyme]<Thyme>System>CSIM01>spSolve.mesa
-- Last editted:
-- by Barth, July 11, 1983 1:42 PM
-- Chen, April 22, 1983 10:56 AM. Changed dvdtmax from 1E13 to 1E15.
-- Wilhelm April 27, 1982 4:04 PM, reformated by Barth and stored under
-- [Cherry]<Barth>Thyme>1.97> .
DIRECTORY spGlobals, spModelDefs, AltoDefs, Real, CWF, RealFns, StringDefs;
spSolve: PROGRAM
IMPORTS spGlobals, Real, CWF, RF: RealFns, StringDefs
EXPORTS spGlobals, spModelDefs =
BEGIN
OPEN spGlobals;
intNodeList: PUBLIC nodePtr ← NIL;
intNodeModFunc: PUBLIC modFuncPtr ← NIL;
otherModFunc: PUBLIC modFuncPtr ← NIL;
solveError: SIGNAL[s: STRING] = CODE;
solveTest: SIGNAL[s: STRING] = CODE;
Retreat: PUBLIC SIGNAL[cause: STRING] = CODE;
Failure: PUBLIC SIGNAL[errorNum: CARDINAL] = CODE;
worstNode: nodePtr;
maxLog: CARDINAL = 16;
curLog: CARDINAL ← 0;
worstNodeLog: ARRAY[0..maxLog) OF RECORD[node: nodePtr,
t: REAL,
v,
dvdt: REAL];
icsSet: BOOLEAN ← FALSE;
checkPoint: BOOLEAN ← FALSE;
canned: BOOLEAN ← FALSE;
numGoodSteps: CARDINAL ← 0;
modelsSkipped: CARDINAL ← 0;
dvdtTrunc: REAL;
t: REAL ← 0.0;
dT: REAL;
AMcount: CARDINAL ← 0;
goodNodeCount: CARDINAL ← 0;
numberOfSteps: LONG CARDINAL;
highStepCount: LONG CARDINAL;
lowStepCount: LONG CARDINAL;
countsPerStep: LONG CARDINAL ← 100000000;
maxCountsPerStep: LONG CARDINAL = 1000000000;
nParms: CARDINAL = 28;
parmNames: ARRAY[0..nParms) OF STRING ← ["maxiter",
"printiter",
"initstep",
"squish",
"tol",
"dt",
"tmax",
"printstep",
"rungekutta",
"tmin",
"printall",
"amiter",
"amup",
"inttol",
"amdown",
"vmin",
"vmax",
"dvdtmax",
"dvdtfac",
"tbreak",
"rkerr",
"amerr",
"amdelay",
"floor",
"minratio",
"rkfac",
"retratio",
""];
Vmin: REAL ← -100.0;
Vmax: REAL ← 100.0;
dvdtmax: REAL ← 1.0e15;
dfltInts: CARDINAL ← 200;
retRatio: REAL ← 0.5;
RKfac: REAL ← 1.0;
minRatio: REAL ← 0.8;
floor: REAL ← 0.01;
AMerr: REAL ← 0.01;
RKerr: REAL ← 0.01;
tBreak: REAL ← 0.0;
dvdtFac: REAL ← 5.0;
squish: REAL ← 0.2;
initStep: REAL ← 0.1;
tol: REAL ← 0.001;
intTol: REAL ← 0.001;
AMup: REAL ← 1.25;
AMdown: REAL ← 0.8;
printdT: REAL ← 0.0;
tMin: REAL ← 0.0;
tMax: REAL ← 0.0;
maxIter: CARDINAL ← 20;
AMiter: CARDINAL ← 2;
AMdelay: CARDINAL ← 5;
printIter: BOOLEAN ← FALSE;
printStep: BOOLEAN ← FALSE;
doRungeKutta: BOOLEAN ← FALSE;
printAll: BOOLEAN ← FALSE;
printMsg: PROCEDURE[fmt: STRING, a1, a2, a3, a4: LONG POINTER ← NIL] =
BEGIN
s: STRING = [128];
CWF.SWF[s, fmt, a1, a2, a3, a4];
printSysWindow[s]
END;
updateFunctions: PROCEDURE[f: modFuncPtr] =
BEGIN
v: REAL;
mb: modBrPtr;
source: argSource;
args, res: argList;
allInactive: BOOLEAN;
modelsSkipped ← 0;
UNTIL f = NIL DO
source ← f↑.arguments;
args ← f↑.argVector;
allInactive ← TRUE;
FOR i: CARDINAL IN [0..LENGTH[source]) DO
args[i] ← source[i]↑.nHist.y;
allInactive ← allInactive AND source[i]↑.marked
ENDLOOP;
WITH mf: f↑ SELECT FROM
fn =>
BEGIN
v ← mf.functionProc[t, args, mf.parmVector];
mf.branch↑.comVal ← v
END;
mod =>
IF ~allInactive THEN
BEGIN
res ← mf.modelResults;
mf.modelProc[args, mf.parmVector, res];
mb ← mf.modelBranches;
UNTIL mb = NIL DO
v ← mb↑.b↑.comVal ← res[mb↑.b↑.modelIndex];
mb ← mb↑.nextBranch
ENDLOOP
END
ELSE modelsSkipped ← modelsSkipped + 1
ENDCASE;
f ← f↑.nextFunction
ENDLOOP;
END;
GSiteration: PROCEDURE[dT: REAL,
print, firstTime, mark: BOOLEAN ← FALSE]
RETURNS[nIter: CARDINAL ← 0] =
BEGIN
E, newE, dvdt, I, newI, cond, csum, otherV, maxI, denom: REAL;
n: nodePtr;
b: branchPtr;
vCur: vSourcePtr;
bLink: branchLinkPtr;
plus, nodeOK, OK: BOOLEAN;
updateFunctions[intNodeModFunc];
IF print THEN printNodeEqn[];
UNTIL nIter = maxIter + maxIter DO
OK ← TRUE;
IF print THEN
CWF.WF2["t = %10.4f, iteration %2d --*n", @t, @nIter];
IF otherModFunc # NIL THEN updateFunctions[otherModFunc];
IF firstTime THEN updateFunctions[intNodeModFunc];
FOR n ← nodeList, n↑.nextNode UNTIL n = NIL DO
IF n = gndNode THEN LOOP;
IF nIter > 1 AND n↑.converged THEN LOOP;
IF print THEN CWF.WF1["Node %s: ", n↑.nodeName↑.name];
bLink ← n↑.branches;
IF n↑.integrate THEN
BEGIN
maxI ← csum ← I ← 0.0;
UNTIL bLink = NIL DO
b ← bLink↑.branch;
plus ← bLink↑.pos;
WITH b↑ SELECT FROM
resistor =>
newI ← (bLink↑.otherNode↑.nHist.y -
n↑.nHist.y)/comVal;
capacitor =>
BEGIN
csum ← csum + comVal;
newI ← comVal*bLink↑.otherNode↑.nHist.f0
END;
inductor =>
newI ← (IF plus THEN -iHist.y
ELSE iHist.y)/comVal;
vSource =>
newI ← IF plus THEN -vsCurrent ELSE vsCurrent;
iSource =>
newI ← IF plus THEN -comVal ELSE comVal
ENDCASE;
maxI ← MAX[maxI, ABS[newI]];
I ← I + newI;
bLink ← bLink↑.nextLink
ENDLOOP;
E ← n↑.nHist.f0;
dvdt ← n↑.nHist.f0 ← I/csum;
IF mark THEN n↑.marked ← ABS[I] <= dvdtTrunc*maxI;
denom ← MAX[ABS[E], ABS[dvdt]];
nodeOK ← IF denom <= Real.SmallestNormalizedNumber THEN TRUE
ELSE tol >= ABS[(dvdt - E)/denom];
IF print THEN
CWF.WF3["Old v' = %9.2f New v' = %9.2f C = %9.2f.*n",
@E, @dvdt, @csum]
END
ELSE
BEGIN
E ← n↑.nHist.y;
IF n↑.curPtr # NIL THEN
BEGIN
vCur ← n↑.curPtr;
newE ← IF vCur↑.posNode = n
THEN vCur↑.negNode↑.nHist.y + vCur↑.comVal
ELSE vCur↑.posNode↑.nHist.y - vCur↑.comVal;
I ← vCur↑.vsCurrent;
newI ← 0.0;
UNTIL bLink = NIL DO
b ← bLink↑.branch;
otherV ← bLink↑.otherNode↑.nHist.y;
plus ← bLink↑.pos;
WITH b↑ SELECT FROM
resistor =>
newI ← newI + (otherV - newE)/comVal;
capacitor =>
newI ← newI + comVal*
(bLink↑.otherNode↑.nHist.f0 -
n↑.nHist.f0);
inductor =>
newI ← newI + (IF plus THEN -iHist.y
ELSE iHist.y)/comVal;
vSource =>
newI ← IF plus THEN newI - vsCurrent
ELSE newI + vsCurrent;
iSource =>
newI ← IF plus THEN newI - comVal
ELSE newI + comVal
ENDCASE;
bLink ← bLink↑.nextLink
ENDLOOP;
newI ← IF n = vCur↑.posNode THEN I + newI
ELSE I - newI;
denom ← MAX[ABS[I], ABS[newI]];
nodeOK ← IF denom <= Real.SmallestNormalizedNumber THEN TRUE
ELSE tol >= ABS[(newI - I)/denom];
IF firstTime THEN n↑.nHist.f0 ← 0.0
ELSE
n↑.nHist.f0 ← squish*n↑.nHist.f1 + (1.0-squish)*
(newE-n↑.nHist.oldy)/dT +
(IF vCur↑.posNode = n
THEN vCur↑.negNode↑.nHist.f0
ELSE vCur↑.posNode↑.nHist.f0);
vCur↑.vsCurrent ← newI;
IF print THEN
CWF.WF2["Old I = %9.2f New I = %9.2f", @I, @newI]
END
ELSE
BEGIN
cond ← newE ← 0.0;
UNTIL bLink = NIL DO
b ← bLink↑.branch;
plus ← bLink↑.pos;
WITH b↑ SELECT FROM
resistor =>
BEGIN
cond ← cond + 1.0/comVal;
newE ← newE +
bLink↑.otherNode↑.nHist.y/comVal
END;
inductor =>
newE ← newE + (IF plus THEN -iHist.y
ELSE iHist.y)/comVal;
vSource =>
newE ← IF plus THEN newE - vsCurrent
ELSE newE + vsCurrent;
iSource =>
newE ← IF plus THEN newE - comVal
ELSE newE + comVal
ENDCASE;
bLink ← bLink↑.nextLink
ENDLOOP;
newE ← newE/cond;
denom ← MAX[ABS[E], ABS[newE]];
nodeOK ← IF denom <= Real.SmallestNormalizedNumber THEN TRUE
ELSE tol >= ABS[(newE - E)/denom]
END;
n↑.nHist.y ← newE;
IF print THEN
CWF.WF2["Old E = %9.2f New E = %9.2f.*n", @E, @newE]
END;
OK ← OK AND nodeOK;
n↑.converged ← nodeOK;
IF ~nodeOK THEN
FOR bLink ← n↑.branches, bLink↑.nextLink UNTIL bLink=NIL DO
bLink↑.otherNode↑.converged ← FALSE
ENDLOOP
ENDLOOP;
nIter ← nIter + 1;
IF nIter = maxIter THEN
printMsg["t = %9.3f, iteration count exceeded", @t];
IF OK THEN EXIT
ENDLOOP;
IF nIter >= maxIter + maxIter THEN
SIGNAL solveError["GS did not converge"]
END;
printNodeEqn: PROCEDURE =
BEGIN
n: nodePtr;
bLink: branchLinkPtr;
branch: branchPtr;
plus: BOOLEAN;
FOR n ← intNodeList, n↑.nextIntNode UNTIL n = NIL DO
CWF.WF3["*n%s -- %13.6f %13.6f*n", n↑.nodeName↑.name,
@n↑.nHist.y, @n↑.nHist.oldy];
FOR bLink ← n↑.branches, bLink↑.nextLink UNTIL bLink = NIL DO
branch ← bLink↑.branch;
plus ← bLink↑.pos;
WITH b: branch↑ SELECT FROM
resistor =>
CWF.WF3["(%13.6f - %13.6f)/%13.6f",
@bLink↑.otherNode↑.nHist.y,
@n↑.nHist.y,@b.comVal];
capacitor =>
CWF.WF2["%13.6fx%13.6f",
@b.comVal, @bLink↑.otherNode↑.nHist.f0];
inductor =>
IF plus THEN
CWF.WF2["-1x%13.6f/%13.6f",@b.iHist.y,@b.comVal]
ELSE CWF.WF2["%13.6f/%13.6f", @b.iHist.y, @b.comVal];
vSource =>
IF plus THEN CWF.WF1["-1x%13.6f", @b.vsCurrent]
ELSE CWF.WF1["%13.6f", @b.vsCurrent];
iSource =>
IF plus THEN CWF.WF1["-1x%13.6f", @b.comVal]
ELSE CWF.WF1["%13.6f", @b.comVal]
ENDCASE;
CWF.WF3[" %s(%13.6f, %13.6f)*n",
bLink↑.otherNode↑.nodeName↑.name,
@bLink↑.otherNode↑.nHist.y,
@bLink↑.otherNode↑.nHist.f0]
ENDLOOP
ENDLOOP
END;
parmValue: PROCEDURE[bool: BOOLEAN]
RETURNS[r: REAL, b: BOOLEAN, c: CARDINAL] =
BEGIN
IF bool THEN
IF item = name THEN
BEGIN
IF LongEqualStrings[newString, "FALSE"] THEN b ← FALSE
ELSE
IF LongEqualStrings[newString, "TRUE"] THEN b ← TRUE
ELSE error[719];
next[]
END
ELSE error[704]
ELSE
BEGIN
r ← getSignedNumber[];
IF r >= 0.0 AND r < 65535.0 THEN c ← Real.FixC[r]
END
END;
runParms: PROCEDURE =
BEGIN
nameIndex: CARDINAL;
WHILE item = name DO
FOR i: CARDINAL IN [0..nParms) DO
nameIndex ← i;
IF newString.length # parmNames[nameIndex].length THEN LOOP;
FOR j: CARDINAL IN [0..newString.length) DO
IF parmNames[i][j] # StringDefs.LowerCase[newString[j]] THEN EXIT;
REPEAT
FINISHED => EXIT;
ENDLOOP;
ENDLOOP;
next[];
IF item # leftArrow THEN error[703] ELSE next[];
SELECT nameIndex FROM
0 => [,, maxIter] ← parmValue[FALSE];
1 => [, printIter,] ← parmValue[TRUE];
2 => [initStep,,] ← parmValue[FALSE];
3 => [initStep,,] ← parmValue[FALSE];
4 => [tol,,] ← parmValue[FALSE];
5 => [printdT,,] ← parmValue[FALSE];
6 => [tMax,,] ← parmValue[FALSE];
7 => [, printStep,] ← parmValue[TRUE];
8 => [, doRungeKutta,] ← parmValue[TRUE];
9 => [tMin,,] ← parmValue[FALSE];
10 => [, printAll,] ← parmValue[TRUE];
11 => [,, AMiter] ← parmValue[FALSE];
12 => [AMup,,] ← parmValue[FALSE];
13 => [intTol,,] ← parmValue[FALSE];
14 => [AMdown,,] ← parmValue[FALSE];
15 => [Vmin,,] ← parmValue[FALSE];
16 => [Vmax,,] ← parmValue[FALSE];
17 => [dvdtmax,,] ← parmValue[FALSE];
18 => [dvdtFac,,] ← parmValue[FALSE];
19 => [tBreak,,] ← parmValue[FALSE];
20 => [RKerr,,] ← parmValue[FALSE];
21 => [AMerr,,] ← parmValue[FALSE];
22 => [,, AMdelay] ← parmValue[FALSE];
23 => [floor,,] ← parmValue[FALSE];
24 => [minRatio,,] ← parmValue[FALSE];
25 => [RKfac,,] ← parmValue[FALSE];
26 => [retRatio,,] ← parmValue[FALSE];
ENDCASE => error[790];
IF item = comma THEN next[]
ENDLOOP
END;
predict: PROCEDURE[h: REAL] =
BEGIN
nodes: nodePtr ← intNodeList;
inds: inductorPtr ← inductorList;
incr: REAL;
h ← h/24.0;
UNTIL nodes = NIL DO
IF ~nodes↑.marked THEN
BEGIN
OPEN nodes↑.nHist;
incr ← h*(55.0*f1 - 59.0*f2 + 37.0*f3 - 9.0*f4);
y ← oldy + incr
END;
nodes ← nodes↑.nextIntNode
ENDLOOP;
UNTIL inds = NIL DO
OPEN inds↑.iHist;
y ← oldy + h*(55.0*f1 - 59.0*f2 + 37.0*f3 - 9.0*f4);
inds ← inds↑.nextInductor
ENDLOOP;
END;
firstCorrector: PROCEDURE[h: REAL] =
BEGIN
nodes: nodePtr ← intNodeList;
inds: inductorPtr ← inductorList;
f: REAL;
h ← h/24.0;
UNTIL nodes = NIL DO
IF ~nodes↑.marked THEN
BEGIN
OPEN nodes↑.nHist;
y ← oldy + h*(9.0*f0 + 19.0*f1 - 5.0*f2 + f3);
oldf ← f0
END;
nodes ← nodes↑.nextIntNode
ENDLOOP;
UNTIL inds = NIL DO
OPEN inds↑.iHist;
f ← inds↑.posNode↑.nHist.y - inds↑.negNode↑.nHist.y;
y ← oldy + h*(9.0*f + 19.0*f1 - 5.0*f2 + f3);
oldf ← f;
inds ← inds↑.nextInductor
ENDLOOP;
END;
iteratingCorrector: PROCEDURE[h: REAL]
RETURNS[newRatio: REAL ← 4.0,
OK: BOOLEAN ← TRUE] =
BEGIN
nodes: nodePtr ← intNodeList;
inds: inductorPtr ← inductorList;
f, deltay, newy, err, denom: REAL;
h ← 3.0*h/8.0;
UNTIL nodes = NIL DO
IF ~nodes↑.marked THEN
BEGIN
OPEN nodes↑.nHist;
deltay ← h*(f0 - oldf);
newy ← y + deltay;
denom ← MAX[ABS[newy], ABS[y]];
IF denom > floor THEN
BEGIN
err ← AMerr*ABS[deltay/denom];
IF intTol < err THEN
BEGIN
OK ← FALSE;
newRatio ← MIN[newRatio, intTol/err];
worstNode ← nodes
END
END;
y ← newy;
oldf ← f0
END;
nodes ← nodes↑.nextIntNode
ENDLOOP;
UNTIL inds = NIL DO
OPEN inds↑.iHist;
f ← inds↑.posNode↑.nHist.y - inds↑.negNode↑.nHist.y;
deltay ← h*(f - oldf);
newy ← y + deltay;
denom ← MAX[ABS[newy], ABS[y]];
IF denom # 0.0 THEN
BEGIN
err ← AMerr*ABS[deltay/denom];
IF intTol < err THEN
BEGIN
OK ← FALSE;
newRatio ← MIN[newRatio, intTol/err]
END
END;
y ← newy;
oldf ← f;
inds ← inds↑.nextInductor
ENDLOOP;
END;
AdamsMoulton: PROCEDURE[dT: REAL]
RETURNS[newRatio: REAL, nGS: CARDINAL, ok: BOOLEAN] =
BEGIN
nAM: CARDINAL;
predict[dT];
nGS ← GSiteration[dT, printAll];
firstCorrector[dT];
nGS ← GSiteration[dT, printAll] + nGS;
FOR nAM ← 2, nAM + 1 UNTIL nAM = AMiter DO
[] ← iteratingCorrector[dT];
nGS ← GSiteration[dT, printAll] + nGS
ENDLOOP;
[newRatio, ok] ← iteratingCorrector[dT];
END;
RKupdate: PROCEDURE[c, h: REAL, last: BOOLEAN]
RETURNS[newRatio: REAL ← 1.0, ok: BOOLEAN ← TRUE] =
BEGIN
nodes: nodePtr ← intNodeList;
inds: inductorPtr ← inductorList;
k, incr, y, err, denom: REAL;
UNTIL nodes = NIL DO
IF ~nodes↑.marked THEN
BEGIN
nodes↑.nHist.sumk ← c*nodes↑.nHist.f0 + nodes↑.nHist.sumk;
IF last THEN
BEGIN
incr ← h*nodes↑.nHist.sumk;
y ← nodes↑.nHist.oldy + incr;
denom ← MAX[ABS[y], ABS[nodes↑.nHist.y]];
IF denom > floor THEN
BEGIN
err ← RKerr*ABS[(y - nodes↑.nHist.y)/denom];
IF intTol < err THEN
BEGIN
ok ← FALSE;
newRatio ← MIN[newRatio, intTol/err]
END
END;
nodes↑.nHist.y ← y;
nodes↑.nHist.sumk ← 0.0
END
ELSE
BEGIN
incr ← h*nodes↑.nHist.f0;
nodes↑.nHist.y ← nodes↑.nHist.oldy + incr
END
END;
nodes ← nodes↑.nextIntNode
ENDLOOP;
UNTIL inds = NIL DO
k ← inds↑.posNode↑.nHist.y - inds↑.negNode↑.nHist.y;
inds↑.iHist.sumk ← c*k + inds↑.iHist.sumk;
IF last THEN
BEGIN
inds↑.iHist.y ← inds↑.iHist.oldy + h*inds↑.iHist.sumk;
inds↑.iHist.sumk ← 0.0
END
ELSE inds↑.iHist.y ← inds↑.iHist.oldy + h*k;
inds ← inds↑.nextInductor
ENDLOOP;
END;
RungeKutta: PROCEDURE[dT: REAL]
RETURNS[newRatio: REAL, nGS: CARDINAL, ok: BOOLEAN] =
BEGIN
time: REAL;
time ← t;
t ← t - 0.5*dT;
[] ← RKupdate[1.0, 0.5*dT, FALSE];
nGS ← GSiteration[0.5*dT, printAll];
[] ← RKupdate[2.0, 0.5*dT, FALSE];
nGS ← GSiteration[0.5*dT, printAll] + nGS;
t ← time;
[] ← RKupdate[2.0, dT, FALSE];
nGS ← GSiteration[dT, printAll] + nGS;
[newRatio, ok] ← RKupdate[1.0, 0.1666667*dT, TRUE]
END;
changeStepSize: PROCEDURE[ratio: REAL, RK: BOOLEAN] =
BEGIN
int: CHARACTER;
IF ratio < 1.0 OR countsPerStep < maxCountsPerStep THEN
BEGIN
numGoodSteps ← 1;
IF worstNode # NIL THEN
BEGIN
worstNodeLog[curLog].node ← worstNode;
worstNodeLog[curLog].t ← t;
worstNodeLog[curLog].v ← worstNode↑.nHist.y;
worstNodeLog[curLog].dvdt ← worstNode↑.nHist.f0
END
ELSE worstNodeLog[curLog] ← [NIL, 0.0, 0.0, 0.0];
curLog ← (curLog + 1) MOD maxLog;
IF ratio > Real.Float[maxCountsPerStep/countsPerStep] THEN
countsPerStep ← maxCountsPerStep
ELSE
countsPerStep ← MAX[1, Real.RoundLI[countsPerStep*ratio]];
int ← IF RK THEN 'R ELSE 'A;
printMsg["t = %10.3f, step = %10lu%c",
@t, @countsPerStep, @int];
IF countsPerStep = 1 THEN SIGNAL solveError["Step too small"]
END
END;
saveState: PROCEDURE =
BEGIN
nodes: nodePtr ← nodeList;
inds: inductorPtr ← inductorList;
volts: vSourcePtr ← vSourceList;
UNTIL nodes = NIL DO
OPEN nodes↑.nHist;
oldy ← y;
f4 ← f3; f3 ← f2; f2 ← f1; f1 ← f0;
IF y < Vmin THEN SIGNAL solveTest["Node below Vmin"]
ELSE
IF y > Vmax THEN SIGNAL solveTest["Node above Vmax"]
ELSE
IF ABS[f0] > dvdtmax THEN
SIGNAL solveTest["Rate above max"];
nodes ← nodes↑.nextNode
ENDLOOP;
UNTIL inds = NIL DO
OPEN inds↑.iHist;
oldy ← y;
f4 ← f3; f3 ← f2; f2 ← f1;
f1 ← inds↑.posNode↑.nHist.y - inds↑.negNode↑.nHist.y;
inds ← inds↑.nextInductor
ENDLOOP;
UNTIL volts = NIL DO
volts↑.oldCurrent ← volts↑.vsCurrent;
volts ← volts↑.nextvSource
ENDLOOP
END;
restoreState: PROCEDURE =
BEGIN
nodes: nodePtr ← nodeList;
inds: inductorPtr ← inductorList;
volts: vSourcePtr ← vSourceList;
UNTIL nodes = NIL DO
nodes↑.nHist.y ← nodes↑.nHist.oldy;
nodes↑.nHist.f0 ← nodes↑.nHist.f1;
nodes ← nodes↑.nextNode
ENDLOOP;
UNTIL inds = NIL DO
inds↑.iHist.y ← inds↑.iHist.oldy;
inds ← inds↑.nextInductor
ENDLOOP;
UNTIL volts = NIL DO
volts↑.vsCurrent ← volts↑.oldCurrent;
volts ← volts↑.nextvSource
ENDLOOP
END;
integration: PROCEDURE[RK: BOOLEAN] =
BEGIN
nGS: CARDINAL;
out, ok: BOOLEAN;
ratio, oldTime, DT: REAL;
oldLowStepcount, oldHighStepCount: LONG CARDINAL;
saveState[];
oldTime ← t;
oldLowStepcount ← lowStepCount;
oldHighStepCount ← highStepCount;
DO
lowStepCount ← lowStepCount + countsPerStep;
IF lowStepCount > maxCountsPerStep THEN
BEGIN
highStepCount ← highStepCount + 1;
lowStepCount ← lowStepCount - maxCountsPerStep;
out ← TRUE
END
ELSE out ← printStep;
t ← tMin + dT*lowStepCount + (dT*maxCountsPerStep)*highStepCount;
DT ← dT*countsPerStep;
BEGIN
ENABLE Retreat =>
BEGIN
printSysWindow[cause];
ok ← FALSE;
ratio ← retRatio;
GOTO FlushIt
END;
IF numGoodSteps < 4 THEN
BEGIN
[ratio, nGS, ok] ← RungeKutta[DT];
IF ok THEN nGS ← GSiteration[DT, printAll,, TRUE] + nGS;
ratio ← RF.SqRt[ratio]*RKfac;
IF ok AND ~RK THEN numGoodSteps ← numGoodSteps + 1
END
ELSE
BEGIN
AMcount ← AMcount + 1;
[ratio, nGS, ok] ← AdamsMoulton[DT];
IF ok THEN
BEGIN
nGS ← GSiteration[DT, printAll,, TRUE] + nGS;
IF ratio > AMup AND AMcount >= AMdelay THEN
BEGIN
ratio ← RF.SqRt[RF.SqRt[ratio]];
changeStepSize[ratio, FALSE];
AMcount ← 0
END
END
ELSE ratio ← RF.SqRt[ratio]
END
EXITS
FlushIt => NULL
END;
IF ok THEN EXIT;
restoreState[];
t ← oldTime;
highStepCount ← oldHighStepCount;
lowStepCount ← oldLowStepcount;
changeStepSize[MIN[ratio, minRatio], numGoodSteps < 4];
AMcount ← 0
ENDLOOP;
advanceCursor[];
IF out THEN printFromList[nGS, t, printIter];
IF out THEN plotFromList[t];
WHILE tBreak > 0.0 AND t >= tBreak DO oldTime ← t ENDLOOP
END;
zeroVars: PROCEDURE =
BEGIN
n: nodePtr ← nodeList;
inds: inductorPtr ← inductorList;
vSs: vSourcePtr ← vSourceList;
gndNode↑.marked ← TRUE;
UNTIL n = NIL DO
n↑.nHist ← initialHistory;
n↑.marked ← FALSE;
n ← n↑.nextNode
ENDLOOP;
UNTIL inds = NIL DO
inds↑.iHist ← initialHistory;
inds ← inds↑.nextInductor
ENDLOOP;
UNTIL vSs = NIL DO
vSs↑.vsCurrent ← vSs↑.oldCurrent ← 0.0;
IF vSs↑.controller = NIL THEN
IF vSs↑.posNode = gndNode THEN vSs↑.negNode↑.marked ← TRUE
ELSE
IF vSs↑.negNode = gndNode THEN vSs↑.posNode↑.marked ← TRUE;
vSs ← vSs↑.nextvSource
ENDLOOP
END;
getStepSize: PROCEDURE RETURNS[stepOK: BOOLEAN] =
BEGIN
vMax, dvdtMax, frac: REAL ← 0.0;
inds: inductorPtr ← inductorList;
nodes: nodePtr ← nodeList;
lowStepCount ← 0;
highStepCount ← 0;
numberOfSteps ← Real.FixC[tMax/printdT];
worstNode ← NIL;
UNTIL nodes = NIL DO
vMax ← MAX[ABS[nodes↑.nHist.y], vMax];
IF ABS[nodes↑.nHist.f0] > dvdtMax THEN
BEGIN
dvdtMax ← ABS[nodes↑.nHist.f0];
worstNode ← nodes
END;
nodes ← nodes↑.nextNode
ENDLOOP;
frac ← IF dvdtMax > 0.0 AND vMax > 0.0 THEN intTol*vMax/(dvdtMax*printdT)
ELSE initStep;
dT ← printdT/maxCountsPerStep;
dvdtTrunc ← dvdtFac*intTol;
countsPerStep ← Real.RoundLI[maxCountsPerStep*MIN[frac, initStep]];
countsPerStep ← MAX[countsPerStep, 1];
IF countsPerStep = 1 THEN
BEGIN
stepOK ← FALSE;
printMsg["Initial step is too small!"]
END
ELSE
BEGIN
stepOK ← TRUE;
printMsg["Initial step is %ld.", @countsPerStep]
END
END;
setICs: PUBLIC PROCEDURE =
BEGIN
n: nodePtr;
b: branchPtr;
v: REAL;
allNodes: BOOLEAN;
zeroVars[];
IF item # leftB THEN error[700] ELSE next[];
IF item = number THEN
BEGIN
tMin ← value;
next[];
IF item = comma THEN next[] ELSE error[702]
END;
DO
allNodes ← item = star;
IF allNodes THEN next[] ELSE [n, b] ← findNodeOrBranch[];
IF item = leftArrow THEN next[] ELSE error[703];
v ← getSignedNumber[];
IF allNodes THEN
FOR n ← nodeList, n↑.nextNode UNTIL n = NIL DO
n↑.nHist.y ← v
ENDLOOP
ELSE
IF n # NIL THEN n↑.nHist.y ← v
ELSE
IF b # NIL THEN
WITH b↑ SELECT FROM
inductor => iHist.y ← v
ENDCASE => error[721]
ELSE error[720];
IF item = comma THEN next[] ELSE EXIT
ENDLOOP;
IF item # rightB THEN error[701, TRUE] ELSE next[];
icsSet ← TRUE
END;
runIt: PUBLIC PROCEDURE RETURNS[REAL] =
BEGIN
n: CARDINAL;
IF item # leftB THEN error[700] ELSE next[];
runParms[];
IF item # rightB THEN error[701] ELSE next[];
IF tMax <= tMin THEN error[731];
IF ~AnyErrors[] THEN
BEGIN
ENABLE Failure => {error[errorNum]; GOTO quit};
IF ~icsSet THEN zeroVars[];
initPlot[tMin, tMin + tMax, TRUE];
t ← tMin;
IF printdT <= 0.0 THEN printdT ← tMax/dfltInts;
dT ← printdT;
n ← GSiteration[0.0, printAll, TRUE, TRUE];
printFromList[n, t, printIter];
plotFromList[t];
IF ~getStepSize[] THEN RETURN[tMin];
numGoodSteps ← 1;
WHILE highStepCount <= numberOfSteps DO
integration[doRungeKutta];
IF canned THEN EXIT;
IF checkPoint THEN
BEGIN
checkPoint ← FALSE;
dumpAll[t]
END
ENDLOOP;
icsSet ← FALSE;
RETURN[t]
EXITS
quit => RETURN[t];
END
ELSE RETURN[tMin]
END;
canIt: PUBLIC PROCEDURE =
BEGIN
canned ← TRUE
END;
checkIt: PUBLIC PROCEDURE =
BEGIN
checkPoint ← TRUE
END;
END.