Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Hopf tori with Python
import os
from math import cos, sin, sqrt, pi
import numpy as np
import pyvista as pv
from planegeometry.geometry import Triangle
A = 0.44
n = 3
def Gamma(t):
alpha = np.pi/2 - (np.pi/2-A)*np.cos(n*t)
beta = t + A*np.sin(2*n*t)
return np.array([
np.sin(alpha) * np.cos(beta),
np.sin(alpha) * np.sin(beta),
np.cos(alpha)
])
def HopfInverse(p, phi):
return np.array([
(1+p[2])*np.cos(phi),
p[0]*np.sin(phi) - p[1]*np.cos(phi),
p[0]*np.cos(phi) + p[1]*np.sin(phi),
(1+p[2])*np.sin(phi)
]) / np.sqrt(2*(1+p[2]))
def Stereo(q):
return 2*q[0:3] / (1-q[3])
def F(t, phi):
return Stereo(HopfInverse(Gamma(t), phi))
def HTmesh(nu=400, nv=200):
angle_u = np.linspace(-np.pi, np.pi, nu)
angle_v = np.linspace(0, np.pi, nv)
u, v = np.meshgrid(angle_u, angle_v)
z,x,y = F(u, v)
grid = pv.StructuredGrid(x, y, z)
mesh = grid.extract_geometry().clean(tolerance=1e-6)
return mesh
mesh = HTmesh()
# -----------------------------------------------------------------------------
def Hexlet(Center, Radius):
s = Radius # side length of the hexagon
Coef = 2/3
a = Coef*(Radius+s/2)/sin(pi/2-2*pi/6)
I = np.array([a, 0]) # inversion pole
## ------------------------------------------------------------------ ////
O1 = np.array([2*a, 0, 0])
# interior sphere
def inversion(M, RRadius):
II = np.array([Coef*(RRadius+RRadius/2)/sin(pi/2-2*pi/6), 0])
S = Coef*(RRadius+RRadius/2) * np.array([cos(2*pi/6),sin(2*pi/6)])
k = np.vdot(S-II,S-II) # negated inversion constant
M = np.asarray(M, dtype=float)
IM = M-II
return II - k/np.vdot(IM,IM)*IM
SmallRadius = Coef*(Radius-s/2)
p1 = inversion((SmallRadius,0), Radius)
p2 = inversion((0,SmallRadius), Radius)
p3 = inversion((-SmallRadius,0), Radius)
tr = Triangle(p1, p2, p3)
cs = tr.circumcircle()
shift = pi/90
for frame_number in range(1, 181):
pltr = pv.Plotter(window_size=[512,512], off_screen=True)
pltr.set_background("#363940")
i = 1
while i<= 6:
beta = i*pi/3 - frame_number*shift; # frame from 1 to 180
ccenter = Coef*Radius*np.array([cos(beta),sin(beta)])
p1 = inversion((0,Coef*Radius/2)+ccenter,Radius)
p2 = inversion((Coef*Radius/2,0)+ccenter,Radius)
p3 = inversion((0,-Coef*Radius/2)+ccenter,Radius)
tr = Triangle(p1, p2, p3)
cs = tr.circumcircle()
center = np.array([cs.center[0], cs.center[1], 0])
r = cs.radius
mesh = HTmesh().copy()
mesh.scale(r*0.05)
mesh.rotate_z(2*frame_number)
mesh.translate(center-O1)
pltr.add_mesh(mesh, color="orangered", specular=20, smooth_shading=True)
i += 1
pltr.set_focus([0,0,0])
pltr.set_position([2, 0, 8])
pltr.show(screenshot="pic_htso%03d.png" % frame_number)
Hexlet((0,0,0),2)
os.system(
"magick convert -dispose previous -loop 0 -delay 8 pic_htso*.png HopfToriSteinerOrbit.gif"
)
View raw

(Sorry about that, but we can’t show files that are this big right now.)

View raw

(Sorry about that, but we can’t show files that are this big right now.)

View raw

(Sorry about that, but we can’t show files that are this big right now.)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment