MODULE ReDifEvolver2D EXPORTS Main;

IMPORT Texture2DGray AS Tx;
IMPORT ParseParams, Process, Wr, Thread, Text, Fmt, FileWr, OSError;
IMPORT Random, Math, LR2;
FROM Stdio IMPORT stderr;

CONST
  Debug = FALSE;

TYPE
  LONG = LONGREAL;
  LONGS = ARRAY OF LONG;

  BOOL = BOOLEAN;
  NAT = CARDINAL;

  TimeDeriv = LONG;
  TimeDerivs = ARRAY OF TimeDeriv;

TYPE
  SpaceDeriv = RECORD
      dx, dy: LONG;
      dxx, dxy, dyy: LONG;
      lap: LONG
    END;

  Nhood = Tx.TexelNeighborhood;
  Nhoods = ARRAY OF Nhood;

  System = OBJECT
      NL: NAT;           (* Number of components (layers) *)
      pixMin: REF LONGS; (* The nominal min texel value per layer *)
      pixMax: REF LONGS; (* The nominal max texel value per layer *)
    METHODS
      print(wr: Wr.T);
        (*
          Prints a description of the system to "wr". *)

      diff(x, y: NAT; t: LONG; READONLY n: Nhoods; VAR d: TimeDerivs);
        (*
          Computes the time derivative "d[i]" for each layer "i",
          given the corresponding neighborhoods "n[i]". *)
    END;

  Options = RECORD
      random: REAL;     (* Amount of random noise to add to input *)
      inName: TEXT;     (* Input filename prefix, or "" if zero texture. *)
      outName: TEXT;    (* Output filename prefix. *)
      UX, VY: NAT;      (* Texture size *)
      VX: NAT;          (* Brick-row displacement *)
      NX, NY: NAT;      (* Output replication factors *)
      system: System;   (* The reaction-diffusion equations *)
      nSteps: NAT;      (* Integration steps *)
      dt: LONG;         (* Integration time step *)
    END;

CONST
  Pi = 3.1415926d0;
  TwoPi = 2.0d0 * Pi;
  CompName = ARRAY OF TEXT{"A", "B", "C", "D", "E", "F", "G", "H"};

PROCEDURE Main() =
  <* FATAL Thread.Alerted, Wr.Failure, OSError.E *>
  VAR ctIter: NAT;
      changed: BOOL;
      t: LONG;
  BEGIN
    WITH
      o = GetOptions(),
      NL = o.system.NL,
      A = NEW(REF ARRAY OF Tx.T, o.system.NL)^,
      DA = NEW(REF ARRAY OF Tx.T, o.system.NL)^,
      rnd = NEW(Random.Default).init(fixed := TRUE),
      zeroTexel = Tx.LongRealTexel{0.0d0}
    DO
      (* Read initial state and allocate working textures: *)
      FOR i := 0 TO NL-1 DO
        IF Text.Equal(o.inName, "") THEN
          A[i] := Tx.New(o.UX, o.VX, o.VY, zeroTexel, zeroTexel);
          Tx.Zero(A[i]);
        ELSE
          WITH
            fName = o.inName & "-" & CompName[i] & ".pgm",
            a = Tx.LongRealTexel{0.5d0 * (o.system.pixMin[i] + o.system.pixMax[i])},
            s = Tx.LongRealTexel{0.5d0 * (o.system.pixMax[i] - o.system.pixMin[i])}
          DO
            A[i] := Tx.Read(fName, o.VX, a, s);
            <* ASSERT A[i].UX = o.UX *>
            <* ASSERT A[i].VX = o.VX *>
            <* ASSERT A[i].VY = o.VY *>
          END
        END;
        IF o.random # 0.0 THEN
          Tx.URandom(A[i], rnd,
            min := Tx.LongRealTexel{FLOAT(o.random, LONGREAL) * o.system.pixMin[i]},
            max := Tx.LongRealTexel{FLOAT(o.random, LONGREAL) * o.system.pixMax[i]},
            clear := FALSE
          )
        END;
        DA[i] := Tx.New(o.UX, o.VX, o.VY, zeroTexel, zeroTexel);
      END;
      (* Integrate system: *)
      ctIter := 0;
      changed := TRUE;
      t := 0.0d0;
      WHILE ctIter < o.nSteps AND changed DO
        changed := Iterate(A, DA, o.system, t, o.dt);
        t := t + o.dt;
        INC(ctIter);
      END;
      (* Write final state and derivatives: *)
      WITH
        ns = Fmt.Pad(Fmt.Int(o.nSteps), 5, '0'),
        fName = o.outName & "-" & ns,
        wr = FileWr.Open(fName & ".txt")
      DO
        o.system.print(wr);
        Wr.PutText(wr, "\n");
        Wr.PutText(wr, "  nSteps = " & Fmt.Int(o.nSteps) & "\n");
        Wr.PutText(wr, "  dt = " & FLR(o.dt, 6) & "\n");
        Wr.PutText(wr, "\n");
        Wr.Flush(wr);
        FOR i := 0 TO NL-1 DO
          WITH
            pixMin = o.system.pixMin[i],
            pixMax = o.system.pixMax[i]
          DO
            Wr.PutText(wr, CompName[i] & " range = [");
            Wr.PutText(wr, FLR(pixMin, 4) & " _ " & FLR(pixMax, 4) & "]");
            Wr.Flush(wr);
            Tx.Write(
              fName & "-" & CompName[i] & ".pgm",
              A[i],
              a := Tx.LongRealTexel{0.5d0 * (pixMin + pixMax)},
              s := Tx.LongRealTexel{0.5d0 * (pixMax - pixMin)},
              NX := o.NX, NY := o.NY
            );
          END;
          WITH
            dev = Tx.InfNorm(DA[i]) + 0.00001d0,
            pixMin = - dev,
            pixMax = + dev
          DO
            Wr.PutText(wr, "d" & CompName[i] & " range = [");
            Wr.PutText(wr, FLR(pixMin, 4) & " _ " & FLR(pixMax, 4) & "]");
            Tx.Write(
              fName & "-d" & CompName[i] & ".pgm",
              DA[i],
              a := zeroTexel,
              s := Tx.LongRealTexel{dev},
              NX := o.NX, NY := o.NY
            );
          END
        END
      END;
    END
  END Main;

PROCEDURE Iterate(
    READONLY A: ARRAY OF Tx.T;
    READONLY DA: ARRAY OF Tx.T;
    system: System;
    t: LONG;
    dt: LONG;
  ): BOOLEAN =
  (*
    Applies one integration step from current texture "A",
    computing the temporal derivatives "DA" and the new texture
    "B = A + dt * DA".

    Returns TRUE if the texture changed.
  *)
  BEGIN
    WITH
      NL = system.NL,
      d = NEW(REF TimeDerivs, NL)^
    DO
      <* ASSERT NUMBER(A) = NL *>
      <* ASSERT NUMBER(DA) = NL *>

      PROCEDURE ComputeDerivatives (* : Tx.MultNeighborhoodVisitProc *) (
          x, y: NAT;
          READONLY n: Nhoods
        ) =
        BEGIN
          system.diff(x, y, t, n, d);
          IF Debug AND x = 1 AND y = 1 THEN DebugDiff(x, y, n, d) END;
          FOR i := 0 TO NL-1 DO
            DA[i].st[y, x] := FLOAT(d[i])
          END;
        END ComputeDerivatives;

      VAR changed: BOOL := FALSE;
      BEGIN
        FOR i := 0 TO NL-1 DO Tx.ZeroBorder(DA[i]) END;
        Tx.MultEnumerateNeighborhoods(A, ComputeDerivatives);
        FOR i := 0 TO NL-1 DO
          Tx.FoldBorder(DA[i]);
          WITH c = Tx.Modify(A[i], dt, DA[i]) DO
            changed := changed OR c
          END
        END;
        RETURN changed
      END
    END;
  END Iterate;

PROCEDURE GetOptions (): Options =
  <* FATAL Thread.Alerted, Wr.Failure *>
  VAR o: Options;
  BEGIN
    WITH pp = NEW(ParseParams.T).init(stderr) DO
      TRY

        IF pp.keywordPresent("-random") THEN
          o.random := pp.getNextReal(0.0, 1000.0)
        ELSE
          o.random := 0.0
        END;
        IF pp.keywordPresent("-inName") THEN
          o.inName := pp.getNext();
        ELSE
          o.inName := "";
        END;

        pp.getKeyword("-outName");
        o.outName := pp.getNext();

        pp.getKeyword("-size");
        o.UX := pp.getNextInt(1, 4096);
        o.VY := pp.getNextInt(1, 4096);

        IF pp.keywordPresent("-offset") THEN
          o.VX := pp.getNextInt(0, o.UX-1)
        ELSE
          o.VX := 0
        END;

        IF pp.keywordPresent("-repeat") THEN
          o.NX := pp.getNextInt(1, 100);
          o.NY := pp.getNextInt(1, 100);
        ELSE
          o.NX := 1;
          o.NY := 1
        END;

        IF pp.keywordPresent("-nSteps") THEN
          o.nSteps := pp.getNextInt(1, 100000)
        ELSE
          o.nSteps := 2000
        END;

        IF pp.keywordPresent("-dt") THEN
          o.dt := pp.getNextLongReal(0.00001d0, 10000.0d0)
        ELSE
          o.dt := 0.01d0
        END;

        pp.getKeyword("-system");
        WITH sysName = pp.getNext() DO
          IF Text.Equal(sysName, "Spots") THEN
            o.system := NEW(SpotsSystem,
              NL := 1,
              pixMin := NewLongs(LONGS{-1.5d0}),
              pixMax := NewLongs(LONGS{+1.5d0}),
              K := pp.getNextLongReal(0.0d0, 1000.0d0)
            )
          ELSIF Text.Equal(sysName, "Dirs") THEN
            o.system := NEW(DirsSystem,
              NL := 2,
              pixMin := NewLongs(LONGS{-1.5d0, -1.5d0}),
              pixMax := NewLongs(LONGS{+1.5d0, +1.5d0}),
              K := pp.getNextLongReal(0.0d0, 1000.0d0)
            )
          ELSIF Text.Equal(sysName, "Lotka") THEN
            o.system := NEW(LotkaSystem,
              NL := 2,
              pixMin := NewLongs(LONGS{-1.0d0, -1.0d0}),
              pixMax := NewLongs(LONGS{+3.0d0, +3.0d0}),
              K := pp.getNextLongReal(0.0d0, 1000.0d0)
            )
          ELSIF Text.Equal(sysName, "Ovals") THEN
            WITH
              K = pp.getNextLongReal(0.0d0, 1000.0d0),
              theta = pp.getNextLongReal(-TwoPi, +TwoPi),
              alpha = pp.getNextLongReal(0.0001d0, 1000.0d0),
              u = LR2.T{Math.cos(theta), Math.sin(theta)}
            DO
              o.system := NEW(OvalsSystem,
                NL := 1,
                pixMin := NewLongs(LONGS{-1.5d0}),
                pixMax := NewLongs(LONGS{+1.5d0}),
                K := K,
                theta := theta,
                alpha := alpha,
                m := SquashedLapCoeffs(u, alpha)
              )
            END
          ELSIF Text.Equal(sysName, "Mods") THEN
            o.system := NEW(ModsSystem,
              NL := 2,
              pixMin := NewLongs(LONGS{-1.5d0, -1.5d0}),
              pixMax := NewLongs(LONGS{+1.5d0, +1.5d0}),
              beta  := pp.getNextLongReal(0.0d0, 1000.0d0),
              gamma := pp.getNextLongReal(0.0d0, 1000.0d0)
            )
          ELSIF Text.Equal(sysName, "Fibers") THEN
            o.system := NEW(FibersSystem,
              NL := 3,
              pixMin := NewLongs(LONGS{-1.5d0, -1.5d0, -1.5d0}),
              pixMax := NewLongs(LONGS{+1.5d0, +1.5d0, +1.5d0}),
              beta  := pp.getNextLongReal(0.0d0, 1000.0d0),
              gamma := pp.getNextLongReal(0.0001d0, 1000.0d0)
            )
          ELSIF Text.Equal(sysName, "Folds") THEN
            o.system := NEW(FoldsSystem,
              NL := 1,
              pixMin := NewLongs(LONGS{-1.5d0}),
              pixMax := NewLongs(LONGS{+1.5d0}),
              alpha := pp.getNextLongReal(0.0d0, 1000.0d0),
              beta  := pp.getNextLongReal(0.0d0, 1000.0d0)
            )
          ELSE
            pp.error("invalid system name = " & sysName)
          END
        END;

        pp.finish();
      EXCEPT
      | ParseParams.Error =>
          Wr.PutText(stderr, "Usage: ReDifEvolver2D\\\n");
          Wr.PutText(stderr, "  [ -inName <name> ] [ -random <amt> ]\\\n");
          Wr.PutText(stderr, "  -outName <name> \\\n");
          Wr.PutText(stderr, "  -size <width> <height> [ -offset <offset> ] \\\n");
          Wr.PutText(stderr, "  -system { Spots <K> | Dirs <K> | ... } \\\n");
          Wr.PutText(stderr, "  [ -nSteps <num> ] [ -dt <timestep> ] \n");
          Process.Exit (1);
      END;
    END;
    RETURN o
  END GetOptions;

PROCEDURE FLR(x: LONG; prec: NAT := 6): TEXT =
  BEGIN
    RETURN Fmt.Pad(Fmt.LongReal(x, prec := prec, style := Fmt.Style.Fix), prec + 6)
  END FLR;

PROCEDURE FR(x: REAL; prec: NAT := 6): TEXT =
  BEGIN
    RETURN Fmt.Pad(Fmt.Real(x, prec := prec, style := Fmt.Style.Fix), prec + 6)
  END FR;

PROCEDURE NewLongs(READONLY a: LONGS): REF LONGS =
  (*
    Modula-3 is so silly....
  *)
  BEGIN
    WITH n = NEW(REF LONGS, NUMBER(a)) DO n^ := a; RETURN n END
  END NewLongs;

PROCEDURE DebugDiff(
    <*UNUSED*> x, y: NAT;
    READONLY n: Nhoods;
    VAR d: TimeDerivs
  ) =
  <* FATAL Thread.Alerted, Wr.Failure *>
  BEGIN
    FOR i := 0 TO LAST(n) DO
      WITH
        nA = n[i], dA = d[i]
      DO
        Wr.PutText(stderr, "--- " & Fmt.Int(i) & " -----------------------\n");
        Wr.PutText(stderr, "d = " & FLR(dA, 6) & "\n");
        Wr.PutText(stderr, "    | " & FR(nA[+1,-1], 4) & FR(nA[+1,00], 4) & FR(nA[+1,+1], 4) & " |\n");
        Wr.PutText(stderr, "n = | " & FR(nA[00,-1], 4) & FR(nA[00,00], 4) & FR(nA[00,+1], 4) & " |\n");
        Wr.PutText(stderr, "    | " & FR(nA[-1,-1], 4) & FR(nA[-1,00], 4) & FR(nA[-1,+1], 4) & " |\n");
      END
    END;
    Wr.PutText(stderr, "-----------------------------\n");
    Wr.Flush(stderr)
  END DebugDiff;

PROCEDURE SpatialDiff(READONLY n: Nhood): SpaceDeriv =
  BEGIN
    WITH
      dx = 0.5d0 * FLOAT(n[00,+1] - n[00,-1], LONG),
      dy = 0.5d0 * FLOAT(n[+1,00] - n[-1,00], LONG),
      dxx = FLOAT(n[00,-1] - n[00,00] - n[00,00] + n[00,+1], LONG),
      dxy = 0.25d0 * FLOAT(n[+1,+1] - n[+1,-1] - n[-1,+1] + n[-1,-1], LONG),
      dyy = FLOAT(n[-1,00] - n[00,00] - n[00,00] + n[+1,00], LONG)
    DO
      RETURN SpaceDeriv{
        dx := dx, dy := dy,
        dxx := dxx, dxy := dxy, dyy := dyy,
        lap := dxx + dyy
      }
    END
  END SpatialDiff;

TYPE
  GenLapCoeffs = RECORD xx, xy, yy: LONG END;

PROCEDURE SquashedLapCoeffs(READONLY u: LR2.T; alpha: LONG): GenLapCoeffs =
  (*
    Computes the generalized Laplacian coefficients that
    will result in a texture squeezed by "alpha" in the
    direction "u" (assumed to be unitary).
  *)
  BEGIN
    WITH
      aa = alpha*alpha
    DO
      RETURN GenLapCoeffs{
        xx := aa*u[0]*u[0] + u[1]*u[1],
        xy := (2.0d0*aa - 2.0d0)*u[0]*u[1],
        yy := u[0]*u[0] + aa*u[1]*u[1]
      }
    END
  END SquashedLapCoeffs;

(* INTERESTING EQUATION SYSTEMS *)

TYPE
  SpotsSystem = System OBJECT
      K: LONG;
    OVERRIDES
      print := SpotsPrint;
      diff := SpotsDiff;
    END;

PROCEDURE SpotsPrint(s: SpotsSystem; wr: Wr.T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr, "Spots: bilevel spots\n");
    Wr.PutText(wr, "  d_t A = A*(1-A^2) + K * Lapl(A)\n");
    Wr.PutText(wr, "  K = " & FLR(s.K, 6) & "\n");
    Wr.Flush(wr)
  END SpotsPrint;

PROCEDURE SpotsDiff(
    s: SpotsSystem;
    <*UNUSED*> x, y: NAT;
    <*UNUSED*> t: LONG;
    READONLY n: Nhoods;
    VAR d: TimeDerivs
  ) =
  BEGIN
    WITH
      nA = n[0], dA = d[0], A = FLOAT(nA[0,0], LONG), gA = SpatialDiff(nA)
    DO
      dA := A*(1.0d0 - A*A) + s.K * gA.lap
    END
  END SpotsDiff;

TYPE
  DirsSystem = System OBJECT
      K: LONG;
    OVERRIDES
      print := DirsPrint;
      diff := DirsDiff
    END;

PROCEDURE DirsPrint(s: DirsSystem; wr: Wr.T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr, "Dirs: direction field\n");
    Wr.PutText(wr, "  d_t A = A*(1-A^2-B^2) + K * Lapl(A)\n");
    Wr.PutText(wr, "  d_t B = B*(1-A^2-B^2) + K * Lapl(B)\n");
    Wr.PutText(wr, "  K = " & FLR(s.K, 6) & "\n");
    Wr.Flush(wr)
  END DirsPrint;

PROCEDURE DirsDiff(
    s: DirsSystem;
    <*UNUSED*> x, y: NAT;
    <*UNUSED*> t: LONG;
    READONLY n: Nhoods;
    VAR d: TimeDerivs
  ) =
  BEGIN
    WITH
      nA = n[0], dA = d[0], A = FLOAT(nA[0,0], LONG), gA = SpatialDiff(nA),
      nB = n[1], dB = d[1], B = FLOAT(nB[0,0], LONG), gB = SpatialDiff(nB)
    DO
      dA := A*(1.0d0 - A*A - B*B) + s.K * gA.lap;
      dB := B*(1.0d0 - A*A - B*B) + s.K * gB.lap
    END
  END DirsDiff;

TYPE
  LotkaSystem = System OBJECT
      K: LONG;
    OVERRIDES
      print := LotkaPrint;
      diff := LotkaDiff
    END;

PROCEDURE LotkaPrint(s: LotkaSystem; wr: Wr.T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr, "Lotka: Lotka-Volterra predator/prey system\n");
    Wr.PutText(wr, "  d_t A = A - A*B + K * Lapl(A)\n");
    Wr.PutText(wr, "  d_t B = A*B - B + K * Lapl(B)\n");
    Wr.PutText(wr, "  K = " & FLR(s.K, 6) & "\n");
    Wr.Flush(wr)
  END LotkaPrint;

PROCEDURE LotkaDiff(
    s: LotkaSystem;
    <*UNUSED*> x, y: NAT;
    <*UNUSED*> t: LONG;
    READONLY n: Nhoods;
    VAR d: TimeDerivs
  ) =
  BEGIN
    WITH
      nA = n[0], dA = d[0], A = FLOAT(nA[0,0], LONG), gA = SpatialDiff(nA),
      nB = n[1], dB = d[1], B = FLOAT(nB[0,0], LONG), gB = SpatialDiff(nB)
    DO
      dA := A*(1.0d0 - B) + s.K * gA.lap;
      dB := B*(A - 1.0d0) + s.K * gB.lap
    END
  END LotkaDiff;

TYPE
  OvalsSystem = System OBJECT
      K: LONG;
      theta: LONG;      (* For documentation only *)
      alpha: LONG;      (* For documentation only *)
      m: GenLapCoeffs;
    OVERRIDES
      print := OvalsPrint;
      diff := OvalsDiff
    END;

PROCEDURE OvalsPrint(s: OvalsSystem; wr: Wr.T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr, "Ovals: squashed bilevel spots\n");
    Wr.PutText(wr, "  d_t A = A*(1-A^2) + K * MLapl(theta, alpha)(A)\n");
    Wr.PutText(wr, "  K =     " & FLR(s.K, 6) & "\n");
    Wr.PutText(wr, "  theta = " & FLR(s.theta, 6) & "\n");
    Wr.PutText(wr, "  alpha = " & FLR(s.alpha, 6) & "\n");
    Wr.Flush(wr)
  END OvalsPrint;

PROCEDURE OvalsDiff(
    s: OvalsSystem;
    <*UNUSED*> x, y: NAT;
    <*UNUSED*> t: LONG;
    READONLY n: Nhoods;
    VAR d: TimeDerivs
  ) =
  BEGIN
    WITH
      nA = n[0], dA = d[0], A = FLOAT(nA[0,0], LONG), gA = SpatialDiff(nA),
      mLapA = s.m.xx * gA.dxx + s.m.xy * gA.dxy + s.m.yy * gA.dyy
    DO
      dA := A*(1.0d0 - A*A) + s.K * mLapA
    END
  END OvalsDiff;

TYPE
  ModsSystem = System OBJECT
      beta: LONG;
      gamma: LONG;
    OVERRIDES
      print := ModsPrint;
      diff := ModsDiff
    END;

PROCEDURE ModsPrint(s: ModsSystem; wr: Wr.T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr, "Mods: spots modulated by permeability\n");
    Wr.PutText(wr, "  d_t A = 0\n");
    Wr.PutText(wr, "  d_t B = B*(1-B^2) + (beta + gama * A) * Lapl(B)\n");
    Wr.PutText(wr, "  beta =  " & FLR(s.beta, 6) & "\n");
    Wr.PutText(wr, "  gamma = " & FLR(s.gamma, 6) & "\n");
    Wr.Flush(wr)
  END ModsPrint;

PROCEDURE ModsDiff(
    s: ModsSystem;
    <*UNUSED*> x, y: NAT;
    <*UNUSED*> t: LONG;
    READONLY n: Nhoods;
    VAR d: TimeDerivs
  ) =
  BEGIN
    WITH
      nA = n[0], dA = d[0], A = FLOAT(nA[0,0], LONG),
      nB = n[1], dB = d[1], B = FLOAT(nB[0,0], LONG), gB = SpatialDiff(nB)
    DO
      dA := 0.0d0;
      dB := B*(1.0d0 - B*B) + MAX(0.0d0, (s.beta + s.gamma * A)) * gB.lap
    END
  END ModsDiff;

TYPE
  FibersSystem = System OBJECT
      beta: LONG;
      gamma: LONG;
    OVERRIDES
      print := FibersPrint;
      diff := FibersDiff
    END;

PROCEDURE FibersPrint(s: FibersSystem; wr: Wr.T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr, "Fibers: self-orienting fibers\n");
    Wr.PutText(wr, "  d_t A = 0\n");
    Wr.PutText(wr, "  d_t B = 0\n");
    Wr.PutText(wr, "  d_t C = C*(1-C^2) + beta * MLapl(C, (A,B), gamma)\n");
    Wr.PutText(wr, "  beta =  " & FLR(s.beta, 6) & "\n");
    Wr.PutText(wr, "  gamma = " & FLR(s.gamma, 6) & "\n");
    Wr.Flush(wr)
  END FibersPrint;

PROCEDURE FibersDiff(
    s: FibersSystem;
    <*UNUSED*> x, y: NAT;
    <*UNUSED*> t: LONG;
    READONLY n: Nhoods;
    VAR d: TimeDerivs
  ) =
  BEGIN
    WITH
      nA = n[0], dA = d[0], A = FLOAT(nA[0,0], LONG),
      nB = n[1], dB = d[1], B = FLOAT(nB[0,0], LONG),
      nC = n[2], dC = d[2], C = FLOAT(nC[0,0], LONG), gC = SpatialDiff(nC),
      m = SquashedLapCoeffs(LR2.T{A, B}, s.gamma),
      mLapC = m.xx * gC.dxx + m.xy * gC.dxy + m.yy * gC.dyy
    DO
      dA := 0.0d0;
      dB := 0.0d0;
      dC := C*(1.0d0 - C*C) + MAX(0.0d0, s.beta) * mLapC
    END
  END FibersDiff;

TYPE
  FoldsSystem = System OBJECT
      alpha : LONG;
      beta : LONG;
    OVERRIDES
      print := FoldsPrint;
      diff := FoldsDiff
    END;

PROCEDURE FoldsPrint(s: FoldsSystem; wr: Wr.T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr, "Folds: spontaneous folds\n");
    Wr.PutText(wr, "  d_t A = A*(1-A^2) + Lapl(A)*(alpha + beta * Bend(A))\n");
    Wr.PutText(wr, "  Bend(A) = (d_xx A)*(d_yy A) - (d_xy A)^2\n");
    Wr.PutText(wr, "  alpha = " & FLR(s.alpha, 6) & "\n");
    Wr.PutText(wr, "  beta =  " & FLR(s.beta, 6) & "\n");
    Wr.Flush(wr)
  END FoldsPrint;

PROCEDURE FoldsDiff(
    s: FoldsSystem;
    <*UNUSED*> x, y: NAT;
    <*UNUSED*> t: LONG;
    READONLY n: Nhoods;
    VAR d: TimeDerivs
  ) =
  BEGIN
    WITH
      nA = n[0], dA = d[0], A = FLOAT(nA[0,0], LONG), gA = SpatialDiff(nA),
      bendA = gA.dxx * gA.dyy - gA.dxy * gA.dxy
    DO
      dA := A*(1.0d0 - A*A) + gA.lap * (s.alpha + s.beta * bendA);
    END
  END FoldsDiff;

(*
  Equation System Template:

TYPE
  ~System = System OBJECT
      ~: LONG;
      ~: LONG;
    OVERRIDES
      print := ~Print;
      diff := ~Diff
    END;

PROCEDURE ~Print(s: ~System; wr: Wr.T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr, "~:\n");
    Wr.PutText(wr, "  ~\n");
    Wr.PutText(wr, "  ~\n");
    Wr.PutText(wr, "  ~\n");
    Wr.PutText(wr, "  ~ = " & FLR(s.~, 6) & "\n");
    Wr.PutText(wr, "  ~ = " & FLR(s.~, 6) & "\n");
    Wr.Flush(wr)
  END ~Print;

PROCEDURE ~Diff(
    s: ~System;
    <*UNUSED*> x, y: NAT;
    READONLY n: Nhoods;
    VAR d: TimeDerivs
  ) =
  BEGIN
    WITH
      nA = n[0], dA = d[0], A = FLOAT(nA[0,0], LONG), gA = SpatialDiff(nA),
      nB = n[1], dB = d[1], B = FLOAT(nB[0,0], LONG), gB = SpatialDiff(nB),
      nC = n[2], dC = d[2], C = FLOAT(nC[0,0], LONG), gC = SpatialDiff(nC),
      ~
    DO
      dA := A*(1.0d0 - A*A - B*B) + s.alpha * gA.lap;
      dB := B*(1.0d0 - A*A - B*B) + s.alpha * gB.lap;
      dC := C*(1.0d0 - C*C) + MAX(0.0d0, s.beta) * mLapC
    END
  END ~Diff;

*)

BEGIN
  Main();
END ReDifEvolver2D.