1 /*
2 * $Id: plato.i,v 1.1 2005-09-18 22:06:16 dhmunro Exp $
3 */
4 /* Copyright (c) 2005, The Regents of the University of California.
5 * All rights reserved.
6 * This file is part of yorick (http://yorick.sourceforge.net).
7 * Read the accompanying LICENSE file for details.
8 */
9
10 local plato;
11 /* DOCUMENT plato.i
12 Contains routines to generate points related to Platonic
13 solids and other pleasing or simple 3D geometrical objects.
14
15 pt_tet All these return points of the solid imbedded in
16 pt_cube a pleasing way inside the cube with corners +-1.
17 pt_oct With a non-zero parameter, the points are instead
18 pt_dodec normalized to unit length.
19 pt_ico
20
21 bucky return points or faces of geodesic dome-like
22 solids and their solid angles
23 */
24
pt_tet(norm)25 func pt_tet(norm)
26 {
27 /* half of corners of a cube, no particular order */
28 p= [[1,1,1],[-1,-1,1],[1,-1,-1],[-1,1,-1]];
29 /* use other choice of corners if norm<0 */
30 if (norm) p*= sqrt(1./3.)*sign(norm);
31 return p;
32 }
33
pt_cube(norm)34 func pt_cube(norm)
35 {
36 /* points at corners of cube, hence 2x2x2 organization
37 * where indices are x, y, z directions */
38 p= array(-1,3,2,2,2);
39 p(1,2,,)= p(2,,2,)= p(3,,,2)= 1;
40 if (norm) p*= sqrt(1./3.);
41 return p;
42 }
43
pt_oct(norm)44 func pt_oct(norm)
45 {
46 /* one point on each face of a cube, hence 2x3 organization
47 * where 1st index selects + or - face of cube
48 * 2nd index selects xyz face of cube */
49 p= array(0,3,2,3);
50 p(1,,1)= p(2,,2)= p(3,,3)= [-1,1];
51 return p;
52 }
53
pt_ico(norm)54 func pt_ico(norm)
55 {
56 /* two points on each face of a cube, hence 2x2x3 organization,
57 * where 1st index selects + or - point on face of cube,
58 * 2nd index selects + or - face of cube
59 * 3rd index selects xyz face of cube */
60 g= 0.5*(sqrt(5.)-1.); /* reciprocal golden ratio */
61 p= [[[0,g,1],[0,-g,1]],[[0,g,-1],[0,-g,-1]]];
62 p= [p, roll(p,[1,0,0]), roll(p,[2,0,0])];
63 if (norm) p/= abs(p(1,..),p(2,..),p(3,..))(-,..);
64 return p;
65 }
66
pt_dodec(norm)67 func pt_dodec(norm)
68 {
69 /* two points on each face of a cube, followed by corners of cube */
70 g= 0.5*(sqrt(5.)-1.); /* reciprocal golden ratio */
71 g2= 1.-g; /* equals g*g */
72 p= [[[0,g2,1],[0,-g2,1]],[[0,g2,-1],[0,-g2,-1]]];
73 p= [p, roll(p,[1,0,0]), roll(p,[2,0,0])];
74 p= grow(p(,*),g*pt_cube()(,*));
75 if (norm) p/= abs(p(1,..),p(2,..),p(3,..))(-,..);
76 return p;
77 }
78
79 func bucky(n,faces,&domega)
80 {
81 /* two rings of five plus two points at poles
82 * rings have cos(theta) = +- g (reciprocal golden ratio),
83 * sin(theta) = sqrt(g)
84 *
85 * The points are organized into five spiral strips running
86 * southeast from the north pole to the south pole. Each strip
87 * consists of a 3x2 array of points which bound a strip of four
88 * triangular faces; the diagonals run from (1,1) to (2,2) and
89 * from (2,1) to (3,2). The point (1,2) is always the north pole,
90 * and the point (3,1) is always the south pole. The points
91 * (2,2) and (3,2) on the eastern edge of one strip are the same
92 * as the points (1,1) and (3,2), respectively, on the western
93 * edge of the strip immediately to the east.
94 * Hence, the dimensionality of the returned array of points
95 * is 3 by 3x2x5. (30 - 4 duplications of north pole
96 * - 4 duplications of south pole - 5*2 other duplicates = 12)
97 *
98 * When you specify n (default 0), bucky will halve each of the
99 * initial 3x2 strips n times, to produce a (2^(n+1)+1)x(2^n+1)
100 * array of equally spaced points in place of the 3x2 points of
101 * the n=0 pattern.
102 * n=0 --> 12 points 20 faces
103 * n=1 --> 42 points 80 faces
104 * n=2 --> 162 points 320 faces
105 * n=4 --> 642 points 1280 faces
106 * n=5 --> 2562 points 5120 faces
107 *
108 * When radius==1, icosahedron apothem==sqrt((2+g)/(2-g)/3),
109 * suggesting a worst case area ratio of (2+g)/(2-g)/3 = 0.631476.
110 * However, by renormalizing the points after each doubling, this
111 * ratio is considerably improved; the worst case in the limit of
112 * many doublings is a little under 0.769.
113 */
114 g= 0.5*(sqrt(5.)-1.); /* reciprocal golden ratio */
115 p= array(0., 3,3,2,5);
116 p(,1,2,)= [0,0,1]; /* north pole */
117 p(,3,1,)= [0,0,-1]; /* south pole */
118 c36= 0.5*(1.+g); /* equals 0.5/g */
119 s36= 0.5*sqrt(2.-g);
120 c72= 0.5*g;
121 s72= 0.5*sqrt(3.+g);
122 g= 1./(1.+2.*g); /* cos of ring angle */
123 s= 2.*g; /* sin is twice cos for this magic angle */
124 ringn= s*[[1,0],[c72,s72],[-c36,s36],[-c36,-s36],[c72,-s72]];
125 rings= s*[[c36,s36],[-c72,s72],[-1,0],[-c72,-s72],[c36,-s36]];
126 p(3,1,1,)= p(3,2,2,)= g;
127 p(1:2,1,1,)= ringn;
128 p(1:2,2,2,)= roll(ringn,[0,-1]);
129 p(3,2,1,)= p(3,3,2,)= -g;
130 p(1:2,2,1,)= rings;
131 p(1:2,3,2,)= roll(rings,[0,-1]);
132 if (!n) n= 0;
133 while (n--) {
134 dims= dimsof(p);
135 dims(3:4)= 2*dims(3:4)-1;
136 q= array(0., dims);
137 q(,1:0:2,1:0:2,)= p;
138 q(,2:-1:2,1:0:2,)= p(,zcen,,);
139 q(,1:0:2,2:-1:2,)= p(,,zcen,);
140 q(,2:-1:2,2:-1:2,)= p(,2:0,2:0,)+p(,1:-1,1:-1,); /* *0.5 */
141 p= q/abs(q(1,..),q(2,..),q(3,..))(-,..);
142 q= x= [];
143 }
144 if (faces) {
145 dims= dimsof(p);
146 q= array(0.,3,2,dims(3)-1,dims(4)-1,5);
147 llur= (pb=p(,1:-1,1:-1,)) + (pc=p(,2:0,2:0,));
148 q(,1,,,)= llur + (pa=p(,1:-1,2:0,)); /* a-b-c */
149 q(,2,,,)= llur + (pd=p(,2:0,1:-1,)); /* d-c-b */
150 domega= q(1,..);
151 p= q/abs(domega,q(2,..),q(3,..))(-,..);
152 domega(1,..)= pt_solid2(pa,pb,pc);
153 domega(2,..)= pt_solid2(pd,pc,pb);
154 }
155 return p;
156 }
157
pt_cross(a,b)158 func pt_cross(a,b)
159 {
160 i = [2,3,1];
161 j = [3,1,2];
162 return a(i,..)*b(j,..) - a(j,..)*b(i,..);
163 }
164
pt_solid(a,b,c)165 func pt_solid(a,b,c)
166 {
167 vab= pt_cross(a,b);
168 sab= sqrt((vab*vab)(sum,..));
169 vbc= pt_cross(b,c);
170 sbc= sqrt((vbc*vbc)(sum,..));
171 vca= pt_cross(c,a);
172 sca= sqrt((vca*vca)(sum,..));
173 cosa= -(vab*vca)(sum,..)/(sab*sca);
174 cosb= -(vbc*vab)(sum,..)/(sbc*sab);
175 cosc= -(vca*vbc)(sum,..)/(sca*sbc);
176 /* this formula is simple, but inaccurate for small triangles */
177 return acos(cosa)+acos(cosb)+acos(cosc)-pi;
178 }
179
pt_solid2(a,b,c)180 func pt_solid2(a,b,c)
181 {
182 vab= pt_cross(a,b);
183 sab= sqrt((vab*vab)(sum,..));
184 vbc= pt_cross(b,c);
185 sbc= sqrt((vbc*vbc)(sum,..));
186 vca= pt_cross(c,a);
187 sca= sqrt((vca*vca)(sum,..));
188 da= 1./(sab*sca);
189 db= 1./(sbc*sab);
190 dc= 1./(sca*sbc);
191 cosa= -(vab*vca)(sum,..)*da;
192 cosb= -(vbc*vab)(sum,..)*db;
193 cosc= -(vca*vbc)(sum,..)*dc;
194 sina= pt_cross(vab,vca);
195 sina= sqrt((sina*sina)(sum,..))*da;
196 sinb= pt_cross(vab,vca);
197 sinb= sqrt((sinb*sinb)(sum,..))*db;
198 sinc= pt_cross(vab,vca);
199 sinc= sqrt((sinc*sinc)(sum,..))*dc;
200 sabc= abs(sina*(cosb*cosc-sinb*sinc)+cosa*(sinb*cosc+cosb*sinc));
201 cabc= 1. - cosa*(cosb*cosc-sinb*sinc) + sina*(sinb*cosc+cosb*sinc);
202 /* this formula is good for small triangles,
203 * bad for triangles with area near 2*pi */
204 return 4.*asin(0.5*sabc/sqrt(cabc*(1.+sqrt(0.5*cabc))));
205 }
206