(* PathG[p,dist]:
	Returns a list of {u_i} as below
	representing a geodesic path starting at
	p+I and moving towards p distance dist.
*)

(* Speed controls step size in hyperbolic length *)
Speed=0.2

(* Basic object is u={z,p,g}:
	z in H
	p in R-infinity
	g in GL_2(R)
Representing the geodesic from z to p; the image
is g(F(z)).  *)

(* Put z in I *)
Fund[{z_,p_,g_}]:= Which[
	Re[z]>1, 
		Fund[{IRef1[z],IRef1[p],g.MRef1}],
	Re[z]<0, 
		Fund[{IRef2[z],IRef2[p],g.MRef2}],
	Abs[z-1/2]<1/2,
		Fund[{IRef0[z],IRef0[p],g.MRef0}],
	True,    {z,p,g}
	]

(* Advance z along geodesic. *)
Flow[{z_,p_,g_}]:= Block[{d,w},
	w = (z-p)^-1;
	w = w + (Exp[Speed]-1) I Im[w];
	travel += Speed;
	Fund[{w^-1+p,p,g}]
	]

(* Develop path *)
PathG[p_,dist_]:=Block[{prec,travel,u,us},
	prec=20+dist/Log[10]//N;
	u = {N[p + I], 
	     N[p,prec], 
	     {{1,0,0},{0,1,0},{0,0,1}}
	};
	u = Fund[u];
	travel=0;
	us={u};
	While[travel<dist,
		u=Flow[u];
		us=Append[us,u]
		];
	us]

(* Utilities for rendering path and tiles *)
(* -------------------------------------- *)
Dev[{z_,p_,g_}]:= Aff[g,F[z]]
Tile[g_]:={Aff[g,1],
	   Aff[g,Omega],
	   Aff[g,Omega^2],
	   Aff[g,1]}

GPath[path_]:=Line[Map[Dev,path]]
GTile[gs_]:=Map[Line[Tile[#]]&,gs]

(* Horocycle paths *)
(* --------------- *)
(* PathH[p,dist]:
	Returns a list of {u_i} as below
	representing a horocycle path starting at
	p+I and moving towards p distance dist.
*)

(* Speed controls step size in hyperbolic length *)
SpeedH=0.1

(* Advance z along a horocycle. *)
(* Use Det[g] to keep orientation right *)
FlowH[{z_,p_,g_}]:= Block[{d,w},
	w = (p-z)^-1;
	w = w + Det[g]*(Exp[SpeedH]-1) Im[w];
	travel += SpeedH;
	Fund[{p-1/w,p,g}]
	]

(* Develop path *)
PathH[p_,dist_] := Path[p,dist,1]

PathH[p_,dist_,r_]:=Block[{prec,travel,u,us},
	prec=20+dist/Log[10]//N;
	u = {N[p + r*I], 
	     N[p,prec], 
	     {{1,0,0},{0,1,0},{0,0,1}}
	};
	u = Fund[u];
	travel=0;
	us={u};
	While[travel<dist,
		u=FlowH[u];
		us=Append[us,u]
		];
	us]
