Duoprism with Julia

duoprism.jl

using Printf
using Makie

A = 4
B = 32

# construction of the vertices
vertices = Array{Float64,3}(undef, A, B, 4)
for i = 1:A
  v1 = [cospi(2 * i / A), sinpi(2 * i / A)]
  for j = 1:B
    v2 = [cospi(2 * j / B), sinpi(2 * j / B)]
    vertices[i, j, :] = vcat(v1, v2)
  end
end

# construction of the edges
function dominates(c1, c2)
  c2[1] > c1[1] || (c2[1] == c1[1] && c2[2] > c1[2])
end

function modulo(x, y)
  ((x % y) + y) % y
end

function getEdges()
  edges = Array{Int64,3}(undef, 2, 2, 2 * A * B)
  counter = 1
  for i = 0:(A-1)
    for j = 0:(B-1)
      c1 = [i, j]
      candidate = [i, modulo(j - 1, B)]
      if dominates(c1, candidate)
        edges[:, :, counter] = hcat(c1, candidate) .+ 1
        counter = counter + 1
      end
      candidate = [i, modulo(j + 1, B)]
      if dominates(c1, candidate)
        edges[:, :, counter] = hcat(c1, candidate) .+ 1
        counter = counter + 1
      end
      candidate = [modulo(i - 1, A), j]
      if dominates(c1, candidate)
        edges[:, :, counter] = hcat(c1, candidate) .+ 1
        counter = counter + 1
      end
      candidate = [modulo(i + 1, A), j]
      if dominates(c1, candidate)
        edges[:, :, counter] = hcat(c1, candidate) .+ 1
        counter = counter + 1
      end
    end
  end
  edges
end

edges = getEdges()

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)
  Meshes.Point(v[1:3] / (sqrt(2) - v[4]))
end

function duoprism(xi)
  # rotated and projected vertices
  vs = [stereog(rotate4d(pi/2, 0, xi, vertices[i, j, :])) for i = 1:A, j = 1:B]
  ## edges
  function cylinder(k)
    p1 = vs[edges[1, 1, k], edges[2, 1, k]]
    p2 = vs[edges[1, 2, k], edges[2, 2, k]]
    tubularmesh([p1, p2], 0.06, 20)
  end
  tubes = [cylinder(k) for k = 1:(2*A*B)]
  function sphere(i, j)
    p = vs[i, j]
    Meshes.Sphere(p, 0.1)
  end
  mesh = reduce(merge, tubes)
  spheres = [Meshes.discretize(sphere(i, j), Meshes.RegularDiscretization(20)) for i in 1:A, j in 1:B]
  merge(mesh, reduce(merge, spheres))
end

nplots = 60
meshes = [duoprism(xi) for xi in LinRange(0, 2*pi/4, nplots+1)[1:nplots]]

function draw(i)
  MeshViz.viz!(meshes[i]; color = :green)
end
for i in 1:nplots
  figg, axx, pltt = MeshViz.viz(Meshes.Box(Meshes.Point(-2.5, -2.5, -0.5), Meshes.Point(2.5, 2.5, 1.5)); alpha = 0)
  axx.show_axis = false
  draw(i)
  scale!(axx.scene, 1.6, 1.6, 1.6)
  Makie.rotate!(axx.scene, 1,0,0)
  png = @sprintf "duoprism%02d.png" i
  Makie.save(png, figg)
end

comm = @cmd "convert -layers OptimizePlus -delay 1x10 'duoprism*.png' duoprism.gif"
# run(comm)

duoprism2.jl

using Printf
import Makie
import Meshes 
import MeshViz
using GLMakie

A = 30 # number of vertices of the first polygon
B = 30 # number of vertices of the second polygon

# vertices of the duoprism
duoprismVertices = Array{Float64,3}(undef, A, B, 4)
for i = 1:A
  v1 = [cospi(2 * i / A), sinpi(2 * i / A)]
  for j = 1:B
    v2 = [cospi(2 * j / B), sinpi(2 * j / B)]
    duoprismVertices[i, j, :] = vcat(v1, v2)
  end
end

# construction of the edges of the duoprism
function dominates(c1, c2)
  c2[1] > c1[1] || (c2[1] == c1[1] && c2[2] > c1[2])
end

function modulo(x, y)
  ((x % y) + y) % y
end

function getEdges()
  edges = Array{Int64,3}(undef, 2, 2, 2 * A * B)
  counter = 1
  for i = 0:(A-1)
    for j = 0:(B-1)
      c1 = [i, j]
      candidate = [i, modulo(j - 1, B)]
      if dominates(c1, candidate)
        edges[:, :, counter] = hcat(c1, candidate) .+ 1
        counter = counter + 1
      end
      candidate = [i, modulo(j + 1, B)]
      if dominates(c1, candidate)
        edges[:, :, counter] = hcat(c1, candidate) .+ 1
        counter = counter + 1
      end
      candidate = [modulo(i - 1, A), j]
      if dominates(c1, candidate)
        edges[:, :, counter] = hcat(c1, candidate) .+ 1
        counter = counter + 1
      end
      candidate = [modulo(i + 1, A), j]
      if dominates(c1, candidate)
        edges[:, :, counter] = hcat(c1, candidate) .+ 1
        counter = counter + 1
      end
    end
  end
  edges
end

Edges = getEdges()

# 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] / (sqrt(2) - v[4])
end

# duoprism
function duoprism(xi)
  # rotated and projected vertices
  vs = [
    stereog(rotate4d(pi/2, 0, xi, duoprismVertices[i, j, :])) 
    for i = 1:A, j = 1:B
  ]
  # edge k
  function tube(k)
    p1 = vs[Edges[1, 1, k], Edges[2, 1, k]]
    p2 = vs[Edges[1, 2, k], Edges[2, 2, k]]
    Meshes.CylinderSurface(tuple(p1...), tuple(p2...), 0.06)
  end
  # all tubular edges
  tubes = [tube(k) for k = 1:(2*A*B)]
  # sphere i,j (vertex)
  function sphere(i, j)
    p = vs[i, j]
    Meshes.Sphere(tuple(p...), 0.1)
  end
  # all spherical vertices
  spheres = [
    Meshes.discretize(sphere(i, j), Meshes.RegularDiscretization(20)) 
    for i in 1:A, j in 1:B
  ]
  # merge everything in a single mesh
  mesh = Meshes.merge(
    reduce(Meshes.merge, spheres), reduce(Meshes.merge, tubes)
  )
end

# frames of the animation
nplots = 60
meshes = [duoprism(xi) for xi in LinRange(0, 2*pi/3, nplots+1)[1:nplots]]

function draw(i)
  MeshViz.viz!(meshes[i]; color = :yellow)
end
for i in 1:nplots
  fig, ax, plt = MeshViz.viz(
    Meshes.Box(Meshes.Point(-2.5, -2.5, -0.5), Meshes.Point(2.5, 2.5, 1.5)); 
    alpha = 0
  )
  ax.show_axis = false
  draw(i)
  Makie.scale!(ax.scene, 1.6, 1.6, 1.6)
  Makie.rotate!(ax.scene, 1.5, 0, 0)
  png = @sprintf "duoprism%02d.png" i
  Makie.save(png, fig)
end

# if you use ImageMagick
comm = @cmd "magick convert -layers OptimizePlus -delay 1x10 'duoprism*.png' duoprism.gif"
# run(comm)

duoprism3.jl

using Printf
import Makie
import Meshes 
import MeshViz
using GLMakie

# 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] / (sqrt(2) - v[4])
end

# helpers for the construction of the edges of the duoprism
function dominates(c1, c2)
  c2[1] > c1[1] || (c2[1] == c1[1] && c2[2] > c1[2])
end

function modulo(x, y)
  ((x % y) + y) % y
end

  
# A: number of vertices of the first polygon
# B: number of vertices of the second polygon

# vertices of the duoprism
function getVertices(A, B) 
   duoprismVertices = Array{Float64,3}(undef, A, B, 4)
  for i = 1:A
    v1 = [cospi(2 * i / A), sinpi(2 * i / A)]
    for j = 1:B
      v2 = [cospi(2 * j / B), sinpi(2 * j / B)]
      duoprismVertices[i, j, :] = vcat(v1, v2)
    end
  end
  duoprismVertices
end

# edges of the duoprism
function getEdges(A, B)
  edges = Array{Int64,3}(undef, 2, 2, 2 * A * B)
  counter = 1
  for i = 0:(A-1)
    for j = 0:(B-1)
      c1 = [i, j]
      candidate = [i, modulo(j - 1, B)]
      if dominates(c1, candidate)
        edges[:, :, counter] = hcat(c1, candidate) .+ 1
        counter = counter + 1
      end
      candidate = [i, modulo(j + 1, B)]
      if dominates(c1, candidate)
        edges[:, :, counter] = hcat(c1, candidate) .+ 1
        counter = counter + 1
      end
      candidate = [modulo(i - 1, A), j]
      if dominates(c1, candidate)
        edges[:, :, counter] = hcat(c1, candidate) .+ 1
        counter = counter + 1
      end
      candidate = [modulo(i + 1, A), j]
      if dominates(c1, candidate)
        edges[:, :, counter] = hcat(c1, candidate) .+ 1
        counter = counter + 1
      end
    end
  end
  edges
end


# duoprism for given angle `xi`
function duoprism(xi, vertices, edges)
  A, B, _ = size(vertices)
  # rotated and projected vertices
  vs = [
    stereog(rotate4d(pi/2, 0, xi, vertices[i, j, :])) 
    for i = 1:A, j = 1:B
  ]
  # edge k
  function tube(k)
    p1 = vs[edges[1, 1, k], edges[2, 1, k]]
    p2 = vs[edges[1, 2, k], edges[2, 2, k]]
    Meshes.CylinderSurface(tuple(p1...), tuple(p2...), 0.06)
  end
  # all tubular edges
  tubes = [tube(k) for k = 1:(2*A*B)]
  # sphere i,j (vertex)
  function sphere(i, j)
    p = vs[i, j]
    Meshes.Sphere(tuple(p...), 0.1)
  end
  # all spherical vertices
  spheres = [
    Meshes.discretize(sphere(i, j), Meshes.RegularDiscretization(20)) 
    for i in 1:A, j in 1:B
  ]
  # merge everything in a single mesh
  mesh = Meshes.merge(
    reduce(Meshes.merge, spheres), reduce(Meshes.merge, tubes)
  )
end

# make duoprism 
A = 3
B = 30
Vertices = getVertices(A, B)
Edges = getEdges(A, B)
function duoprismMesh(angle, nplots) 
  meshes = [
    duoprism(xi, Vertices, Edges) 
    for xi in LinRange(0, angle, nplots+1)[1:nplots]
  ]
end

# frames of the animation
angle = 2*pi/3
nplots = 60
meshes = duoprismMesh(angle, nplots)

# draw one frame
function drawFrame(i)
  MeshViz.viz!(meshes[i]; color = :yellow)
end

function drawDuoprism(alpha1, alpha2)
  for i in 1:nplots
    fig, ax, plt = MeshViz.viz(
      Meshes.Box(Meshes.Point(-2.5, -2.5, -0.5), Meshes.Point(2.5, 2.5, 1.5)); 
      alpha = 0
    )
    ax.show_axis = false
    drawFrame(i)
    Makie.scale!(ax.scene, 1.6, 1.6, 1.6)
    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

# if you use ImageMagick
comm = @cmd "magick convert -layers OptimizePlus -delay 1x10 'zzpic*.png' duoprism.gif"
# run(comm)