MODULE SPPlot;

FROM Math IMPORT cos, sin, sqrt;
IMPORT Wr, Thread, SPRange, SPTriang, PSPlot, HLR3, LR3, LR3Extras, LR4;
FROM SPTriang IMPORT Arc, Org, Dest;
FROM SPQuad IMPORT NullArc;
FROM PSPlot IMPORT Color, Invisible, Blue, Black;
FROM Stdio IMPORT stderr;

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

PROCEDURE Open(name: TEXT; eps: BOOLEAN): File =
  VAR f: PSPlot.File;
  BEGIN
    IF eps THEN
      f :=  NEW(PSPlot.EPSFile).open(name & ".eps");
    ELSE
      WITH fp = NEW(PSPlot.PSFile).open(name & ".ps") DO
        fp.beginPage();
        fp.beginDrawing();
        f := fp
      END
    END;
    (* Set plotting scale. *)
    (* Assume that the unit sphere will project onto the unit circle. *)
    f.setScale(PSPlot.Axis.X, -1.3D0, 1.3D0);
    f.setScale(PSPlot.Axis.Y, -1.3D0, 1.3D0);
    RETURN f
  END Open;
  
PROCEDURE Close(f: PSPlot.File) =
  BEGIN
    TYPECASE f OF
    | PSPlot.EPSFile => 
        (* Nothing to do *)
    | PSPlot.PSFile(fp) =>
        fp.endDrawing();
        fp.endPage();
    ELSE
      <* ASSERT FALSE *>
    END;
    f.close()
  END Close;
  
PROCEDURE PerspMap(READONLY obs, upp: HLR3.Point): HLR3.PMap =
  BEGIN
    (* 
      Choose the image plane so that the sphere's 
      outline will have unit radius. *)
    WITH 
      w = obs.c[0], x = obs.c[1], y = obs.c[2], z = obs.c[3],
      rr = x*x + y*y + z*z,
      den = rr + sqrt(rr * (rr - w*w))
    DO
      <* ASSERT w >= 0.0d0 *>
      <* ASSERT rr > w*w *>
      WITH
        foc = HLR3.Point{c := LR4.T{den, w*x, w*y, w*z}}
      DO
        RETURN HLR3.PerspMap(obs := obs, foc := foc, upp := upp)
      END
    END
  END PerspMap;

PROCEDURE Sphere(f: File) =
  BEGIN
    (*
      Assumes the perspective is such that the sphere's outline
      becomes the unit circle. *)
    f.ellipse(0.0d0, 0.0d0, 1.0d0, 0.0d0, 0.0d0, 1.0d0)
  END Sphere;
  
PROCEDURE Axis(
    f: File; 
    READONLY map: PMap; 
    dir: LR3.T;
    length: LONG;
    sense: SIGN := 0;
  ) =
  CONST 
    N = 20;
    TipLength = 0.150d0;
    TipRadius = 0.015d0;
    TwoPi = 2.0d0 * 3.1415926535897932384626433833d0;
  BEGIN
    (* Head of arrow *)
    WITH
      c = LR3.Dir(dir),     (* Point where arrow exits sphere. *)
      cm = Project(c, map)
    DO
      IF sense = 0
      OR sense = -1 AND NOT PointIsVisible(cm)
      OR sense = +1 AND PointIsVisible(cm)
      THEN
        WITH
          t = LR3.Scale(length, dir), (* Tip of arrow. *)
          tm = Project(t, map),
          b = LR3.Scale(length - TipLength, dir), (* End of shaft. *)
          (* Two orthogonal directions: *)
          u = FindOrthoDir(dir),
          v = LR3Extras.Cross(dir, u)
        DO
          f.segment(cm[0], cm[1], tm[0], tm[1]);
          FOR i := 1 TO N DO
            WITH
              alpha = TwoPi*FLOAT(i, LONG)/FLOAT(N, LONG),
              c = cos(alpha),
              s = sin(alpha),
              p = LR3.Add(b, LR3.Mix(c*TipRadius, u, s*TipRadius, v)),
              pm = Project(p, map)
            DO
              f.segment(tm[0], tm[1], pm[0], pm[1])
            END
          END
        END
      END
    END;
    (* Tail of arrow *)
    WITH
      c = LR3.Dir(LR3.Neg(dir)),     (* Point where arrow enters sphere. *)
      cm = Project(c, map)
    DO
      IF sense = 0
      OR sense = -1 AND NOT PointIsVisible(cm)
      OR sense = +1 AND PointIsVisible(cm)
      THEN
        WITH
          t = LR3.Scale(-length, dir), (* Tail of arrow. *)
          tm = Project(t, map)
        DO
          f.segment(cm[0], cm[1], tm[0], tm[1]);
        END
      END
    END
  END Axis;
  
PROCEDURE CoordAxes(
    f: File; 
    READONLY map: PMap; 
    sense: SIGN;
  ) =
  BEGIN
    Axis(f, map, dir := LR3.T{1.0d0, 0.0d0, 0.0d0}, length := 1.25d0, sense := sense);
    Axis(f, map, dir := LR3.T{0.0d0, 1.0d0, 0.0d0}, length := 1.25d0, sense := sense);
    Axis(f, map, dir := LR3.T{0.0d0, 0.0d0, 1.0d0}, length := 1.25d0, sense := sense);
  END CoordAxes; 

PROCEDURE FindOrthoDir(d: LR3.T): LR3.T =
  VAR iMax, iMed, iMin: CARDINAL;
      u: LR3.T;
  BEGIN
    (* Find the max, med, and min coordinates of "d": *)
    IF ABS(d[0]) >= ABS(d[1]) THEN
      iMax := 0; iMed := 1
    ELSE
      iMax := 1; iMed := 0
    END;
    IF ABS(d[2]) <= ABS(d[iMed]) THEN
      iMin := 2
    ELSIF ABS(d[2]) <= ABS(d[iMax]) THEN
      iMin := iMed; iMed := 2;
    ELSE
      iMin := iMed; iMed := iMax; iMax := 2
    END;
    (* Build result: *)
    u[iMax] := -d[iMed]; u[iMed] := d[iMax]; u[iMin] := 0.0d0;
    RETURN LR3.Dir(u)
  END FindOrthoDir;

PROCEDURE Segment(
    f: File; 
    READONLY map: PMap; 
    p, q: Coords; 
    N: CARDINAL := 20;
  ) =
  VAR am: LR3.T := Project(p, map);
  BEGIN
    FOR i := 1 TO N DO
      WITH
        s = FLOAT(i, LONG)/FLOAT(N, LONG),
        b = LR3.Dir(LR3.Mix(1.0d0 - s, p, s, q)),
        bm = Project(b, map)
      DO
        IF PointIsVisible(am) AND PointIsVisible(bm) THEN
          f.segment(am[0], am[1], bm[0], bm[1])
        END;
        am := bm
      END
    END
  END Segment;

PROCEDURE Triangulation(
    f: File;
    READONLY map: PMap;
    READONLY tri: SPTriang.T;
    N: CARDINAL := 20;
  ) =

  PROCEDURE PlotDelaunayEdge(e:Arc) =
    BEGIN 
      WITH 
        op = Org(e).c, 
        dp = Dest(e).c 
      DO
        Segment(f, map, op, dp, N);
      END;
    END PlotDelaunayEdge;
    
  BEGIN
    WITH arc = tri.arc^ DO
      FOR i := 0 TO LAST(arc) BY 2 DO
        IF arc[i] # NullArc THEN PlotDelaunayEdge(arc[i]) END
      END
    END
  END Triangulation;

PROCEDURE Isolines(
    f: File;
    READONLY map: PMap;
    func: ScalarField;
    M, N: CARDINAL;
    fStep: REAL;
    fStart: REAL := 0.0;
    sign: SIGN;
  )=
      
  PROCEDURE ProcessTriangle(
      READONLY P: Coords; fP: LONG;
      READONLY Q: Coords; fQ: LONG;
      READONLY R: Coords; fR: LONG
    ) =
     VAR kMax, kMin: INTEGER;
     
    BEGIN (* ProcessTriangle *)
      WITH
        fStep = FLOAT(fStep, LONG),
        fStart = FLOAT(fStart, LONG),
        fMin = MIN(MIN(fP, fQ), fR),
        fMax = MAX(MAX(fP, fQ), fR),
        Pm = Project(P, map),
        Qm = Project(Q, map),
        Rm = Project(R, map)
      DO
        IF NOT TriangIsVisible(Pm, Qm, Rm) THEN RETURN END;
        
        IF sign = -1 THEN
          kMin := CEILING((fMin-fStart)/fStep);
          kMax := CEILING((MIN(0.0d0, fMax)-fStart)/fStep) - 1;
        ELSIF  sign = 1 THEN
          kMin := FLOOR((MAX(0.0d0, fMin)-fStart)/fStep) + 1;
          kMax := FLOOR((fMax-fStart)/fStep);
        END;
        FOR k := kMin TO kMax DO
          WITH
            fLevel = fStart + FLOAT(k, LONG) * fStep
          DO
            PlotZerosInTriangle(
              f,
              Pm[0], Pm[1], fP - fLevel,
              Qm[0], Qm[1], fQ - fLevel,
              Rm[0], Rm[1], fR - fLevel
            )
          END
        END;
      END;
    END ProcessTriangle; 
  
  BEGIN
    WITH 
      prevLon = NEW(REF ARRAY OF LONG, M),
      thisLon = NEW(REF ARRAY OF LONG, M),
      prevFun = NEW(REF ARRAY OF LONG, M),
      thisFun = NEW(REF ARRAY OF LONG, M)
    DO 
    
      PROCEDURE Initialize() =
        (* Plots isolines around the south pole. *)
        VAR  P, Q : Coords;
        BEGIN
          WITH
            pi = 3.1415926535897932384626433833d0,
            poleC = Coords{0.0d0, 0.0d0, -1.0d0},
            poleF = func(poleC),
            lat = FLOAT(1,LONG) * pi/FLOAT(N,LONG) - pi/2.0d0
          DO
            prevLon[0] := 0.0d0;
            P := XYZFromLatLon(lat, prevLon[0]);
            prevFun[0] := func(P);
            FOR j:= 1 TO M DO 
              prevLon[j MOD M]:= prevLon[j-1] + 2.0d0*pi/FLOAT(M,LONG);
              Q := XYZFromLatLon(lat, prevLon[j MOD M]);
              prevFun[j MOD M] := func(Q);
              ProcessTriangle(poleC, poleF, P, prevFun[j-1], Q, prevFun[j MOD M]);
              P := Q;
            END;
          END;
        END Initialize;

      PROCEDURE MainLoop() =
        (* 
          Plots isolines zone by zone, from south to north,
          except near poles. Each zone is a row of north-pointing
          triangles, interleaved witha row of down-pointing triangles.
        *)
        VAR O, P, Q, R : Coords;
            fO, fR : LONG;
            longo: LONG;
        BEGIN
          FOR i := 2 TO N-1 DO
            WITH
              pi = 3.1415926535897932384626433833d0,
              prevLat = FLOAT(i-1,LONG)*pi/FLOAT(N,LONG) - pi/2.0d0,
              latat = FLOAT(i,LONG)*pi/FLOAT(N,LONG) - pi/2.0d0
            DO 
              longo := 2.0d0*pi*0.5d0*FLOAT(i,LONG)/FLOAT(M,LONG);
              O := XYZFromLatLon(latat, longo);
              fO := func(O);
              P := XYZFromLatLon(prevLat, prevLon[0]);
              FOR  j:= 1 TO M  DO
                WITH 
                  longr = longo + 2.0d0*pi/FLOAT(M,LONG)
                DO 
                  Q := XYZFromLatLon(prevLat, prevLon[j MOD M]);
                  R := XYZFromLatLon(latat, longr);
                  fR := func(R);
                  thisLon[j-1] := longo;
                  thisFun[j-1] := fO;
                  ProcessTriangle(O, fO, P, prevFun[j-1], Q, prevFun[j MOD M]);
                  ProcessTriangle(O, fO, Q, prevFun[j MOD M], R, fR);
                  fO := fR;
                  P := Q;
                  O := R;
                  longo := longr;
                END; 
              END; 
              FOR j := 0 TO M-1 DO 
                prevLon[j] := thisLon[j];
                prevFun[j] := thisFun[j]
              END;
            END;
          END;
        END MainLoop;

      PROCEDURE Finalize() =
        (* Plots isolines around north poleC. *)
        VAR P, Q : Coords;
        BEGIN
          WITH
            pi = 3.1415926535897932384626433833d0,
            poleC = Coords{0.0d0, 0.0d0, 1.0d0},
            poleF = func(poleC),
            lat = FLOAT(N-1,LONG)*pi/FLOAT(N,LONG) - pi/2.0d0
          DO
            P := XYZFromLatLon(lat, thisLon[0]);
            FOR j:= 1 TO M DO
              Q := XYZFromLatLon(lat, thisLon[j MOD M]);
              ProcessTriangle(poleC, poleF, P,  thisFun[j-1],
                            Q,  thisFun[j MOD M]);
              P := Q;              
            END;
          END;
        END Finalize;
     
      BEGIN
        Initialize();
        MainLoop();
        Finalize();
      END;
    END;
  END Isolines;
  
PROCEDURE PaintValues(
    f: File;
    READONLY map: PMap;     (* Perspective projection matrix. *)
    func: ScalarField;      (* Function to plot. *)
    M, N: CARDINAL;         (* Grid size. *)
    fMin: LONG;             (* Minimum function value. *)
    cMin: Color;            (* Color to use for "fMin". *)
    clipMin: BOOLEAN;       (* TRUE omits parts below "fMin". *)
    fMax: LONG;             (* Maximum function VALUE. *)
    cMax: Color;            (* Color to use for "fMax". *)
    clipMax: BOOLEAN;       (* TRUE omits parts above "fMax". *)
    dLight: LR3.T;          (* Direction towards main light source. *)
    shadow: REAL := 0.1;    (* Amount of darkening by shadow. *)
  ) =
  PROCEDURE ProcessTriangle(
      READONLY P: Coords; fP: LONG;
      READONLY Q: Coords; fQ: LONG;
      READONLY R: Coords; fR: LONG
    ) =

    BEGIN (* ProcessTriangle *)
      WITH
        Pm = Project(P, map),
        Qm = Project(Q, map),
        Rm = Project(R, map)
      DO
        IF NOT TriangIsVisible(Pm, Qm, Rm) THEN RETURN END;
        WITH
          fMed = (fP + fQ + fR) / 3.0d0,
          tc = InterpolateColor(fMed, fMin, cMin, clipMin, fMax, cMax, clipMax)
        DO
          IF tc = PSPlot.Invisible THEN RETURN END;
          WITH
            dMed = LR3.Dir(LR3.Add(LR3.Add(P, Q), R)),
            illum = 1.0 + shadow * FLOAT(LR3.Cos(dMed, dLight), REAL),
            ac = ClipColor(Color{illum*tc[0], illum*tc[1], illum*tc[2]})
          DO
            f.setFillColor(ac);
            f.triangle(Pm[0], Pm[1], Qm[0], Qm[1], Rm[0], Rm[1])
          END
        END;
      END;
    END ProcessTriangle; 
  
  BEGIN
    WITH 
      prevLon = NEW(REF ARRAY OF LONG, M),
      thisLon = NEW(REF ARRAY OF LONG, M),
      prevFun = NEW(REF ARRAY OF LONG, M),
      thisFun = NEW(REF ARRAY OF LONG, M)
    DO 
    
      PROCEDURE Initialize() =
        (* Plots isolines around the south pole. *)
        VAR  P, Q : Coords;
        BEGIN
          WITH
            pi = 3.1415926535897932384626433833d0,
            poleC = Coords{0.0d0, 0.0d0, -1.0d0},
            poleF = func(poleC),
            lat = FLOAT(1,LONG) * pi/FLOAT(N,LONG) - pi/2.0d0
          DO
            prevLon[0] := 0.0d0;
            P := XYZFromLatLon(lat, prevLon[0]);
            prevFun[0] := func(P);
            FOR j:= 1 TO M DO 
              prevLon[j MOD M]:= prevLon[j-1] + 2.0d0*pi/FLOAT(M,LONG);
              Q := XYZFromLatLon(lat, prevLon[j MOD M]);
              prevFun[j MOD M] := func(Q);
              ProcessTriangle(poleC, poleF, P, prevFun[j-1], Q, prevFun[j MOD M]);
              P := Q;
            END;
          END;
        END Initialize;

      PROCEDURE MainLoop() =
        (* 
          Plots isolines zone by zone, from south to north,
          except near poles. Each zone is a row of north-pointing
          triangles, interleaved witha row of down-pointing triangles.
        *)
        VAR O, P, Q, R : Coords;
            fO, fR : LONG;
            longo: LONG;
        BEGIN
          FOR i := 2 TO N-1 DO
            WITH
              pi = 3.1415926535897932384626433833d0,
              prevLat = FLOAT(i-1,LONG)*pi/FLOAT(N,LONG) - pi/2.0d0,
              latat = FLOAT(i,LONG)*pi/FLOAT(N,LONG) - pi/2.0d0
            DO 
              longo := 2.0d0*pi*0.5d0*FLOAT(i,LONG)/FLOAT(M,LONG);
              O := XYZFromLatLon(latat, longo);
              fO := func(O);
              P := XYZFromLatLon(prevLat, prevLon[0]);
              FOR  j:= 1 TO M  DO
                WITH 
                  longr = longo + 2.0d0*pi/FLOAT(M,LONG)
                DO 
                  Q := XYZFromLatLon(prevLat, prevLon[j MOD M]);
                  R := XYZFromLatLon(latat, longr);
                  fR := func(R);
                  thisLon[j-1] := longo;
                  thisFun[j-1] := fO;
                  ProcessTriangle(O, fO, P, prevFun[j-1], Q, prevFun[j MOD M]);
                  ProcessTriangle(O, fO, Q, prevFun[j MOD M], R, fR);
                  fO := fR;
                  P := Q;
                  O := R;
                  longo := longr;
                END; 
              END; 
              FOR j := 0 TO M-1 DO 
                prevLon[j] := thisLon[j];
                prevFun[j] := thisFun[j]
              END;
            END;
          END;
        END MainLoop;

      PROCEDURE Finalize() =
        (* Plots isolines around north poleC. *)
        VAR P, Q : Coords;
        BEGIN
          WITH
            pi = 3.1415926535897932384626433833d0,
            poleC = Coords{0.0d0, 0.0d0, 1.0d0},
            poleF = func(poleC),
            lat = FLOAT(N-1,LONG)*pi/FLOAT(N,LONG) - pi/2.0d0
          DO
            P := XYZFromLatLon(lat, thisLon[0]);
            FOR j:= 1 TO M DO
              Q := XYZFromLatLon(lat, thisLon[j MOD M]);
              ProcessTriangle(poleC, poleF, P,  thisFun[j-1],
                            Q,  thisFun[j MOD M]);
              P := Q;              
            END;
          END;
        END Finalize;
     
      BEGIN
        Initialize();
        MainLoop();
        Finalize();
      END;
    END;
  END PaintValues;
  
PROCEDURE InterpolateColor(
    f: LONG;
    fMin: LONG;             (* Minimum function value. *)
    cMin: Color;            (* Color to use for "fMin". *)
    clipMin: BOOLEAN;       (* TRUE omits parts below "fMin". *)
    fMax: LONG;             (* Maximum function VALUE. *)
    cMax: Color;            (* Color to use for "fMax". *)
    clipMax: BOOLEAN;       (* TRUE omits parts above "fMax". *)
  ): Color =
  BEGIN
    IF f < fMin THEN
      IF clipMin THEN RETURN PSPlot.Invisible ELSE RETURN cMin END
    ELSIF f > fMax THEN
      IF clipMax THEN RETURN PSPlot.Invisible ELSE RETURN cMax END
    ELSE
      WITH
        s = FLOAT((f - fMin)/(fMax - fMin), REAL),
        t = 1.0 - s
      DO
        RETURN Color{
          s*cMax[0] + t*cMin[0],
          s*cMax[1] + t*cMin[1],
          s*cMax[2] + t*cMin[2]
        }
      END
    END
  END InterpolateColor;
  
PROCEDURE ClipColor(c: Color): Color =
  BEGIN
    WITH
      m = MAX(MAX(c[0], c[1]), c[2])
    DO
      IF m <= 1.0 THEN 
        RETURN c 
      ELSE
        WITH
          a = Color{c[0]/m, c[1]/m, c[2]/m},
          yc = 0.3*c[0] + 0.6*c[1] + 0.1*c[2],
          ya = yc/m,
          s = MIN(1.0, yc) - ya
        DO
          RETURN Color{
            a[0] + s*(1.0 - a[0]),
            a[1] + s*(1.0 - a[1]),
            a[2] + s*(1.0 - a[2])
          }
        END
      END
    END
  END ClipColor;

PROCEDURE XYZFromLatLon(lat, lon: LONG): Coords =
  BEGIN
    WITH 
      cLat = cos(lat), sLat = sin(lat), 
      cLon = cos(lon), sLon = sin(lon)
    DO
      RETURN Coords{cLat*cLon, cLat*sLon, sLat}
    END
  END XYZFromLatLon;

PROCEDURE PlotZerosInTriangle( 
    f: File;
    Px, Py: LONG;  Pf: LONG;
    Qx, Qy: LONG;  Qf: LONG;
    Rx, Ry: LONG;  Rf: LONG;
  ) =
  (*
    Given a triangle "P,Q,R" on the plane,
    and function values at its vertices, interpolates a 
    linear function through that data, and plots the set of 
    points where that function is zero. 
  *)
  BEGIN
    IF (Pf = 0.0d0) AND (Qf = 0.0d0) AND (Rf = 0.0d0) THEN
      (* Ignore triangle. *)
      RETURN
      (*
        f.segment(Px, Py, Qx, Qy);
        f.segment(Qx, Qy, Rx, Ry);
        f.segment(Rx, Ry, Px, Py);
        f.setFillColor(PSPlot.Black);
        f.triangle(Px, Py, Qx, Qy, Rx, Ry);
      *)
    END;
    IF (Pf = 0.0d0) AND (Qf = 0.0d0) THEN 
      f.segment(Px, Py, Qx, Qy);
    END;
    IF (Qf = 0.0d0) AND (Rf = 0.0d0) THEN 
      f.segment(Qx, Qy, Rx, Ry);
    END;
    IF (Pf = 0.0d0) AND (Rf = 0.0d0) THEN
      f.segment(Px, Py, Rx, Ry);
    END;
    IF (ORD(Pf = 0.0d0) + ORD(Qf = 0.0d0) + ORD(Rf = 0.0d0)) >= 2 THEN
      RETURN; 
    END;
    IF (Pf >= 0.0d0) AND (Qf >= 0.0d0) AND (Rf >= 0.0d0) THEN 
      RETURN;
    END;
    IF (Pf <= 0.0d0) AND (Qf <= 0.0d0) AND (Rf <= 0.0d0) THEN 
      RETURN;
    END;
    (* Now there is at least one positive corner and one negative corner *)
    (* Permute corners so that Rf <= Pf <= Qf: *)
    IF Rf > Pf THEN
      VAR temp: LONG;
      BEGIN
        temp := Px; Px := Rx; Rx := temp;
        temp := Py; Py := Ry; Ry := temp;
        temp := Pf; Pf := Rf; Rf := temp;
      END;
    END;
    IF Pf > Qf THEN
      VAR temp: LONG;
      BEGIN
        temp:= Px; Px := Qx; Qx := temp;
        temp:= Py; Py := Qy; Qy := temp;
        temp:= Pf; Pf := Qf; Qf := temp; 
      END;
    END;
    IF Rf > Pf THEN
      VAR temp: LONG;
      BEGIN
        temp := Px; Px := Rx; Rx := temp;
        temp := Py; Py := Ry; Ry := temp;
        temp := Pf; Pf := Rf; Rf := temp;
      END;
    END;
    <* ASSERT Rf <= Pf *>
    <* ASSERT Pf <= Qf *>
    <* ASSERT Rf < 0.0d0 *>
    <* ASSERT Qf > 0.0d0 *>
    (* Now swap q,r so that sign(Pf) # sign(Rf): *)
    IF Pf < 0.0d0 THEN
      VAR temp: LONG;
      BEGIN
        temp := Qx; Qx := Rx; Rx := temp;
        temp := Qy; Qy := Ry; Ry := temp;
        temp := Qf; Qf := Rf; Rf := temp; 
      END;
    END;
    (*now Qf # 0.0d0, Rf # 0.0d0 *)
    (* also Pf = 0.0d0, OR (sign(Pf) # sig(Qf) AND sig(Pf) # sig(Rf)) *)
    <* ASSERT Qf # 0.0d0 *>
    <* ASSERT Rf # 0.0d0 *>
    <* ASSERT (Qf > 0.0d0) # (Rf > 0.0d0) *>
    <* ASSERT 
      (Pf = 0.0d0) OR ((Pf > 0.0d0) = (Qf > 0.0d0) AND (Pf > 0.0d0) # (Rf > 0.0d0))
    *>
    IF (Pf # 0.0d0) THEN 
      WITH
        SPR = Rf/(Rf - Pf),
        SRP = Pf/(Pf - Rf),
        PRx = Px*SPR + Rx*SRP,
        PRy = Py*SPR + Ry*SRP,
        SQR = Rf/(Rf - Qf),
        SRQ = Qf/(Qf - Rf),
        QRx = Qx*SQR + Rx*SRQ,
        QRy = Qy*SQR + Ry*SRQ
      DO
        f.segment(PRx, PRy, QRx, QRy);
      END
    ELSE
      WITH
        SQR = Rf/(Rf - Qf),
        SRQ = Qf/(Qf - Rf),
        QRx = Qx*SQR + Rx*SRQ,
        QRy = Qy*SQR + Ry*SRQ
      DO
        <* ASSERT (Qf > 0.0d0) # (Rf > 0.0d0) *>
        f.segment(Px, Py, QRx, QRy)
      END;
    END;
  END PlotZerosInTriangle;

PROCEDURE Vectors(
    f: File;
    READONLY map: PMap;
    fld: VectorField;
    N: CARDINAL;
    scale: LONG;
  ) =
      
  PROCEDURE ProcessPoint(READONLY p: Coords) =

    BEGIN (* ProcessPoint *)
      WITH pm = Project(p, map) DO
        IF NOT PointIsVisible(pm) THEN RETURN END;
        WITH 
          q = LR3.Add(p, LR3.Scale(scale, fld(p))), 
          qm = Project(q, map)
        DO
          f.circle(pm[0], pm[1], 0.2);
          f.segment(pm[0], pm[1], qm[0], qm[1])
        END
      END
    END ProcessPoint; 
  
  VAR r: LR3.T;
  BEGIN
    FOR k := 0 TO 2 DO
      FOR s := -1 TO +1 BY 2 DO
        FOR i := 0 TO N-1 DO
          FOR j := 0 TO N-1 DO
            WITH
              u = 2.0d0 * (FLOAT(i, LONG) + 0.5d0)/FLOAT(N, LONG) - 1.0d0,
              v = 2.0d0 * (FLOAT(j, LONG) + 0.5d0)/FLOAT(N, LONG) - 1.0d0
            DO
              r[k] := FLOAT(s, LONG);
              r[(k+1) MOD 3] := u;
              r[(k+2) MOD 3] := v;
              r := LR3.Dir(r);
              ProcessPoint(r)
            END
          END
        END
      END
    END;
  END Vectors;

PROCEDURE Project(READONLY p: Coords; READONLY map: HLR3.PMap): LR3.T =
  BEGIN
    WITH
      q = HLR3.Point{c := LR4.T{1.0d0, p[0], p[1], p[2]}},
      g = HLR3.MapPoint(q, map),
      w = g.c[0]
    DO
      RETURN LR3.T{g.c[1]/w, g.c[2]/w, g.c[3]/w}
    END;
  END Project;

PROCEDURE TriangIsVisible(READONLY Pm, Qm, Rm: LR3.T): BOOLEAN =
  (* Given perspective-mapped corners, decide whether triangle is visible *)
  BEGIN
    WITH
      plane = HLR3.PlaneFromThreePoints(
        HLR3.FromCartesian(Pm),
        HLR3.FromCartesian(Qm),
        HLR3.FromCartesian(Rm)
      ),
      sCenter = HLR3.Side(HLR3.Point{c := LR4.T{1.0d0, 0.0d0, 0.0d0, 0.0d0}}, plane),
      sObs = HLR3.Side(HLR3.Point{c := LR4.T{0.0d0, 0.0d0, 0.0d0, 1.0d0}}, plane)
    DO
      RETURN sCenter # sObs
    END
  END TriangIsVisible;                   

PROCEDURE PointIsVisible(READONLY pm: Coords): BOOLEAN =
  (* 
    Given a perspective-mapped point "pm", tests if it is visible.
    More precisely, tests whether it it lies in front of or behind
    the plane that contains the sphere's horizon. *)
  BEGIN
    RETURN pm[2] >= 0.0d0
  END PointIsVisible;                   

PROCEDURE Everything(
    func: ScalarField;      (* The function to plot *)
    tri: REF SPTriang.T;    (* Triangulation to draw, or NIL *)
    fMax: LONG;             (* Nominal maximum absolute function value *)
    fStep: REAL;            (* Isoline spacing. *)
    M, N: CARDINAL;         (* Grid size. *)
    obs: HLR3.Point;        (* Observer's position *)
    upp: HLR3.Point;        (* Camera vertical reference *)
    dLight: LR3.T;          (* Direction towards main light source. *)
    outName: TEXT;          (* Output file name, minus extension *)
    eps: BOOL;              (* TRUE for Encapsulated Postscript (".eps") format *)
    verbose: BOOL := TRUE;  (* TRUE to print messages along the way. *)
  ) =
  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};
    LineWidth = 0.2;
  BEGIN
    WITH
      map = PerspMap(obs, upp),
      psFile = Open(outName, eps)
    DO
      IF verbose THEN Wr.PutText(stderr, "Back axes...\n") END;
      psFile.setLineColor(Black);
      psFile.setFillColor(Invisible);
      psFile.setLineWidth(3.0 * LineWidth);
      CoordAxes(psFile, map, -1);

      psFile.setLineColor(Black);
      psFile.setFillColor(Invisible);
      psFile.setLineWidth(1.5 * LineWidth);
      Sphere(psFile);

      IF verbose THEN Wr.PutText(stderr, "Low shading...\n") END;
      psFile.setLineColor(Invisible);
      PaintValues(
        psFile, map,
        func, 
        M, N,
        fMin := -fMax, 
        cMin := FMinColor,
        clipMin := FALSE,
        fMax := 0.0d0,
        cMax := FMedColor,
        clipMax := TRUE,
        dLight := dLight,
        shadow := 0.1
      ); 
      IF verbose THEN Wr.PutText(stderr, "High shading...\n") END;
      PaintValues(
        psFile, map,
        func, 
        M := M, N := N,
        fMin := 0.0d0,
        cMin := FMedColor,
        clipMin := TRUE,
        fMax := fMax,
        cMax := FMaxColor,
        clipMax := FALSE,
        dLight := dLight,
        shadow := 0.1
      ); 

      IF tri # NIL THEN
        IF verbose THEN Wr.PutText(stderr, "Triangulation...\n") END;
        psFile.setLineWidth(LineWidth);
        psFile.setLineColor(Black);
        psFile.setFillColor(Invisible);
        Triangulation(psFile, map, tri^);
      END;

      IF verbose THEN Wr.PutText(stderr, "Negative isolines...\n") END;
      psFile.setLineColor(Blue);
      psFile.setLineWidth(0.75 * LineWidth);
      Isolines(
        psFile, map,
        func, 
        M, N,
        fStep := fStep,
        fStart := -0.5 * fStep,
        sign := -1
      ); 
      IF verbose THEN Wr.PutText(stderr, "Positive isolines...\n") END;
      psFile.setLineColor(Color{0.5, 0.0, 0.0});
      Isolines(
        psFile, map,
        func, 
        M, N,
        fStep := fStep,
        fStart := +0.5 * fStep,
        sign := +1
      ); 

      IF verbose THEN Wr.PutText(stderr, "Front axes...\n") END;
      psFile.setLineColor(Black);
      psFile.setFillColor(Invisible);
      psFile.setLineWidth(3.0 * LineWidth);
      CoordAxes(psFile, map, +1);

      Close(psFile);
    END
  END Everything;

PROCEDURE BothSides(
    func: ScalarField;      (* The function to plot. *)
    tri: REF SPTriang.T;    (* Triangulation to draw, or NIL *)
    fMax: LONG;             (* Nominal MAX(ABS(fun)), or 0 to compute it. *)
    fStep: REAL;            (* Isoline spacing. *)
    M, N: CARDINAL;         (* Grid size. *)
    obs: HLR3.Point;        (* Observer's position *)
    upp: HLR3.Point;        (* Camera vertical reference *)
    dLight: LR3.T;          (* Direction towards main light source. *)
    outName: TEXT;          (* Output file name, minus extension *)
    eps: BOOL;              (* TRUE for Encapsulated Postscript (".eps") format *)
    verbose: BOOL := TRUE;  (* TRUE to print messages along the way. *)
  ) =
  VAR tag: TEXT;
  BEGIN
    IF fMax = 0.0d0 THEN
      fMax := RoundMax(SPRange.OnSphere(func, 40))
    END;
    
    FOR side := -1 TO +1 BY 2 DO
      IF side = +1 THEN tag :=  "v" ELSE tag := "f" END;
      Everything(
        func := func,
        tri := tri,
        outName := outName & tag,
        fMax := fMax,
        fStep := fStep,
        M := M, N := N,
        obs := obs,
        upp := upp,
        dLight := dLight,
        eps := eps,
        verbose := verbose
      );
      (* Reverse position of observer: *)
      FOR k := 1 TO 3 DO obs.c[k] := - obs.c[k] END;
    END
  END BothSides;

PROCEDURE RoundMax(x: LONG): LONG =
  VAR z: LONG := 1.0d0;
  BEGIN
    (* Allow for rounding errors in the algorithm below: *)
    x := 0.9999999d0 * x;
    (* If too small, don't bother: *)
    IF x < 1.0d-30 THEN RETURN 1.0d-30 END;
    (* Finds the smallest power of 10 that is no less than "x": *)
    WHILE z/10.0d0 > x DO z := z/10.0d0 END;
    WHILE z <= x DO z := z * 10.0d0 * z END;
    (* Finds a nice multiple of that power of 10: *)
    IF x < 0.25d0 * z THEN
      RETURN 0.25d0 * z 
    ELSIF x < 0.5d0 * z THEN 
      RETURN 0.5d0 * z
    ELSE
      RETURN z
    END
  END RoundMax;

BEGIN
END SPPlot.