BeginPackage["HPolygons`", {"HDraw`"}]

(* 
  
  This package defines some routines for calculating and visualizing
  Euclidean and hyperbolic geometry in the unit disk.
  
  Requires Mathematica version 2.0 or greater.
  
  The newest version of this package is available for anonymous ftp at
  nic.funet.fi, directory "pub/math/riemann/mathematica".
  
  Version 0.1. Last modified: Nov 1st, 1994.
  
  Send comments and bug reports to 
  
    Juha Haataja (e-mail: jhaataja@csc.fi)
    Center for Scientific Computing
    PO Box 405
    FIN-02101 Espoo, Finland
  
*)

(*----------------------------------------------------------------------*)

Unprotect[UnitPolygon, RegularPolygon, FindCenter, 
  CircleAroundPolygon, Diagonals, Bisectors, FilterPolygons, 
  HypInvertPolygon, PolygonInvGroup];

(*----------------------------------------------------------------------*)

(* Descriptions of functions *)

UnitPolygon::usage = 
"UnitPolygon[n] forms a unit polygon of given number of sides.";

RegularPolygon::usage = 
"RegularPolygon[n,angle] constructs a regular polygon with a given
number of sides and a given angle between the sides.";

FindCenter::usage =
"FindCenter[Polygon[pts]]";

CenterPoint::usage =
"CenterPoint[Polygon[pts]]";

CircleAroundPolygon::usage =
"CircleAroundPolygon[Polygon[pts]]";

Diagonals::usage =
"Diagonals[Polygon[pts]]";

Bisectors::usage =
"Bisectors[Polygons[pts]]";

FilterPolygons::usage =
"FilterPolygons[polygonlist]";

PolygonLayer::usage =
"";

AddPolygonLayer::usage =
"";

HypInvertPolygon::usage = 
"HypInvertPolygon[polygon, {ind1, ind2, ind3}]. Invert a set of
points (see also Polygon) across a given edge or a set of edges.";

PolygonInversion::usage = 
"PolygonInversion[Polygon[{p1,p2,p3}]]";

PolygonInvGroup::usage = 
"PolygonInvGroup[Polygon[{p1,p2,p3}], level]. Make an inversion group
for a polygon.";

(*----------------------------------------------------------------------*)

Begin["`Private`"]

(*----------------------------------------------------------------------*)

(* Routines for numerical ordering of points at the center of polygons *)

NPointLessQ[{_, Point[{p1x_,p1y_}]}, {_, Point[{p2x_,p2y_}]}] :=
  If[p1x < p2x - GeomEps, True, 
    If[p1x > p2x + GeomEps, False,
      If[p1y < p2y - GeomEps, True, 
        False
  ]]];

NPointLessQ[Point[{p1x_,p1y_}], Point[{p2x_,p2y_}]] :=
  If[p1x < p2x - GeomEps, True, 
    If[p1x > p2x + GeomEps, False,
      If[p1y < p2y - GeomEps, True, 
        False
  ]]];

(*----------------------------------------------------------------------*)

(* Construct a polygon of wanted size *)

UnitPolygon[sides_, rotation_:0, scaling_:1] :=
  Polygon[Table[Point[{scaling*Sin[rotation + i*2*Pi/sides],
    scaling*Cos[rotation + i*2*Pi/sides]}], {i,0,sides-1}]];

(* Construct a regular polygon *)

RegularPolygon[sides_, angle_, rotation_:0]:=
  Module[{beta = Pi/sides, dist, i, alfa},
    alfa = If[SameQ[angle,Indeterminate], Pi/4*(1-6/sides), angle/2];
    dist = Hyp2EuclDist[Abs[ArcCosh[Cot[alfa]*Cot[beta]]]];
    UnitPolygon[sides, rotation, dist]];

(*----------------------------------------------------------------------*)

(* Find the hyperbolic center of a regular polygon *)

FindCenter[Polygon[{p1_,p2_,p3_,___}]] := BSolve[p1,p2,p3];

CenterPoint[pg:Polygon[{p1_,p2_,p3_,___}]] := 
  CenterPoint[pg] = Point[BSolve[p1,p2,p3]];

(*----------------------------------------------------------------------*)

CircleAroundPolygon[poly:Polygon[pts_List]] :=
  Module[{cp, pt = First[pts[[1]]]},
    cp = FindCenter[poly];
    Circle[cp, EuclDist[cp,pt]]
  ];

(* Make the arcs from polygon corners to the centre *)

Diagonals[poly:Polygon[pts_List]] :=
  Module[{c, i, len = Length[pts]},
    If[Mod[len,2]==0, 
      Table[Line[pts[[i]],pts[[i+len/2]]], {i, 1, len/2}],
      c = CenterPoint[poly];
      Table[Line[pts[[i]], c], {i, 1, len}]]]

(* Make the arcs from middle points of polygon edges to the centre *)

Bisectors[poly:Polygon[pts_List]] :=
  Module[{sides, centers, i, len = Length[pts]},
    sides = Map[Apply[Line,#]&, Thread[{pts,RotateLeft[pts]}]];
    centers = Map[NHypMiddle,sides];
    If[Mod[len,2]==0,
      Table[Line[centers[[i]],centers[[i+len/2]]], {i, 1, len/2}],
      Table[Line[centers[[i]], pts[[(1+Mod[i-1+Quotient[len,2],len])]]],
        {i,1,len}]
    ]
  ];

(*----------------------------------------------------------------------*)

(* Invert a polygon across an given edge *)

HypInvertPolygon[Polygon[pts_List], index_Integer] :=
  Module[{pt = RotateLeft[pts, index - 1], arc, inv},
    arc = Take[pt,2]; 
    inv = Inversion[Line @@ arc]; 
    Polygon[Join[arc, Map[inv, Drop[pt,2]]]]
  ];

(* Invert a polygon across some of its edges *)

HypInvertPolygon[poly:Polygon[pts_List], indexes_List] :=
    Map[HypInvertPolygon[poly,#]&, indexes];

(* Invert a polygon across all its edges *)

HypInvertPolygon[poly:Polygon[pts_List]] :=
  Map[HypInvertPolygon[poly,#]&, Table[i,{i,1,Length[pts]}]];

(* Functions for inversions *)

PolygonInversion[Polygon[pts_List], index_Integer] :=
  Module[{pt = RotateLeft[pts, index - 1], arc},
    arc = Take[pt,2]; 
    Inversion[Line @@ arc]
  ];

PolygonInversion[poly:Polygon[pts_List], indexes_List] :=
    Map[PolygonInversion[poly,#]&, indexes];

PolygonInversion[poly:Polygon[pts_List]] :=
  Map[PolygonInversion[poly,#]&, Table[i,{i,1,Length[pts]}]];

(*----------------------------------------------------------------------*)

(* Eliminate duplicates from a list of polygons  *)

FilterPolygons[pls_List] :=
  Module[{p = Flatten[pls], plist},
    plist = Thread[{p,Map[CenterPoint, p]}];
    plist = Sort[plist, NPointLessQ];
    plist = FixedPoint[Replace[#, {H___,poly1:{Polygon[__], Point[p1_]}, 
      {Polygon[__], Point[p2_]}, E___} /; 
        SamePtQ[p1,p2] :> {H, poly1, E}]&, plist];
    Map[#[[1]]&, plist]
  ];

(* Make an inversion group for a polygon *)

PolygonLayer[polygons_List] :=
  FilterPolygons[Map[HypInvertPolygon, polygons]]

PolygonLayer[p:Polygon[__]] :=
  FilterPolygons[HypInvertPolygon[p]]

AddPolygonLayer[pgl_List, index_] :=
  FilterPolygons[{pgl,Map[HypInvertPolygon[#,index]&, pgl]}]

AddPolygonLayer[pg:Polygon[__], index_] :=
  FilterPolygons[{pg,HypInvertPolygon[pg,index]}]

PolygonInvGroup[pg:Polygon[pts_List], level_Integer:1] :=
  If[level<=0, pg,
    Module[{index = Table[i,{i,1,Length[pts]}], res},
      res = Nest[AddPolygonLayer[#,index]&, pg, level]
    ]
  ];

PolygonInvGroup[pgl_List, level_Integer:1] :=
  If[level<=0, FilterPolygons[pgl],
      FilterPolygons[Map[PolygonInvGroup[#,level]&, Flatten[pgl]]]
  ];

(*----------------------------------------------------------------------*)

End[];

(*----------------------------------------------------------------------*)

Protect[UnitPolygon, RegularPolygon, FindCenter, 
  CircleAroundPolygon, Diagonals, Bisectors, FilterPolygons, 
  HypInvertPolygon, PolygonInvGroup];

(*----------------------------------------------------------------------*)

EndPackage[];

(*----------------------------------------------------------------------*)


