MODULE SMap;
(* 
  A map will be represented by an object with two fields: "f" and "c"
  where "f" is a face describing the map and "c" is a corner of the
  map whose left face is "f".  In fact, "f" is not any face of the
  map. To provide auxiliary information for point location, we will
  assume the existence of a reference point which will be defined as
  the North Pole (the point whose homogeneous coordinate is
  "[1,0,0,1]" and whose Plucker coefficients is
  "<0,0,1,0,0,0>"). Thus, a map is described by its face containing
  the reference point. In the particular case, when the reference
  point is on a border of the map, that is, when it is on an edge or
  it is a vertex, we will store in the field "c", the corner representing
  that edge or that vertex and "f" is the left face of that
  corner. It's worth mentioning that if the reference point lies on a
  face (i.e, it's not on an edge neither on a vertex) the field
  "c" will be set to "NIL".

  In the procedures involving intersection between an edge "e" and map
  (or some borders of it), we use a parameter "IntersectAndBreak" or
  "IntersectButNoBreak" to speficy whether we must only compute the
  intersection points or must break the edges also. In both cases, we
  store the intersection points in a list, where each element consists
  of the "coefficients" of the intersection point (we need store these
  coefficients because sometimes we don't break the existent edges)
  and either the "new" corner generated in breaking an edge of the map
  (when using "IntersectAndBreak") or an "old" corner of the map whose
  intersection with the edge has generated that point (when using
  "IntersectButNoBreak"). It is important to say that the intersection
  points are always computed as "(new edge) \meet (existent edge)".
  
  Besides, when the intersection point is generated by an isolated
  vertex, we store in the field "ipsubst" the corner that substitutes
  the isolated vertex after connecting it (via "Splice"). On the other
  hand, if the intersection point wasn't generated by an isolated
  vertex then "ivsubst" will be "NIL".

  Created in Sep. 16, 1997 by Marcus Vinicius A. Andrade
  Last edited in May 15, 1998 by Marcus Vinicius A. Andrade
*)
  
FROM SMGeo IMPORT Point, Circle, Arc, Sign;
FROM SMC IMPORT Face, Border, Edge, Vertex, Corner; 
FROM PSPlot IMPORT Color;

IMPORT SMGeo, SMC, SMBorder, SMPoint, SMCircle, SMDraw, Fmt, FileFmt, NGet, NPut, FGet, FPut, FileRd, FileWr, Lex, Rd, Wr, Text, Stdio, Extras;

CONST 
  (* ===== Constants for file format ===== *)
  SMVersion = "97-10-20";
  
  (* ===== Constant defining the kind of perturbation ==== *)
  CPerturb = +1;
  NotPerturb = 0;

  (* ===== Keys to control the walk on map, faces and borders === *)
  NotGetFaces = FALSE;
  IntersectAndBreak = TRUE;
  IntersectButNoBreak = FALSE;

REVEAL 
  (* ====== list of intersection points ====== *) 
  InterPoint = BRANDED REF RECORD 
    p      : Point;       (* coordinates of the intersection point  *)
    which  : Corner;      (* the (onext) corner where the new edge leaves "p" *)
    coming : Corner;      (* the (oprev) corner where the new edge comes to "p"*)
    next   : InterPoint;  (* next intersection point                *)
  END;  
   
  ListOfInterPoints = BRANDED REF RECORD  
    n : CARDINAL;
    ult: InterPoint;
    l : InterPoint;
  END;

  (* ========== Definition of Map ================ *)  
  T = PublicMap BRANDED OBJECT 
    ValidNumber : BOOLEAN:=TRUE; (* TRUE => map's numeration is valid          *)
    NV, NE      : CARDINAL;      (* number of Vertices and Edges               *)
  OVERRIDES                
    init := CreateMap;
    makeVertex := MakeVertex;
    delIsolatedVertex := DelIsolatedVertex;
    makeEdge := MakeEdge;
    delEdge := DelEdge;
    whereTo := WhereTo;
    whichEdge := WhichEdge;
    pointLocation := PointLocation;
    intersectEdgeAndElement := IntersectEdgeAndElement;
    numerate := NumerateMap;
    read := ReadMap;
    write := WriteMap;
    draw := DrawMap;
  END; 

VAR RefPoint : InterPoint;

(* ====================================================================== *)
(*   procedures to manipulate the geometry of vertices, edges and arcs    *)
(* ====================================================================== *)

PROCEDURE SetVertexGeometry(v: Corner; p: Point) = 
  BEGIN
    v.org().setGeom(p)
  END SetVertexGeometry;

PROCEDURE GetVertexGeometry(v: Corner) : Point = 
  BEGIN
    RETURN v.org().getGeom()
  END GetVertexGeometry;

PROCEDURE SetEdgeGeometry(e: Corner; sc: Circle; p: Point:=NIL; q: Point:=NIL) = 
  BEGIN
    IF p#NIL AND NOT e.isOval() THEN  
      SetVertexGeometry(e,p)
    END;

    e.halfEdge().initGeom(sc,+1);          (* --- sets the geometry of "e" --- *) 

    WITH es = e.sym() DO 
      (* --- "e.edg" and "es.edg" use a same edge's feature record; 
              but "e" and "es" have opposite direction given by "d" --- *)  
      es.setHalfEdge(e.halfEdge().edge(),-1);   

      IF q # NIL AND NOT es.isOval() THEN 
        IF SMGeo.AreSamePoints(p,q) THEN
          (* --- since "p = q" then "e.org" and "es.org" use a  
                 same vertex's feature record --- *) 
          es.setVertex(e.org())      
        ELSE 
          SetVertexGeometry(es,q);   
        END
      END;
    END
  END SetEdgeGeometry;

PROCEDURE GetEdgeGeometry(e: Corner) : Circle = 
  BEGIN
    RETURN e.halfEdge().getGeom()
  END GetEdgeGeometry;

PROCEDURE GetArcGeometry(e: Corner) : Arc = 
  BEGIN
    RETURN e.getGeoOfArc(e.sym())
  END GetArcGeometry;

PROCEDURE CreateCornerAndFeatures(sc: Circle; p,q: Point) : Corner = 
  BEGIN
    IF p=NIL AND q=NIL THEN (* --- creates a oval --- *)
      WITH e=NEW(Corner).createOval() DO 
        SetEdgeGeometry(e,sc);
        RETURN e
      END
    ELSIF sc=NIL AND p#NIL THEN (* --- creates an isolated vertex --- *)
      WITH v=NEW(Corner).createIsolatedVertex() DO 
        SetVertexGeometry(v,p);
        RETURN v
      END
    ELSE (* --- creates an arc --- *)
      WITH e = NEW(Corner).createEdge() DO

        (* --- if "p = q", sets the vertex of "e.sym()" to the same as "e" --- *)
        IF SMGeo.AreSamePoints(p,q) THEN 
          e.sym().setVertex(e.org())
        END;

        SetEdgeGeometry(e,sc,p,q);

        RETURN e
      END
    END
  END CreateCornerAndFeatures;

(* ================================================= *)
(*   Manipulation of the List of Intersection Points *)
(* ================================================= *)
PROCEDURE InitializeLIP() : ListOfInterPoints = 
  BEGIN
    WITH lip = NEW(ListOfInterPoints) DO
      lip^.n:=0;
      lip^.l:=NIL;
      lip^.ult:=NIL;
      RETURN lip
    END
  END InitializeLIP;

PROCEDURE InsNode(lip: ListOfInterPoints; p: Point; wh,come: Corner) =
  BEGIN
    WITH node = NEW(InterPoint) DO
      node^.p := p;
      node^.which := wh; 
      node^.coming := come;
      node^.next := NIL;
      IF lip^.n > 0 THEN (* list has at least one element *)
        lip^.ult^.next := node;
      ELSE  (* node is the first element *)
        lip^.l := node
      END;  
      lip^.ult:=node;
      INC(lip^.n)
    END;
  END InsNode;

<* UNUSED *> PROCEDURE MergeList(l1,l2: ListOfInterPoints) : ListOfInterPoints =
  BEGIN
    IF l2^.l = NIL THEN 
      RETURN l1
    ELSIF l1^.l = NIL THEN 
      RETURN l2
    ELSE
      l1^.ult^.next := l2^.l; 
      l1^.ult:=l2^.ult;
      l1^.n:=l1^.n+l2^.n;
      RETURN l1
    END
  END MergeList;  

(* ================================================= *)
(*  Manipulation of the Stack of Vertices and Edges  *)
(* ================================================= *)
TYPE 
  VStack = REF ARRAY OF Vertex;
  EStack = REF ARRAY OF Edge;
  BStack = REF ARRAY OF Border;
  FStack = REF ARRAY OF Face;
  
CONST 
  StackSize=127;

PROCEDURE IsOutOfStackRange(stack: REFANY; i: CARDINAL) : BOOLEAN = 
  BEGIN
    TYPECASE stack OF 
    | VStack(vStack) => RETURN i > NUMBER(vStack^)-1;
    | EStack(eStack) => RETURN i > NUMBER(eStack^)-1;
    | BStack(bStack) => RETURN i > NUMBER(bStack^)-1;
    | FStack(fStack) => RETURN i > NUMBER(fStack^)-1;
    ELSE 
      RETURN TRUE;
    END
  END IsOutOfStackRange;

PROCEDURE DoubleStackSize(stack: REFANY) : REFANY = 
  BEGIN
    TYPECASE stack OF 
    | VStack(vStack) => WITH
                          ss = NUMBER(vStack^),
                          ssn = 2*ss,
                          newstack = NEW(VStack,ssn)
                        DO
                          SUBARRAY(newstack^,0,ss) := vStack^;
                          RETURN newstack
                        END;
    | EStack(eStack) => WITH
                          ss = NUMBER(eStack^),
                          ssn = 2*ss,
                          newstack = NEW(EStack,ssn)
                        DO
                          SUBARRAY(newstack^,0,ss) := eStack^;
                          RETURN newstack;
                        END;

    | BStack(bStack) => WITH
                          ss = NUMBER(bStack^),
                          ssn = 2*ss,
                          newstack = NEW(BStack,ssn)
                        DO
                          SUBARRAY(newstack^,0,ss) := bStack^;
                          RETURN newstack
                        END;
    | FStack(fStack) => WITH
                          ss = NUMBER(fStack^),
                          ssn = 2*ss,
                          newstack = NEW(FStack,ssn)
                        DO
                          SUBARRAY(newstack^,0,ss) := fStack^;
                          RETURN newstack
                        END;
    ELSE
       RETURN stack;
    END;
  END DoubleStackSize;

(* ================================================= *)
(*          General Interna Procedures               *)
(* ================================================= *)
PROCEDURE SetRefFace(m: T; f: Face) = 
  VAR c: Corner;
  BEGIN 
    m.setRefFace(f);
    c := m.refCorner();
    IF c#NIL AND SMGeo.AreSamePoints(RefPoint^.p,GetVertexGeometry(c)) THEN 
      WHILE c.left() # f DO 
        c:=c.oNext();
      END;
      m.setRefCorner(c)
    END    
  END SetRefFace;

PROCEDURE NextPointOnCircle(ip: InterPoint; sc: Circle; lip: ListOfInterPoints) : InterPoint RAISES ANY =

  (* --------------------------------------------------------------------- *)
  (* Given a list of points on a circle "sc" and a point "p", also on that *)
  (* edge, returns the first point after "p" in that list considering      *)
  (* the orientation of "e"                                                *)
  (* --------------------------------------------------------------------- *)
  VAR f,r,br,bf  : InterPoint;  (* bf=before "fap"; br=before "r" *)
  BEGIN
    
    IF lip^.l = NIL THEN   (* if list is empty *)
      RETURN NIL
    ELSE
      bf:=NIL;
      f:=lip^.l;         (* "f" starts as the first point in the list *)
      br:=f;
      r:=f^.next;
      
      (* --- checks if the first point in "lip" is equal to "ip" --- *)
      IF NOT SMGeo.AreSamePoints(ip^.p,f^.p) THEN  
        WHILE r # NIL DO 
          WITH s = SMGeo.CircOrd3PointsOnCircle(ip^.p,r^.p,f^.p,sc) DO
            IF s > 0 THEN   (* "r" is between "p" and "fap" *)
              bf:=br;
              f:=r;
            ELSIF s=0 THEN  (* "r" is equal to "p" or to "fap" *)
              IF SMGeo.AreSamePoints(ip^.p,r^.p) THEN  
                (* "r" is equal to the initial point then we can stop the search *)
                bf:=br;
                f:=r;
                EXIT 
              ELSIF SMGeo.AreSamePoints(ip^.p,f^.p) THEN 
                EXIT 
              END
            END
          END;
          br:=r;
          r:=r^.next;
        END
      END;
      (* ====== remove "f" from the list "lip" ====== *)
      IF bf = NIL THEN 
        (* "f" is the first point in the list *)
        lip^.l:=lip^.l^.next
      ELSE
        bf^.next:=f^.next  
      END;
      f^.next:=NIL;
      lip^.n:=lip^.n-1;
      RETURN f
    END
  END NextPointOnCircle;

PROCEDURE SortPointsOnEdge(lip: ListOfInterPoints; e: Corner) : ListOfInterPoints RAISES ANY = 
  VAR slip : ListOfInterPoints;  (* sorted list *)
      ip,ult  : InterPoint;
      sc : Circle;
  BEGIN
    sc:=GetEdgeGeometry(e);
    IF lip^.n = 1 THEN 
      WITH auxslip = ClearMultipleOccurrences(lip,sc) DO 
        RETURN auxslip
      END 
    ELSE 
      slip := InitializeLIP();
      (* ---------------------------------------------------------- *)
      (* Creates an InterPoint record for an initial point on edge; *)
      (* if "e" is an oval, it's taken the first point in "lip"     *)
      (* ---------------------------------------------------------- *)
      ip := NEW(InterPoint);
      IF e.isOval() AND lip^.l # NIL THEN  
        ip^.p := lip^.l^.p
      ELSE
        ip^.p := GetVertexGeometry(e)
      END;
      ip^.which := e;
      ip^.coming := NIL;
      ip^.next:=NIL;

      (* ====== Inserts the first point in the list "slip" ====== *)
      WITH ipaux = NextPointOnCircle(ip,sc,lip) DO 
        InsNode(slip,GetVertexGeometry(ipaux^.which),ipaux^.which,ipaux^.coming)
      END;

      (* ====== Sorts the points in "lip" and puts them in "slip" ====== *)
      ult:= slip^.l;
      ip := lip^.l;
      WHILE ip # NIL DO 
        ip := NextPointOnCircle(ult,sc,lip);
        IF ip # NIL THEN 
          InsNode(slip,GetVertexGeometry(ip^.which),ip^.which,ip^.coming);
          ult:=ip
        END
      END;

      WITH auxslip = ClearMultipleOccurrences(slip,sc) DO 
        RETURN auxslip
      END;
    END
  END SortPointsOnEdge;

PROCEDURE UpdateLeavingAndComing(ip: InterPoint; sc: Circle; VAR leave, come: Corner) RAISES ANY = 
  BEGIN
    IF NumCornersAtSameVertex(ip^.which) > 1 THEN 
      EVAL InternalWhereTo(ip^.which,sc,leave);
      EVAL InternalWhereTo(ip^.which,SMGeo.InvCircle(sc),come)
    ELSE 
      leave:=ip^.which;
      come:=NIL;
    END;
  END UpdateLeavingAndComing;

PROCEDURE ClearMultipleOccurrences(slip: ListOfInterPoints; sc: Circle) : ListOfInterPoints RAISES ANY =
  (* In the list of intersection points, each point must occur exactly
     one time; so, after sorting the points on the circle "sc" we have
     to check if any point appear more than one time (because of many
     corner with a same origin that was intersected by the circle
     "sc") and, if so, we have to determine which corner on
     the vertex must be preserved; the other corners must be
     eliminated. *)
  VAR auxlip : ListOfInterPoints;
      fp,p   : InterPoint;
      leave,come : Corner;
      numoccur: CARDINAL;
  BEGIN
    IF slip^.n < 1 THEN  
      RETURN slip 
    ELSIF slip^.n = 1 THEN
      fp:=slip^.l;
      UpdateLeavingAndComing(fp,sc,leave,come);
      fp^.which:=leave;
      fp^.coming:=come;
      RETURN slip
    ELSE 
      auxlip:=InitializeLIP();
      fp:=slip^.l;
      p:=slip^.l^.next;
      WHILE fp#NIL DO 
        numoccur := 1;
        WHILE p#NIL AND SMGeo.AreSamePoints(GetVertexGeometry(fp^.which),
                                            GetVertexGeometry(p^.which)) DO 
          p:=p^.next;
          INC(numoccur);
        END;
          
        UpdateLeavingAndComing(fp,sc,leave,come);
        
        InsNode(auxlip,fp^.p,leave,come);
        fp:=p;
        IF p#NIL THEN    (* not finish yet *)
          p:=p^.next
        END
      END;
      RETURN auxlip
    END;
  END ClearMultipleOccurrences;

PROCEDURE GeneratePointOnBorder(b: Border) : Point = 
  (* ------------------------------------------------------------- *)
  (* returns a point on the border "b" and the corner where it is  *)
  (* ------------------------------------------------------------- *)
  VAR p: Point;
  BEGIN
    WITH c = b.borderDesc() DO
      IF NOT c.isOval() THEN 
        p:=GetVertexGeometry(c)
      ELSE
        p:=SMGeo.GeneratePointOnCircle(GetEdgeGeometry(c))
      END;
    END;
    RETURN p
  END GeneratePointOnBorder;

<* UNUSED *> PROCEDURE GetCrossingCornerAtVertex(crn: Corner; c: Circle) : Corner = 
  (* returns the next corner on the border "crn.borderOf()" or 
     "crn.sym.borderOf()" such that its supporting circle is different of "c" *)
  VAR d: Corner;
  BEGIN
    IF crn.isIsolatedVertex() OR 
          SMGeo.AreSameCircles(c,GetEdgeGeometry(crn)) = 0 THEN 
      RETURN crn
    END;
    
    WITH 
      bc = crn.borderOf(),
      bcs = crn.sym().borderOf()
    DO 
      d:=crn;
      REPEAT 
        d:=d.oPrev();
      UNTIL d=crn OR 
           (d.borderOf() = bc OR d.borderOf() = bcs 
              AND SMGeo.AreSameCircles(c,GetEdgeGeometry(d)) = 0);
    END;
    IF d=crn THEN 
      RETURN NIL 
    ELSE 
      RETURN d
    END
  END GetCrossingCornerAtVertex;

PROCEDURE SideOfPointRelativeToFirstCrossingCorner(p: Point; circ: Circle; lip: ListOfInterPoints; OnlyOneBorder: BOOLEAN:=FALSE; bc: Border:=NIL) : Sign RAISES ANY = 
  VAR pp,cp : InterPoint;
      q  := NEW(Point);
      d1 : Corner;
      aretang : BOOLEAN;
      s     : Sign;
  BEGIN 
    pp:=NEW(InterPoint);
    pp^.p:=p;
    pp^.which:=NIL;
    pp^.coming := NIL;
    pp^.next:=NIL;
   
    REPEAT 
      cp:=NextPointOnCircle(pp,circ,lip);

      IF (cp=NIL OR cp^.which.isIsolatedVertex()) THEN 
        (* since the edge didn't cross any edge of the border or the 
           first crossing was at an isolated vertex then the point is 
           on the left side of the border *)
        s:=+1;
      ELSIF SMGeo.IsPointInArc(p,GetArcGeometry(cp^.which),s) THEN         
         (* if point "p" is on that edge, breaks repeat-until and returns 0 *)
          RETURN 0
      ELSE 
        (* if the intersection point is a vertex  *) 
        IF NOT cp^.which.isOval() 
             AND SMGeo.AreSamePoints(cp^.p,GetVertexGeometry(cp^.which)) THEN 

          IF NumCornersAtSameVertex(cp^.which) > 1 THEN 
            WITH 
              cinv = SMGeo.InvCircle(circ)
            DO 
              EVAL InternalWhereTo(cp^.which,cinv,d1,OnlyOneBorder,bc);

              IF d1 = NIL THEN 
                (* there isn't a crossing corner at that vertex *)
                s:=1
              ELSE 
                (* since "WhereTo" returns the previous corner relative 
                   to "cinv" then "cinv" is entering on the positive 
                   side of the border "d1.border"; so, we only need to check 
                   whether "d1" is on the border "bc" or not. *)
                IF d1.borderOf() = bc THEN 
                  s := +1
                ELSE 
                  s := -1
                END 
              END
            END
          ELSE (* there is only one corner at that vertex *)
            s:=1
          END
        ELSE   (* the intersection point is not a vertex *)
          EVAL SMCircle.Intersection(circ^,GetEdgeGeometry(cp^.which)^,q^,aretang);

          IF SMGeo.AreSamePoints(cp^.p,q) THEN 
            (* edge whose supporting circle is "circ" is entering in the border *)
            s:=+1
          ELSE
            s:=-1
          END
        END
      END
    UNTIL s # 0;
    
    RETURN s
  END SideOfPointRelativeToFirstCrossingCorner;

PROCEDURE StraightLocalization(m: T; p,pb: Point; b: Border; c: Circle) : Sign RAISES ANY =
  (* since there is a rational circle "c" connecting the points "p" and "q" 
     and since "pb" is a point on the border "b", we can determine the 
     position of the point "p" relative to "b" intersecting the arc "(p,pb,c)" 
     with the border "b" and checking the position of $p$ relative to the 
     "first" intersection point. We CANNOT use the parity of the number of 
     intersection because in a vertex we can have an indefinite number of edges; 
     not only two as in a simple polygon.*)

  VAR lip : ListOfInterPoints;
  BEGIN 
    WITH 
      e = CreateCornerAndFeatures(c,p,pb)
    DO 
      lip:=m.intersectEdgeAndElement(e,IntersectButNoBreak,b);

      IF lip^.l = NIL THEN 
        (* the edge "e" connecting the searched point "p" and the generated 
           point on the border didn't cross any edge of the border; maybe, 
           because "e" is tangent to the borders's edge at the generated 
           point; so, since we take "e \meet border's edge" then "p" is 
           on the positive side of the border *)
        RETURN 1
      ELSE
        RETURN SideOfPointRelativeToFirstCrossingCorner(p,c,lip,TRUE,b);  
      END
    END
  END StraightLocalization;

PROCEDURE IndirectLocalization(m : T; p,pb: Point; b: Border) : Sign RAISES ANY = 
  (* since there isn't a rational circle connecting the point "p" and 
     a point on the border "b" then we will determine the position 
     of the point "p" relative to the border "b" indirectly; firstly, 
     we will compute the position of the "RefPoint" relative to the 
     border "b"; secondly, we will determine the intersection points 
     of the arc "(p,RefPoint)" and use these intersection points to  
     determine if the relative position of the point "p" is the 
     same as "RefPoint" *)
  VAR sref   : Sign;
      c      : Circle;
      lip    : ListOfInterPoints;
  BEGIN
  
    (* locates RefPoint relative to the border "b" *)
    EVAL SMGeo.CircleByTwoPoints(RefPoint^.p,pb,c);
    sref:=StraightLocalization(m,RefPoint^.p,pb,b,c);
    
    IF SMGeo.AreSamePoints(RefPoint^.p,p) THEN 
      RETURN sref
    END;

    (* relative position of RefPoint and "p" *)
    EVAL SMGeo.CircleByTwoPoints(RefPoint^.p,p,c);
    WITH 
      e = CreateCornerAndFeatures(c,RefPoint^.p,p)
    DO 
      lip:=m.intersectEdgeAndElement(e,IntersectButNoBreak,b);

      IF lip^.l = NIL THEN 
        (* since there isn't any intersection point then the position of 
           "p" relative to the border "b" is the same as "RefPoint" *)
        RETURN sref
      ELSE
        RETURN SideOfPointRelativeToFirstCrossingCorner(p,c,lip,TRUE,b);  

      END
    END
  END IndirectLocalization;

PROCEDURE LocatePointRelativeToBorder(m : T; p: Point; b: Border) : Sign RAISES ANY =
  (* ------------------------------------------------------------------ *) 
  (* Returns +1 if the point "p" is on the left side of the border "b"; *)
  (* that is, if "p" is "inside" the border "b". If "p" is on the       *)
  (* border "b" then returns 0 and if "p" is on the right side of "b"   *)
  (* then returns -1.                                                   *)
  (* ------------------------------------------------------------------ *) 
  VAR pb : Point;
      c : Circle;
  BEGIN
    WITH crn = b.borderDesc() DO 
      IF crn.isOval() THEN 
        RETURN SMGeo.SidePointCircle(p,GetEdgeGeometry(crn))
      END
    END;

    pb:=GeneratePointOnBorder(b);
    IF SMGeo.AreSamePoints(p,pb) THEN 
      RETURN 0
    ELSE 
      IF SMGeo.CircleByTwoPoints(p,pb,c) THEN 
        (* there is a rational circle connecting the points "p" and "pb" *)
        RETURN StraightLocalization(m,p,pb,b,c)
      ELSE
        RETURN IndirectLocalization(m,p,pb,b)
      END;
    END
  END LocatePointRelativeToBorder;


PROCEDURE MergeBorders(m: T; f1,f2 : Face) : Face =
   (* ---------------------------------------------------------- *)
   (* Returns a "new" face with all the borders of "f1" and "f2" *)
   (* ---------------------------------------------------------- *)
   VAR sf,nf  : Face;  
      sb     : SMBorder.T;
      b      : Border;
  BEGIN
    (* ------------------------------------------------------------- *)
    (* choose the "smallest" face, i.e. the face with fewest borders *) 
    (* ------------------------------------------------------------- *)
    IF f1.genus() <= f2.genus() THEN 
      sf := f1;
      nf := f2;
    ELSE
      sf := f2;
      nf := f1;
    END;
    (* ------------------------------------------------------------- *)
    (* transfer the borders of the "smallest" face to the "new face  *)
    (* ------------------------------------------------------------- *)
    WITH it = sf.bordersOf().iterate() DO
      WHILE it.next(sb) DO 
        b:=NARROW(sb,Border);
        nf.transferBorder(b)
      END;
    END;
    (* ------------------------------------------------------------- *)
    (* if the "smallest" face is the face containing the reference   *)
    (* point then updates the record describing the map              *)
    (* ------------------------------------------------------------- *)
    IF m.refFace() = sf THEN  
      SetRefFace(m,nf)
    END;
    RETURN nf
  END MergeBorders;

PROCEDURE DistributeBorders(m: T; f,fc,fd : Face) RAISES ANY = 
  (* 
    Given the set of borders "f", distributes its borders among the
    faces "fc" and "fd". At the begining, each face "fc" and "fd" have
    only one border, respectively "bc" and "bd", and any border of "f"
    doesn't intersect the initial borders of "fc" and "fd". So, we
    need only to determine in which face each border must be
    included. To realize this operation, we generate a point "p" on an
    initial border (by definition, the border "bc") and for each
    border "b" of "f" we generate a point "pb" on "b"; create an edge
    connecting "p" to "pb"; intersect this edge with the borders
    "bc" and "bd"; get the next point on that edge (sorted from the
    point "pb") and finally, the corner at that point gives us in
    which face the border "b" must be included *)
    
  VAR bc,b : Border;        
      sb : SMBorder.T; 
      s: Sign;
      
  BEGIN
    (* ----------------------------------------------------- *)
    (* Distributes the borders of "f" between "fc" and "fd"  *)
    (* ----------------------------------------------------- *)
    IF fc.genus() > 0  THEN (* "fc" has (at least) one border *)
      IF fc.bordersOf().iterate().next(sb) THEN 
        bc:=NARROW(sb,Border);
      END
    ELSE 
      RAISE ValidationError(" **** ERROR in DistributeBorders --- one face is empty");
    END;

    WITH itf = f.bordersOf().iterate() DO
      WHILE itf.next(sb) DO 
        b:=NARROW(sb,Border);
        WITH 
          p = GeneratePointOnBorder(b)
        DO 
          s := LocatePointRelativeToBorder(m,p,bc);
          IF s > 0 THEN 
            fc.transferBorder(b)
          ELSIF s < 0 THEN 
            fd.transferBorder(b)
          ELSE 
            RAISE ValidationError(" **** ERROR in DistributeBorders --- intersection between distinc borders")
          END
        END  
      END
    END;

    (* ------------------------------------------------------------ *)
    (* If "f" (face to be distributed) is the face containing the   *)
    (* reference point, then determine which face will receive it   *)
    (* ------------------------------------------------------------ *)
    IF f = m.refFace() THEN 
      s := LocatePointRelativeToBorder(m,RefPoint^.p,bc);
      IF s >= 0 THEN 
        SetRefFace(m,fc);
      ELSE 
        SetRefFace(m,fd);
      END
    END;
  END DistributeBorders;

PROCEDURE BreakEdge(m: T; p: Point; e: Corner) : Corner RAISES ANY = 
(* breaks the edge "e" inserting on it the point "p"; that is, 
   creates a pair of new corners "c" and "c.sym", such that 
   "c" is inserted just after "e" and "c.sym" after "e.sym". 
   Returns the corner "c". *)

  VAR wasoval : BOOLEAN;
      es,caux : Corner;
      s       : Sign;
  BEGIN
    wasoval := e.isOval();
    IF NOT wasoval THEN 
      IF SMGeo.AreSamePoints(p,GetVertexGeometry(e)) THEN 
        RETURN e
      END;
      IF SMGeo.AreSamePoints(p,GetVertexGeometry(e.sym())) THEN 
        RETURN e.sym()
      END
    END;
    
    es:=e.sym();
    WITH 
      c  = SMC.DropOnEdge(e),
      con = c.oNext()
    DO
      (* --- set the vertex's geometry of the new corners --- *)
      c.org().setGeom(p);
      con.org().setGeom(p);
      
      IF NOT wasoval THEN 
        (* --- updates the edge's geometry of the corners --- *)
        
        WITH 
          fwd = e.halfEdge().fwd()
        DO 
          e.sym().halfEdge().setGeom(e.halfEdge().edge(),-fwd);
          c.halfEdge().initGeom(e.halfEdge().getGeom(),+1);
          c.sym().halfEdge().setGeom(c.halfEdge().edge(),-1);
        END;

        (* ---------------------------------------------------- *)
        (* if "e" OR "es" is the Reference Corner, updates it;  *) 
        (* we have to check the case when the Reference point   *)
        (* is a vertex because in this case it is on (at least) *)
        (* two edges but not is the origin of these two corners *)
        (* ---------------------------------------------------- *)
        caux:=NIL;
        IF e = m.refCorner() THEN 
          caux:=c;
        ELSIF es = m.refCorner() THEN 
          caux:=con;
        END;
        IF caux # NIL THEN 
          IF SMGeo.IsPointInArc(RefPoint^.p,GetArcGeometry(caux),s) THEN 
            IF s >= 0 THEN (* "s>=0" implies that "p" is the origin or is on edge *)
              m.setRefCorner(caux)
            END
          END
        END;
      END;
      m.ValidNumber:=FALSE;
      RETURN c
    END;
  END BreakEdge;

PROCEDURE MakeOval(m: T; e: Corner; f: Face) RAISES ANY =
  VAR s : Sign;
  BEGIN
    WITH
      b = e.borderOf(),
      fc = b.faceOf(),
      fcs = e.sym().left()
    DO
      DistributeBorders(m,f,fc,fcs);
    END;
    IF m.refCorner() = NIL 
       AND SMGeo.IsPointInArc(RefPoint^.p,GetArcGeometry(e),s) THEN
      IF e.left() = m.refFace() THEN 
        m.setRefCorner(e)
      ELSIF e.sym().left() = m.refFace() THEN 
        m.setRefCorner(e.sym())
      ELSE 
        <* NOWARN *> Wr.PutText(Stdio.stdout,"**** ERRO em MakeOval --- inconsistencia na atualzacao de RefCorner n");
        <* NOWARN *> Wr.Flush(Stdio.stdout);
        WITH c=0 DO EVAL 1 DIV c END
      END
    END;
    m.ValidNumber:=FALSE;
  END MakeOval;

<* UNUSED *> PROCEDURE DelOval(m: T; e: Corner) : Face RAISES {ValidationError} = 
  BEGIN
    IF NOT e.isOval() THEN 
      RAISE ValidationError("In DelOval: the corner MUST be a Oval")
    ELSE
      WITH 
        be = e.borderOf(),
        fe = be.faceOf(),
        es = e.sym(),
        bes = es.borderOf(),
        fes = bes.faceOf()
      DO
        be.removeCorner(e);
        bes.removeCorner(es);
        m.ValidNumber:=FALSE;
        RETURN MergeBorders(m,fe,fes)
      END;
    END
  END DelOval;

PROCEDURE CheckConnectionAlreadyExists(c : Corner; e: Corner) : BOOLEAN =
  (* checks if there is an edge with same geometry as "e" already 
     connecting the corner "c" *)
  VAR sc   : Circle;
      d    : Corner;
      exist: BOOLEAN;
  BEGIN 
    sc := GetEdgeGeometry(e);
    exist:=FALSE;
    IF NOT c.isIsolatedVertex() THEN 
      d:=c;
      REPEAT 
        WITH 
          s=SMGeo.AreSameCircles(GetEdgeGeometry(d),sc)
        DO 
          exist := (s=1)
        END;
        d:=d.oNext()
      UNTIL exist OR (c=d);
    END;
    RETURN exist
  END CheckConnectionAlreadyExists;

PROCEDURE Connect(m: T; v: Corner; e: Corner) RAISES ANY = 
  VAR f  : Face;
  BEGIN
    e.setVertex(v.org());
    IF NOT CheckConnectionAlreadyExists(v,e) THEN 
      WITH b = e.borderOf() DO
        f := v.left();
        IF f # b.faceOf() THEN 
          f.transferBorder(b)
        END;

        IF v.isIsolatedVertex() THEN
          (* to replace the isolated vertex with the edge just remove the vertex  *)
          m.delIsolatedVertex(v);
          IF SMGeo.AreSamePoints(RefPoint^.p,GetVertexGeometry(e)) THEN 
            m.setRefCorner(e)
          ELSIF SMGeo.AreSamePoints(RefPoint^.p,GetVertexGeometry(e.sym())) THEN 
            m.setRefCorner(e.sym())
          END;
        ELSE 
          v.splice(e); 
          WITH 
            nfe = e.left(),
            nfv = v.left()
          DO 
            (* ----------------------------------------------------------------- *)
            (* After "Splice", if "v" and "e" are incident to a same face then   *)
            (* the borders of the "old" face "f" must be transfered to the "new" *)
            (* "nfe"; otherwise, the borders of "f" must be distribuited between *)
            (* the two different "new" faces "nfv" and "nfe"                     *)
            (* ----------------------------------------------------------------- *)
            IF nfe = nfv THEN 
              nfe.mergeFaces(f);   
              IF m.refFace() = f THEN 
                SetRefFace(m,nfe)
              END;
            ELSE 
              DistributeBorders(m,f,nfe,nfv);
            END;
          END;
        END
      END;
    END
  END Connect;

(* ================================================= *)
(*          Implementation of the methods            *)     
(* ================================================= *)

PROCEDURE CreateMap(self: T): T =
  BEGIN
    EVAL NARROW(self, SMC.Map).init();
    self.ValidNumber:= TRUE;
    self.NV:=0;
    self.NE:=0;
    RETURN self
  END CreateMap;

PROCEDURE MakeVertex(self: T; p: Point) : Corner RAISES ANY = 
  VAR v : Corner;
  BEGIN
    WITH elem = self.pointLocation(p) DO 
      TYPECASE elem OF 
      | Face(f) => BEGIN
                     (* --- creates a new isolated vertex in a dummy face --- *)
                     v := CreateCornerAndFeatures(NIL,p,NIL);

                     (* --- transfers the vertex from its dummy face to "f" --- *)
                     f.transferBorder(v.borderOf());
                   END   
      | Corner(e) => v:=BreakEdge(self,p,e);
      
      ELSE
      END;
      self.ValidNumber := FALSE; 
    END;
    (* --- if the new vertex lies on the reference point AND 
           there isn't any reference corner yet; set it   --- *) 
    IF self.refCorner() = NIL AND SMGeo.AreSamePoints(RefPoint^.p,p) THEN 
      self.setRefCorner(v)
    END;
    RETURN v
  END MakeVertex;

PROCEDURE DelIsolatedVertex(self: T; v: Corner) RAISES {ValidationError} =
  BEGIN
    IF NOT v.isIsolatedVertex() THEN 
      RAISE ValidationError("In DelIsolatedVertex: the corner MUST be an isolated vertex")
    ELSE
      v.borderOf().removeCorner(v);
      (* --- If the isolated vertex "v" lies on the reference point, 
             sets the field "c" to NIL --- *)
      IF self.refCorner() # NIL AND 
            SMGeo.AreSamePoints(RefPoint^.p,GetVertexGeometry(v)) THEN 
        self.setRefCorner(NIL)
      END;
    END;
    self.ValidNumber:=FALSE;
  END DelIsolatedVertex;

PROCEDURE NumCornersAtSameVertex(c: Corner) : CARDINAL = 
  VAR num : CARDINAL;
      d : Corner;
  BEGIN
    num:=1;
    d:=c.oNext();
    WHILE c # d DO 
      d:=d.oNext();
      INC(num)
    END;
    RETURN num
  END NumCornersAtSameVertex;

PROCEDURE CheckConnectionAtVertexWithOnlyOneCorner(ip: InterPoint; leaving: BOOLEAN; e: Corner) : Corner = 
  (* if the vertex to be connect has only one corner (for example, an
     isolated vertex) and if it was already connect, get the corner
     that replaced it or that have connected it; on the other hand, if
     the vertex will be first time connected just keeps track which
     corner will replace or connect it. The parameter "leaving"
     indicates whether we are coming or leaving at that point. *)
  VAR crn : Corner;
  BEGIN
    IF ip^.which.isIsolatedVertex() THEN 
      (* if it is an isolated vextex, so at the second connection it
         wasn't an isolated anymore and we need to return the corner
         that connected the "old" isolated vertex *)
      IF ip^.coming = NIL THEN 
      (* the vertex was NOT connect yet; so stores the edge that will
         connect the vertex (or replace it, in case of isolated
         vextex); so, the vertex will NOT be an isolated vertex
         anymore. *)

        ip^.coming := e;
        RETURN ip^.which
      ELSE 
        RETURN ip^.coming
      END
    ELSIF NumCornersAtSameVertex(ip^.which) = 1  THEN 
      (* it is a vertex with only one edge and we will connect at
         this vertex; so, the next connection need to consider the "new" 
         edge that was connected the vertex --- we cannot join this case 
         with the case above because after the first connection we will not 
         have only one corner at the vertex but we need to connect the edge 
         used in the first connection *)
      crn:=ip^.which;
      IF ip^.coming = NIL THEN
        ip^.coming := e;  (* for use in case of "mustclose" *)
        ip^.which := e    (* for use in the second connection *)
      END;
      RETURN crn
    ELSIF leaving THEN 
      RETURN ip^.which
    ELSE 
      RETURN ip^.coming
    END
  END CheckConnectionAtVertexWithOnlyOneCorner;

PROCEDURE MakeEdge(self: T; c: Circle; p: Point:=NIL; q: Point:=NIL) : Corner RAISES ANY = 
  VAR ip1,ip2      : InterPoint; 
      lip,slip     : ListOfInterPoints;
      u,v,crn      : Corner;
      e,en         : Corner;
      f1           : Face;               (* origin's face *)
      mustclose    : BOOLEAN;
      s            : Sign;
      p1           : Point;              (* point to search *)
  BEGIN
    e := CreateCornerAndFeatures(c,p,q);
    IF p#NIL THEN     (* --- edge to be created is not a oval --- *)
      p1 := p;
      v:=self.makeVertex(p);
    
      IF q#NIL AND NOT SMGeo.AreSamePoints(p,q) THEN 
        u:=self.makeVertex(q);
      ELSE 
        u:=v;
      END;
      IF NumCornersAtSameVertex(v) > 1 THEN 
        EVAL InternalWhereTo(v,c,en);
        f1:=en.left()
      ELSE 
        f1:=v.left()
      END
    ELSE 
      p1 := SMGeo.GeneratePointOnCircle(c);

      TYPECASE self.pointLocation(p1) OF 
      | Corner(d) => BEGIN
                       IF NOT d.isOval() AND SMGeo.AreSamePoints(p1,GetVertexGeometry(d)) THEN 
                         EVAL InternalWhereTo(d,c,en);
                         f1:=en.left()
                       ELSE 
                         f1:=d.left()
                       END
                     END;
      | Face(f) => f1:=f
      ELSE 
      END;
    END;

    lip := self.intersectEdgeAndElement(e,IntersectAndBreak,self,f1);
    
    (* --- inserts the ending points of the edge "e" --- *)
    IF NOT e.isOval() THEN 
      IF NOT v.isIsolatedVertex() THEN 
        InsNode(lip,GetVertexGeometry(v),v,NIL)
      END;
      IF NOT u.isIsolatedVertex() THEN 
        InsNode(lip,GetVertexGeometry(u),u,NIL)
      END;
    END;
    
    IF (lip^.l = NIL) THEN 
      IF (p=NIL) AND (q=NIL) THEN  (* --- creates a oval --- *)
        MakeOval(self,e,f1)
      END;
    ELSE
      slip := SortPointsOnEdge(lip,e);
      (* since, in case of edges, we first insert its extremes, then
         after sorting, if we get only one intersection point then we
         are inserting a loop and the intersection point stored in
         "slip" is an isolated vertex. Thus, "mustclose" will be true. *)

      mustclose:=e.isOval() OR SMGeo.AreSamePoints(p,q);
      
      (* -------------------------------------------------------------------- *)
      (* since the point "ip^.p" can lie on a existing vertex, we need to do  *)
      (* a "WhereTo" on "v1" to decide which corner on "v1" will be connected *)
      (* In fact, since the intersection points are computed by               *)
      (* "new-edge \meet existent-edge", all "ip^.which" follows the rule:    *)
      (*                                                                      *)
      (*                   /|\ f3                                             *)
      (*                <---|---- ("ip2^.which)                               *)
      (*                    | f2                                              *)
      (*                <---|---- ("ip1^.which)                               *)
      (*                    |  f1                                             *)
      (*               ("new" edge)                                           *)
      (*                                                                      *)
      (* So, in each step in the edge creation we must execute a "WhereTo" on *)
      (* the "new" vertex in "ip1^.which" to find which corner must be used   *)
      (* in the connection with "ip2^.which" itself. To create the next       *)
      (* subedge, repeat the process with the pair ("ip2";"ip3") ...          *)
      (* -------------------------------------------------------------------- *)

      IF slip^.n > 1 THEN 
        ip1:=slip^.l;
        ip2:=ip1^.next;
        WHILE ip2 # NIL DO 
          en := NEW(Corner).createEdge();
          SetEdgeGeometry(en,c);

          crn:=CheckConnectionAtVertexWithOnlyOneCorner(ip1,TRUE,en);

          Connect(self,crn,en);

          crn:=CheckConnectionAtVertexWithOnlyOneCorner(ip2,FALSE,en.sym());

          Connect(self,crn,en.sym());

          IF self.refCorner() = NIL 
              AND SMGeo.IsPointInArc(RefPoint^.p,GetArcGeometry(en),s) THEN 
            IF en.left() = self.refFace() THEN 
              self.setRefCorner(en)
            ELSIF en.sym().left() = self.refFace() THEN 
              self.setRefCorner(en.sym())
            END   
          END;
          ip1:=ip2;
          ip2:=ip2^.next;
        END
      END;

      (* As we saw above, if "slip^.n=1" then "mustclose = TRUE" and 
         so, this case is treated bellow *)

      IF mustclose THEN
        en := NEW(Corner).createEdge();
        SetEdgeGeometry(en,c);
        ip1:=slip^.ult;
        ip2:=slip^.l;
        crn:=CheckConnectionAtVertexWithOnlyOneCorner(ip1,TRUE,en);

        (* the insertion of an oval or loop tangent to an edge, we need 
           to do a "WhereTo" at the vertex to be connected to determine 
           which corner will be connected *)
           
        IF SMGeo.AreSamePoints(GetVertexGeometry(crn),ip1^.p) 
            AND NumCornersAtSameVertex(crn) > 1 THEN 
          EVAL InternalWhereTo(crn,c,u);
          crn:=u
        END;

        Connect(self,crn,en);
        crn:=CheckConnectionAtVertexWithOnlyOneCorner(ip2,FALSE,en.sym());

        (* idem above *)
        IF SMGeo.AreSamePoints(GetVertexGeometry(crn),ip2^.p) 
            AND NumCornersAtSameVertex(crn) > 1 THEN 
          EVAL InternalWhereTo(crn,SMGeo.InvCircle(c),u);
          crn:=u
        END;

        Connect(self,crn,en.sym())
      END;
    END;
    self.ValidNumber:=FALSE;
    RETURN en;
  END MakeEdge;

PROCEDURE DelEdge(<* UNUSED *>self: T; <* UNUSED *>e: Corner) = 
  BEGIN
  END DelEdge;

PROCEDURE PointLocation(self: T; p: Point) : REFANY RAISES ANY =
  (* starting from the reference point, we walk around the borders of
     the reference face looking for intersecting points; we construct
     an oriented edge "e" from the reference point to the point "p"
     and since all the edges "c" of the borders of the reference face 
     have the reference face on its left side then we can check any
     crossing by computing "e \meet c" *)

  VAR sc    : Circle;
      ip,cp : InterPoint;
      lip   : ListOfInterPoints;
      startingface : Face;
      crn,d : Corner;
      s     : Sign;
  BEGIN
    IF SMGeo.AreSamePoints(RefPoint^.p,p)  THEN 
      IF self.refCorner() # NIL THEN (* --- that vertex already exists --- *)
        RETURN self.refCorner()
      ELSE 
        RETURN self.refFace()
      END
    END;
      
    EVAL SMGeo.CircleByTwoPoints(RefPoint^.p,p,sc);

    WITH e = CreateCornerAndFeatures(sc,RefPoint^.p,p) DO
      WITH rc = self.refCorner() DO 
        IF rc # NIL AND NOT rc.isOval() AND 
           SMGeo.AreSamePoints(RefPoint^.p,GetVertexGeometry(rc)) THEN 
           (* there is a corner whose vertex is the RefPoint *)
          EVAL InternalWhereTo(rc,sc,d);
          startingface := d.left()
        ELSE
          startingface := self.refFace()
        END
      END;
      
      lip:=self.intersectEdgeAndElement(e,IntersectButNoBreak,self,startingface);

      IF lip^.l = NIL THEN 
        (* -----------------------------------------------------------*)
        (* if the edge connecting "p" to RefPoint doesn't intersect   *)
        (* the border "b" then "p" and RefPoint are on a same face.   *)
        (* ---------------------------------------------------------- *)
        RETURN self.refFace();
      ELSE
        ip:=NEW(InterPoint);
        ip^.p:=p;
        ip^.which:=NIL;
        ip^.coming := NIL;
        ip^.next:=NIL;

        cp:=NextPointOnCircle(ip,SMGeo.InvCircle(sc),lip);

        crn := cp^.which;
        IF crn.isIsolatedVertex() THEN 
          IF SMGeo.AreSamePoints(GetVertexGeometry(crn),p) THEN 
            RETURN crn
          ELSE
            RETURN crn.left()
          END
        ELSE
          IF NOT crn.isOval() THEN 
            IF SMGeo.AreSamePoints(cp^.p,GetVertexGeometry(crn)) 
                  AND NumCornersAtSameVertex(crn) > 1 THEN 
               EVAL InternalWhereTo(crn,sc,d);
               crn:=d
            END;
            (* checks the very particular case when the point "p" lies on 
               the ending vertex of an existent dangling edge *)
            IF SMGeo.AreSamePoints(p,GetVertexGeometry(crn.sym())) THEN 
              crn:=crn.sym()
            END
          END;
          
          IF SMGeo.IsPointInArc(p,GetArcGeometry(crn),s) THEN 
            (* "p" is on the arc defined by the corner "crn" *)
            RETURN crn 
          ELSE 
            s := SMGeo.SidePointCircle(p,GetEdgeGeometry(crn));
            (* since "crn" is the closest corner to the point "p" so we need 
               to check the position of "p" relative to the border of "crn" *)
            WITH 
              sconfirm = LocatePointRelativeToBorder(self,p,crn.borderOf())
            DO 
              IF sconfirm = 1 THEN 
                (* "p" is on the positive side of the border *)
                RETURN crn.left()
              ELSE 
                (* "p" is on the negative side of the border *)
                RETURN crn.sym().left()
              END
            END
          END
        END
      END
    END
  END PointLocation;
 
PROCEDURE InternalWhereTo(v: Corner; c: Circle; VAR e: Corner; OnlyOneBorder: BOOLEAN:=FALSE; b: Border:=NIL): REFANY RAISES ANY = 
  (* in general case, computes the function "WhereTo" as specified in the interface; 
     but, if the parameter "OnlyOneBorder" is defined as value "TRUE" so consider only 
     the corners on the border "b" *)
  VAR a1,a2      : Arc;
      v1,v2      : Corner;
      aretang : BOOLEAN;
      s : Sign;
      q : Point;
  BEGIN
    IF SMGeo.SidePointCircle(GetVertexGeometry(v),c) # 0 THEN 
      RAISE ValidationError("In WhereTo -- Circle doesn't pass through point ")
    ELSIF v.isIsolatedVertex() THEN 
      e:=v;
      RETURN v.left()
    ELSIF NumCornersAtSameVertex(v) = 1 THEN 
      e:=v;
      RETURN e
    ELSE 
      WITH 
        p = GetVertexGeometry(v),
        sarc = SMGeo.MakeArc(p,p,c)  
      DO 
        v1:=v;
        v2:=v.oNext();
        
        IF OnlyOneBorder THEN 
          WHILE v2.borderOf() # b AND v2.sym().borderOf() # b AND v2 # v1 DO 
            v2:=v2.oNext()
          END;
          IF v2=v1 THEN 
            e:=v1;
            RETURN e
          END 
        END;

        LOOP 
          a1 := SMGeo.MakeArc(p,p,GetEdgeGeometry(v1));
          a2 := SMGeo.MakeArc(p,p,GetEdgeGeometry(v2));
          s := SMGeo.CircOrdOf3ArcsAroundVertex(a1,sarc,a2,p);
          IF s = 1 THEN 
            (* --- in counterclockwise direction, "sarc" is after "a1" 
                   and before "a2" --- *)
            e:=v1;
            RETURN v1.left()
          ELSIF s = 0 THEN 
            (* --- the supporting circles of two arcs are coincident --- *)
            q:=NEW(Point);
            EVAL SMCircle.Intersection(c^,SMGeo.ScircleOfArc(a1)^,q^,aretang);
            IF SMPoint.IsCpointNull(q^) THEN 
              e:=v1;
              RETURN v1
            ELSE 
              EVAL SMCircle.Intersection(c^,SMGeo.ScircleOfArc(a2)^,q^,aretang);
              IF SMPoint.IsCpointNull(q^) THEN 
                e:=v2;
                RETURN v2
              END
            END;
          END;        
          v1:=v2;
          v2:=v2.oNext();
          IF OnlyOneBorder THEN 
            WHILE v2.borderOf() # b AND v2.sym().borderOf() # b AND v2 # v1 DO 
              v2:=v2.oNext()
            END;
            IF v2=v1 THEN 
              RETURN v1
            END 
          END;
        END
      END
    END 
  END InternalWhereTo;

PROCEDURE WhereTo(<* UNUSED *>self: T; v: Corner; c: Circle; VAR e: Corner): REFANY RAISES ANY = 
  BEGIN
    RETURN InternalWhereTo(v,c,e)  
  END WhereTo;
  
PROCEDURE WhichEdge(self: T; v: Corner; c: Circle): Corner RAISES ANY =
  VAR ip  := NEW(InterPoint);
      lip : ListOfInterPoints;
      f   : Face;   (* starting face where the searching will start *)
      en  : Corner;
  BEGIN
    IF NOT v.isOval() THEN 
      WITH
        p = GetVertexGeometry(v),
        e = CreateCornerAndFeatures(c,p,p)
      DO
        WITH elem = self.pointLocation(p) DO 
          TYPECASE elem OF 
          | Corner(e) =>  BEGIN
                            IF SMGeo.AreSamePoints(p,GetVertexGeometry(e)) THEN 
                              EVAL InternalWhereTo(e,c,en);
                              f:=en.left()
                            ELSE 
                              f:=e.left()
                            END
                          END;
          | Face(faux) => f:=faux
          ELSE 
          END        
        END;
        lip := self.intersectEdgeAndElement(e,IntersectButNoBreak,self,f);
        IF lip^.l = NIL THEN  (* === Didn't meet any edge === *)
          RETURN v           
        ELSE
          ip^.p:=p;
          ip^.which:=v;
          ip^.coming := NIL;
          ip^.next:=NIL;
          ip := NextPointOnCircle(ip,c,lip);
          WHILE lip^.l # NIL AND ip^.which.isIsolatedVertex() DO 
            ip := NextPointOnCircle(ip,c,lip);
          END;
          IF lip^.l = NIL THEN (* --- didn't meet an edge --- *) 
            RETURN ip^.which;
          ELSE 
            RETURN NIL 
          END  
        END
      END
    ELSE  (* v is an oval *)
      RETURN v;
    END
  END WhichEdge;


PROCEDURE IntersectEdgeAndElement(self: T; e: Corner; AlsoBreak: BOOLEAN;
                       elem: REFANY; startingface:Face:=NIL) : ListOfInterPoints RAISES ANY =
  (* Computes the points of intersection between the edge "e" and an
     element of the map (i.e the whole map, a face or a border); the
     parameter "AlsoBreak" indicates whether the elements intersecting
     the edge "e" will be broken. If "elem" is the whole map then
     "startingface" speficies a face where to start the searching.  *)

  VAR lip : ListOfInterPoints;
      earc,earcinv : Arc;
      ecirc : Circle;
      aretangent : BOOLEAN; 
      
  PROCEDURE IntersectTwoCorners(c: Corner) RAISES ANY = 
    VAR sp : Sign;

    PROCEDURE IntersectAndInsert(narc: Arc) RAISES ANY = 
      VAR p  : Point;
          nc : Corner;
          d1,d2 : Corner;
          ok : BOOLEAN;
      BEGIN
        WITH 
          carc = GetArcGeometry(c)
        DO 
          IF AlsoBreak THEN 
            ok := SMGeo.IntersectionOfTwoArcs(narc,carc,NotPerturb,aretangent,p)
          ELSE 
            ok := SMGeo.IntersectionOfTwoArcs(narc,carc,CPerturb,aretangent,p)
          END;
          IF ok AND aretangent 
              AND (NOT c.isOval() AND SMGeo.AreSamePoints(p,GetVertexGeometry(c)))
              AND NumCornersAtSameVertex(c) > 1 THEN 
            EVAL InternalWhereTo(c,ecirc,d1);
            EVAL InternalWhereTo(c,SMGeo.InvCircle(ecirc),d2);
            IF d1 = d2 AND NOT AlsoBreak THEN 
              (* the tangent circle "ecirc" is on a same side, that is, 
                 it doesn't crosses between two corners; so, if "AlsoBreak" is 
                 FALSE (in case of "PointLocation"), we don't consider the point 
                 "q" as an intersection *)
              ok:=FALSE
            END 
          END;

          IF ok THEN 
            IF AlsoBreak AND 
               (c.isOval() OR NOT SMGeo.AreSamePoints(p,GetVertexGeometry(c))) THEN 
              nc := BreakEdge(self,p,c);
            ELSE
              nc := c
            END;
            InsNode(lip,p,nc,NIL);
          END
        END
      END IntersectAndInsert;

    BEGIN   (* ---- IntersectTwoCorners ---- *)
      IF c.isIsolatedVertex() THEN 
         (* WITH v = c.org(), vn = v.getNum() DO *)
          WITH p = GetVertexGeometry(c) DO 
            IF SMGeo.IsPointInArc(p,earc,sp) THEN 
              InsNode(lip,p,c,NIL);
            END;
          (* END  *)
        END
      ELSE 
        IntersectAndInsert(earc);
        IF NOT AlsoBreak THEN 
          IntersectAndInsert(earcinv);
        END;
      END
    END IntersectTwoCorners;
    
  PROCEDURE IntersectEdgeWithMap() RAISES ANY = 
    VAR fStack : FStack;
        ip: InterPoint;
        wc: Corner;
    BEGIN
      WITH cnt=self.getTCounter() DO 
        fStack:=NEW(FStack,cnt.NF+1);
      END;

      (* --- computes the intersection between "e" and the starting face "f" --- *)
      startingface.walkOn(IntersectTwoCorners,NotGetFaces); 
      fStack[startingface.getNum()]:=startingface;

      ip:=lip^.l;
      WHILE ip # NIL DO
        (* if the edge "e" passes throw a vertex and there are more
           than two corners at this vertex then we have to do a
           "whereTo" on that vertex to decide in which face to continue 
           the search for intersection points. *)
        
        IF NumCornersAtSameVertex(ip^.which) > 2 AND
             SMGeo.AreSamePoints(ip^.p,GetVertexGeometry(ip^.which)) THEN 
          EVAL InternalWhereTo(ip^.which,GetEdgeGeometry(e),wc)
        ELSE
          wc:=ip^.which.sym()
        END;
                
        WITH
          f=wc.left(),
          fn=f.getNum()
        DO 
          IF fStack[fn] = NIL THEN   (* face was not visited yet *)
            f.walkOn(IntersectTwoCorners,NotGetFaces); 
            fStack[fn]:=f
          END
        END;
        ip:=ip^.next
      END
    END IntersectEdgeWithMap;
    
  BEGIN  (* ---- IntersectEdgeAndElement ---- *)
    (* --- initializes the list of intersection points --- *)
    lip := InitializeLIP();

    earc := GetArcGeometry(e);
    ecirc := GetEdgeGeometry(e);

    (* for point location we have to compute the intersection points 
       with the opposite arc *)
    IF NOT AlsoBreak THEN   
      earcinv := SMGeo.InvArc(earc)
    END;

    IF NOT self.ValidNumber THEN 
      self.numerate()
    END;
    TYPECASE elem OF 
    | T(m) => IntersectEdgeWithMap();         
    | Face(f) => f.walkOn(IntersectTwoCorners, NotGetFaces); 
    | Border(b) => b.walkOn(IntersectTwoCorners, NotGetFaces); 
    ELSE
    END;

    RETURN lip
  END IntersectEdgeAndElement;

PROCEDURE NumerateMap(self: T) RAISES ANY =
  VAR cnum,vnum,enum,bnum,fnum : CARDINAL;
      cnt : SMC.TCounter;
      vStack : VStack;
      eStack : EStack;
      bStack : BStack;
      fStack : FStack;
  
  PROCEDURE VisitAndNumerate(ic: Corner) = 
    BEGIN
      (* --- numerates the corner --- *) 

      INC(cnum);
      ic.setNum(cnum);

      (* --- marks and numerates vertex, edge and face --- *)
      IF NOT ic.isOval() THEN 
        WITH 
          v = ic.org(),
          vn = v.getNum() 
        DO 
          WHILE IsOutOfStackRange(vStack,vn) DO 
            vStack := DoubleStackSize(vStack) 
          END;
          (* --- checks if the vertex's number is up-to-date --- *)
          IF vn = 0 OR vStack[vn] # v THEN
            INC(vnum);
            IF IsOutOfStackRange(vStack,vnum) THEN vStack:=DoubleStackSize(vStack) END;
            v.setNum(vnum);
            vStack[vnum]:=v
          END
        END
      END; 
      IF NOT ic.isIsolatedVertex() THEN 
        WITH 
          sc = ic.halfEdge().edge(),
          scn = sc.getNum()
        DO 
          (* --- checks if the edge's number is up-to-date --- *)
          WHILE IsOutOfStackRange(eStack,scn) DO 
            eStack:=DoubleStackSize(eStack) 
          END;
          IF scn = 0 OR eStack[scn] # sc THEN
            INC(enum);
            
            IF IsOutOfStackRange(eStack,enum) THEN 
              eStack:=DoubleStackSize(eStack);
            END;
            sc.setNum(enum);
            eStack[enum]:=sc
          END
        END;
      END; 
      WITH 
        b = ic.borderOf(),
        bn = b.getNum()
      DO 
        (* --- checks if the border's number is up-to-date --- *)
        WHILE IsOutOfStackRange(bStack,bn) DO 
          bStack:=DoubleStackSize(bStack) 
        END;
        IF bn = 0 OR bStack[bn] # b THEN
          INC(bnum);
          IF IsOutOfStackRange(bStack,bnum) THEN bStack:=DoubleStackSize(bStack) END;
          b.setNum(bnum);
          bStack[bnum]:=b
        END
      END;

      WITH 
        f = ic.left(),
        fn = f.getNum()
      DO 
        (* --- checks if the face's number is up-to-date --- *)
        WHILE IsOutOfStackRange(fStack,fn) DO 
          fStack:=DoubleStackSize(fStack) 
        END;
        IF fn = 0 OR fStack[fn] # f THEN
          INC(fnum);
          IF IsOutOfStackRange(fStack,fnum) THEN fStack:=DoubleStackSize(fStack) END;
          f.setNum(fnum);
          fStack[fnum]:=f
        END
      END;
    END VisitAndNumerate;

  BEGIN
    vStack := NEW(VStack, StackSize);
    eStack := NEW(EStack, StackSize);
    bStack := NEW(BStack, StackSize);
    fStack := NEW(FStack, StackSize);
    
    cnum:=0; vnum:=0; enum:=0;  bnum:=0;  fnum:=0;
    self.walkOn(VisitAndNumerate);
    cnt.NC:=cnum; cnt.NB:=bnum; cnt.NF:=fnum;
    self.setTCounter(cnt);
    self.NV:=vnum; self.NE:=enum; 
    
    self.ValidNumber:=TRUE;    
  END NumerateMap;

PROCEDURE ReadMap(self: T; name: TEXT) RAISES ANY =
  CONST dig = SET OF CHAR{'0'..'9'};
  VAR rd : FileRd.T;
      vStack : VStack;
      eStack : EStack;
      cnt    : SMC.TCounter;
      
  PROCEDURE ReadVertices() RAISES ANY = 
    BEGIN
      FileFmt.ReadHeader(rd,"Vertices Data","");
      EVAL FileFmt.ReadComment(rd,'#');
      FGet.Skip(rd);
      WHILE Rd.GetChar(rd) IN dig DO 
        Rd.UnGetChar(rd);
        WITH  vnum = FGet.Int(rd) DO 
          IF vStack[vnum] # NIL THEN
            vStack[vnum].setNum(vnum);
            vStack[vnum].setStyle(FGet.Int(rd));
            vStack[vnum].setGeom(SMGeo.ReadPoint(rd))
          END;
        END;

        (* --- ignores the equivalent Spoint --- *)
        FGet.Skip(rd);
        FGet.Match(rd,"<=>");
        EVAL SMPoint.ReadSpoint(rd);

        FGet.EOL(rd);
        FGet.Skip(rd)
      END;
      Rd.UnGetChar(rd);
      FileFmt.ReadFooter(rd,"Vertices Data")
    END ReadVertices;

  PROCEDURE ReadEdges() RAISES ANY = 
    BEGIN
      FileFmt.ReadHeader(rd,"Edges Data","");
      EVAL FileFmt.ReadComment(rd,'#');
      Lex.Skip(rd);
      WHILE Rd.GetChar(rd) IN dig DO 
        Rd.UnGetChar(rd);
        WITH  enum = FGet.Int(rd) DO 
          IF eStack[enum] # NIL THEN
            eStack[enum].setNum(enum);
            eStack[enum].setStyle(FGet.Int(rd));
            eStack[enum].setGeom(SMGeo.ReadCircle(rd));
          END;
        END;

        FGet.EOL(rd);
        Lex.Skip(rd)
      END;
      Rd.UnGetChar(rd);
      FileFmt.ReadFooter(rd,"Edges Data")
    END ReadEdges;

  BEGIN
    rd:=FileRd.Open(name);

    EVAL FileFmt.ReadComment(rd,'#');  (* ---  discards the title   --- *)
    EVAL FileFmt.ReadComment(rd,'#');  (* --- discards the comments --- *)
    
    (* --- reads the element's counter  --- *) 
    cnt.NC := NGet.Int(rd,"Corners");  FGet.EOL(rd);
    self.NV:= NGet.Int(rd,"Vertices"); FGet.EOL(rd); 
    self.NE:= NGet.Int(rd,"Edges");   FGet.EOL(rd);
    cnt.NB := NGet.Int(rd,"Borders");  FGet.EOL(rd);
    cnt.NF := NGet.Int(rd,"Faces");  FGet.EOL(rd);
    EVAL FileFmt.ReadComment(rd,'#');  (* --- discards a blank line --- *)
    
    (* --- initializes the topological counter's --- *)
    self.setTCounter(cnt);

    (* --- initializes stack in --- *)
    vStack:=NEW(VStack,self.NV+1);
    eStack:=NEW(EStack,self.NE+1);
    
    (* --- reads the topology --- *)
    FileFmt.ReadHeader(rd,"Topology","");
    
    NARROW(self, SMC.Map).read(rd, vStack, eStack);
    
    EVAL FileFmt.ReadComment(rd,'#');  (* ---  discards the blank line --- *)
    FileFmt.ReadHeader(rd,"Features","");
    
    ReadVertices();
    ReadEdges();

    FileFmt.ReadFooter(rd,"Features");
  END ReadMap;


(* =================================================== *)
(*       Internal procedures used by WriteMap          *)         
(* =================================================== *)
PROCEDURE NumDigits(n: CARDINAL) : CARDINAL = 
  VAR w: CARDINAL:=1;
  BEGIN
    WHILE n>9 DO INC(w); n:=n DIV 10 END;
    RETURN w
  END NumDigits;

PROCEDURE WriteMap(self: T; name: TEXT; comment: TEXT) RAISES ANY = 
  CONST 
    HeaderVertex = "number,   style,    geometry";
    HeaderEdge   = "number,   style,    geometry";
  VAR wr  : FileWr.T;
      vWidth : CARDINAL := NumDigits(self.NV);
      eWidth : CARDINAL := NumDigits(self.NE);
      
      vStack: VStack;
      eStack: EStack;
  
  BEGIN
    wr:=FileWr.Open(name);
    IF NOT self.ValidNumber THEN 
      self.numerate()
    END;

    (* --- initializes the stack of vertices and edges --- *)
    vStack:=NEW(VStack, self.NV+1);
    eStack:=NEW(EStack, self.NE+1);

    (* --- writes title and comments --- *)
    FileFmt.WriteComment(wr,"<<<< Spherical Map - Version "&SMVersion&" >>>>",'#');
    IF NOT Text.Empty(comment) THEN 
      FileFmt.WriteComment(wr,comment,'#');
    END;

    (* --- writes the element's counter --- *) 
    WITH cnt = self.getTCounter() DO 
      NPut.Int(wr,Fmt.Pad("Corners ",8),cnt.NC);  FPut.EOL(wr);
      NPut.Int(wr,Fmt.Pad("Vertices",8),self.NV);  FPut.EOL(wr); 
      NPut.Int(wr,Fmt.Pad("Edges   ",8),self.NE);  FPut.EOL(wr);
      NPut.Int(wr,Fmt.Pad("Borders ",8),cnt.NB);  FPut.EOL(wr);
      NPut.Int(wr,Fmt.Pad("Faces   ",8),cnt.NF);  FPut.EOL(wr);
    END;
    
    (* --- writes topology --- *)
    NARROW(self, SMC.Map).write(wr, vStack, eStack);
    
    (* --- writes features of vertices and edges --- *)
    FileFmt.WriteComment(wr," ",'#');    (* --- writes a blank line --- *)
    FileFmt.WriteHeader(wr,"Features","");
    FileFmt.WriteHeader(wr,"Vertices Data","");
    FileFmt.WriteComment(wr,HeaderVertex,'#');

    (* --- writes the vertices --- *)
    FOR i:=0 TO self.NV DO 
      WITH v = vStack[i] DO  
        IF v # NIL THEN 
          <* NOWARN *> Wr.PutText(wr,"    "&Fmt.Pad(Fmt.Int(i),vWidth)&"    ");
          <* NOWARN *> Wr.PutText(wr,"    "&Fmt.Int(v.getStyle())&"      ");
          
          (* writes Cpoint's Plucker coefficients and its equivalent Spoint *)
          WITH p = v.getGeom() DO 
            SMPoint.PrintCpoint(wr,p^);
            <* NOWARN *> Wr.PutText(wr,"   <=>   ");
            SMPoint.PrintSpoint(wr,SMPoint.CToS(p^),8);
          END;
          FPut.EOL(wr)
        END
      END;
    END;
    FileFmt.WriteFooter(wr,"Vertices Data");

    (* --- writes the edges --- *)
    FileFmt.WriteHeader(wr,"Edges Data","");
    FileFmt.WriteComment(wr,HeaderEdge,'#');
    FOR i:=0 TO self.NE DO 
      WITH e = eStack[i] DO  
        IF e # NIL THEN 
          <* NOWARN *> Wr.PutText(wr,"   "&Fmt.Pad(Fmt.Int(i),eWidth)&"    ");
          <* NOWARN *> Wr.PutText(wr,"   "&Fmt.Int(e.getStyle())&"     ");
          WITH sc = e.getGeom() DO SMCircle.Print(wr,sc^) END;
          FPut.EOL(wr)
        END
      END;
    END;
    FileFmt.WriteFooter(wr,"Edges Data");
    
    FileFmt.WriteFooter(wr,"Features");
    <* NOWARN *> Wr.Close(wr);
  END WriteMap;

PROCEDURE DrawMap(self: T; name: TEXT) RAISES ANY = 
  VAR sd: SMDraw.T;
      params : SMDraw.Atributes;

      vStack: VStack;
      eStack: EStack;

  PROCEDURE DrawCorner(c: Corner) = 
    BEGIN
      IF c.isIsolatedVertex() THEN 
        WITH vn = c.org().getNum() DO 
          (* --- checks if that vertex wasn't used yet --- *)
          IF vStack[vn] = NIL THEN 
            WITH sp=SMPoint.CToS(GetVertexGeometry(c)^) DO 
              sd.drawPoint(sp,params);
            END;
            vStack[vn]:=c.org()
          END
        END
      ELSE 
        WITH en = c.halfEdge().edge().getNum() DO 
          (* --- checks if that edge wasn't used yet --- *)
          IF eStack[en] = NIL THEN 
            IF c.isOval() THEN 
              WITH circ=GetEdgeGeometry(c) DO 
                sd.drawCircle(circ^,params);
              END;
            ELSE
              sd.drawSegment(GetArcGeometry(c),params);
            END;
            eStack[en]:=c.halfEdge().edge()
          END
        END
      END
    END DrawCorner;
  
  BEGIN
    sd:=NEW(SMDraw.T).init(name);

    IF NOT self.ValidNumber THEN 
      self.numerate()
    END;

    (* --- initializes the stack of vertices and edges --- *)
    vStack:=NEW(VStack, self.NV+1);
    eStack:=NEW(EStack, self.NE+1);

    (* --- initializes params; must be replaced by customer's parametrization --- *)
    params.color:=Color{0.0,0.0,0.0};
    params.width:=0.3;
    params.colorfill:=Color{0.0,0.0,0.0};
    params.radius:=0.5;

    (* --- writes corners and builds stack with vertices, edges and faces --- *)
    self.walkOn(DrawCorner);
    
    sd.close();
  END DrawMap;

PROCEDURE Overlay(filename1,filename2 : TEXT) : T RAISES ANY = 
  CONST dig = SET OF CHAR{'0'..'9'};
  TYPE reccrn = RECORD 
                  orgnum,enum,onext,lnext : INTEGER;
                  fwd : Sign;
                  used : BOOLEAN;
                END;
  VAR m1,m2 : T;
      mr,ms : T;
      fr,fs : TEXT;
      CList : REF ARRAY OF reccrn;
      VList : REF ARRAY OF Point;  (* vertices of the smallest map *)
      EList : REF ARRAY OF Circle; (* supporting circles of the smallest map *)
      NC,NC1,NC2,vn,en : INTEGER;

  PROCEDURE GetHeader(fn: TEXT; m:T; VAR nc: INTEGER) RAISES ANY = 
    VAR rd: FileRd.T; 
    BEGIN 
      rd:=FileRd.Open(fn);

      EVAL FileFmt.ReadComment(rd,'#');  (* ---  discards the title   --- *)
      EVAL FileFmt.ReadComment(rd,'#');  (* --- discards the comments --- *)

      (* --- reads the element's counter  --- *) 
      nc := NGet.Int(rd,"Corners");  FGet.EOL(rd);
      m.NV:= NGet.Int(rd,"Vertices"); FGet.EOL(rd); 
      m.NE:= NGet.Int(rd,"Edges");   FGet.EOL(rd);
      EVAL NGet.Int(rd,"Borders");  FGet.EOL(rd);
      EVAL NGet.Int(rd,"Faces");  FGet.EOL(rd);
      EVAL FileFmt.ReadComment(rd,'#');  (* --- discards a blank line --- *)
    END GetHeader; 

  PROCEDURE GetVertexAndEdges(fn: TEXT; <* UNUSED *>m: T) RAISES ANY = 
    VAR rd: FileRd.T; 
    BEGIN 
      rd:=FileRd.Open(fn);

      EVAL FileFmt.ReadComment(rd,'#');  (* ---  discards the title   --- *)
      EVAL FileFmt.ReadComment(rd,'#');  (* --- discards the comments --- *)

      (* --- discards element's counter  --- *) 
      EVAL NGet.Int(rd,"Corners");  FGet.EOL(rd);
      EVAL NGet.Int(rd,"Vertices"); FGet.EOL(rd); 
      EVAL NGet.Int(rd,"Edges");   FGet.EOL(rd);
      EVAL NGet.Int(rd,"Borders");  FGet.EOL(rd);
      EVAL NGet.Int(rd,"Faces");  FGet.EOL(rd);
      EVAL FileFmt.ReadComment(rd,'#');  (* --- discards a blank line --- *)

      (* --- discards the number of the referece face and corner  --- *) 
      EVAL FileFmt.ReadComment(rd,'#');  (* ---  discards the blank line --- *)
      FileFmt.ReadHeader(rd,"Topology",""); (* ---  discards header line --- *)
      EVAL FileFmt.ReadComment(rd,'#');  (* ---  discards the blank line --- *)
      EVAL NGet.Int(rd,"Reference Face");  FGet.EOL(rd);
      EVAL NGet.Int(rd,"Reference Corner");  FGet.EOL(rd);
      EVAL FileFmt.ReadComment(rd,'#');  (* ---  discards the blank line --- *)
      
      (* --- reads input file and stores it in stackin --- *)
      FileFmt.ReadHeader(rd,"Corners Data","");
      EVAL FileFmt.ReadComment(rd,'#');
      FGet.Skip(rd);
      WHILE Rd.GetChar(rd) IN dig DO 
        Rd.UnGetChar(rd);
        WITH cnum=FGet.Int(rd) DO
          IF Extras.ReadNum(rd,vn) THEN 
            CList[cnum].orgnum:=vn
          ELSE 
            CList[cnum].orgnum:=0 
          END;
          IF Extras.ReadNum(rd,en) THEN 
            CList[cnum].enum:=en 
          ELSE 
            CList[cnum].enum:=0 
          END;
          CList[cnum].fwd:=FGet.Int(rd);
          EVAL FGet.Int(rd); (* --- skips border number --- *)
          CList[cnum].onext:=FGet.Int(rd);
          CList[cnum].lnext:=FGet.Int(rd);
          CList[cnum].used := FALSE;
        END;
        FGet.EOL(rd);
        FGet.Skip(rd)
      END;
      Rd.UnGetChar(rd);
      FileFmt.ReadFooter(rd,"Corners Data"); 

      (* --- discards borders data --- *)
      FileFmt.ReadHeader(rd,"Borders Data","");
      EVAL FileFmt.ReadComment(rd,'#');
      FGet.Skip(rd);
      WHILE Rd.GetChar(rd) IN dig DO 
        EVAL Rd.GetLine(rd);
        FGet.Skip(rd)
      END;
      Rd.UnGetChar(rd);
      FileFmt.ReadFooter(rd,"Borders Data"); 

      (* --- discards borders data --- *)
      FileFmt.ReadHeader(rd,"Faces Data","");
      EVAL FileFmt.ReadComment(rd,'#');
      FGet.Skip(rd);
      WHILE Rd.GetChar(rd) IN dig DO 
        EVAL Rd.GetLine(rd);
        FGet.Skip(rd)
      END;
      Rd.UnGetChar(rd);
      FileFmt.ReadFooter(rd,"Faces Data"); 
      FileFmt.ReadFooter(rd,"Topology"); 
      EVAL FileFmt.ReadComment(rd,'#');

      (* --- reads vertices's geometry --- *)
      FileFmt.ReadHeader(rd,"Features","");
      FileFmt.ReadHeader(rd,"Vertices Data","");
      EVAL FileFmt.ReadComment(rd,'#');
      FGet.Skip(rd);

      WHILE Rd.GetChar(rd) IN dig DO 
        Rd.UnGetChar(rd);
        WITH  vnum = FGet.Int(rd) DO 
          
          EVAL FGet.Int(rd);  (* --- discards style --- *)
          IF VList[vnum] = NIL THEN
            VList[vnum] := SMGeo.ReadPoint(rd)
          ELSE 
            <* NOWARN *> Wr.PutText(Stdio.stdout," ********* ERRO lendo geometria do vertice vnum = "&Fmt.Int(vnum)&" \n");
            <* NOWARN *> Wr.Flush(Stdio.stdout);
            WITH c=0 DO EVAL 1 DIV c END 
          END;
        END;

        (* --- ignores the equivalent Spoint --- *)
        EVAL Rd.GetLine(rd);
        FGet.Skip(rd)
      END;
      Rd.UnGetChar(rd);
      FileFmt.ReadFooter(rd,"Vertices Data");

      (* --- reads edge's geometry --- *)
      FileFmt.ReadHeader(rd,"Edges Data","");
      
      EVAL FileFmt.ReadComment(rd,'#');
      FGet.Skip(rd);
      WHILE Rd.GetChar(rd) IN dig DO 
        Rd.UnGetChar(rd);
        WITH  enum = FGet.Int(rd) DO 
          EVAL FGet.Int(rd);  (* --- discards style --- *)
          IF EList[enum] = NIL THEN
            EList[enum] := SMGeo.ReadCircle(rd)
          ELSE 
            <* NOWARN *> Wr.PutText(Stdio.stdout," ********* ERRO lendo geometria da aresta enum = "&Fmt.Int(enum)&" \n");
            <* NOWARN *> Wr.Flush(Stdio.stdout);
            WITH c=0 DO EVAL 1 DIV c END 
          END;
        END;

        FGet.EOL(rd);
        Lex.Skip(rd)
      END;
      Rd.UnGetChar(rd);
      FileFmt.ReadFooter(rd,"Edges Data");
      FileFmt.ReadFooter(rd,"Features")
    END GetVertexAndEdges;
    
    PROCEDURE MakeIsolatedVerticesAndArcs() RAISES ANY = 
      VAR circ : Circle;
      BEGIN
        FOR i:=1 TO NC DO 
          IF NOT CList[i].used THEN
            (* --- this corner was not created yet --- *)
            IF CList[i].fwd = 0 THEN 
              (* -- it is an isolated vertex --- *)
              EVAL mr.makeVertex(VList[CList[i].orgnum]);
            ELSE 
              WITH 
                vn = CList[i].orgnum,
                en = CList[i].enum,
                symnum = CList[CList[i].lnext].onext
              DO 
                IF CList[i].fwd > 0 THEN 
                  circ := EList[en]
                ELSE 
                  circ := SMGeo.InvCircle(EList[en])
                END;
                
                IF vn = 0 THEN (* --- it is an oval --- *)
                  EVAL mr.makeEdge(circ);
                ELSE (* --- it is an arc --- *)
                  WITH 
                    org = VList[CList[i].orgnum],
                    dst = VList[CList[symnum].orgnum]
                  DO 
                    EVAL mr.makeEdge(circ,org,dst);
                  END;
                END;
                CList[symnum].used:=TRUE
              END 
            END;
            CList[i].used:=TRUE
          END
        END
      END MakeIsolatedVerticesAndArcs;
  
  BEGIN
    m1:=NEW(T).init();
    m2:=NEW(T).init();
    
    (* --- get the maps's counter --- *)
    GetHeader(filename1,m1,NC1);
    GetHeader(filename2,m2,NC2);
    
    IF NC1 >= NC2 THEN  
      (* map "m1" has more corners than "m2" *)
      mr := m1;
      fr := filename1;
      ms := m2;
      fs := filename2;
      NC := NC2;
    ELSE 
      (* map "m2" has more corners than "m1" *)
      mr := m2;
      fr := filename2;
      ms := m1;
      fs := filename1;
      NC := NC1
    END;

    (* --- builds the biggest map --- *)
    mr.read(fr);
    
    (* --- gets isolated vertex and arcs of the smallest map --- *)
    CList:=NEW(REF ARRAY OF reccrn,NC+1);
    VList:=NEW(REF ARRAY OF Point,ms.NV+1);
    EList:=NEW(REF ARRAY OF Circle,ms.NE+1);
    
    GetVertexAndEdges(fs,ms);
    
    MakeIsolatedVerticesAndArcs();
    
    RETURN mr;
   
  END Overlay;


BEGIN
  RefPoint := NEW(InterPoint);
  RefPoint^.p:=SMGeo.NorthPole; 
  RefPoint^.which:=NIL;
  RefPoint^.coming:=NIL 
END SMap.