MODULE SPPWDraw EXPORTS Main;

(* Plots level curves of a piecewise spherical function *)

IMPORT Wr, Rd, Fmt, Text, Thread, Process, ParseParams;
IMPORT SPPWFunction, SPTriang, SPPlot, HLR3, LR3, LR4, PSPlot; 
FROM Stdio IMPORT stderr, stdin;
FROM SPFile IMPORT OpenRead;
FROM PSPlot IMPORT Color, Invisible, Blue, Black;

TYPE Point = SPPWFunction.Point;

TYPE 
  Options = RECORD
      triName: TEXT;        (* Triang file (minus ".tri" extension) or "-" *)
      sfnName: TEXT;        (* Function file (minus ".sfn" extension) or "-" *)
      outName: TEXT;        (* Output file (minus ".ps" extension) *)
      eps: BOOLEAN;         (* TRUE to generate ".eps" instead of ".ps" *)
      obs: HLR3.Point;      (* Position of observer *)
      upp: HLR3.Point;      (* Approximate "up" direction vector *)
      gridSize: CARDINAL;   (* Number of grid steps in Lat & Lon *)
      axes: BOOLEAN;        (* TRUE to draw the coordinate axes *)
      (* Parameters for IsoLines: *)
      fStep: REAL;          (* Value difference between successive level curves *)
      fStart: REAL;         (* Value of first level curve *)
      lineWidth: REAL;      (* Line width in millimeters *)
      (* Parameters for PaintValues: *)
      fMin, fMax: LONGREAL; (* Nominal min and max function values. *)
      light: LR3.T;         (* Light source direction *)
      shadow: REAL;         (* Relative shadow darkening factor *)
    END;

<* FATAL Wr.Failure, Thread.Alerted *>

(* CONST LightBlue = Color{0.8, 1.0, 1.0}; *)

PROCEDURE Main() =
  VAR obs: HLR3.Point;
      map: SPPlot.PMap;
  BEGIN
    WITH 
      o = GetOptions(),
      triRd = OpenTriReader(o.triName),
      tri = MakeTriRef(SPTriang.Read(triRd)),
      sfnRd = OpenSFNReader(o.sfnName),
      fn = NEW(SPPWFunction.T),
      psFile = SPPlot.Open(o.outName, o.eps),
      fMed = (o.fMin + o.fMax) / 2.0d0
    DO  
      fn.tri := tri;
      fn.read(sfnRd);
      IF LR4.Norm(o.obs.c) = 0.0d0 THEN
        obs := HLR3.Point{c := LR4.T{0.0d0, fn.supp[1], fn.supp[2], fn.supp[3]}};
        Wr.PutText(stderr, "using obs = " & FHLR3(obs) & "\n");
      ELSE
        obs := o.obs;
      END;
      map := SPPlot.PerspMap(obs, o.upp);
      
      PROCEDURE EvalFunction(READONLY p: Point): LONGREAL =
        BEGIN
          RETURN fn.eval(p)
        END EvalFunction;
        
      CONST
        FMaxColor = Color{1.0, 0.5, 0.4};
        FMedColor = Color{0.9, 0.9, 0.9};
        FMinColor = Color{0.1, 0.9, 0.9};

      BEGIN
        IF o.axes THEN
          psFile.setLineColor(Black);
          psFile.setFillColor(Invisible);
          psFile.setLineWidth(3.0 * o.lineWidth);
          SPPlot.CoordAxes(psFile, map, -1);
        END;

        psFile.setLineColor(Black);
        psFile.setFillColor(Invisible);
        psFile.setLineWidth(1.5 * o.lineWidth);
        SPPlot.Sphere(psFile);

        Wr.PutText(stderr, "Low shading...\n");
        psFile.setLineColor(Invisible);
        SPPlot.PaintValues(
          psFile, map,
          EvalFunction, 
          M := o.gridSize, N := o.gridSize,
          fMin := o.fMin,
          cMin := FMinColor,
          clipMin := FALSE,
          fMax := fMed,
          cMax := FMedColor,
          clipMax := TRUE,
          dLight := o.light,
          shadow := o.shadow
        ); 
        Wr.PutText(stderr, "High shading...\n");
        SPPlot.PaintValues(
          psFile, map,
          EvalFunction, 
          M := o.gridSize, N := o.gridSize,
          fMin := fMed,
          cMin := FMedColor,
          clipMin := TRUE,
          fMax := o.fMax,
          cMax := FMaxColor,
          clipMax := FALSE,
          dLight := o.light,
          shadow := o.shadow
        ); 
        
        Wr.PutText(stderr, "Triangulation...\n");
        psFile.setLineWidth(o.lineWidth);
        psFile.setLineColor(Black);
        psFile.setFillColor(Invisible);
        SPPlot.Triangulation(psFile, map, tri^);
        
        Wr.PutText(stderr, "Negative isolines...\n");
        psFile.setLineColor(Blue);
        psFile.setLineWidth(0.75 * o.lineWidth);
        SPPlot.Isolines(
          psFile, map,
          EvalFunction, 
          M := o.gridSize, N := o.gridSize,
          fStep := o.fStep,
          fStart := o.fStart,
          sign := -1
        ); 
        Wr.PutText(stderr, "Positive isolines...\n");
        psFile.setLineColor(Color{0.5, 0.0, 0.0});
        SPPlot.Isolines(
          psFile, map,
          EvalFunction, 
          M := o.gridSize, N := o.gridSize,
          fStep := o.fStep,
          fStart := o.fStart,
          sign := 1
        ); 

        IF o.axes THEN
          psFile.setLineColor(Black);
          psFile.setFillColor(Invisible);
          psFile.setLineWidth(3.0 * o.lineWidth);
          SPPlot.CoordAxes(psFile, map, +1);
        END;

        SPPlot.Close(psFile);
      END
    END
  END Main;
  
PROCEDURE MakeTriRef(READONLY tri: SPTriang.T): REF SPTriang.T = 
  BEGIN
    WITH t = NEW(REF SPTriang.T) DO t^ := tri; RETURN t END
  END MakeTriRef;
  
PROCEDURE OpenSFNReader(name: TEXT): Rd.T =
  BEGIN
    IF Text.Equal(name, "-") THEN 
      RETURN stdin
    ELSE
      RETURN OpenRead(name & ".sfn")
    END
  END OpenSFNReader;
  
PROCEDURE OpenTriReader(name: TEXT): Rd.T =
  BEGIN
    IF Text.Equal(name, "-") THEN 
      RETURN stdin
    ELSE
      RETURN OpenRead(name & ".tri")
    END
  END OpenTriReader;

PROCEDURE GetOptions (): Options =
  VAR o: Options;
  BEGIN
    TRY
      WITH  pp = NEW(ParseParams.T).init(stderr) DO 
          
        pp.getKeyword("-triName");                           
        o.triName := pp.getNext();  
      
        IF pp.keywordPresent("-sfnName") THEN                              
          o.sfnName := pp.getNext()
        ELSE
          o.sfnName := "-"
        END;  
      
        pp.getKeyword("-outName");                               
        o.outName := pp.getNext();  
        
        IF pp.keywordPresent("-eps") THEN
          o.eps := TRUE
        ELSIF pp.keywordPresent("-ps") THEN
          o.eps := FALSE
        ELSE
          o.eps := FALSE
        END;
      
        o.axes := pp.keywordPresent("-axes");
      
        IF pp.keywordPresent("-obs") THEN
          o.obs := ParsePointOption(pp)
        ELSE
          (* Default is all zeros, meaning rotate sphere for best view. *)
          o.obs := HLR3.Point{c := LR4.T{0.0d0, 0.0d0, 0.0d0, 0.0d0}}
        END;

        IF pp.keywordPresent("-up") THEN
          o.upp := ParsePointOption(pp)
        ELSE
          o.upp := HLR3.Point{c := LR4.T{0.0d0, 0.0d0, 0.0d0, 1.0d0}}
        END;
                  
        IF pp.keywordPresent("-gridSize") THEN
          o.gridSize := pp.getNextInt(5, 1000)
        ELSE
          o.gridSize := 100
        END;

        IF pp.keywordPresent("-fStep") THEN
          o.fStep := pp.getNextReal()
        ELSE
          o.fStep := 0.1
        END;
      
        IF pp.keywordPresent("-fStart") THEN
          o.fStart := pp.getNextReal()
        ELSE
          o.fStart := 0.0
        END;
        
        IF pp.keywordPresent("-lineWidth") THEN
          o.lineWidth := pp.getNextReal()
        ELSE
          o.lineWidth := 0.2
        END;

        IF pp.keywordPresent("-fMin") THEN
          o.fMin := pp.getNextLongReal()
        ELSE
          o.fMin := -1.0d0
        END;
      
        IF pp.keywordPresent("-fMax") THEN
          o.fMax := pp.getNextLongReal(o.fMin, LAST(LONGREAL))
        ELSE
          o.fMax := +1.0d0
        END;
        
        IF pp.keywordPresent("-light") THEN
          o.light := ParseDirOption(pp)
        ELSE
          o.light := LR3.Dir(LR3.T{0.1d0, 0.2d0, 1.0d0})
        END;
                  
        IF pp.keywordPresent("-shadow") THEN
          o.shadow := pp.getNextReal(0.0, 1.0)
        ELSE
          o.shadow := 0.1
        END;
        
        pp.finish();
      END
    EXCEPT                                                            
    | ParseParams.Error =>                                              
        Wr.PutText(stderr, "Usage: SPPWDraw \\\n");
        Wr.PutText(stderr, "  -triName <name> ] \\\n");
        Wr.PutText(stderr, "  [ -sfnName <name> ] \\\n");
        Wr.PutText(stderr, "  -outName <name> [ -eps | -ps ] \\\n");
        Wr.PutText(stderr, "  [ -obs OW OX OY OZ ] \\\n");
        Wr.PutText(stderr, "  [ -up UW UX UY UZ ] \\\n");
        Wr.PutText(stderr, "  [ -gridSize NUM ] \\\n");
        Wr.PutText(stderr, "  [ -fStep NUM ] [ -fStart NUM ] \\\n");
        Wr.PutText(stderr, "  [ -lineWidth MILLIMETERS ] \\\n");
        Wr.PutText(stderr, "  [ -fMin NUM ] [ -fMax NUM ] \\\n");
        Wr.PutText(stderr, "  [ -light DX DY DZ ] [ -shadow NUM ]\n");
        Process.Exit (1);
    END;
    RETURN o
  END GetOptions;
  
PROCEDURE ParsePointOption(pp: ParseParams.T): HLR3.Point RAISES {ParseParams.Error} =
  VAR p: HLR3.Point;
  BEGIN
    FOR i := 0 TO 3 DO
      p.c[i] := pp.getNextLongReal(-1.0d20, +1.0d20)
    END;
    RETURN p
  END ParsePointOption;  
  
PROCEDURE ParseDirOption(pp: ParseParams.T): LR3.T RAISES {ParseParams.Error} =
  VAR d: LR3.T;
  BEGIN
    FOR i := 0 TO 2 DO
      d[i] := pp.getNextLongReal(-1.0d0, +1.0d0)
    END;
    RETURN LR3.Dir(d)
  END ParseDirOption;  
  
PROCEDURE FHLR3(p: HLR3.Point): TEXT =
  BEGIN
    RETURN
      "{ " &
      Fmt.LongReal(p.c[0], Fmt.Style.Fix, 2) & 
      ", " &
      Fmt.LongReal(p.c[1], Fmt.Style.Fix, 2) & 
      ", " &
      Fmt.LongReal(p.c[2], Fmt.Style.Fix, 2) & 
      ", " &
      Fmt.LongReal(p.c[3], Fmt.Style.Fix, 2) & 
      " }"
  END FHLR3;

BEGIN
  Main();
END SPPWDraw.