kleinian_double_1_15_cusp_group_an.py

"""
in trace_a s d=[0,0] n=1
in trace_b s d=[0,0] n=1
in colorize s d=1 n=2
in reverse s d=0 n=2
out centers s
out radii s
out hue s
"""

# source: http://bugman123.com/Fractals/index.html Kleinian Double 1/15 Cusp Group Pov-Ray Code
import cmath
import numpy as np
# from mathutils import Vector
from animation_nodes.data_structures import Vector3DList, EdgeIndicesList


ta = trace_a[0] + trace_a[1]*1.0j
tb = trace_b[0] + trace_b[1]*1.0j
# ta = 1.958591030 - 0.011278560 * 1.0j
# tb = 2.0 + 0.0j
# tab = (ta*tb - cmath.sqrt(ta**2*tb**2 - 4*(ta**2 + tb**2))) / 2.0
tab = (ta*tb + cmath.sqrt((ta**2)*(tb**2) - 4*(ta**2 + tb**2))) / 2.0
z0 = (tab-2.0)*tb / (tb*tab - 2*ta + 2*tab*1.0j)

b = np.array([[(tb - 2.0j)/2.0, tb/2.0], [tb/2.0, (tb + 2.0j)/2.0]])
B = np.linalg.inv(b)
a = np.array([[tab, (tab-2.0)/z0], [(tab+2.0)*z0, tab]]).dot(B)
# a = np.array([[ta/2.0, (ta*tab - 2*tb + 4.0j)/(2*tab + 4)/z0], [(ta*tab - 2*tb - 4.0j)*z0/(2*tab - 4), ta/2.0]])
A = np.linalg.inv(a)

def Fix(X):
    return (X[0][0] - X[1][1] \
         - cmath.sqrt(4.0*X[0][1]*X[1][0] + (X[0][0] - X[1][1])**2)) / (2.0*X[1][0])

def ToMatrix(z, r):
    Z = np.array([[z, r**2 - z*z.conjugate()], [1.0, -z.conjugate()]])
    return (1.0j/r) * Z

def MotherCircle(M1, M2, M3):
    z1, z2, z3 = Fix(M1), Fix(M2), Fix(M3)
    x1, x2, x3 = z1.real, z2.real, z3.real
    y1, y2, y3 = z1.imag, z2.imag, z3.imag
    zx = (x3**2)*(y1 - y2) + ((x1**2) + (y1 - y2)*(y1-y3))*(y2-y3) + (x2**2)*(y3-y1)
    zy = -(x2**2)*x3 + (x1**2)*(x3 - x2) + x3*(y1 - y2)*(y1 + y2) \
        + x1*(x2**2 - x3**2 + y2**2 - y3**2) + x2*(x3**2-y1**2+y3**2) \
        / (2*(x3*(y1 - y2) + x1*(y2 - y3) + x2*(y3 - y1)))
    z = zx + (zy * 1.0j)
    return ToMatrix(z, abs(z1 - z))

# C1 = MotherCircle(b, np.dot(a, np.dot(b, A)), np.dot(a, np.dot(b, np.dot(A, B))))
# C2 = MotherCircle(np.dot(b, np.linalg.matrix_power(a,15)), np.dot(a, np.dot(b, np.linalg.matrix_power(a,14))), np.dot(a, np.dot(b, np.dot(A,B))));
C1 = MotherCircle(b, a.dot(b.dot(A)), a.dot(b.dot(A.dot(B))))
C2 = MotherCircle(b.dot(np.linalg.matrix_power(a,15)), a.dot(b.dot(np.linalg.matrix_power(a,14))), a.dot(b.dot(A.dot(B))));

def Reflect(C, M):
    return M.dot(C.dot(np.linalg.inv( \
            np.array([[M[0][0].conjugate(), M[0][1].conjugate()], \
                    [M[1][0].conjugate(), M[1][1].conjugate()]]))))

def zcen(Z):
    return (Z[0][0] / Z[1][0])

def rad(Z):
    return (1.0j/Z[1][0]).real

def spiral(C0, M, n, colorize_, reverse_):
    c = []
    r = []
    h = []
    C = C0
    i = 0
    while i < (n - 1):
        if colorize_:
            hue = (18 + (1 if reverse else -1)*i)/15 % 1
            h.append(hue)
        z = zcen(C)
        c.append([float(z.real), float(z.imag), 0.0])
        r.append(float(rad(C)))
        i = i + 1
        C = Reflect(C, M)
    return c, r, h

centers_np = np.array([])
radii_np = np.array([])
hue_np = np.array([])

c, r, h = spiral(C1, a, 83, colorize, reverse)
centers_np = np.append(centers_np, c)
radii_np = np.append(radii_np, r)
hue_np = np.append(hue_np, h)
c, r, h = spiral(C1, A, 83, colorize, reverse)
centers_np = np.append(centers_np, c)
radii_np = np.append(radii_np, r)
hue_np = np.append(hue_np, h)
c, r, h = spiral(C2, a, 91, colorize, reverse)
centers_np = np.append(centers_np, c)
radii_np = np.append(radii_np, r)
hue_np = np.append(hue_np, h)
c, r, h = spiral(C2, A, 76, colorize, reverse)
centers_np = np.append(centers_np, c)
radii_np = np.append(radii_np, r)
hue_np = np.append(hue_np, h)

centers_list = centers_np.reshape(-1,3).tolist()
centers = Vector3DList()
for center in centers_list:
    centers.append(Vector((center[0], center[1], center[2])))
    
radii = radii_np.tolist()
hue = hue_np.tolist()
# spiral(C1, a, 83, 1, 0)
# spiral(C1, A, 83, 1, 1)
# spiral(C2, a, 91, 0, 0)
# spiral(C2, A, 76, 0, 1)