Linked cyclides circle by circle with Julia

LinkedCyclidesCbyC.jl

using Meshes
using Rotations
import TransformsBase as TB
using MeshViz
using GLMakie

function HopfFiber(q, t)
  [
    q[1]*cos(t) + q[2]*sin(t),
    sin(t)*(1 + q[3]),
    cos(t)*(1 + q[3]),
    q[1]*sin(t) - q[2]*cos(t)
  ] / sqrt(2*(1 + q[3]))
end

# rotation in 4d (right-isoclinic); xi is the angle of rotation
function rotate4d(alpha, beta, xi, vec)
  a = cos(xi)
  b = sin(alpha) * cos(beta) * sin(xi)
  c = sin(alpha) * sin(beta) * sin(xi)
  d = cos(alpha) * sin(xi)
  x, y, z, w = vec
  [
    a*x - b*y - c*z - d*w,
    a*y + b*x + c*w - d*z,
    a*z - b*w + c*x + d*y,
    a*w + b*z - c*y + d*x
  ]
end
  
# stereographic projection
function stereog(v)
  v[1:3] / (1 - v[4])
end

#
nCirclesByCyclide = 100
theta_ = LinRange(0, 2*pi, nCirclesByCyclide+1)[1:nCirclesByCyclide]
phi = 1 # -pi/2 < phi < pi/2; close to pi/2 <=> big hole

function Circle(theta)
  circle4d3pts = [
    rotate4d(
      pi/2, 0, 1, 
      HopfFiber([cos(theta)*cos(phi), sin(theta)*cos(phi), sin(phi)], t)
    )
    for t in (0, 2, 4)
  ]
  p1, p2, p3 = [stereog(p) for p in circle4d3pts]
  torus = Torus(
    Meshes.Point(tuple(p1...)), 
    Meshes.Point(tuple(p2...)), 
    Meshes.Point(tuple(p3...)), 
    0.2
  )
  discretize(torus)
end

circles = [
  Circle(theta) for theta in theta_
]

mesh1 = reduce(merge, circles)
trsfrm = Rotate(AngleAxis(2*pi/3, 0.0, 1.0, 0.0))
mesh2, _ = TB.apply(trsfrm, mesh1)
mesh3, _ = TB.apply(trsfrm, mesh2)

function drawLC()
  fig, ax, plt = MeshViz.viz(mesh1; color = :red)
  ax.show_axis = false
  MeshViz.viz!(mesh2; color = :green)
  MeshViz.viz!(mesh3; color = :blue)
  Makie.scale!(ax.scene, 1.6, 1.6, 1.6)
  Makie.rotate!(fig.scene, Meshes.Vec3f(1, 0, 0), pi/4)
  Makie.rotate!(ax.scene, Meshes.Vec3f(0, 0, 1), -pi/4)
  Makie.save("LinkedCyclides.png", fig)
end