Cyclides and pseudo-cyclides with 'rgl' (R)

.gitignore

zzpic*

cyclideWithTexture.R

library(rgl)

mesh <- MeshesOperations::cyclideMesh(
  a = 97, c = 32, mu = 57, nu = 414L, nv = 414L
)

library(magick)
img <- image_read("ElectricityTexture.jpg")
rst <- as.raster(img)
mesh[["material"]] <- list(color = c(rst))

open3d(windowRect = c(50, 50, 562, 562))
bg3d(
  sphere = FALSE, texture = "SpaceBackground.png", col = "white"
)
view3d(0, 0, zoom = 0.8)
shade3d(mesh, shininess = 128)

# animation ####
M <- par3d("userMatrix")
movie3d(
  par3dinterp(
    time = seq(0, 1, len = 9),
    userMatrix = list(
      M,
      rotate3d(M, pi, 1, 0, 0),
      rotate3d(M, pi, 1, 1, 0),
      rotate3d(M, pi, 1, 1, 1),
      rotate3d(M, pi, 0, 1, 1),
      rotate3d(M, pi, 0, 1, 0),
      rotate3d(M, pi, 1, 0, 1),
      rotate3d(M, pi, 0, 0, 1),
      M
    )
  ),
  fps = 120,
  duration = 1,
  dir = ".",
  movie = "zzpic",
  convert = FALSE,
  webshot = FALSE
)

library(gifski)
gifski(
  png_files = Sys.glob("zzpic*.png"), 
  gif_file = "cyclideWithTexture.gif",
  width = 512,
  height = 512,
  delay = 1/12
)

pngfiles <- Sys.glob("zzpic*.png")
file.remove(pngfiles)

ElectricityTexture.jpg

SpaceBackground.png

TorturedCyclide.R

library(rgl)

scos <- function(x,alpha) sign(cos(x)) * abs(cos(x))^alpha
ssin <- function(x,alpha) sign(sin(x)) * abs(sin(x))^alpha

twistedTorusMesh <- function(alpha, beta, ntwists, R, r, S, s, arc, ...){
  vertices <- matrix(NA_real_, nrow = 3L, ncol = S*s)
  colors <- matrix(NA_character_, nrow = 4L, ncol = S*s)
  full <- arc == 2*pi
  SS <- ifelse(full, S, S-1)
  indices <- matrix(NA_integer_, nrow = 4L, ncol = SS*s)
  u_ <- if(full){
    seq(0, 2*pi, length.out = S+1)[-1L]
  }else{
    seq(0, arc, length.out = S)
  }
  v_ <- seq(0, 2*pi, length.out = s+1)[-1L]
  for(i in 1:S){
    cos_ui <- scos(u_[i], alpha)
    sin_ui <- ssin(u_[i], alpha)
    cos_2ui <- cos(ntwists*u_[i])
    sin_2ui <- sin(ntwists*u_[i])
    cx <- R * cos_ui
    cy <- R * sin_ui
    for(j in 1:s){
      rcos3_vj <- r * scos(v_[j], beta)
      rsin3_vj <- r * ssin(v_[j], beta)
      vertices[, (i-1)*s+j] <- c(
        cx + cos_ui * (rcos3_vj*cos_2ui + rsin3_vj*sin_2ui),
        cy + sin_ui * (rcos3_vj*cos_2ui + rsin3_vj*sin_2ui),
        rsin3_vj*cos_2ui - rcos3_vj*sin_2ui
      )
      colors[, (i-1)*s+j] <- 
        ifelse(floor(8*v_[j]/pi) %% 2 == 0, "#440154", "#FDE725")
    }
  }
  # quads
  s <- as.integer(s)
  for(i in 1L:SS){
    ip1 <- ifelse(i==S, 1L, i+1L)
    for(j in 1L:s){
      jp1 <- ifelse(j==s, 1L, j+1L)
      indices[,(i-1)*s+j] <- 
        c((i-1L)*s+j, (i-1L)*s+jp1, (ip1-1L)*s+jp1, (ip1-1L)*s+j)   
    }
  }
  out <- list(
    vb = rbind(vertices, 1),
    ib = indices,
    primitivetype = "quad",
    material = list(color=colors, ...),
    normals = NULL 
  )
  class(out) <- c("mesh3d", "shape3d")
  out
}

# only for mu>c (ring cyclide)
twistedCyclideMesh <- function(
    mu, a, c, ntwists=2, alpha=3, beta=2, S=128, s=64, arc=2*pi, ...
){
  b2 <- a * a - c * c;
  bb <- sqrt(b2 * (mu * mu - c * c))
  omega <- (a * mu + bb) / c
  Omega0 <- c(omega, 0, 0)
  inversion <- function(M) {
    OmegaM <- M - Omega0
    k <- c(crossprod(OmegaM))
    OmegaM / k + Omega0
  }
  h <- (c * c) / ((a - c) * (mu - c) + bb)
  r <- (h * (mu - c)) / ((a + c) * (mu - c) + bb)
  R <- (h * (a - c)) / ((a - c) * (mu + c) + bb)
  bb2 <- b2 * (mu * mu - c * c)
  denb1 <- c * (a*c - mu*c + c*c - a*mu - bb)
  b1 <- (a*mu*(c-mu)*(a+c) - bb2 + c*c + bb*(c*(a-mu+c) - 2*a*mu))/denb1
  denb2 <- c * (a*c - mu*c - c*c + a*mu + bb)
  b2 <- (a*mu*(c+mu)*(a-c) + bb2 - c*c + bb*(c*(a-mu-c) + 2*a*mu))/denb2
  omegaT <- (b1 + b2)/2
  OmegaT <- c(omegaT, 0, 0)
  tmesh <- twistedTorusMesh(alpha, beta, ntwists, R, r, S, s, arc, ...)
  tmesh$vb[1:3,] <- apply(tmesh$vb[1:3,] + OmegaT, 2, inversion)
  addNormals(tmesh)
}

# draw ####
# twisted torus
mesh <- twistedTorusMesh(alpha=1, ntwists=2, R=3, r=1, S=150, s=150, arc=2*pi)
open3d(windowRect = c(50,50,500,500))
bg3d(rgb(54,57,64, maxColorValue = 255))
view3d(0,-55)
shade3d(addNormals(mesh))

# twisted cyclide
a = 70; mu = 26; c = 14
mesh <- twistedCyclideMesh(
  mu, a, c, arc = 2*pi, S = 300, s = 150, alpha = 0.5, beta = 0.5, ntwists = 3
)

open3d(windowRect = c(50, 50, 562, 562))
bg3d(
  sphere = FALSE, texture = "SpaceBackground.png", col = "white"
)
view3d(0, -55, zoom = 0.75)
shade3d(mesh, meshColor = "legacy")

# animation ####
movie3d(spin3d(axis = c(0, 0, 1), rpm = 10),
        duration = 6, fps = 20,
        movie = "zzpic", dir = ".",
        convert = FALSE, webshot = FALSE,
        startTime = 1/20)

library(gifski)
gifski(
  png_files = Sys.glob("zzpic*.png"), 
  gif_file = "torturedCyclide.gif",
  width = 512,
  height = 512,
  delay = 1/13
)

pngfiles <- Sys.glob("zzpic*.png")
file.remove(pngfiles)

TwistedCyclide.R

library(rgl)

scos <- function(x,alpha) sign(cos(x)) * abs(cos(x))^alpha
ssin <- function(x,alpha) sign(sin(x)) * abs(sin(x))^alpha

twistedTorusMesh <- function(alpha, ntwists, R, r, S, s, arc, ...){
  vertices <- matrix(NA_real_, nrow = 3L, ncol = S*s)
  colors <- matrix(NA_character_, nrow = 4L, ncol = S*s)
  full <- arc == 2*pi
  SS <- ifelse(full, S, S-1)
  indices <- matrix(NA_integer_, nrow = 4L, ncol = SS*s)
  u_ <- if(full){
    seq(0, 2*pi, length.out = S+1)[-1L]
  }else{
    seq(0, arc, length.out = S)
  }
  v_ <- seq(0, 2*pi, length.out = s+1)[-1L]
  for(i in 1:S){
    cos_ui <- cos(u_[i])
    sin_ui <- sin(u_[i])
    cos_2ui <- cos(ntwists*u_[i])
    sin_2ui <- sin(ntwists*u_[i])
    cx <- R * cos_ui
    cy <- R * sin_ui
    for(j in 1:s){
      rcos3_vj <- r * scos(v_[j], alpha)
      rsin3_vj <- r * ssin(v_[j], alpha)
      vertices[, (i-1)*s+j] <- c(
        cx + cos_ui * (rcos3_vj*cos_2ui + rsin3_vj*sin_2ui),
        cy + sin_ui * (rcos3_vj*cos_2ui + rsin3_vj*sin_2ui),
        rsin3_vj*cos_2ui - rcos3_vj*sin_2ui
      )
      colors[, (i-1)*s+j] <- 
        ifelse(floor(8*v_[j]/pi) %% 2 == 0, "#440154", "#FDE725")
    }
  }
  # quads
  s <- as.integer(s)
  for(i in 1L:SS){
    ip1 <- ifelse(i==S, 1L, i+1L)
    for(j in 1L:s){
      jp1 <- ifelse(j==s, 1L, j+1L)
      indices[,(i-1)*s+j] <- 
        c((i-1L)*s+j, (i-1L)*s+jp1, (ip1-1L)*s+jp1, (ip1-1L)*s+j)   
    }
  }
  out <- list(
    vb = rbind(vertices, 1),
    ib = indices,
    primitivetype = "quad",
    material = list(color=colors, ...),
    normals = NULL 
  )
  class(out) <- c("mesh3d", "shape3d")
  out
}

# only for mu>c (ring cyclide)
twistedCyclideMesh <- function(
    mu, a, c, ntwists=2, alpha=3, S=128, s=64, arc=2*pi, ...
){
  b2 <- a * a - c * c;
  bb <- sqrt(b2 * (mu * mu - c * c))
  omega <- (a * mu + bb) / c
  Omega0 <- c(omega, 0, 0)
  inversion <- function(M) {
    OmegaM <- M - Omega0
    k <- c(crossprod(OmegaM))
    OmegaM / k + Omega0
  }
  h <- (c * c) / ((a - c) * (mu - c) + bb)
  r <- (h * (mu - c)) / ((a + c) * (mu - c) + bb)
  R <- (h * (a - c)) / ((a - c) * (mu + c) + bb)
  bb2 <- b2 * (mu * mu - c * c)
  denb1 <- c * (a*c - mu*c + c*c - a*mu - bb)
  b1 <- (a*mu*(c-mu)*(a+c) - bb2 + c*c + bb*(c*(a-mu+c) - 2*a*mu))/denb1
  denb2 <- c * (a*c - mu*c - c*c + a*mu + bb)
  b2 <- (a*mu*(c+mu)*(a-c) + bb2 - c*c + bb*(c*(a-mu-c) + 2*a*mu))/denb2
  omegaT <- (b1 + b2)/2
  OmegaT <- c(omegaT, 0, 0)
  tmesh <- twistedTorusMesh(alpha, ntwists, R, r, S, s, arc, ...)
  tmesh$vb[1:3,] <- apply(tmesh$vb[1:3,] + OmegaT, 2, inversion)
  addNormals(tmesh)
}

# draw ####
# twisted torus
mesh <- twistedTorusMesh(alpha=1, ntwists=2, R=3, r=1, S=150, s=150, arc=2*pi)
open3d(windowRect = c(50,50,500,500))
bg3d(rgb(54,57,64, maxColorValue = 255))
view3d(0,-55)
shade3d(addNormals(mesh))

# twisted cyclide
a = 94; mu = 56; c = 34
mesh <- twistedCyclideMesh(
  mu, a, c, arc = 2*pi, S = 300, s = 150, alpha = 3, ntwists = 2
)

open3d(windowRect = c(50, 50, 562, 562))
bg3d(
  sphere = FALSE, texture = "SpaceBackground.png", col = "white"
)
view3d(0, -55, zoom = 0.75)
shade3d(mesh, meshColor = "legacy")

# animation ####
movie3d(spin3d(axis = c(0, 0, 1), rpm = 10),
        duration = 6, fps = 20,
        movie = "zzpic", dir = ".",
        convert = FALSE, webshot = FALSE,
        startTime = 1/20)

library(gifski)
gifski(
  png_files = Sys.glob("zzpic*.png"), 
  gif_file = "twistedCyclide.gif",
  width = 512,
  height = 512,
  delay = 1/11
)

pngfiles <- Sys.glob("zzpic*.png")
file.remove(pngfiles)

cyclide_alpha-beta.R

cyclide_alpha-beta_0.5-0.5.gif

cyclide_alpha-beta_3-1.gif

cyclideWithTexture.gif

torturedCyclide.gif

twistedCyclide.gif