Creating Calabi Yau Manifold in python - Kukan Kogei -空間工芸-

f:id:asahidari:20200608145413p:plain

Overview

I got interested in manifold (in math)¹ after creating some minimal surfaces and implementing 4d rotation in my previous post. especially, Calabi-Yau manifold² looks beautiful. After finding some websites to modeling this, I tried to make this 3d model.

The way to create Calabi Yau manifold

I found some useful websites³ to make the surface with javascript. In this entry, I'll make it with python and Blender. (Actually I do not fully understand the meanings of this equations, so I just implemented this following the provided information.)

Mathematical equations describing the manifold are here. This is merely the outline. So please visit the link pages to know more details. Visualization of this surface can be done with projection of this equation onto 3d space.

$z_1^n+z_2^n=1$

Parameter z1 and z2 are expressed in the following way.

$\begin{eqnarray} \left\{ \begin{array}{1} z_1=\mathrm{e}^{\mathrm{i}\phi_1}[cos(x+\mathrm{i}y)]^\left(2/n\right) \\ z_2=\mathrm{e}^{\mathrm{i}\phi_2}[sin(x+\mathrm{i}y)]^\left(2/n\right) \end{array} \right. \end{eqnarray}$

$\phi_1=\frac{2\pi\mathrm{k}_1}{\mathrm{n}}\qquad (0\leq \mathrm{k}_1 < \mathrm{n}) \qquad \phi_2=\frac{2\pi\mathrm{k}_2}{\mathrm{n}}\qquad (0\leq \mathrm{k}_2 < \mathrm{n})$

Parameter k1 and k2 individually take Integer values from 0 to n - 1, and results in n x n parts of the manifold. And parameter x and y take the values in this range.

$0 \le x < \frac{\pi}{2}\qquad -\frac{\pi}{2} \le y < \frac{\pi}{2}$

Coordinates in 3d space are described here.

$(\mathrm{Re}(z_1), \mathrm{Re}(z_2), \mathrm{Im}(z_1)\cos(a)+\mathrm{Im}(z_2)\sin(a))$

This Fourth value is ignored.

$\mathrm{Im}(z_1)\sin(a)-\mathrm{Im}(z_2)\cos(a))$

Implementing with numpy

Here is the code to make the manifold. To use complex numbers, You have to import cmath library.

import numpy as np
import math, cmath

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.ticker as plticker
from matplotlib import cm

from scipy.spatial import Delaunay
import surf2stl_python.surf2stl as s2s

N = 9
a = 0.4
row, col = 30, 30
writeSTL = False

def calcZ1(x, y, k, n):
    return cmath.exp(1j*(2*cmath.pi*k/n)) * (cmath.cos(x+y*1j)**(2/n))

def calcZ2(x, y, k, n):
    return cmath.exp(1j*(2*cmath.pi*k/n)) * (cmath.sin(x+y*1j)**(2/n))

def calcZ1Real(x, y, k, n):
    return (calcZ1(x, y, k, n)).real

def calcZ2Real(x, y, k, n):
    return (calcZ2(x, y, k, n)).real

def calcZ(x, y, k1_, k2_, n, a_):
    z1 = calcZ1(x, y, k1, n)
    z2 = calcZ2(x, y, k2, n)
    return z1.imag * math.cos(a_) + z2.imag*math.sin(a_)

# set param range
x = np.linspace(0, math.pi/2, col)
y = np.linspace(-math.pi/2, math.pi/2, row)
x, y = np.meshgrid(x, y)

# init graph
fig = plt.figure(figsize=(18,8))

for n in range(2, N):
    ax = fig.add_subplot(2, 4, n - 1, projection='3d')
    ax.view_init(elev=15, azim=15)
    ax.set_title("n=%d" % n)
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    loc = plticker.MultipleLocator(base=1.0) # this locator puts ticks at regular intervals
    ax.xaxis.set_major_locator(loc)
    ax.yaxis.set_major_locator(loc)
    ax.zaxis.set_major_locator(loc)

    count = 0
    for k1 in range(n):
        for k2 in range(n):
            # calc X, Y, Z values
            X = np.frompyfunc(calcZ1Real, 4, 1)(x, y, k1, n).astype('float32')
            Y = np.frompyfunc(calcZ2Real, 4, 1)(x, y, k2, n).astype('float32')
            Z = np.frompyfunc(calcZ, 6, 1)(x, y, k1, k2, n, a).astype('float32')

            ax.plot_surface(X, Y, Z, cmap=cm.ocean, linewidth=0)

            # write to a STL file
            if writeSTL:
                X_ = X.flatten()
                Y_ = Y.flatten()
                Z_ = Z.flatten()
                delaunay_tri = Delaunay(np.array([x_, y_]).T)
                s2s.tri_write('output_calabi-yau_n%d_%d.stl'% (n, count), X_, Y_, Z_, delaunay_tri)

            count += 1

If successfully done, these graphs are displayed. f:id:asahidari:20200608180626p:plain This code is uploaded in my gist.⁴

Implementing with sympy

This code is sympy version. I found sympy's symbols can not be reused in loop processing, so in this code, symbols are generated in each loop.

import sympy as sp

import numpy as np
import math
from sympy.utilities.lambdify import lambdify

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

from scipy.spatial import Delaunay
import surf2stl_python.surf2stl as s2s

# function to bind sympy plotting to backend figure and ax
# https://stackoverflow.com/questions/60325325/putting-together-plots-of-matplotlib-and-sympy
def move_sympyplot_to_axes(p, ax):
    backend = p.backend(p)
    backend.ax = ax
    # Fix for > sympy v1.5
    backend._process_series(backend.parent._series, ax, backend.parent)
    backend.ax.spines['right'].set_color('none')
    backend.ax.spines['bottom'].set_position('zero')
    backend.ax.spines['top'].set_color('none')
    plt.close(backend.fig)

    
N = 9
A = 0.4
writeSTL = False

# setup plotting boxes
fig = plt.figure(figsize=(16,8))
axes = []
for g in range(7):
    axes.append(fig.add_subplot(2,4,g+1, projection='3d'))
    axes[g].view_init(elev=30, azim=45)
    axes[g].set_title("n=%d" % (g+2))
    axes[g].set_xlabel('X')
    axes[g].set_ylabel('Y')


a = sp.symbols('a')
for index in range(2, N):
    
    n = sp.symbols('n_%d'%index)
    p_array = []
    count = 0
    for i in range(index):
        k1 = sp.symbols('k1_%d_%d'%(index, i))
        for j in range(index):
            k2 = sp.symbols('k2_%d_%d'%(index, j))

            # declare symbols
            z1 = sp.symbols('z1_%d_%d'%(index, count))
            z2 = sp.symbols('z2_%d_%d'%(index, count))
            x = sp.symbols('x_%d_%d'%(index, count))
            y = sp.symbols('y_%d_%d'%(index, count))
            X = sp.symbols('X_%d_%d'%(index, count))
            Y = sp.symbols('Y_%d_%d'%(index, count))
            Z = sp.symbols('Z_%d_%d'%(index, count))

            # prepare z1, z2 equations
            z1 = sp.exp(sp.I*((2*sp.pi*k1)/n)) * (sp.cos(x+y*sp.I))**(2/n)
            z2 = sp.exp(sp.I*((2*sp.pi*k2)/n)) * (sp.sin(x+y*sp.I))**(2/n)
            z1 = z1.subs([(n, index), (k1, i)])
            z2 = z2.subs([(n, index), (k2, j)])

            # set equations to symbols
            X = sp.re(z1)
            Y = sp.re(z2)
            Z = sp.im(z1)*sp.cos(a) + sp.im(z2)*sp.sin(a)
            Z = Z.subs(a, A)

            # draw each parts
            p = sp.plotting.plot3d_parametric_surface(X, Y, Z, (x, 0, sp.pi/2), (y, -sp.pi/2, sp.pi/2), show=False)            
            p_array.append(p)

            # write a part into a STL file
            if writeSTL:
                x0 = np.linspace(0, math.pi/2, 30)
                y0 = np.linspace(-math.pi/2, math.pi/2, 30)
                x_, y_ = np.meshgrid(x0, y0)
                x_, y_ = x_.flatten(), y_.flatten()

                X_ = lambdify((x, y), X, "numpy")
                Y_ = lambdify((x, y), Y, "numpy")
                Z_ = lambdify((x, y), Z, "numpy")

                X_npf = np.frompyfunc(X_, 2, 1)
                Y_npf = np.frompyfunc(Y_, 2, 1)
                Z_npf = np.frompyfunc(Z_, 2, 1)

                X_np = X_npf(x_, y_).astype(float)
                Y_np = Y_npf(x_, y_).astype(float)
                Z_np = Z_npf(x_, y_).astype(float)

                delaunay_tri = Delaunay(np.array([x_, y_]).T)
                s2s.tri_write('output_calabi-yau_n%d_%d.stl'%(index, count), X_np, Y_np, Z_np, delaunay_tri)

            # inclement symbol number
            count += 1

    # show all parts
    for k in range(1, index*index):
        p_array[0].append(p_array[k][0])

    move_sympyplot_to_axes(p_array[0], axes[index-2])
    # print(p_array[0])
    # p_array[0].show()
    
plt.show()

These graphs will be displayed. f:id:asahidari:20200608181717p:plain

This is also uploaded in my gist.⁵

Exporting the manifold parts into STL files

In both ways (using numpy and sympy), these codes output each parts of the surface into STL files. Surf2stl_python script is a light weight script, and you can get the code from my github⁶.

Importing STL files and Combining the parts in Blender

After exporting STL files, You need to combine all of the parts into one object. You can use various 3d modeling software, but In this document, I'll use Blender.

To import the STL files, select menu->files->import->stl, and import all the parts. After loading the files, join them together with Ctrl+J. After successfully united, in edit mode, mesh->clean up->merge by distances, and you may also have to fill some holes. f:id:asahidari:20200608180024p:plain Done! Thanks for reading. f:id:asahidari:20200608180053p:plain