Home Manual Packages Global Index Keywords Quick Reference ``` /* * plato.i -- \$Id\$ * * Copyright (c) 1998. See accompanying LEGAL file for details. */ local plato ; /* DOCUMENT plato.i Contains routines to generate points related to Platonic solids and other pleasing or simple 3D geometrical objects. pt_tet All these return points of the solid imbedded in pt_cube a pleasing way inside the cube with corners +-1. pt_oct With a non-zero parameter, the points are instead pt_dodec normalized to unit length. pt_ico bucky return points or faces of geodesic dome-like solids and their solid angles */ func pt_tet (norm) { /* half of corners of a cube, no particular order */ p= [[1,1,1],[-1,-1,1],[1,-1,-1],[-1,1,-1]]; /* use other choice of corners if norm<0 */ if (norm) p*= sqrt(1./3.)*sign(norm); return p; } func pt_cube (norm) { /* points at corners of cube, hence 2x2x2 organization * where indices are x, y, z directions */ p= array(-1,3,2,2,2); p(1,2,,)= p(2,,2,)= p(3,,,2)= 1; if (norm) p*= sqrt(1./3.); return p; } func pt_oct (norm) { /* one point on each face of a cube, hence 2x3 organization * where 1st index selects + or - face of cube * 2nd index selects xyz face of cube */ p= array(0,3,2,3); p(1,,1)= p(2,,2)= p(3,,3)= [-1,1]; return p; } func pt_ico (norm) { /* two points on each face of a cube, hence 2x2x3 organization, * where 1st index selects + or - point on face of cube, * 2nd index selects + or - face of cube * 3rd index selects xyz face of cube */ g= 0.5*(sqrt(5.)-1.); /* reciprocal golden ratio */ p= [[[0,g,1],[0,-g,1]],[[0,g,-1],[0,-g,-1]]]; p= [p, roll(p,[1,0,0]), roll(p,[2,0,0])]; if (norm) p/= abs(p(1,..),p(2,..),p(3,..))(-,..); return p; } func pt_dodec (norm) { /* two points on each face of a cube, followed by corners of cube */ g= 0.5*(sqrt(5.)-1.); /* reciprocal golden ratio */ g2= 1.-g; /* equals g*g */ p= [[[0,g2,1],[0,-g2,1]],[[0,g2,-1],[0,-g2,-1]]]; p= [p, roll(p,[1,0,0]), roll(p,[2,0,0])]; p= grow(p(,*),g*pt_cube()(,*)); if (norm) p/= abs(p(1,..),p(2,..),p(3,..))(-,..); return p; } func bucky (n,faces,&domega) { /* two rings of five plus two points at poles * rings have cos(theta) = +- g (reciprocal golden ratio), * sin(theta) = sqrt(g) * * The points are organized into five spiral strips running * southeast from the north pole to the south pole. Each strip * consists of a 3x2 array of points which bound a strip of four * triangular faces; the diagonals run from (1,1) to (2,2) and * from (2,1) to (3,2). The point (1,2) is always the north pole, * and the point (3,1) is always the south pole. The points * (2,2) and (3,2) on the eastern edge of one strip are the same * as the points (1,1) and (3,2), respectively, on the western * edge of the strip immediately to the east. * Hence, the dimensionality of the returned array of points * is 3 by 3x2x5. (30 - 4 duplications of north pole * - 4 duplications of south pole - 5*2 other duplicates = 12) * * When you specify n (default 0), bucky will halve each of the * initial 3x2 strips n times, to produce a (2^(n+1)+1)x(2^n+1) * array of equally spaced points in place of the 3x2 points of * the n=0 pattern. * n=0 --> 12 points 20 faces * n=1 --> 42 points 80 faces * n=2 --> 162 points 320 faces * n=4 --> 642 points 1280 faces * n=5 --> 2562 points 5120 faces * * When radius==1, icosahedron apothem==sqrt((2+g)/(2-g)/3), * suggesting a worst case area ratio of (2+g)/(2-g)/3 = 0.631476. * However, by renormalizing the points after each doubling, this * ratio is considerably improved; the worst case in the limit of * many doublings is a little under 0.769. */ g= 0.5*(sqrt(5.)-1.); /* reciprocal golden ratio */ p= array(0., 3,3,2,5); p(,1,2,)= [0,0,1]; /* north pole */ p(,3,1,)= [0,0,-1]; /* south pole */ c36= 0.5*(1.+g); /* equals 0.5/g */ s36= 0.5*sqrt(2.-g); c72= 0.5*g; s72= 0.5*sqrt(3.+g); g= 1./(1.+2.*g); /* cos of ring angle */ s= 2.*g; /* sin is twice cos for this magic angle */ ringn= s*[[1,0],[c72,s72],[-c36,s36],[-c36,-s36],[c72,-s72]]; rings= s*[[c36,s36],[-c72,s72],[-1,0],[-c72,-s72],[c36,-s36]]; p(3,1,1,)= p(3,2,2,)= g; p(1:2,1,1,)= ringn; p(1:2,2,2,)= roll(ringn,[0,-1]); p(3,2,1,)= p(3,3,2,)= -g; p(1:2,2,1,)= rings; p(1:2,3,2,)= roll(rings,[0,-1]); if (!n) n= 0; while (n--) { dims= dimsof(p); dims(3:4)= 2*dims(3:4)-1; q= array(0., dims); q(,1:0:2,1:0:2,)= p; q(,2:-1:2,1:0:2,)= p(,zcen,,); q(,1:0:2,2:-1:2,)= p(,,zcen,); q(,2:-1:2,2:-1:2,)= p(,2:0,2:0,)+p(,1:-1,1:-1,); /* *0.5 */ p= q/abs(q(1,..),q(2,..),q(3,..))(-,..); q= x= []; } if (faces) { dims= dimsof(p); q= array(0.,3,2,dims(3)-1,dims(4)-1,5); llur= (pb=p(,1:-1,1:-1,)) + (pc=p(,2:0,2:0,)); q(,1,,,)= llur + (pa=p(,1:-1,2:0,)); /* a-b-c */ q(,2,,,)= llur + (pd=p(,2:0,1:-1,)); /* d-c-b */ domega= q(1,..); p= q/abs(domega,q(2,..),q(3,..))(-,..); domega(1,..)= pt_solid2(pa,pb,pc); domega(2,..)= pt_solid2(pd,pc,pb); } return p; } func pt_cross (a,b) { c= 0.*a; c(1,..)= a(2,..)*b(3,..)-a(3,..)*b(2,..); c(2,..)= a(3,..)*b(1,..)-a(1,..)*b(3,..); c(3,..)= a(1,..)*b(2,..)-a(2,..)*b(1,..); return c; } func pt_solid (a,b,c) { vab= pt_cross(a,b); sab= sqrt((vab*vab)(sum,..)); vbc= pt_cross(b,c); sbc= sqrt((vbc*vbc)(sum,..)); vca= pt_cross(c,a); sca= sqrt((vca*vca)(sum,..)); cosa= -(vab*vca)(sum,..)/(sab*sca); cosb= -(vbc*vab)(sum,..)/(sbc*sab); cosc= -(vca*vbc)(sum,..)/(sca*sbc); /* this formula is simple, but inaccurate for small triangles */ return acos(cosa)+acos(cosb)+acos(cosc)-pi; } func pt_solid2 (a,b,c) { vab= pt_cross(a,b); sab= sqrt((vab*vab)(sum,..)); vbc= pt_cross(b,c); sbc= sqrt((vbc*vbc)(sum,..)); vca= pt_cross(c,a); sca= sqrt((vca*vca)(sum,..)); da= 1./(sab*sca); db= 1./(sbc*sab); dc= 1./(sca*sbc); cosa= -(vab*vca)(sum,..)*da; cosb= -(vbc*vab)(sum,..)*db; cosc= -(vca*vbc)(sum,..)*dc; sina= pt_cross(vab,vca); sina= sqrt((sina*sina)(sum,..))*da; sinb= pt_cross(vab,vca); sinb= sqrt((sinb*sinb)(sum,..))*db; sinc= pt_cross(vab,vca); sinc= sqrt((sinc*sinc)(sum,..))*dc; sabc= abs(sina*(cosb*cosc-sinb*sinc)+cosa*(sinb*cosc+cosb*sinc)); cabc= 1. - cosa*(cosb*cosc-sinb*sinc) + sina*(sinb*cosc+cosb*sinc); /* this formula is good for small triangles, * bad for triangles with area near 2*pi */ return 4.*asin(0.5*sabc/sqrt(cabc*(1.+sqrt(0.5*cabc)))); } ```