(* Copyright 1993 DCC-IMECC, Universidade Estadual de Campinas *)
MODULE VTGDirSpots EXPORTS Main;

IMPORT Wr, ParseParams, Random, Process, Thread;
IMPORT R1State, VTGR1, VTGR1Extra, R3State, VTGR3;
FROM Stdio IMPORT stderr;

<*UNUSED*>
CONST
  Comment =
    "Defining equation:\n" &
    "\n" &
    "  da/dt = a(1 - a^2) + K U_alpha^v(a) \n" &
    "\n" &
    "The operator U is a Laplacian stretched by alpha in\n" &
    "the direction of the unit vector field v, which is provided\n" &
    "by three input VGM files.";

TYPE
  Options = RECORD
      name: TEXT;                   (* File name prefix for output *)
      field: ARRAY [0..2] OF TEXT;  (* File names for direction field *)
      nx, ny, nz: CARDINAL;         (* Texture size *)
      dxdy, dxdz, dydz: CARDINAL;   (* Texture offsets *)
      dt: REAL;                     (* Time step *)
      tol: REAL;                    (* Cutoff value of time derivative *)
      fill: REAL;                   (* Fraction of positive/negative seed vals *)
      latDiff: REAL;                (* Transversal diffusion constant. *)
      lngDiff: REAL;                (* Longitudinal diffusion constant. *)
      nit: CARDINAL;                (* Max interations *)
    END;

PROCEDURE Main() =
  BEGIN
    WITH
      o = GetOptions(),
      dir = NEW(REF VTGR3.StateArray, o.nz, o.ny, o.nx)^,
      isoCoef = o.latDiff,
      dirCoef = o.lngDiff - o.latDiff
    DO

      VAR ds: VTGR1Extra.Derivatives;
      
      PROCEDURE DirSpotsOp(
          READONLY ts: VTGR1.LocalState;
          x, y, z: CARDINAL
        ): R1State.T =
        VAR d: R1State.T;
        BEGIN
          VTGR1Extra.ComputeDerivatives(ts, ds);
	  WITH
            f = ds.f,
	    u = dir[z,y,x],
            ux = u[0],
	    uy = u[1],
	    uz = u[2],
            f2c = 1.0 - f*f,
	    isoLap = ds.fxx + ds.fyy + ds.fzz,
	    dirLap =
	      (ds.fxx * ux + ds.fxy * uy + ds.fxz * uz) * ux +
	      (ds.fxy * ux + ds.fyy * uy + ds.fyz * uz) * uy +
	      (ds.fxz * ux + ds.fyz * uy + ds.fzz * uz) * uz
	  DO
	    d := f * f2c + isoCoef * isoLap + dirCoef * dirLap;
	    (* Ensure convergence: *)
	    WITH fnew = f + 2.0 * d * o.dt DO
	      IF ABS(fnew) > 1.2 THEN
		WITH ffix = MIN(1.2, MAX(-1.2, fnew)) DO
		  d := (ffix - f)/(2.0 * o.dt)
		END
	      END
	    END
	  END;
          RETURN d
        END DirSpotsOp;

      BEGIN
        WITH
          s = NEW(REF VTGR1.StateArray, o.nz, o.ny, o.nx)^,
          d = NEW(REF VTGR1.StateArray, o.nz, o.ny, o.nx)^,
	  rand = NEW(Random.Default).init(fixed := TRUE)
        DO
          FOR i := 0 TO 2 DO
	    ReadDirComponent(o.field[i], dir, i, o.dxdy, o.dxdz, o.dydz)
	  END;
	  WITH
	    lo = MAX(-1.0, 2.0 * o.fill - 2.0),
	    hi = MIN(+1.0, 2.0 * o.fill)
	  DO
	    VTGR1.Randomize(s, lo, hi, rand);
	  END;
          VTGR1.Evolve(
            s,
            o.dxdy, o.dxdz, o.dydz,
            d,
            o.dt, o.tol, o.nit,
            DirSpotsOp
          );
          WriteVal(o.name, s, o.dxdy, o.dxdz, o.dydz);
        END
      END
    END
  END Main;

PROCEDURE ReadDirComponent(
    name: TEXT;
    VAR dir: VTGR3.StateArray;
    which: [0..2];
    dxdy, dxdz, dydz: INTEGER;
  ) =

  PROCEDURE PutComp(VAR s: R3State.T; v: REAL) =
    BEGIN
      s[which] := 2.0 * v - 1.0
    END PutComp;

  BEGIN
    VTGR3.Read(name, dir, dxdy, dxdz, dydz, put := PutComp)
  END ReadDirComponent;

PROCEDURE WriteVal(
    name: TEXT;
    READONLY s: VTGR1.StateArray;
    dxdy, dxdz, dydz: INTEGER;
  ) =

  PROCEDURE GetComp(READONLY st: R1State.T): REAL =
    BEGIN
     RETURN 0.5 + 0.5 * st
    END GetComp;

  BEGIN
    VTGR1.Write(name, s, dxdy, dxdz, dydz, get := GetComp);
  END WriteVal;

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

	pp.getKeyword("-field");
	FOR i := 0 TO 2 DO
	  o.field[i] := pp.getNext();
	END;

	IF pp.keywordPresent("-size") THEN
	  o.nx := pp.getNextInt(1, 100);
	  o.ny := pp.getNextInt(1, 100);
	  o.nz := pp.getNextInt(1, 100);
	ELSE
	  o.nx := 40;
	  o.ny := 40;
	  o.nz := 40;
	END;

	IF pp.keywordPresent("-offsets") THEN
	  o.dxdy := pp.getNextInt(0, o.nx-1);
	  o.dxdz := pp.getNextInt(0, o.nx-1);
	  o.dydz := pp.getNextInt(0, o.ny-1);
	ELSE
	  o.dxdy := 0;
	  o.dxdz := 0;
	  o.dydz := 0;
	END;

	IF pp.keywordPresent("-fill") THEN
	  o.fill := pp.getNextReal(0.0, 1.0)
	ELSE
	  o.fill := 0.5
	END;

	IF pp.keywordPresent("-lngDiff") THEN
	  o.lngDiff := pp.getNextReal(0.0, 100.0)
	ELSE
	  o.lngDiff := 1.0
	END;

	IF pp.keywordPresent("-latDiff") THEN
	  o.latDiff := pp.getNextReal(0.0, 100.0)
	ELSE
	  o.latDiff := 0.5
	END;

	IF pp.keywordPresent("-dt") THEN
	  o.dt := pp.getNextReal(0.0001, 1000.0)
	ELSE
	  o.dt := 0.05
	END;

	IF pp.keywordPresent("-tol") THEN
	  o.tol := pp.getNextReal(0.0001, 1000.0)
	ELSE
	  o.tol := 0.001
	END;

	IF pp.keywordPresent("-nit") THEN
	  o.nit := pp.getNextInt(0, 10000)
	ELSE
	  o.nit := 200
	END;

	pp.finish();
      EXCEPT
      | ParseParams.Error =>
	  Wr.PutText(stderr, "Usage: VTGDirSpots\\\n");
	  Wr.PutText(stderr, "  -name <name>\\\n");
	  Wr.PutText(stderr, "  -field <name> <name> <name>\\\n");
	  Wr.PutText(stderr, "  [ -size <nx> <ny> <nz> ]\\\n");
	  Wr.PutText(stderr, "  [ -offsets <dxdy>, <dxdz>, <dydz> ]\\\n");
	  Wr.PutText(stderr, "  [ -fill <num> ]\\\n");
	  Wr.PutText(stderr, "  [ -latDiff <num> ] [-lngDiff <num>]\\\n");
	  Wr.PutText(stderr, "  [ -tol <num> ] [ -dt <num> ] [ -nit <num> ]\n");
	  Process.Exit (1);
      END;
    END;
    RETURN o
  END GetOptions;

BEGIN
  Main()
END VTGDirSpots.