GENERIC INTERFACE RDS3(State);

(*
  WHERE
    /\ State.Dim: CONST CARDINAL;
    /\ State.T: TYPE;
    /\ State.Throw: PROCEDURE(lo, hi: REAL): State.T;
    /\ State.Update: 
        PROCEDURE(VAR s: State.T; READONLY d: State.T; dt: REAL; VAR change: REAL);
    /\ State.Select: PROCEDURE(READONLY s: State.T; which: [0..State.Dim-1]): REAL;
  *)

TYPE
  LocalState = ARRAY [-1..+1] OF ARRAY [-1..+1] OF ARRAY [-1..+1] OF State.T;
  (* 
    Texel states in the neighborhood of a grid point.  
    Note: indices are "[dx,dy,dz]" --- opposite to the standard VGM order! *)
  
  Operator = PROCEDURE (READONLY ts: LocalState; x, y, z: CARDINAL): State.T;
  (* 
    Assumed to compute the temporal derivatives "ds/dt" of texel
    with coordinates "(x,y,z)", given its current state "s" and the 
    states of neighboring texels. *)

PROCEDURE Randomize(
    VAR s: ARRAY OF ARRAY OF ARRAY OF State.T;
    lo, hi: REAL;
  );
  (*
    Fills the elements of "s" by calling "State.Throw(lo, hi)". *)
    
PROCEDURE Evolve(
    VAR s: ARRAY OF ARRAY OF ARRAY OF State.T;
    dxdy, dxdz, dydz: CARDINAL;
    VAR d: ARRAY OF ARRAY OF ARRAY OF State.T;
    dt, tol: REAL;
    nit: CARDINAL;
    op: Operator
  );
  (*
    Modifies a discrete periodic 3D texture "s" according to the
    rule "s := s + dt * d", where "d = ds/dt = op(s, x, y, z)".
    
    Note: the texel with coordinates "(x,y,z)" is stored in 
    eleemnt "s[z,y,x]" of the array!
    
    The texture will have the array "s" as its fundamental region, 
    and period vectors "(nx,0,0)", "(dxdy, ny, 0)", "(dxdz, dydz, nz)".
    
    The evolution is stopped once the time derivative "ds/dt"
    becomes less than "tol", or after "nit" iterations. *)

TYPE 
  MapProc = PROCEDURE (READONLY st: State.T): REAL;
    (*
      Should map the state "st" to the range [0__1]. *)

PROCEDURE Write(
    name: TEXT;
    READONLY s: ARRAY OF ARRAY OF ARRAY OF State.T;
    dxdy, dxdz, dydz: CARDINAL;
    map: MapProc;
  );
  (*
    Writes some component of the texture as a VGM file
    called "name.vgm".  
    
    Scaleslinearly the range "[0__1]" returned by the "map" 
    procedure to "[0..255]". *)

(*** BUILDING BLOCKS ***)

PROCEDURE GenerationStep(
    nx, ny, nz: CARDINAL;
    dxdy, dxdz, dydz: CARDINAL;
    VAR s: ARRAY OF ARRAY OF ARRAY OF State.T;
    VAR d: ARRAY OF ARRAY OF ARRAY OF State.T; (* Workspace *)
    op: Operator;
    dt: REAL;
    VAR maxChange, avgChange: REAL;
  );
  (*
    Performs one step of the system's evolution. Namely, computes 
    the derivatives "d := ds/dt" by calling "op", then updates 
    "s := s + d * dt", and computes the maximum and average 
    absolute change in the texel elements. *)

END RDS3.