drawing_klein_group_maskit_sv.py

"""
in u s d=[0,2] n=1
in r s d=1 n=2
in epsilon s d=0.005 n=2
in levmax s d=50 n=2
out verts_out v
out edges_out s
"""

import numpy as np
import cmath
import math

# mobius transformation for a point
def mobiusOnPoint(mat, z):
    a, b, c, d = mat[0][0], mat[0][1], \
                    mat[1][0], mat[1][1]
    if (c * z + d) == 0:
        return np.inf + 0.j
    return (a * z + b) / (c * z + d)

# Calc fix point (plus) for a matrix
def fixPointPlus(mat):
    if mat[1][0] == 0:
        return np.inf
    a, b, c, d = mat[0][0], mat[0][1], \
                    mat[1][0], mat[1][1]
    return ((a-d) + cmath.sqrt((a+d)**2 - 4.)) / (2.*c)

def go_forward(lev, tags, gens, word):
    lev = lev + 1
    tags[lev] = (tags[lev-1] + 1) % 4
    word[lev] = word[lev-1].dot(gens[tags[lev]])
    return lev
    
def go_backward(lev):
    return lev-1

def available_turn(lev, tags):
    return ((tags[lev+1]-1) % 4) != ((tags[lev]+2) % 4)

def turn_and_go_forward(lev, tags, gens, word):
    tags[lev+1] = ((tags[lev+1]-1) % 4)
    if lev < 0:
        word[0] = gens[tags[0]]
    else:
        word[lev+1] = word[lev].dot(gens[tags[lev+1]])
    return lev+1

def branch_termination(lev, levmax, tags, word, fp, endpt):
    # new_z = mobiusOnPoint(word[lev], endpt[tags[lev]])
    zz = np.array([])
    for i in range(0, 3):
        zz_c = mobiusOnPoint(word[lev], fp[tags[lev]][i])
        zz = np.append(zz, np.array([zz_c.real, zz_c.imag]))
    zz = zz.reshape(3, 2)
    
    # new_point = np.array([new_z.real, new_z.imag])
    # if np.linalg.norm(new_point - old_point) < epsilon \
    if (lev == levmax - 1) \
        or (math.fabs(np.linalg.norm(zz[1] - zz[0])) < epsilon \
        and math.fabs(np.linalg.norm(zz[2] - zz[1])) < epsilon):
        return True, zz
    return False, np.array([0, 0, 0])
    
    
# ======= Main program ======
verts = []
edges = []
mat_id = np.array([[1., 0.], [0., 1.]])

m = u[0] + u[1] * 1.0j
a = np.array([[-m * 1.0j, -1.0j], [-1.0j, 0]])
b = np.array([[1, 2], [0, 1]])

A = np.linalg.inv(a)
B = np.linalg.inv(b)

gens = np.array([a, b, A, B])

begpt, endpt = np.array([]), np.array([])
for k in range(0, 4):
    begpt = np.append(begpt, fixPointPlus( gens[(k+1)%4] * gens[(k+2)%4] * gens[(k+3)%4] * gens[k] ))
    endpt = np.append(endpt, fixPointPlus( gens[(k-1)%4] * gens[(k-2)%4] * gens[(k-3)%4] * gens[k] ))
    
# repet = np.array([[b*A*B*a, a, B*A*b*a], \
#                 [A*B*a*b, b, a*B*A*b], \
#                 [B*a*b*A, A, b*a*B*A], \
#                 [a*b*A*B, B, A*b*a*B]])
repet = np.array([[b.dot(A.dot(B.dot(a))), a, B.dot(A.dot(b.dot(a)))], \
                [A.dot(B.dot(a.dot(b))), b, a.dot(B.dot(A.dot(b)))], \
                [B.dot(a.dot(b.dot(A))), A, b.dot(a.dot(B.dot(A)))], \
                [a.dot(b.dot(A.dot(B))), B, A.dot(b.dot(a.dot(B)))]])
fp = np.array([[fixPointPlus(p) for p in points] for points in repet])
print("fp:", fp)

lev = 0
word = np.array([[1 + 0.j, 0+0.j, 0+0.j, 1+0.j]  \
                for i in range(0, levmax)]).reshape(levmax, 2, 2)
word[0] = gens[0]
tags = np.zeros(levmax, dtype=int)

# oldpoint = np.array([begpt[3].real, begpt[3].imag])
# verts_out.append([oldpoint[0], oldpoint[1], 0])

count = 10000

while not (lev < 0 and tags[0] == 1):
    
    count2 = 100
    while count2 > 0:
        result = branch_termination(lev, levmax, tags, word, fp, endpt)
        if result[0]:
            # print("New point detected!")
            # newpoint = result[1]
            # oldpoint = result[1]
            #verts_out.append([newpoint[0], newpoint[1], 0])
            for zz in result[1]:
                verts.append([zz[0], zz[1], 0])
            # if len(verts) > 3:
            #     edges.append([len(verts)-4, len(verts)-3])
            edges.append([len(verts)-3, len(verts)-2])
            edges.append([len(verts)-2, len(verts)-1])
            break
        
        if lev >= 0:
            lev = go_forward(lev, tags, gens, word)
            
        count2 = count2 - 1
    
    # print("count2:", count2)
    lev = go_backward(lev)
    count3 = 100
    while lev >= 0 and not available_turn(lev, tags):
        lev = go_backward(lev)
        count3 = count3 - 1
        if count3 < 0:
            break
    # print("count3:", count3)
    
    lev = turn_and_go_forward(lev, tags, gens, word)

    if lev == 0 and tags[0] == 0:
        print("lev == 0 and tags[0] == 0")
        break
    
    count = count - 1
    if count < 0:
        print("count < 0")
        break

verts_out = [verts]
edges_out = [edges]

print("len(verts)", len(verts))
print("count:", count)