MODULE SMGeo;

IMPORT LR3,I3,I4,I4Extras,I6,HI3,SMPoint,SMCircle,Rd,Wr,Stdio;

REVEAL Arc = BRANDED REF RECORD
               u,v : Plane;     (* ending points *) 
               c   : Circle;    (* scircle connecting the ending points *)
             END;               

(* ------------------------------------------------------------ *)
(*   Internal procedure to check a Arc, tha is to verify        *)
(*   if the origin and destination are on the supporting circle *)
(* ------------------------------------------------------------ *)

PROCEDURE CheckArc(READONLY s : Arc) RAISES {ValidationError} =
  BEGIN
    WITH 
      p = GetDstOfArc(s),
      q = GetOrgOfArc(s)
    DO 
      IF SidePointCircle(p,s^.c) # 0 THEN 
        RAISE ValidationError("Origem nao pertence ao circulo de suporte"); 
      END;
      IF SidePointCircle(q,s^.c) # 0 THEN 
        RAISE ValidationError("Destino nao pertence ao circulo de suporte"); 
      END
    END;
  END CheckArc;

PROCEDURE AreSamePoints(p,q: Point) : BOOLEAN =
  VAR pa,qa   : SMPoint.AT;
  BEGIN
    IF p=NIL OR q=NIL THEN 
      RETURN FALSE
    ELSE 
      WITH
        okp = SMPoint.IsApoint(p^,pa),
        okq = SMPoint.IsApoint(q^,qa)
      DO
        IF okp # okq THEN 
          RETURN FALSE
        ELSIF okp THEN 
          RETURN (I3.Equal(pa.a,qa.a) = 1)
        ELSE
          RETURN (I6.Equal(p^.c,q^.c) = 1)
        END
      END
    END
  END AreSamePoints;

PROCEDURE AreSameCircles(c1,c2 : Circle) : Sign = 
  BEGIN
    RETURN SMCircle.AreCirclesEqual(c1^,c2^)  
  END AreSameCircles;
  
PROCEDURE InvStabbingLine(p: Point) : Point =
  BEGIN
    RETURN NEW(Point,c:=I6.Neg(p^.c))
  END InvStabbingLine;
  
PROCEDURE InvCircle(c: Circle) : Circle =
  BEGIN
    WITH rc = SMCircle.InvCircle(c^) DO
      RETURN NEW(Circle,s:=rc.s)
    END
  END InvCircle;

PROCEDURE AllocArc() : Arc =
  VAR s := NEW(Arc);
  BEGIN
    s^.u := NEW(Plane);
    s^.v := NEW(Plane);
    s^.c := NEW(Circle);
    RETURN s
  END AllocArc;

PROCEDURE SidePointCircle(p: Point; c: Circle) : Sign =
  BEGIN
    WITH
      pc = SMCircle.PlaneFromCircle(c^),
      lp = SMPoint.LineFromCpoint(p^),
      q  = HI3.PointFromLineAndPlane(lp,pc)
    DO
      IF HI3.IsInvalidPoint(q) THEN 
        RETURN 0
      ELSE
        WITH 
          sq = HI3.SideOfSphere(q),
          sm = HI3.Side(HI3.MiddlePoint(lp),pc),
          sd = HI3.Side(HI3.LineDir(lp),pc)
        DO 
          CASE sq OF 
          | -1 => RETURN sd
          |  0 => IF sm = sd THEN RETURN sm
                  ELSE RETURN 0 END
          | +1 => RETURN sm
          END
        END
      END 
    END
  END SidePointCircle;

PROCEDURE MakeArc(p,q : Point; c: Circle) : Arc RAISES ANY =
  BEGIN
    WITH s = NEW(Arc) DO
      IF p= NIL AND q = NIL THEN (* --- is a ring --- *)
        s^.u:=NIL;
        s^.v:=NIL
      ELSE 
        s^.u:=PlaneDefPointOnCircle(p,c);
        IF p=q THEN (* --- is a loop --- *)
          s^.v:=s^.u
        ELSE 
          s^.v:=PlaneDefPointOnCircle(q,c);
        END;
      END;
      s^.c:=c;
      RETURN s
    END;
  END MakeArc;

PROCEDURE GetPointOfArc(p: Plane; c: Circle) : Point =
  (* Returns the point corresponding to the intersection of the 
     plane "p" with the circle "c". *)
  BEGIN
    WITH
      r  = HI3.LineFromTwoPlanes(SMCircle.PlaneFromCircle(c^),p^),
      np = NEW(Point)
    DO
      np^:=SMPoint.Stab(SMPoint.CpointFromLine(r));
      RETURN np
    END 
  END GetPointOfArc;

PROCEDURE GetOrgOfArc(s: Arc) : Point =
  BEGIN
    IF s^.u = NIL THEN 
      RETURN NIL
    ELSE
      RETURN GetPointOfArc(s^.u,s^.c)
    END
  END GetOrgOfArc;

PROCEDURE GetDstOfArc(s: Arc) : Point =
  BEGIN
    (* in a loop where "origin=destination", "s^.v" can be NIL and 
       so the destination is given by "s^.u" *)
    IF s^.u = NIL THEN 
      RETURN NIL
    ELSIF s^.v = NIL THEN 
      RETURN GetPointOfArc(s^.u,s^.c)
    ELSE       
      RETURN GetPointOfArc(s^.v,s^.c)
    END
  END GetDstOfArc;

PROCEDURE ScircleOfArc(s: Arc) : Circle = 
  BEGIN
    RETURN s^.c
  END ScircleOfArc;

PROCEDURE InvArc(s: Arc) : Arc RAISES ANY =
  BEGIN
    RETURN MakeArc(GetDstOfArc(s),GetOrgOfArc(s),InvCircle(s^.c))
  END InvArc;

PROCEDURE PlaneDefPointOnCircle(p: Point; c: Circle) : Plane RAISES ANY =
  VAR np := NEW(Plane);
  BEGIN
    WITH q=BringStabLineToCircle(p,c) DO 
      np^ := HI3.PlaneOrthogToPlaneThroughLine(SMCircle.PlaneFromCircle(c^),
                                               SMPoint.LineFromCpoint(q^));
    END;
    RETURN np
  END PlaneDefPointOnCircle;

PROCEDURE BringStabLineToCircle(p: Point; c: Circle): Point RAISES {ValidationError} =
  (* Given a point "p" on the \Scircle "c", returns the same point 
     represented by a stabbing line that lies on the plane "splane(c)" *)
  VAR pa : SMPoint.AT;
  BEGIN
    IF HI3.IsLineOnPlane(SMPoint.LineFromCpoint(p^),
                         SMCircle.PlaneFromCircle(c^)) THEN
      RETURN p
    ELSE
      IF SMPoint.IsApoint(p^,pa) THEN 
        WITH
          hpa = SMPoint.ApointToHomogPoint(pa),
          lp  = HI3.LineFromTwoPoints(SMCircle.DCenter(c^),hpa),
          np = NEW(Point)
        DO 
          np^:=SMPoint.CpointFromLine(lp);
          RETURN np
        END
      ELSE
        RAISE ValidationError("BringStabLineToScircle: line is on circle and point is not Apoint")
      END
    END
  END BringStabLineToCircle;

PROCEDURE CircOrd3PointsOnCircle(p,q,r: Point; c: Circle) : Sign RAISES ANY = 

  PROCEDURE Ahead2Cpoints(u,v : Point; c: Circle) : Sign RAISES ANY =
    VAR dm,dc : Sign;
        cb,cdl,cbref := NEW(Circle);
    BEGIN
      WITH
        l  = SMPoint.LineFromCpoint(u^),        
        b  = HI3.PlaneFromLineAndPoint(l,SMCircle.Snormal(c^)),
        O = HI3.Origin
      DO
        cb^:=SMCircle.CircleFromPlane(b);

        dc:=HI3.Side(O,b);

        IF dc = 0 THEN 
          RETURN -1*SidePointCircle(v,cb)
        END;
        CASE dc * SidePointCircle(v,cb) OF 
        | -1 => dm:=+1;
        |  0 => IF AreSamePoints(u,v) THEN dm:=0
                ELSE dm:=+1 END;
        | +1 => BEGIN
                  cdl^:=SMCircle.CircleFromPlane(HI3.SPolarOfPoint(HI3.LineDir(l)));

                  IF SidePointCircle(v,cdl) = -1 THEN 
                    cbref^:=SMCircle.CircleFromPlane(HI3.ReflectPlaneAcrossOrg(b));
                    dm:=dc*SidePointCircle(v,cbref)
                  ELSE
                    dm:=-1
                  END
                END
        ELSE 
        END;
        RETURN dc*dm
      END
    END Ahead2Cpoints;
    
  BEGIN
    WITH
      spq = Ahead2Cpoints(p,q,c),
      sqr = Ahead2Cpoints(q,r,c),
      srp = Ahead2Cpoints(r,p,c),
      s   = spq+sqr+srp
    DO
      IF s = 0 THEN 
        RETURN 0
      ELSE 
        RETURN (s DIV ABS(s))
      END
    END
  END CircOrd3PointsOnCircle;

PROCEDURE IsPointInArc(p: Point; sarc: Arc; VAR sp: Sign) : BOOLEAN RAISES ANY = 
  BEGIN
    IF p = NIL THEN 
      RETURN FALSE
    ELSE 
      IF sarc^.c = NIL THEN (* arc is an isolated vertex *)
        sp:=+1;
        RETURN AreSamePoints(p,GetOrgOfArc(sarc))
      ELSE 
        sp := SidePointCircle(p,sarc^.c);
        IF sp # 0 THEN  
          (* "p" isn't on the supporting circle *) 
          RETURN FALSE 
        ELSIF sarc^.u = NIL AND sarc^.v = NIL THEN (* arc is an entire circle *)
          RETURN TRUE 
        ELSE
          WITH
            org = GetOrgOfArc(sarc),
            dst = GetDstOfArc(sarc),
            s = CircOrd3PointsOnCircle(org,p,dst,sarc^.c)
          DO 
            IF s = 0 THEN 
              IF AreSamePoints(p,org) THEN 
                sp:=+1 
              ELSIF AreSamePoints(p,dst) THEN 
                sp:=-1 
              ELSE 
                sp:=0
              END
            END;
            RETURN s >= 0 
          END
        END        
      END
    END
  END IsPointInArc;

PROCEDURE IntersectionOfTwoArcs(A,B: Arc; perturb: Sign; VAR tangent: BOOLEAN; 
                                          VAR p: Point) : BOOLEAN RAISES ANY =
  (*    
    To deal with degenerated cases (such as, overlapping and
    intersection with extremities) we assume that the arc "A" is
    perturbed by \epsilon around its origin. The direction
    (counterclockwise or clockwise) of this perturbation is defined by
    the parameter "perturb"; in fact, "pertub \in {-1,0,1}" where
        . "pertub = -1" means clockwise pertubation, 
        . "pertub = +1" means counterclockwise pertubation,
        . "pertub =  0" means "DON'T DO PERTUBATION"
    
    The parameter "tangent" returns whether the two arcs are tangent
    or not. *)
    
                                          
  VAR sa, sb: Sign;
      q : Point;
      BAUX : Arc;
  PROCEDURE CheckPerturbation(q: Point; aa,bb: Arc) : BOOLEAN RAISES ANY = 
    (* -------------------------------------------------------------------- *)
    (* checks if the point "q = aa \meet bb" lies on the arc "aa" perturbed *)
    (* (rotated) by  "\epsilon" around its origin in the direction given by *)
    (* "perturb" being that "pertub=+1" implies a counterclockwise rotation,*)
    (* "-1" clockwise and "0" implies DO NOT consider pertubation           *)
    (* -------------------------------------------------------------------- *)
    BEGIN
      IF perturb = 0 THEN 
        RETURN TRUE
      ELSE 
        (* --------------------------------------------------------- *)
        (* breaks the arc "aa=(u,v,c)" at point "q" creating arcs    *)
        (* "sa1=(q,v,c)" and "sa2=(q,u,\neg c)". So, checks if the   *)
        (* circular order of "sa1","bb","sa2" at point "q" coincides *)
        (* with the pertubation's direction. See                     *)
        (*                                                           *)
        (*                         bb                                *)
        (*                    u    |    v                            *)
        (*  (original aa)     .----*----.                            *)
        (*                   sa2   q   sa1                           *)
        (* --------------------------------------------------------- *)
        WITH 
          org = PlaneDefPointOnCircle(q,aa^.c),
          sa1 = NEW(Arc,u:=org,v:=aa^.v,c:=aa^.c),
          sa2 = NEW(Arc,u:=org,v:=aa^.u,c:=InvCircle(aa^.c)),
          sord  = CircOrdOf3ArcsAroundVertex(sa1,bb,sa2,q)
        DO 
          RETURN sord = perturb
        END      
      END
    END CheckPerturbation;

  BEGIN
    q := NEW(Point);
    IF SMCircle.Intersection(A^.c^,B^.c^,q^,tangent) THEN 
      IF SMPoint.IsCpointNull(q^) THEN  
        (* --------------------------------------------------------------- *) 
        (* if "q" is "<0,0,0,0,0,0>" then "A.c" and "B.c" are coincident.  *) 
        (* By the pertubation adopted (rotation of "A" around its origin)  *)
        (* there is a intersection iff "A.org" is on the arc "B"           *)
        (* --------------------------------------------------------------- *) 
        WITH Aorg = GetOrgOfArc(A) DO 
          IF IsPointInArc(Aorg,B,sb) THEN 
            p := Aorg;
            RETURN TRUE 
          END;
        END;
      ELSIF IsPointInArc(q,A,sa) AND IsPointInArc(q,B,sb) THEN 
        IF sb = 0  THEN 
          (* "q" is internal to "B"; so any pertubation on "A" 
             DOESN'T modify the intersection *)
          p := q;   
          RETURN TRUE 
        ELSIF sb = 1 THEN (* "q" is the origin of "B" *)
          BAUX:=B
        ELSE (* "q" is the destination of "B" *)
          BAUX:=InvArc(B);
        END;
        IF CheckPerturbation(q,A,BAUX) THEN  
          p:=q;
          RETURN TRUE 
        END
      END;
    END;
    RETURN FALSE 
  END IntersectionOfTwoArcs;
    
PROCEDURE GeneratePointOnCircle(c: Circle) : Point =
  VAR p := NEW(Point);
  BEGIN
    WITH
      dc = SMCircle.DCenter(c^),
      q = HI3.Point{c:=I4Extras.GenerateOrtho(c^.s)}
    DO
      p^ := SMPoint.Stab(SMPoint.CpointFromLine(HI3.LineFromTwoPoints(dc,q)));
      RETURN p
    END 
  END GeneratePointOnCircle;

PROCEDURE CirclePassingThroughPoint(p: Point) : Circle RAISES ANY =
  VAR a,b := NEW(Plane);
      c := NEW(Circle);
  BEGIN
    HI3.TwoPlanesByLine(SMPoint.LineFromCpoint(p^),a^,b^);
    c^ := SMCircle.CircleFromPlane(a^);
    RETURN c
  END CirclePassingThroughPoint;

PROCEDURE CircleByTwoPoints(p,q: Point; VAR sc : Circle) : BOOLEAN RAISES ANY = 
  VAR ap : SMPoint.AT;
  BEGIN
    sc := NEW(Circle);
    WITH
      lp = SMPoint.LineFromCpoint(p^),
      lq = SMPoint.LineFromCpoint(q^)
    DO
      IF I6.Equal(lp.k,lq.k) # 0 THEN 
        (* "lp" and "lq" are coincident or opposite *)
        sc:=CirclePassingThroughPoint(p);
        RETURN TRUE 
      END;
      IF SMPoint.IsApoint(q^,ap) THEN 
        WITH b = HI3.PlaneFromLineAndPoint(lp,SMPoint.ApointToHomogPoint(ap)) DO 
          IF NOT I4.IsAllZero(b.f) THEN 
            sc^:=SMCircle.CircleFromPlane(b);
            RETURN TRUE 
          END 
        END
      END;
      IF SMPoint.IsApoint(p^,ap) THEN 
        WITH a = HI3.PlaneFromLineAndPoint(lq,SMPoint.ApointToHomogPoint(ap)) DO 
          IF NOT I4.IsAllZero(a.f) THEN 
            sc^:= SMCircle.CircleFromPlane(a);
            RETURN TRUE 
          END 
        END
      END;
      
      (* --- "p" and "q" are not A-points OR "p" is on Line(q) OR vice-versa --- *)
      IF HI3.AreLinesParallel(lp,lq) # 0 THEN 
        (* "lp" and "lq" are parallel or coparallel but not coincident *)

        sc^:=SMCircle.CircleFromPlane(
                 HI3.PlaneFromLineAndPoint(lp,HI3.MiddlePoint(lq)))
      ELSE
        WITH pq = HI3.FrontMeetOfTwoLines(lp,lq) DO
          IF NOT HI3.IsInvalidPoint(pq) THEN (* "lp" meets "lq" *)
            IF HI3.SideOfSphere(pq) >= 0 THEN (* "pq" isn't inside S2 *)
              sc^:=SMCircle.CircleFromPlane(
                       HI3.PlaneFromLineAndPoint(lp,HI3.MiddlePoint(lq)))
            ELSE
              sc^:=SMCircle.CircleFromPlane(
                      HI3.PlaneFromLineAndPoint(lp,HI3.LineDir(lq)))
            END
          ELSE
            (* there isn't a rational circle passing through these points *)
            sc:=NIL;
            RETURN FALSE 
          END
        END
      END;
    END;
    RETURN TRUE 
  END CircleByTwoPoints;

PROCEDURE CircOrdOf3ArcsAroundVertex(A,B,C : Arc; v: Point) : Sign RAISES ANY =
  VAR pa : SMPoint.AT;
  BEGIN
    WITH
      a = SMCircle.PlaneFromCircle(A^.c^),
      b = SMCircle.PlaneFromCircle(B^.c^),
      c = SMCircle.PlaneFromCircle(C^.c^),
      l = SMPoint.LineFromCpoint(v^)
    DO
      IF NOT SMPoint.IsApoint(v^,pa) THEN 
        (* The line "l" is not tangent and is common to the planes "a","b","c"  *)
        RETURN HI3.CircOrdOf3PlanesAroundALine(a,b,c,l)
      ELSE           
        (* The tangent plane "v*" is rational *)
        WITH
          v_ast = HI3.SPolarOfPoint(SMPoint.ApointToHomogPoint(pa)),
          ta    = HI3.LineFromTwoPlanes(a,v_ast),
          tb    = HI3.LineFromTwoPlanes(b,v_ast),
          tc    = HI3.LineFromTwoPlanes(c,v_ast),
          eqab  = I6.Equal(ta.k,tb.k),
          eqbc  = I6.Equal(tb.k,tc.k),
          eqca  = I6.Equal(tc.k,ta.k),
          dta   = HI3.LineDir(ta),
          dtb   = HI3.LineDir(tb),
          dtc   = HI3.LineDir(tc)
        DO
          IF  eqab # 1 AND eqbc # 1 AND eqca # 1 THEN 
            (* The three tangent lines are different *)
            WITH
              Omega = HI3.Plane{f:=I4.Axis(0)},
              l_inf = HI3.LineFromTwoPlanes(Omega,v_ast)
            DO
              RETURN HI3.CircOrdOf3PointsOnALine(dta,dtb,dtc,l_inf)
            END 
          ELSIF eqab = 1 AND eqbc = 1 AND eqca = 1 THEN 
            (* The three tangent lines are equal *)
            RETURN HI3.CircOrdOf3PlanesAroundALine(a,b,c,ta)
          ELSE
            (* Only two tangent lines are equal (or opposite) *)
            IF eqab  = 1 THEN (* ta = tb and tb # tc *)
              RETURN HI3.Ahead2PlanesAroundLine(a,b,dta)
            ELSIF eqbc = 1 THEN (* tb = tc and ta # tc *) 
              RETURN HI3.Ahead2PlanesAroundLine(b,c,dtb) 
            ELSE (* ta = tc and ta # tb *) 
              RETURN HI3.Ahead2PlanesAroundLine(c,a,dtc) 
            END
          END
        END
      END 
    END 
  END CircOrdOf3ArcsAroundVertex;

PROCEDURE ReadPoint(rd: Rd.T) : Point RAISES ANY = 
  BEGIN
    WITH p=NEW(Point) DO
      p^:=SMPoint.ReadCpoint(rd);
      RETURN p
    END   
  END ReadPoint;

PROCEDURE ReadCircle(rd: Rd.T) : Circle RAISES ANY = 
  BEGIN
    WITH sc=NEW(Circle) DO
      sc^:=SMCircle.Read(rd);
      RETURN sc
    END   
  END ReadCircle;
  
PROCEDURE ReadArc(rd: Rd.T; eps: LONGREAL) : Arc RAISES ANY = 
  VAR s  : Arc;
      sep: TEXT;
      bp : SMPoint.CT;
  BEGIN
     s := AllocArc();
     IF rd = Stdio.stdin THEN 
       Wr.PutText(Stdio.stdout,"Entre com o scircle de suporte = ");
       Wr.Flush(Stdio.stdout);
     END;
     s^.c^:=SMCircle.Read(rd);
     WITH
       pc = SMCircle.PlaneFromCircle(s^.c^),
       pcd= HI3.PolarComplementOfPlane(pc)
     DO
       IF rd = Stdio.stdin THEN 
         Wr.PutText(Stdio.stdout,"Entre com XYZ (reais) do ponto inicial = ");
         Wr.Flush(Stdio.stdout);
       END;
       bp:=SMPoint.SToC(SMPoint.ReadSpoint(rd),eps);
       s^.u^:=HI3.PlaneFromLineAndPoint(SMPoint.LineFromCpoint(bp),pcd);
       IF rd = Stdio.stdin THEN 
         Wr.PutText(Stdio.stdout,"Entre com XYZ (reais) do ponto final   = ");
         Wr.Flush(Stdio.stdout);
       END;
       bp:=SMPoint.SToC(SMPoint.ReadSpoint(rd),eps);
       s^.v^:=HI3.PlaneFromLineAndPoint(SMPoint.LineFromCpoint(bp),pcd);
       IF rd # Stdio.stdin THEN 
         sep:=Rd.GetLine(rd)
       END;
     END;
     RETURN s
  END ReadArc;

PROCEDURE PrintArc(wr: Wr.T; s: Arc) RAISES ANY =  
  VAR a,b : Point;
  BEGIN
    Wr.PutText(wr,"Arco => ");
    a:=GetOrgOfArc(s);
    b:=GetDstOfArc(s);
    IF a#NIL AND b#NIL THEN (* --- if is not a ring --- *)
      WITH
        pa=SMPoint.CToS(a^),
        pb=SMPoint.CToS(b^)
      DO
        LR3.Print(wr,pa.s,"A = (",",",")   ");
        LR3.Print(wr,pb.s,"B = (",",",")   ");
      END;    
    END;
    SMCircle.Print(wr,s^.c^);
    Wr.Flush(wr);
  END PrintArc;

BEGIN
  (* ------------------------------ *)
  (* sets Cpoint to "<0,0,1,0,0,0>" *)
  (* ------------------------------ *)
  NorthPole:=NEW(Point);
  NorthPole^:=SMPoint.CT{c:=I6.Axis(2)}
END SMGeo.