-- File: [Cherry]<Thyme>System>C03>spSolve.mesa
-- Last editted:
-- S. Chen, February 12, 1984  8:13 PM
-- 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
        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
                mf.branch^.comVal _ mf.functionProc[t, args, mf.parmVector];
              END;
            mod =>
              IF ~allInactive THEN
                BEGIN
                  oldArgs: argList;
                  oldArgs_ mf.oldArgVector;
                  FOR i: CARDINAL IN [0..LENGTH[source]) DO
                    oldArgs[i]_ source[i]^.nHist.oldy;
                  ENDLOOP;
                  res _ mf.modelResults;
                  mf.modelProc[args, oldArgs, mf.parmVector, res];
                  FOR i: CARDINAL IN [0..LENGTH[source]) DO
                    source[i]^.nHist.y_ args[i];
                  ENDLOOP;
                  mb _ mf.modelBranches;
                  UNTIL mb = NIL DO
                    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.
2/12/84:-
  original: File: [Cherry]<Thyme>System>CSIM02>spSolve.mesa
  modified to support oldArgVector