(* Created 1993 by Rober M. Rosi *)

MODULE Triang;

IMPORT 
  Oct, R3, LR3, LR3Extras, FileRd, FileWr, Random, Math, Thread, 
  Rd, Wr, Fmt, Lex, FloatMode, OSError, Text, TextRd, TextWr;
FROM Oct IMPORT 
  RBits, RotBits, Sym, Onext, Lnext, Rprev, Rot, Tor, Flip, Oprev, 
  Splice;
  
REVEAL
  Edge = PublicEdge BRANDED OBJECT
    org: ARRAY RBits OF Node;
  OVERRIDES
    init := EdgeInit;
  END;

REVEAL 
  Node = PublicNode BRANDED OBJECT 
    END;

REVEAL     
  Vertex = PublicVertex BRANDED OBJECT 
    END;

REVEAL 
  Face = PublicFace BRANDED OBJECT
    END;
    
(* === INIT METHODS === *)

PROCEDURE EdgeInit(e: Edge): Edge=
  BEGIN
    EVAL NARROW(e,Oct.Edge).init(); 
    e.org[0] := NIL; (* For now *)
    e.org[1] := NIL; (* For now *)
    e.org[2] := NIL; (* For now *)
    e.org[3] := NIL; (* For now *)
    RETURN  e;
  END EdgeInit;

(* === ELEMENT CREATION === *)
  
PROCEDURE MakeEdge(): Arc =
  BEGIN
    WITH e = NEW(Edge).init() DO
      RETURN Arc{edge := e, bits := 0};
    END;
  END MakeEdge;
  
PROCEDURE MakeVertex(): Vertex =
  BEGIN RETURN NEW(Vertex) END MakeVertex;

PROCEDURE MakeFace(): Face =
  BEGIN RETURN NEW(Face) END MakeFace;

(* === ARC PROPERTIES === *)

PROCEDURE Org(e: Arc): Node =
  BEGIN
    WITH c = NARROW(e.edge, Edge).org[RotBits(e)] DO
      RETURN c;
    END;
  END Org;

PROCEDURE SetOrg(e: Arc; n: Node) =
  BEGIN
    WITH c = NARROW(e.edge, Edge).org[RotBits(e)] DO
      c := n
    END;
  END SetOrg;

PROCEDURE SetAllOrgs(a: Arc; n: Node) =
  VAR t: Arc := a;
  BEGIN
    REPEAT
      SetOrg(t, n);
      t := Onext(t)
    UNTIL t = a
  END SetAllOrgs;

PROCEDURE Left(a: Arc): Node =
  BEGIN
    RETURN Org(Tor(a))
  END Left;
   
PROCEDURE SetLeft(e: Arc; n: Node) =
  BEGIN
    SetOrg(Tor(e), n)
  END SetLeft;

PROCEDURE SetAllLefts(a: Arc; n: Node) =
  VAR t: Arc := a;
  BEGIN
    REPEAT
      SetLeft(t, n);
      t := Lnext(t)
    UNTIL t = a
  END SetAllLefts;
  
PROCEDURE OrgV(a: Arc): Vertex =
  BEGIN
    RETURN NARROW(Org(a), Vertex)
  END OrgV;
   
PROCEDURE LeftF(a: Arc): Face =
  BEGIN
    RETURN NARROW(Org(Tor(a)), Face)
  END LeftF;
   
(* === CONSTRUCTION TOOLS === *)

PROCEDURE Connect(a, b: Arc): Arc =
  VAR e: Arc;
  BEGIN
    e := MakeEdge();
    Splice(e, a);
    Splice(Sym(e), b);
    RETURN e;
  END Connect;
  
PROCEDURE AddSpokes(a: Arc): Arc =
  VAR e, b, c, d: Arc;
  BEGIN
    WITH v = MakeVertex() DO 
      e := MakeEdge();
      Splice(e, a); 
      SetOrg(Sym(e), v);
      SetOrg(e, Org(a));

      b := a;
      d := e;
      c := Lnext(a);
      WHILE c # e DO
        WITH
          r = Connect(c, Sym(d)),
          f = MakeFace()
        DO
          SetOrg(r, Org(c));
          SetOrg(Sym(r), v);
          SetLeft(r, f);
          SetLeft(Sym(d), f);
          SetLeft(b, f);
          b := c;
          d := r;
          c := Lnext(c)
        END
      END;

      WITH
        f = MakeFace()
      DO
        SetLeft(e, f);
        SetLeft(Sym(d), f);
        SetLeft(b, f);
      END;

      RETURN e
    END
  END AddSpokes;

PROCEDURE MakeGrid(nx, ny: CARDINAL): ARRAY [0..7] OF Arc =
      
  VAR EdgeCount: CARDINAL := 0;
  
  PROCEDURE MakeCell(s: Arc): Arc =
    (* 
      Adds four triangles nested in the corner between "s" 
      (on the west side, pointing north) and "Oprev(s)" 
      (on the south side, pointing east).
      Returns the east edge of the new cell, pointing north. *)
    VAR b, c, t: Arc;
        v: Vertex;
    BEGIN
      v := MakeVertex();
      t := Sym(Oprev(s));

      b := MakeEdge(); b.edge.num := EdgeCount; INC(EdgeCount);
      Splice(Oprev(t), b);
      SetOrg(b, Org(t));
      SetOrg(Sym(b), v);
      
      c := MakeEdge(); c.edge.num := EdgeCount; INC(EdgeCount);
      Splice(c, Sym(b));
      Splice(Sym(s), Sym(c));
      SetOrg(c, v);
      SetOrg(Sym(c), Org(Sym(s)));

      EVAL AddSpokes(b);

      RETURN b
    END MakeCell;
  
  PROCEDURE MakeCellRow(a: Arc): Arc =
  (* 
    Adds a new row of cells, bounded on the west side by arc "a"
    (assumed to be pointing north).  Returns the east edge of the 
    easternmost cell, pointing north. *)
    BEGIN
      FOR col := 0 TO nx-1 DO
        a := MakeCell(a);
      END;
      RETURN a
    END MakeCellRow;

  VAR b, t: Arc;
      ca: ARRAY [0..7] OF Arc;
  BEGIN
    (* Create bottom row of edges: *)  
    FOR col := 0 TO nx-1 DO
      WITH c = MakeEdge() DO
        c.edge.num := EdgeCount; INC(EdgeCount);
        IF col = 0 THEN
          WITH v = MakeVertex() DO
            SetOrg(c, v);
          END;
          t := c;
        ELSE
          Splice(Sym(b), c);
          SetOrg(c, Org(Sym(b)));
        END;
        b := c;
        WITH w = MakeVertex() DO
          SetOrg(Sym(b), w)
        END;
      END
    END;

    ca[1] := Flip(t);
    ca[2] := Sym(b);
    FOR row := 0 TO ny-1 DO
      (* "t" is the first edge on the top of the last row of cells *)
      (* "b" is the last edge on the top of the last row of cells *)
      (* Both are pointing from left to right *)
      (* Add another row of cells, beginning just above "t": *)
      WITH a = MakeEdge() DO
        a.edge.num := EdgeCount; INC(EdgeCount);
        IF row = 0 THEN ca[0] := a END;
        IF row = ny-1 THEN ca[7] := Flip(Sym(a)) END;

        Splice(a, t);
        SetOrg(a, Org(t));
        WITH u = MakeVertex() DO
          SetOrg(Sym(a), u)
        END;
        WITH f = MakeCellRow(a) DO
          t := Oprev(Sym(a));
          b := Sym(Onext(Sym(f)));

          IF row = 0 THEN ca[3] := Flip(f) END;
          IF row = ny-1 THEN ca[4] := Sym(f) END;
        END
      END
    END;
    ca[6] := t;
    ca[5] := Flip(Sym(b));
    RETURN ca
  END MakeGrid;
  
PROCEDURE GridCell(a: Arc; READONLY x, y: CARDINAL): Arc =
  VAR c: Arc;
  BEGIN
    c := a;
    FOR i := 0 TO x-1 DO
      c := Oprev(Oprev(Oprev(Oprev(Sym(c)))))
    END;
    FOR j := 0 TO y-1 DO 
      c := Onext(Onext(Sym(Onext(Onext(c)))))
    END;
    RETURN c;
  END GridCell;
  
PROCEDURE EnumGridVertices(
    a: Arc; nx, ny: CARDINAL;
    visit: GridVisitProc;
  ) =
  VAR e: Arc;
  BEGIN
    FOR iy := 0 TO ny-1 DO 
      e := a;
      FOR ix := 0 TO nx-1 DO
        visit(e, 2*ix, 2*iy);
        visit(Sym(Onext(e)), 2*ix+1, 2*iy+1);
        IF iy = ny-1 THEN
          visit(Sym(Onext(Onext(e))), 2*ix, 2*iy+2)
        END;
        IF ix = nx-1 THEN
          visit(Sym(e), 2*ix+2, 2*iy)
        END;
        IF ix = nx-1 AND iy = ny-1 THEN 
          visit(Sym(Oprev(Oprev(Sym(e)))), 2*ix+2, 2*iy+2)
        END;
        IF ix < nx-1 THEN
          e := Oprev(Oprev(Oprev(Oprev(Sym(e)))))
        END
      END;
      IF iy < ny-1 THEN
        a := Onext(Onext(Sym(Onext(Onext(a)))))
      END
    END
  END EnumGridVertices;

PROCEDURE SetQuarterGridProperties(
    c: Arc; 
    order: CARDINAL; 
    vertexColor: R3.T;
    vertexRadius: REAL;
    edgeColor: R3.T;
    edgeRadius: REAL;
    facePatch: CARDINAL;
    faceColor: R3.T;
    faceTransp: R3.T;
  ) =
   
  PROCEDURE SetCornerVertex(v: Vertex) =
    BEGIN
      v.color := vertexColor;
      v.radius := vertexRadius;
    END SetCornerVertex;
  
  PROCEDURE SetQuarterEdge(a: Arc) =
    BEGIN
      FOR i := 0 TO order DO 
        IF i # 0 THEN 
          WITH v = OrgV(a) DO
            v.color := edgeColor;
            v.radius := edgeRadius;
          END;
        END;
        WITH e = NARROW(a.edge, Edge) DO 
          e.color := edgeColor;
          e.radius := edgeRadius
        END;
        a := Onext(Onext(Sym(a)));
        IF i MOD 2 # 0 THEN a := Onext(Onext(a)) END;
      END;
    END SetQuarterEdge;
    
  PROCEDURE SetLeftTriangle(t: Face) =
    BEGIN
      t.patch := facePatch;
      t.color := faceColor;
      t.transp := faceTransp;
    END SetLeftTriangle;
  
  PROCEDURE SetRowPatches(d: Arc; row: CARDINAL) =
    BEGIN
      d := Sym(d);
      SetLeftTriangle(Left(d));
      d := Onext(d);
      FOR col := 1 TO order-row-1 DO 
        d := Onext(d);
        SetLeftTriangle(Left(d));
        d := Sym(d);
        SetLeftTriangle(Left(d));
        d := Onext(Flip(d));
      END
    END SetRowPatches;
    
  BEGIN
    SetCornerVertex(Org(c));
    SetQuarterEdge(Oprev(c));
    FOR row := 0 TO order-1 DO 
      SetRowPatches(c, row);
      c := Oprev(Sym(c))
    END
  END SetQuarterGridProperties;
 
PROCEDURE Glue(a, b: Arc; n: CARDINAL): Arc =
  VAR ta, tb, pa: Arc;
  BEGIN
    (* Sanity check: *)
    <* ASSERT n >= 1 *>
    ta := a; tb := b;
    FOR i := 2 TO n DO 
      ta := Lnext(ta); tb := Rprev(tb);
      <* ASSERT ta # a *>
      <* ASSERT tb # b *>
    END;
    
    Splice(a, Oprev(b));
    SetAllOrgs(a, Org(a));
    
    FOR i := 1 TO n DO
      ta := Lnext(a); tb := Rprev(b);
      Splice(ta, Oprev(tb));
      <* ASSERT Onext(a) = b *>
      <* ASSERT Onext(Sym(b)) = Sym(a) *>
      SetAllOrgs(ta, Org(ta));
      (* Disconnect b: *)
      WITH f = Left(b) DO
        Splice(b, Oprev(b));
        Splice(Sym(b), Oprev(Sym(b)));
        SetLeft(a, f)
      END;
      pa := a; 
      a := ta; b := tb
    END;
    
    RETURN pa
  END Glue;

(* === GLOBAL PROCEDURES === *)

PROCEDURE NumberVertices(a: Arc): CARDINAL = 
  VAR n: CARDINAL := 0;
  
  PROCEDURE Visit(c: Arc) =
    VAR cn: Arc;  
    BEGIN
      WITH v = Org(c) DO
        cn := c;
        REPEAT
          WITH vn = Org(cn) DO
            <* ASSERT vn = v *>
            vn.num := n
          END;
          cn := Onext(cn);
        UNTIL (cn = c);
      END;
      INC(n);
    END Visit;
    
  BEGIN
    Oct.EnumVertices(a, Visit);
    RETURN n
  END NumberVertices;

PROCEDURE MakeTopology(a: Arc): Topology =
  VAR top: Topology;
  
  BEGIN
    top.NV := NumberVertices(a); 
    top.NF := NumberVertices(Rot(a));

    top.edge := Oct.NumberEdges(ARRAY OF Arc{a});
    top.NE := NUMBER(top.edge^);
    
    top.out := NEW(REF ARRAY OF Arc, top.NV); 
    top.vertex := NEW(REF ARRAY OF Vertex, top.NV);
    top.face := NEW(REF ARRAY OF Face, top.NF);
    top.side := NEW(REF ARRAY OF Arc, top.NF); 

    FOR ei := 0 TO top.NE-1 DO
      VAR c: Arc := top.edge[ei];
      BEGIN
        FOR k := 0 TO 1 DO
          WITH
            v = Org(c), vi = v.num,
            f = Left(c), fi = f.num 
          DO
            top.vertex[vi] := v; top.out[vi] := c;
            top.face[fi] := f; top.side[fi] := Tor(c);
          END;
          c := Sym(c)
        END
      END
    END;
    RETURN top
  END MakeTopology;

PROCEDURE MakeAdjacencyMatrix(READONLY top: Topology): REF AdjacencyMatrix =
  VAR m := NEW(REF ARRAY OF ARRAY OF BOOLEAN, top.NV, top.NV);
  BEGIN
    FOR i := 0 TO top.NV-1 DO 
      WITH mi = m[i] DO 
        FOR j := 0 TO top.NV-1 DO 
          mi[j] := FALSE
        END
      END
    END;
    FOR ei := 0 TO top.NE-1 DO
      WITH a = top.edge[ei] DO
        WITH i = Org(a).num, j = Org(Sym(a)).num DO
          m[i,j] := TRUE;
          m[j,i] := TRUE;
        END
      END
    END;
    RETURN m
  END MakeAdjacencyMatrix;

PROCEDURE CheckOutAndSide(READONLY top: Topology) =
  BEGIN
    FOR i := 0 TO top.NV-1 DO 
      WITH
        v = top.vertex[i],
        e = top.out[i]
      DO
        <* ASSERT v.num = i *>
        <* ASSERT Org(e) = v *>
      END
    END;
    FOR i := 0 TO top.NF-1 DO
      WITH
        f = top.face[i],
        e = top.side[i]
      DO
        <* ASSERT f.num = i *>
        <* ASSERT Org(e) = f *>
      END
    END
  END CheckOutAndSide;
  
PROCEDURE GetVariableVertices(READONLY top: Topology; VAR vr: ARRAY OF BOOLEAN) =
  BEGIN
    FOR i := 0 TO top.NV-1 DO
      WITH v = top.vertex[i] DO 
        vr[i] := v.exists AND (NOT v.fixed) 
      END;
    END;
  END GetVariableVertices;

(* === GEOMETRIC TOOLS === *)

PROCEDURE InitCoords(
    coins: Random.T; 
    VAR c: Coords; 
    radius: REAL := 1.0
  ) =
  BEGIN
    WITH r = FLOAT(radius, LONGREAL) DO
      FOR i := 0 TO LAST(c) DO
        c[i] := LR3.T{
          coins.longreal(-r, r),
          coins.longreal(-r, r),
          coins.longreal(-r, r) 
        }
      END
    END
  END InitCoords;

PROCEDURE Barycenter(READONLY top: Topology; READONLY c: Coords): LR3.T =
  VAR B: LR3.T := LR3.T{0.0d0, ..};
      N: CARDINAL := 0;
  BEGIN
    FOR i := 0 TO LAST(c) DO 
      WITH v = top.vertex[i] DO
        IF v.exists THEN B := LR3.Add(B, c[i]); INC(N) END;
      END;
    END;
    RETURN LR3.Scale(1.0d0/FLOAT(N, LONGREAL), B)
  END Barycenter;

PROCEDURE MeanVertexNorm(READONLY top: Topology; READONLY c: Coords): LONGREAL =
  VAR S: LONGREAL := 0.0d0;
      N: CARDINAL := 0;
  BEGIN
    FOR i := 0 TO LAST(c) DO 
      WITH v = top.vertex[i] DO
        IF v.exists THEN 
          S := S + LR3.NormSqr(c[i]);
          INC(N)
        END;
      END;
    END;
    RETURN Math.sqrt(S/FLOAT(N,LONGREAL))
  END MeanVertexNorm;

PROCEDURE MeanEdgeLength(READONLY top: Topology; READONLY c: Coords): LONGREAL =
  VAR S: LONGREAL := 0.0d0;
      N: CARDINAL := 0;
  BEGIN
    FOR i := 0 TO top.NE-1 DO 
      WITH a = top.edge[i], e = NARROW(a.edge, Edge) DO
        IF e.exists THEN 
          WITH o = OrgV(a).num, d = OrgV(Sym(a)).num DO 
            S := S + LR3.DistSqr(c[o], c[d])
          END;
          INC(N)
        END;
      END;
    END;
    RETURN Math.sqrt(S/FLOAT(N,LONGREAL))
  END MeanEdgeLength;

PROCEDURE Displace(READONLY top: Topology; d: LR3.T; VAR c: Coords) =
  BEGIN
    FOR i := 0 TO LAST(c) DO 
      IF top.vertex[i].exists THEN
        WITH vc = c[i] DO vc := LR3.Add(vc, d) END
      END
    END
  END Displace;

PROCEDURE Scale(READONLY top: Topology; s: LONGREAL; VAR c: Coords) =
  BEGIN
    FOR i := 0 TO LAST(c) DO 
      IF top.vertex[i].exists THEN
        WITH vc = c[i] DO vc := LR3.Scale(s, vc) END;
      END
    END
  END Scale;

PROCEDURE NormalizeVertexNorms(READONLY top: Topology; VAR c: Coords) =
  BEGIN
    WITH b = Barycenter(top, c) DO Displace(top, LR3.Neg(b), c) END;
    WITH s = MeanVertexNorm(top, c) DO Scale(top, 1.0d0/s, c) END;
  END NormalizeVertexNorms;

PROCEDURE NormalizeEdgeLengths(READONLY top: Topology; VAR c: Coords) =
  BEGIN
    WITH b = Barycenter(top, c) DO Displace(top, LR3.Neg(b), c) END;
    WITH s = MeanEdgeLength(top, c) DO Scale(top, 1.0d0/s, c) END;
  END NormalizeEdgeLengths;

PROCEDURE FaceCross(a: Arc; READONLY c: Coords): LR3.T =
  VAR be: Arc; bo, n: LR3.T; bn: CARDINAL; 
  BEGIN
    IF NOT LeftF(a).exists THEN
      RETURN LR3.T{1.0d0, 0.0d0, 0.0d0}
    ELSE
      WITH
        on = OrgV(a).num,
        o = c[on],
        an = OrgV(Sym(a)).num 
      DO
        be := a;
        bn := an;
        bo := LR3.Sub(c[bn], o);
        LOOP
          WITH
            ce = Lnext(be),
            cn = OrgV(Sym(be)).num
          DO
            IF cn = on THEN EXIT END;
            WITH 
              co = LR3.Sub(c[cn], o),
              r = LR3Extras.Cross(bo, co)
            DO
              IF bn = an THEN n := r ELSE n := LR3.Add(n, r) END;
              bn := cn; bo := co; be := ce
            END
          END
        END
      END;
      RETURN n
    END
  END FaceCross;

PROCEDURE FaceNormal(a: Arc; READONLY c: Coords): LR3.T =
  BEGIN
    WITH
      n = FaceCross(a, c),
      s = LR3.Norm(n)
    DO
      IF s = 0.0d0 THEN 
        RETURN LR3.T{1.0d0, 0.0d0, 0.0d0}
      ELSE
        RETURN LR3.Scale(1.0d0/s, n)
      END
    END
  END FaceNormal;

PROCEDURE FaceBarycenter(a: Arc; READONLY c: Coords): LR3.T =
  BEGIN
    WITH
      e = Lnext(a),
      i = Lnext(e),
      ac = c[OrgV(a).num],
      ec = c[OrgV(e).num],
      ic = c[OrgV(i).num],
      bar = LR3.Scale(1.0d0/3.0d0, LR3.Add(LR3.Add(ac, ec), ic))
    DO
      RETURN bar
    END
  END FaceBarycenter;

PROCEDURE VertexCross(a: Arc; READONLY c: Coords): LR3.T =
  VAR ao: Arc; 
      sum: LR3.T := LR3.T{0.0d0, ..};
  BEGIN
    WITH
      uv = OrgV(a),
      u = c[uv.num]
    DO
      IF NOT uv.exists THEN
        RETURN LR3.T{1.0d0, 0.0d0, 0.0d0}
      ELSE
        ao := a;
        REPEAT
          WITH
            an = Onext(ao),
            ov = OrgV(Sym(ao)), (* DestV(ao) *)  
            nv = OrgV(Sym(an))  (* DestV(an) *)
          DO
            IF ov.exists AND nv.exists THEN
              WITH 
                o = c[ov.num],
                n = c[nv.num],
                r = LR3Extras.Cross(LR3.Sub(o, u), LR3.Sub(n, u))
              DO
                IF ao = a THEN sum := r ELSE sum := LR3.Add(sum, r) END
              END
            END;
            ao := an
          END;
        UNTIL (ao = a);
        RETURN sum
      END
    END
  END VertexCross;

PROCEDURE VertexNormal(a: Arc; READONLY c: Coords): LR3.T =
  BEGIN
    WITH
      n = VertexCross(a, c),
      s = LR3.Norm(n)
    DO
     IF s = 0.0d0 THEN 
        RETURN LR3.T{1.0d0, 0.0d0, 0.0d0}
      ELSE
        RETURN LR3.Scale(1.0d0/s, n)
      END
    END
  END VertexNormal;

PROCEDURE NeighborBarycenter(a: Arc; READONLY c: Coords): LR3.T =
  VAR an: Arc;
      n: CARDINAL := 0;
      b := LR3.T{0.0d0, ..};
  BEGIN
    an := a;
    REPEAT
      WITH v = OrgV(Sym(an)) DO
        IF v.exists THEN b := LR3.Add(b, c[v.num]); INC(n) END
      END;
      an := Onext(an)
    UNTIL (an = a);
    IF n = 0 THEN
      RETURN b
    ELSE
      RETURN LR3.Scale(1.0d0/FLOAT(n, LONGREAL), b)
    END;
  END NeighborBarycenter;

PROCEDURE ComputeAllVertexNormals(
    READONLY top: Topology; 
    READONLY c: Coords;
  ): REF ARRAY OF LR3.T =
  BEGIN
    WITH
      rvn = NEW(REF ARRAY OF LR3.T, top.NV),
      vn = rvn^
    DO
      FOR i := 0 TO top.NV-1 DO 
        vn[i] := VertexNormal(top.out[i], c)
      END;
      RETURN rvn
    END;
  END ComputeAllVertexNormals; 

(* === INPUT/OUTPUT === *)

CONST
  Boole = ARRAY BOOLEAN OF CHAR {'F', 'T'};
  AlphaChars = SET OF CHAR{'\t', '\n', '\r', '\f', ' ', 'A'..'Z', 'a'..'z'};
  
PROCEDURE Read(name: TEXT): Configuration =
 
  PROCEDURE ReadBool(rd: Rd.T): BOOLEAN RAISES {Lex.Error} =
    <* FATAL Rd.Failure, Rd.EndOfFile, Thread.Alerted *>
    BEGIN
      WITH c = Rd.GetChar(rd) DO
        IF c = 'T' THEN 
          RETURN TRUE
        ELSIF c = 'F' THEN
          RETURN FALSE
        ELSE 
          Rd.UnGetChar(rd);
          RAISE Lex.Error
        END
      END
    END ReadBool;

  <* FATAL Rd.Failure, Rd.EndOfFile, Thread.Alerted, FloatMode.Trap, Lex.Error *>
  <* FATAL OSError.E *>
  VAR top: Topology; 
      c: REF Coords;
      comments: TEXT;
  BEGIN
    WITH 
      rd = FileRd.Open(name & ".top")
    DO
      (* Header comments: *)
      comments := ReadComments(rd);

      (* Element counts: *)
      Lex.Skip(rd, cs := AlphaChars);
      top.NV := Lex.Int(rd);
      Lex.Skip(rd, cs := AlphaChars);
      top.NF := Lex.Int(rd);
      Lex.Skip(rd, cs := AlphaChars);
      top.NE := Lex.Int(rd);
      Lex.Skip(rd);

      WITH
        map = NEW(REF ARRAY OF Oct.Edge, top.NE)^
      DO
        (* Create edges and build edge map: *)
        top.edge := NEW(REF ARRAY OF Arc, top.NE);
        FOR i := 0 TO top.NE-1 DO 
          top.edge[i] := MakeEdge();
          top.edge[i].edge.num := i;
          map[i] := top.edge[i].edge
        END;

        (* Create vertex records: *)
        top.vertex := NEW(REF ARRAY OF Vertex, top.NV);
        top.out := NEW(REF ARRAY OF Arc, top.NV);
        c := NEW(REF Coords, top.NV);
        FOR i := 0 TO top.NV-1 DO top.vertex[i] := MakeVertex(); END;

        (* Create face records: *)
        top.face := NEW(REF ARRAY OF Face, top.NF);
        top.side := NEW(REF ARRAY OF Arc, top.NF);
        FOR i := 0 TO top.NF-1 DO top.face[i] := MakeFace() END;


        (* Read vertex data: *)
        EVAL ReadComments(rd);
        FOR j := 0 TO top.NV-1 DO 
          Lex.Skip(rd);
          WITH 
            nv = Lex.Int(rd),
            v = top.vertex[nv]
          DO 
            <* ASSERT nv = j *>
            v.num := nv;
            Lex.Skip(rd); 
            v.exists := ReadBool(rd);
            Lex.Skip(rd);
            v.fixed := ReadBool(rd);
            Lex.Skip(rd);
            WITH cv = c[nv] DO
              cv[0] := Lex.LongReal(rd);
              Lex.Skip(rd);
              cv[1] := Lex.LongReal(rd);
              Lex.Skip(rd);
              cv[2] := Lex.LongReal(rd);
            END;
            Lex.Skip(rd);
            WITH cc = v.color DO
              cc[0] := Lex.Real(rd);
              Lex.Skip(rd);
              cc[1] := Lex.Real(rd);
              Lex.Skip(rd);
              cc[2] := Lex.Real(rd); 
            END;
            Lex.Skip(rd);
            v.radius := Lex.Real(rd);
          END;      
        END;
        Lex.Skip(rd);

        (* Read face data: *)
        EVAL ReadComments(rd);
        FOR j := 0 TO top.NF-1 DO 
          Lex.Skip(rd);
          WITH 
            nf = Lex.Int(rd),
            f = top.face[nf]
          DO
            <* ASSERT nf = j *>
            f.num := nf;
            Lex.Skip(rd);
            f.exists := ReadBool(rd);
            Lex.Skip(rd);
            WITH cc = f.color DO
              cc[0] := Lex.Real(rd);
              Lex.Skip(rd);
              cc[1] := Lex.Real(rd);
              Lex.Skip(rd);
              cc[2] := Lex.Real(rd); 
            END;
            Lex.Skip(rd);
            WITH cc = f.transp DO
              cc[0] := Lex.Real(rd);
              Lex.Skip(rd);
              cc[1] := Lex.Real(rd);
              Lex.Skip(rd);
              cc[2] := Lex.Real(rd); 
            END;
            Lex.Skip(rd);
            f.patch := Lex.Int(rd);
          END;
        END;
        Lex.Skip(rd);

        (* Read edge data: *)
        EVAL ReadComments(rd);
        FOR j := 0 TO top.NE-1 DO 
          Lex.Skip(rd);
          WITH 
            ne = Lex.Int(rd),
            e = NARROW(top.edge[ne].edge, Edge)
          DO
            <* ASSERT ne = j *>
            <* ASSERT top.edge[ne].bits = 0 *>
            Oct.ReadEdge(rd, e, map);
            FOR k := 0 TO 3 DO 
              Lex.Skip(rd);
              WITH 
                n = Lex.Int(rd),
                vf = e.org[k]
              DO 
                IF Rd.GetChar(rd) = 'v' THEN 
                  vf := top.vertex[n];
                  top.out[n] := Arc{edge := e, bits := 2*k};
                ELSE
                  vf := top.face[n];
                  top.side[n] := Arc{edge := e, bits := 2*k};                  
                END;
              END;
            END;
            Lex.Skip(rd);
            e.exists := ReadBool(rd);
            Lex.Skip(rd);
            e.spring := ReadBool(rd);
            Lex.Skip(rd);
            WITH ec = e.color DO
              ec[0] := Lex.Real(rd);
              Lex.Skip(rd);
              ec[1] := Lex.Real(rd);
              Lex.Skip(rd);
              ec[2] := Lex.Real(rd)
            END;
            Lex.Skip(rd);
            e.radius := Lex.Real(rd); 
          END;
        END;
        CheckOutAndSide(top);
        Rd.Close(rd); 
        RETURN Configuration{top, c, comments}
      END;
    END;  
  END Read;

PROCEDURE Write(
    name: TEXT; 
    READONLY top: Topology; 
    READONLY c: Coords;
    comments: TEXT := "";
  ) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  <* FATAL OSError.E *>
  BEGIN
    WITH 
      wr = FileWr.Open(name & ".top"),
      vWidth = NumDigits(top.NV - 1),
      fWidth = NumDigits(top.NF - 1),
      eWidth = NumDigits(top.NE - 1)
    DO 

      PROCEDURE WriteCoord(x: LONGREAL) =
        BEGIN
          Wr.PutText(wr, Fmt.Pad(Fmt.LongReal(x, Fmt.Style.Sci, prec := 7), 14));
        END WriteCoord;
        
      PROCEDURE WritePoint(READONLY c: LR3.T) =
        BEGIN
          WriteCoord(c[0]);
          Wr.PutText(wr, " ");
          WriteCoord(c[1]);
          Wr.PutText(wr, " ");
          WriteCoord(c[2]);
        END WritePoint;
        
      PROCEDURE WriteIntensity(r: REAL) =
        BEGIN
          Wr.PutText(wr, Fmt.Real(r, Fmt.Style.Fix, prec := 4));
        END WriteIntensity;

      PROCEDURE WriteColor(READONLY c: R3.T) =
        BEGIN
          WriteIntensity(c[0]);
          Wr.PutText(wr, " ");                     
          WriteIntensity(c[1]);
          Wr.PutText(wr, " ");                     
          WriteIntensity(c[2]);
        END WriteColor;
        
      PROCEDURE WriteRadius(r: REAL) =
        BEGIN
          Wr.PutText(wr, Fmt.Real(r, prec := 4));
        END WriteRadius;
        
      BEGIN
        IF NOT Text.Empty(comments) THEN
          WriteComments(wr, comments & "\n")
        END;
        Wr.PutText(wr, "vertices "); Wr.PutText(wr, Fmt.Int(top.NV) & "\n");  
        Wr.PutText(wr, "faces "); Wr.PutText(wr, Fmt.Int(top.NF) & "\n"); 
        Wr.PutText(wr, "edges "); Wr.PutText(wr, Fmt.Int(top.NE) & "\n");

        WriteComments(wr, "\nVertex data:\n");
        FOR i := 0 TO top.NV-1 DO
          WITH v = top.vertex[i] DO 
            Wr.PutText(wr, Fmt.Pad(Fmt.Int(v.num), vWidth)); 
            Wr.PutText(wr, "  ");
            Wr.PutText(wr, Fmt.Char(Boole[v.exists]));
            Wr.PutText(wr, " ");
            Wr.PutText(wr, Fmt.Char(Boole[v.fixed]));
            Wr.PutText(wr, "  ");
            WritePoint(c[v.num]);
            Wr.PutText(wr, "  ");
            WriteColor(v.color);
            Wr.PutText(wr, "  ");
            WriteRadius(v.radius);
            Wr.PutText(wr, "\n")
          END       
        END;

        WriteComments(wr, "\nFace data:\n");
        FOR i := 0 TO top.NF-1 DO
          WITH f = top.face[i] DO 
            Wr.PutText(wr, Fmt.Pad(Fmt.Int(f.num), fWidth)); 
            Wr.PutText(wr, "  ");
            Wr.PutText(wr, Fmt.Char(Boole[f.exists]));
            Wr.PutText(wr, "  ");
            WriteColor(f.color);
            Wr.PutText(wr, "  ");
            WriteColor(f.transp);
            Wr.PutText(wr, "  ");
            Wr.PutText(wr, Fmt.Pad(Fmt.Int(f.patch), fWidth));
            Wr.PutText(wr, "\n")
          END
        END;

        WriteComments(wr, "\nEdge data:\n");
        FOR i := 0 TO top.NE-1 DO
          WITH e = NARROW(top.edge[i].edge, Edge) DO 
            Wr.PutText(wr, Fmt.Pad(Fmt.Int(e.num), eWidth)); 
            Wr.PutText(wr, "  ");
            Oct.PrintEdge(wr, e, eWidth);
            Wr.PutText(wr, "  ");
            FOR j := 0 TO 3 DO
              WITH n = e.org[j] DO 
                TYPECASE n OF 
                | NULL => <* ASSERT FALSE *>
                | Vertex(v) => Wr.PutText(wr, Fmt.Pad(Fmt.Int(v.num), vWidth) & "v ");
                | Face(f) => Wr.PutText(wr, Fmt.Pad(Fmt.Int(f.num), fWidth) & "f ");
                ELSE <* ASSERT FALSE *>
                END;
              END;
            END;
            Wr.PutText(wr, " ");
            Wr.PutText(wr, Fmt.Char(Boole[e.exists]));
            Wr.PutText(wr, " ");
            Wr.PutText(wr, Fmt.Char(Boole[e.spring]));
            Wr.PutText(wr, "  ");
            WriteColor(e.color);
            Wr.PutText(wr, "  ");
            WriteRadius(e.radius);
            Wr.PutText(wr, "\n");
          END;
        END;       
      END;
      Wr.Close(wr); 
    END;     
  END Write;

PROCEDURE NumDigits(n: CARDINAL): CARDINAL =
  (* Width of "n" when printed. *)
  VAR w: CARDINAL := 1;
  BEGIN
    WHILE n > 9 DO INC(w); n := n DIV 10 END;
    RETURN w
  END NumDigits;

PROCEDURE WritePOV(
    wr: Wr.T;
    READONLY top: Topology; 
    READONLY c: Coords;
    smooth: BOOLEAN;
  ) =
  <* FATAL Wr.Failure, Thread.Alerted *>

  PROCEDURE WritePOVCoord(c: LONGREAL) =
    BEGIN
      Wr.PutText(wr, Fmt.LongReal(c, Fmt.Style.Fix, prec := 4))
    END WritePOVCoord;

  PROCEDURE WritePOVPoint(READONLY p: LR3.T) =
    BEGIN
      Wr.PutText(wr, "<");
      WritePOVCoord(p[1]); 
      Wr.PutText(wr, ",");
      WritePOVCoord(p[2]); 
      Wr.PutText(wr, ",");
      WritePOVCoord(-p[0]);
      Wr.PutText(wr, ">");
    END WritePOVPoint;

  PROCEDURE WritePOVColor(READONLY cr: R3.T) =
    BEGIN
      Wr.PutText(wr, "rgb <");
      Wr.PutText(wr, Fmt.Real(cr[0], Fmt.Style.Fix, 4));
      Wr.PutText(wr, ",");
      Wr.PutText(wr, Fmt.Real(cr[1], Fmt.Style.Fix, 4));
      Wr.PutText(wr, ",");
      Wr.PutText(wr, Fmt.Real(cr[2], Fmt.Style.Fix, 4));
      Wr.PutText(wr, ">");
    END WritePOVColor;

  PROCEDURE WritePOVColorTransp(READONLY cr: R3.T; tr: REAL) =
    BEGIN
      Wr.PutText(wr, "rgbf <");
      Wr.PutText(wr, Fmt.Real(cr[0]));
      Wr.PutText(wr, ",");
      Wr.PutText(wr, Fmt.Real(cr[1]));
      Wr.PutText(wr, ",");
      Wr.PutText(wr, Fmt.Real(cr[2]));
      Wr.PutText(wr, ",");
      Wr.PutText(wr, Fmt.Real(tr));
      Wr.PutText(wr, ">");
    END WritePOVColorTransp;

  PROCEDURE WritePOVCylinder(READONLY o, d: LR3.T; radius: REAL; READONLY cr: R3.T) =
    BEGIN
      Wr.PutText(wr, "  cylinder {\n");
      Wr.PutText(wr, "    "); WritePOVPoint(o); Wr.PutText(wr, ",\n");
      Wr.PutText(wr, "    "); WritePOVPoint(d); Wr.PutText(wr, ",\n");
      Wr.PutText(wr, "    " & Fmt.Real(radius) & "\n");
      Wr.PutText(wr, "    open\n");
      Wr.PutText(wr, "    texture { edgetexture } \n");
      Wr.PutText(wr, "    pigment { color "); WritePOVColor(cr); Wr.PutText(wr, " }\n");
      Wr.PutText(wr, "  }\n");
      Wr.PutText(wr, "\n");
      Wr.Flush(wr);
    END WritePOVCylinder;

  PROCEDURE WritePOVSphere(READONLY p: LR3.T; radius: REAL; READONLY cr: R3.T) =
    BEGIN
      Wr.PutText(wr, "  sphere {\n");
      Wr.PutText(wr, "    "); WritePOVPoint(p);
      Wr.PutText(wr, ", " & Fmt.Real(radius) & "\n");
      Wr.PutText(wr, "    texture { edgetexture } \n");  
      Wr.PutText(wr, "    pigment { color "); WritePOVColor(cr); Wr.PutText(wr, " }\n");
      Wr.PutText(wr, "  } \n");
      Wr.PutText(wr, "\n");
      Wr.Flush(wr);
    END WritePOVSphere;

  PROCEDURE WritePOVTriangle(READONLY a, b, c: LR3.T; READONLY cr: R3.T; tr: REAL) =
    BEGIN
      Wr.PutText(wr, "  triangle {");
      Wr.PutText(wr, "    "); WritePOVPoint(a); Wr.PutText(wr, ",\n");
      Wr.PutText(wr, "    "); WritePOVPoint(b); Wr.PutText(wr, ",\n");
      Wr.PutText(wr, "    "); WritePOVPoint(c); Wr.PutText(wr, "\n");
      Wr.PutText(wr, "    texture { facetexture } \n");
      Wr.PutText(wr, "    pigment { color "); 
        WritePOVColorTransp(cr, tr); 
        Wr.PutText(wr, " }\n");
      Wr.PutText(wr, "  } \n");
      Wr.PutText(wr, "\n");
      Wr.Flush(wr);
    END WritePOVTriangle;
    
  PROCEDURE WritePOVSmoothTriangle(
      READONLY a, an: LR3.T;
      READONLY b, bn: LR3.T;
      READONLY c, cn: LR3.T; 
      READONLY cr: R3.T; 
      tr: REAL
    ) =
    BEGIN
      Wr.PutText(wr, "  smooth_triangle {\n");
      Wr.PutText(wr, "    "); WritePOVPoint(a); Wr.PutText(wr, ", ");
        WritePOVPoint(an); Wr.PutText(wr, ",\n");
      Wr.PutText(wr, "    "); WritePOVPoint(b); Wr.PutText(wr, ", ");
        WritePOVPoint(bn); Wr.PutText(wr, ",\n");
      Wr.PutText(wr, "    "); WritePOVPoint(c); Wr.PutText(wr, ", ");
        WritePOVPoint(cn); Wr.PutText(wr, "\n");
      Wr.PutText(wr, "    texture { facetexture } \n");  
      Wr.PutText(wr, "    pigment { color "); 
        WritePOVColorTransp(cr, tr); 
        Wr.PutText(wr, " }\n");
      Wr.PutText(wr, "  } \n");
      Wr.PutText(wr, "\n");
      Wr.Flush(wr);
    END WritePOVSmoothTriangle;
  
  VAR  rvn: REF ARRAY OF LR3.T;
  BEGIN
    FOR i := 0 TO top.NE-1 DO 
      WITH 
        e = top.edge[i],
        er = NARROW(e.edge, Edge)
      DO 
        IF er.exists AND er.radius > 0.0 THEN 
          WITH o = OrgV(e).num, d = OrgV(Sym(e)).num DO    
            WritePOVCylinder(c[o], c[d], er.radius, er.color)
          END
        END;
      END;
    END;
    
    FOR i := 0 TO top.NV-1 DO 
      WITH v = top.vertex[i] DO 
        IF v.exists AND v.radius > 0.0 THEN 
          WritePOVSphere(c[i], v.radius, v.color)
        END
      END
    END; 
    IF smooth THEN rvn := ComputeAllVertexNormals(top, c) END;
    FOR i := 0 TO top.NF-1 DO 
      WITH f = top.face[i] DO
        IF f.exists THEN 
          WITH 
            ea = Rot(top.side[i]),  an = OrgV(ea).num,
            eb = Lnext(ea),         bn = OrgV(eb).num,
            ec = Lnext(eb),         cn = OrgV(ec).num,
            t3 = f.transp,
            transp = (t3[0] + t3[1] + t3[2]) / 3.0
          DO    
            IF bn # an AND cn # an AND Lnext(ec) = ea THEN
              IF smooth THEN
                WITH
                  vn = rvn^
                DO
                  (* Face is a triangle, we can show it: *)
                  WritePOVSmoothTriangle(
                    c[an], vn[an],
                    c[bn], vn[bn],
                    c[cn], vn[cn],
                    f.color, transp
                  )
                END
              ELSE
                WritePOVTriangle(
                  c[an], c[bn], c[cn],
                  f.color, transp
                )
              END
            END
          END
        END
      END
    END;
    Wr.Flush(wr);
  END WritePOV;
  
PROCEDURE WriteX3D(
    wr: Wr.T; 
    READONLY top: Topology; 
    READONLY c: Coords;
    all: BOOLEAN;
  ) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  
  PROCEDURE WriteX3DCoord(x: LONGREAL) =
    BEGIN
      Wr.PutText(wr, Fmt.LongReal(x*1000.0d0, prec := 5))
    END WriteX3DCoord;
  
  PROCEDURE WriteX3DColor(READONLY cr: R3.T) =
    BEGIN
      Wr.PutText(wr, Fmt.Int(ROUND(255.0*cr[0])));
      Wr.PutText(wr, " ");
      Wr.PutText(wr, Fmt.Int(ROUND(255.0*cr[1])));
      Wr.PutText(wr, " ");
      Wr.PutText(wr, Fmt.Int(ROUND(255.0*cr[2])));
    END WriteX3DColor;
    
  PROCEDURE WriteX3DPoint(READONLY p: LR3.T) =
    BEGIN
      WriteX3DCoord(p[0]);
      Wr.PutText(wr, " ");
      WriteX3DCoord(p[1]);
      Wr.PutText(wr, " ");
      WriteX3DCoord(p[2]);
    END WriteX3DPoint;
    
  BEGIN
    FOR i := 0 TO top.NF-1 DO 
      IF all OR top.face[i].exists THEN 
        WITH
          ec = Rot(top.side[i]),
          ed = Lnext(ec), 
          ee = Lnext(ed), 
          k = OrgV(ec).num,
          l = OrgV(ed).num,
          m = OrgV(ee).num,
          color = LeftF(ec).color 
        DO
          IF ed # ec AND ee # ec AND Lnext(ee) = ec THEN
            (* Face is a triangle, we can show it: *)
            Wr.PutText(wr, "# face " & Fmt.Int(i) & "\n");
            Wr.PutText(wr, "  "); WriteX3DColor(color); Wr.PutText(wr, "# color \n\n");
            Wr.PutText(wr, "  "); WriteX3DPoint(c[k]); Wr.PutText(wr, "\n");
            Wr.PutText(wr, "  "); WriteX3DPoint(c[l]); Wr.PutText(wr, "\n");
            Wr.PutText(wr, "  "); WriteX3DPoint(c[m]); Wr.PutText(wr, "\n");
            Wr.PutText(wr, "\n");
            Wr.Flush(wr);
          END
        END;
      END;
    END;
  END WriteX3D;
  
<*OBSOLETE*>
PROCEDURE ReadVTX(inFile: TEXT): REF Coords =
  VAR NV: CARDINAL;
      rc: REF Coords;
  <* FATAL Rd.Failure, Thread.Alerted, OSError.E *>
  <* FATAL Lex.Error, FloatMode.Trap *>
  BEGIN
    WITH 
      rd = FileRd.Open(inFile & ".vtx")
    DO
      (* Element counts: *)
      Lex.Skip(rd, Lex.Blanks);
      Lex.Match(rd, "vertices"); 
      NV := Lex.Int(rd);
      rc := NEW(REF Coords, NV);
      
      WITH c = rc^ DO
        FOR i := 0 TO NV-1 DO 
          WITH ci = c[i] DO
            Lex.Skip(rd);
            ci[0] := Lex.LongReal(rd);
            Lex.Skip(rd);
            ci[1] := Lex.LongReal(rd);
            Lex.Skip(rd);
            ci[2] := Lex.LongReal(rd);
          END
        END
      END;
      
      Rd.Close(rd);
      RETURN rc
    END;  
  END ReadVTX;
  
PROCEDURE WriteRLuisSt(
    wr: Wr.T;
    READONLY c: Coords;
    READONLY v: Coords;
    prec: CARDINAL := 4;
    comments: TEXT := "";
  ) =
    
  PROCEDURE WriteRLuisCoord(x: LONGREAL) =
    <* FATAL Wr.Failure, Thread.Alerted *>    
    BEGIN
      Wr.PutText(wr, Fmt.LongReal(x, Fmt.Style.Sci, prec := prec))
    END WriteRLuisCoord;

  PROCEDURE WriteRLuisPoint(READONLY p: LR3.T) =   
    <* FATAL Wr.Failure, Thread.Alerted *>    
    BEGIN
      WriteRLuisCoord(p[1]);
      Wr.PutText(wr, " ");
      WriteRLuisCoord(p[0]);
      Wr.PutText(wr, " ");
      WriteRLuisCoord(p[2]);
    END WriteRLuisPoint;

  <* FATAL Wr.Failure, Thread.Alerted *>    
  BEGIN
    WITH
      NV = NUMBER(c),
      d = NumDigits(NV-1)
    DO
      Wr.PutText(wr, "begin state (format of 95-07-30)\n");
      WriteComments(wr, "\n" & comments & "\n");
      Wr.PutText(wr, "vertices = " & Fmt.Int(NV) & "\n");
      Wr.PutText(wr, "t = 0\n");

      FOR i := 0 TO NV-1 DO 
        Wr.PutText(wr, Fmt.Pad(Fmt.Int(i), d) & ": "); 
        WriteRLuisPoint(c[i]); 
        Wr.PutText(wr, "  ");
        IF v[i][0] = 0.0d0 AND v[i][1] = 0.0d0 AND v[i][2] = 0.0d0 THEN
          Wr.PutText(wr, "0 0 0")
        ELSE
          WriteRLuisPoint(v[i])
        END; 
        Wr.PutText(wr, "\n");
      END;

      (* Write to file "name.st" *)
      Wr.PutText(wr, "end state\n");
      Wr.Flush(wr);
    END;
  END WriteRLuisSt;
  
PROCEDURE WriteRLuisSuMa(
    outFile: TEXT;
    READONLY top: Topology;    
    all: BOOLEAN;
    comments: TEXT := "";
  ) =
  <* FATAL Wr.Failure, Thread.Alerted, OSError.E *>
  
  BEGIN
    WITH
      su = FileWr.Open(outFile & ".su"),
      ma = FileWr.Open(outFile & ".ma"),
      NV = top.NV,
      NF = top.NF,
      vd = NumDigits(NV-1),
      fd = NumDigits(NF-1)
    DO
    
      PROCEDURE WriteRLuisIntensity(r: REAL) =
        BEGIN
          Wr.PutText(ma, Fmt.Real(r, Fmt.Style.Fix, prec := 3));
        END WriteRLuisIntensity;

      PROCEDURE WriteRLuisColor(READONLY c: R3.T) =
        BEGIN
          WriteRLuisIntensity(c[0]);
          Wr.PutText(ma, " ");                     
          WriteRLuisIntensity(c[1]);
          Wr.PutText(ma, " ");                     
          WriteRLuisIntensity(c[2]);
        END WriteRLuisColor;

      BEGIN    
        Wr.PutText(su, "begin surface (format of 95-07-30)\n");
        Wr.PutText(su, "triangles = " & Fmt.Int(top.NF) & "\n");
        WriteComments(su, "\n" & comments & "\n");

        Wr.PutText(ma, "begin materials (format of 95-07-30)\n");
        Wr.PutText(ma, "materials = " & Fmt.Int(top.NF) & "\n");
        WriteComments(ma, "\n" & comments & "\n");

        FOR i := 0 TO NF-1 DO 
          WITH f = top.face[i], s = top.side[i] DO
            IF (all OR f.exists) AND Oct.Degree(s) = 3 THEN 
              WITH
                ec = Rot(s),
                ed = Lnext(ec), 
                ee = Lnext(ed), 
                k = OrgV(ec).num,
                l = OrgV(ed).num,
                m = OrgV(ee).num
              DO
                Wr.PutText(su, Fmt.Pad(Fmt.Int(k), vd) & " ");
                Wr.PutText(su, Fmt.Pad(Fmt.Int(l), vd) & " ");
                Wr.PutText(su, Fmt.Pad(Fmt.Int(m), vd) & "  ");
                Wr.PutText(su, Fmt.Pad(Fmt.Int(i), fd) & "  ");  (* material code *)
                Wr.PutText(su, "0 0 0");           
                Wr.PutText(su, "\n");

                Wr.PutText(ma, "0.150 0.150 0.150");             (* Ambient color *)
                Wr.PutText(ma, "  "); WriteRLuisColor(f.color);  (* Diffuse color *)
                Wr.PutText(ma, "  0.300 0.300 0.300");           (* Specular color *)
                Wr.PutText(ma, "  "); WriteRLuisColor(f.transp); (* Transparency *)
                Wr.PutText(ma, "  1");                           (* Index of refraction *)
                Wr.PutText(ma, "  100");                         (* Gloss exponent? *)
                Wr.PutText(ma, "\n");
              END;
            END;
          END
        END;
        Wr.PutText(su, "end surface\n");
        Wr.Close(su);

        Wr.PutText(ma, "end materials\n");
        Wr.Close(ma);
      END
    END
  END WriteRLuisSuMa;

CONST CommentPrefix: CHAR = '|';
EXCEPTION MissingFinalNewLine;

PROCEDURE WriteComments(wr: Wr.T; comments: TEXT) =
  (* 
    Writes the given "comments" text to "wr", with a CommentPrefix
    in front of every line.  Supplies a final '\n' if the text is 
    non-empty but does not end with newline. *)
  
  VAR rd: Rd.T := TextRd.New(comments);
  
  PROCEDURE CopyLine() RAISES {Rd.EndOfFile} =
  (*
    Copy one line from "rd" to "wr", prefixed by CommentPrefix. 
    Supplies a final '\n' if the next line exists but does NOT
    end with newline. Raises Rd.EndOfFile if there are no more
    lines in "rd". *)
    
    <* FATAL Rd.Failure, Wr.Failure, Thread.Alerted *>
    VAR c: CHAR;
    BEGIN
      c := Rd.GetChar(rd); (* If EOF here, propagate to caller *)
      Wr.PutChar(wr, CommentPrefix);
      Wr.PutChar(wr, c);
      WHILE c # '\n' DO
        TRY c := Rd.GetChar(rd) EXCEPT Rd.EndOfFile => c := '\n' END;
        Wr.PutChar(wr, c)
      END
    END CopyLine;

  BEGIN
    TRY LOOP CopyLine() END EXCEPT Rd.EndOfFile => (* Ok *) END;
  END WriteComments;

PROCEDURE ReadComments(rd: Rd.T): TEXT =
  (*
    Reads zero or more lines from "rd" that begin with CommentPrefix, strips the
    leading CommentPrefix of each line, and returns the concatenation of those lines
    as a single TEXT with embedded and terminating newline chars. *)

  VAR wr: Wr.T := TextWr.New();

  PROCEDURE CopyLine() RAISES {Rd.EndOfFile} =
  (*
    Copy one comment line from "rd" to "wr", removing the CommentPrefix
    but leaving the final (mandatory) newline.
    Raises Rd.EndOfFile if "rd" is exhausted or the next char is 
    not CommentPrefix. *)
    <* FATAL Rd.Failure, Wr.Failure, Thread.Alerted *>
    <* FATAL MissingFinalNewLine *>
    VAR c: CHAR;
    BEGIN
      c := Rd.GetChar(rd); (* If EOF here, propagate to caller *)
      IF c # CommentPrefix THEN Rd.UnGetChar(rd); RAISE Rd.EndOfFile END;
      REPEAT
        TRY c := Rd.GetChar(rd) EXCEPT Rd.EndOfFile => RAISE MissingFinalNewLine END;
        Wr.PutChar(wr, c)
      UNTIL c = '\n'
    END CopyLine;

  VAR twr: TextWr.T;
      txt: TEXT;
  BEGIN
    TRY LOOP CopyLine() END EXCEPT Rd.EndOfFile => (* Ok *) END;
    twr := wr;
    txt := TextWr.ToText(twr);
    RETURN txt
  END ReadComments;  
 
<*UNUSED*>
PROCEDURE LR3Print(wr: Wr.T; w: LR3.T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr,
      "( " & Fmt.LongReal(w[0], prec := 7) & " "
           & Fmt.LongReal(w[1], prec := 7) & " "
           & Fmt.LongReal(w[2], prec := 7) & " )"
    );
  END LR3Print;

BEGIN
END Triang.