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, FromRefText, Length, ROPE],
spGlobals USING [Aborted, argList, argNodes, branchLinkPtr, branchPtr, Canned, CheckPoint, dumpAll, error, findNodeOrBranch, getSignedNumber, Handle, initPlot, InitPP, initPrint, makeStringNB, maxCountsPerStep, maxLog, modBrPtr, modFuncPtr, next, NextStage, nodePtr, plotFromList, printFromList, PutMsgLine, RefFuncBody, RefModBody, RefC, RefI, RefL, RefR, RefV, ResetCheckPoint, ShowDetails, UpdatePP];

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= 35;
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", "tUnit", "iUnit", "vUnit", "ymin", "ymax", ""];

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=> [squish, , ]_ 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];
29=> [tUnit, , ]_ parmValue[notBoolean];
30=> [iUnit, , ]_ parmValue[notBoolean];
31=> [vUnit, , ]_ parmValue[notBoolean];
32=> [yMin, , ]_ parmValue[notBoolean];
33 => [yMax, , ]_ parmValue[notBoolean];
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];
plotFromList[handle, handle.vars.t];
UpdatePP[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];
IF handle.vars.item = string THEN {
handle.vars.title _ Rope.FromRefText[handle.vars.newString];
IF handle.vars.title = NIL THEN handle.vars.title _ "<Circuit>";
next[handle];
IF handle.vars.item = comma THEN next[handle] ELSE error[handle, 603]
};
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];
InitPP[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];
UpdatePP[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
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.
Chen, July 22, 1985 8:09:51 pm PDT, => Cedar6.0.


ξFile: spSolve.mesa
Copyright (C) 1985 by Xerox Corporation.  All rights reserved.
Last Edited by:
Sweetsun Chen, November 5, 1985 1:26:59 am PST
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]];
ΚΫ˜J™Jš
œ>™>™J™.—J˜š
Ο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˜—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š
œA˜AJš
œ$˜$Jš
œ ˜ J˜—šœœ%˜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šœœœ˜Fšœœ˜#Jš
œ<˜<Jšœœ
œ!˜@Jš
œ
˜
Jšœœœ˜EJš
œ˜—J˜Jšœœœ˜GJšœ&œ˜@šœœ˜#Jšœ&œ˜7Jšœœ˜-J˜%J˜'J˜HJ˜FJ˜ šœ˜"J˜;—J˜$Jšœ2œœ˜>J˜?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š
œ·
˜·
—J–
"mycedar" stylešœœ˜0J˜——…—dΎƒ‡