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

(* 
  
  This package defines some routines for calculating and visualizing
  Euclidean and hyperbolic geometry in the unit disk.
  
  This package was developed for Mathematica version 2.2.
  
  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
  
*)

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

(* Descriptions of functions *)

GeomPrecision::usage =
"The internal precision of floating-point calculations.";

GeomEps::usage =
"The internal threshold for regarding two floating-point numbers as 
equal."; 

NGP::usage =
"NGP[args] converts the arguments into floating-point numbers using
GeomPrecision digits.";

ZToReIm::usage = 
"ZToReIm[z] converts a complex number z to a list of real and imaginary
parts.";

ReImToZ::usage = 
"ReImToZ[{x,y}] converts a list of two numbers to a complex number.";

EuclDist::usage = 
"EuclDist[Point[p1],Point[p2]] calculates the Euclidean distance of
two points.";

NEuclDist::usage = 
"EuclDist[Point[p1],Point[p2]] calculates numerically the Euclidean
distance of two points.";

HypDist::usage = 
"HypDist[Point[p1],Point[p2]] calculates the hyperbolic distance of
two points in the unit disk model.";

NHypDist::usage = 
"HypDist[Point[p1],Point[p2]] calculates numerically the hyperbolic
distance of two points in the unit disk model.";

Hyp2EuclDist::usage = 
"Hyp2EuclDist[dist] converts the hyperbolic distance to Euclidean
distance from the origin in the unit disk model.";

MovePoint::usage = 
"MovePoint[dist, theta]. Move a complex number a given hyperbolic
distance from origin to a given direction.";

ParallelQ::usage =
"ParallelQ[p1,p2] tests whether points p1 and p2 lie on the
same line through the origin.";

SamePtQ::usage =
"SamePtQ[p1,p2] tests whether the floating-point values of point p1
and point p2 are approximately equal.";

Angle::usage = 
"Angle[z1] finds the angle of a complex number in relation to
the origin.";

RotatePoint::usage = 
"RotatePoint[z1, angle]. Rotate a complex number around origin.";

ZPoint::usage =
"ZPoint[x+I*y] results in Point[{x,y}], which can be used as an
argument to the routines Line, Polygon etc.";

Line::usage =
"Line[Point[p1],Point[p2]] defines a hyperbolic or Euclidean line
segment between Point[p1] and Point[p2].";

Geodesic::usage = 
"Geodesic[Point[p1],Point[p2]] defines a hyperbolic geodesic going
through the two given points from the unit disk to the unit disk.";

UnitCircle::usage = 
"Example: HDraw[UnitCircle]. Defines the unit circle for calculations 
and drawing.";

Polygon::usage =
"Polygon[{Point[p1],Point[p2],Point[p3]}] defines a polygon (in the
example a triangle, which can be drawn as a hyperbolic or Euclidean
polygon.";

FilledPolygon::usage =
"PlainFilledPolygon defines for drawing a filled and bordered polygon.
See also Polygon.";

PlainFilledPolygon::usage =
"PlainFilledPolygon defines for drawing a filled polygon.
See also Polygon.";

HypMiddle::usage =
"HypMiddle[Line[p1,p2]] finds out the hyperbolic middle point of line
segment from point p1 to point p2. See also NHypMiddle.";

NHypMiddle::usage =
"NHypMiddle[Line[p1,p2]] finds out numerically the hyperbolic middle
point of line segment from point p1 to point p2.";

Invert::usage = 
"Invert[Circle[{cx,cy},r], Point[p1]] inverts the point p1 accross
the given circle.";

EuclReflect::usage = 
"Eucl[Line[p1,p2], Point[p1]] inverts the point p1 accross
the given Euclidean line.";

Inversion::usage =
"Inversion[Line[p1,p2]] generates of function for inverting points
across geodesics and lines.";

HSolve::usage =
"HSolve[line1,line2] or HSolve[line, circle] or HSolve[circle,circle]
finds the common points of lines and circles.";

HypSolve::usage =
"HypSolve is an auxilliary routine for the HSolve command.";

BSolve::usage = 
"BSolve[p1, p2, p3] finds out the center of the circle which goes
through points p1, p2 and p3.";

HDraw::usage =
"HDraw[{objects}, options] draws the given geometric objects using the
unit disk model of hyperbolic geometry. See also Options[HDraw].";

HypDrawRules::usage =
"HypDrawRules is the rule set of the command HDraw for drawing
points, lines and geodesics in the hyperbolic unit disk model.
See also EuclDrawRules.";

EuclDrawRules::usage =
"EuclDrawRules is the rule set of the command HDraw for drawing
points, lines and geodesics in the Euclidean model. See also
HypDrawRules.";

Make::usage = 
"Make[obj] converts lines, polygons and circles for drawing.";

FilledArc::usage =
"FilledArc[p1,p2] calculates a set of successive points on the
geodesic arc segment from the point p1 to the point p2.";

ArcSegments::usage =
"ArcSegments[p1,p2] is an auxilliary routine for the FilledArc routine.";

PointLabel::usage =
"PointLabel[Point[p], text] labels a point. See also Options[PointLabel].";

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

Options[HDraw] = {DrawUnitCircle -> True, Hyperbolic -> True};

Options[PointLabel] = {DefaultFont -> {"Times-Roman", 9}};

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

(* Descriptions of options for functions *)

DrawUnitCircle::usage =
"DrawUnitCircle is an option of HDraw. If True, the unit circle is drawn; 
otherwise not.";

Hyperbolic::usage =
"Hyperbolic is an option of HDraw. If True, the hyperbolic geometry model
is used for drawing the figures.";

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

Begin["`Private`"]

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

(* General utilities *)

(* Tools for filtering options *)

FilterOptions[command_Symbol, opts___] :=
  Module[{keywords = First /@ Options[command]},
    Sequence @@ Select[{opts}, MemberQ[keywords, First[#]]&]]

(* Flatten a list to two nested levels *)

FlattenObjectList[lst_List] :=
  FixedPoint[Replace[#1, {A___, {X___List}, B___} :>
    {A, X, B}] &, lst]

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

(* The precision of floating-point calculations *)

GeomPrecision = $MachineEpsilon;  
GeomEps = 1000*GeomPrecision;
GeomDigits = Floor[1-Log[10,GeomPrecision]];

NGP[arg_] := N[arg, GeomDigits];
NCGP[arg_] := Chop[N[arg, GeomDigits], 10^(-7)];

SetAttributes[NGP, Listable];
SetAttributes[NCGP, Listable];

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

(* Utilities *)

(* Complex number conversions *)

ZToReIm[Point[z_]] := {Re[z],Im[z]};
ZToReIm[z_] := {Re[z],Im[z]};
ReImToZ[{x_,y_}] := x + I*y;

ZPoint[z_] := Point[ZToReIm[z]]

ConjI[{a_,b_}] := {-b,a};  

(* Coordinate transformations *)

RotatePoint[p:Point[{x_,y_}], angle_] :=
  Point[{x*Cos[angle] - y*Sin[angle], y*Cos[angle] + x*Sin[angle]}]

Angle[Point[{x_,y_}]] := Arg[x+y*I];

PolarCoord[Point[a_]] := PolarCoord[a];  

PolarCoord[{x_,y_}] :=
  Module[{r=Sqrt[x^2+y^2], t, sgn, ang},
    t = NGP[ArcSin[y/r]];
    sgn = Map[(Sign[NGP[{x,y}]]==#)&, 
      {{1,1}, {-1,1}, {-1,-1}, {1,-1}, {0,-1}, {0,1}, {-1,0}, {1,0}}];
    ang = {t,Pi-t,Pi-t,2*Pi+t,t,t,Pi,t};
    t = ang[[Flatten[Position[sgn,True]]]];
    Flatten[{r,t}]
  ];  

ParallelQ[Point[p1_], Point[p2_]] := ParallelQ[p1,p2];

ParallelQ[{x_,y_},{z_,w_}] := Abs[NGP[x*w-y*z]] < GeomEps;   

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

(* Distance functions *)

EuclDist[Point[{a1_,a2_}],Point[{b1_,b2_}]] := 
  Simplify[Sqrt[Expand[(a1-b1)^2+(a2-b2)^2]]];

EuclDist[{a1_,a2_},{b1_,b2_}] := 
  Simplify[Sqrt[Expand[(a1-b1)^2+(a2-b2)^2]]];

NEuclDist[Point[{a1_,a2_}],Point[{b1_,b2_}]] := 
  Sqrt[NGP[(a1-b1)^2+(a2-b2)^2]];

NEuclDist[{a1_,a2_},{b1_,b2_}] := Sqrt[NGP[(a1-b1)^2+(a2-b2)^2]];

HypDist[Point[{a1_,a2_}],Point[{b1_,b2_}]] := 
  Module[{r1=a1^2+a2^2, r2=b1^2+b2^2, r3=a1*b1+a2*b2},
    Abs[ArcCosh[((1+r1)*(1+r2)-4*r3)/((1-r1)*(1-r2))]]
  ];

NHypDist[Point[{a1_,a2_}],Point[{b1_,b2_}]] := 
  Module[{r1=a1^2+a2^2, r2=b1^2+b2^2, r3=a1*b1+a2*b2},
    Abs[NGP[ArcCosh[((1+r1)*(1+r2)-4*r3)/((1-r1)*(1-r2))]]]
  ];

(* Numerical comparison of points *)

SamePtQ[Point[{a1_,a2_}], Point[{b1_, b2_}]] := 
  NGP[(a1-b1)^2+(a2-b2)^2] < GeomEps;

SamePtQ[{a1_,a2_}, {b1_,b2_}] := 
  NGP[(a1-b1)^2+(a2-b2)^2] < GeomEps;

(* Move a complex number hyperbolic distance dist from origin to
direction angle *)

MovePoint[dist_,angle_:0] :=
  Point[Hyp2EuclDist[dist]*{Cos[angle], Sin[angle]}];

Hyp2EuclDist[dist_] :=
  If[NGP[dist]==Infinity, 1, (Exp[dist] - 1)/(Exp[dist] + 1)];

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

(* Find the hyperbolic center of an arc segment *)

NHypMiddle[arc:Line[Point[p1:{p1x_,p1y_}],Point[p2:{p2x_,p2y_}]]] :=
  HypMiddle[NGP[arc]]

HypMiddle[Point[p1_], Point[p2_]] :=
  HypMiddle[Line[Point[p1],Point[p2]]];

NHypMiddle[Point[p1_], Point[p2_]] :=
  HypMiddle[Line[Point[NGP[p1]],Point[NGP[p2]]]];

HypMiddle[arc:Line[Point[p1:{p1x_,p1y_}],Point[p2:{p2x_,p2y_}]]] :=
  Module[{m,n,s,t,x,y,x0,y0,r0,eq,sol,res,k,b},
    {m,n} = HypSolve[Make[arc],Circle[{0,0},1]];
    s = EuclDist[p1,n]/EuclDist[p1,m];
    t = EuclDist[p2,m]/EuclDist[p2,n];
    If[!ParallelQ[p1,p2],
      (* a geodesic circle *)
      {x0,y0} = First[Make[arc]]; 
      r0 = EuclDist[{x0,y0},p1]; 
      eq = (x - x0)^2 + (y - y0)^2 == r0^2;
      sol = Solve[{s*((x - m[[1]])^2 + (y - m[[2]])^2) ==
        t*((x - n[[1]])^2 + (y - n[[2]])^2), eq}, {x, y}],
      (* line through origin *)
      If[Abs[NGP[p1x-p2x]] > GeomEps,
        k = (p1y-p2y)/(p1x-p2x); b = p1y - k*p1x;
        y = k*x + b;
        sol = Solve[s*((x - m[[1]])^2 + (y - m[[2]])^2) ==
          t*((x - n[[1]])^2 + (y - n[[2]])^2), x],
        sol = Solve[s*((p1x - m[[1]])^2 + (y - m[[2]])^2) ==
          t*((p1x - n[[1]])^2 + (y - n[[2]])^2), y]]
    ];
    res = Map[({x,y} /. #)&, sol];
    Simplify[Point @@ Select[res, (NGP[#[[1]]^2+#[[2]]^2] <= 1)&]]
  ];

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

(* Geometric Transformations *)

BSolve[Point[p1_], Point[p2_], Point[p3_]] :=
  BSolve[p1, p2, p3];

BSolve[p1:{p1x_,p1y_}, p2:{p2x_,p2y_}, p3:{p3x_,p3y_}] :=
  Module[{b1, b2, vars = {x,y}, c},
    b1 = If[NEuclDist[p1,p2] < GeomEps, 
      p1x*y - p1y*x == 0,
      (-p1x^2-p1y^2+p2x^2+p2y^2)/2 + (p1x-p2x)*x + (p1y-p2y)*y == 0];
    b2 = If[NEuclDist[p2,p3] < GeomEps, 
      p2x*y - p2y*x == 0,
      (-p2x^2-p2y^2+p3x^2+p3y^2)/2 + (p2x-p3x)*x + (p2y-p3y)*y == 0];
    {c} = vars /. Solve[{b1,b2}, vars];
    c
  ];

Invert[{x_,y_}] := {x,y}/(x^2+y^2);             

Invert[Circle[{x_,y_},r_], Point[{px_,py_}]] := 
  Point[Invert[Circle[{x,y},r], {px,py}]]

Invert[cr:Circle[{x_,y_},r_], l:Line[Point[p1_],Point[p2_]]] := 
  Map[Invert[cr,#]&, l, {2}]

Invert[cr:Circle[{x_,y_},r_], g:Geodesic[Point[p1_],Point[p2_]]] := 
  Map[Invert[cr,#]&, g, {2}]

Invert[cr:Circle[{x_,y_},r_], pl:Polygon[pts_List]] := 
  Polygon[Map[Invert[cr,#]&, pts]]

Invert[Circle[{x_,y_},r_], {px_,py_}] := 
  Module[{denom = px^2 + py^2 - 2*px*x + x^2 - 2*py*y + y^2},
    {px*r^2 - 2*px*x^2 + x^3 + x*(px^2 + py^2 - r^2 - 2*py*y + y^2), 
      py*r^2 + (px^2 + py^2 - r^2)*y - 2*px*x*y + x^2*y - 2*py*y^2 + y^3}/
      denom
  ];   

EuclReflect[Line[Point[p_],Point[q_]], obj_]:=
  EuclReflect[Line[{p,q}], obj];

EuclReflect[l:Line[lst_List], Point[p_]]:= Point[EuclReflect[l, p]];

EuclReflect[l:Line[lst_List], Line[p1_,p2_]]:= 
  Line[EuclReflect[l,p1], EuclReflect[l,p2]];

EuclReflect[l:Line[lst_List], Geodesic[p1_,p2_]]:= 
  Geodesic[EuclReflect[l,p1], EuclReflect[l,p2]];

EuclReflect[l:Line[lst_List], Polygon[pts_List]]:= 
  Polygon[Map[EuclReflect[l,#]&, pts]];

EuclReflect[Line[{{px_,py_},{qx_,qy_}}], {p1_,p2_}]:=
  Module[{tmp, denom, s1, s2},
    denom = (-px + qx)^2 + (py - qy)^2;
    tmp = (-(py*p1) + px*p2 + py*qx - p2*qx - px*qy + p1*qy)/denom;
    s1 = p1 + 2*(py - qy)*tmp;
    s2 = p2 + 2*(-px + qx)*tmp;
    {s1,s2}
  ];

Inversion[Geodesic[Point[a_],Point[b_]]] := 
  Inversion[Line[Point[a],Point[b]]]
  
Inversion[Line[Point[p1_],Point[p2_]]] :=  
  Module[{vars = {x,y}, b1, b2, c, d, r},
    If[ParallelQ[p1,p2],
      Return[Function[obj, EuclReflect[Line[{p1,p2}], obj]]]
    ];
    c = BSolve[p1, Invert[p1], p2];
    r = First[PolarCoord[p1-c]];
    Function[obj, Invert[Circle[c,r], obj]]
  ];

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

(* Solve circles and lines *)

HypSolve[Circle[p1:{p1x_,p1y_},r1_,___], Circle[p2:{p2x_,p2y_},r2_,___]] :=
  Module[{dx=p2x-p1x, dy=p2y-p1y, len, t0, t1},
    len = dx^2 + dy^2;
    t0 = len+r1^2-r2^2;
    t1 = 4*len*r1^2-t0^2;
    If[NGP[t1] < 0, Return[{}]];
    t1 = Sqrt[t1];
    len = 2*len;
    {{p1x + (dx*t0 - dy*t1)/len, p1y + (dy*t0 + dx*t1)/len}, 
     {p1x + (dx*t0 + dy*t1)/len, p1y + (dy*t0 - dx*t1)/len}}
  ];

HypSolve[Circle[c_,r_,th___],Line[{a_,b_}]] := 
  HypSolve[Line[{a,b}], Circle[c,r,th]];

HypSolve[Line[{p1:{p1x_,p1y_},p2:{p2x_,p2y_}}], Circle[c:{cx_,cy_},r_,___]] := 
  Module[{dx=p2x-p1x, dy=p2y-p1y, ex=cx-p1x, ey=cy-p1y, len, t0, t1},
    len = dx^2 + dy^2;
    t0 = {(dx*(dx*ex + dy*ey))/len, (dy*(dx*ex + dy*ey))/len};
    t1 = r^2 -(ex-dx)^2 - (ey-dy)^2;
    If[NGP[t1] < 0, Return[{}]];    
    t1 = Sqrt[t1/len];
    {{dx + p1x + dx*t1, dy + p1y + dy*t1}, 
     {dx + p1x - dx*t1, dy + p1y - dy*t1}}
  ];

HypSolve[Line[{p1:{p1x_,p1y_},p2:{p2x_,p2y_}}], 
         Line[{q1:{q1x_,q1y_},q2:{q2x_,q2y_}}]] :=
  Module[{dpx=p2x-p1x, dpy=p2y-p1y, dqx=q2x-q1x, dqy=q2y-q1y, d, n1, n2},
    d = dqy*(p1x - p2x) + dqx*(-p1y + p2y);
    n1 = p2x*(-(dqx*p1y) + q1y*q2x - q1x*q2y) + 
      p1x*(dqx*p2y - q1y*q2x + q1x*q2y);
    n2 = -(dqy*p1y*p2x) + dqy*p1x*p2y + p2y*(q1y*q2x - q1x*q2y) + 
      p1y*(-(q1y*q2x) + q1x*q2y);
    {n1/d, n2/d}
  ];

HSolveRules = {
  UnitCircle -> Circle[{0,0},1],
  Line[Point[p1_List],Point[p2_List]] :> Line[{p1,p2}],
  Circle[Point[p1_List],Point[p2_List]] :> 
    Circle[p1, EuclDist[p1, p2]]
};

HSolve[obj1_,obj2_] := 
  Module[{t1, t2, pts},
    t1 = (obj1 /. HSolveRules);
    t2 = (obj2  /. HSolveRules);
    pts = HypSolve[t1, t2];
    First[Map[Point, Select[pts, (NGP[Dot[#,#]] <= 1)&]]]
  ];

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

(* Conversions to graphics primitives *)

CommonDrawRules = {UnitCircle -> Circle[{0,0},1]};

HypDrawRules = Join[CommonDrawRules, {
    Line[Point[p1_List],Point[p2_List]] :> 
      Make[Line[Point[p1],Point[p2]]],
    Geodesic[Point[p1_List],Point[p2_List]] :> 
      Make[Geodesic[Point[p1],Point[p2]]],
    Circle[Point[p1_List],Point[p2_List]] :> 
      Make[Circle[Point[p1],Point[p2]]],
    Circle[c_,r_,o___] :> Circle[c,r,o],
    Polygon[pts_List] :> Make[Polygon[pts]],
    FilledPolygon[pts_List] :> Make[FilledPolygon[pts]],
    PlainFilledPolygon[pts_List] :> Make[PlainFilledPolygon[pts]]
  }];

EuclDrawRules = Join[CommonDrawRules, {
    Line[Point[p1_List],Point[p2_List]] :> Line[{p1,p2}],
    Circle[Point[p1_List],Point[p2_List]] :> 
               Circle[p1, NEuclDist[p1, p2]],
    Circle[c_,r_,o___] :> Circle[c,r,o],
    Polygon[pts_List] :> EuclPolygon[pts],
    FilledPolygon[pts_List] :> {FilledEuclPolygon[pts], 
               GrayLevel[0], EuclPolygon[pts]},
    PlainFilledPolygon[pts_List] :> FilledEuclPolygon[pts]
  }];

Make[Line[Point[p1:{p1x_,p1y_}], Point[p2:{p2x_,p2y_}]]] := 
  Module[{pt1, pt2, b1, b2, c, d, r, tmp1, tmp2},
    If[ParallelQ[p1,p2], Return[Line[{p1,p2}]]];
    pt1 = NGP[If[NGP[p1x^2+p1y^2] < 1.0, 
      p1, p1*(1-1000*GeomEps)]];
    pt2 = NGP[If[NGP[p2x^2+p2y^2] < 1.0, 
      p2, p2*(1-1000*GeomEps)]];
    c = BSolve[pt1, Invert[pt1], pt2];
    {r,tmp1} = PolarCoord[pt1-c];
    tmp2 = (PolarCoord[pt2-c])[[2]];
    {tmp1,tmp2} = Sort[{tmp1,tmp2}, OrderedQ[{NCGP[#1],NCGP[#2]}]&];
    If[NCGP[Pi-(tmp2-tmp1)] < 0, {tmp1,tmp2} = {tmp2-2*Pi,tmp1}];
    Circle[NGP[c], NGP[r], NCGP[{tmp1,tmp2}]]
  ];

Make[Geodesic[Point[p1:{p1x_,p1y_}], Point[p2:{p2x_,p2y_}]]] := 
  Module[{vars = {x,y}, b1, b2, c1, c2, d, r2, r, tmp1, tmp2, alpha, beta},
    If[ParallelQ[p1,p2], t = (p2-p1)/Sqrt[(p1x-p2x)^2+(p1y-p2y)^2];
      Return[Line[{-t,t}]]];
    {c1,c2} = BSolve[p1, Invert[p1], p2];
    r2 = (p1x-c1)^2+(p1y-c2)^2;
    {alpha,beta} = HypSolve[Circle[{0,0},1], Circle[{c1,c2},Sqrt[r2]]];
    {r,tmp1} = PolarCoord[alpha-{c1,c2}];
    tmp2 = (PolarCoord[beta-{c1,c2}])[[2]];
    {tmp1,tmp2} = Sort[{tmp1,tmp2}, OrderedQ[{NCGP[#1],NCGP[#2]}]&];
    If[NCGP[Pi-(tmp2-tmp1)] < 0, {tmp1,tmp2} = {tmp2-2Pi,tmp1}];
    Circle[NGP[{c1,c2}], NGP[r], NCGP[{tmp1,tmp2}]]
  ];

EuclPolygon[pts_List] := 
  Module[{pl = Flatten[pts]},
    Line[Join[pl,{First[pl]}] /. Point[p_] :> p]
  ];

FilledEuclPolygon[pts_List] := Polygon[pts] /. Point[p_] :> p

FilledArc[{p1_, p2_}] := 
  ((Make[Line[p1, p2]] /. Circle -> ArcSegments) /. Line[{pts__}] :> {pts})

ArcSegments[c_List, r_, {phi1_, phi2_}] := 
  Table[NCGP[{r*Cos[phi], r*Sin[phi]} + c], 
    {phi, NGP[phi1], NGP[phi2], Sign[NGP[phi2-phi1]]*
                       Abs[NGP[phi2-phi1]]/Floor[10+5*r]}]

CombineSegments[l1_List, l2_List] :=
  If[SamePtQ[Last[l1],First[l2]], FlattenObjectList[{l1,l2}],
    If[SamePtQ[First[l1],First[l2]], FlattenObjectList[{Reverse[l2], l1}],
      If[SamePtQ[First[l1],Last[l2]], FlattenObjectList[{l2,l1}],
         FlattenObjectList[{l1,Reverse[l2]}]]]]

Make[Polygon[pts_List]] :=
  Module[{l = Flatten[pts], arcs, t},
    sides = Map[Apply[Line,#]&, Thread[{l,RotateLeft[l]}]];
    Flatten[Map[Make,sides],1]
  ];

Make[FilledPolygon[pts_List]] :=
    {Make[PlainFilledPolygon[pts]], GrayLevel[0], Make[Polygon[pts]]}

Make[PlainFilledPolygon[pts_List]] :=
  Module[{l = Flatten[pts], arcs, t},
    l = Join[l,{First[l]}];
    arcs = Partition[l, 2, 1];
    l = Map[FilledArc, arcs];
    t = Fold[CombineSegments, First[l], Rest[l]];
    Polygon[t]
  ];

Make[Circle[Point[p1:{p1x_,p1y_}], Point[p2:{p2x_,p2y_}]]] := 
  Module[{hl=Make[Line[Point[p1],Point[p2]]], 
          vars = {x,y}, tmp, c, b1, b2},
    If[SameQ[Line,Head[hl]], (* then *)
      If[SamePtQ[p1,{0,0}], Return[Circle[p1,NEuclDist[p1,p2]]]];
      c = Simplify[Expand[(p1+Invert[p1])/2]];
      tmp = Invert[Circle[c,NEuclDist[p1,c]]][p2];
      c = Simplify[Expand[(p2+tmp)/2]],
      c = hl[[1]]; (* else *)
      b1 = Expand[(vars-p2).(p2-c)] == 0;
      b2 = Expand[vars.ConjI[p1]] == 0;
      {c} = vars /. Solve[{b1,b2}, vars]
    ];
    Circle[c,NEuclDist[c,p2]]
  ];   

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

(* Label functions *)

PointLabel[Point[p_], label_, opts___Rule] := 
  Module[{ofs, style},
    style = DefaultFont /. {opts} /. Options[PointLabel];
    Text[FontForm[label, style], p]
  ];

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

(* Drawing Function *)    

HDraw[geom_, opts___Rule] := 
  Module[{hypq, drawrules, ucir, obj, pltrange},
    hypq = Hyperbolic /. {opts} /. Options[HDraw];
    drawrules = If[hypq, HypDrawRules, EuclDrawRules];
    ucir = DrawUnitCircle /. {opts} /. Options[HDraw];
    obj = If[ucir, NGP[{geom, Thickness[0.004], UnitCircle}] /. drawrules,
                   NGP[geom] /. drawrules];
    pltrange = PlotRange /. {opts} /. {PlotRange -> All};
    Show[Graphics[obj], AspectRatio -> Automatic, PlotRange -> pltrange,
      FilterOptions[Graphics, opts]]
  ];

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

End[];

EndPackage[];

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