Site icon R-bloggers

Drawing a stereographic duoprism

[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.

In this post, I’ll show how to draw a stereographic duoprism using R, Asymptote and POV-Ray.

With R

library(rgl)

A <- 8L # number of sides of the first polygon
B <- 4L # number of sides of the second polygon

# construction of the vertices
vertices <- array(NA_real_, dim = c(A,B,4L))
for(i in 1L:A){
  v1 <- c(cos(i/A*2*pi), sin(i/A*2*pi))
  for(j in 1L:B){
    v2 <- c(cos(j/B*2*pi), sin(j/B*2*pi))
    vertices[i,j,] <- c(v1,v2)
  }
}

# construction of the edges
edges <- array(NA_integer_, dim = c(2L,2L,2L*A*B))
dominates <- function(c1, c2){
  c2[1L]>c1[1L] || (c2[1L]==c1[1L] && c2[2L]>c1[2L])
}
counter <- 1L
for(i in seq_len(A)-1L){
  for(j in seq_len(B)-1L){
    c1 <- c(i,j)
    candidate <- c(i, (j-1L)%%B)
    if(dominates(c1, candidate)){
      edges[,,counter] <- cbind(c1, candidate) + 1L
      counter <- counter + 1L
    }
    candidate <- c(i, (j+1L)%%B)
    if(dominates(c1, candidate)){
      edges[,,counter] <- cbind(c1, candidate) + 1L
      counter <- counter + 1L
    }
    candidate <- c((i-1L)%%A, j)
    if(dominates(c1, candidate)){
      edges[,,counter] <- cbind(c1, candidate) + 1L
      counter <- counter + 1L
    }
    candidate <- c((i+1L)%%A, j)
    if(dominates(c1, candidate)){
      edges[,,counter] <- cbind(c1, candidate) + 1L
      counter <- counter + 1L
    }
  }
}

# stereographic projection
stereog <- function(v){
  v[1L:3L] / (sqrt(2) - v[4L])
}

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

# stereographic edge
stereoEdge <- function(verts, v1, v2){
  P <- verts[v1[1L], v1[2L], ]
  Q <- verts[v2[1L], v2[2L], ]
  PQ <- sphericalSegment(P, Q, 100L)
  pq <- t(apply(PQ, 1L, stereog))
  dists <- sqrt(apply(pq, 1L, crossprod))
  cylinder3d(pq, radius = dists/15, sides = 60)
}

# projected vertices
vs <- apply(vertices, c(1L,2L), stereog)

####~~~~ 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:(2L*A*B)){
  v1 <- edges[, 1L, k]
  v2 <- edges[, 2L, k]
  edge <- stereoEdge(vertices, v1, v2)
  shade3d(edge, color = "gold")
}
## plot the vertices
for(i in 1L:A){
  for(j in 1L:B){
    v <- vs[,i,j]
    spheres3d(v, radius = sqrt(c(crossprod(v)))/10 , color = "gold2")
  }
}

With Asymptote

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

currentprojection = orthographic(4,4,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 = {
    "DP000", "DP001", "DP002", "DP003", "DP004", "DP005",
    "DP006", "DP007", "DP008", "DP009", "DP010", "DP011",
    "DP012", "DP013", "DP014", "DP015", "DP016", "DP017",
    "DP018", "DP019", "DP020", "DP021", "DP022", "DP023",
    "DP024", "DP025", "DP026", "DP027", "DP028", "DP029",
    "DP030", "DP031", "DP032", "DP033", "DP034", "DP035",
    "DP036", "DP037", "DP038", "DP039", "DP040", "DP041",
    "DP042", "DP043", "DP044", "DP045", "DP046", "DP047",
    "DP048", "DP049", "DP050", "DP051", "DP052", "DP053",
    "DP054", "DP055", "DP056", "DP057", "DP058", "DP059",
    "DP060", "DP061", "DP062", "DP063", "DP064", "DP065",
    "DP066", "DP067", "DP068", "DP069", "DP070", "DP071",
    "DP072", "DP073", "DP074", "DP075", "DP076", "DP077",
    "DP078", "DP079", "DP080", "DP081", "DP082", "DP083",
    "DP084", "DP085", "DP086", "DP087", "DP088", "DP089",
    "DP090", "DP091", "DP092", "DP093", "DP094", "DP095",
    "DP096", "DP097", "DP098", "DP099", "DP100", "DP101",
    "DP102", "DP103", "DP104", "DP105", "DP106", "DP107",
    "DP108", "DP109", "DP110", "DP111", "DP112", "DP113",
    "DP114", "DP115", "DP116", "DP117", "DP118", "DP119",
    "DP120", "DP121", "DP122", "DP123", "DP124", "DP125",
    "DP126", "DP127", "DP128", "DP129", "DP130", "DP131",
    "DP132", "DP133", "DP134", "DP135", "DP136", "DP137",
    "DP138", "DP139", "DP140", "DP141", "DP142", "DP143",
    "DP144", "DP145", "DP146", "DP147", "DP148", "DP149",
    "DP150", "DP151", "DP152", "DP153", "DP154", "DP155",
    "DP156", "DP157", "DP158", "DP159", "DP160", "DP161",
    "DP162", "DP163", "DP164", "DP165", "DP166", "DP167",
    "DP168", "DP169", "DP170", "DP171", "DP172", "DP173",
    "DP174", "DP175", "DP176", "DP177", "DP178", "DP179"};

// lexicographic order ---------------------------------------------------------
bool dominates(int[] e1, int[] e2){
    return e2[0]>e1[0] || (e2[0]==e1[0] && e2[1]>e1[1]);
}

// vertices --------------------------------------------------------------------
int A = 8;
int B = 4;

struct quadruple {
    real x;
    real y;
    real z;
    real t;
}

real[][] poly1 = new real[A][2];
for(int i = 0; i < A; ++i){
    poly1[i][0] = cos(i/A*2pi);
    poly1[i][1] = sin(i/A*2pi);
}
real[][] poly2 = new real[B][2];
for(int i = 0; i < B; ++i){
    poly2[i][0] = cos(pi/B+i/B*2pi);
    poly2[i][1] = sin(pi/B+i/B*2pi);
}
quadruple[][] vertices = new quadruple[A][B];
for(int i = 0; i < A; ++i){
    for(int j = 0; j < B; ++j){
        quadruple v;
        v.x = poly1[i][0]; v.y = poly1[i][1]; 
        v.z = poly2[j][0]; v.t = poly2[j][1];
        vertices[i][j] = v;
    }
}

// edges -----------------------------------------------------------------------
int[][][] edges;
for(int i = 0; i < A; ++i){
    for(int j = 0; j < B; ++j){
        int[] e = {i,j};
        int[] candidate = {i,(j-1)%B};
        if(dominates(e,candidate)){
            int[][] edge = {e,candidate}; 
            edges.push(edge);
        }
        int[] candidate = {i,(j+1)%B};
        if(dominates(e,candidate)){
            int[][] edge = {e,candidate}; 
            edges.push(edge);
        }
        int[] candidate = {(i-1)%A,j};
        if(dominates(e,candidate)){
            int[][] edge = {e,candidate}; 
            edges.push(edge);
        }
        int[] candidate = {(i+1)%A,j};
        if(dominates(e,candidate)){
            int[][] edge = {e,candidate}; 
            edges.push(edge);
        }
    }
}

// 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); 
}

// bounding box ----------------------------------------------------------------
real f=3, h = 4.5, g = 1.5;
path3 boundingbox = (-h,0,-f)--(-h,0,g)--(h,0,f)--(h,0,-g)--cycle;

// draw the duoprism -----------------------------------------------------------
int n = 100;
real r = sqrt(2);
real alpha = pi/2, beta = 0;

for(int file = 0; file < 180; ++file){
    real xi = 2*file*pi/180;
    picture pic;
    // draw bounding box
    draw(pic, boundingbox, rgb("363940ff")+opacity(0));
    // draw edges
    for(int k = 0; k < 2*A*B; ++k){
        quadruple A = vertices[edges[k][0][0]][edges[k][0][1]];
        quadruple B = vertices[edges[k][1][0]][edges[k][1][1]];
        path3 p3 = 
            stereoPath(rotate4d(alpha, beta, xi, A), 
                       rotate4d(alpha, beta, xi, 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 < A; ++i){
        for(int j = 0; j < B; ++j){
            triple Asg = 
                stereog(rotate4d(alpha, beta, xi, vertices[i][j]), r);
            draw(pic, shift(Asg)*scale3(length(Asg)/10)*unitsphere, purple);
        }
    }
    // add and save picture
    add(pic);
    shipout(files[file], bbox(rgb("363940ff"), FillDraw(rgb("363940ff"))));
    erase();
}

/* to do the animation
gs -dSAFER -dBATCH -dNOPAUSE -dEPSCrop -sDEVICE=png16m -r600 -sOutputFile=zpic%03d.png DP*.eps
mogrify -resize 512x zpic*.png
gifski --fps 12 zpic*.png -o DuoprismStereo.gif
*/

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
  right x*image_width/image_height
}

// sun -------------------------------------------------------------------------
light_source {< 4000,6000,-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 A = 4;
#declare B = 30;

#declare poly1 = array[A];
#for(i,0,A-1)
  #declare poly1[i] = array[2] {cos(i/A*2*pi), sin(i/A*2*pi)};
#end
#declare poly2 = array[B];
#for(i,0,B-1)
  #declare poly2[i] = array[2] {cos(i/B*2*pi), sin(i/B*2*pi)};
#end
#declare vertices = array[A][B];
#for(i,0,A-1)
  #for(j,0,B-1)
    #declare vertices[i][j] = 
      < poly1[i][0], poly1[i][1], poly2[j][0], poly2[j][1] >;
  #end
#end

/* ----- edges ----- */
#macro dominates(e1,e2)
  (e2[0]>e1[0]) | ((e2[0]=e1[0]) & (e2[1]>e1[1]))
#end
#declare nedges = 2*A*B;
#declare edges = array[nedges];
#declare k=0;
#for(i,0,A-1)
  #for(j,0,B-1)
    #local e = array[2] {i,j};
    #local candidate = array[2] {i,mod(mod(j-1,B)+B,B)};
    #if(dominates(e,candidate))
      #local edge = array[2] {e,candidate};
      #declare edges[k] = edge;
      #declare k = k+1;
    #end
    #local candidate = array[2] {i,mod(mod(j+1,B)+B,B)};
    #if(dominates(e,candidate))
      #local edge = array[2] {e,candidate};
      #declare edges[k] = edge;
      #declare k = k+1;
    #end
    #local candidate = array[2] {mod(mod(i-1,A)+A,A),j};
    #if(dominates(e,candidate))
      #local edge = array[2] {e,candidate};
      #declare edges[k] = edge;
      #declare k = k+1;
    #end
    #local candidate = array[2] {mod(mod(i+1,A)+A,A),j};
    #if(dominates(e,candidate))
      #local edge = array[2] {e,candidate};
      #declare edges[k] = edge;
      #declare k = k+1;
    #end
  #end
#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)
  <q.x,q.y,q.z> / (sqrt(2)-q.t)
#end

/* rotated and projected vertices */
#macro ProjectedVertices(theta,phi,xi)
  #local out = array[A][B];
  #for(i,0,A-1)
    #for(j,0,B-1)
      #local out[i][j] = StereographicProjection(
                            rotate4d(theta,phi,xi,vertices[i][j])
                          );
    #end
  #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, n)
  #local out = array[n+1];
  #for(i, 0, n)
    #local pt = P + (i/n)*(Q-P);
    #local out[i] = sqrt(2)/vlength4(pt) * pt;
  #end
  out
#end

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

/*-----------------------------------------*/
/*-----      draw the duoprism       ------*/
/*-----------------------------------------*/
#declare theta = pi/2;
#declare phi = 0;
#declare xi = 2*frame_number*pi/180; 

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

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

object {
  union {
    /* draw edges */
    #for(i, 0, 2*A*B-1)
      Edge(vertices, edges[i][0], edges[i][1], 
                    theta, phi, xi, edgeTexture)
    #end
    /* draw vertices */
    #for(i,0,A-1)
      #for(j,0,B-1)
        sphere {
          vs[i][j], vlength(vs[i][j])/10
          texture { edgeTexture }
        }
      #end
    #end  
  }
  translate <-3, 6, -15>
  scale 0.8
}

/* ini file 
Width = 512
Height = 512
Antialias = On
Antialias_Threshold = 0.3
Input_File_Name = DuoprismStereographic.pov
Initial_Clock = 0
Final_Clock = 1
Initial_Frame = 0
Final_Frame = 179
Subset_Start_Frame = 0
Cyclic_Animation = on
*/


Here is another one. This is a hexagonal duoprism with a cell colored in red.

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

// camera ----------------------------------------------------------------------
camera {
  location <0, 0,-10>
  look_at 0
  angle 45
  right x*image_width/image_height
}

// light sources ---------------------------------------------------------------
light_source { <0,0,-100> White shadowless } 
light_source { <100,0,-100> White shadowless } 

// moon ------------------------------------------------------------------------
light_source { 
  <-1000, 800, 3000> 
  color White
  shadowless
  looks_like { 
    sphere { 
      <0,0,0>, 300 
      texture { 
        pigment { 
          color Yellow
        }
        normal { 
          bumps 0.5
          scale 50
        }
        finish { 
          emission 0.8   
          diffuse 0.2
          phong 1
        }                                        
      } 
    } 
  } 
} 

// sky -------------------------------------------------------------------------
plane { 
  <0,1,0>, 1 hollow  
  texture { 
    pigment { 
      color rgb <0.01, 0.01, 0.2> 
    }
    finish { 
      emission 0.5 
      diffuse 0.5 
    } 
  }  
  scale 10000
}

// the clouds ------------------------------------------------------------------
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.00 rgb <0.5,0.5,0.5>        ]
      }
      scale 2.5
      translate <0,1,0>
    }
    finish { 
      emission 0.25 
      diffuse 0
    } 
  }      
  scale 5000
}

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

// sea -------------------------------------------------------------------------
plane { 
  <0,1,0>, -1 hollow
  texture{
    pigment{
      rgb <.098,.098,.439>
    }
    finish {
      ambient 0.15
      diffuse 0.55
      brilliance 6.0
      phong 0.8
      phong_size 120
      reflection 0.2
    }
    normal { 
      bumps 0.95
      turbulence .05 
      scale <1,0.25,1> 
    }
  }
}

// vertices --------------------------------------------------------------------
#declare a = sqrt(3) / 2;
#declare vertices = array[36]
  {
    <a, 0.5, a, 0.5>,
    <a, 0.5, 0.0, 1.0>,
    <a, 0.5, -a, 0.5>,
    <a, 0.5, -a, -0.5>,
    <a, 0.5, 0.0, -1.0>,
    <a, 0.5, a, -0.5>,
    <0.0, 1.0, a, 0.5>,
    <0.0, 1.0, 0.0, 1.0>,
    <0.0, 1.0, -a, 0.5>,
    <0.0, 1.0, -a, -0.5>,
    <0.0, 1.0, 0.0, -1.0>,
    <0.0, 1.0, a, -0.5>,
    <-a, 0.5, a, 0.5>,
    <-a, 0.5, 0.0, 1.0>,
    <-a, 0.5, -a, 0.5>,
    <-a, 0.5, -a, -0.5>,
    <-a, 0.5, 0.0, -1.0>,
    <-a, 0.5, a, -0.5>,
    <-a, -0.5, a, 0.5>,
    <-a, -0.5, 0.0, 1.0>,
    <-a, -0.5, -a, 0.5>,
    <-a, -0.5, -a, -0.5>,
    <-a, -0.5, 0.0, -1.0>,
    <-a, -0.5, a, -0.5>,
    <0.0, -1.0, a, 0.5>,
    <0.0, -1.0, 0.0, 1.0>,
    <0.0, -1.0, -a, 0.5>,
    <0.0, -1.0, -a, -0.5>,
    <0.0, -1.0, 0.0, -1.0>,
    <0.0, -1.0, a, -0.5>,
    <a, -0.5, a, 0.5>,
    <a, -0.5, 0.0, 1.0>,
    <a, -0.5, -a, 0.5>,
    <a, -0.5, -a, -0.5>,
    <a, -0.5, 0.0, -1.0>,
    <a, -0.5, a, -0.5>
  };
#declare facetVertices = array[12] {0,5,6,30,11,35,12,17,18,23,24,29};
#declare otherVertices = array[24] 
  {1,2,3,4,7,8,
  9,10,13,14,15,16,
  19,20,21,22,25,26,
  27,28,31,32,33,34};

// edges -------------------------------------------------------------------
#declare facetEdges = array[18][2]
  {
    {0, 5},
    {0, 6},
    {0, 30},
    {5, 11},
    {5, 35},
    {6, 11},
    {6, 12},
    {11, 17},
    {12, 17},
    {12, 18},
    {17, 23},
    {18, 23},
    {18, 24},
    {23, 29},
    {24, 29},
    {24, 30},
    {29, 35},
    {30, 35}
  };
#declare otherEdges = array[54][2]
  {
    {0, 1},
    {1, 2},
    {1, 7},
    {1, 31},
    {2, 3},
    {2, 8},
    {2, 32},
    {3, 4},
    {3, 9},
    {3, 33},
    {4, 5},
    {4, 10},
    {4, 34},
    {6, 7},
    {7, 8},
    {7, 13},
    {8, 9},
    {8, 14},
    {9, 10},
    {9, 15},
    {10, 11},
    {10, 16},
    {12, 13},
    {13, 14},
    {13, 19},
    {14, 15},
    {14, 20},
    {15, 16},
    {15, 21},
    {16, 17},
    {16, 22},
    {18, 19},
    {19, 20},
    {19, 25},
    {20, 21},
    {20, 26},
    {21, 22},
    {21, 27},
    {22, 23},
    {22, 28},
    {24, 25},
    {25, 26},
    {25, 31},
    {26, 27},
    {26, 32},
    {27, 28},
    {27, 33},
    {28, 29},
    {28, 34},
    {30, 31},
    {31, 32},
    {32, 33},
    {33, 34},
    {34, 35}
  };

// macros ----------------------------------------------------------------------
#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, n)
  #local out = array[n+1];
  #for(i, 0, n)
    #local pt = P + (i/n)*(Q-P);
    #local out[i] = sqrt(2)/vlength4(pt) * pt;
  #end
  out
#end

#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

#macro StereographicProjection(q) 
  acos(q.t/sqrt(2))/sqrt(2-q.t*q.t) * <q.x,q.y,q.z>
#end

#macro ProjectedFacetVertices(theta, phi, xi)
  #local out = array[12];
  #for(i, 0, 11)
    #local out[i] = 
      StereographicProjection(
        rotate4d(theta, phi, xi, vertices[facetVertices[i]])
      );
  #end
  out
#end

#macro ProjectedOtherVertices(theta, phi, xi)
  #local out = array[24];
  #for(i, 0, 23)
    #local out[i] = 
      StereographicProjection(
        rotate4d(theta, phi, xi, vertices[otherVertices[i]])
      );
  #end
  out
#end

// texture ---------------------------------------------------------------------
#declare edgeTexture1 = 
  texture {
    New_Penny
    finish {
      ambient 0.01
      diffuse 2
      reflection 0
      brilliance 8
      specular 0.1
      roughness 0.1
    }
  };

#declare edgeTexture2 = 
  texture {
    pigment { Red }
    finish {
      ambient 0.01
      diffuse 2
      reflection 0
      brilliance 8
      specular 0.1
      roughness 0.1
    }
  };

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

// draw ------------------------------------------------------------------------
#declare theta = pi/2;
#declare phi = 0;
#declare xi = 2*frame_number*pi/180;
#declare vsFacet = ProjectedFacetVertices(theta, phi, xi);
#declare vsOther = ProjectedOtherVertices(theta, phi, xi);
object {
  union {
    #for(i, 0, 53)
      Edge(vertices, otherEdges[i][0], otherEdges[i][1], 
                            theta, phi, xi, edgeTexture1)
    #end
    #for(i, 0, 17)
      Edge(vertices, facetEdges[i][0], facetEdges[i][1], 
                            theta, phi, xi, edgeTexture2)
    #end
    #for(i, 0, 23)
      sphere {
        vsOther[i], vlength(vsOther[i])/10
        texture { edgeTexture1 }
      }
    #end  
    #for(i, 0, 11)
      sphere {
        vsFacet[i], vlength(vsFacet[i])/10
        texture { edgeTexture2 }
      }
    #end  
  }
  scale 0.5
  rotate <60, 0, 0>
  translate <0, 0.5, -2>
}

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.