球体を裏返す

例えば， このYouTubeのビデオ のように，球体を連続的に裏返す方法について話しているビデオがいくつかネット上にあります．Adam BednorzとWitold Bednorzによる 素晴らしい論文, Differential and its Applications 64, 59 (2019) ( arXiv にも掲載)のおかげで，このいわゆる球の反転(表面をつまんだり破ったりせずに，言い換えればトポロジーを保ったまま球を裏返すこと)をプロットすることができます．

数学が絡むとちょっと，いや，絡んでいるように見えるかもしれませんね．重要なのは，このアニメーションで見られる全体的なプロセスです．まず，球体を波型にして少し引き伸ばし，滑らかな変換のためのスペースを確保してから，ローブ同士をねじり，球体を展開するためにこのプロセスを逆にします．この変換が本当に滑らかであることは明らかではありませんが，これはBednorzとBednorzによる論文で証明されています．

import numpy as np

import pyvista as pv

# define some parameters
n_steps = 30  # number of steps for a given stage of the animation
Q = 2 / 3  # arbitrary < 1
w = 2  # arbitrary > 0
n = 2  # arbitrary integer > 1, the number of "lobes"
beta = 1  # arbitrary > 1
alpha_final = 1  # arbitrary > 1
eta_final = 2  # arbitrary > 1
kappa = (n - 1) / (2 * n)

外転の実装に必要なマッピングの連鎖を定義してみましょう:

def sphere_to_cylinder(theta, phi):
    """Map from spherical polar coordinates to cylindrical ones.

    Input theta goes from -pi/2 to pi/2, phi goes from -pi to pi.
    Output h goes from -infinity to infinity, phi goes from -phi to phi.

    """
    h = w * np.sin(theta) / np.cos(theta) ** n
    # phi is unchanged
    return h, phi


def cylinder_to_wormhole(h, phi, t, p, q):
    """Map from a cylinder to an open wormhole using Eq. (4).

    Input h goes from -infinity to infinity, phi goes from -phi to phi.
    Output is an (x, y, z) point embedded in 3d space.

    The parameters t, p, q vary during the eversion process.
    Start from |t| > 1 (fixed), p = 1 and q = 0. End at p = 0, qt = +-1.

    """
    x = t * np.cos(phi) + p * np.sin((n - 1) * phi) - h * np.sin(phi)
    y = t * np.sin(phi) + p * np.cos((n - 1) * phi) + h * np.cos(phi)
    z = h * np.sin(n * phi) - t / n * np.cos(n * phi) - q * t * h
    return x, y, z


def close_wormhole(x0, y0, z0, eta, xi, alpha):
    """Close the wormhole using Eqs. (7)-(8).

    Input is an (x0, y0, z0) point embedded in 3d space.
    Output is an (x2, y2, z2) == (x'', y'', z'') point embedded in 3d space.

    The parameters eta, xi, alpha vary during the eversion process.

    """
    # Eq. (7): (x, y, z) -> (x', y', z')
    denominator = xi + eta * (x**2 + y**2)
    x1 = x0 / (denominator**kappa)
    y1 = y0 / (denominator**kappa)
    z1 = z0 / denominator

    gamma = 2 * np.sqrt(alpha * beta)
    # singular case, Eq (9):
    if np.isclose(gamma, 0):
        denominator = x1**2 + y1**2
        x2 = x1 / denominator
        y2 = y1 / denominator
        z2 = -z1
        return x2, y2, z2

    # Eq. (8): (x', y', z') -> (x'', y'', z'')
    exponential = np.exp(gamma * z1)
    numerator = alpha - beta * (x1**2 + y1**2)
    denominator = alpha + beta * (x1**2 + y1**2)
    x2 = x1 * exponential / denominator
    y2 = y1 * exponential / denominator
    z2 = numerator / denominator * exponential / gamma - (alpha - beta) / (alpha + beta) / gamma
    return x2, y2, z2


def unfold_sphere(theta, phi, t, q, eta, lamda):
    """Unfold the sphere using Eqs. (12), (15), (10).

    Input is a (theta, phi) point in spherical coordinates.
    Output is an (x, y, z) point embedded in 3d space.

    The parameter lamda varies. Lamda = 1 is the final stage of the
    wormhole closing, and lamda = 0 is the recovered sphere.

    """
    # apply Eqs. (12), (15)
    # fmt: off
    x = (
        t * (1 - lamda + lamda * np.cos(theta)**n) * np.cos(phi)
        - lamda * w * np.sin(theta) * np.sin(phi)
    )
    x /= np.cos(theta)**n
    y = (
        t * (1 - lamda + lamda * np.cos(theta)**n) * np.sin(phi)
        + lamda * w * np.sin(theta) * np.cos(phi)
    )
    y /= np.cos(theta) ** n
    z = (
        lamda * (
            (w * np.sin(theta) * (np.sin(n * phi) - q * t)) / np.cos(theta)**n
            - t / n * np.cos(n * phi)
        )
        - (1 - lamda) * eta**(1 + kappa) * t * abs(t)**(2 * kappa)
            * np.sin(theta) / np.cos(theta)**(2 * n)  # noqa
    )
    # fmt: on

    # apply Eq. (10)
    denominator = x**2 + y**2
    x2 = x * eta**kappa / denominator ** (1 - kappa)
    y2 = y * eta**kappa / denominator ** (1 - kappa)
    z2 = -z / eta / denominator
    return x2, y2, z2

ここで，論文の表1のプロセスを実行することによって，関数を連鎖させます．下から順に t = -1/Q とし，行ごとに変化するパラメータを線形に変化させながらステップアップを続け，一番上で t = -1/Q から t = 1/Q へ，そして上から下へと戻っていきます．各フレームを GIF に保存します．

pyvista.Plotter.add_mesh() のキーワードパラメータ backface_params をうまく利用して，内側と外側で異なる色でプロットできるようにしました．

# plot options to use for each frame
opts = dict(
    color='aquamarine',
    specular=1.0,
    specular_power=50.0,
    backface_params=dict(color='forestgreen'),
    smooth_shading=True,
    reset_camera=True,
)

# use a small figure window to reduce the size of the GIF
plotter = pv.Plotter(window_size=(300, 300))
plotter.open_gif('sphere_eversion.gif')


def save_frame(x, y, z):
    """Helper to generate and store a frame of the eversion."""
    plotter.clear()
    plotter.enable_lightkit()
    plotter.add_mesh(pv.StructuredGrid(x, y, z), **opts)
    plotter.write_frame()


# initial parameters, will be updated
t = -1 / Q
q = Q
p = xi = alpha = 0
eta = 1

# sphere -> inverted wormhole
theta, phi = np.mgrid[-np.pi / 2 : np.pi / 2 : 200j, -np.pi : np.pi : 400j]
h, phi = sphere_to_cylinder(theta, phi)
for lamda in np.linspace(0, 1, n_steps, endpoint=False):
    x2, y2, z2 = unfold_sphere(theta, phi, t, q, eta, lamda)
    save_frame(x2, y2, z2)

# inverted wormhole -> unfolded wormhole
x, y, z = cylinder_to_wormhole(h, phi, t, p, q)
xis = np.linspace(0, 1, n_steps)
alphas = np.linspace(0, alpha_final, n_steps)
etas = np.linspace(1, eta_final, n_steps)
for xi, alpha, eta in zip(xis, alphas, etas):
    x2, y2, z2 = close_wormhole(x, y, z, eta, xi, alpha)
    save_frame(x2, y2, z2)

# unfolded wormhole -> closed wormhole
for q in np.linspace(Q, 0, n_steps):
    p = 1 - abs(q * t)
    x, y, z = cylinder_to_wormhole(h, phi, t, p, q)
    x2, y2, z2 = close_wormhole(x, y, z, eta, xi, alpha)
    save_frame(x2, y2, z2)

# closed wormhole turned inside out (flip sign of time)
# unfolded wormhole -> closed wormhole
for t in np.linspace(-1 / Q, 1 / Q, n_steps):
    p = 1 - abs(q * t)
    x, y, z = cylinder_to_wormhole(h, phi, t, p, q)
    x2, y2, z2 = close_wormhole(x, y, z, eta, xi, alpha)
    save_frame(x2, y2, z2)

# closed wormhole -> unfolded wormhole
for q in np.linspace(0, Q, n_steps + 1)[1:]:
    p = 1 - abs(q * t)
    x, y, z = cylinder_to_wormhole(h, phi, t, p, q)
    x2, y2, z2 = close_wormhole(x, y, z, eta, xi, alpha)
    save_frame(x2, y2, z2)

# unfolded wormhole -> inverted wormhole
x, y, z = cylinder_to_wormhole(h, phi, t, p, q)
xis = np.linspace(1, 0, n_steps + 1)[1:]
alphas = np.linspace(alpha_final, 0, n_steps + 1)[1:]
etas = np.linspace(eta_final, 1, n_steps + 1)[1:]
for xi, alpha in zip(xis, alphas):
    x2, y2, z2 = close_wormhole(x, y, z, eta, xi, alpha)
    save_frame(x2, y2, z2)

# inverted wormhole -> sphere
for lamda in np.linspace(1, 0, n_steps + 1)[1:]:
    x2, y2, z2 = unfold_sphere(theta, phi, t, q, eta, lamda)
    save_frame(x2, y2, z2)

plotter.close()

t = 0 の中間状態の静止画を見ると，球の内側と外側の2つのローブが見える，見事な対称配置になっていることがわかります．

t = q = 0
xi = p = 1
eta = eta_final
alpha = alpha_final

x, y, z = cylinder_to_wormhole(h, phi, t, p, q)
x2, y2, z2 = close_wormhole(x, y, z, eta, xi, alpha)

plotter = pv.Plotter(window_size=(512, 512))
plotter.add_mesh(pv.StructuredGrid(x2, y2, z2), **opts)
plotter.show()