16
$\begingroup$

Creating a kaleidoscope in Mathematica is not a new topic at all, examples can be found from the links like Wolfram reference and Wolfram demonstrations.

My question is how to create Poincaré disk type of kaleidoscopes in Mathematica like the guys do in Java. Any idea and suggestion would be gratefully appreciated.

enter image description here enter image description here enter image description here

$\endgroup$
24
$\begingroup$

Hyperbolic lines in the Poincare disk are arcs of circles which, if extended, would meet the unit circle boundary orthogonally. Tiling the disk requires reflection of a tile's vertices and circular edges in those same circular edges. This corresponds to inversion in a circle.

Begin with a central polygon having p edges, with q such polygons meeting at each vertex. Then reflect this polygon in each of its circular edges to form p new polygons. Nest this iteration to the depth required.

Here is some code (perhaps too much). With these basics you can add variations like colouring hyperbolic polygons (not straight-edged polygons), or marking tiles with lines at the midpoints (forming the dual tiling), or extending to 3D, or truncating the tile corners to form additional shapes.

The function HyperbolicLine draws hyperbolic lines between points P and Q in the Poincare disk.

HyperbolicLine[{{Px_, Py_}, {Qx_, Qy_}}] :=
   If[N[Chop[Px Qy - Py Qx]] =!= 0.,
      Circle[OrthoCentre[{{Px,Py},{Qx,Qy}}],OrthoRadius[{{Px, Py},{Qx, Qy}}],
             OrthoAngles[{{Px, Py}, {Qx, Qy}}]],
      Line[{{Px, Py}, {Qx, Qy}}]]

The three functions within are as follows.

OrthoCentre[{{Px_, Py_}, {Qx_, Qy_}}] := 
   With[{d = 2 Px Qy - 2 Py Qx, p = 1 + Px^2, q = 1 + Qx^2 + Qy^2}, 
      If[N[d] =!= 0.,
         {p Qy + Py^2 Qy - Py q, -p Qx - Py^2 Qx + Px q}/d, 
         ComplexInfinity]]

OrthoRadius[{{Px_, Py_}, {Qx_, Qy_}}] := 
   If[N[Chop[Px Qy - Py Qx]] =!= 0., 
      Sqrt[Total[OrthoCentre[{{Px, Py}, {Qx, Qy}}]^2] - 1],
      Infinity]

OrthoAngles[{{Px_, Py_}, {Qx_, Qy_}}] := 
   Block[{a, b, c = OrthoCentre[{{Px, Py}, {Qx, Qy}}]},
      If[(a = N[Apply[ArcTan, {Px, Py} - c]]) < 0., a = a + 2 π]; 
      If[(b = N[Apply[ArcTan, {Qx, Qy} - c]]) < 0., 
      b = b + 2 π]; {a, b} = Sort[{a, b}]; 
      If[b - a > π, {b, a + 2 π}, {a, b}]]

The Inversion function inverts points, lines, and polygons in a circle.

Inversion[Circle[{Cx_, Cy_}, r_], {Px_, Py_}] := {Cx, Cy} + 
   r^2 {Px - Cx, Py - Cy}/((Cx - Px)^2 + (Cy - Py)^2)
Inversion[Circle[{Cx_, Cy_}, r_, {a_, b_}], {Px_, Py_}] := {Cx, Cy} + 
   r^2 {Px - Cx, Py - Cy}/((Cx - Px)^2 + (Cy - Py)^2)

Inversion[Circle[{Cx_, Cy_}, r_, {a_, b_}], p_Line] := 
   Map[Inversion[Circle[{Cx, Cy}, r], #] &, p, {2}]

Inversion[Circle[{Cx_, Cy_}, r_, {a_, b_}], p_Polygon] := 
   Map[Inversion[Circle[{Cx, Cy}, r], #] &, p, {2}]

Inversion[Line[{{Px_, Py_}, {Qx_, Qy_}}], {Ux_, Uy_}] := 
   With[{u = Px - Qx, v = Qy - Py},
      {-Ux (v^2 - u^2) - 2 u v Uy, Uy (v^2 - u^2) - 2 u v Ux}/(u^2 + v^2)]
Inversion[Line[{{Px_, Py_}, {Qx_, Qy_}}], p_Polygon] :=
   Map[Inversion[Line[{{Px, Py}, {Qx, Qy}}], #] &, p, {2}]

Inversion[Circle[{Cx_, Cy_}, r_], c_List] := 
   Map[Inversion[Circle[{Cx, Cy}, r], #] &, c]

Invert points of a polygon or list of polygons in a hyperbolic line (which is a line or a circle).

PolygonInvert[p_Polygon] := 
   Map[Inversion[HyperbolicLine[#], p] &, 
       Partition[Join[p[[1]], {p[[1, 1]]}], 2, 1]]
PolygonInvert[p_List] := Flatten[Map[PolygonInvert[#] &, p]]

Join the points of a hyperbolic polygon with straight or hyperbolic lines.

LineRule = Polygon[x_] :> Line[Join[x, {x[[1]]}]];
HyperbolicLineRule = 
   Polygon[x_] :> 
      Map[HyperbolicLine, Partition[Join[x, {x[[1]]}], 2, 1]];

The function CentralPolygon finds the initial polygon to begin the iterations. The first point of the polygon (with p sides) is in the first quadrant at angle π/p with respect to the horizontal. Each vertex is surrounded by q such polygons. The polygon may be rotated counter-clockwise through ϕ radians by using the optional third argument.

CentralPolygon[p_Integer, q_Integer, ϕ_: 0] := 
   With[{r = (Cot[π/p] Cot[π/q] - 1)/Sqrt[Cot[π/p]^2 Cot[π/q]^2 - 1],
        θ = π Range[1, 2 p - 1, 2]/p}, 
        r Map[{{Cos[ϕ], -Sin[ϕ]}, {Sin[ϕ], Cos[ϕ]}}.# &,
             Transpose[{Cos[θ], Sin[θ]}]]]

Delete duplicate polygons, allowing a tolerance.

PolygonUnion[p_Polygon, tol_: 10.^-10] := p
PolygonUnion[p_List, tol_: 10.^-10] := 
   With[{q = p /. Polygon[x_] :> N[Polygon[Round[x, 10.^-10]]]}, 
        DeleteDuplicates[q]]

HyperbolicTessellation tiles polygons with p sides, q of them meeting at each vertex, nesting down to level k. The entire tessellation is rotated counter-clockwise by ϕ radians.

HyperbolicTessellation[p_Integer, q_Integer, ϕ_, k_Integer, t_: 10.^-10] := 
   Map[PolygonUnion[#, t] &, 
       NestList[PolygonInvert, Polygon[CentralPolygon[p,q,ϕ]],k][[{-2, -1}]]] /; k > 0

HyperbolicTessellation[p_Integer, q_Integer, ϕ_, k_Integer, t_: 10.^-10] := 
   Polygon[CentralPolygon[p, q, ϕ]] /; k == 0

And finally,

HyperbolicTessellationGraphics[p_Integer, q_Integer, ϕ_, k_Integer,
      rule_RuleDelayed, opts___] := 
   Graphics[{Circle[{0, 0}, 1], 
      HyperbolicTessellation[p, q, ϕ, k, 10.^-10] /. rule}, opts]

For example,

HyperbolicTessellationGraphics[7, 3, 0., 3, HyperbolicLineRule, PlotLabel -> "{7,3}"]

Heptagons, 3 at each vertex

Graphics[{EdgeForm[{Thin, Black}], 
   Flatten[HyperbolicTessellation[5,4,0.,4]] /. Polygon[x_] :>
      {ColorData["DarkRainbow", Norm[Mean[x]]], Polygon[x]}}]

Pentagons, 4 at each vertex

Block[{poly = Flatten[HyperbolicTessellation[7, 3, 0., 4]]}, 
   Graphics[{EdgeForm[{Thin, Black}], 
      poly /. Polygon[x_] :> Polygon[x, VertexColors -> 
       Map[ColorData["SouthwestColors", Norm[#]^1.2 + RandomReal[{0.,0.3}]]&,x]]}]]

7,3 vertex colouring

And an example of truncation of a heptagon tiling, left as an exercise...

7,3 truncated

$\endgroup$
  • 3
    $\begingroup$ This would make for a very nice Demonstration (as in demonstrations.wolfram.com). $\endgroup$ – Daniel Lichtblau Nov 9 '14 at 16:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.