1 /*****************************************************************************
2 * #ident "Id: main.c,v 3.27 2002-01-06 16:23:01+02 rl Exp "
5 * Kaleidoscopic construction of uniform polyhedra
6 * Copyright (c) 1991-2002 Dr. Zvi Har'El <rl@math.technion.ac.il>
8 * Redistribution and use in source and binary forms, with or without
9 * modification, are permitted provided that the following conditions
12 * 1. Redistributions of source code must retain the above copyright
13 * notice, this list of conditions and the following disclaimer.
15 * 2. Redistributions in binary form must reproduce the above copyright
16 * notice, this list of conditions and the following disclaimer in
17 * the documentation and/or other materials provided with the
20 * 3. The end-user documentation included with the redistribution,
21 * if any, must include the following acknowledgment:
22 * "This product includes software developed by
23 * Dr. Zvi Har'El (http://www.math.technion.ac.il/~rl/)."
24 * Alternately, this acknowledgment may appear in the software itself,
25 * if and wherever such third-party acknowledgments normally appear.
27 * This software is provided 'as-is', without any express or implied
28 * warranty. In no event will the author be held liable for any
29 * damages arising from the use of this software.
33 * Deptartment of Mathematics,
34 * Technion, Israel Institue of Technology,
35 * Haifa 32000, Israel.
36 * E-Mail: rl@math.technion.ac.il
38 * ftp://ftp.math.technion.ac.il/kaleido/
39 * http://www.mathconsult.ch/showroom/unipoly/
41 * Adapted for xscreensaver by Jamie Zawinski <jwz@jwz.org> 25-Apr-2004
43 *****************************************************************************
57 #include "polyhedra.h"
59 extern const char *progname;
62 #define MAXLONG 0x7FFFFFFF
65 #define MAXDIGITS 10 /* (int)log10((double)MAXLONG) + 1 */
69 #define DBL_EPSILON 2.2204460492503131e-16
71 #define BIG_EPSILON 3e-2
72 #define AZ M_PI/7 /* axis azimuth */
73 #define EL M_PI/17 /* axis elevation */
76 fprintf (stderr, "%s: %s\n", progname, (x)); \
80 #define Free(lvalue) do {\
82 free((char*) lvalue);\
87 #define Matfree(lvalue,n) do {\
89 matfree((char*) lvalue, n);\
93 #define Malloc(lvalue,n,type) do {\
94 if (!(lvalue = (type*) calloc((n), sizeof(type)))) \
98 #define Realloc(lvalue,n,type) do {\
99 if (!(lvalue = (type*) realloc(lvalue, (n) * sizeof(type)))) \
103 #define Calloc(lvalue,n,type) do {\
104 if (!(lvalue = (type*) calloc(n, sizeof(type))))\
108 #define Matalloc(lvalue,n,m,type) do {\
109 if (!(lvalue = (type**) matalloc(n, (m) * sizeof(type))))\
113 #define Sprintfrac(lvalue,x) do {\
114 if (!(lvalue=sprintfrac(x)))\
118 #define numerator(x) (frac(x), frax.n)
119 #define denominator(x) (frac(x), frax.d)
120 #define compl(x) (frac(x), (double) frax.n / (frax.n-frax.d))
127 /* NOTE: some of the int's can be replaced by short's, char's,
128 or even bit fields, at the expense of readability!!!*/
129 int index; /* index to the standard list, the array uniform[] */
130 int N; /* number of faces types (atmost 5)*/
131 int M; /* vertex valency (may be big for dihedral polyhedra) */
132 int V; /* vertex count */
133 int E; /* edge count */
134 int F; /* face count */
136 int chi; /* Euler characteristic */
137 int g; /* order of symmetry group */
138 int K; /* symmetry type: D=2, T=3, O=4, I=5 */
139 int hemi;/* flag hemi polyhedron */
140 int onesided;/* flag onesided polyhedron */
141 int even; /* removed face in pqr| */
142 int *Fi; /* face counts by type (array N)*/
143 int *rot; /* vertex configuration (array M of 0..N-1) */
144 int *snub; /* snub triangle configuration (array M of 0..1) */
145 int *firstrot; /* temporary for vertex generation (array V) */
146 int *anti; /* temporary for direction of ideal vertices (array E) */
147 int *ftype; /* face types (array F) */
148 int **e; /* edges (matrix 2 x E of 0..V-1)*/
149 int **dual_e; /* dual edges (matrix 2 x E of 0..F-1)*/
150 int **incid; /* vertex-face incidence (matrix M x V of 0..F-1)*/
151 int **adj; /* vertex-vertex adjacency (matrix M x V of 0..V-1)*/
152 double p[4]; /* p, q and r; |=0 */
153 double minr; /* smallest nonzero inradius */
154 double gon; /* basis type for dihedral polyhedra */
155 double *n; /* number of side of a face of each type (array N) */
156 double *m; /* number of faces at a vertex of each type (array N) */
157 double *gamma; /* fundamental angles in radians (array N) */
158 char *polyform; /* printable Wythoff symbol */
159 char *config; /* printable vertex configuration */
160 char *group; /* printable group name */
161 char *name; /* name, standard or manifuctured */
162 char *dual_name; /* dual name, standard or manifuctured */
165 Vector *v; /* vertex coordinates (array V) */
166 Vector *f; /* face coordinates (array F)*/
173 static Polyhedron *polyalloc(void);
174 static Vector rotate(Vector vertex, Vector axis, double angle);
176 static Vector sum3(Vector a, Vector b, Vector c);
177 static Vector scale(double k, Vector a);
178 static Vector sum(Vector a, Vector b);
179 static Vector diff(Vector a, Vector b);
180 static Vector pole (double r, Vector a, Vector b, Vector c);
181 static Vector cross(Vector a, Vector b);
182 static double dot(Vector a, Vector b);
183 static int same(Vector a, Vector b, double epsilon);
185 static char *sprintfrac(double x);
187 static void frac(double x);
188 static void matfree(void *mat, int rows);
189 static void *matalloc(int rows, int row_size);
191 static Fraction frax;
194 static const struct {
195 char *Wythoff, *name, *dual, *group, *class, *dual_class;
196 short Coxeter, Wenninger;
199 /****************************************************************************
200 * Dihedral Schwarz Triangles (D5 only)
201 ***************************************************************************/
203 /* 0 */ {"2 5|2", "Pentagonal Prism",
204 "Pentagonal Dipyramid",
210 /* 2 */ {"|2 2 5", "Pentagonal Antiprism",
211 "Pentagonal Deltohedron",
216 /* (2 2 5/2) (D2/5) */
217 /* 4 */ {"2 5/2|2", "Pentagrammic Prism",
218 "Pentagrammic Dipyramid",
224 /* 6 */ {"|2 2 5/2", "Pentagrammic Antiprism",
225 "Pentagrammic Deltohedron",
230 /* (5/3 2 2) (D3/5) */
232 /* 8 */ {"|2 2 5/3", "Pentagrammic Crossed Antiprism",
233 "Pentagrammic Concave Deltohedron",
239 /****************************************************************************
241 ***************************************************************************/
244 /* 10 */ {"3|2 3", "Tetrahedron",
246 "Tetrahedral (T[1])",
251 /* 12 */ {"2 3|3", "Truncated Tetrahedron",
252 "Triakistetrahedron",
253 "Tetrahedral (T[1])",
258 /* 14 */ {"3/2 3|3", "Octahemioctahedron",
260 "Tetrahedral (T[2])",
266 /* 16 */ {"3/2 3|2", "Tetrahemihexahedron",
268 "Tetrahedral (T[3])",
273 /****************************************************************************
275 ***************************************************************************/
278 /* 18 */ {"4|2 3", "Octahedron",
285 /* 20 */ {"3|2 4", "Cube",
292 /* 22 */ {"2|3 4", "Cuboctahedron",
293 "Rhombic Dodecahedron",
299 /* 24 */ {"2 4|3", "Truncated Octahedron",
300 "Tetrakishexahedron",
306 /* 26 */ {"2 3|4", "Truncated Cube",
313 /* 28 */ {"3 4|2", "Rhombicuboctahedron",
314 "Deltoidal Icositetrahedron",
320 /* 30 */ {"2 3 4|", "Truncated Cuboctahedron",
321 "Disdyakisdodecahedron",
327 /* 32, 33, 66, and 67 are chiral, existing in both left and right handed
328 (enantiomeric) forms, so it would make sense to display both versions.
331 /* 32 */ {"|2 3 4", "Snub Cube",
332 "Pentagonal Icositetrahedron",
333 "Octahedral (O[1]), Chiral",
337 /* (3/2 4 4) (O2b) */
339 /* 34 */ {"3/2 4|4", "Small Cubicuboctahedron",
340 "Small Hexacronic Icositetrahedron",
341 "Octahedral (O[2b])",
347 /* 36 */ {"3 4|4/3", "Great Cubicuboctahedron",
348 "Great Hexacronic Icositetrahedron",
354 /* 38 */ {"4/3 4|3", "Cubohemioctahedron",
361 /* 40 */ {"4/3 3 4|", "Cubitruncated Cuboctahedron",
362 "Tetradyakishexahedron",
369 /* 42 */ {"3/2 4|2", "Great Rhombicuboctahedron",
370 "Great Deltoidal Icositetrahedron",
376 /* 44 */ {"3/2 2 4|", "Small Rhombihexahedron",
377 "Small Rhombihexacron",
384 /* 46 */ {"2 3|4/3", "Stellated Truncated Hexahedron",
385 "Great Triakisoctahedron",
391 /* 48 */ {"4/3 2 3|", "Great Truncated Cuboctahedron",
392 "Great Disdyakisdodecahedron",
397 /* (4/3 3/2 2) (O11) */
399 /* 50 */ {"4/3 3/2 2|", "Great Rhombihexahedron",
400 "Great Rhombihexacron",
401 "Octahedral (O[11])",
406 /****************************************************************************
408 ***************************************************************************/
411 /* 52 */ {"5|2 3", "Icosahedron",
413 "Icosahedral (I[1])",
418 /* 54 */ {"3|2 5", "Dodecahedron",
420 "Icosahedral (I[1])",
425 /* 56 */ {"2|3 5", "Icosidodecahedron",
426 "Rhombic Triacontahedron",
427 "Icosahedral (I[1])",
432 /* 58 */ {"2 5|3", "Truncated Icosahedron",
433 "Pentakisdodecahedron",
434 "Icosahedral (I[1])",
439 /* 60 */ {"2 3|5", "Truncated Dodecahedron",
440 "Triakisicosahedron",
441 "Icosahedral (I[1])",
446 /* 62 */ {"3 5|2", "Rhombicosidodecahedron",
447 "Deltoidal Hexecontahedron",
448 "Icosahedral (I[1])",
453 /* 64 */ {"2 3 5|", "Truncated Icosidodecahedron",
454 "Disdyakistriacontahedron",
455 "Icosahedral (I[1])",
460 /* 32, 33, 66, and 67 are chiral, existing in both left and right handed
461 (enantiomeric) forms, so it would make sense to display both versions.
464 /* 66 */ {"|2 3 5", "Snub Dodecahedron",
465 "Pentagonal Hexecontahedron",
466 "Icosahedral (I[1]), Chiral",
470 /* (5/2 3 3) (I2a) */
472 /* 68 */ {"3|5/2 3", "Small Ditrigonal Icosidodecahedron",
473 "Small Triambic Icosahedron",
474 "Icosahedral (I[2a])",
479 /* 70 */ {"5/2 3|3", "Small Icosicosidodecahedron",
480 "Small Icosacronic Hexecontahedron",
481 "Icosahedral (I[2a])",
486 /* 72 */ {"|5/2 3 3", "Small Snub Icosicosidodecahedron",
487 "Small Hexagonal Hexecontahedron",
488 "Icosahedral (I[2a])",
492 /* (3/2 5 5) (I2b) */
494 /* 74 */ {"3/2 5|5", "Small Dodecicosidodecahedron",
495 "Small Dodecacronic Hexecontahedron",
496 "Icosahedral (I[2b])",
502 /* 76 */ {"5|2 5/2", "Small Stellated Dodecahedron",
503 "Great Dodecahedron",
504 "Icosahedral (I[3])",
505 "Kepler-Poinsot Solid",
506 "Kepler-Poinsot Solid",
509 /* 78 */ {"5/2|2 5", "Great Dodecahedron",
510 "Small Stellated Dodecahedron",
511 "Icosahedral (I[3])",
512 "Kepler-Poinsot Solid",
513 "Kepler-Poinsot Solid",
516 /* 80 */ {"2|5/2 5", "Great Dodecadodecahedron",
517 "Medial Rhombic Triacontahedron",
518 "Icosahedral (I[3])",
523 /* 82 */ {"2 5/2|5", "Truncated Great Dodecahedron",
524 "Small Stellapentakisdodecahedron",
525 "Icosahedral (I[3])",
530 /* 84 */ {"5/2 5|2", "Rhombidodecadodecahedron",
531 "Medial Deltoidal Hexecontahedron",
532 "Icosahedral (I[3])",
537 /* 86 */ {"2 5/2 5|", "Small Rhombidodecahedron",
538 "Small Rhombidodecacron",
539 "Icosahedral (I[3])",
544 /* 88 */ {"|2 5/2 5", "Snub Dodecadodecahedron",
545 "Medial Pentagonal Hexecontahedron",
546 "Icosahedral (I[3])",
552 /* 90 */ {"3|5/3 5", "Ditrigonal Dodecadodecahedron",
553 "Medial Triambic Icosahedron",
554 "Icosahedral (I[4])",
559 /* 92 */ {"3 5|5/3", "Great Ditrigonal Dodecicosidodecahedron",
560 "Great Ditrigonal Dodecacronic Hexecontahedron",
561 "Icosahedral (I[4])",
566 /* 94 */ {"5/3 3|5", "Small Ditrigonal Dodecicosidodecahedron",
567 "Small Ditrigonal Dodecacronic Hexecontahedron",
568 "Icosahedral (I[4])",
573 /* 96 */ {"5/3 5|3", "Icosidodecadodecahedron",
574 "Medial Icosacronic Hexecontahedron",
575 "Icosahedral (I[4])",
580 /* 98 */ {"5/3 3 5|", "Icositruncated Dodecadodecahedron",
581 "Tridyakisicosahedron",
582 "Icosahedral (I[4])",
587 /* 100 */ {"|5/3 3 5", "Snub Icosidodecadodecahedron",
588 "Medial Hexagonal Hexecontahedron",
589 "Icosahedral (I[4])",
593 /* (3/2 3 5) (I6b) */
595 /* 102 */ {"3/2|3 5", "Great Ditrigonal Icosidodecahedron",
596 "Great Triambic Icosahedron",
597 "Icosahedral (I[6b])",
602 /* 104 */ {"3/2 5|3", "Great Icosicosidodecahedron",
603 "Great Icosacronic Hexecontahedron",
604 "Icosahedral (I[6b])",
609 /* 106 */ {"3/2 3|5", "Small Icosihemidodecahedron",
610 "Small Icosihemidodecacron",
611 "Icosahedral (I[6b])",
616 /* 108 */ {"3/2 3 5|", "Small Dodecicosahedron",
617 "Small Dodecicosacron",
618 "Icosahedral (I[6b])",
622 /* (5/4 5 5) (I6c) */
624 /* 110 */ {"5/4 5|5", "Small Dodecahemidodecahedron",
625 "Small Dodecahemidodecacron",
626 "Icosahedral (I[6c])",
632 /* 112 */ {"3|2 5/2", "Great Stellated Dodecahedron",
634 "Icosahedral (I[7])",
635 "Kepler-Poinsot Solid",
636 "Kepler-Poinsot Solid",
639 /* 114 */ {"5/2|2 3", "Great Icosahedron",
640 "Great Stellated Dodecahedron",
641 "Icosahedral (I[7])",
642 "Kepler-Poinsot Solid",
643 "Kepler-Poinsot Solid",
646 /* 116 */ {"2|5/2 3", "Great Icosidodecahedron",
647 "Great Rhombic Triacontahedron",
648 "Icosahedral (I[7])",
653 /* 118 */ {"2 5/2|3", "Great Truncated Icosahedron",
654 "Great Stellapentakisdodecahedron",
655 "Icosahedral (I[7])",
660 /* 120 */ {"2 5/2 3|", "Rhombicosahedron",
662 "Icosahedral (I[7])",
667 /* 122 */ {"|2 5/2 3", "Great Snub Icosidodecahedron",
668 "Great Pentagonal Hexecontahedron",
669 "Icosahedral (I[7])",
675 /* 124 */ {"2 5|5/3", "Small Stellated Truncated Dodecahedron",
676 "Great Pentakisdodecahedron",
677 "Icosahedral (I[9])",
682 /* 126 */ {"5/3 2 5|", "Truncated Dodecadodecahedron",
683 "Medial Disdyakistriacontahedron",
684 "Icosahedral (I[9])",
689 /* 128 */ {"|5/3 2 5", "Inverted Snub Dodecadodecahedron",
690 "Medial Inverted Pentagonal Hexecontahedron",
691 "Icosahedral (I[9])",
695 /* (5/3 5/2 3) (I10a) */
697 /* 130 */ {"5/2 3|5/3", "Great Dodecicosidodecahedron",
698 "Great Dodecacronic Hexecontahedron",
699 "Icosahedral (I[10a])",
704 /* 132 */ {"5/3 5/2|3", "Small Dodecahemicosahedron",
705 "Small Dodecahemicosacron",
706 "Icosahedral (I[10a])",
711 /* 134 */ {"5/3 5/2 3|", "Great Dodecicosahedron",
712 "Great Dodecicosacron",
713 "Icosahedral (I[10a])",
718 /* 136 */ {"|5/3 5/2 3", "Great Snub Dodecicosidodecahedron",
719 "Great Hexagonal Hexecontahedron",
720 "Icosahedral (I[10a])",
724 /* (5/4 3 5) (I10b) */
726 /* 138 */ {"5/4 5|3", "Great Dodecahemicosahedron",
727 "Great Dodecahemicosacron",
728 "Icosahedral (I[10b])",
732 /* (5/3 2 3) (I13) */
734 /* 140 */ {"2 3|5/3", "Great Stellated Truncated Dodecahedron",
735 "Great Triakisicosahedron",
736 "Icosahedral (I[13])",
741 /* 142 */ {"5/3 3|2", "Great Rhombicosidodecahedron",
742 "Great Deltoidal Hexecontahedron",
743 "Icosahedral (I[13])",
748 /* 144 */ {"5/3 2 3|", "Great Truncated Icosidodecahedron",
749 "Great Disdyakistriacontahedron",
750 "Icosahedral (I[13])",
755 /* 146 */ {"|5/3 2 3", "Great Inverted Snub Icosidodecahedron",
756 "Great Inverted Pentagonal Hexecontahedron",
757 "Icosahedral (I[13])",
761 /* (5/3 5/3 5/2) (I18a) */
763 /* 148 */ {"5/3 5/2|5/3", "Great Dodecahemidodecahedron",
764 "Great Dodecahemidodecacron",
765 "Icosahedral (I[18a])",
769 /* (3/2 5/3 3) (I18b) */
771 /* 150 */ {"3/2 3|5/3", "Great Icosihemidodecahedron",
772 "Great Icosihemidodecacron",
773 "Icosahedral (I[18b])",
777 /* (3/2 3/2 5/3) (I22) */
779 /* 152 */ {"|3/2 3/2 5/2","Small Retrosnub Icosicosidodecahedron",
780 "Small Hexagrammic Hexecontahedron",
781 "Icosahedral (I[22])",
785 /* (3/2 5/3 2) (I23) */
787 /* 154 */ {"3/2 5/3 2|", "Great Rhombidodecahedron",
788 "Great Rhombidodecacron",
789 "Icosahedral (I[23])",
794 /* 156 */ {"|3/2 5/3 2", "Great Retrosnub Icosidodecahedron",
795 "Great Pentagrammic Hexecontahedron",
796 "Icosahedral (I[23])",
801 /****************************************************************************
803 ***************************************************************************/
805 /* 158 */ {"3/2 5/3 3 5/2", "Great Dirhombicosidodecahedron",
806 "Great Dirhombicosidodecacron",
813 static int last_uniform = sizeof (uniform) / sizeof (uniform[0]);
817 static int unpacksym(const char *sym, Polyhedron *P);
818 static int moebius(Polyhedron *P);
819 static int decompose(Polyhedron *P);
820 static int guessname(Polyhedron *P);
821 static int newton(Polyhedron *P, int need_approx);
822 static int exceptions(Polyhedron *P);
823 static int count(Polyhedron *P);
824 static int configuration(Polyhedron *P);
825 static int vertices(Polyhedron *P);
826 static int faces(Polyhedron *P);
827 static int edgelist(Polyhedron *P);
830 kaleido(const char *sym,
831 int need_coordinates, int need_edgelist, int need_approx,
836 * Allocate a Polyhedron structure P.
838 if (!(P = polyalloc()))
841 * Unpack input symbol into P.
843 if (!unpacksym(sym, P))
846 * Find Mebius triangle, its density and Euler characteristic.
851 * Decompose Schwarz triangle.
856 * Find the names of the polyhedron and its dual.
863 * Solve Fundamental triangles, optionally printing approximations.
865 if (!newton(P,need_approx))
868 * Deal with exceptional polyhedra.
873 * Count edges and faces, update density and characteristic if needed.
878 * Generate printable vertex configuration.
880 if (!configuration(P))
883 * Compute coordinates.
885 if (!need_coordinates && !need_edgelist)
902 * Allocate a blank Polyhedron structure and initialize some of its nonblank
905 * Array and matrix field are allocated when needed.
911 Calloc(P, 1, Polyhedron);
919 * Free the struture allocated by polyalloc(), as well as all the array and
923 polyfree(Polyhedron *P)
943 Matfree(P->dual_e, 2);
944 Matfree(P->incid, P->M);
945 Matfree(P->adj, P->M);
950 matalloc(int rows, int row_size)
954 if (!(mat = malloc(rows * sizeof (void *))))
956 while ((mat[i] = malloc(row_size)) && ++i < rows)
967 matfree(void *mat, int rows)
970 free(((void **)mat)[rows]);
975 * compute the mathematical modulus function.
980 return (i%=j)>=0?i:j<0?i-j:i+j;
985 * Find the numerator and the denominator using the Euclidean algorithm.
990 static const Fraction zero = {0,1}, inf = {1,0};
997 if (fabs(s) > (double) MAXLONG)
999 f = (long) floor (s);
1002 frax.n = frax.n * f + r0.n;
1003 frax.d = frax.d * f + r0.d;
1004 if (x == (double)frax.n/(double)frax.d)
1012 * Unpack input symbol: Wythoff symbol or an index to uniform[]. The symbol is
1013 * a # followed by a number, or a three fractions and a bar in some order. We
1014 * allow no bars only if it result from the input symbol #80.
1017 unpacksym(const char *sym, Polyhedron *P)
1019 int i = 0, n, d, bars = 0;
1021 while ((c = *sym++) && isspace(c))
1023 if (!c) Err("no data");
1025 while ((c = *sym++) && isspace(c))
1028 Err("no digit after #");
1032 while ((c = *sym++) && isdigit(c))
1033 n = n * 10 + c - '0';
1036 if (n > last_uniform)
1037 Err("index too big");
1039 while ((c = *sym++) && isspace(c))
1042 Err("data exceeded");
1043 sym = uniform[P->index = n - 1].Wythoff;
1048 while ((c = *sym++) && isspace(c))
1051 if (i == 4 && (bars || P->index == last_uniform - 1))
1055 Err("not enough fractions");
1058 Err("data exceeded");
1061 Err("too many bars");
1068 while ((c = *sym++) && isdigit(c))
1069 n = n * 10 + c - '0';
1070 if (c && isspace (c))
1071 while ((c = *sym++) && isspace(c))
1075 if ((P->p[i++] = n) <= 1)
1079 while ((c = *sym++) && isspace(c))
1081 if (!c || !isdigit(c))
1084 while ((c = *sym++) && isdigit(c))
1085 d = d * 10 + c - '0';
1087 Err("zero denominator");
1089 if ((P->p[i++] = (double) n / d) <= 1)
1095 * Using Wythoff symbol (p|qr, pq|r, pqr| or |pqr), find the Moebius triangle
1096 * (2 3 K) (or (2 2 n)) of the Schwarz triangle (pqr), the order g of its
1097 * symmetry group, its Euler characteristic chi, and its covering density D.
1098 * g is the number of copies of (2 3 K) covering the sphere, i.e.,
1100 * g * pi * (1/2 + 1/3 + 1/K - 1) = 4 * pi
1102 * D is the number of times g copies of (pqr) cover the sphere, i.e.
1104 * D * 4 * pi = g * pi * (1/p + 1/q + 1/r - 1)
1106 * chi is V - E + F, where F = g is the number of triangles, E = 3*g/2 is the
1107 * number of triangle edges, and V = Vp+ Vq+ Vr, with Vp = g/(2*np) being the
1108 * number of vertices with angle pi/p (np is the numerator of p).
1111 moebius(Polyhedron *P)
1113 int twos = 0, j, len = 1;
1115 * Arrange Wythoff symbol in a presentable form. In the same time check the
1116 * restrictions on the three fractions: They all have to be greater then one,
1117 * and the numerators 4 or 5 cannot occur together. We count the ocurrences
1118 * of 2 in `two', and save the largest numerator in `P->K', since they
1119 * reflect on the symmetry group.
1122 if (P->index == last_uniform - 1) {
1123 Malloc(P->polyform, ++len, char);
1124 strcpy(P->polyform, "|");
1126 Calloc(P->polyform, len, char);
1127 for (j = 0; j < 4; j++) {
1130 Sprintfrac(s, P->p[j]);
1131 if (j && P->p[j-1]) {
1132 Realloc(P->polyform, len += strlen (s) + 1, char);
1133 strcat(P->polyform, " ");
1135 Realloc (P->polyform, len += strlen (s), char);
1136 strcat(P->polyform, s);
1140 if ((k = numerator (P->p[j])) > P->K) {
1144 } else if (k < P->K && k == 4)
1149 Realloc(P->polyform, ++len, char);
1150 strcat(P->polyform, "|");
1154 * Find the symmetry group P->K (where 2, 3, 4, 5 represent the dihedral,
1155 * tetrahedral, octahedral and icosahedral groups, respectively), and its
1158 if (twos >= 2) {/* dihedral */
1163 Err("numerator too large");
1164 P->g = 24 * P->K / (6 - P->K);
1167 * Compute the nominal density P->D and Euler characteristic P->chi.
1168 * In few exceptional cases, these values will be modified later.
1170 if (P->index != last_uniform - 1) {
1172 P->D = P->chi = - P->g;
1173 for (j = 0; j < 4; j++) if (P->p[j]) {
1174 P->chi += i = P->g / numerator(P->p[j]);
1175 P->D += i * denominator(P->p[j]);
1180 Err("nonpositive density");
1186 * Decompose Schwarz triangle into N right triangles and compute the vertex
1187 * count V and the vertex valency M. V is computed from the number g of
1188 * Schwarz triangles in the cover, divided by the number of triangles which
1189 * share a vertex. It is halved for one-sided polyhedra, because the
1190 * kaleidoscopic construction really produces a double orientable covering of
1191 * such polyhedra. All q' q|r are of the "hemi" type, i.e. have equatorial {2r}
1192 * faces, and therefore are (except 3/2 3|3 and the dihedra 2 2|r) one-sided. A
1193 * well known example is 3/2 3|4, the "one-sided heptahedron". Also, all p q r|
1194 * with one even denominator have a crossed parallelogram as a vertex figure,
1195 * and thus are one-sided as well.
1198 decompose(Polyhedron *P)
1201 if (!P->p[1]) { /* p|q r */
1203 P->M = 2 * numerator(P->p[0]);
1205 Malloc(P->n, P->N, double);
1206 Malloc(P->m, P->N, double);
1207 Malloc(P->rot, P->M, int);
1209 for (j = 0; j < 2; j++) {
1210 P->n[j] = P->p[j+2];
1213 for (j = P->M / 2; j--;) {
1217 } else if (!P->p[2]) { /* p q|r */
1221 Malloc(P->n, P->N, double);
1222 Malloc(P->m, P->N, double);
1223 Malloc(P->rot, P->M, int);
1225 P->n[0] = 2 * P->p[3];
1227 for (j = 1; j < 3; j++) {
1228 P->n[j] = P->p[j-1];
1233 if (fabs(P->p[0] - compl (P->p[1])) < DBL_EPSILON) {/* p = q' */
1234 /* P->p[0]==compl(P->p[1]) should work. However, MSDOS
1235 * yeilds a 7e-17 difference! Reported by Jim Buddenhagen
1236 * <jb1556@daditz.sbc.com> */
1239 if (P->p[0] != 2 && !(P->p[3] == 3 && (P->p[0] == 3 ||
1246 } else if (!P->p[3]) { /* p q r| */
1249 Malloc(P->n, P->N, double);
1250 Malloc(P->m, P->N, double);
1251 Malloc(P->rot, P->M, int);
1253 for (j = 0; j < 3; j++) {
1254 if (!(denominator(P->p[j]) % 2)) {
1255 /* what happens if there is more then one even denominator? */
1256 if (P->p[(j+1)%3] != P->p[(j+2)%3]) { /* needs postprocessing */
1257 P->even = j;/* memorize the removed face */
1258 P->chi -= P->g / numerator(P->p[j]) / 2;
1261 } else {/* for p = q we get a double 2 2r|p */
1262 /* noted by Roman Maeder <maeder@inf.ethz.ch> for 4 4 3/2| */
1263 /* Euler characteristic is still wrong */
1268 P->n[j] = 2 * P->p[j];
1272 } else { /* |p q r - snub polyhedron */
1275 P->V = P->g / 2;/* Only "white" triangles carry a vertex */
1276 Malloc(P->n, P->N, double);
1277 Malloc(P->m, P->N, double);
1278 Malloc(P->rot, P->M, int);
1279 Malloc(P->snub, P->M, int);
1282 P->m[0] = P->n[0] = 3;
1283 for (j = 1; j < 4; j++) {
1293 * Sort the fundamental triangles (using bubble sort) according to decreasing
1294 * n[i], while pushing the trivial triangles (n[i] = 2) to the end.
1301 for (j = 0; j < last; j++) {
1302 if ((P->n[j] < P->n[j+1] || P->n[j] == 2) && P->n[j+1] != 2) {
1306 P->n[j] = P->n[j+1];
1309 P->m[j] = P->m[j+1];
1311 for (i = 0; i < P->M; i++) {
1314 else if (P->rot[i] == j+1)
1317 if (P->even != -1) {
1320 else if (P->even == j+1)
1328 * Get rid of repeated triangles.
1330 for (J = 0; J < P->N && P->n[J] != 2;J++) {
1332 for (j = J+1; j < P->N && P->n[j]==P->n[J]; j++)
1336 for (i = j; i < P->N; i++) {
1337 P->n[i - k] = P->n[i];
1338 P->m[i - k] = P->m[i];
1341 for (i = 0; i < P->M; i++) {
1344 else if (P->rot[i] > J)
1352 * Get rid of trivial triangles.
1355 J = 1; /* hosohedron */
1359 for (i = 0; i < P->M; i++) {
1360 if (P->rot[i] >= P->N) {
1361 for (j = i + 1; j < P->M; j++) {
1362 P->rot[j-1] = P->rot[j];
1364 P->snub[j-1] = P->snub[j];
1373 Realloc(P->n, P->N, double);
1374 Realloc(P->m, P->N, double);
1375 Realloc(P->rot, P->M, int);
1377 Realloc(P->snub, P->M, int);
1382 static int dihedral(Polyhedron *P, const char *name, const char *dual_name);
1386 * Get the polyhedron name, using standard list or guesswork. Ideally, we
1387 * should try to locate the Wythoff symbol in the standard list (unless, of
1388 * course, it is dihedral), after doing few normalizations, such as sorting
1389 * angles and splitting isoceles triangles.
1392 guessname(Polyhedron *P)
1394 if (P->index != -1) {/* tabulated */
1395 P->name = uniform[P->index].name;
1396 P->dual_name = uniform[P->index].dual;
1397 P->group = uniform[P->index].group;
1398 P->class = uniform[P->index].class;
1399 P->dual_class = uniform[P->index].dual_class;
1401 } else if (P->K == 2) {/* dihedral nontabulated */
1404 Malloc(P->name, sizeof ("Octahedron"), char);
1405 Malloc(P->dual_name, sizeof ("Cube"), char);
1406 strcpy(P->name, "Octahedron");
1407 strcpy(P->dual_name, "Cube");
1410 P->gon = P->n[0] == 3 ? P->n[1] : P->n[0];
1412 return dihedral(P, "Antiprism", "Deltohedron");
1414 return dihedral(P, "Crossed Antiprism", "Concave Deltohedron");
1415 } else if (!P->p[3] ||
1419 Malloc(P->name, sizeof("Cube"), char);
1420 Malloc(P->dual_name, sizeof("Octahedron"), char);
1421 strcpy(P->name, "Cube");
1422 strcpy(P->dual_name, "Octahedron");
1425 P->gon = P->n[0] == 4 ? P->n[1] : P->n[0];
1426 return dihedral(P, "Prism", "Dipyramid");
1427 } else if (!P->p[1] && P->p[0] != 2) {
1429 return dihedral(P, "Hosohedron", "Dihedron");
1432 return dihedral(P, "Dihedron", "Hosohedron");
1434 } else {/* other nontabulated */
1435 static const char *pre[] = {"Tetr", "Oct", "Icos"};
1436 Malloc(P->name, 50, char);
1437 Malloc(P->dual_name, 50, char);
1438 sprintf(P->name, "%sahedral ", pre[P->K - 3]);
1440 strcat (P->name, "One-Sided ");
1442 strcat(P->name, "Convex ");
1444 strcat(P->name, "Nonconvex ");
1445 strcpy(P->dual_name, P->name);
1446 strcat(P->name, "Isogonal Polyhedron");
1447 strcat(P->dual_name, "Isohedral Polyhedron");
1448 Realloc(P->name, strlen (P->name) + 1, char);
1449 Realloc(P->dual_name, strlen (P->dual_name) + 1, char);
1455 dihedral(Polyhedron *P, const char *name, const char *dual_name)
1459 Sprintfrac(s, P->gon < 2 ? compl (P->gon) : P->gon);
1460 i = strlen(s) + sizeof ("-gonal ");
1461 Malloc(P->name, i + strlen (name), char);
1462 Malloc(P->dual_name, i + strlen (dual_name), char);
1463 sprintf(P->name, "%s-gonal %s", s, name);
1464 sprintf(P->dual_name, "%s-gonal %s", s, dual_name);
1470 * Solve the fundamental right spherical triangles.
1471 * If need_approx is set, print iterations on standard error.
1474 newton(Polyhedron *P, int need_approx)
1477 * First, we find initial approximations.
1481 Malloc(P->gamma, P->N, double);
1483 P->gamma[0] = M_PI / P->m[0];
1486 for (j = 0; j < P->N; j++)
1487 P->gamma[j] = M_PI / 2 - M_PI / P->n[j];
1488 errno = 0; /* may be non-zero from some reason */
1490 * Next, iteratively find closer approximations for gamma[0] and compute
1491 * other gamma[j]'s from Napier's equations.
1494 fprintf(stderr, "Solving %s\n", P->polyform);
1496 double delta = M_PI, sigma = 0;
1497 for (j = 0; j < P->N; j++) {
1499 fprintf(stderr, "%-20.15f", P->gamma[j]);
1500 delta -= P->m[j] * P->gamma[j];
1503 printf("(%g)\n", delta);
1504 if (fabs(delta) < 11 * DBL_EPSILON)
1506 /* On a RS/6000, fabs(delta)/DBL_EPSILON may occilate between 8 and
1507 * 10. Reported by David W. Sanderson <dws@ssec.wisc.edu> */
1508 for (j = 0; j < P->N; j++)
1509 sigma += P->m[j] * tan(P->gamma[j]);
1510 P->gamma[0] += delta * tan(P->gamma[0]) / sigma;
1511 if (P->gamma[0] < 0 || P->gamma[0] > M_PI)
1512 Err("gamma out of bounds");
1513 cosa = cos(M_PI / P->n[0]) / sin(P->gamma[0]);
1514 for (j = 1; j < P->N; j++)
1515 P->gamma[j] = asin(cos(M_PI / P->n[j]) / cosa);
1517 Err(strerror(errno));
1522 * Postprocess pqr| where r has an even denominator (cf. Coxeter &al. Sec.9).
1523 * Remove the {2r} and add a retrograde {2p} and retrograde {2q}.
1526 exceptions(Polyhedron *P)
1529 if (P->even != -1) {
1531 Realloc(P->n, P->N, double);
1532 Realloc(P->m, P->N, double);
1533 Realloc(P->gamma, P->N, double);
1534 Realloc(P->rot, P->M, int);
1535 for (j = P->even + 1; j < 3; j++) {
1536 P->n[j-1] = P->n[j];
1537 P->gamma[j-1] = P->gamma[j];
1539 P->n[2] = compl(P->n[1]);
1540 P->gamma[2] = - P->gamma[1];
1541 P->n[3] = compl(P->n[0]);
1543 P->gamma[3] = - P->gamma[0];
1551 * Postprocess the last polyhedron |3/2 5/3 3 5/2 by taking a |5/3 3 5/2,
1552 * replacing the three snub triangles by four equatorial squares and adding
1553 * the missing {3/2} (retrograde triangle, cf. Coxeter &al. Sec. 11).
1555 if (P->index == last_uniform - 1) {
1558 Realloc(P->n, P->N, double);
1559 Realloc(P->m, P->N, double);
1560 Realloc(P->gamma, P->N, double);
1561 Realloc(P->rot, P->M, int);
1562 Realloc(P->snub, P->M, int);
1565 for (j = 3; j; j--) {
1567 P->n[j] = P->n[j-1];
1568 P->gamma[j] = P->gamma[j-1];
1570 P->m[0] = P->n[0] = 4;
1571 P->gamma[0] = M_PI / 2;
1573 P->n[4] = compl(P->n[1]);
1574 P->gamma[4] = - P->gamma[1];
1575 for (j = 1; j < 6; j += 2) P->rot[j]++;
1585 * Compute edge and face counts, and update D and chi. Update D in the few
1586 * cases the density of the polyhedron is meaningful but different than the
1587 * density of the corresponding Schwarz triangle (cf. Coxeter &al., p. 418 and
1589 * In these cases, spherical faces of one type are concave (bigger than a
1590 * hemisphere), and the actual density is the number of these faces less the
1591 * computed density. Note that if j != 0, the assignment gamma[j] = asin(...)
1592 * implies gamma[j] cannot be obtuse. Also, compute chi for the only
1593 * non-Wythoffian polyhedron.
1596 count(Polyhedron *P)
1599 Malloc(P->Fi, P->N, int);
1600 for (j = 0; j < P->N; j++) {
1601 P->E += temp = P->V * numerator(P->m[j]);
1602 P->F += P->Fi[j] = temp / numerator(P->n[j]);
1605 if (P->D && P->gamma[0] > M_PI / 2)
1606 P->D = P->Fi[0] - P->D;
1607 if (P->index == last_uniform - 1)
1608 P->chi = P->V - P->E + P->F;
1613 * Generate a printable vertex configuration symbol.
1616 configuration(Polyhedron *P)
1619 for (j = 0; j < P->M; j++) {
1621 Sprintfrac(s, P->n[P->rot[j]]);
1622 len += strlen (s) + 2;
1624 Malloc(P->config, len, char);
1625 /* strcpy(P->config, "(");*/
1626 strcpy(P->config, "");
1628 Realloc(P->config, len, char);
1629 strcat(P->config, ", ");
1631 strcat(P->config, s);
1634 /* strcat (P->config, ")");*/
1635 if ((j = denominator (P->m[0])) != 1) {
1636 char s[MAXDIGITS + 2];
1637 sprintf(s, "/%d", j);
1638 Realloc(P->config, len + strlen (s), char);
1639 strcat(P->config, s);
1645 * Compute polyhedron vertices and vertex adjecency lists.
1646 * The vertices adjacent to v[i] are v[adj[0][i], v[adj[1][i], ...
1647 * v[adj[M-1][i], ordered counterclockwise. The algorith is a BFS on the
1648 * vertices, in such a way that the vetices adjacent to a givem vertex are
1649 * obtained from its BFS parent by a cyclic sequence of rotations. firstrot[i]
1650 * points to the first rotaion in the sequence when applied to v[i]. Note that
1651 * for non-snub polyhedra, the rotations at a child are opposite in sense when
1652 * compared to the rotations at the parent. Thus, we fill adj[*][i] from the
1653 * end to signify clockwise rotations. The firstrot[] array is not needed for
1654 * display thus it is freed after being used for face computations below.
1657 vertices(Polyhedron *P)
1661 Malloc(P->v, P->V, Vector);
1662 Matalloc(P->adj, P->M, P->V, int);
1663 Malloc(P->firstrot, P->V, int); /* temporary , put in Polyhedron
1664 structure so that may be freed on
1666 cosa = cos(M_PI / P->n[0]) / sin(P->gamma[0]);
1672 P->v[1].x = 2 * cosa * sqrt(1 - cosa * cosa);
1674 P->v[1].z = 2 * cosa * cosa - 1;
1677 P->adj[0][1] = -1;/* start the other side */
1678 P->adj[P->M-1][1] = 0;
1680 P->firstrot[1] = P->snub[P->M-1] ? 0 : P->M-1 ;
1683 for (i = 0; i < newV; i++) {
1685 int last, one, start, limit;
1686 if (P->adj[0][i] == -1) {
1687 one = -1; start = P->M-2; limit = -1;
1689 one = 1; start = 1; limit = P->M;
1692 for (j = start; j != limit; j += one) {
1695 temp = rotate (P->v[P->adj[j-one][i]], P->v[i],
1696 one * 2 * P->gamma[P->rot[k]]);
1697 for (J=0; J<newV && !same(P->v[J],temp,BIG_EPSILON); J++)
1703 if (J == newV) { /* new vertex */
1704 if (newV == P->V) Err ("too many vertices");
1705 P->v[newV++] = temp;
1710 P->adj[P->M-1][J] = i;
1715 P->firstrot[J] = !P->snub[last] ? last :
1716 !P->snub[k] ? (k+1)%P->M : k ;
1726 * Compute polyhedron faces (dual vertices) and incidence matrices.
1727 * For orientable polyhedra, we can distinguish between the two faces meeting
1728 * at a given directed edge and identify the face on the left and the face on
1729 * the right, as seen from the outside. For one-sided polyhedra, the vertex
1730 * figure is a papillon (in Coxeter &al. terminology, a crossed parallelogram)
1731 * and the two faces meeting at an edge can be identified as the side face
1732 * (n[1] or n[2]) and the diagonal face (n[0] or n[3]).
1735 faces(Polyhedron *P)
1738 Malloc (P->f, P->F, Vector);
1739 Malloc (P->ftype, P->F, int);
1740 Matalloc (P->incid, P->M, P->V, int);
1741 P->minr = 1 / fabs (tan (M_PI / P->n[P->hemi]) * tan (P->gamma[P->hemi]));
1742 for (i = P->M; --i>=0;) {
1744 for (j = P->V; --j>=0;)
1745 P->incid[i][j] = -1;
1747 for (i = 0; i < P->V; i++) {
1749 for (j = 0; j < P->M; j++) {
1751 int pap=0;/* papillon edge type */
1752 if (P->incid[j][i] != -1)
1754 P->incid[j][i] = newF;
1756 Err("too many faces");
1757 P->f[newF] = pole(P->minr, P->v[i], P->v[P->adj[j][i]],
1758 P->v[P->adj[mod(j + 1, P->M)][i]]);
1759 P->ftype[newF] = P->rot[mod(P->firstrot[i] + ((P->adj[0][i] <
1760 P->adj[P->M - 1][i])
1765 pap = (P->firstrot[i] + j) % 2;
1771 if ((i0 = P->adj[J][k]) == i) break;
1772 for (J = 0; J < P->M && P->adj[J][i0] != k; J++)
1775 Err("too many faces");
1776 if (P->onesided && (J + P->firstrot[i0]) % 2 == pap) {
1777 P->incid [J][i0] = newF;
1783 P->incid [J][i0] = newF;
1796 * Compute edge list and graph polyhedron and dual.
1797 * If the polyhedron is of the "hemi" type, each edge has one finite vertex and
1798 * one ideal vertex. We make sure the latter is always the out-vertex, so that
1799 * the edge becomes a ray (half-line). Each ideal vertex is represented by a
1800 * unit Vector, and the direction of the ray is either parallel or
1801 * anti-parallel this Vector. We flag this in the array P->anti[E].
1804 edgelist(Polyhedron *P)
1806 int i, j, *s, *t, *u;
1807 Matalloc(P->e, 2, P->E, int);
1808 Matalloc(P->dual_e, 2, P->E, int);
1811 for (i = 0; i < P->V; i++)
1812 for (j = 0; j < P->M; j++)
1813 if (i < P->adj[j][i]) {
1815 *t++ = P->adj[j][i];
1822 Malloc(P->anti, P->E, int);
1824 for (i = 0; i < P->V; i++)
1825 for (j = 0; j < P->M; j++)
1826 if (i < P->adj[j][i])
1829 *s++ = P->incid[mod(j-1,P->M)][i];
1830 *t++ = P->incid[j][i];
1832 if (P->ftype[P->incid[j][i]]) {
1833 *s = P->incid[j][i];
1834 *t = P->incid[mod(j-1,P->M)][i];
1836 *s = P->incid[mod(j-1,P->M)][i];
1837 *t = P->incid[j][i];
1839 *u++ = dot(P->f[*s++], P->f[*t++]) > 0;
1847 sprintfrac(double x)
1852 Malloc(s, sizeof ("infinity"), char);
1853 strcpy(s, "infinity");
1854 } else if (frax.d == 1) {
1855 char n[MAXDIGITS + 1];
1856 sprintf(n, "%ld", frax.n);
1857 Malloc(s, strlen (n) + 1, char);
1860 char n[MAXDIGITS + 1], d[MAXDIGITS + 1];
1861 sprintf(n, "%ld", frax.n);
1862 sprintf(d, "%ld", frax.d);
1863 Malloc(s, strlen (n) + strlen (d) + 2, char);
1864 sprintf(s, "%s/%s", n, d);
1870 dot(Vector a, Vector b)
1872 return a.x * b.x + a.y * b.y + a.z * b.z;
1876 scale(double k, Vector a)
1885 diff(Vector a, Vector b)
1894 cross(Vector a, Vector b)
1897 p.x = a.y * b.z - a.z * b.y;
1898 p.y = a.z * b.x - a.x * b.z;
1899 p.z = a.x * b.y - a.y * b.x;
1904 sum(Vector a, Vector b)
1913 sum3(Vector a, Vector b, Vector c)
1922 rotate(Vector vertex, Vector axis, double angle)
1925 p = scale(dot (axis, vertex), axis);
1926 return sum3(p, scale(cos(angle), diff(vertex, p)),
1927 scale(sin(angle), cross(axis, vertex)));
1930 static Vector x, y, z;
1933 * rotate the standard frame
1936 rotframe(double azimuth, double elevation, double angle)
1938 static const Vector X = {1,0,0}, Y = {0,1,0}, Z = {0,0,1};
1941 axis = rotate(rotate (X, Y, elevation), Z, azimuth);
1942 x = rotate(X, axis, angle);
1943 y = rotate(Y, axis, angle);
1944 z = rotate(Z, axis, angle);
1948 * rotate an array of n Vectors
1951 rotarray(Vector *new, Vector *old, int n)
1954 *new++ = sum3(scale(old->x, x), scale(old->y, y), scale(old->z, z));
1960 same(Vector a, Vector b, double epsilon)
1962 return fabs(a.x - b.x) < epsilon && fabs(a.y - b.y) < epsilon
1963 && fabs(a.z - b.z) < epsilon;
1967 * Compute the polar reciprocal of the plane containing a, b and c:
1969 * If this plane does not contain the origin, return p such that
1970 * dot(p,a) = dot(p,b) = dot(p,b) = r.
1972 * Otherwise, return p such that
1973 * dot(p,a) = dot(p,b) = dot(p,c) = 0
1978 pole(double r, Vector a, Vector b, Vector c)
1982 p = cross(diff(b, a), diff(c, a));
1985 return scale(1 / sqrt(dot(p, p)), p);
1987 return scale(r/ k , p);
1996 static void rotframe(double azimuth, double elevation, double angle);
1997 static void rotarray(Vector *new, Vector *old, int n);
1998 static int mod (int i, int j);
2002 push_point (polyhedron *p, Vector v)
2004 p->points[p->npoints].x = v.x;
2005 p->points[p->npoints].y = v.y;
2006 p->points[p->npoints].z = v.z;
2011 push_face3 (polyhedron *p, int x, int y, int z)
2013 p->faces[p->nfaces].npoints = 3;
2014 Malloc (p->faces[p->nfaces].points, 3, int);
2015 p->faces[p->nfaces].points[0] = x;
2016 p->faces[p->nfaces].points[1] = y;
2017 p->faces[p->nfaces].points[2] = z;
2022 push_face4 (polyhedron *p, int x, int y, int z, int w)
2024 p->faces[p->nfaces].npoints = 4;
2025 Malloc (p->faces[p->nfaces].points, 4, int);
2026 p->faces[p->nfaces].points[0] = x;
2027 p->faces[p->nfaces].points[1] = y;
2028 p->faces[p->nfaces].points[2] = z;
2029 p->faces[p->nfaces].points[3] = w;
2037 construct_polyhedron (Polyhedron *P, Vector *v, int V, Vector *f, int F,
2038 const char *name, const char *dual,
2039 const char *class, const char *star,
2040 double azimuth, double elevation, double freeze)
2042 int i, j, k=0, l, ll, ii, *hit=0, facelets;
2047 Malloc (result, 1, polyhedron);
2048 memset (result, 0, sizeof(*result));
2053 rotframe(azimuth, elevation, freeze);
2054 Malloc(temp, V, Vector);
2055 rotarray(temp, v, V);
2057 Malloc(temp, F, Vector);
2058 rotarray(temp, f, F);
2061 result->number = P->index + 1;
2062 result->name = strdup (name);
2063 result->dual = strdup (dual);
2064 result->wythoff = strdup (P->polyform);
2065 result->config = strdup (P->config);
2066 result->group = strdup (P->group);
2067 result->class = strdup (class);
2072 Malloc (result->points, V + F * 13, point);
2073 result->npoints = 0;
2075 result->nedges = P->E;
2076 result->logical_faces = F;
2077 result->logical_vertices = V;
2078 result->density = P->D;
2079 result->chi = P->chi;
2081 for (i = 0; i < V; i++)
2082 push_point (result, v[i]);
2085 * Auxiliary vertices (needed because current VRML browsers cannot handle
2086 * non-simple polygons, i.e., ploygons with self intersections): Each
2087 * non-simple face is assigned an auxiliary vertex. By connecting it to the
2088 * rest of the vertices the face is triangulated. The circum-center is used
2089 * for the regular star faces of uniform polyhedra. The in-center is used for
2090 * the pentagram (#79) and hexagram (#77) of the high-density snub duals, and
2091 * for the pentagrams (#40, #58) and hexagram (#52) of the stellated duals
2092 * with configuration (....)/2. Finally, the self-intersection of the crossed
2093 * parallelogram is used for duals with form p q r| with an even denominator.
2095 * This method do not work for the hemi-duals, whose faces are not
2096 * star-shaped and have two self-intersections each.
2098 * Thus, for each face we need six auxiliary vertices: The self intersections
2099 * and the terminal points of the truncations of the infinite edges. The
2100 * ideal vertices are listed, but are not used by the face-list.
2102 * Note that the face of the last dual (#80) is octagonal, and constists of
2103 * two quadrilaterals of the infinite type.
2106 if (*star && P->even != -1)
2107 Malloc(hit, F, int);
2108 for (i = 0; i < F; i++)
2110 (frac(P->n[P->ftype[i]]), frax.d != 1 && frax.d != frax.n - 1)) ||
2114 denominator (P->m[0]) != 1))) {
2115 /* find the center of the face */
2117 if (!*star && P->hemi && !P->ftype[i])
2120 h = P->minr / dot(f[i],f[i]);
2121 push_point(result, scale (h, f[i]));
2123 } else if (*star && P->even != -1) {
2124 /* find the self-intersection of a crossed parallelogram.
2125 * hit is set if v0v1 intersects v2v3*/
2126 Vector v0, v1, v2, v3, c0, c1, p;
2128 v0 = v[P->incid[0][i]];
2129 v1 = v[P->incid[1][i]];
2130 v2 = v[P->incid[2][i]];
2131 v3 = v[P->incid[3][i]];
2132 d0 = sqrt(dot(diff(v0, v2), diff(v0, v2)));
2133 d1 = sqrt(dot (diff(v1, v3), diff(v1, v3)));
2134 c0 = scale(d1, sum(v0, v2));
2135 c1 = scale(d0, sum(v1, v3));
2136 p = scale(0.5 / (d0 + d1), sum(c0, c1));
2137 push_point (result, p);
2138 p = cross(diff(p, v2), diff(p, v3));
2139 hit[i] = (dot(p, p) < 1e-6);
2140 } else if (*star && P->hemi && P->index != last_uniform - 1) {
2141 /* find the terminal points of the truncation and the
2142 * self-intersections.
2149 Vector v0, v1, v2, v3, v01, v03, v21, v23, v0123, v0321 ;
2151 double t = 1.5;/* truncation adjustment factor */
2152 j = !P->ftype[P->incid[0][i]];
2153 v0 = v[P->incid[j][i]];/* real vertex */
2154 v1 = v[P->incid[j+1][i]];/* ideal vertex (unit vector) */
2155 v2 = v[P->incid[j+2][i]];/* real */
2156 v3 = v[P->incid[(j+3)%4][i]];/* ideal */
2157 /* compute intersections
2158 * this uses the following linear algebra:
2159 * v0123 = v0 + a v1 = v2 + b v3
2160 * v0 x v3 + a (v1 x v3) = v2 x v3
2161 * a (v1 x v3) = (v2 - v0) x v3
2162 * a (v1 x v3) . (v1 x v3) = (v2 - v0) x v3 . (v1 x v3)
2165 v0123 = sum(v0, scale(dot(cross(diff(v2, v0), v3), u) / dot(u,u),
2167 v0321 = sum(v0, scale(dot(cross(diff(v0, v2), v1), u) / dot(u,u),
2169 /* compute truncations */
2170 v01 = sum(v0 , scale(t, diff(v0123, v0)));
2171 v23 = sum(v2 , scale(t, diff(v0123, v2)));
2172 v03 = sum(v0 , scale(t, diff(v0321, v0)));
2173 v21 = sum(v2 , scale(t, diff(v0321, v2)));
2175 push_point(result, v01);
2176 push_point(result, v23);
2177 push_point(result, v0123);
2178 push_point(result, v03);
2179 push_point(result, v21);
2180 push_point(result, v0321);
2182 } else if (*star && P->index == last_uniform - 1) {
2183 /* find the terminal points of the truncation and the
2184 * self-intersections.
2197 Vector v0, v1, v2, v3, v4, v5, v6, v7, v01, v07, v21, v23;
2198 Vector v43, v45, v65, v67, v0123, v0721, v4365, v4567;
2199 double t = 1.5;/* truncation adjustment factor */
2201 for (j = 0; j < 8; j++)
2202 if (P->ftype[P->incid[j][i]] == 3)
2204 v0 = v[P->incid[j][i]];/* real {5/3} */
2205 v1 = v[P->incid[(j+1)%8][i]];/* ideal */
2206 v2 = v[P->incid[(j+2)%8][i]];/* real {3} */
2207 v3 = v[P->incid[(j+3)%8][i]];/* ideal */
2208 v4 = v[P->incid[(j+4)%8][i]];/* real {5/2} */
2209 v5 = v[P->incid[(j+5)%8][i]];/* ideal */
2210 v6 = v[P->incid[(j+6)%8][i]];/* real {3/2} */
2211 v7 = v[P->incid[(j+7)%8][i]];/* ideal */
2212 /* compute intersections */
2214 v0123 = sum(v0, scale(dot(cross(diff(v2, v0), v3), u) / dot(u,u),
2217 v0721 = sum(v0, scale(dot(cross(diff(v2, v0), v1), u) / dot(u,u),
2220 v4567 = sum(v4, scale(dot(cross(diff(v6, v4), v7), u) / dot(u,u),
2223 v4365 = sum(v4, scale(dot(cross(diff(v6, v4), v5), u) / dot(u,u),
2225 /* compute truncations */
2226 v01 = sum(v0 , scale(t, diff(v0123, v0)));
2227 v23 = sum(v2 , scale(t, diff(v0123, v2)));
2228 v07 = sum(v0 , scale(t, diff(v0721, v0)));
2229 v21 = sum(v2 , scale(t, diff(v0721, v2)));
2230 v45 = sum(v4 , scale(t, diff(v4567, v4)));
2231 v67 = sum(v6 , scale(t, diff(v4567, v6)));
2232 v43 = sum(v4 , scale(t, diff(v4365, v4)));
2233 v65 = sum(v6 , scale(t, diff(v4365, v6)));
2235 push_point(result, v01);
2236 push_point(result, v23);
2237 push_point(result, v0123);
2238 push_point(result, v07);
2239 push_point(result, v21);
2240 push_point(result, v0721);
2241 push_point(result, v45);
2242 push_point(result, v67);
2243 push_point(result, v4567);
2244 push_point(result, v43);
2245 push_point(result, v65);
2246 push_point(result, v4365);
2251 * Each face is printed in a separate line, by listing the indices of its
2252 * vertices. In the non-simple case, the polygon is represented by the
2253 * triangulation, each triangle consists of two polyhedron vertices and one
2256 Malloc (result->faces, F * 10, face);
2261 for (i = 0; i < F; i++) {
2265 denominator (P->m[0]) != 1)) {
2266 for (j = 0; j < P->M - 1; j++) {
2267 push_face3 (result, P->incid[j][i], P->incid[j+1][i], ii);
2271 push_face3 (result, P->incid[j][i], P->incid[0][i], ii++);
2274 } else if (P->even != -1) {
2275 if (hit && hit[i]) {
2276 push_face3 (result, P->incid[3][i], P->incid[0][i], ii);
2277 push_face3 (result, P->incid[1][i], P->incid[2][i], ii);
2279 push_face3 (result, P->incid[0][i], P->incid[1][i], ii);
2280 push_face3 (result, P->incid[2][i], P->incid[3][i], ii);
2285 } else if (P->hemi && P->index != last_uniform - 1) {
2286 j = !P->ftype[P->incid[0][i]];
2288 push_face3 (result, ii, ii + 1, ii + 2);
2289 push_face4 (result, P->incid[j][i], ii + 2, P->incid[j+2][i], ii + 5);
2290 push_face3 (result, ii + 3, ii + 4, ii + 5);
2293 } else if (P->index == last_uniform - 1) {
2294 for (j = 0; j < 8; j++)
2295 if (P->ftype[P->incid[j][i]] == 3)
2297 push_face3 (result, ii, ii + 1, ii + 2);
2299 P->incid[j][i], ii + 2, P->incid[(j+2)%8][i], ii + 5);
2300 push_face3 (result, ii + 3, ii + 4, ii + 5);
2302 push_face3 (result, ii + 6, ii + 7, ii + 8);
2304 P->incid[(j+4)%8][i], ii + 8, P->incid[(j+6)%8][i],
2306 push_face3 (result, ii + 9, ii + 10, ii + 11);
2311 result->faces[result->nfaces].npoints = P->M;
2312 Malloc (result->faces[result->nfaces].points, P->M, int);
2313 for (j = 0; j < P->M; j++)
2314 result->faces[result->nfaces].points[j] = P->incid[j][i];
2319 int split = (frac(P->n[P->ftype[i]]),
2320 frax.d != 1 && frax.d != frax.n - 1);
2321 for (j = 0; j < V; j++) {
2322 for (k = 0; k < P->M; k++)
2323 if (P->incid[k][j] == i)
2330 for (l = P->adj[k][j]; l != j; l = P->adj[k][l]) {
2331 for (k = 0; k < P->M; k++)
2332 if (P->incid[k][l] == i)
2334 if (P->adj[k][l] == ll)
2335 k = mod(k + 1 , P->M);
2336 push_face3 (result, ll, l, ii);
2340 push_face3 (result, ll, j, ii++);
2347 Malloc (pp, 100, int);
2351 for (l = P->adj[k][j]; l != j; l = P->adj[k][l]) {
2352 for (k = 0; k < P->M; k++)
2353 if (P->incid[k][l] == i)
2355 if (P->adj[k][l] == ll)
2356 k = mod(k + 1 , P->M);
2360 result->faces[result->nfaces].npoints = pi;
2361 result->faces[result->nfaces].points = pp;
2369 * Face color indices - for polyhedra with multiple face types
2370 * For non-simple faces, the index is repeated as many times as needed by the
2375 if (!*star && P->N != 1) {
2376 for (i = 0; i < F; i++)
2377 if (frac(P->n[P->ftype[i]]), frax.d == 1 || frax.d == frax.n - 1)
2378 result->faces[ff++].color = P->ftype[i];
2380 for (j = 0; j < frax.n; j++)
2381 result->faces[ff++].color = P->ftype[i];
2383 for (i = 0; i < facelets; i++)
2384 result->faces[ff++].color = 0;
2388 if (*star && P->even != -1)
2398 /* External interface (jwz)
2402 free_polyhedron (polyhedron *p)
2412 for (i = 0; i < p->nfaces; i++)
2413 Free (p->faces[i].points);
2421 construct_polyhedra (polyhedron ***polyhedra_ret)
2424 double azimuth = AZ;
2425 double elevation = EL;
2429 polyhedron **result;
2430 Malloc (result, last_uniform * 2 + 3, polyhedron*);
2432 while (index < last_uniform) {
2436 sprintf(sym, "#%d", index + 1);
2437 if (!(P = kaleido(sym, 1, 0, 0, 0))) {
2438 Err (strerror(errno));
2441 result[count++] = construct_polyhedron (P, P->v, P->V, P->f, P->F,
2442 P->name, P->dual_name,
2444 azimuth, elevation, freeze);
2446 result[count++] = construct_polyhedron (P, P->f, P->F, P->v, P->V,
2447 P->dual_name, P->name,
2449 azimuth, elevation, freeze);
2454 *polyhedra_ret = result;
2455 count++; /* leave room for teapot */