--- title: "Drawing a duoprism" author: "Stéphane Laurent" date: '2019-07-08' tags: R, maths, geometry, graphics, rgl, asymptote, povray output: md_document: variant: markdown preserve_yaml: true html_document: highlight: kate keep_md: no prettify: yes linenums: yes prettifycss: minimal highlighter: pandoc-solarized --- In geometry of 4 dimensions, a *duoprism* is a polytope resulting from the Cartesian product of two polygons. We show how to draw a duoprism using R, Asymptote and POV-Ray, when the two polygons are regular. # Drawing a duoprism with R ```r A <- 15 B <- 23 # construction of the vertices vertices <- array(NA_real_, dim = c(A,B,4)) for(i in 1:A){ v1 <- c(cos(i/A*2*pi), sin(i/A*2*pi)) for(j in 1: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(2,2,2*A*B)) dominates <- function(c1, c2){ c2[1]>c1[1] || (c2[1]==c1[1] && c2[2]>c1[2]) } counter <- 1 for(i in seq_len(A)-1L){ for(j in seq_len(B)-1L){ c1 <- c(i,j) candidate <- c(i, (j-1)%%B) if(dominates(c1, candidate)){ edges[,,counter] <- cbind(c1, candidate) + 1L counter <- counter + 1 } candidate <- c(i, (j+1)%%B) if(dominates(c1, candidate)){ edges[,,counter] <- cbind(c1, candidate) + 1L counter <- counter + 1 } candidate <- c((i-1)%%A, j) if(dominates(c1, candidate)){ edges[,,counter] <- cbind(c1, candidate) + 1L counter <- counter + 1 } candidate <- c((i+1)%%A, j) if(dominates(c1, candidate)){ edges[,,counter] <- cbind(c1, candidate) + 1L counter <- counter + 1 } } } # stereographic projection stereog <- function(v){ v[1:3] / (sqrt(2) - v[4]) } # projected vertices vs <- apply(vertices, c(1,2), stereog) # plot library(rgl) ## plot the edges for(k in 1:(2*A*B)){ p1 <- vs[, edges[1L,1L,k], edges[2L,1L,k]] p2 <- vs[, edges[1L,2L,k], edges[2L,2L,k]] shade3d(cylinder3d(rbind(p1,p2), radius = 0.03, sides = 30), color = "gold") } ## plot the vertices for(i in 1:A){ for(j in 1:B){ spheres3d(vs[,i,j], radius = 0.05, color = "gold2") } } ``` ![](figures/Duoprism_R.png) # Drawing a duoprism with Asymptote ```cpp settings.render = 4; settings.outformat = "pdf"; import solids; 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"); // lexicographic order --------------------------------------------------------- bool dominates(int[] e1, int[] e2){ return e2[0]>e1[0] || (e2[0]==e1[0] && e2[1]>e1[1]); } // vertices -------------------------------------------------------------------- int A = 30; int B = 4; 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(i/B*2pi); poly2[i][1] = sin(i/B*2pi); } real[][][] vertices = new real[A][B][4]; for(int i = 0; i < A; ++i){ for(int j = 0; j < B; ++j){ vertices[i][j] = concat(poly1[i],poly2[j]); } } // 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); } } } // stereographic projection ---------------------------------------------------- triple stereog(real[] x, real r){ return (x[0],x[1],x[2])/(r-x[3]); } // projected vertices ---------------------------------------------------------- triple[][] vs = new triple[A][B]; for(int i = 0; i < A; ++i){ for(int j = 0; j < B; ++j){ vs[i][j] = stereog(vertices[i][j], sqrt(2)); } } // draw the duoprism ----------------------------------------------------------- // draw edges for(int k = 0; k < edges.length; ++k){ path3 p = vs[edges[k][0][0]][edges[k][0][1]]..vs[edges[k][1][0]][edges[k][1][1]]; draw(tube(p, scale(0.05)*unitcircle), material(rgb(255,215,0), emissivepen=gray(0.1))); } // draw vertices for(int i = 0; i < A; ++i){ for(int j = 0; j < B; ++j){ draw(shift(vs[i][j])*scale3(0.1)*unitsphere, rgb(255,215,0)); } } ``` ![](figures/Duoprism_Asy.png) # Drawing a duoprism with POV-Ray ```povray #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 } // 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 White } 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 = 23; #declare B = 29; #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 // macros ---------------------------------------------------------------------- #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) * #end #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 // texture --------------------------------------------------------------------- #declare edgeTexture = texture { New_Penny finish { ambient 0.01 diffuse 2 reflection 0 brilliance 8 specular 0.1 roughness 0.1 } }; // draw an edge ---------------------------------------------------------------- #macro Edge(verts,pair1,pair2) cylinder { verts[pair1[0]][pair1[1]] verts[pair2[0]][pair2[1]], 0.05 texture { edgeTexture } } #end // draw ------------------------------------------------------------------------ #declare vs = ProjectedVertices(0, 0, 2*frame_number*pi/180); object { union { #for(i,0,nedges-1) Edge(vs,edges[i][0],edges[i][1]) #end #for(i,0,A-1) #for(j,0,B-1) sphere { vs[i][j], 0.06 texture { edgeTexture } } #end #end } scale 0.5 rotate <90, 0, 0> translate <0, 0.5, -2> } /* ini file Input_File_Name = Duoprism.pov Initial_Clock = 0 Final_Clock = 1 Initial_Frame = 0 Final_Frame = 179 Cyclic_Animation = on */ ``` ![](figures/Duoprism_POVRAY.gif)