hyperbolic_order5square_generators.py

# START
from mpmath import mp
import numpy as np
import matplotlib.pyplot as plt
import cmath

# Set the desired precision
mp.dps = 50

# Define the angles of the fundamental domain
alpha = mp.pi / 2 # Internal angle of the hyperbolic square
beta = mp.pi / 2 # Internal angle of the hyperbolic square
gamma = 2 * mp.pi / 5 # Angle at the vertex where five squares meet
# Angles for testing triangles instead:
#alpha = mp.pi / 2
#beta = mp.pi / 3
#gamma = mp.pi / 7


# Compute the side lengths of the fundamental domain using the hyperbolic law of cosines
a = mp.acosh((mp.cos(alpha) + mp.cos(beta) * mp.cos(gamma)) / (mp.sin(beta) * mp.sin(gamma)))
b = mp.acosh((mp.cos(beta) + mp.cos(gamma) * mp.cos(alpha)) / (mp.sin(gamma) * mp.sin(alpha)))
c = mp.acosh((mp.cos(gamma) + mp.cos(alpha) * mp.cos(beta)) / (mp.sin(alpha) * mp.sin(beta)))

# Compute the Euclidean representation of the vertices in the Poincaré disk model
z1 = mp.mpc(0) + mp.mpc(0, 1) * mp.sinh(a / 2)
z2 = z1 + mp.exp(mp.mpc(0, 1) * gamma) * mp.tanh(b / 2)
z3 = z1 + mp.tanh(c / 2)

# Define the Möbius transformation given two fixed points in matrix form
def mobius_matrix(z1, z2):
    return mp.matrix([[z1 - z2 + 1, -z1], [z2, z1 * z2 - z2 + 1]])

# Compute the generators in matrix form
g1_matrix = mobius_matrix(z1, z2)  # Reflection across z1-z2
g2_matrix = mobius_matrix(z2, z3)  # Reflection across z2-z3
g3_matrix = mobius_matrix(z1, z3)  # Reflection across z1-z3

# Define a function to compute the Möbius transformation from a 2x2 matrix
def mobius_from_matrix(matrix, z):
    return (matrix[0, 0] * z + matrix[0, 1]) / (matrix[1, 0] * z + matrix[1, 1])

# Check if the generators satisfy the relations
def is_identity_matrix(matrix, tol=1e-10):
    for i in range(2):
        for j in range(2):
            if i == j:
                if not mp.almosteq(matrix[i, j], 1, tol):
                    return False
            else:
                if not mp.almosteq(matrix[i, j], 0, tol):
                    return False
    return True

# Define the inverse of a Möbius transformation
def mobius_inverse(matrix):
    return mp.matrix([[matrix[1, 1], -matrix[0, 1]], [-matrix[1, 0], matrix[0, 0]]])

# Check if the generators satisfy the relations
identity = mp.matrix([[1, 0], [0, 1]])
check_relation1 = is_identity_matrix(g1_matrix * g1_matrix)
check_relation2 = is_identity_matrix(g2_matrix * g2_matrix)
check_relation3 = is_identity_matrix(g3_matrix * g3_matrix)

g123 = g1_matrix * g2_matrix * g3_matrix
check_relation4 = is_identity_matrix(g123 * g123 * g123 * g123 * g123)

print("Relation 1:", check_relation1)
print("Relation 2:", check_relation2)
print("Relation 3:", check_relation3)
print("Relation 4:", check_relation4)


# Setup for drawing a test plot

def draw_line_segment(z1, z2, ax, color="black"):
    z1 = complex(z1)
    z2 = complex(z2)
    k = z1.conjugate() * z2
    center = (z2 - z1 * k) / (1 - k)
    radius = abs(center - z1)
    angle1 = cmath.phase(z2 - center)
    angle2 = cmath.phase(z1 - center)
    
    if angle1 < 0:
        angle1 += 2 * np.pi
    if angle2 < 0:
        angle2 += 2 * np.pi
    
    if abs(angle1 - angle2) < 1e-5:
        return  # Do not draw the arc if the angles are too close
    
    center_x, center_y = center.real, center.imag
    arc = plt.Circle((center_x, center_y), radius, color=color, fill=False, lw=1.5)
    ax.add_patch(arc)
    ax.set_xlim(-1.5, 1.5)
    ax.set_ylim(-1.5, 1.5)
    arc.set_clip_box(ax.bbox)


def generate_transformed_points(vertices, generator_matrices, depth=0, max_depth=4):
    if depth >= max_depth:
        return vertices
    
    transformed_vertices = []
    for vertex in vertices:
        for generator_matrix in generator_matrices:
            transformed_vertex = mobius_from_matrix(generator_matrix, vertex)
            transformed_vertices.append(transformed_vertex)
    
    return vertices + generate_transformed_points(transformed_vertices, generator_matrices, depth + 1, max_depth)

vertices = generate_transformed_points([z1, z2, z3], [g1_matrix, g2_matrix, g3_matrix], max_depth=6)

fig, ax = plt.subplots(figsize=(10, 10))
ax.axis("equal")

# Draw the boundary of the Poincaré disk
disk_boundary = plt.Circle((0, 0), 1, color="black", fill=False)
ax.add_patch(disk_boundary)

# Draw the edges of the fundamental domain
draw_line_segment(z1, z2, ax, color="red")
draw_line_segment(z2, z3, ax, color="red")
draw_line_segment(z3, z1, ax, color="red")

# Draw the tiling
for i in range(len(vertices)):
    for j in range(i + 1, len(vertices)):
        if np.abs(np.abs(vertices[i] - vertices[j]) - 1) < 1e-6:
            draw_line_segment(vertices[i], vertices[j], ax)

plt.show()
# END