INTERFACE Integrator;

(*
  This interface defines the abstract data type "Integrator.T",
  an object that numerically solves ordinary differential equations
  of the form 
  |
  |    ds/dt = F(t, s)
  |
  where "s" is an unknown function from time "t" to some 
  'state space' "R^n", and "F" is a client-provided function.
  
  Created 95/06 by R.L.Liesenfeld, C.H.Q.Forster, and J.Stolfi
*)

TYPE
  Time = LONGREAL;            (* A time value *)
  Coord = LONGREAL;           (* A state coordinate *)
  Speed = LONGREAL;           (* Derivative of "Coord" relative to "Time" *)
  State = ARRAY OF Coord;     (* The system's state. *)
  Error = ARRAY OF Coord;     (* Error estimate for a "State". *)
  Velocity = ARRAY OF Speed;  (* Derivative of "State" relative to "Time". *)

  T = OBJECT METHODS

      setSize(n: CARDINAL): T;
        (*
          Allocates internal working storage for a system with "n"
          state variables. *)
          
      integrate(
          t: Time;
          VAR s: State;
          diff: DiffProc;
          dt: Time;
          checkStep: CheckStepProc;
          checkState: CheckStateProc := NIL;
        ): Time;
        (*
          Integrates the equation "ds/dt = diff(t, s)" from time "t"
          and state "s", until a limit set by the "checkState" and
          "checkStep" procedures below.
          
          Upon entry, the vector "s" should contain the initial state
          at time "t". Upon exit, "s" will contain the final state.
          The corresponding time is returned as result of the call.
          
          Integration will proceed by a sequence of chained steps.
          The length of the first step will be "dt"; that of later
          steps is set by the client during the integration, 
          
          
          The integrator will call "diff(t, s, v)" several times per
          step to compute the derivative "v = ds/dt" for a state "s"
          and time "t".  Beware that those times "t" are not always
          increasing, and that the states "s" may be widely off the
          true path.
          
          If the "checkState" procedure is not NIL, the integrator
          will call "checkState(t, s, v)", right after "diff(t, s,
          v)", for every valid state in the path (not including the
          'tentative' states mentioned above).
          
          The "checkState" procedure is expected to return the
          approximate upper time bount "tf" for stopping the
          integration.  If "tf" is equal to "t" (or earlier), the
          integration will stop right there.  If "tf" is later than
          "t", the next step will be taken and will not extend beyond
          "tf".  (If "checkState" is NIL, the integration will 
          continue anyway.)
          
          Note that "checkState" is always applied to the first state
          on the path; and is not applied to the final state, except
          when "checkState" itself terminates the integration by
          returning "t" as result.
          
          If the "checkStep" procedure is not NIL, the integrator will
          call "checkStep(ta, sa, va, tb, sb, er)" at every step to
          decide whether the integration should be truncated
          prematurely.  The call says that the last step went from
          time "ta" and state "sa" to time "tb" and state "sb"; that
          the derivative at "ta" was "va"; and that the estimated
          error vector at "tb" is "er".  The "checkStep" procedure may
          then return
            
            (1) "ta"  (or earlier) to stop the integration at "ta"; or
            
            (2) "tb"  to stop the integration at "tb"; or
            
            (3) any value "tc" in the interval "(ta __ tb)", to
                redo the integration from "ta" to "t";
                 
            (4) any value "tc" greater than "tb", to continue the 
                integration beyond that point.
                 
          In case (3), after the integration from "ta" to "tc"
          is redone, the "checkStep" procedure will be called again
          to inspect the truncated step, and perhaps refine the estimate
          of the stopping time "tc".  The client should make sure that
          this terminates after a finite numebr of steps.
          
          The "checkStep" procedure is urged to use the "computeStep"
          method below to compute the appropriate 
          
          In cases (2) and (4), the "checkStep" procedure is allowed
          to change the state "s" before returning to the integrator.
          For instance, the procedure may adjust "s" so as to satisy
          any constraints or conservation laws which are implied by
          the original differential equations but could be violated in
          a long integration due to the accumulation of error.
        *)
        
      computeStep(dt: Time; error, tol: Coord): Time;
        (*
          Given that a time step of size "dt" generated the specified 
          integration "error", returns the time step that should 
          make the integration error approximately "tol".
        *)
        
      name(): TEXT;
        (*
          Returns the integrator's name, possibly with
          any internal parameters. *)
    END;

  DiffProc = PROCEDURE (t: Time; READONLY s: State; VAR v: Velocity);
    (*
      Called by the integrator to evaluate the right-hand side of the
      differential equation. *)

  CheckStateProc = PROCEDURE (t: Time; READONLY s: State; READONLY v: Velocity): Time;
    (*
      Called by the integrator for every new state along the 
      path.  Returns estimated stopping time. *)

  CheckStepProc = PROCEDURE (
      ta: Time;
      READONLY sa: State;
      READONLY va: Velocity; 
      tb: Time;
      VAR sb: State;
      READONLY er: Error;
    ): Time;
    (*
      Called by the integrator to decide whether to take the 
      step from state "sa" at time "ta" to state "sb" at time "tb". *)

END Integrator.