library(rgl)

# vertices ####
phi <- (1+sqrt(5))/2
a <- 1/sqrt(3)
b <- a/phi
c <- a*phi
vertices <-
  rbind(
    c( a,  a,  a),
    c( a,  a, -a),
    c( a, -a,  a),
    c(-a, -a,  a),
    c(-a,  a, -a),
    c(-a,  a,  a),
    c( 0,  b, -c),
    c( 0, -b, -c),
    c( 0, -b,  c),
    c( c,  0, -b),
    c(-c,  0, -b),
    c(-c,  0,  b),
    c( b,  c,  0),
    c( b, -c,  0),
    c(-b, -c,  0),
    c(-b,  c,  0),
    c( 0,  b,  c),
    c( a, -a, -a),
    c( c,  0,  b),
    c(-a, -a, -a)
  )

# tetrahedra vertices indices ####
alltetraIdxs <-
  rbind(
    c( 1,13,16),c(13,10,16),c(13, 1,10),c( 1,16,10),
    c( 3, 0,17),c( 0, 4,17),c( 3,17, 4),c( 0, 3, 4),
    c(14, 5,18),c( 5, 6,18),c( 6,14,18),c( 5,14, 6),
    c(11,12, 2),c(12, 7, 2),c(12,11, 7),c( 2, 7,11),
    c( 9,15,19),c(15, 8,19),c(15, 9, 8),c( 8, 9,19),
    c( 8,12,17),c( 8,10,12),c(10,17,12),c(10, 8,17),
    c( 1, 2, 5),c( 1,19, 2),c(19, 1, 5),c(19, 5, 2),
    c(18, 3,15),c(18, 7, 3),c(18,15, 7),c( 7,15, 3),
    c( 0,13,11),c( 0, 6,13),c( 6,11,13),c(11, 6, 0),
    c( 4,14, 9),c( 4,16,14),c(16, 9,14),c( 4, 9,16)
  ) + 1
tetraIdxs <- lapply(split(1:40, gl(10,4)),
                    function(is) alltetraIdxs[is,])

# tetrahedra vertices ####
tetrahedra <- lapply(tetraIdxs,
                     function(t) plyr::alply(t, 1,
                                             function(idxs) vertices[idxs,]))

# rgl - compound of tetrahedra ####
colors <- rainbow(10)
lapply(1:10,
       function(i) lapply(tetrahedra[[i]],
                          function(face) triangles3d(face, color=colors[i])))


# edges vertices indices ####
alledgesIdxs <-
  rbind(
    c( 1,13),c(13,16),c(16, 1),c(13,10),c(10,16),c( 1,10),
    c( 3, 0),c( 0,17),c(17, 3),c( 0, 4),c( 4,17),c( 4, 3),
    c(14, 5),c( 5,18),c(18,14),c( 5, 6),c( 6,18),c( 6,14),
    c(11,12),c(12, 2),c( 2,11),c(12, 7),c( 7, 2),c(11, 7),
    c( 9,15),c(15,19),c(19, 9),c(15, 8),c( 8,19),c( 9, 8),
    c( 8,12),c(12,17),c(17, 8),c( 8,10),c(10,12),c(10,17),
    c( 1, 2),c( 2, 5),c( 5, 1),c( 1,19),c(19, 2),c( 5,19),
    c(18, 3),c( 3,15),c(15,18),c(18, 7),c( 7, 3),c(15, 7),
    c( 0,13),c(13,11),c(11, 0),c( 0, 6),c( 6,13),c( 6,11),
    c( 4,14),c(14, 9),c( 9, 4),c( 4,16),c(16,14),c(16, 9)
  ) + 1
edgesIdxs <- lapply(split(1:60, gl(10,6)),
                    function(is) alledgesIdxs[is,])

# edges vertices ####
edges <- lapply(edgesIdxs,
                function(vs) plyr::alply(vs, 1,
                                         function(idxs) vertices[idxs,]))

# vertices per tetrahedron ####
verticesIdxs <- lapply(edgesIdxs, function(idxs) unique(c(idxs)))
vs <- lapply(verticesIdxs, function(idxs) vertices[idxs,])

# rgl - interlocked tetrahedra ####
nicesegment3d <- function(p1p2, radius=0.05, col="black"){
  shade3d(cylinder3d(p1p2, radius=radius, sides=30), color=col)
}
# edges
lapply(1:10,
       function(i) lapply(edges[[i]],
                          function(v1v2) nicesegment3d(v1v2, col=colors[i])))
# vertices
lapply(1:10,
       function(i) spheres3d(vs[[i]], radius=0.05, col=colors[i]))