MODULE SpherePlot;

IMPORT Wr, Thread, Math, H2, R3, R3Extras, R3x3, R3x3Extras, PSPlot;  
FROM PSPlot IMPORT LetterXSize, LetterYSize, Color;

REVEAL 
  T = Public BRANDED OBJECT
      f: PSPlot.PSFile;
      m: R3x3.T;     (* View rotation matrix *)
      d: LONGREAL;   (* Viewing distance (1.0 = on the sphere). *)
    OVERRIDES
      frontSide := T_frontSide;
      point := T_Point;
      seg := T_Seg;
      close := T_Close;
    END;
    
PROCEDURE T_frontSide(t: T; a: H2.Point): BOOLEAN =
  BEGIN
    WITH u = R3.Dir(R3x3.MapRow(a.c, t.m)) DO
      RETURN FLOAT(u[0], LONGREAL) > 1.0d0/t.d
    END
  END T_frontSide;

PROCEDURE T_Point(t: T; a: H2.Point; radius: REAL) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  CONST 
    N = 8;
  VAR u, v, vAnt: R3.T;
  BEGIN
    u := R3.Dir(R3x3.MapRow(a.c, t.m));
    (* Check if point is visible: *)
    t.f.setLineColor(Color{0.0, 0.0, 0.0});
    t.f.setLineSolid();
    IF FLOAT(u[0], LONGREAL) > 1.0d0/t.d THEN
      t.f.setLineWidth(0.3);
      radius := 0.5*radius;
      (* Should fill polygon with black *)
    ELSE
      t.f.setLineWidth(0.1);
      (* Should fill polygon with white *)
    END;
    
    (* Compute a vector "v" orthogonal to "u": *)
    WITH
      u0 = u[0]*u[0],
      u1 = u[1]*u[1],
      u2 = u[2]*u[2]
    DO
      IF u0 >= u1 AND u0 >= u2 THEN
        v := R3.T{-u[1], u[0], 0.0}
      ELSIF u1 >= u2 THEN
        v := R3.T{0.0, -u[2], u[1]}
      ELSE
        v := R3.T{u[2], 0.0, -u[0]}
      END;
    END;
    v := R3.Dir(R3Extras.Cross(u, v));
    
    (* Now draw circle: *)
    WITH
      r = FLOAT(radius, LONGREAL),
      s = Math.sqrt(1.00d0 - r*r),
      ang = (2.0d0 * 3.1415926535897932384626433833d0)/FLOAT(N, LONGREAL),
      rm = R3x3Extras.Rotation(ang, u, 1.0d0)
    DO
      v := R3.Dir(R3.Add(R3.Scale(r, v), R3.Scale(s, u)));
      FOR i := 0 TO N-1 DO
        vAnt := v;
        v := R3x3.MapRow(v, rm);
        DrawShortSeg(t, vAnt, v)
      END;
    END
  END T_Point;
  
PROCEDURE T_Seg(t: T; a, b: H2.Point; width: REAL) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  CONST 
    Eps = 0.01d0;
    DashSteps = 2;
    StopRadius = 2.0d0 * Eps;
  VAR u, v, uPrev: R3.T; 
      i: CARDINAL;
  BEGIN
    t.f.setLineWidth(width);
    t.f.setLineColor(Color{0.0, 0.0, 0.0});
    t.f.setLineSolid();
    u := R3.Dir(R3x3.MapRow(a.c, t.m));
    v := R3.Dir(R3x3.MapRow(b.c, t.m));
    IF R3.Dist(u, v) < 0.001d0 THEN RETURN END;
    WITH
      w = R3.Dir(R3Extras.Cross(u, v)),
      r = R3x3Extras.Rotation(Eps, w)
    DO
      i := 0;
      WHILE R3.Dist(u, v) > StopRadius DO
        uPrev := u;
        u := R3.Dir(R3x3.MapRow(u, r));
        IF FLOAT(u[0], LONGREAL) > 1.0d0/t.d
        OR (i DIV DashSteps) MOD 2 = 0 THEN
          DrawShortSeg(t, uPrev, u)
        END;
        INC(i)
      END;
      DrawShortSeg(t, u, v);
    END;
  END T_Seg;
  
PROCEDURE DrawShortSeg(t: T; READONLY a, b: R3.T) =
  BEGIN
    WITH
      wa = FLOAT(a[0], LONGREAL),
      sa = t.d/(t.d - wa),

      wb = FLOAT(b[0], LONGREAL),
      sb = t.d/(t.d - wb),

      xa = sa*FLOAT(a[1], LONGREAL),
      ya = sa*FLOAT(a[2], LONGREAL),
      
      xb = sb*FLOAT(b[1], LONGREAL),
      yb = sb*FLOAT(b[2], LONGREAL)
    DO
      t.f.segment(xa, ya, xb, yb)
    END
  END DrawShortSeg;
  
PROCEDURE T_Close(t: T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    WITH
      s = t.f
    DO
      s.endDrawing();
      s.endPage();
      s.close()
    END  
  END T_Close;
  
PROCEDURE New(name: TEXT): T =
  <* FATAL Wr.Failure, Thread.Alerted *>
  CONST 
    Turn = -0.10d0;  (* West-to-east turn. *)
    Tilt = 0.35d0;   (* Top-to-viewer tilt. *)
  VAR t : T;
  BEGIN
    t := NEW(T);
    WITH
      s  = NEW(PSPlot.PSFile).open(name)
    DO
      s.beginPage();
      s.beginDrawing(
          xSize := 160.0d0, ySize := 160.0d0,
          xCenter := LetterXSize/2.0d0, yCenter := LetterYSize/2.0d0
        );
      s.setLineColor(Color{0.0, 0.0, 0.0});
      s.setLineWidth(0.5);
      s.setLineSolid();
      s.setScale(PSPlot.Axis.X, -1.0d0,1.0d0);
      s.setScale(PSPlot.Axis.Y, -1.0d0,1.0d0);
      s.caption(name);
      t.m := R3x3.Mul(
               R3x3Extras.Rotation(Turn, R3.T{0.0, 0.0, 1.0}),
               R3x3Extras.Rotation(Tilt, R3.T{0.0, 1.0, 0.0})
             );
      t.d := 1.0001d0;
      s.setFillColor(Color{1.0, 1.0, 1.0});
      s.circle(0.0d0, 0.0d0, FLOAT(t.d/Math.sqrt(t.d*t.d-1.0d0),REAL));
      s.setFillColor(Color{0.0, 0.0, 0.0});
      t.f := s
    END;
    RETURN t
  END New;

BEGIN
  
END SpherePlot.