drawing_klein_group_sv.py

"""
in tr_al s d=[2,0] n=1
in tr_bl s d=[2,0] 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.]])

# make complex trace a and b
tr_a = tr_al[0] + tr_al[1] * 1.0j
tr_b = tr_bl[0] + tr_bl[1] * 1.0j

# create tr_ab, za, a and b
tr_ab = ((tr_a * tr_b) + ((-1.)**r) * \
        cmath.sqrt(tr_a**2. * tr_b**2. - 4. * (tr_a**2. + tr_b**2.)))/2.0
# print("tr_ab:", tr_ab)
# z0 = (tr_ab - 2) * tr_b / (tr_b * tr_ab - 2 * tr_a + 2 * tr_ab * 1.0j)
# a = np.array([[tr_a/2.0, (tr_a*tr_ab - 2*tr_b + 4.0j)/(2.0*tr_ab - 4.0) /z0], \
#         [(tr_a*tr_ab - 2*tr_b - 4.0j)*z0/(2.0*tr_ab - 4.0), tr_a/2.0]])
# b = np.array([[(tr_b-2.0j)/2.0, tr_b/2.0], [tr_b/2.0, (tr_b+2.0j)/2.0]])
z0 = (tr_ab - 2.) * tr_b / (tr_b * tr_ab - 2. * tr_a + tr_ab * 2.0j)
a = np.array([[tr_a/2.0, (tr_a*tr_ab - 2.*tr_b + 4.0j)/(2.0*tr_ab + 4.0)/z0 ], \
        [(tr_a*tr_ab - 2.*tr_b - 4.0j)*z0/(2.0*tr_ab - 4.0), tr_a/2.0]])
b = np.array([[(tr_b-2.0j)/2.0, tr_b/2.0], [tr_b/2.0, (tr_b+2.0j)/2.0]])

print("z0:", z0)
print("a:", a)
print("b:", b)
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)