DIRECTORY
BasicTime USING [Now],
IO USING [char, int, PutF, PutFR, real, rope],
Labels USING [Set],
Real USING [Fix, FixC, FixI, Float, RoundLI, SmallestNormalizedNumber],
RealFns USING [SqRt],
RefText USING [Equal],
Rope USING [Equal, Length, ROPE],
spGlobals USING [Aborted, argList, argNodes, branchLinkPtr, branchPtr, Canned, CheckPoint, dumpAll, error, findNodeOrBranch, getSignedNumber, Handle, initPlot, initPrint, makeStringNB, maxCountsPerStep, maxLog, modBrPtr, modFuncPtr, next, NextStage, nodePtr, plotFromList, printFromList, PutMsgLine, RefFuncBody, RefModBody, RefC, RefI, RefL, RefR, RefV, ResetCheckPoint, ShowDetails];

spSolve: CEDAR PROGRAM
IMPORTS BasicTime, IO, Labels, Real, RealFns, RefText, spGlobals
EXPORTS spGlobals=

BEGIN OPEN spGlobals;

Retreat: PUBLIC SIGNAL[cause: Rope.ROPE]= CODE;
Failure: PUBLIC SIGNAL[errorNum: NAT]= CODE;
solveError: SIGNAL[s: Rope.ROPE]= CODE;
solveTest: SIGNAL[s: Rope.ROPE]= CODE;

updateFunctions: PROC[handle: Handle, f: modFuncPtr]= {
mb: modBrPtr;
source: argNodes;
args, res: argList;
allInactive: BOOL;

handle.vars.modelsSkipped_ 0;
UNTIL f=NIL DO
source_ f.arguments;
args_ f.argVector;
allInactive_ TRUE;
IF source # NIL THEN
FOR i: NAT IN [0..source.size) DO
args[i]_ source[i].nHist.y;
allInactive_ (allInactive AND source[i].marked);
ENDLOOP;
WITH f.body SELECT FROM
fcn: RefFuncBody => {
fcn.branch.comVal_ fcn.functionProc[handle.vars.t, args, f.parmVector];
};
mod: RefModBody =>
IF ~allInactive THEN {
oldArgs: argList_ mod.oldArgVector;
FOR i: NAT IN [0..source.size) DO
oldArgs[i]_ source[i].nHist.oldy;
ENDLOOP;
res_ mod.modelResults;
mod.modelProc[args, oldArgs, f.parmVector, res];
FOR i: NAT IN [0..source.size) DO
source[i].nHist.y_ args[i];
ENDLOOP;
mb_ mod.modelBranches;
UNTIL mb=NIL DO
mb.b.comVal_ res[mb.b.modelIndex];
mb_ mb.nextBranch;
ENDLOOP;
}
ELSE handle.vars.modelsSkipped_ handle.vars.modelsSkipped + 1;
ENDCASE;
f_ f.nextMFPtr;
ENDLOOP;
}; -- updateFunctions

GSiteration: PROC[handle: Handle, dT: REAL, print, firstTime, mark: BOOL_ FALSE]
RETURNS[nIter: INT_ 0]= {
E, newE, dvdt, I, newI, cond, csum, otherV, maxI, denom, tempDvDt: REAL;
n: nodePtr;
b: branchPtr;
vCur: branchPtr;
vCurBody: RefV;
bLink: branchLinkPtr;
plus, nodeOK, OK: BOOL;
realMaxIter: INT = IF handle.vars.maxIter >= LAST[INT]/2 THEN LAST[INT]
ELSE handle.vars.maxIter + handle.vars.maxIter;

updateFunctions[handle, handle.vars.intNodeModFunc];
IF print THEN printNodeEqn[handle];
UNTIL nIter > realMaxIter DO
OK_ TRUE;
IF print THEN handle.msgStream.PutF["t= %10.4f, iteration %2d --\n",
IO.real[handle.vars.t], IO.int[nIter]];
IF handle.vars.otherModFunc # NIL THEN
updateFunctions[handle, handle.vars.otherModFunc];
IF firstTime THEN updateFunctions[handle, handle.vars.intNodeModFunc];
FOR n_ handle.vars.nodeList, n.nextNode UNTIL n=NIL DO
IF n=handle.vars.gndNode THEN LOOP;
IF nIter > 1 AND n.converged THEN LOOP;
IF print THEN handle.msgStream.PutF["Node %g: ",
IO.rope[n.nodeName.name]];
bLink_ n.branches;
IF n.integrate THEN {
maxI_ csum_ I_ 0.0;
UNTIL bLink=NIL DO
b _ bLink.branch;
plus _ bLink.pos;
WITH b.body SELECT FROM
r: RefR => newI_ (bLink.otherNode.nHist.y - n.nHist.y)/b.comVal;
c: RefC => {
csum_ csum + b.comVal;
newI_ b.comVal*bLink.otherNode.nHist.f0;
};
l: RefL => newI_ (IF plus THEN -l.iHist.y ELSE l.iHist.y)/b.comVal;
v: RefV => newI_ IF plus THEN -v.vsCurrent ELSE v.vsCurrent;
i: RefI => newI_ IF plus THEN -b.comVal ELSE b.comVal
ENDCASE;
maxI_ MAX[maxI, ABS[newI]];
I_ I + newI;
bLink_ bLink.nextLink;
ENDLOOP;
E_ n.nHist.f0;
tempDvDt_ I/csum;
dvdt_ n.nHist.f0_ IF ABS[tempDvDt] < handle.vars.ignorDvDtBelow THEN 0.0 ELSE tempDvDt;
IF mark THEN n.marked_ (ABS[I] <= handle.vars.dvdtTrunc*maxI);
denom_ MAX[ABS[E], ABS[dvdt]];
nodeOK_ IF denom <= Real.SmallestNormalizedNumber THEN TRUE
ELSE (handle.vars.tol >= ABS[(dvdt - E)/denom]);
IF print THEN handle.msgStream.PutF["Old v'= %9.2f New v'= %9.2f C= %9.2f.\n",
IO.real[E], IO.real[dvdt], IO.real[csum]];
}
ELSE {
E_ n.nHist.y;
IF n.curPtr # NIL THEN {
vCur_ n.curPtr;
vCurBody_ NARROW[vCur.body];
newE_ IF vCur.posNode=n THEN vCur.negNode.nHist.y + vCur.comVal
ELSE vCur.posNode.nHist.y - vCur.comVal;
I_ vCurBody.vsCurrent;
newI_ 0.0;
UNTIL bLink=NIL DO
b _ bLink.branch;
otherV_ bLink.otherNode.nHist.y;
plus _ bLink.pos;
WITH b.body SELECT FROM
r: RefR => newI_ newI + (otherV - newE)/b.comVal;
c: RefC => newI_ newI + b.comVal*(bLink.otherNode.nHist.f0 - n.nHist.f0);
l: RefL => newI_ newI + (IF plus THEN -l.iHist.y ELSE l.iHist.y)/b.comVal;
v: RefV => newI_ IF plus THEN newI - v.vsCurrent ELSE newI + v.vsCurrent;
i: RefI => newI_ IF plus THEN newI - b.comVal ELSE newI + b.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 (handle.vars.tol >= ABS[(newI - I)/denom]);
n.nHist.f0_ 0.0;
IF NOT firstTime THEN {
tempDvDt _ handle.vars.squish*n.nHist.f1 + 
(1.0-handle.vars.squish)*(newE-n.nHist.oldy)/dT +
(IF vCur.posNode=n THEN vCur.negNode.nHist.f0 ELSE vCur.posNode.nHist.f0);
IF ABS[tempDvDt] > handle.vars.ignorDvDtBelow THEN n.nHist.f0_ tempDvDt;
};
vCurBody.vsCurrent_ newI;
IF print THEN handle.msgStream.PutF["Old I= %9.2f New I= %9.2f",
IO.real[I], IO.real[newI]];
}
ELSE {
cond_ newE_ 0.0;
UNTIL bLink=NIL DO
b _ bLink.branch;
plus _ bLink.pos;
WITH b.body SELECT FROM
r: RefR => {
cond_ cond + 1.0/b.comVal;
newE_ newE + bLink.otherNode.nHist.y/b.comVal;
};
l: RefL => newE_ newE + (IF plus THEN -l.iHist.y ELSE l.iHist.y)/b.comVal;
v: RefV => newE_ IF plus THEN newE - v.vsCurrent
ELSE newE + v.vsCurrent;
i: RefI => newE_ IF plus THEN newE - b.comVal
ELSE newE + b.comVal
ENDCASE;
bLink_ bLink.nextLink
ENDLOOP;
newE_ newE/cond;
denom_ MAX[ABS[E], ABS[newE]];
nodeOK_ IF denom <= Real.SmallestNormalizedNumber THEN TRUE
ELSE (handle.vars.tol >= ABS[(newE - E)/denom]);
};
n.nHist.y_ newE;
IF print THEN handle.msgStream.PutF["Old E= %9.2f New E= %9.2f.\n",
IO.real[E], IO.real[newE]];
};
OK_ OK AND nodeOK;
n.converged_ nodeOK;
IF ~nodeOK THEN
IF nIter = realMaxIter THEN
SIGNAL solveError[IO.PutFR["Node %g GS did not converge.",
IO.rope[spGlobals.makeStringNB[handle, n, NIL]]]];
FOR bLink_ n.branches, bLink.nextLink UNTIL bLink=NIL DO
bLink.otherNode.converged_ FALSE;
ENDLOOP;
ENDLOOP; -- n
nIter_ nIter + 1;
IF handle.vars.canned AND NOT handle.vars.forcedDump THEN SIGNAL Aborted;
IF nIter=handle.vars.maxIter THEN
handle.msgStream.PutF["t= %9.3f, iteration count exceeded",
IO.real[handle.vars.t]];
IF OK THEN EXIT;
ENDLOOP; -- nIter
}; -- GSiteration

printNodeEqn: PROC[handle: Handle]= {
n: nodePtr;
bLink: branchLinkPtr;
branch: branchPtr;
plus: BOOL;

FOR n_ handle.vars.intNodeList, n.nextIntNode UNTIL n=NIL DO
handle.msgStream.PutF["\n%g -- %13.6f %13.6f\n",
IO.rope[n.nodeName.name],
IO.real[n.nHist.y],
IO.real[n.nHist.oldy]];
FOR bLink_ n.branches, bLink.nextLink UNTIL bLink=NIL DO
branch_ bLink.branch;
plus _ bLink.pos;
WITH branch.body SELECT FROM
r: RefR => handle.msgStream.PutF["(%13.6f - %13.6f)/%13.6f",
IO.real[bLink.otherNode.nHist.y],
IO.real[n.nHist.y],
IO.real[branch.comVal] ];
c: RefC => handle.msgStream.PutF["%13.6fx%13.6f",
IO.real[branch.comVal],
IO.real[bLink.otherNode.nHist.f0] ];
l: RefL =>
IF plus THEN handle.msgStream.PutF["-1x%13.6f/%13.6f",
IO.real[l.iHist.y],
IO.real[branch.comVal]]
ELSE handle.msgStream.PutF["%13.6f/%13.6f",
IO.real[l.iHist.y],
IO.real[branch.comVal]];
v: RefV =>
IF plus THEN handle.msgStream.PutF["-1x%13.6f",
IO.real[v.vsCurrent]]
ELSE handle.msgStream.PutF["%13.6f", IO.real[v.vsCurrent]];
i: RefI =>
IF plus THEN handle.msgStream.PutF["-1x%13.6f",
IO.real[branch.comVal]]
ELSE handle.msgStream.PutF["%13.6f",IO.real[branch.comVal]];
ENDCASE;
handle.msgStream.PutF["   %s(%13.6f, %13.6f)\n",
IO.rope[bLink.otherNode.nodeName.name],
IO.real[bLink.otherNode.nHist.y],
IO.real[bLink.otherNode.nHist.f0]];
ENDLOOP;
ENDLOOP;
}; -- printNodeEqn

runParms: PROC[handle: Handle]= {
OPEN handle.vars;

parmValue: PROC[booleanOrNot: {boolean, notBoolean}]
RETURNS[r: REAL, b: BOOL, n: NAT]= {
IF booleanOrNot=boolean THEN
IF item=name THEN {
IF RefText.Equal[newString, "FALSE", FALSE] THEN b_ FALSE
ELSE
IF RefText.Equal[newString, "TRUE", FALSE] THEN b_ TRUE
ELSE error[handle, 719];
next[handle]
}
ELSE error[handle, 704]
ELSE {
r_ getSignedNumber[handle];
IF r >= 0.0 AND r < Real.Float[LAST[NAT]] THEN n_ Real.FixC[r]; -- else?
};
}; -- parmValue

nParms: NAT= 30;
parmNames: ARRAY[0..nParms) OF REF TEXT= ["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", "ignorDvDtBelow", "showRetreat", ""];

nameIndex: NAT;
tempR: REAL;
WHILE item=name DO
FOR i: NAT IN [0..nParms) DO
nameIndex_ i;
IF RefText.Equal[newString, parmNames[nameIndex], FALSE] THEN EXIT;
ENDLOOP;
next[handle];
IF item # leftArrow THEN error[handle, 703] ELSE next[handle];
SELECT nameIndex FROM
0=> {
limit: INT_ LAST[INT];
[tempR, , ]_ parmValue[notBoolean];
IF tempR >= Real.Float[limit] THEN maxIter_ limit
ELSE IF tempR > 0 THEN maxIter_ Real.Fix[tempR];
};
1=> [ , printIter, ]_ parmValue[boolean];
2=> [initStep, , ]_ parmValue[notBoolean];
3=> [initStep, , ]_ parmValue[notBoolean];
4=> [tol, , ]_ parmValue[notBoolean];
5=> [printdT, , ]_ parmValue[notBoolean];
6=> [tMax, , ]_ parmValue[notBoolean];
7=> [ , printStep, ]_ parmValue[boolean];
8=> [ , doRungeKutta, ]_ parmValue[boolean];
9=> [tMin, , ]_ parmValue[notBoolean];
10=> [ , printAll, ]_ parmValue[boolean];
11=> [ , , AMiter]_ parmValue[notBoolean];
12=> [AMup, , ]_ parmValue[notBoolean];
13=> [intTol, , ]_ parmValue[notBoolean];
14=> [AMdown, , ]_ parmValue[notBoolean];
15=> [Vmin, , ]_ parmValue[notBoolean];
16=> [Vmax, , ]_ parmValue[notBoolean];
17=> [dvdtmax, , ]_ parmValue[notBoolean];
18=> [dvdtFac, , ]_ parmValue[notBoolean];
19=> [tBreak, , ]_ parmValue[notBoolean];
20=> [RKerr, , ]_ parmValue[notBoolean];
21=> [AMerr, , ]_ parmValue[notBoolean];
22=> [ , , AMdelay]_ parmValue[notBoolean];
23=> [floor, , ]_ parmValue[notBoolean];
24=> [minRatio, , ]_ parmValue[notBoolean];
25=> [RKfac, , ]_ parmValue[notBoolean];
26=> [retRatio, , ]_ parmValue[notBoolean];
27=> [ignorDvDtBelow, , ]_ parmValue[notBoolean];
28 => [ , showRetreat, ]_ parmValue[boolean];
ENDCASE=> error[handle, 790];
IF item=comma THEN next[handle];
ENDLOOP;
}; -- runParms

predict: PROC[handle: Handle, h: REAL]= {
nodes: nodePtr_ handle.vars.intNodeList;
inds: branchPtr_ handle.vars.inductorList;
incr: REAL;

h_ h/24.0;

UNTIL nodes=NIL DO
IF ~nodes.marked THEN {
OPEN nodes.nHist;
incr_ h*(55.0*f1 - 59.0*f2 + 37.0*f3 - 9.0*f4);
y_ oldy + incr;
};
nodes_ nodes.nextIntNode;
ENDLOOP;

UNTIL inds=NIL DO
indsBody: RefL_ NARROW[inds.body];
{OPEN indsBody.iHist;
y_ oldy + h*(55.0*f1 - 59.0*f2 + 37.0*f3 - 9.0*f4)};
inds_ indsBody.nextInductor;
ENDLOOP;
}; -- predict

firstCorrector: PROC[handle: Handle, h: REAL]= {
nodes: nodePtr_ handle.vars.intNodeList;
inds: branchPtr_ handle.vars.inductorList;
f: REAL;

h_ h/24.0;

UNTIL nodes=NIL DO
IF ~nodes.marked THEN {
OPEN nodes.nHist;
y_ oldy + h*(9.0*f0 + 19.0*f1 - 5.0*f2 + f3);
oldf_ f0;
};
nodes_ nodes.nextIntNode;
ENDLOOP;

UNTIL inds=NIL DO
indsBody: RefL_ NARROW[inds.body];
{OPEN indsBody.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_ indsBody.nextInductor;
ENDLOOP;
}; -- firstCorrector

iteratingCorrector: PROC[handle: Handle, h: REAL] RETURNS[newRatio: REAL_ 4.0,
OK: BOOL_ TRUE]= {
nodes: nodePtr_ handle.vars.intNodeList;
inds: branchPtr_ handle.vars.inductorList;
f, deltay, newy, err, denom: REAL;

h_ 3.0*h/8.0;
UNTIL nodes=NIL DO
IF ~nodes.marked THEN {
OPEN nodes.nHist;
deltay_ h*(f0 - oldf);
newy_ y + deltay;
denom_ MAX[ABS[newy], ABS[y]];
IF denom > handle.vars.floor AND deltay > 1.0E-30 THEN {
err_ handle.vars.AMerr*ABS[deltay/denom];
IF handle.vars.intTol < err THEN {
OK_ FALSE;
newRatio_ MIN[newRatio, handle.vars.intTol/err];
handle.vars.worstNode_ nodes;
};
};
y_ newy;
oldf_ f0;
};
nodes_ nodes.nextIntNode;
ENDLOOP;

UNTIL inds=NIL DO
indsBody: RefL_ NARROW[inds.body];
{ OPEN indsBody.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 > handle.vars.floor AND deltay > 1.0E-30 THEN {
err_ handle.vars.AMerr*ABS[deltay/denom];
IF handle.vars.intTol < err THEN {
OK_ FALSE;
newRatio_ MIN[newRatio, handle.vars.intTol/err];
};
};
y_ newy;
oldf_ f;
};
inds_ indsBody.nextInductor;
ENDLOOP;
}; -- iteratingCorrector

AdamsMoulton: PROC[handle: Handle, dT: REAL]
RETURNS[newRatio: REAL, nGS: INT, ok: BOOL]= {
nAM: NAT;

predict[handle, dT];
nGS_ GSiteration[handle, dT, handle.vars.printAll];
firstCorrector[handle, dT];
nGS_ GSiteration[handle, dT, handle.vars.printAll] + nGS;
FOR nAM_ 2, nAM + 1 UNTIL nAM=handle.vars.AMiter DO
[]_ iteratingCorrector[handle, dT];
nGS_ GSiteration[handle, dT, handle.vars.printAll] + nGS
ENDLOOP;
[newRatio, ok]_ iteratingCorrector[handle, dT];
}; -- AdamsMoulton

RKupdate: PROC[handle: Handle, c, h: REAL, last: BOOL]
RETURNS[newRatio: REAL_ 1.0, ok: BOOL_ TRUE]= {
nodes: nodePtr_ handle.vars.intNodeList;
inds: branchPtr_ handle.vars.inductorList;
k, incr, y, err, denom: REAL;

UNTIL nodes=NIL DO
IF ~nodes.marked THEN {
nodes.nHist.sumk_ c*nodes.nHist.f0 + nodes.nHist.sumk;
IF last THEN {
incr_ h*nodes.nHist.sumk;
y_ nodes.nHist.oldy + incr;
denom_ MAX[ABS[y], ABS[nodes.nHist.y]];
IF denom > handle.vars.floor THEN {
err_ handle.vars.RKerr*ABS[(y - nodes.nHist.y)/denom];
IF handle.vars.intTol < err THEN {
ok_ FALSE;
IF handle.vars.countsPerStep=1 THEN
SIGNAL solveError[IO.PutFR["\nNode %g convergence problems.",
IO.rope[spGlobals.makeStringNB[handle, nodes, NIL]]]];
newRatio_ MIN[newRatio, handle.vars.intTol/err]
};
};
nodes.nHist.y_ y;
nodes.nHist.sumk_ 0.0;
}
ELSE {
incr_ h*nodes.nHist.f0;
nodes.nHist.y_ nodes.nHist.oldy + incr;
};
};
nodes_ nodes.nextIntNode;
ENDLOOP;

UNTIL inds=NIL DO
indsBody: RefL_ NARROW[inds.body];
k_ inds.posNode.nHist.y - inds.negNode.nHist.y;
indsBody.iHist.sumk_ c*k + indsBody.iHist.sumk;
IF last THEN {
indsBody.iHist.y_ indsBody.iHist.oldy + h*indsBody.iHist.sumk;
indsBody.iHist.sumk_ 0.0
}
ELSE indsBody.iHist.y_ indsBody.iHist.oldy + h*k;
inds_ indsBody.nextInductor;
ENDLOOP;
}; -- RKupdate

RungeKutta: PROC[handle: Handle, dT: REAL]
RETURNS[newRatio: REAL, nGS: NAT, ok: BOOL]= {
time: REAL;

time_ handle.vars.t;
handle.vars.t_ handle.vars.t - 0.5*dT;
[] _ RKupdate[handle, 1.0, 0.5*dT, FALSE];
nGS_ GSiteration[handle, 0.5*dT, handle.vars.printAll];
[] _ RKupdate[handle, 2.0, 0.5*dT, FALSE];
nGS_ GSiteration[handle, 0.5*dT, handle.vars.printAll] + nGS;
handle.vars.t_ time;
[] _ RKupdate[handle, 2.0, dT, FALSE];
nGS_ GSiteration[handle, dT, handle.vars.printAll] + nGS;
[newRatio, ok]_ RKupdate[handle, 1.0, 0.1666667*dT, TRUE]
}; -- RungeKutta

changeStepSize: PROC[handle: Handle, ratio: REAL, RK: BOOL]= {
int: CHAR;

IF ratio < 1.0 OR handle.vars.countsPerStep < maxCountsPerStep THEN {
handle.vars.numGoodSteps_ 1;

handle.vars.worstNodeLog[handle.vars.curLog]_
IF handle.vars.worstNode # NIL THEN 
[visited: FALSE,
node: handle.vars.worstNode,
t: handle.vars.t,
v: handle.vars.worstNode.nHist.y,
dvdt: handle.vars.worstNode.nHist.f0]
ELSE [FALSE, NIL, 0.0, 0.0, 0.0];

handle.vars.curLog_ (handle.vars.curLog + 1) MOD maxLog;

IF handle.vars.countsPerStep=1 THEN SIGNAL solveError["Step too small"];

handle.vars.countsPerStep_
IF ratio > Real.Float[maxCountsPerStep/handle.vars.countsPerStep] THEN
maxCountsPerStep
ELSE MAX[1, Real.RoundLI[handle.vars.countsPerStep*ratio]];

int_ IF RK THEN 'R ELSE 'A;
IF ShowDetails[handle] THEN {
Labels.Set[handle.time,
IO.PutFR["%12.6e", IO.real[handle.vars.t]]];
Labels.Set[handle.step,
IO.PutFR["%10d %g", IO.int[handle.vars.countsPerStep], IO.char[int]]];
};
};
}; -- changeStepSize

saveState: PROC[handle: Handle]= {
inds: branchPtr_ handle.vars.inductorList;
volts: branchPtr_ handle.vars.vSourceList;
indsBody: RefL;
voltsBody: RefV;

FOR nodes: nodePtr_ handle.vars.nodeList, nodes.nextNode UNTIL nodes=NIL DO
OPEN nodes.nHist;

SolveTrouble: PROC [explanation: Rope.ROPE] = {
SIGNAL solveTest[IO.PutFR["Node %g %g", IO.rope[spGlobals.makeStringNB[handle, nodes, NIL]], IO.rope[explanation]]];
};

oldy_ y;
f4_ f3; f3_ f2; f2_ f1; f1_ f0;
SELECT TRUE FROM
y < handle.vars.Vmin => SolveTrouble["below Vmin"];
y > handle.vars.Vmax => SolveTrouble["above Vmax"];
ABS[f0] > handle.vars.dvdtmax => SolveTrouble["change rate too large"];
ENDCASE => NULL;
ENDLOOP;

UNTIL inds=NIL DO
indsBody_ NARROW[inds.body];
{ OPEN indsBody.iHist;
oldy_ y;
f4_ f3; f3_ f2; f2_ f1;
f1_ inds.posNode.nHist.y - inds.negNode.nHist.y;
};
inds_ indsBody.nextInductor;
ENDLOOP;

UNTIL volts=NIL DO
voltsBody_ NARROW[volts.body];
voltsBody.oldCurrent_ voltsBody.vsCurrent;
volts_ voltsBody.nextvSource;
ENDLOOP;
}; -- saveState

restoreState: PROC[handle: Handle]= {
nodes: nodePtr_ handle.vars.nodeList;
inds: branchPtr_ handle.vars.inductorList;
volts: branchPtr_ handle.vars.vSourceList;
indsBody: RefL;
voltsBody: RefV;

UNTIL nodes=NIL DO
nodes.nHist.y_ nodes.nHist.oldy;
nodes.nHist.f0_ nodes.nHist.f1;
IF ABS[nodes.nHist.f0] < handle.vars.ignorDvDtBelow THEN nodes.nHist.f0_ 0.0;
nodes_ nodes.nextNode;
ENDLOOP;

UNTIL inds=NIL DO
indsBody_ NARROW[inds.body];
indsBody.iHist.y_ indsBody.iHist.oldy;
inds_ indsBody.nextInductor;
ENDLOOP;

UNTIL volts=NIL DO
voltsBody_ NARROW[volts.body];
voltsBody.vsCurrent_ voltsBody.oldCurrent;
volts_ voltsBody.nextvSource;
ENDLOOP;
}; -- restoreState

integration: PROC[handle: Handle, RK: BOOL]= {
nGS: INT;
out, ok: BOOL_ FALSE; 
ratio, oldTime, DT: REAL;
oldLowStepcount, oldHighStepCount: INT;

saveState[handle];
oldTime_ handle.vars.t;
oldLowStepcount_ handle.vars.lowStepCount;
oldHighStepCount_ handle.vars.highStepCount;
DO
handle.vars.lowStepCount_ handle.vars.lowStepCount + handle.vars.countsPerStep;
IF handle.vars.lowStepCount > maxCountsPerStep THEN {
handle.vars.highStepCount_ handle.vars.highStepCount + 1;
handle.vars.lowStepCount_ handle.vars.lowStepCount - maxCountsPerStep;
out_ TRUE;
}
ELSE out_ handle.vars.printStep;

handle.vars.t_ handle.vars.tMin + handle.vars.dT*handle.vars.lowStepCount +
(handle.vars.dT*maxCountsPerStep)*handle.vars.highStepCount;
DT_ handle.vars.dT*handle.vars.countsPerStep;

{ ENABLE Retreat=> {
IF ~Rope.Equal[cause, handle.vars.lastRetreatCause] THEN {
IF Rope.Length[cause] > 0 AND handle.vars.showRetreat THEN
handle.msgStream.PutF["%g at t= %12.6e.\n",
IO.rope[cause], IO.real[handle.vars.t]];
handle.vars.lastRetreatCause_ cause;
IF handle.vars.countsPerStep = 1 THEN REJECT;
};
ok_ FALSE;
ratio_ handle.vars.retRatio;
GOTO FlushIt;
};

IF handle.vars.numGoodSteps < 4 THEN {
[ratio, nGS, ok]_ RungeKutta[handle, DT];
IF ok THEN nGS_ nGS +
GSiteration[handle, DT, handle.vars.printAll,, TRUE];
ratio_ RealFns.SqRt[ratio]*handle.vars.RKfac;
IF ok AND ~RK THEN
handle.vars.numGoodSteps_ 1 + handle.vars.numGoodSteps;
}
ELSE {
handle.vars.AMcount_ 1 + handle.vars.AMcount;
[ratio, nGS, ok]_ AdamsMoulton[handle, DT];
IF ok THEN {
nGS_ nGS + GSiteration[handle, DT, handle.vars.printAll,, TRUE];
IF ratio > handle.vars.AMup AND
handle.vars.AMcount >= handle.vars.AMdelay THEN {
ratio_ RealFns.SqRt[RealFns.SqRt[ratio]];
changeStepSize[handle, ratio, FALSE];
handle.vars.AMcount_ 0;
};
}
ELSE ratio_ RealFns.SqRt[ratio];
};
EXITS
FlushIt=> NULL;
};
IF ok THEN EXIT;
restoreState[handle];
handle.vars.t_ oldTime;
handle.vars.highStepCount_ oldHighStepCount;
handle.vars.lowStepCount_ oldLowStepcount;
changeStepSize[handle, MIN[ratio, handle.vars.minRatio],
(handle.vars.numGoodSteps < 4)];
handle.vars.AMcount_ 0
ENDLOOP;

NextStage[handle];
IF out THEN printFromList[handle, nGS, handle.vars.t, handle.vars.printIter];
IF out THEN plotFromList[handle, handle.vars.t];
WHILE handle.vars.tBreak > 0.0 AND handle.vars.t >= handle.vars.tBreak DO
oldTime_ handle.vars.t ENDLOOP;
}; -- integration

zeroVars: PROC[handle: Handle]= {
n: nodePtr_ handle.vars.nodeList;
inds: branchPtr_ handle.vars.inductorList;
vSs: branchPtr_ handle.vars.vSourceList;
indsBody: RefL;
vSsBody: RefV;

handle.vars.gndNode.marked_ TRUE;

UNTIL n=NIL DO
n.nHist_ [];
n.marked_ FALSE;
n_ n.nextNode;
ENDLOOP;

UNTIL inds=NIL DO
indsBody_ NARROW[inds.body];
indsBody.iHist_ [];
inds_ indsBody.nextInductor;
ENDLOOP;

UNTIL vSs=NIL DO
vSsBody_ NARROW[vSs.body];
vSsBody.vsCurrent_ vSsBody.oldCurrent_ 0.0;
IF vSs.controller=NIL THEN
IF vSs.posNode=handle.vars.gndNode THEN vSs.negNode.marked_ TRUE
ELSE
IF vSs.negNode=handle.vars.gndNode THEN vSs.posNode.marked_ TRUE;
vSs_ vSsBody.nextvSource;
ENDLOOP;
}; -- zeroVars

getStepSize: PROC[handle: Handle] RETURNS[stepOK: BOOL]= {
vMax, dvdtMax, frac: REAL_ 0.0;
inds: branchPtr_ handle.vars.inductorList;
nodes: nodePtr_ handle.vars.nodeList;

handle.vars.lowStepCount_ 0;
handle.vars.highStepCount_ 0;
handle.vars.numberOfSteps_ Real.FixI[handle.vars.tMax/handle.vars.printdT];
handle.vars.worstNode_ NIL;

UNTIL nodes=NIL DO
vMax_ MAX[ABS[nodes.nHist.y], vMax];
IF ABS[nodes.nHist.f0] > dvdtMax THEN {
dvdtMax_ ABS[nodes.nHist.f0];
handle.vars.worstNode_ nodes;
};
nodes_ nodes.nextNode;
ENDLOOP;

frac_ IF dvdtMax > 0.0 AND vMax > 0.0 THEN
handle.vars.intTol*vMax/(dvdtMax*handle.vars.printdT)
ELSE handle.vars.initStep;
handle.vars.dT_ handle.vars.printdT/maxCountsPerStep;
handle.vars.dvdtTrunc_ handle.vars.dvdtFac*handle.vars.intTol;
handle.vars.countsPerStep_
Real.RoundLI[maxCountsPerStep*MIN[frac, handle.vars.initStep]];
handle.vars.countsPerStep_
MAX[handle.vars.countsPerStep, 1];
IF handle.vars.countsPerStep=1 THEN {
stepOK_ FALSE;
PutMsgLine[handle, "Initial step is too small!"];
}
ELSE {
stepOK_ TRUE;
};
}; -- getStepSize

setICs: PUBLIC PROC[handle: Handle]= {
n: nodePtr;
b: branchPtr;
v: REAL;
allNodes: BOOL;

zeroVars[handle];
IF handle.vars.item # leftB THEN error[handle, 700] ELSE next[handle];
IF handle.vars.item=number THEN {
handle.vars.tMin_ handle.vars.value;
next[handle];
IF handle.vars.item=comma THEN next[handle] ELSE error[handle, 702];
};
DO
allNodes_ (handle.vars.item=star);
IF allNodes THEN next[handle] ELSE [n, b]_ findNodeOrBranch[handle];
IF handle.vars.item=leftArrow THEN next[handle] ELSE error[handle, 703];
v_ getSignedNumber[handle];
IF allNodes THEN
FOR n_ handle.vars.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.body SELECT FROM
l: RefL => l.iHist.y_ v;
ENDCASE=> error[handle, 721]
ELSE error[handle, 720];
IF handle.vars.item=comma THEN next[handle] ELSE EXIT;
ENDLOOP;
IF handle.vars.item # rightB THEN error[handle, 701, TRUE] ELSE next[handle];
handle.vars.icsSet_ TRUE;
}; -- setICs

runIt: PUBLIC PROC[handle: Handle] RETURNS[REAL]= {
n: INT;

IF handle.vars.item # leftB THEN error[handle, 700] ELSE next[handle];
runParms[handle];
IF handle.vars.item # rightB THEN error[handle, 701] ELSE next[handle];
IF handle.vars.tMax <= handle.vars.tMin THEN error[handle, 731];
IF handle.vars.errCount <= 0 THEN {
ENABLE Failure => {error[handle, errorNum]; GOTO quit};
IF ~handle.vars.icsSet THEN zeroVars[handle];
handle.vars.simTime_ BasicTime.Now[];
initPrint[handle, handle.vars.simTime];
initPlot[handle, handle.vars.tMin, handle.vars.tMin + handle.vars.tMax];
handle.vars.t_ handle.vars.tMin;
IF handle.vars.printdT <= 0.0 THEN
handle.vars.printdT_ handle.vars.tMax/handle.vars.dfltInts;
handle.vars.dT_ handle.vars.printdT;
n_ GSiteration[handle, 0.0, handle.vars.printAll, TRUE, TRUE];
printFromList[handle, n, handle.vars.t, handle.vars.printIter];
plotFromList[handle, handle.vars.t];
IF ~getStepSize[handle] THEN RETURN[handle.vars.tMin];
handle.vars.numGoodSteps_ 1;
WHILE handle.vars.highStepCount <= handle.vars.numberOfSteps DO
integration[handle, handle.vars.doRungeKutta];
IF Canned[handle] THEN EXIT;
IF CheckPoint[handle] THEN {
ResetCheckPoint[handle];
dumpAll[handle, handle.vars.t];
};
ENDLOOP;
handle.vars.icsSet_ FALSE;
RETURN[handle.vars.t];
EXITS
quit=> RETURN[handle.vars.t];
}
ELSE RETURN[handle.vars.tMin]
}; -- runIt

END.

CHANGE LOG
$File: [Cherry]<Thyme>Cedar5.2>System>spSolve.mesa
Last Edited by: SChen, May 28, 1985 6:54:27 pm PDT
since ~allInactive, from above we know that source is not Nil
n.nHist.f0_ IF firstTime THEN 0.0 ELSE
handle.vars.squish*n.nHist.f1 + 
(1.0-handle.vars.squish)*(newE-n.nHist.oldy)/dT +
(IF vCur.posNode=n THEN vCur.negNode.nHist.f0
ELSE vCur.posNode.nHist.f0);
handle.msgStream.PutF["Initial step is %d.\n",
IO.int[handle.vars.countsPerStep]];
Wilhelm, April 27, 1982  4:04 PM
Barth, 7-May-82 10:32:52 PDT
Chen, April 22, 1983  10:56 AM. Changed dvdtmax from 1E13 to 1E15.
Barth, July 11, 1983  1:42 PM
Chen, February 12, 1984  8:13 PM modified to support oldArgVector.
Chen, June 12, 1984 6:10:21 pm PDT, cedarized.
Barth, February 8, 1985 3:17:34 pm PST
McCreight, April 4, 1985 1:57:23 pm PST
Chen, May 9, 1985 6:40:28 pm PDT:
1. RealOps.SqRt => RealFns.SqRt; 
2. moved lastRetreatCause from a global variable to be one on the handle.vars.
3. in order to get around the concurrency problem in floating point software, added a new run-parameter ignorDvDtBelow to ignor dvdt when its magnitude is smaller than ignorDvDtBelow.
4. added another new run-parameter, showRetreat, so that if one really wants to see the reasons of retreats during "integration", he may do so by setting it true in the run statement.
�Ê:˜J™1J™2J˜š
Ïk	˜	Jšœ
�œ˜Jš�œ
�œ&˜.Jšœ�œ˜Jšœ�œ=˜GJšœ�œ˜Jšœ�œ	˜Jšœ�œ�œ˜!Jšœ
�œò˜�J˜—šœ	�œ
�˜Jš�œ�œ+˜@Jš�œ˜J˜—Jš�œ
�œ˜J˜Jš	œ	�œ
�œ
�œ�œ
˜/Jš	œ	�œ
�œ�œ�œ
˜,Jšœ�œ	�œ�œ
˜'Jšœ�œ	�œ�œ
˜&J˜šœ�œ"˜7J˜
J˜J˜Jšœ
�œ
˜J˜J˜š�œ�œ
�˜J˜J˜Jšœ
�œ
˜š�œ
�œ
�˜š�œ�œ
�œ�˜!J˜Jšœ�œ˜0Jš�œ
˜——š�œ�œ
�˜˜J˜GJ˜—˜š�œ�œ˜Jš
œ=™=J˜#š�œ�œ
�œ�˜!J˜!—Jš�œ
˜J˜J˜0š�œ�œ
�œ�˜!J˜—Jš�œ
˜J˜š�œ�œ
�˜J˜"J˜—Jš�œ
˜J˜
—Jš�œ:˜>—Jš�œ
˜—J˜Jš�œ
˜—JšœÏc˜J˜—š
Ïnœ�œ�œ�œ�œ
˜PJš�œ�œ˜Jš�
œ�
œ3�œ
˜HJ˜J˜
J˜J˜J˜Jšœ�œ�œ
˜šœ
�œ�œ�œ
�œ�œ
�œ
�œ
˜GJš�œ+˜/—J˜J˜4Jš�œ�œ˜#š�œ�˜Jš�œ�œ
˜	š�œ�œ7˜DJš�œ�œ
˜'—š�œ�œ
�˜&J˜2—Jš�œ�œ5˜Fš�œ%�œ�œ
�˜6Jš�œ�œ
�œ
˜#Jš�œ�œ
�œ
�œ
˜'š�œ�œ#˜0Jš�œ˜—J˜š�œ
�œ˜Jšœ�
œ˜š�œ�œ
�˜J˜J˜š�œ�œ
�˜J˜@˜J˜J˜(J˜—Jšœ�œ�œ�œ˜CJšœ�œ�œ�œ
˜<Jšœ�œ�œ�œ	˜5Jš�œ
˜—Jšœ�œ�œ˜Jš�
œ�
œ˜J˜Jš�œ
˜—Jš�
œ
˜Jšœ
�
œ˜Jš	œ�œ
�œ(�œ�œ
˜WJš�œ�œ�œ
�
œ!˜>Jš	œ�œ
�œ
�
œ�œ˜šœ�œ(�œ
�˜;Jš�œ�œ	�
œ
˜0—š�œ�œA˜NJš�œ�
œ�œ
�œ
˜*—J˜
—š�œ˜Jš�
œ˜
š�œ�œ
�œ˜J˜Jšœ
�œ˜šœ�œ�œ#˜?Jš�œ$˜(—Jš�
œ˜J˜
š�œ�œ
�˜J˜J˜ J˜š�œ�œ
�˜J˜1J˜IJšœ�œ�œ�œ˜JJšœ�œ�œ�œ˜IJšœ�œ�œ�œ˜BJš�œ
˜—J˜Jš�œ
˜—Jšœ�œ�œ
�
œ�œ
�
œ˜4Jš	œ�œ
�œ
�
œ�œ˜šœ�œ(�œ
�˜;Jš�œ�œ	�
œ
˜0—šœ�œ�œ�™&™ J™1šœ
�œ�œ™-Jš�œ™———Jš
œ˜š�œ
�œ�œ˜˜+J˜1Jšœ
�œ�œ�œ˜J—Jš�œ
�œ(�œ˜HJ˜—J˜š�œ�œ3˜@Jš�œ�
œ�œ
˜—J˜
—š�œ˜J˜š�œ�œ
�˜J˜J˜š�œ�œ
�˜˜J˜J˜.J˜—Jšœ�œ�œ�œ˜Jšœ�œ�œ˜0Jš�œ˜—šœ�œ�œ˜-Jš�œ˜—Jš�œ
˜—J˜Jš�œ
˜—J˜Jš	œ�œ
�œ
�
œ�œ˜šœ�œ(�œ
�˜;Jš�œ�œ	�
œ
˜0—J˜—J˜š�œ�œ6˜CJš�œ�
œ�œ
˜—J˜—Jš�œ�œ
�œ˜J˜š�œ	�˜š�œ�˜š�œ�œ&˜:Jš�œ(�œ˜2——š�œ#�œ�œ
�˜8Jšœ�œ
˜!Jš�œ
˜——Jš�œž˜
—J˜Jš
�œ�œ
�œ�œ
�œ	˜Iš�œ�˜!š
œ;˜;Jš�œ˜——Jš�œ
�œ
�œ
�œ
˜Jš�œž˜—Jšœž˜J˜—šœ�œ˜%J˜J˜J˜Jšœ�œ
˜J˜š�œ+�œ�œ
�˜<š
œ0˜0Jš�œ˜Jš�œ˜Jš�œ˜—š�œ#�œ�œ
�˜8J˜J˜š�œ
�œ
�˜š
œ<˜<Jš�œ˜!Jš�œ˜Jš�œ˜—š
œ1˜1Jš�œ˜Jš�œ"˜$—˜
š�œ�œ*˜6Jš�œ˜Jš�œ˜—š�œ'˜+Jš�œ˜Jš�œ˜——˜
š�œ�œ#˜/Jš�œ˜—Jš�œ!�œ˜;—˜
š�œ�œ#˜/Jš�œ˜—Jš�œ �œ˜<—Jš�œ
˜—š
œ0˜0Jš�œ%˜'Jš�œ˜!Jš�œ!˜#—Jš�œ
˜—Jš�œ
˜—Jšœž˜J˜—šœ
�œ˜!Jš�œ
˜J˜šœ�œ%˜4Jš�œ�œ�œ�œ˜$š�œ�˜š�œ�œ˜Jš�œ#�œ�œ�˜9š
�˜Jš�œ"�œ�œ�˜7Jš�œ˜—J˜J˜
—Jš�œ˜—š�œ˜J˜Jš�œ
�œ�œ
�œ�œž	˜HJ˜—Jšœž˜J˜—Jšœ�œ˜Jš	œ�œ�œ
�œ
�œ·˜ÞJ˜Jšœ�œ
˜Jšœ�œ
˜š�œ�˜š�œ�œ
�œ
�˜J˜
Jš�œ0�œ�œ
�œ
˜CJš�œ
˜—J˜
Jš�œ�œ�œ˜>š�œ�˜˜Jšœ�œ�œ
�œ˜J˜#Jš�œ�œ˜1Jš�œ
�œ�œ˜0J˜—J˜)J˜*J˜*J˜%J˜)J˜&J˜)J˜,J˜&J˜)J˜*J˜'J˜)J˜)J˜'J˜'J˜*J˜*J˜)J˜(J˜(J˜+J˜(J˜+J˜(J˜+J˜1Jš
œ-˜-Jš�œ˜—Jš�œ�œ˜ Jš�œ
˜—Jšœž˜J˜—šœ	�œ�œ˜)J˜(J˜*Jšœ�œ
˜J˜J˜
J˜š�œ�œ
�˜š�œ�œ˜Jš�œ
˜J˜/J˜J˜—J˜Jš�œ
˜—J˜š�œ�œ
�˜Jšœ�œ˜"šœ
�œ˜J˜4—J˜Jš�œ
˜—Jšœž
˜
J˜—šœ�œ�œ˜0J˜(J˜*Jšœ�œ
˜J˜J˜
J˜š�œ�œ
�˜š�œ�œ˜Jš�œ
˜J˜-J˜	J˜—J˜Jš�œ
˜—J˜š�œ�œ
�˜Jšœ�œ˜"šœ
�œ˜J˜/J˜,J˜	—J˜Jš�œ
˜—Jšœž˜J˜—š	œ�œ�œ�œ�œ˜NJš�œ�œ�œ˜J˜(J˜*Jšœ�œ
˜"J˜J˜
š�œ�œ
�˜š�œ�œ˜Jš�œ
˜J˜J˜Jšœ�œ
�œ�œ˜š�œ�œ�œ˜8Jšœ�œ˜)š�œ�œ˜"Jš�œ�œ
˜
Jšœ
�œ#˜0J˜J˜—J˜—J˜J˜	J˜—J˜Jš�œ
˜—J˜š�œ�œ
�˜Jšœ�œ˜"šœ�œ˜J˜/J˜J˜Jšœ�œ
�œ�œ˜š�œ�œ�œ˜8Jšœ�œ˜)š�œ�œ˜"Jš�œ�œ
˜
Jšœ
�œ#˜0J˜—J˜—J˜J˜J˜—J˜Jš�œ
˜—Jšœž˜J˜—šŸœ�œ�œ
˜,Jš�œ�œ�œ�œ˜.Jšœ�œ
˜	J˜J˜J˜3J˜J˜9š�œ�œ�˜3J˜#J˜8Jš�œ
˜—J˜/Jšœž˜J˜—šŸœ�œ�œ�œ
˜6Jš�œ�œ�œ�œ˜/J˜(J˜*Jšœ�œ
˜J˜š�œ�œ
�˜š�œ�œ˜J˜6š�œ�œ˜J˜J˜Jšœ�œ
�œ�œ˜'š�œ�œ˜#Jšœ�œ˜6š�œ�œ˜"Jšœ�œ
˜
š�œ�˜#š�œ�œ)˜=Jš�œ,�œ˜6——Jšœ
�œ"˜/J˜—J˜—J˜J˜J˜
—š�œ˜J˜J˜'J˜—J˜—J˜Jš�œ
˜—J˜š�œ�œ
�˜Jšœ�œ˜"J˜/J˜/š�œ�œ˜J˜>J˜J˜
—Jš�œ-˜1J˜Jš�œ
˜—Jšœž˜J˜—šŸ
œ�œ�œ
˜*Jš�œ�œ�œ�œ˜.Jšœ�œ
˜J˜J˜J˜&Jšœ#�œ˜*J˜7Jšœ#�œ˜*J˜=J˜Jšœ�œ˜&J˜9Jšœ4�œ
˜9Jšœž
˜J˜—š	œ�œ�œ�œ�œ˜>Jšœ�œ
˜
J˜š�œ
�œ.�œ˜EJ˜J˜˜-š�œ�œ
�œ
˜$Jšœ
�œx˜‡
—Jš�œ�œ
�œ˜!J˜—Jšœ-�œ˜8J˜Jš�œ�œ
�œ˜HJ˜˜š�œ@�˜FJ˜—Jš�œ
�œ3˜;J˜—Jš	œ�œ
�œ
�œ�œ˜š�œ�œ˜˜Jš�œ�œ˜,—˜Jš�œ�œ!�œ
˜F—J˜—J˜—Jšœž˜J˜—šœ�œ˜"J˜*J˜*J˜J˜J˜š�œ6�œ�œ
�˜KJš�œ
˜J˜šŸœ�œ�œ˜/Jš
�œ�œ�œ,�œ�œ˜tJ˜J˜—J˜J˜š�œ
�œ
�˜Jš
œ3˜3Jš
œ3˜3Jš�œD˜GJš�œ�œ
˜—Jš�œ
˜—J˜š�œ�œ
�˜Jšœ
�œ˜šœ�œ˜J˜J˜J˜0J˜—J˜Jš�œ
˜—J˜š�œ�œ
�˜Jšœ�œ
˜J˜*J˜Jš�œ
˜—Jšœž˜J˜—šœ�œ˜%J˜%J˜*J˜*J˜J˜J˜š�œ�œ
�˜J˜ J˜Jš�œ
�œ.�œ˜MJ˜Jš�œ
˜—J˜š�œ�œ
�˜Jšœ
�œ˜J˜&J˜Jš�œ
˜—J˜š�œ�œ
�˜Jšœ�œ
˜J˜*J˜Jš�œ
˜—Jšœž˜J˜—šœ
�œ�œ�œ˜.Jšœ�œ
˜	Jšœ	�œ˜Jšœ�œ�œ
˜Jšœ#�œ
˜'J˜J˜J˜J˜*J˜,š
�˜J˜Oš�œ-�œ˜5J˜9J˜FJšœ�œ
˜
J˜
—Jš�œ˜ J˜˜KJ˜<—Jš�œ+˜-J˜šœ�œ˜š�œ2�œ˜:š�œ�œ�˜:˜+Jš�œ�œ˜(——Jš
œ$˜$Jš�œ�œ
�œ
˜-Jš
œ˜—Jšœ�œ
˜
J˜Jš�œ	˜
J˜J˜š�œ�œ˜&Jšœ%�œ˜)š�œ�œ˜Jšœ�œ�œ˜5—J˜-š�œ�œ�œ
�˜J˜7—J˜
—š�œ˜J˜-Jšœ'�œ˜+š�œ�œ˜Jšœ�œ�œ˜@š�œ�˜Jšœ+�œ˜1J˜)Jšœ�œ˜%J˜J˜—J˜
—Jš�œ˜ J˜——š
�˜Jšœ
�œ
˜—J˜Jš�œ�œ
�œ
˜J˜J˜J˜,J˜*šœ�œ˜8J˜ —J˜—Jš�œ
˜J˜J˜Jš�œ�œB˜MJš�œ�œ%˜0š�œ�œ%�˜IJšœ�œ
˜—Jšœž˜J˜—šœ
�œ˜!J˜!J˜*J˜(J˜J˜J˜Jšœ�œ
˜!J˜š�œ�œ
�˜J˜Jšœ
�œ
˜J˜Jš�œ
˜—J˜š�œ�œ
�˜Jšœ
�œ˜J˜J˜Jš�œ
˜—J˜š�œ�œ
�˜Jšœ	�œ˜J˜+š�œ�œ
�˜Jš�œ!�œ�˜@š
�˜Jš�œ!�œ�œ
˜A——J˜Jš�œ
˜—Jšœž˜J˜—šœ
�œ�œ	�œ˜:Jšœ�œ˜J˜*J˜%J˜J˜J˜J˜KJšœ�œ
˜J˜š�œ�œ
�˜Jšœ�œ
�œ˜$š�œ
�œ�œ˜'Jšœ	�œ˜J˜J˜—J˜Jš�œ
˜—J˜šœ�œ�œ�˜*J˜5Jš�œ˜—J˜5J˜>˜Jšœ�œ˜?—˜Jš�œ˜"—š�œ�œ˜%Jšœ�œ
˜J˜1J˜
—š�œ˜Jšœ�œ
˜
š
œ.™.Jš
œ#™#—J˜—Jšœž˜J˜—šœ�œ
�œ˜&J˜J˜
Jšœ�œ
˜Jšœ
�œ
˜J˜J˜Jš�œ�œ�œ˜Fš�œ�œ˜!J˜$J˜
Jš�œ�œ�œ˜DJ˜—š
�˜J˜"Jš�œ
�œ�œ"˜DJš�œ�œ�œ˜HJ˜š�œ
�˜š�œ%�œ�œ
�˜6J˜
Jš
�˜——š
�˜Jš�œ�œ
�œ
˜š
�˜š�œ�œ
�˜š�œ�œ
�˜J˜Jš�œ˜——Jš�œ˜——Jš�œ�œ�œ
�œ
˜6Jš�œ
˜—Jš�œ�œ�œ�œ˜MJšœ�œ
˜Jšœž	˜J˜—š	œ�œ
�œ�œ
�œ˜3Jšœ�œ
˜J˜Jš�œ�œ�œ˜FJ˜Jš�œ�œ�œ˜GJš�œ&�œ˜@š�œ�œ˜#Jš�œ&�œ˜7Jš�œ�œ˜-J˜%J˜'J˜HJ˜ š�œ�˜"J˜;—J˜$Jšœ2�œ�œ˜>J˜?J˜$Jš�œ�œ
�œ˜6J˜š�œ8�˜?J˜.Jš�œ�œ
�œ
˜š�œ�œ˜J˜J˜J˜—Jš�œ
˜—Jšœ�œ
˜Jš�œ˜š
�˜Jšœ�œ˜—J˜
—Jš�œ
�œ˜Jšœž˜J˜—Jš�œ
˜J˜š�œ
�˜
J™ Jš
œ™J™BJ™J™BJ™.J™&J™'™!J™!Jš
œN™NJ™·
Jš
œ·
™·
———…—^šø