TYPE
  Attractor = RECORD 
      nCharges: CARDINAL;               (* Number of discrete charges *)
      pt: REF ARRAY OF StateInvariants; (* Position of charges *)
      wt: REF ARRAY OF REAL;            (* Magnitude of charges *)
      totCharge: LONGREAL;              (* Sum of "wt[r]" for all "r" *)
      radius: LONGREAL;                 (* Nominal radius of each charge *)
    END;
    (*
      An "Attractor" is a weighted set of invariant vectors.
      Think of each element as a fuzzy charge at point "q[i]"
      of five-dimensional space, with approximate radius "r"
      and magnitude "w". *)
      
PROCEDURE WriteAttractor(name: TEXT; READONLY a: Attractor);
  (*
    Writes the attarctos to a file "name.plot".  Each line
    contains "2 + 4*State.N" real numbers, separated by spaces:
    the radius, the weight, and the attractor coordinates. *)

PROCEDURE LocalEnergyFromNeighborhood(
    READONLY n: StateNeighborhood;        (* The full neighborhood of a texel *)
    READONLY a: Attractor;                (* The reference charges *)
    scale: LONGREAL;                      (* Energy scaling factor *)
    VAR (*OUT*) E: LONGREAL;              (* Energy associated with this texel *)
    VAR (*OUT*) dEdn: StateNeighborhood;  (* Partials "dE/dn[y,x,k]" *)
  );
  (*
    Computes the energy term "E" at some texel of some texture "t",
    from the local neighborhood "n[-1..+1][-1..+1][k]" of the texture.
    Also computes the partial derivatives of "E" relative to the
    components of "n". *)

PROCEDURE LocalEnergyFromDerivatives(
    READONLY v: StateDerivatives;       (* Spatial derivatives of texture *)
    READONLY a: Attractor;              (* Reference charges *)
    scale: LONGREAL;                    (* Energy scaling factor *)
    VAR (*OUT*) E: LONGREAL;            (* Energy of texture *)
    VAR (*OUT*) dEdv: StateDerivatives; (* Partials dE/dv[k][i]  *)
  );
  (*
    Computes the energy term "E" at some texel of some texture "t",
    from the local spatial derivatives "v[k][0..5]" of the texture.
    Also computes the partial derivatives of "E" relative to the
    elements of "v". *)
    
PROCEDURE LocalEnergyFromInvariants(
    READONLY p: StateInvariants;        (* Local invariants of a texel *)
    READONLY a: Attractor;              (* The reference charges *)
    scale: LONGREAL;                    (* Energy scaling factor *)
    VAR (*OUT*) E: LONGREAL;            (* Energy associated with this texel *)
    VAR (*OUT*) dEdp: StateInvariants;  (* Partials "dE/dp[k][i]" *)
  );
  (*
    Computes the energy term "E" at some texel of some texture "t",
    from the local isotropic invariants "p[k][0..4]" of the texture.
    Also computes the partial derivatives of "E" relative to the
    elements of "p". *)

PROCEDURE Add(VAR s; READONLY t: T);
  (*
    Performs "s := s+t". *)

PROCEDURE Modify(VAR s: T; alpha: LONGREAL; READONLY t: T): BOOLEAN;
  (*
    Performs "s := s + alpha * t"; returns TRUE if s changed. *)

PROCEDURE Mix(alpha: REAL; VAR s: T; beta: LONGREAL; READONLY t: T);
  (*
    Performs "s := alpha * s + beta * t". *)

PROCEDURE InfNorm(READONLY s: T): LONGREAL;
  (*
    Returns the max absolute value of the components of "s". *)

PROCEDURE Mix(alpha: REAL; READONLY s: T; beta: REAL; READONLY t: T): T =
  BEGIN RETURN alpha*s + beta*t END Mix;
  
PROCEDURE Add(READONLY s, t: T): T =
  BEGIN RETURN s+t END Add;
  
PROCEDURE InfNorm(READONLY s: T): REAL =
  BEGIN RETURN ABS(s) END InfNorm;
  
              Pixel.GetCoords(tyx, tc);
              Pixel.getCoords(fyx, fc);
              FOR k := 0 TO Pixel.N-1 DO 
                WITH tck = tc[k], fck = fc[k] DO
                  pak[0] := sa * pak[0] + sb * pbk[0];
                  pak[1] := sa * pak[1] + sb * pbk[1];
                  pak[2] := sa * pak[2] + sb * pbk[2];
                  pak[3] := sa * pak[3] + sb * pbk[3];
                  pak[4] := sa * pak[4] + sb * pbk[4];
                END
              END;
              nt = Pixel.Mix(1.0, ot, dt, fy[x])
            DO
              IF nt # ot THEN changed := TRUE; ot := nt END


(*
  print neighborhod:
      PutText(stderr, "    n   =\n");
      FOR y := +1 TO -1 BY -1 DO 
        PutText(stderr, "      ");
        FOR x := -1 TO +1 BY 1 DO
          PutText(stderr, Fmt.Real(n[y,x]) & " ");
        END;
        PutText(stderr, "\n");
      END;
    PutText(stderr, "\n");

    PutText(stderr, "    ck = " & Fmt.Real(ck));
    PutText(stderr, "    d00 = " & Fmt.Real(d00) & "\n");

PROCEDURE Eval(READONLY t: T; x, y: REAL): Pixel.RealT =
  BEGIN
    Reduce(t, x, y);
    WITH 
      ix = FLOOR(x),
      iy = FLOOR(y),
      sx = x - FLOAT(ix),
      sy = y - FLOAT(iy)
    DO
      IF sx = 0.0 AND sy = 0.0 THEN
        (* Fetch and return texel: *)
    RETURN t.v[j,i]
  END Eval;

*)