MODULE SpherePlot;

IMPORT Wr, Thread, FileStream, Math, H2, R3, R3x3, PSPlot;

REVEAL 
  T = Public BRANDED OBJECT
      f: Wr.T;
      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: *)
    IF FLOAT(u[0], LONGREAL) > 1.0d0/t.d THEN
      PSPlot.SetPen(t.f, 0.0, 0.3, 0.0, 0.0);
      radius := 0.5*radius;
      (* Should fill polygon with black *)
    ELSE
      PSPlot.SetPen(t.f, 0.0, 0.1, 0.0, 0.0);
      (* 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(R3.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 = R3x3.Rotation(ang, u)
    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.f, vAnt, v, t.d)
      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
    PSPlot.SetPen(t.f, 0.0, width, 0.0, 0.0);
    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(R3.Cross(u, v)),
      r = R3x3.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.f, uPrev, u, t.d)
        END;
        INC(i)
      END;
      DrawShortSeg(t.f, u, v, t.d);
    END;
  END T_Seg;
  
PROCEDURE DrawShortSeg(f: Wr.T; READONLY a, b: R3.T; d: LONGREAL) =
  BEGIN
    WITH
      wa = FLOAT(a[0], LONGREAL),
      sa = d/(d - wa),

      wb = FLOAT(b[0], LONGREAL),
      sb = d/(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
      PSPlot.DrawSegment(f, xa, ya, xb, yb)
    END
  END DrawShortSeg;
  
PROCEDURE T_Close(t: T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    PSPlot.DrawFrame(t.f);
    PSPlot.EndPage(t.f);
    PSPlot.EndFile(t.f);
    Wr.Close(t.f);
  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. *)
  BEGIN
    WITH 
      t = NEW(T)
    DO
      t.f := FileStream.OpenWrite(name);
      PSPlot.BeginFile(t.f);
      PSPlot.BeginPage(t.f, 1, -1.1D0, 1.1D0, -1.1D0, 1.1D0, 1, 1);
      PSPlot.AddCaption(t.f, name);
      t.m := R3x3.Mul(
          R3x3.Rotation(Turn, R3.T{0.0, 0.0, 1.0}),
          R3x3.Rotation(Tilt, R3.T{0.0, 1.0, 0.0})
        );
      t.d := 5.0d0;
      PSPlot.DrawCircle(t.f, 0.0d0, 0.0d0, t.d/Math.sqrt(t.d*t.d-1.0d0));
      RETURN t
    END;
  END New;

BEGIN
  
END SpherePlot.