Several Tilings of the Poincare Disk

In Euclidean geometry, there are only three ways to tile the plane using congruent regular polygons so that edges match up (using equilateral triangles, squares, or regular hexagons). There is much more freedom in hyperbolic geometry and a hyperbolic plane may be tiled in infinitely many ways. Here we illustrate this by showing two distinct tilings of the Poincare disk with equilateral triangles and two distinct tilings by regular, right angled polygons.

In the Poincare disk model of hyperbolic geometry, points are defined to be points within the interior of some fixed circle and lines are defined to be open diameters of that circle or circular arcs perpendicular to the fixed circle. Here are the tilings (the function PTiling is defined below. 

trianglePics = {PTiling[3,4,3], PTiling[3,5,3]};
rightAngledPics = {PTiling[5, 2 5 - 4], PTiling[6, 2 6 - 4]};
Show[GraphicsArray[{
	{Show[GraphicsArray[trianglePics],
		PlotLabel -> "Distinct Tilings by Equilateral Triangles"]},
	{Show[GraphicsArray[rightAngledPics],
		PlotLabel -> "Distinct Tilings by Right Angled Polygons"]}}]]
[Groovy Tilings of the Poincare Disk]

These images were generated with the following Mathematica code. 

(* *************************************************

   PTiling[n,k]  generates a groovy hyperbolic 
   picture when 2 < n < k.
   
   PTiling[3,k]  generates a tiling of the Poincare
   disk consisting of equilateral triangles with 
   angle (n-2) Pi/k.
   
   PTiling[n, 2n-4]    generates a tiling of the Poincare
   disk consisting of regular, right-angled n-gons,
   when n >= 5.
   
   The optional third argument, depth, should probably
   be left at 2.  3 sometimes works.

*************************************************** *)


PTiling[n_Integer, k_Integer, depth_:2] := Module[
	{aH, a, q, r, idealPoints, z1Ideal, z2Ideal,
	init, initCircs, PGamma,
	
	(* Auxilliary Function Declarations *)
	IdealPLine, RotaterMatrix, Rotate, Invert, 
	Iter, Reflect, arc, toR2, toC, mag},
	
	(* Auxilliary Function Definitions *)
	
	IdealPLine[{z1Given_, z2Given_}] := Module[
	{center, z1, z2},
	
		z1 = z1Given // N;  z2 = z2Given // N;
		center = Exp[I (Arg[z2] + Arg[z1])/2] / 
			Cos[(Arg[z2] - Arg[z1])/2] // N;
		arc[{z1,z2}, center]
	];
	
	arc[{z1_, z2_}, z0_] := Module[{theta1, theta2},
		theta1 = Arg[z1 - z0]; 
		theta2 = Arg[z2 - z0];
	
		If[Abs[N[theta1 - theta2]] > N[Pi], 
			If[N[theta1] < N[theta2], 
				theta1 = N[theta1 + 2Pi], 
				theta2 = N[theta2 + 2Pi]]];
	
		Circle[toR2[z0], Abs[z1-z0],
			Sort[{theta1, theta2}, N[#1] < N[#2] &]
		] // N
	];
	
	Invert[Circle[c_, r1_, ___], z_] := 
		r1^2/Conjugate[z-toC[c]] + toC[c] // N;
	Reflect[circ1_Circle, Circle[c2_, r2_, thetaList_]] := 
		IdealPLine[Invert[circ1, toC[c2 + r2 
			{Cos[#],Sin[#]}]]& /@ thetaList] // N;

	Iter[PLineList_List] := Union[PLineList, Select[Flatten[
		Outer[Reflect, PLineList, PLineList]], (# =!= Null &)]];

	RotaterMatrix[theta_] := 
		{{Cos[theta], -Sin[theta]}, {Sin[theta], Cos[theta]}};
	Rotate[Circle[c_, r_, thetaList_], theta_] := 
		Circle[RotaterMatrix[theta].c, r, theta + thetaList]  // N;
		
	toR2 = N[{Re[#],Im[#]}]&;
	toC  = N[#[[1]] + #[[2]] I] &;
	mag = N[Sqrt[#.#]]&;

	(* Basic Definitions *)
	aH = ArcCosh[(Cos[Pi/n] Cos[Pi/2] + 
		Cos[(n-2)Pi/(2k)])/(Sin[Pi/n] Sin[Pi/2])];
	a = (Exp[aH] - 1)/(Exp[aH] + 1) // N;
	q = (a + 1/a)/2 // N;
	r = q-a // N;
	PGamma = {Thickness[.008], GrayLevel[.3],
		Circle[{0,0},1]};
	
	(* The Action *)
	idealPoints = {x,y} /. NSolve[{x^2 + y^2 == 1, 
		(x-q)^2 + y^2 == r^2}, {x,y}]; 
	{z1Ideal, z2Ideal} = toC /@ idealPoints;
	
	init = IdealPLine[{z1Ideal, z2Ideal}];
	initCircs = NestList[Rotate[#, 2Pi/n]&, init, n-1];
	
	Show[Graphics[{Nest[Iter, initCircs, depth], PGamma}],
		AspectRatio -> Automatic]];