planes = ComputeFacePlanes(t, coords)^,
      
PROCEDURE ComputeFacePlanes(
    READONLY t: Topology;
    READONLY coords: Coords;
  ): REF ARRAY OF H3.Plane =
  BEGIN
    WITH
      rpl = NEW(REF ARRAY OF H3.Plane, t.NF),
      pl = rpl^
    DO
      FOR i := 0 TO t.NF-1 DO
        pl[i] := FacePlane(Rot(t.side[i]), coords)
      END;
      RETURN rpl
    END;
  END ComputeFacePlanes;
  
PROCEDURE FacePlane(a: Arc; READONLY coords: Coords): H3.Plane =
  BEGIN
    WITH
      b = Onext(a),
      vo = H3.FromCartesian(coords[NOrg(a)]),
      va = H3.FromCartesian(coords[NOrg(Sym(a))]),
      vb = H3.FromCartesian(coords[NOrg(Sym(b))]),
      plane =  H3.Join(vo, va, vb)
    DO
      RETURN plane
    END
  END FacePlane;



PROCEDURE ZeroVal(sc: T; org: REFANY);
  (*
    Sets to zero the "val" fields of all triangles of 
    vertices of "sc" that 
    are incident to specific original elements.
    
    Specifically, sets to zero the "val" fields of any vertex 
    incident to
    
      every vertex "v" with "v.org" in "{vorg,eorg,torg}" , and
      
      every edge "e" with "e.org" in "{eorg,torg}", and
      
      every triangle "t" with "t.org = torg".
      
  *)
  BEGIN
    WHILE sc.v # NIL DO
      WITH v = sc.v.v^ DO
        IF v.org = vorg OR v.org = eorg OR v.org = torg THEN
          v.val := 0.0d0
        END
      END;
      sc.v := sc.v.rest
    END;
    WHILE sc.e # NIL DO
      WITH e = sc.e.e^ DO
        IF e.org = eorg OR e.org = torg THEN
          FOR i := 0 TO 1 DO e.vert[i].val := 0.0d0 END
        END
      END;
      sc.e := sc.e.rest
    END;
    WHILE sc.t # NIL DO
      WITH t = sc.t.t^ DO
        IF t.org = eorg OR t.org = torg THEN
          FOR j := 0 TO 2 DO 
            FOR i := 0 TO 1 DO t.edge[j].vert[i].val := 0.0d0 END
          END
        END
      END;
      sc.t := sc.t.rest
    END;
  END ZeroVal;

  
  FrgVertexList = REF RECORD
      v: REF Vertex;
      rest: VertexList
    END;

  FrgEdgeList = REF RECORD
      v: REF FrgEdge;
      rest: FrgEdgeList
    END;
 
  FrgTriangleList = REF RECORD
      v: REF FrgTriangle;
      rest: FrgTriangleList
    END;

  Vertex = RECORD        (* Input vertex *)
      vno: VertexNo;         (* Original vertex id *)
      icz: H3.Point;         (* ICS coordinates of vertex *)
    END;

  VertexNo = CARDINAL;

  Edge = RECORD          (* Input edge *)
      eno: EdgeNo;           (* Original edge id *)
      vert: VertexPair;      (* Edge endpoints *)
    END;

  VertexPair = ARRAY [0..2] OF REF Vertex;

  Triangle = RECORD      (* Input triangle *)
      no: TriangleNo;       (* Original triangle id *)
      edge: EdgeTriplet;        (* Edges *)
      edir: EDirTriplet;        (* Edge directions *)
      xmin, xmax: REAL;         (* Range of ICS X coordinate *)
      ymin, ymax: REAL;         (* Range of ICS Y coordinate *)
      zmin, zmax: REAL;         (* Range of ICS Z coordinate *)
      norm: R3.T;               (* Normal direction in WCS *)
      colr: Colors;             (* Color fractions by light source *)
      visible: BOOLEAN;         (* (Work) TRUE if partly visible *)
    END;

  EdgeTriplet = ARRAY [0..2] OF REF Edge;

  EDirTriplet = ARRAY [0..2] OF BITS 8 FOR [0..1];
    (* 
      Vertex "i" of triangle "t" is "t.edge[i].vert[t.edir[i]]".
    *)

(* FRAGMENTS *)

        neg := ConsFragment(, neg)
      ELSIF sb < 0.0d0 THEN
        WITH
          u = CutFragSide(a^, sa, b^, sb),
          v = CutFragSide(a^, sa, c^, sc)
        DO
          pos := ConsFragment(NewFragment(whole, a, u, v, flip), pos);
          neg := ConsFragment(NewFragment(whole, u, b, c, flip), neg);
          neg := ConsFragment(NewFragment(whole, v, u, c, flip), neg);
        END
      ELSIF sb = 0.0d0 AND sc < 0.0d0 THEN
        WITH
          v = CutFragSide(a^, sa, c^, sc)
        DO
          pos := ConsFragment(NewFragment(whole, a, b, v, flip), pos);
          neg := ConsFragment(NewFragment(whole, v, b, c, flip), neg);
        END
      ELSIF sc < 0.0d0 THEN
        WITH
          u = CutFragSide(b^, sb, c^, sc),
          v = CutFragSide(a^, sa, c^, sc)
        DO
          pos := ConsFragment(NewFragment(whole, a, b, v, flip), pos);
          pos := ConsFragment(NewFragment(whole, b, u, v, flip), pos);
          neg := ConsFragment(NewFragment(whole, u, c, v, flip), neg);
        END
      ELSE
        pos := ConsFragment(, pos);
      END
    END

PROCEDURE SortVertices(
    VAR a: REF R3.T; VAR sa: LONGREAL;
    VAR b: REF R3.T; VAR sb: LONGREAL;
    VAR c: REF R3.T; VAR sc: LONGREAL;
    VAR flip: BOOLEAN;
  ) =
  VAR t: R3.T; st: LONGREAL;
  BEGIN
    flip := FALSE;
    IF sa < sb THEN
      t := a; st := sa;
      a := b; sa := sb;
      b := t; sb := st;
      flip := NOT flip;
    END;
    IF sb < sc THEN
      t := b; st := sb;
      b := c; sb := sc;
      c := t; sc := st;
      flip := NOT flip;
    END;
    IF sa < sb THEN
      t := a; st := sa;
      a := b; sa := sb;
      b := t; sb := st;
      flip := NOT flip;
    END;
  END SortVertices;

PROCEDURE NewFragment(
    whole: REF Triangle;
    READONLY a, b, c: REF R3.T; 
    flip: BOOLEAN;
  ): REF Fragment =
  BEGIN
    WITH 
      r = NEW(REF Fragment),
       = r^
    DO
      .whole := whole;
      .a := a;
      IF flip THEN 
        .b := c; .c := b
      ELSE
        .b := b; .c := c
      END;
      RETURN r
    END;
  END NewFragment;
      
PROCEDURE ConsFragment(
    : REF Fragment;
    rest: FragmentList;
  ): FragmentList =
  BEGIN
    WITH 
      lst = NEW(FragmentList)
    DO
      lst. := ;
      lst.rest := rest;
      RETURN lst
    END
  END ConsFragment;


    (* Initialize lists for region = whole space. *)
    (* Beware of aliasing between "s" and "in"/"bd"/"ot": *)
    VAR t: Scene := s;
    BEGIN 
      bd := Scene{NIL,NIL,NIL};
      ot := Scene{NIL,NIL,NIL};
      ls := Scene{NIL,NIL,NIL}; 
      in := t
    END;

PROCEDURE Visible(...) =
  BEGIN
    (* Take each triangle in "op", and filter scene triangles against it: *)
    WHILE op.t # NIL DO
      WITH obstacle = op.t DO
        (* Remove manually the pieces of "sc" that are coplanar with "obstacle": *)
        self := EGScene.T{NIL, NIL, NIL};
        othr := EGScene.T{NIL, NIL, NIL};
        
        (* Pieces of op.t => selfS, rest => othrS *)
        
        VAR s: VertexList := sc.v; r: VertexList;
        BEGIN
          sc.v := NIL;
          WHILE s # NIL DO
            r := s.rest;
            IF s.v.org = obstacle.org THEN
              s.rest := self.v; self.v := s
            ELSE
              s.rest := othr.v; othr.v := s
            END;
            s := r
          END;
        END;
        
        VAR s: EdgeList := sc.e; r: EdgeList;
        BEGIN
          sc.e := NIL;
          WHILE s # NIL DO
            r := s.rest;
            IF s.e.org = obstacle.org THEN
              s.rest := self.e; self.e := s
            ELSE
              s.rest := othr.e; othr.e := s
            END;
            s := r
          END;
        END;
        
        VAR s: TriangleList := sc.t; r: TriangleList;
        BEGIN
          sc.t := NIL;
          WHILE s # NIL DO
            r := s.rest;
            IF s.t.org = obstacle.org THEN
              s.rest := self.t; self.t := s
            ELSE
              s.rest := othr.t; othr.t := s
            END;
            s := r
          END;
        END;
        
        (* Clip the region by the shadow volume: *)
        otS := Scene{NIL, NIL, NIL};
        TriShadow(othr, self, obstacle, obs, inS, bdS, otS, inF := FALSE);
        sc := otS; TriangClip.CloseScene(sc)
      END;
      
    END
  END Visible;