(* This package computes the function
	F : I -> T
   between the ideal triangle I spanning 0,1,inf in the
   upper half-plane, and the Euclidean triangle
   T spanning the cube roots of unity.

   F : (0,1,inf) -> (Omega,Omega^2,1)
*)
(* These parameters set the accuracy of computations:
	NumTermsT for Tri : D -> T
	NumTermsL for Lambda : I -> D
   Lev controls height at which F is rounded
   to a cusp.
*)
NumTermsT = 15
NumTermsL = 5
Lev = 10.0

(* H = upper half-plane
   I = ideal triangle in H
   T = equilaterial triangle
   D = unit disk
   L = lattice with fundamental domain T
*)

(* Tri : D -> T *)
(*----------------------------------------------------*)
(* Vertices of T are fixed by Tri *)

Omega=(-1/2)+(I Sqrt[3]/2)//N

(* h1 : D -> {z in D | Re[z]>0} *)
m1[z_]:=-Omega*(z-Omega^2)/(z-Omega)
m2[z_]:=-I*(1-2/(z+1))
h1[z_]:=m2[Sqrt[m1[z]]]

(* h2 : D -> Hexagon *)
series = Series[(1+z^6)^(-1/3),{z,0,NumTermsT}]
series = Normal[series]
series = N[Integrate[series,z]]
series = (1.5 series)/(series/.z->1)
h2[z_]=series

(* h3 : D -> T if z is near 1 *)
h3[z_]:=1-h2[h1[z]^(1/3)]

Tri[z_]:=Which[
	Abs[z-1]<1.01,h3[z],
	Abs[z-Omega]<1.01,Omega*h3[z/Omega],
	Abs[z-Omega^2]<1.01,h3[Omega*z]/Omega
	]

(*----------------------------------------------------*)
(* Lambda : I -> H *)
(* Lambda : 0,1,Inf -> 1,Inf,0 *)
(* Ahlfors p. 279 *)

LambdaFn[t_]:=Block[{e1,e2,e3},
	e1 = 1 + 2 Sum[Cos[Pi n t]^(-2),
		{n,1,NumTermsL}];
	e2 = 2 Sum[Sin[Pi (n+1/2) t]^(-2),
		{n,0,NumTermsL}];
	e3 = 2 Sum[Cos[Pi (n+1/2) t]^(-2),
		{n,0,NumTermsL}];
	Return[(e3-e2)/(e1-e2)]
]

(* Use modular group to get in domain
   of rapid convergence *)

Lambda[t_]:= Which[
	Re[t]> 1/2, 1/(1-1/Lambda[t-1]),
	Re[t]<-1/2, 1/(1-1/Lambda[t+1]),
	Abs[t]<1,  1-Lambda[-1/t],
	True, LambdaFn[t]
	]

(*----------------------------------------------------*)
(* U : H -> D *)
(* 0,1,inf -> 1,Omega,Omega^2 *)
U[z_]:=Omega^2*(z+Omega^2)/(z+Omega)

(*----------------------------------------------------*)
(* F : I -> T   *)
(*   : (0,1,inf) -> (Omega,Omega^2,1) *)
F[z_]:= Which[
	Im[z//N]>Lev,
		1.0,
	Im[-1/z//N]>Lev,
		Omega//N,
	Im[1/(1-z)//N]>Lev,
		Omega^2//N,
	True,
		Tri[N[U[N[Lambda[N[z]]]]]]
	]
