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