Site icon R-bloggers

Stereographic truncated tesseract

[This article was first published on Saturn Elephant, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

The stereographic truncated tesseract

We show how to draw a stereographic truncated tesseract with rgl (R), Asymptote, and POV-Ray.

The truncated tesseract is a uniform polychoron. Among its cells, there are sixteen tetrahedra, and we fill the faces of these tetrahedra.

Using Asymptote and POV-Ray, we include a 4D rotation to animate the figure.

We call it stereographic because we map each edge to the 3-sphere before applying the stereographic projection.

Drawing with rgl (R)

library(rgl)
library(cxhull)
library(abind)

# vertices #####################################################################
sqr2p1 <- sqrt(2) + 1
vertices <- rbind(
  c(-1, -sqr2p1, -sqr2p1, -sqr2p1),
  c(-1, -sqr2p1, -sqr2p1, sqr2p1),
  c(-1, -sqr2p1, sqr2p1, -sqr2p1),
  c(-1, -sqr2p1, sqr2p1, sqr2p1),
  c(-1, sqr2p1, -sqr2p1, -sqr2p1),
  c(-1, sqr2p1, -sqr2p1, sqr2p1),
  c(-1, sqr2p1, sqr2p1, -sqr2p1),
  c(-1, sqr2p1, sqr2p1, sqr2p1),
  c(1, -sqr2p1, -sqr2p1, -sqr2p1),
  c(1, -sqr2p1, -sqr2p1, sqr2p1),
  c(1, -sqr2p1, sqr2p1, -sqr2p1),
  c(1, -sqr2p1, sqr2p1, sqr2p1),
  c(1, sqr2p1, -sqr2p1, -sqr2p1),
  c(1, sqr2p1, -sqr2p1, sqr2p1),
  c(1, sqr2p1, sqr2p1, -sqr2p1),
  c(1, sqr2p1, sqr2p1, sqr2p1),
  c(-sqr2p1, -1, -sqr2p1, -sqr2p1),
  c(-sqr2p1, -1, -sqr2p1, sqr2p1),
  c(-sqr2p1, -1, sqr2p1, -sqr2p1),
  c(-sqr2p1, -1, sqr2p1, sqr2p1),
  c(-sqr2p1, 1, -sqr2p1, -sqr2p1),
  c(-sqr2p1, 1, -sqr2p1, sqr2p1),
  c(-sqr2p1, 1, sqr2p1, -sqr2p1),
  c(-sqr2p1, 1, sqr2p1, sqr2p1),
  c(sqr2p1, -1, -sqr2p1, -sqr2p1),
  c(sqr2p1, -1, -sqr2p1, sqr2p1),
  c(sqr2p1, -1, sqr2p1, -sqr2p1),
  c(sqr2p1, -1, sqr2p1, sqr2p1),
  c(sqr2p1, 1, -sqr2p1, -sqr2p1),
  c(sqr2p1, 1, -sqr2p1, sqr2p1),
  c(sqr2p1, 1, sqr2p1, -sqr2p1),
  c(sqr2p1, 1, sqr2p1, sqr2p1),
  c(-sqr2p1, -sqr2p1, -1, -sqr2p1),
  c(-sqr2p1, -sqr2p1, -1, sqr2p1),
  c(-sqr2p1, -sqr2p1, 1, -sqr2p1),
  c(-sqr2p1, -sqr2p1, 1, sqr2p1),
  c(-sqr2p1, sqr2p1, -1, -sqr2p1),
  c(-sqr2p1, sqr2p1, -1, sqr2p1),
  c(-sqr2p1, sqr2p1, 1, -sqr2p1),
  c(-sqr2p1, sqr2p1, 1, sqr2p1),
  c(sqr2p1, -sqr2p1, -1, -sqr2p1),
  c(sqr2p1, -sqr2p1, -1, sqr2p1),
  c(sqr2p1, -sqr2p1, 1, -sqr2p1),
  c(sqr2p1, -sqr2p1, 1, sqr2p1),
  c(sqr2p1, sqr2p1, -1, -sqr2p1),
  c(sqr2p1, sqr2p1, -1, sqr2p1),
  c(sqr2p1, sqr2p1, 1, -sqr2p1),
  c(sqr2p1, sqr2p1, 1, sqr2p1),
  c(-sqr2p1, -sqr2p1, -sqr2p1, -1),
  c(-sqr2p1, -sqr2p1, -sqr2p1, 1),
  c(-sqr2p1, -sqr2p1, sqr2p1, -1),
  c(-sqr2p1, -sqr2p1, sqr2p1, 1),
  c(-sqr2p1, sqr2p1, -sqr2p1, -1),
  c(-sqr2p1, sqr2p1, -sqr2p1, 1),
  c(-sqr2p1, sqr2p1, sqr2p1, -1),
  c(-sqr2p1, sqr2p1, sqr2p1, 1),
  c(sqr2p1, -sqr2p1, -sqr2p1, -1),
  c(sqr2p1, -sqr2p1, -sqr2p1, 1),
  c(sqr2p1, -sqr2p1, sqr2p1, -1),
  c(sqr2p1, -sqr2p1, sqr2p1, 1),
  c(sqr2p1, sqr2p1, -sqr2p1, -1),
  c(sqr2p1, sqr2p1, -sqr2p1, 1),
  c(sqr2p1, sqr2p1, sqr2p1, -1),
  c(sqr2p1, sqr2p1, sqr2p1, 1)
)

# convex hull is the polytope itself, since it is convex #######################
hull <- cxhull(vertices)

# edges ########################################################################
edges <- hull$edges

# (triangular) faces of the tetrahedra ######################################### 
ridgeSizes <- sapply(hull$ridges, function(ridge) length(ridge$vertices))
triangles <- t(sapply(hull$ridges[which(ridgeSizes==3)], 
                      function(ridge) ridge$vertices))

# stereographic projection #####################################################
sproj <- function(p, r){
  c(p[1], p[2], p[3])/(r-p[4])
}

# spherical segment ############################################################
sphericalSegment <- function(P, Q, r, n){
  out <- matrix(NA_real_, nrow = n+1, ncol = 4)
  for(i in 0:n){
    pt <- P + (i/n)*(Q-P)
    out[i+1L, ] <- r/sqrt(c(crossprod(pt))) * pt
  }
  out
}

# stereographic edge ###########################################################
stereoEdge <- function(verts, edge, r, n){
  P <- verts[edge[1], ]
  Q <- verts[edge[2], ]
  PQ <- sphericalSegment(P, Q, r, n)
  pq <- t(apply(PQ, 1, function(M){sproj(M, r)}))
  dists <- sqrt(apply(pq, 1, crossprod))
  cylinder3d(pq, radius = log1p(dists/4)/4, sides = 60)
}

# stereographic subdivision (to fill the triangles) ############################
midpoint4 <- function(A, B, r){
  M <- (A + B) / 2
  lg <- sqrt(c(crossprod(M))) / r
  M / lg
} 

subdiv0 <- function(A, B, C, r){
  mAB <- midpoint4(A, B, r)
  mAC <- midpoint4(A, C, r)
  mBC <- midpoint4(B, C, r)
  out <- array(NA_real_, dim = c(4L, 3L, 4L))
  out[1,,] <- rbind(A, mAB, mAC)
  out[2,,] <- rbind(B, mBC, mAB)
  out[3,,] <- rbind(C, mAC, mBC)
  out[4,,] <- rbind(mAB, mBC, mAC)
  out
}

subdiv <- function(n, A, B, C, r){
  if(n == 1){
    return(subdiv0(A, B, C, r))
  }
  triangles <- subdiv(n-1, A, B, C, r)
  out <- array(NA_real_, dim = c(0L, 3L, 4L))
  for(i in 1:(4^(n-1))){
    trgl <- triangles[i,,]
    out <- abind(out, subdiv0(trgl[1L,], trgl[2L,], trgl[3L,], r), along = 1L)
  }
  out
}

# mesh maker ###################################################################
stereoTriangle <- function(n, A, B, C, r){
  triangles <- subdiv(n, A, B, C, r)
  ntriangles <- dim(triangles)[1L]
  indices <- matrix(1L:(3L*ntriangles), nrow = 3L, ncol = ntriangles)
  vertices <- matrix(NA_real_, nrow = 3L, ncol = 3L*ntriangles)
  for(i in 1L:ntriangles){
    trgl <- triangles[i,,]
    vertices[, 3L*(i-1L)+1L] <- sproj(trgl[1L,], r)
    vertices[, 3L*(i-1L)+2L] <- sproj(trgl[2L,], r)
    vertices[, 3L*(i-1L)+3L] <- sproj(trgl[3L,], r)
  }
  mesh <- tmesh3d(vertices, indices, homogeneous = FALSE)
  Rvcg::vcgClean(mesh, sel = c(0,7), silent = TRUE)
}

# projected vertices ###########################################################
r <- sqrt(1+3*sqr2p1^2)
vs <- t(apply(vertices, 1L, function(M){sproj(M, r)}))

####~~~~ plot ~~~~##############################################################
open3d(windowRect = c(50, 50, 562, 562), zoom = 0.9)
bg3d(rgb(54, 57, 64, maxColorValue = 255))
## plot the edges
for(k in 1L:nrow(edges)){
  edge <- stereoEdge(vertices, edges[k,], r, n = 100)
  shade3d(edge, color = "gold")
}
## plot the vertices
for(i in 1L:nrow(vs)){
  v <- vs[i,]
  spheres3d(v, radius = log1p(sqrt(c(crossprod(v)))/4)/3, color = "gold2")
}
## plot the filled triangles
for(i in 1L:nrow(triangles)){
  trgl <- triangles[i,]
  A <- vertices[trgl[1L], ]
  B <- vertices[trgl[2L], ]
  C <- vertices[trgl[3L], ]
  mesh <- stereoTriangle(6, A, B, C, r)
  shade3d(mesh, color = "red")
}

Drawing with Asymptote

settings.render = 4;
settings.outformat = "eps";
import tube;
size(200,0);

currentprojection = orthographic(1, 2, 4);
currentlight = light(gray(0.85), ambient=black, specularfactor=3,
                     (100,100,100), specular=gray(0.9), viewport=false);
currentlight.background = rgb("363940ff");

// files to be saved -----------------------------------------------------------
string[] files = {
    "TT000", "TT001", "TT002", "TT003", "TT004", "TT005",
    "TT006", "TT007", "TT008", "TT009", "TT010", "TT011",
    "TT012", "TT013", "TT014", "TT015", "TT016", "TT017",
    "TT018", "TT019", "TT020", "TT021", "TT022", "TT023",
    "TT024", "TT025", "TT026", "TT027", "TT028", "TT029",
    "TT030", "TT031", "TT032", "TT033", "TT034", "TT035",
    "TT036", "TT037", "TT038", "TT039", "TT040", "TT041",
    "TT042", "TT043", "TT044", "TT045", "TT046", "TT047",
    "TT048", "TT049", "TT050", "TT051", "TT052", "TT053",
    "TT054", "TT055", "TT056", "TT057", "TT058", "TT059",
    "TT060", "TT061", "TT062", "TT063", "TT064", "TT065",
    "TT066", "TT067", "TT068", "TT069", "TT070", "TT071",
    "TT072", "TT073", "TT074", "TT075", "TT076", "TT077",
    "TT078", "TT079", "TT080", "TT081", "TT082", "TT083",
    "TT084", "TT085", "TT086", "TT087", "TT088", "TT089",
    "TT090", "TT091", "TT092", "TT093", "TT094", "TT095",
    "TT096", "TT097", "TT098", "TT099", "TT100", "TT101",
    "TT102", "TT103", "TT104", "TT105", "TT106", "TT107",
    "TT108", "TT109", "TT110", "TT111", "TT112", "TT113",
    "TT114", "TT115", "TT116", "TT117", "TT118", "TT119",
    "TT120", "TT121", "TT122", "TT123", "TT124", "TT125",
    "TT126", "TT127", "TT128", "TT129", "TT130", "TT131",
    "TT132", "TT133", "TT134", "TT135", "TT136", "TT137",
    "TT138", "TT139", "TT140", "TT141", "TT142", "TT143",
    "TT144", "TT145", "TT146", "TT147", "TT148", "TT149",
    "TT150", "TT151", "TT152", "TT153", "TT154", "TT155",
    "TT156", "TT157", "TT158", "TT159", "TT160", "TT161",
    "TT162", "TT163", "TT164", "TT165", "TT166", "TT167",
    "TT168", "TT169", "TT170", "TT171", "TT172", "TT173",
    "TT174", "TT175", "TT176", "TT177", "TT178", "TT179"};

// vertices --------------------------------------------------------------------
struct quadruple {
    real x;
    real y;
    real z;
    real t;
}

quadruple array2quadruple(real[] v){
    quadruple q;
    q.x = v[0]; q.y = v[1]; q.z = v[2]; q.t = v[3];
    return q;
}

real x = 1 + sqrt(2);
real r = sqrt(1 + 3 * x * x);
real[][] vertices0 = {
    {-1.0, -x, -x, -x},
    {-1.0, -x, -x, x},
    {-1.0, -x, x, -x},
    {-1.0, -x, x, x},
    {-1.0, x, -x, -x},
    {-1.0, x, -x, x},
    {-1.0, x, x, -x},
    {-1.0, x, x, x},
    {1.0, -x, -x, -x},
    {1.0, -x, -x, x},
    {1.0, -x, x, -x},
    {1.0, -x, x, x},
    {1.0, x, -x, -x},
    {1.0, x, -x, x},
    {1.0, x, x, -x},
    {1.0, x, x, x},
    {-x, -1.0, -x, -x},
    {-x, -1.0, -x, x},
    {-x, -1.0, x, -x},
    {-x, -1.0, x, x},
    {-x, 1.0, -x, -x},
    {-x, 1.0, -x, x},
    {-x, 1.0, x, -x},
    {-x, 1.0, x, x},
    {x, -1.0, -x, -x},
    {x, -1.0, -x, x},
    {x, -1.0, x, -x},
    {x, -1.0, x, x},
    {x, 1.0, -x, -x},
    {x, 1.0, -x, x},
    {x, 1.0, x, -x},
    {x, 1.0, x, x},
    {-x, -x, -1.0, -x},
    {-x, -x, -1.0, x},
    {-x, -x, 1.0, -x},
    {-x, -x, 1.0, x},
    {-x, x, -1.0, -x},
    {-x, x, -1.0, x},
    {-x, x, 1.0, -x},
    {-x, x, 1.0, x},
    {x, -x, -1.0, -x},
    {x, -x, -1.0, x},
    {x, -x, 1.0, -x},
    {x, -x, 1.0, x},
    {x, x, -1.0, -x},
    {x, x, -1.0, x},
    {x, x, 1.0, -x},
    {x, x, 1.0, x},
    {-x, -x, -x, -1.0},
    {-x, -x, -x, 1.0},
    {-x, -x, x, -1.0},
    {-x, -x, x, 1.0},
    {-x, x, -x, -1.0},
    {-x, x, -x, 1.0},
    {-x, x, x, -1.0},
    {-x, x, x, 1.0},
    {x, -x, -x, -1.0},
    {x, -x, -x, 1.0},
    {x, -x, x, -1.0},
    {x, -x, x, 1.0},
    {x, x, -x, -1.0},
    {x, x, -x, 1.0},
    {x, x, x, -1.0},
    {x, x, x, 1.0} };
    
quadruple[] vertices = new quadruple[vertices0.length];
for(int i = 0; i < vertices0.length; ++i){
    vertices[i] = array2quadruple(vertices0[i]);
}

// edges -----------------------------------------------------------------------
int[][] edges = {
    {0, 8},
    {0, 16},
    {0, 32},
    {0, 48},
    {1, 9},
    {1, 17},
    {1, 33},
    {1, 49},
    {2, 10},
    {2, 18},
    {2, 34},
    {2, 50},
    {3, 11},
    {3, 19},
    {3, 35},
    {3, 51},
    {4, 12},
    {4, 20},
    {4, 36},
    {4, 52},
    {5, 13},
    {5, 21},
    {5, 37},
    {5, 53},
    {6, 14},
    {6, 22},
    {6, 38},
    {6, 54},
    {7, 15},
    {7, 23},
    {7, 39},
    {7, 55},
    {8, 24},
    {8, 40},
    {8, 56},
    {9, 25},
    {9, 41},
    {9, 57},
    {10, 26},
    {10, 42},
    {10, 58},
    {11, 27},
    {11, 43},
    {11, 59},
    {12, 28},
    {12, 44},
    {12, 60},
    {13, 29},
    {13, 45},
    {13, 61},
    {14, 30},
    {14, 46},
    {14, 62},
    {15, 31},
    {15, 47},
    {15, 63},
    {16, 20},
    {16, 32},
    {16, 48},
    {17, 21},
    {17, 33},
    {17, 49},
    {18, 22},
    {18, 34},
    {18, 50},
    {19, 23},
    {19, 35},
    {19, 51},
    {20, 36},
    {20, 52},
    {21, 37},
    {21, 53},
    {22, 38},
    {22, 54},
    {23, 39},
    {23, 55},
    {24, 28},
    {24, 40},
    {24, 56},
    {25, 29},
    {25, 41},
    {25, 57},
    {26, 30},
    {26, 42},
    {26, 58},
    {27, 31},
    {27, 43},
    {27, 59},
    {28, 44},
    {28, 60},
    {29, 45},
    {29, 61},
    {30, 46},
    {30, 62},
    {31, 47},
    {31, 63},
    {32, 34},
    {32, 48},
    {33, 35},
    {33, 49},
    {34, 50},
    {35, 51},
    {36, 38},
    {36, 52},
    {37, 39},
    {37, 53},
    {38, 54},
    {39, 55},
    {40, 42},
    {40, 56},
    {41, 43},
    {41, 57},
    {42, 58},
    {43, 59},
    {44, 46},
    {44, 60},
    {45, 47},
    {45, 61},
    {46, 62},
    {47, 63},
    {48, 49},
    {50, 51},
    {52, 53},
    {54, 55},
    {56, 57},
    {58, 59},
    {60, 61},
    {62, 63} };

// tetrahedra ------------------------------------------------------------------
int[][] tetrahedra = {
    {0, 16, 32, 48},
    {11, 27, 43, 59},
    {12, 28, 44, 60},
    {8, 24, 40, 56},
    {9, 25, 41, 57},
    {15, 31, 47, 63},
    {13, 29, 45, 61},
    {14, 30, 46, 62},
    {10, 26, 42, 58},
    {3, 19, 35, 51},
    {2, 18, 34, 50},
    {1, 17, 33, 49},
    {4, 20, 36, 52},
    {5, 21, 37, 53},
    {6, 22, 38, 54},
    {7, 23, 39, 55} };

// rotation in 4D space (right-isoclinic) --------------------------------------
quadruple rotate4d(real alpha, real beta, real xi, quadruple vec){
    real a = cos(xi);
    real b = sin(alpha)*cos(beta)*sin(xi);
    real c = sin(alpha)*sin(beta)*sin(xi);
    real d = cos(alpha)*sin(xi);
    real p = vec.x;
    real q = vec.y;
    real r = vec.z;
    real s = vec.t;
    quadruple out;
    out.x = a*p - b*q - c*r - d*s;
    out.y = a*q + b*p + c*s - d*r;
    out.z = a*r - b*s + c*p + d*q;
    out.t = a*s + b*r - c*q + d*p;
    return out;
}


// stereographic projection ----------------------------------------------------
triple stereog(quadruple A, real r){
    return (A.x, A.y, A.z) / (r - A.t);
}

// stereographic path ----------------------------------------------------------
path3 stereoPath(quadruple A, quadruple B, real r, int n){
    path3 out;
    for(int i = 0; i <= n; ++i){
        real t = i/n;
        quadruple M;
        real x = (1-t)*A.x + t*B.x;
        real y = (1-t)*A.y + t*B.y;
        real z = (1-t)*A.z + t*B.z;
        real t = (1-t)*A.t + t*B.t;
        real lg = sqrt(x*x + y*y + z*z + t*t) / r;
        M.x = x / lg; M.y = y / lg; M.z = z / lg; M.t = t / lg;
        out = out .. stereog(M, r);
    }
    return out;
}

// section transformation ------------------------------------------------------
transform T(path3 p3, real t, int n){
    triple M = relpoint(p3, t/(n/4));
    return scale(length(M)/15); 
}

// stereographic subdivision (to fill the tetrahedra) --------------------------
quadruple midpoint4(quadruple A, quadruple B, real r){
    quadruple M;
    real x = (A.x + B.x)/2;
    real y = (A.y + B.y)/2;
    real z = (A.z + B.z)/2;
    real t = (A.t + B.t)/2;
    real lg = sqrt(x*x + y*y + z*z + t*t) / r;
    M.x = x / lg; M.y = y / lg; M.z = z / lg; M.t = t / lg;
    return(M);
} 

quadruple[][] subdiv0(quadruple A, quadruple B, quadruple C, real r){
    quadruple m01 = midpoint4(A, B, r);
    quadruple m02 = midpoint4(A, C, r);
    quadruple m12 = midpoint4(B, C, r);
    return new quadruple[][] {
        {A, m01, m02}, 
        {B, m12, m01}, 
        {C, m02, m12}, 
        {m01, m12, m02}
    };
}

quadruple[][] subdiv(int n, quadruple A, quadruple B, quadruple C, real r){
  if(n == 1){
    return subdiv0(A, B, C, r);
  }
  quadruple[][] triangles = subdiv(n-1, A, B, C, r);
  quadruple[][] out = new quadruple[0][3];
  for(int i = 0; i < 4^(n-1); ++i){
    quadruple[] trgl = triangles[i];
    out.append(subdiv0(trgl[0], trgl[1], trgl[2], r));
  }
  return out;
}

// mesh maker ------------------------------------------------------------------
struct Mesh {
    triple[] vertices;
    int[][] indices;
}
Mesh stereoTriangle(int n, quadruple A, quadruple B, quadruple C, real r){
    quadruple[][] triangles = subdiv(n, A, B, C, r);
    triple[] vertices;
    int[][] indices;
    int faceIndex = 0;
    for(int i = 0; i < triangles.length; ++i){
        quadruple[] triangle = triangles[i];
        vertices.push(stereog(triangle[0], r));
        vertices.push(stereog(triangle[1], r));
        vertices.push(stereog(triangle[2], r));
        int[] x = {faceIndex, faceIndex+1, faceIndex+2};
        indices.push(x);
        faceIndex += 3;
    }
    Mesh out;
    out.vertices = vertices;
    out.indices = indices;
    return out;
}


// "bounding box" (to clip the animation) --------------------------------------
path3 boundingbox = scale(0.5,1,1)*circle(c=O, r=4, normal=Z);

// plot ------------------------------------------------------------------------
int n = 75;
int depth = 6;
real alpha = 0, beta = 0;

for(int f = 0; f < 180; ++f){
    // rotation angle
    real xi = 2*f*pi/180;
    // rotated vertices
    quadruple[] vs = new quadruple[vertices.length];
    for(int i = 0; i < vertices.length; ++i){
        vs[i] = rotate4d(alpha, beta, xi, vertices[i]);
    }
    // new picture
    picture pic;
    // draw boundingbox
    draw(pic, boundingbox, rgb("363940ff")+opacity(0));
    // draw edges
    for(int k = 0; k < edges.length; ++k){
        quadruple A = vs[edges[k][0]];
        quadruple B = vs[edges[k][1]];
        path3 p3 = stereoPath(A, B, r, n);
        transform S(real t){
            return T(p3, t, n);
        }
        draw(pic, tube(p3, unitcircle, S), rgb(139, 0, 139), 
                render(compression=Low, merge=true));
    }
    // draw vertices
    for(int i = 0; i < vertices.length; ++i){
        triple Asg = stereog(vs[i], r);
        draw(pic, shift(Asg)*scale3(length(Asg)/10)*unitsphere, purple);
    }
    // draw tetrahedra
    for(int i = 0; i < tetrahedra.length; ++i){
        int[] t = tetrahedra[i];
        Mesh mesh = stereoTriangle(depth, vs[t[0]], vs[t[1]], vs[t[2]], r);
        draw(pic, mesh.vertices, mesh.indices, m = red);
        Mesh mesh = stereoTriangle(depth, vs[t[0]], vs[t[1]], vs[t[3]], r);
        draw(pic, mesh.vertices, mesh.indices, m = red);
        Mesh mesh = stereoTriangle(depth, vs[t[0]], vs[t[2]], vs[t[3]], r);
        draw(pic, mesh.vertices, mesh.indices, m = red);
        Mesh mesh = stereoTriangle(depth, vs[t[1]], vs[t[2]], vs[t[3]], r);
        draw(pic, mesh.vertices, mesh.indices, m = red);
    }
    // add picture and save
    add(pic);
    shipout(files[f], bbox(rgb("363940ff"), FillDraw(rgb("363940ff"))));
    erase();
}   

Drawing with POV-Ray

#version 3.7;
global_settings { assumed_gamma 1 }
#include "colors.inc"
#include "textures.inc"

/* camera */
camera {
    location <-11, 7,-32>
    look_at 0
    angle 45
}

// sun -------------------------------------------------------------------------
light_source {< 1000,3000,-6000> color rgb<1,1,1>*0.9}             // sun 
light_source {<-11, 7,-32>  color rgb<0.9,0.9,1>*0.1 shadowless}   // flash

// sky -------------------------------------------------------------------------
plane{<0,1,0>,1 hollow  
       texture{ pigment{ bozo turbulence 1.3
                         color_map { [0.00 rgb <0.24, 0.32, 1.0>*0.6]
                                     [0.75 rgb <0.24, 0.32, 1.0>*0.6]
                                     [0.83 rgb <1,1,1>]
                                     [0.95 rgb <0.25,0.25,0.25>]
                                     [1.0 rgb <0.5,0.5,0.5>]}
                        scale<1,1,1>*2.5  translate< 0,0,3>
                       }
                finish {ambient 1 diffuse 0} }      
       scale 10000}

// fog on the ground -----------------------------------------------------------
fog { fog_type   2
      distance   50
      color      Gray10  
      fog_offset 0.1
      fog_alt    1.5
      turbulence 1.8
    }

// ground ----------------------------------------------------------------------
plane { <0,1,0>, 0 
    texture {
        pigment { color rgb <0.95,0.9,0.73>*0.35 }
        normal { bumps 2 scale <0.25,0.25,0.25>*0.5 turbulence 0.5 } 
        finish { phong 0.1 }
    } 
} 


/* vertices */
#declare sqr2p1 = sqrt(2) + 1;
#declare v0 = < -1.0 , -sqr2p1 , -sqr2p1 , -sqr2p1 >;
#declare v1 = < -1.0 , -sqr2p1 , -sqr2p1 , sqr2p1 >;
#declare v2 = < -1.0 , -sqr2p1 , sqr2p1 , -sqr2p1 >;
#declare v3 = < -1.0 , -sqr2p1 , sqr2p1 , sqr2p1 >;
#declare v4 = < -1.0 , sqr2p1 , -sqr2p1 , -sqr2p1 >;
#declare v5 = < -1.0 , sqr2p1 , -sqr2p1 , sqr2p1 >;
#declare v6 = < -1.0 , sqr2p1 , sqr2p1 , -sqr2p1 >;
#declare v7 = < -1.0 , sqr2p1 , sqr2p1 , sqr2p1 >;
#declare v8 = < 1.0 , -sqr2p1 , -sqr2p1 , -sqr2p1 >;
#declare v9 = < 1.0 , -sqr2p1 , -sqr2p1 , sqr2p1 >;
#declare v10 = < 1.0 , -sqr2p1 , sqr2p1 , -sqr2p1 >;
#declare v11 = < 1.0 , -sqr2p1 , sqr2p1 , sqr2p1 >;
#declare v12 = < 1.0 , sqr2p1 , -sqr2p1 , -sqr2p1 >;
#declare v13 = < 1.0 , sqr2p1 , -sqr2p1 , sqr2p1 >;
#declare v14 = < 1.0 , sqr2p1 , sqr2p1 , -sqr2p1 >;
#declare v15 = < 1.0 , sqr2p1 , sqr2p1 , sqr2p1 >;
#declare v16 = < -sqr2p1 , -1.0 , -sqr2p1 , -sqr2p1 >;
#declare v17 = < -sqr2p1 , -1.0 , -sqr2p1 , sqr2p1 >;
#declare v18 = < -sqr2p1 , -1.0 , sqr2p1 , -sqr2p1 >;
#declare v19 = < -sqr2p1 , -1.0 , sqr2p1 , sqr2p1 >;
#declare v20 = < -sqr2p1 , 1.0 , -sqr2p1 , -sqr2p1 >;
#declare v21 = < -sqr2p1 , 1.0 , -sqr2p1 , sqr2p1 >;
#declare v22 = < -sqr2p1 , 1.0 , sqr2p1 , -sqr2p1 >;
#declare v23 = < -sqr2p1 , 1.0 , sqr2p1 , sqr2p1 >;
#declare v24 = < sqr2p1 , -1.0 , -sqr2p1 , -sqr2p1 >;
#declare v25 = < sqr2p1 , -1.0 , -sqr2p1 , sqr2p1 >;
#declare v26 = < sqr2p1 , -1.0 , sqr2p1 , -sqr2p1 >;
#declare v27 = < sqr2p1 , -1.0 , sqr2p1 , sqr2p1 >;
#declare v28 = < sqr2p1 , 1.0 , -sqr2p1 , -sqr2p1 >;
#declare v29 = < sqr2p1 , 1.0 , -sqr2p1 , sqr2p1 >;
#declare v30 = < sqr2p1 , 1.0 , sqr2p1 , -sqr2p1 >;
#declare v31 = < sqr2p1 , 1.0 , sqr2p1 , sqr2p1 >;
#declare v32 = < -sqr2p1 , -sqr2p1 , -1.0 , -sqr2p1 >;
#declare v33 = < -sqr2p1 , -sqr2p1 , -1.0 , sqr2p1 >;
#declare v34 = < -sqr2p1 , -sqr2p1 , 1.0 , -sqr2p1 >;
#declare v35 = < -sqr2p1 , -sqr2p1 , 1.0 , sqr2p1 >;
#declare v36 = < -sqr2p1 , sqr2p1 , -1.0 , -sqr2p1 >;
#declare v37 = < -sqr2p1 , sqr2p1 , -1.0 , sqr2p1 >;
#declare v38 = < -sqr2p1 , sqr2p1 , 1.0 , -sqr2p1 >;
#declare v39 = < -sqr2p1 , sqr2p1 , 1.0 , sqr2p1 >;
#declare v40 = < sqr2p1 , -sqr2p1 , -1.0 , -sqr2p1 >;
#declare v41 = < sqr2p1 , -sqr2p1 , -1.0 , sqr2p1 >;
#declare v42 = < sqr2p1 , -sqr2p1 , 1.0 , -sqr2p1 >;
#declare v43 = < sqr2p1 , -sqr2p1 , 1.0 , sqr2p1 >;
#declare v44 = < sqr2p1 , sqr2p1 , -1.0 , -sqr2p1 >;
#declare v45 = < sqr2p1 , sqr2p1 , -1.0 , sqr2p1 >;
#declare v46 = < sqr2p1 , sqr2p1 , 1.0 , -sqr2p1 >;
#declare v47 = < sqr2p1 , sqr2p1 , 1.0 , sqr2p1 >;
#declare v48 = < -sqr2p1 , -sqr2p1 , -sqr2p1 , -1.0 >;
#declare v49 = < -sqr2p1 , -sqr2p1 , -sqr2p1 , 1.0 >;
#declare v50 = < -sqr2p1 , -sqr2p1 , sqr2p1 , -1.0 >;
#declare v51 = < -sqr2p1 , -sqr2p1 , sqr2p1 , 1.0 >;
#declare v52 = < -sqr2p1 , sqr2p1 , -sqr2p1 , -1.0 >;
#declare v53 = < -sqr2p1 , sqr2p1 , -sqr2p1 , 1.0 >;
#declare v54 = < -sqr2p1 , sqr2p1 , sqr2p1 , -1.0 >;
#declare v55 = < -sqr2p1 , sqr2p1 , sqr2p1 , 1.0 >;
#declare v56 = < sqr2p1 , -sqr2p1 , -sqr2p1 , -1.0 >;
#declare v57 = < sqr2p1 , -sqr2p1 , -sqr2p1 , 1.0 >;
#declare v58 = < sqr2p1 , -sqr2p1 , sqr2p1 , -1.0 >;
#declare v59 = < sqr2p1 , -sqr2p1 , sqr2p1 , 1.0 >;
#declare v60 = < sqr2p1 , sqr2p1 , -sqr2p1 , -1.0 >;
#declare v61 = < sqr2p1 , sqr2p1 , -sqr2p1 , 1.0 >;
#declare v62 = < sqr2p1 , sqr2p1 , sqr2p1 , -1.0 >;
#declare v63 = < sqr2p1 , sqr2p1 , sqr2p1 , 1.0 >;
#declare vertices = array[64] 
    {v0,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,
    v16,v17,v18,v19,v20,v21,v22,v23,v24,v25,v26,v27,v28,v29,v30,v31,
    v32,v33,v34,v35,v36,v37,v38,v39,v40,v41,v42,v43,v44,v45,v46,v47,
    v48,v49,v50,v51,v52,v53,v54,v55,v56,v57,v58,v59,v60,v61,v62,v63};

/* edges */
#declare edges = array[128][2] 
{ { 0 , 8 }
, { 0 , 16 }
, { 0 , 32 }
, { 0 , 48 }
, { 1 , 9 }
, { 1 , 17 }
, { 1 , 33 }
, { 1 , 49 }
, { 2 , 10 }
, { 2 , 18 }
, { 2 , 34 }
, { 2 , 50 }
, { 3 , 11 }
, { 3 , 19 }
, { 3 , 35 }
, { 3 , 51 }
, { 4 , 12 }
, { 4 , 20 }
, { 4 , 36 }
, { 4 , 52 }
, { 5 , 13 }
, { 5 , 21 }
, { 5 , 37 }
, { 5 , 53 }
, { 6 , 14 }
, { 6 , 22 }
, { 6 , 38 }
, { 6 , 54 }
, { 7 , 15 }
, { 7 , 23 }
, { 7 , 39 }
, { 7 , 55 }
, { 8 , 24 }
, { 8 , 40 }
, { 8 , 56 }
, { 9 , 25 }
, { 9 , 41 }
, { 9 , 57 }
, { 10 , 26 }
, { 10 , 42 }
, { 10 , 58 }
, { 11 , 27 }
, { 11 , 43 }
, { 11 , 59 }
, { 12 , 28 }
, { 12 , 44 }
, { 12 , 60 }
, { 13 , 29 }
, { 13 , 45 }
, { 13 , 61 }
, { 14 , 30 }
, { 14 , 46 }
, { 14 , 62 }
, { 15 , 31 }
, { 15 , 47 }
, { 15 , 63 }
, { 16 , 20 }
, { 16 , 32 }
, { 16 , 48 }
, { 17 , 21 }
, { 17 , 33 }
, { 17 , 49 }
, { 18 , 22 }
, { 18 , 34 }
, { 18 , 50 }
, { 19 , 23 }
, { 19 , 35 }
, { 19 , 51 }
, { 20 , 36 }
, { 20 , 52 }
, { 21 , 37 }
, { 21 , 53 }
, { 22 , 38 }
, { 22 , 54 }
, { 23 , 39 }
, { 23 , 55 }
, { 24 , 28 }
, { 24 , 40 }
, { 24 , 56 }
, { 25 , 29 }
, { 25 , 41 }
, { 25 , 57 }
, { 26 , 30 }
, { 26 , 42 }
, { 26 , 58 }
, { 27 , 31 }
, { 27 , 43 }
, { 27 , 59 }
, { 28 , 44 }
, { 28 , 60 }
, { 29 , 45 }
, { 29 , 61 }
, { 30 , 46 }
, { 30 , 62 }
, { 31 , 47 }
, { 31 , 63 }
, { 32 , 34 }
, { 32 , 48 }
, { 33 , 35 }
, { 33 , 49 }
, { 34 , 50 }
, { 35 , 51 }
, { 36 , 38 }
, { 36 , 52 }
, { 37 , 39 }
, { 37 , 53 }
, { 38 , 54 }
, { 39 , 55 }
, { 40 , 42 }
, { 40 , 56 }
, { 41 , 43 }
, { 41 , 57 }
, { 42 , 58 }
, { 43 , 59 }
, { 44 , 46 }
, { 44 , 60 }
, { 45 , 47 }
, { 45 , 61 }
, { 46 , 62 }
, { 47 , 63 }
, { 48 , 49 }
, { 50 , 51 }
, { 52 , 53 }
, { 54 , 55 }
, { 56 , 57 }
, { 58 , 59 }
, { 60 , 61 }
, { 62 , 63 } };

/* tetrahedra */
#declare tetrahedra = array[16][4]
{ { 0 , 16 , 32 , 48 }
, { 11 , 27 , 43 , 59 }
, { 12 , 28 , 44 , 60 }
, { 8 , 24 , 40 , 56 }
, { 9 , 25 , 41 , 57 }
, { 15 , 31 , 47 , 63 }
, { 13 , 29 , 45 , 61 }
, { 14 , 30 , 46 , 62 }
, { 10 , 26 , 42 , 58 }
, { 3 , 19 , 35 , 51 }
, { 2 , 18 , 34 , 50 }
, { 1 , 17 , 33 , 49 }
, { 4 , 20 , 36 , 52 }
, { 5 , 21 , 37 , 53 }
, { 6 , 22 , 38 , 54 }
, { 7 , 23 , 39 , 55 } };

#declare faces = array[64][3];
#for(i, 0, 15)
    #declare faces[4*i][0] = tetrahedra[i][0];
    #declare faces[4*i][1] = tetrahedra[i][1];
    #declare faces[4*i][2] = tetrahedra[i][2];
    #declare faces[4*i+1][0] = tetrahedra[i][0];
    #declare faces[4*i+1][1] = tetrahedra[i][1];
    #declare faces[4*i+1][2] = tetrahedra[i][3];
    #declare faces[4*i+2][0] = tetrahedra[i][0];
    #declare faces[4*i+2][1] = tetrahedra[i][2];
    #declare faces[4*i+2][2] = tetrahedra[i][3];
    #declare faces[4*i+3][0] = tetrahedra[i][1];
    #declare faces[4*i+3][1] = tetrahedra[i][2];
    #declare faces[4*i+3][2] = tetrahedra[i][3];
#end


/* rotation in 4D space */
#macro rotate4d(theta, phi, xi, vec)
    #local a = cos(xi);
    #local b = sin(theta)*cos(phi)*sin(xi);
    #local c = sin(theta)*sin(phi)*sin(xi);
    #local d = cos(theta)*sin(xi);
    #local p = vec.x;
    #local q = vec.y;
    #local r = vec.z;
    #local s = vec.t;
    < a*p - b*q - c*r - d*s
    , a*q + b*p + c*s - d*r
    , a*r - b*s + c*p + d*q
    , a*s + b*r - c*q + d*p >
#end

/* stereographic projection */
#macro StereographicProjection(q, r)
    <q.x,q.y,q.z>/(r-q.t)
#end

/* rotated and projected vertices */
#macro ProjectedVertices(theta, phi, xi, r)
    #local out = array[64];
    #for(i, 0, 63)
        #local out[i] = StereographicProjection(
                            rotate4d(theta,phi,xi,vertices[i]), r
                        );
    #end
    out
#end

/* macro spherical segment */
#macro vlength4(P)
    sqrt(P.x*P.x + P.y*P.y + P.z*P.z + P.t*P.t)
#end

#macro sphericalSegment(P, Q, r, n)
    #local out = array[n+1];
    #for(i, 0, n)
        #local pt = P + (i/n)*(Q-P);
        #local out[i] = r/vlength4(pt) * pt;
    #end
    out
#end

/* macro to draw an edge */
#macro Edge(verts, v1, v2, theta, phi, xi, r, Tex)
    #local P = verts[v1];
    #local Q = verts[v2];
    #local PQ = sphericalSegment(P, Q, r, 100);
    sphere_sweep {
        b_spline 101
        #for(k,0,100)
            #local O = 
                StereographicProjection(rotate4d(theta,phi,xi,PQ[k]), r);
            O log(1+vlength(O)/4)/2
        #end
        texture { 
            Tex
        }
    }
#end

/* stereographic subdivision (to fill the triangular faces) */
#macro midpoint4(A, B, r)
    #local xmid = (A.x + B.x)/2;
    #local ymid = (A.y + B.y)/2;
    #local zmid = (A.z + B.z)/2;
    #local tmid = (A.t + B.t)/2;
    #local lg = sqrt(xmid*xmid + ymid*ymid + zmid*zmid + tmid*tmid) / r;
    < xmid / lg, ymid / lg, zmid / lg, tmid / lg >
#end

#macro subdiv0(A, B, C, r)
    #local mAB = midpoint4(A, B, r);
    #local mAC = midpoint4(A, C, r);
    #local mBC = midpoint4(B, C, r);
    #local trgl1 = array[3] {A, mAB, mAC};
    #local trgl2 = array[3] {B, mAB, mBC};
    #local trgl3 = array[3] {C, mBC, mAC};
    #local trgl4 = array[3] {mAB, mAC, mBC};
    array[4] {trgl1, trgl2, trgl3, trgl4}
#end

#macro subdiv(A, B, C, r, depth)
    #if(depth=1)
        #local out = subdiv0(A, B, C, r);
    #else
        #local triangles = subdiv(A, B, C, r, depth-1);
        #local out = array[pow(4,depth)];
        #for(i, 0, pow(4,depth-1)-1)
            #local trgl = triangles[i];
            #local trgls = subdiv0(trgl[0], trgl[1], trgl[2], r);
            #local out[4*i] = trgls[0];
            #local out[4*i+1] = trgls[1];
            #local out[4*i+2] = trgls[2];
            #local out[4*i+3] = trgls[3];
        #end
    #end
    out
#end

/*-------------------------------------------*/
/*-----      draw the polychoron       ------*/
/*-------------------------------------------*/
#declare theta = pi/2;
#declare phi = 0;
#declare xi = 2*frame_number*pi/180;
#declare r = sqrt(1+3*sqr2p1*sqr2p1); 
#declare depth = 5;

#declare vs = ProjectedVertices(theta, phi, xi, r);

#declare stereoTriangles = array[64];
#for(i, 0, 63)
    #local triangles4 = 
        subdiv(
            vertices[faces[i][0]], 
            vertices[faces[i][1]], 
            vertices[faces[i][2]], r, depth);
    #declare stereoTriangles[i] = array[dimension_size(triangles4,1)][3];
    #for(j, 0, dimension_size(triangles4,1)-1)
        #local trgl4 = triangles4[j];
        #declare stereoTriangles[i][j][0] = 
            StereographicProjection(rotate4d(theta, phi, xi, trgl4[0]), r);
        #declare stereoTriangles[i][j][1] = 
            StereographicProjection(rotate4d(theta, phi, xi, trgl4[1]), r);
        #declare stereoTriangles[i][j][2] = 
            StereographicProjection(rotate4d(theta, phi, xi, trgl4[2]), r);
    #end
#end

#declare edgeTexture = texture {
    pigment { color DarkPurple }
    finish {
      ambient .1
      diffuse .9
      reflection 0
      specular 1
      metallic
    }
};

object {
    union {
        /* draw edges */
        #for(i, 0, 127)
            Edge(vertices, edges[i][0], edges[i][1], 
                    theta, phi, xi, r, edgeTexture)
        #end
        /* draw vertices */ 
        #for(i,0,63)
            sphere {
                vs[i], vlength(vs[i])/15
                texture { edgeTexture }
            }
        #end  
        /* fill triangles */
        mesh {
            #for(i, 0, 63)
                #local trgl = stereoTriangles[i];
                #for(j, 0, dimension_size(trgl,1)-1)
                    triangle {
                        trgl[j][0], trgl[j][1], trgl[j][2]
                    }
                #end
            #end
            texture {
                pigment { Red }
                finish {
                    ambient 0.5
                    reflection 0
                    brilliance 4
                }
            }
        }
     }
     scale 0.5
     translate <-8, 6, -25>
}

To leave a comment for the author, please follow the link and comment on their blog: Saturn Elephant.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.