(* Copyright (c) 1992 The Geometry Center; University of Minnesota
   1300 South Second Street;  Minneapolis, MN  55454, USA;
   
This file is part of CirclePack. CirclePack is free software; you can
redistribute it and/or modify it only under the terms given in the
file COPYING, which you should have received along with this file.
This and other software may be obtained via anonymous ftp from
geom.umn.edu; email: software@geom.umn.edu. *)

(* Authors: Silvio Levy, Oliver Goodman *)
(* $Id: Hyperbolic.m,v 1.4 1992/07/17 16:12:57 oag Exp $ *)

BeginPackage["Hyperbolic`", "Minkowski`", "Normalize`", "Circle3D`"]

(* NOTE: This package is under development.
   The documentation is terse.  Read the whole thing once, to see
   what's available; also look at the enclosed transcript. *)

UnitCircle::usage = "\n  UnitCircle draws a unit circle at the origin."

ModelQ::usage = "\n
  ModelQ[keyword] tests whether the keyword is one of the five\n
  models known to the package: Minkowski, Projective, Hemisphere,\n
  PoincareBall, and UpperHalfSpace."

Projective::usage = "\n
  Projective[x], where x is a list of length n, describes a point in\n
  hyperbolic n-space, in the projective model. Projective[Vector[x,y]]\n
  describes a tangent vector at x, whose tip is at x+y."
  
PoincareBall::usage = "\n
  PoincareBall[x] and PoincareBall[Vector[x,y]] are like their\n
  counterparts with head Projective (q.v.)"
  
UpperHalfSpace::usage = "\n
  UpperHalfSpace[x] and UpperHalfSpace[Vector[x,y]] are like their\n
  counterparts with head Projective (q.v.)"
  
Minkowski::usage = "\n
  Minkowski[x], where x is a list of length n+1, describes a point in\n 
  hyperbolic n-space, in the hyperboloid model.  The conversion routines\n
  assume x has length i, that is, x_1^2+...+x_n^2-x_(n+1)^2=-1.\n
  Minkowski[Vector[x,y]] describes a tangent vector at x, whose tip is\n
  at x+y.  The conversion routines assume x has length i and y is\n
  orthogonal to x, that is, x_1y_1+...+x_n y_n-x_(n+1)y_(n+1)=0."
  
Hemisphere::usage = "\n
  Hemisphere[x] and Hemisphere[Vector[x,y]] are like their\n
  counterparts with head Projective (q.v.)"

Convert::usage = "\n
  Convert[x,(options)], where x is of type Projective, PoincareBall,\n
  etc., converts x to another model, specified by the Model option.\n
  The default Model is Automatic, which specifies no conversion at all.\n
  Convert[x,model] is the same as Convert[x,Model->model]."
 
Vector::usage = "\n
  Vector[x,y] denotes a Euclidean vector at x of length y.  It is\n
  normally used as an argument to Projective, PoincareBall, etc."
  
Distance::usage = "\n
  Distance[model1[x],model2[y]], where x and y are lists, gives the\n
  hyperbolic distance between the two points."

VectorLength::usage = "\n
  VectorLength[model[Vector[x,y]]] gives the hyperbolic length of the\n
  given tangent vector."

Angle::usage = "\n
  Angle[model1[Vector[x,y1]],model2[Vector[x,y2]]] gives the angle\n
  between two tangent vectors at the same point."

Drag::usage = "\n
  Drag[model[Vector[x,y]]] gives the vector obtained by parallel\n
  transport of Vector[x,y] along the geodesic it determines, by a\n
  distance equal to the length of y. Drag[model[Vector[x,y]],d]\n
  drags Vector[x,y] a distance d."

Geodesic::usage = "Geodesic[model[x], model[y]] gives the infinite 
  geodesic through the points x, y."

IntersectionPoints::usage = "IntersectionPoints[object1, object2] gives 
  a list of points in the Projective model for all the intersections 
  between the two objects. The objects may be Geodesics of Circles"

(* WARNING: the syntax for the next three commands is likely to change *)

Frame::usage = "\n
  Minkowski[Frame[x]], where x is an (n+1)\times(n+1) matrix,\n
  represents an orthonormal basis for (n+1)-dimensional Minkowski\n
  space. The last element should have length -1, the others length 1,\n
  and all inner products should be 0.  In hyperbolic space this\n
  corresponds to giving a point and an orthonormal basis for the\n
  tangent space at that point."

CompleteFrame::usage = "\n
  CompleteFrame[model[Vector[x,y]],sign_:1] yields a frame having\n
  x and y (or their counterparts in Minkowski space) in the last\n
  and last-but-one positions.  The optional argument sign gives\n
  the sign of the frame's determinant"
  
Isometry::usage = "\n
  Isometry[frame1,frame2], where frame1 and frame2 are of the form\n
  Minkowski[Frame[x]], gives the unique isometry that takes frame1\n
  to frame2.  Isometry[,] . model[x], where x is a point or a vector,\n
  gives the image of x under the isometry.  Isometries can be composed\n
  using the . operator, and inverted using Inverse[]."

Peel::usage = "\n
  Peel[x] removes the head of the expression x, which should have\n
  only one element."

(* GRAPHICS COMMANDS *)

ResolveHyperbolic::usage = "\n
  ResolveHyperbolic is a set of rules that resolves objects of type\n
  PoincareBall, Projective and UpperHalfSpace when they occur inside\n
  a Point or Line. It is automatically called by Display, Show, etc."

ResolveHyperbolic::toohigh = "Cannot render line in dimension > 3"

(* END OF DOCUMENTATION *)

Options[Convert] = {Model->Automatic}

Begin["`private`"]

Cross[x_,y_] := RotateLeft[RotateLeft[y] x - RotateLeft[x] y]

Peel[_[x_]] := x

UnitCircle=Circle[{0,0},1]

Models = {Minkowski, Projective, Hemisphere, PoincareBall, UpperHalfSpace}
ModelQ[x_] := Or @@ (SameQ[x,#]& /@ Models)
Angle[x_List,y_List] := ArcCos[ x.y / Sqrt[(x.x)(y.y)] ]

Convert[source_?ModelQ[x_?VectorQ],target_?ModelQ] := 
  target[convert[x,source,target]]
Convert[x:(_?ModelQ[_?VectorQ]),opts___] := 
  Block[{m=Model/.{opts}/.Options[Convert]},
    If[SameQ[m,Automatic], x, Convert[x,m]]]
Convert[model_?ModelQ[Frame[m_?MatrixQ]], target_?ModelQ]:= 
   target[Frame[convert[#, model, target]& /@ m]]

convert[x_,target_?ModelQ,target_] := x (* no conversion needed *)
convert[x_,source_?ModelQ,target_?ModelQ] := 
  convert[convert[x,source,Hemisphere],Hemisphere,target] /;
    !SameQ[source,Hemisphere] && !SameQ[target,Hemisphere]

convert[Vector[x_?VectorQ,v_?VectorQ],source_,target_] :=
  Vector[convert[x,source,target],ConvertAt[v,x,source,target]]

convert[x_List,Projective,Hemisphere] := Append[x,Sqrt[1-x.x]]
ConvertAt[v_List,x_List,Projective,Hemisphere] := Append[v,-x.v/Sqrt[1-x.x]]
convert[x_List,Hemisphere,Projective] := Drop[x,-1]
ConvertAt[v_List,x_List,Hemisphere,Projective] := Drop[v,-1]

convert[x_List,Minkowski,Hemisphere] := MapAt[0&,x,{-1}]/x[[-1]] /; 
  Abs[Minkowski[x].Minkowski[x]]<.1 (* i.e. is 0 *)
convert[x_List,Minkowski,Hemisphere] := MapAt[1&,x,{-1}]/x[[-1]]
ConvertAt[v_List,x_List,Minkowski,Hemisphere] :=
  Append[Drop[v,-1]x[[-1]]-Drop[x,-1]v[[-1]],-v[[-1]]]/x[[-1]]^2
convert[x_List,Hemisphere,Minkowski] := Append[Drop[x,-1],1] /; x[[-1]]==0
convert[x_List,Hemisphere,Minkowski] := Append[Drop[x,-1],1]/Abs[x[[-1]]]
ConvertAt[v_List,x_List,Hemisphere,Minkowski] :=
  ConvertAt[v,x,Minkowski,Hemisphere] 

convert[x_List,PoincareBall,Hemisphere] := Block[{p=1+x.x},Append[2x,2-p]/p]
ConvertAt[v_List,x_List,PoincareBall,Hemisphere] :=
  Block[{p=1+x.x},Append[2v/p,0]-(4x.v/p^2)Append[x,1]]
convert[x_List,Hemisphere,PoincareBall] := Drop[x,-1]/(1+x[[-1]])
ConvertAt[v_List,x_List,Hemisphere,PoincareBall] :=
  Drop[v,-1]/(1+x[[-1]])-v[[-1]]Drop[x,-1]/(1+x[[-1]])^2

convert[x_List,UpperHalfSpace,Hemisphere] := 
  Block[{p=1+x.x},Insert[2x,-2+p,-2]/p]
ConvertAt[v_List,x_List,UpperHalfSpace,Hemisphere] :=
  Block[{p=1+x.x},Insert[2v/p,0,-2]-(4x.v/p^2)Insert[x,-1,-2]]
convert[x_List,Hemisphere,UpperHalfSpace] := Drop[x,{-2}]/(1-x[[-2]])
ConvertAt[v_List,x_List,Hemisphere,UpperHalfSpace] :=
  Drop[v,{-2}]/(1-x[[-2]])+v[[-2]]Drop[x,{-2}]/(1-x[[-2]])^2

(* GEOMETRIC QUANTITIES AND CONSTRUCTIONS *)

Distance[x_,y_] := ArcCosh[- Convert[x,Minkowski] . Convert[y,Minkowski]]
Angle[x_,y_] := Block[
  {v1=Last /@ Convert[x,Minkowski], v2=Last /@ Convert[y,Minkowski]},
  ArcCos[v1.v2 / Sqrt[(v1.v1)(v2.v2)]]
] /; x[[1,1]]==y[[1,1]] (* base points must be the same *)

(* visual midpoint in PoincareBall model *)
vmid[x_PoincareBall,y_PoincareBall] := Block[
   {p = Convert[x,Projective], q = Convert[y,Projective]},
   Convert[Projective[(Peel[p]+Peel[q])/2],PoincareBall]]

Mid[x:(m_?ModelQ)[_],y:(_?ModelQ)[_]] := Block[
   {p = Convert[x,Minkowski], q = Convert[y,Minkowski]},
   Convert[Minkowski[(Peel[p]+Peel[q])/Sqrt[2 - 2 p.q]],m]]

(* Abbreviated Quadratic Forms for the equations of circles and geodesics *)
AQForm[Geodesic[Projective[p_], Projective[q_]]]:= 
   Flatten[{{0,0,0},Cross[Append[u,1],Append[v,1]]}, 1]
AQForm[Geodesic[{m1_?ModelQ[p_], m2_?ModelQ[q_]}]]:=
   AQForm[Geodesic[Convert[m1[p], Projective], Convert[m2[q], Projective]]]
AQForm[Circle[m_?ModelQ[c_],r_]]:=
   AQForm[Circle[Convert[m[c], Projective],r]]
AQForm[Circle[Projective[{a_,b_}],r_]]:= 
   Block[{k = Cosh[r]^2(1 - a^2 - b^2)},
   {a^2 + k, 2 a b, b^2 + k, -2 a, -2 b, 1 - k}]
AQForm[UnitCircle] = {1,0,1,0,0,-1}

(* Test for types which can be converted to AQForm *)
AQType[Circle[_,_]] = True
AQType[Geodesic[_,_]] = True
AQType[UnitCircle] = True
AQType[_] = False

(* We can use AQForms to find the intersection points of two objects *)
IntersectionPoints[o1_?AQType, o2_?AQType]:= Block[{x,y}, 
   Select[Projective[{x,y} /. #]& /@ 
      Solve[{{x^2, x*y, y^2, x, y, 1}.AQForm[o1] == 0, 
         {x^2, x*y, y^2, x, y, 1}.AQForm[o2] == 0}, {x,y}],
      (Im[#[[1]]]=={0,0})&]]

(* GRAPHICAL DISPLAY *)

Unprotect[Display]
Display[x_,y_/;Or @@ (!FreeQ[y,#]& /@ Models)] := 
  Display[x,y /. ResolveHyperbolic /. ResolveExtras]
Protect[Display]

ResolveHyperbolic = {
  Point[x:(_?ModelQ[_])] :> Point @@ Convert[x],
  Text[a_,x:(_?ModelQ[_]),c___] :> Text[a,Peel[Convert[x]],c],
  g:Geodesic[_?ModelQ[_], _?ModelQ[_]]:> 
    Line[IntersectionPoints[g, UnitCircle]],
  Line[x:{(_?ModelQ[___])..}] :> Block[{xx=Convert /@ x},
    If[SameQ[Head[xx[[1]]],Projective], Line[Peel /@ xx],
      DoLine[#,False]& /@ Transpose[{Drop[xx,-1],Drop[xx,1]}]]],
  Polygon[x:{(_?ModelQ[___])..}] :> Block[{xx=Convert /@ x},
    If[SameQ[Head[xx[[1]]],Projective], Polygon[Peel /@ xx],
      Polygon[Flatten[
        ((DoLine[#,True] /. ResolveExtras)[[1]])& /@ 
        Transpose[{xx,RotateLeft[xx]}],1]]]],
  ob_?((#===Circle||#===Disk)&)[PoincareBall[v_],r_] :> Block[
    {r0 = 2 ArcTanh[Sqrt[v.v]], a, b},
    a = Tanh[(r0 - r)/2];b = Tanh[(r0 + r)/2];
    ob[If[r0 == 0, v 0, v (a + b)/(2 Sqrt[v.v])],(b - a)/2]]}

(* This set of rules turns Circles and Ellipses into appropriate 
   Line type objects. They need to be part of the display rules. 
*)
ResolveExtras = {
  Ellipse[c_, ax1_, ax2_] :>
    Line[Table[Cos[t] ax1 + Sin[t] ax2 + c, {t,0, N[2 Pi], N[Pi/30]}]], 
  Circle[c_,radius_,{beg_,end_}] :> Block[
    {tmp=end-beg, foo=(MaxBend/.Options[ParametricPlot])Degree//N, n, i},
    n=Ceiling[Abs[tmp]/foo];
    Line[Table[{Cos[beg+i tmp/n],Sin[beg+i tmp/n]} radius + c , {i,0,n}]]]}

mindot[x_,y_] := x.y - 2x[[-1]]y[[-1]]

DoLine[{PoincareBall[x_],PoincareBall[y_]},orderflag_] :=
  Block[{x0=Append[2x,1+x.x],y0=Append[2y,1+y.y],sol},
    sol=NullSpace[{mindot[#,x0],mindot[#,y0],-Last[#]}&/@{x0,y0}];
    If[Length[sol]>1,
      Line[{x,y}],
      Block[{ctr=MapAt[sol[[1,3]]+#&, sol[[1,1]] x0 + sol[[1,2]] y0 ,-1]},
        If[Max[Abs[ctr]] < 10^-6 || Abs[ctr[[-1]] ]*10^6 < Max[Abs[ctr]],
          Line[{x,y}],
          circ[Drop[ctr,-1]/ctr[[-1]],{x,y},orderflag]]]]]

circ[center_,{x_,y_},orderflag_]:=
  Switch[Length[center],
    2, Circle[center,Sqrt[center.center-1],
              fixangles[ArcTan@@(#-center)& /@ {x,y},orderflag]],
    3, Block[{x0=x-center,y0=y-center,ang,yy},
               ang=ArcCos[x0.y0/x0.x0];
               yy=(y0-x0 Cos[ang])/Sin[ang];
               Circle3D[center,{x0,yy},{0,ang}]],
    _, Message[ResolveHyperbolic::toohigh]]

DoLine[{x_UpperHalfSpace,y_},orderflag_] :=
  Block[{x1,x2,y1,y2},
    {{x1,x2},{y1,y2}}=Peel/@ {x,y};
    If[Abs[x1-y1]*10^6<Abs[x1+y1+x2^2-y2^2],Line[{{x1,x2},{y1,y2}}],
      Block[{z=(x1+y1+(x2^2-y2^2)/(x1-y1))/2},
        Circle[{z,0},Sqrt[(z-x1)^2+x2^2],
          fixangles[ArcTan@@(#-{z,0})& /@ {{x1,x2},{y1,y2}},orderflag]]]]]

(* fixangles adjusts c by a multiple of 2Pi such that a and c differ by 
   at most Pi. If orderflag is False then a and c are swapped if necessary 
   so that a-c is positive. 
*)
fixangles[{a_,c_},orderflag_] := If[orderflag,
  Switch[Floor[N[(c-a)/Pi]],1,{a,c-N[2Pi]},0,{a,c},-1,{a,c},-2,{a,c+N[2Pi]}],
  Switch[Floor[N[(c-a)/Pi]],1,{c-N[2Pi],a},0,{a,c},-1,{c,a},-2,{a,c+N[2Pi]}]]

Drag[Minkowski[Vector[x_,v_]],delta_] :=
  Block[{r=Sqrt[Minkowski[v].Minkowski[v]],s=Sinh[delta],c=Cosh[delta]},
    Minkowski[Vector[c x + s v/r, s r x + c v]]]
Drag[model_?ModelQ[v_],delta_] := 
  Convert[Drag[Convert[model[v],Minkowski],delta],model] /; !SameQ[model,Minkowski]
Drag[v_] := Drag[v,VectorLength[v]]

VectorLength[x_] := Block[{v=Last /@ Convert[x,Minkowski]}, Sqrt[(v.v)]]

CompleteFrame[Minkowski[Frame[x_?MatrixQ]],sign_:1] :=
  If[Length[x]==Length[x[[1]]],
    If[Det[x] sign>0,Minkowski[Frame[x]],
      Minkowski[Frame[MapAt[Minus,x,{{1}}]]]],
    CompleteFrame[Minkowski[Frame[Prepend[x,
      Peel[Normalize[Minkowski[First[
        NullSpace[Transpose[MapAt[Minus,Transpose[x],{{-1}}]]]]]]]]]],sign]]

CompleteFrame[Minkowski[Vector[x_?VectorQ,y_?VectorQ]],sign_:1] :=
  CompleteFrame[Minkowski[Frame[{Peel[Normalize[Minkowski[y]]],x}]],sign]

CompleteFrame[model_?ModelQ[x_],sign_:1] :=
  CompleteFrame[Convert[model[x],Minkowski],sign]

Isometry[x:Minkowski[Frame[_?MatrixQ]],y:Minkowski[Frame[_?MatrixQ]]] :=
  Isometry[y] . Inverse[Isometry[x]]

Isometry[Minkowski[Frame[x_?MatrixQ]]] := Isometry[Transpose[x]]

Inverse[Isometry[x_]] ^:=
  Isometry[MapAt[Minus,Transpose[MapAt[Minus,x,{{-1}}]],{{-1}}]]

Isometry[x_?MatrixQ] . Isometry[y_?MatrixQ] ^:= Isometry[x.y]
  
Isometry[M_?MatrixQ] . Minkowski[x_?VectorQ] ^:= Minkowski[M.x]

Isometry[M_?MatrixQ] . Minkowski[Vector[x_?VectorQ,v_?VectorQ]] ^:=
  Minkowski[Vector[M.x,M.v]]

Isometry[M_?MatrixQ] . model_?ModelQ[x_] ^:=
  Convert[Isometry[M].Convert[model[x],Minkowski],model] /;
    !SameQ[model,Minkowski]

End[]
EndPackage[]
