Tesseract with Julia

tesseract.jl

import Meshes
using MeshViz
using GLMakie
using Makie
using Printf

# vectors obtained by changing signs of the given vectors
function changesOfSigns(vecs)
  N = length(vecs[1])
  signs = Iterators.product(ntuple(i -> (-1, 1), N)...)
  return [[sign[i] * vec[i] for i in 1:N]
          for vec in vecs for sign in signs
          if all(vec[i] != 0 || sign[i] == 1 for i in eachindex(vec))]
end

# tesseract vertices
vertices = changesOfSigns([[1, 1, 1, 1]])

# tesseract edges
edges = Vector{Tuple{Int,Int}}(undef, 0)
for i in 1:(length(vertices)-1)
  for j in (i+1):length(vertices)
    if count(==(0), vertices[i] - vertices[j]) == 3
      push!(edges, (i, j))
    end
  end
end

# rotation in 4D (right-isoclinic)
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)
  sv = v[[1, 2, 3]] / (2 - v[4])
  tuple(sv...)
end

# rotated tesseract for given angle `xi`
function tesseract(xi)
  # rotated and projected vertices
  vs = [Meshes.Point(stereog(rotate4d(0, 0, xi, vertices[i]))) for i = 1:16]
  ## edges
  function cylinder(edge)
    P = vs[edge[1]]
    Q = vs[edge[2]]
    Meshes.CylinderSurface(P, Q, 0.08)
  end
  tubes = [cylinder(edge) for edge in edges]
  function sphere(i)
    Meshes.Sphere(vs[i], 0.1)
  end
  mesh = reduce(merge, tubes)
  spheres = [Meshes.discretize(sphere(i), Meshes.RegularDiscretization(20)) for i in 1:16]
  merge(mesh, reduce(merge, spheres))
end

# draw tesseract
function drawTesseract(nplots, alpha1, alpha2)
  meshes = [tesseract(xi) for xi in LinRange(0, pi, nplots+1)[1:nplots]]
  for i in 1:nplots
    fig, ax, plt = MeshViz.viz( # invisible bounding box
      Meshes.Box(Meshes.Point(-2, -2, -2), Meshes.Point(2, 2, 2)); 
      alpha = 0
    )
    ax.show_axis = false
    MeshViz.viz!(meshes[i]; color = :crimson)
    Makie.scale!(ax.scene, 1.65, 1.65, 1.65)
    Makie.rotate!(fig.scene, Meshes.Vec3f(1, 0, 0), alpha1)
    Makie.rotate!(ax.scene, Meshes.Vec3f(0, 1, 0), alpha2)
    png = @sprintf "zzpic%03d.png" i
    Makie.save(png, fig)
  end
end