Rhombic Spirallohedra

plus
 Space-fillings, surfaces, and conchoids;
 plane sections, ruled surfaces,
 the volume of the general k-armed spirallohedron;
 the composition of spirallohedra by rhombic hexahedra,
 and more ...

Russell Towle; version of December, 2003
Box 141, Dutch Flat CA 95714    rtowle@inreach.com

Introduction; examples of spirallohedra

Rhombic Spirallohedra

In the July 2003 update, the special packing of mixed 3- and 6-armed spirallohedra is constructed, and plane sections of individual spirallohedra, or of packed spirallohedra, are found. Also, the “vertical lines” connecting the transformed equatorial faces of the “parent” polar zonohedron are found.

In the October 2002 update of the Rhombic Spirallohedra notebook, the volume of the general k-armed spirallohedron is found, and issues relating to the close-packings of 3- and 4-armed spirallohedra are explored. Spirallohedra arise by subtraction of some of the rhombic hexahedra which close-pack to compose a polar zonohedron; in the last, most recent section, we find these subtracted hexahedra, examine the volumes of hexahedra, and make spirallohedra from hexahedra. The (newest) last section of the notebook contains many subsections.

The discussion below describes what a spirallohedron is, and how these forms were discovered. Since rhombic spirallohedra are closely related to polar zonohedra, we begin there.


Polar zonohedra are bounded by n*(n-1) rhombic faces, where n is the number of directions in which edges occur. There is an axis of n-fold symmetry (herein, always coincident with the z axis) which intersects the zonohedron at each “pole,” where n equal rhombs meet. Simple instances include the cube (n=3) and Kepler’s rhombic dodecahedron (n=4); both of these exist in nature as crystals. For any one n>=3, an infinite variety of polar zonohedra may be constructed, ranging from oblate, squashed forms, to prolate, spindle-like forms; the determining parameter is the angle between edges and a plane orthogonal to the polar zonohedron’s symmetry axis, which I call “pitch,” 0<pitch<90 degrees. Curiously, when pitch equals arc tan Sqrt[1/2] or 35.2643+ degrees, a polar zonohedron is an isometric, orthogonal “shadow” of an n-cube (See §13 of H.S.M. Coxeter's Regular Polytopes, Dover, 1973).

(*A polar projection, and three equatorial projections,
of n=12 polar zonohedra.*)
Show[PolarZonohedron[12,35.2643],
PlotLabel->"n=12, pitch=35.2643",
ViewPoint->{0,0,3}];
Show[PolarZonohedron[12,35.2643],
PlotLabel->"n=12, pitch=35.2643",
ViewPoint->{0,-3,0}];
Show[PolarZonohedron[12,15],(*oblate*)
PlotLabel->"n=12, pitch=15",
ViewPoint->{0,-3,0}];
Show[PolarZonohedron[12,60],(*prolate*)
PlotLabel->"n=12, pitch=60",
ViewPoint->{0,-3,0}];

It so happens that the polar zonohedra can be dissected into the “tetrahedral” number, n-things-taken-three-at-a-time, rhombic hexahedra. These are somewhat like squished cubes, in fact, when pitch=35.2643+ degrees, they arise as projections of the ordinary cubes which bound, in part, an n-cube. When n is even moderately large, say, n=12, this dissection occurs in many different ways. When it is performed in every possible way the polar zonohedron of pitch=35.2643+ degrees is not just a shadow, but a full-fledged orthogonal projection, of the n-cube.

The general zonohedron is determined by a set of three or more vectors spanning the 3-space, and each vector determines a zone of parallel edges, and faces bounded in part by these edges, which “girdle” the zonohedron. Convex zonohedra have central symmetry, and every face, edge and vertex has an equal, opposite, and similarly-situated counter-face, counter-edge, and counter-vertex. Any two zones of faces cross one another at a face and its counter-face. Below, one of the twelve zones of faces of an n=12 polar zonohedron is depicted.

n=12;pitch=35.2643;v=vectors[n,pitch];
a=PolarZonohedron[n,pitch];

(*color the second zone differently*)
Show[colorzone[2,a,v],
LightSources -> 
   {{{-1, 1, 1}, RGBColor[1, 1, 1]},
   {{1, -2, 2}, RGBColor[1, 1, 1]}},
Boxed->False,
Background->GrayLevel[0],
ViewPoint->{0,-.5,2}];

(*show the second zone by itself*)
Show[onezone[2,a,v],
Boxed->False,
Background->GrayLevel[0],
ViewPoint->{0,-.5,2}];

(*show all faces not in the second zone*)
Show[notzone[2,a,v],
Boxed->False,
Background->GrayLevel[0],
ViewPoint->{0,-.5,2}];

The “rhombic spirallohedra” constructed in this notebook were discovered accidentally, in early 1997, while attempting to develop a function which dissected a polar zonohedron into its component hexahedra cyclically, that is, using nested loops in a way completely analogous to my method of creating the rhombic faces of a polar zonohedron. I managed to construct the “outermost” hexahedra in this way, such that all faces of the bounding polar zonohedron were accounted for. To all outward appearances, such an array of hexahedra looked like a polar zonohedron; but I knew the truth: it was hollow. Wondering what the interior might look like, I removed all hexahedra above the “equator” of the array, and looked inside.

(*Make the hemisphere of hexahedra*)
(*n even, n>7*)
n=12;
pitch=35.2643;
hemi=If[EvenQ[n],Floor[n/2]-1,Floor[n/2]];
v=vectors[n,pitch];

q=Table[mark={0,0,0};
	Table[hex[{ v[[j]], v[[1+Mod[j,n]]], v[[1+Mod[i+j,n]]]},mark],
		{i,n-hemi}],{j,1,n-1,2}];
(*Change "{i,n-hemi}" to "{i,n-2}" to make it complete*)


Show[Graphics3D[ Polygon /@ Flatten[q ,2] ],
Boxed->False,
Axes->False,
ViewPoint->{0,-.25, 3}];

I saw a series of spiraling ridges and valleys, and immediately realized they implied the existence of a “hidden polyhedron,” composed by a subset of the rhombic hexahedra, which polyhedron would have spiraling valleys and ridges. I set out to construct this hidden polyhedron by its faces, not by its hexahedra, and after a time, succeeded. I was able to construct what I would now call (n/2)-armed spirallohedra.

Eventually, more types of spirallohedra were recognized, and the common elements in their construction were fused to make a Spirallohedron function. Spirallohedra, as I conceive them, only exist when n>=6, and n not prime. All k-armed spirallohedra can be thought to arise by subtracting sets of rhombic hexahedra from the exterior of a polar zonohedron, said zonohedron having been dissected into hexahedra in a special, cyclic way. This construction is effected in the Volume II section, below.

Rhombic spirallohedra exist in enantiomorphic pairs, depending upon the sense of their spiraling ridges and valleys.

(*An n=12, 3-armed spirallohedron*)
Show[a=Spirallohedron[12,3],
Axes->True,
ViewPoint->2.25*unit[ {0,0,2} ]];

(*The enantiomorphic spirallohedron*)
Show[Graphics3D[ Map[{-1,-1,-1}*#&, a[[1]], {-2}] ],
Axes->True,
ViewPoint->2.25*unit[ {0,0,2} ]];

I have found that for j=n/k, rhombic spirallohedra are bounded by n(n-j) faces, 2n(n-j) edges, and n*(n-j)+2 vertices. Using my parse function (which takes a Graphics3D polyhedron and reduces it to a list of Head-less polygons) and my vunion function (which finds the union of a set of n-vectors) the number of faces and number of vertices can be verified.

n=12;k=3;j=(n/k);
a=Spirallohedron[n,k];
Print["For n=",n,", k=",k,", j=",j,", there are ", Length[parse[a]], " faces."](*number of faces*)
Print["For n=",n,", k=",k,", j=",j, ", there are ",
Length[vunion[Flatten[parse[a],1]]], " vertices."](*number of vertices*)

The complete syntax for the Spirallohedron function is
           Spirallohedron[n, k, pitch, factor, ratio], where
           
1.  n is the (integer) number of directions in which edges occur.
2.  k is a factor of n, such that k>2 and k<n.
3.  pitch is the angle, in degrees (0<pitch<90), between edges and a plane orthogonal to the symmetry axis.
4.  factor is an integer factor of n, 1<=factor<=n. Ordinarily, if this option is used, 1<factor<n.
5.  ratio is a non-negative real number, by which every factor-th vector is multiplied.

By default, the Spirallohedron function uses n=12,  k=3, pitch=35.2643, factor=1, and ratio=1. Hence one may simply write:

n=24;k=4;
Show[Spirallohedron[n,k],
ViewPoint->{0,0,2}];

Non-equilateral equivalents to the k-armed spirallohedra may be constructed using factor and ratio, so that every factor-th vector is multiplied by ratio. Note that the same spirallohedron is shown in top and bottom views; yet the faces surrounding the axis of k-fold symmetry, in the center of each image, are different. Although every face of a spirallohedron has a counter-face, faces and counter-faces are not similarly situated. Thus, spirallohedra probably should not be regarded as zonohedra.

Show[Spirallohedron[12,6,35.2643,2,4],
ViewPoint->{0,0,2.5}];(*top view*)
Show[Spirallohedron[12,6,35.2643,2,4],
ViewPoint->{0,0,-2.5}];(*bottom view*)

One of the forms below (when k equals 6) is what results if, as mentioned above, a shell of rhombic hexahedra is created which accounts for all exterior faces of the polar zonohedron, so that the interior faces of this shell define a “hidden” zonohedron. The routine fails for odd n and when n=4 (no “hidden” zonohedron exists). Perhaps the most striking forms arise when k equals 3; I call these  “3-armed spirallohedra.”

Show[Spirallohedron[12,6],(*6-armed*)
Boxed->False,
ViewPoint->{0,0,3}];
Show[Spirallohedron[12,4],(*4-armed*)
Boxed->False,
ViewPoint->{0,0,3}];
Show[Spirallohedron[12,3],(*3-armed*)
Boxed->False,
ViewPoint->{0,0,3}];

The Spirallohedron function works somewhat as follows. The iterated routine is a table of tables of tables, with complicated and opaque constraints applied to the loops. Within these tables the quad function is invoked, which creates a parallelogram from a list of two vectors. The vectors are created in the vectors function, which accepts four arguments, n, pitch, factor, and ratio. All the equilateral spirallohedra created represent types of “hidden” zonohedra which can arise by taking away certain sets of rhombic hexahedra from a polar zonohedron.

It is quite possible that my Spirallohedron function can be simplified, perhaps reduced to two nested loops, rather than three.

To construct a k-armed spirallohedron, two constraints must be met: k must divide n without remainder, and k must be greater than or equal to 3. Hence n>=6, and n not prime. A method to discover all possible values for k, for various n, is as follows:

t=Table[
a=Table[{i,Mod[n,i]},{i,2,(n/2)-1}];
b=Select[a,#[[2]]==0&];
facs=If[b!={},Transpose[b][[1]],{}];
Join[{n}, Map[n/#&, facs]],(*the list of k's*)
{n, 3, 48, 1}];
t1=Select[t, Length[#]>1&](*discard empty sets, output->{n, k1, k2, ...}*)

Spirallohedra do not have central symmetry, and one doubts whether they should be called zonohedra. Although every face has a counter-face, they are not similarly situated, and although we may pick out a zone of faces easily enough, the lack of symmetry is evident. Later, in the sections on close-packing, we will find that the faces of k-armed spirallohedra are distributed in patches, k(k-1) patches altogether, and that every patch has an equal and opposite counterpatch. It should follow that the zones of faces are distributed in families, such that one family of zones crosses another family at a patch and counter-patch.

n=12;k=3;pitch=35.2643;v=vectors[n,pitch];
a=Spirallohedron[n,k,pitch];

(*color the second zone differently*)
Show[colorzone[2,a,v],
LightSources -> 
   {{{-1, 1, 1}, RGBColor[1, 1, 1]},
   {{1, -2, 2}, RGBColor[1, 1, 1]}},
Boxed->False,
Background->GrayLevel[0],
ViewPoint->{0,0,-2.5}];

(*show the second zone by itself*)
Show[onezone[2,a,v],
Boxed->False,
Background->GrayLevel[0],
ViewPoint->{0,0,2.5}];

(*show all faces not in the second zone*)
Show[notzone[2,a,v],
Boxed->False,
Background->GrayLevel[0],
ViewPoint->{0,0,2.5}];

Initialization

Miscellaneous; Polar zonohedra; Rhombic Spirallohedra

mag[v_]:=Sqrt[v.v](*magnitude of vector 'v'*)

unit[v_]:=v/mag[v](*unit vector in direction of 'v'*)

cent[v_List]:=Plus@@v/Length[v](*centroid of point list*)

mycross[ {ax_, ay_, az_}, {bx_, by_, bz_} ] :=
	{ay bz - az by, az bx - ax bz, ax by - ay bx}(*cross product*)
	
tosegs[gon_]:=Transpose[ {gon, RotateLeft[gon]} ](*polygon to edges*)

(*area of parallelogram*)
sz[f_List] := mag[mycross[ f[[2]] - f[[1]], f[[3]] - f[[2]] ]]

(* treat 4gons specially for efficiency *)
totalarea[v_List] /; Length[v] == 4 :=
	sz[v]

(* treat triangles specially for efficiency *)
totalarea[v_List] /; Length[v] == 3 :=
	sz[v]/2
	
(* general case: use polygon center as additional vertex *)
totalarea[v_List] :=
    Module[{bary = cent[v], circ, t},
        circ = tosegs[v];
        t=Apply[ {#1, #2, bary}&, circ, {1} ];
        Sum[sz[t[[i]]]/2,{i,Length[t]}]
    ]

(*get Headless polygons from Graphics3D object or List*)
parse[x_]:= Cases[ x, {{_,_,_}..}, Infinity]

(*sort polygons by area*)
colorSort[poly_]:=(*~30* faster than old version*)
Block[{q1,q2,p},
p=parse[poly];
q1=Table[{p[[i]],totalarea[ p[[i]] ] },{i,Length[p]}];
q2=Split[Sort[q1,#1[[2]]<#2[[2]]&],
TrueQ[(Abs[#1[[2]]-#2[[2]]]<.0001)]&];
Print["You will need at least ",Length[q2],
" colors in the colorlist below."];(*how many colors?*)
Table[Cases[q2[[i]],{{_,_,_}..},Infinity],{i,Length[q2]}]
]

(*remove coincident faces; by Xah Lee*)
dePair[poly_]:=
Block[{a},
a=Sort[((Chop[Rationalize[(cent@#),10^-3]])@#&/@ (poly)),
OrderedQ[{Head@#1,Head@#2}]&];
Nest[
	(
	(DeleteCases[#,Null,{1}]&) @ (Flatten[#,1]&) @
	((If[SameQ@@(Head/@#),Null,#]&) /@
		Partition[Join[RotateRight[#],{Null,Null,Null}],2])
	)&,a,2]/.({___}[li_List]:>li)
]

(*n vectors*)
vectors[n_Integer,pitch_,factor_Integer,ratio_]:=
Block[{e},
Table[If[Mod[i,factor]==0,e=ratio,e=1];
	N[{e Cos[Degree pitch] Cos[2Pi i/n],
	e Cos[Degree pitch] Sin[2Pi i/n],
	e Sin[Degree pitch]}],
	{i,0,n-1}]]

vectors[n_:12,pitch_:35.2643,factor_:1,ratio_:1]:=
vectors[n, pitch, factor, ratio]
	
(*parallelogram from two vectors, update global 'mark'*)
quad[x_,ref_]:={ref, ref+x[[1]], ref+x[[1]]+x[[2]], mark=ref+x[[2]]}

PolarZonohedron[n_Integer,pitch_,factor_Integer,ratio_]:=
Block[{v=vectors[n,pitch,factor,ratio],h=(Plus @@ vectors[n,pitch,factor,ratio])/2, q},
	q=Table[mark=-h;
		Table[quad[{v[[j]],v[[1+Mod[i+j-1,n]]]},mark],{i,n-1}],{j,n}];
	Graphics3D[ Polygon /@ Flatten[q ,1] ]
]

PolarZonohedron[n_:12,pitch_:35.2643,factor_:1,ratio_:1]:=
PolarZonohedron[n, pitch, factor, ratio]

Spirallohedron[n_Integer,karm_Integer,pitch_,factor_,ratio_]:=
Block[{v=vectors[n,pitch,factor,ratio],
h=((Plus @@ vectors[n,pitch,factor,ratio])[[3]])/2,
klim,q,t},
klim=(n/karm);
q=Table[Table[If[k==1,mark={0,0,-h},mark=({0,0,-h}+Sum[ v[[ 1+Mod[j+(t-2),n] ]],{t,k-1}])];
Table[quad[{v[[1+Mod[j+(k-2),n]]],v[[1+Mod[i+j+(k-2),n]]]},mark],
	{i,klim-(k-1),n-k}],
		{j,1,n-1,klim}],
			{k,klim}];
Graphics3D[ Polygon /@ Flatten[q ,2] ]
]

Spirallohedron[n_:12, karm_:3, pitch_:35.2643, factor_:1, ratio_:1]:=
Spirallohedron[n, karm, pitch, factor, ratio]

(*from Graphics:Shapes: translate a polyhedron, etc.*)
trans[shape_, vec_List] :=
Block[{t = N[vec]},
 shape /. {
  poly:Polygon[_] :> Map[(t + #)&, poly, {2}],
  ln:Line[_]    :> Map[(t + #)&, ln, {2}],
  pt:Point[_]  :> Map[(t + #)&, pt,{1}] }
    ]

vunion::usage = "vunion[v] returns the union of a list of n-vectors.";

vunion[v_]:=
Union[Chop[v],
SameTest->((Sign[#1]==Sign[#2]) && (Chop[#1-#2]==Table[0,{Length[v[[1]]]}]) &)]

(*makes an array of ijk points, translations for close-packing*)
makearray[{inum_,jnum_,knum_},trans_]:=
Block[{ir,jr,kr,f1,f2,f3},
ir=(inum-1)/2;
jr=(jnum-1)/2;
kr=(knum-1)/2;
Flatten[Table[
f1=i*trans[[1]];
f2=j*trans[[2]]+f1;
f3=k*trans[[3]]+f2;
f3,
{i,-ir,ir,1},
{j,-jr,jr,1},
{k,-kr,kr,1}
],2]
]

(*makes an array of ijk points, translations for close-packing*)
make2Darray[{inum_,jnum_},trans_]:=
Block[{ir,jr,f1,f2},
ir=(inum-1)/2;
jr=(jnum-1)/2;
Flatten[Table[
f1=i*trans[[1]];
f2=j*trans[[2]]+f1;
f2,
{i,-ir,ir,1},
{j,-jr,jr,1}
],1]
]

tol=0.000001;
same[v1_, v2_]:=
   And @@ Map[Abs[#]<tol&, v1-v2]
   
facevectors[gon_]:= Map[#[[1]]-#[[2]]&,Partition[Join[gon,{First[gon]}],2,1]]

(*get kth zone from zonohedron 'poly' determined by vectors 'star'*)
onezone[k_Integer,poly_,star_]:= 
Block[{qaz=Map[facevectors[#]&,b=parse[poly]]},
zaq=Map[First,Position[Map[Map[same[#, star[[k]] ]&,# ]&,qaz],True]];
Graphics3D[ Polygon /@ b[[zaq]] ] 
]

(*get all faces not in kth zone from zonohedron 'poly'
determined by vectors 'star'*)
notzone[k_Integer,poly_,star_]:= 
Block[{qaz=Map[facevectors[#]&,b=parse[poly]]},
zaq=Union[Map[First,Position[Map[Map[same[#, star[[k]] ]&,# ]&,qaz],True]]];
xaq=Select[Range[Length[b]],MemberQ[zaq,#]==False&];
Graphics3D[ Polygon /@ b[[xaq]] ] 
]

(*color kth zone from zonohedron 'poly' determined by vectors 'star'*)
colorzone[k_Integer,poly_,star_]:= 
Block[{qaz=Map[facevectors[#]&,b=parse[poly]]},
zaq=Union[Map[First,Position[Map[Map[same[#, star[[k]] ]&,# ]&,qaz],True]]];
xaq=Select[Range[Length[b]],MemberQ[zaq,#]==False&];
Graphics3D[ {{SurfaceColor[RGBColor[1,.3,.3] ],Polygon /@ b[[zaq]]},
{SurfaceColor[RGBColor[.5,.5,1] ],Polygon /@ b[[xaq]]}} ] 
]

(*a function which creates a rhombic hexahedron from three vectors*)
(*it then updates global variable 'mark'*)
hex[x_,ref_]:={
{a=ref, b=ref+x[[1]], c=ref+x[[1]]+x[[2]], d=ref+x[[2]]},
{a,a+x[[3]],b+x[[3]],b},
{b,b+x[[3]],c+x[[3]],c},
{c,c+x[[3]],d+x[[3]],d},
{d,d+x[[3]],a+x[[3]],a},
{mark=a+x[[3]],b+x[[3]],c+x[[3]],d+x[[3]]}
}

For Ruled Surfaces

(*makes lines along ridges of (n/factor)-armed spirallohedron*)
ridges[n_Integer,kk_Integer,pitch_]:=
Block[{v=vectors[n,pitch],h=(Plus @@ vectors[n,pitch])/2, q, factor},
factor=(n/kk);
	q=Table[mark=-h;
		Table[zad[ v[[1+Mod[i+j-1,n]]], mark],{i,0,n-1}],{j,1,n,factor}];
	Graphics3D[ {AbsoluteThickness[2],Line /@ Flatten[q ,1]} ]
]

ridges[n_:12,kk_:1,pitch_:35.2643]:=
ridges[n, kk, pitch]

(*updates position on ridge*)
zad[x_,ref_]:={ref, mark=ref+x}

(*assembles list of vertices of (n/factor)-armed spirallohedron for ruling*)
spiverts[n_Integer,kk_Integer,pitch_]:=
Block[{v=vectors[n,pitch],h=(Plus @@ vectors[n,pitch])/2, q, factor},
factor=(n/kk);
	q=Table[mark=-h;
		Table[zad[ v[[1+Mod[i+j-1,n]]], mark],{i,0,n-1}],{j,1,n,factor}];
	Chop[q]
]

spiverts[n_:12,kk_:1,pitch_:35.2643]:=
spiverts[n, kk, pitch]

RSruled[x_]:=(*uses output of spiverts[]*)
Block[{b,c},
b=Transpose[x];
Graphics3D[
Table[c=Last /@ b[[k]];{RGBColor[0,0,1], Line[ Join[ c, {First@c}] ]},
{k, Length[b]}]
](*end Graphics3D*)
]

RSgons[x_]:=(*uses output of spiverts[]*)
Map[#[[2]]&, Transpose[x], {2}]

Sections, volume

(*divide-by-zero-protected*)
(*find the section, if it exists, of a segment '{p, q}'
by the subspace perpendicular to the unit vector 'vec'
at distance 'd' from the origin*)

section[{p_, q_}, d_, vec_] :=
Block[{t,tempdot,denom},
tempdot=(q - p) . vec;
denom=If[Chop[tempdot]==0,.000000001,tempdot];
t = (d*vec . vec - p . vec)/denom; 
If[t >= 0 && t <= 1, p + t*(q - p)]
]

(*find the distance from the point 'pt' to the plane
(or n-1)-space) with unit normal 'norm' at distance 'dist'
from the origin.*)

ptspace[pt_,norm_,dist_]:=Abs[(Plus @@ (pt*norm))+dist]

hullvol[poly_]:=(*changed, to hopefully handle non-centered objects*)
Block[{a,c,transpoly},
c=Chop[cent[ vunion[Flatten[parse[poly],1]] ]];
transpoly=trans[poly, -c];
(*first divide faces into triangles*)
a=Flatten[Map[toTriangles[#]&, parse[transpoly] ],1];
(*sum one-sixth the volumes of the determined parallelpipeds*)
N[Plus @@ Flatten[Map[(Abs[Minors[#, 3]]/6)&, a]]]
]

(*from Wouter Meussen*)
myOrdering[z_List]:=(Sort[Transpose[{z,Range[Length[z]]}]]//Transpose) [[2]]

(*from Wouter Meussen*)
ang2[{v1_List,v2_List},v2_]:=1.
ang2[{v1_List,v2_List},v1_]:=-1.
ang2[{v1_List,v2_List},v3_List]:=
Block[{a,b},
a=If[mag[v2-v1]==0,.000000001,mag[v2-v1]];
b=If[mag[v3-v2]==0,.000000001,mag[v3-v2]];
(v2-v1).(v3-v2)/a/b
]

(*from Wouter Meussen*)
ang[{v1_List,v2_List},v2_|v1_]:=1.
ang[{v1_List,v2_List},v3_List]:=
Block[{a,b},
a=If[mag[v2-v1]==0,.000000001,mag[v2-v1]];
b=If[mag[v3-v2]==0,.000000001,mag[v3-v2]];
(v2-v1).(v3-v2)/a/b
]

(*from Wouter Meussen*)
(*finds convex polygonal hull of 3D points in plane*)
convexhull[pts_]:=
Block[{cen,dis,idx,xfar,far,xsec,hull},
cen=cent[pts];
dis=mag[#-cen]& /@ pts;
idx=myOrdering[dis];
xfar=Last[idx];
far=pts[[xfar]];
xsec=Sort[{10+  ang2[{cen,far},pts[[#]] ]-10.,
	mag[(far-pts[[#]])] ,#} &/@idx][[-2,-1]];
hull={xfar,xsec};

FixedPoint[Last[AppendTo[hull, 
Sort[ {10+ang[  pts[[Take[hull,-2] ]],
pts[[#]] ]-10.,
mag[(pts[[ Last[hull] ]] -pts[[#]])] ,#} &/@idx] [[-3,-1]] ] ]&,0,
SameTest->(First[hull]==#2 &)];
pts[[ Drop[hull,-1] ]]
]

(*These functions appear elsewhere under different names*)

(* turn n-gon into list of triangles *)

toTriangles[vlist_List] /; Length[vlist] == 3 := {vlist}

(* treat 4gons specially for efficiency *)

toTriangles[vlist_List] /; Length[vlist] == 4 :=
	{vlist[[{1, 2, 3}]], vlist[[{3, 4, 1}]]}

(* general case: use center of gravity as additional vertex *)
(* in this way, some nonconvex polygons can also be rendered correctly *)

toTriangles[vlist_List] :=
    Module[{bary = (Plus @@ vlist)/Length[vlist], circ},
        circ = Partition[ Append[vlist, First[vlist]], 2, 1 ];
        Apply[ {#1, #2, bary}&, circ, {1} ]
    ]

Volume II

(*a function which creates a rhombic hexahedron from three vectors*)
(*it then updates global variable 'mark'*)
newhex[{i_,j_,k_},ref_]:=Graphics3D[
Polygon/@
{
{ref, ref+i, ref+i+j, ref+j},(*face ij*)
{ref+j,ref+j+i,ref+j+i+k,ref+j+k},(*face ik*)

{ref, ref+j, ref+j+k, mark=ref+k},(*face jk and update mark*)
{ref+k, ref+k+j, ref+k+j+i, ref+k+i},(*face ji*)

{ref, ref+k, ref+k+i, ref+i},(*face ki*)
{ref+i, ref+i+k, ref+i+k+j, ref+i+j}(*face kj*)
}]

(*create the vectors which determine a polar zonohedron*)
vectors[n_Integer,pitch_,edge_]:=
Table[
	N[{edge Cos[Degree pitch] Cos[2Pi i/n],
	edge Cos[Degree pitch] Sin[2Pi i/n],
	edge Sin[Degree pitch]}],
	{i,0,n-1}]

(*note: pitch in radians*)
kvolume[n_,k_,pitch_,edge_]:=(edge^3*n*Cos[pitch]^2*Csc[Pi/n]^2*(-(k*Sin[(2*Pi)/k]) + n*Sin[(2*Pi)/n])*Sin[pitch])/(4*k)

(*from Standard Packages, Discrete Math, Combinatorica.nb*)
KSubsets[l_List,0] := { {} }
KSubsets[l_List,1] := Partition[l,1]
KSubsets[l_List,k_Integer?Positive] := {l} /; (k == Length[l])
KSubsets[l_List,k_Integer?Positive] := {}  /; (k > Length[l])

KSubsets[l_List,k_Integer?Positive] :=
	Join[
		Map[(Prepend[#,First[l]])&, KSubsets[Rest[l],k-1]],
		KSubsets[Rest[l],k]
	]

(*makes the face-index of equivalent faces above and below the equator equal*)
cycleface[trip_]:=Map[If[#>hemi,n-#,#]&, trip]

(*take face-indices ijk and form product of sines*)
sinproduct[f_]:=
StringJoin[
If[f[[1]]==1, StringJoin[ "Sin[Pi/",ToString[n],"] " ],
	 StringJoin[ "Sin[",ToString[ f[[1]] ],"Pi/",ToString[n],"] "] ],
If[f[[2]]==1, StringJoin[ "Sin[Pi/",ToString[n],"] " ],
	 StringJoin[ "Sin[",ToString[ f[[2]] ],"Pi/",ToString[n],"] "] ],
"Sin[",ToString[ f[[3]] ],"Pi/",ToString[n],"] "
]

Clear[i,j,k,n,t,edge](*this is how I got symbolic hexahedral volume*)
jj=edge*{{Cos[t]*Cos[2*i*Pi/n],Cos[t]*Sin[2*i*Pi/n],Sin[t]},
{Cos[t]*Cos[2*j*Pi/n],Cos[t]*Sin[2*j*Pi/n],Sin[t]},
{Cos[t]*Cos[2*k*Pi/n],Cos[t]*Sin[2*k*Pi/n],Sin[t]}};

Simplify[  jj[[1]].(Cross[ jj[[2]], jj[[3]] ]) ]

(*how many hexahedra compose an {n,k} spirallohedron?*)
hexacore[n_,k_]:=((-2 + k)*(-1 + k)*n^3)/(6*k^2)

symhexvol[{i_,j_,k_}, n_,t_]:=Abs[
-4*edge^3*Cos[t]^2*Sin[((i - j)*Pi)/n]*Sin[((i - k)*Pi)/n]*Sin[((j - k)*Pi)/n]*Sin[t]
]

(*sort hexahedra by volume*)
volumeSort[trip_]:=
Block[{q1,q3},
q1=Map[{#[[1]], Det[ v[[ #[[1]] ]] ], #[[2]] }&, trip];
q3=Split[Sort[q1,#1[[2]]<#2[[2]]&],
TrueQ[(Abs[#1[[2]]-#2[[2]]]<.0001)]&];
Print["You will need at least ",Length[q3],
" colors in the colorlist below."];
Table[Cases[q3[[i]],Graphics3D[_],Infinity],{i,Length[q3]}]
];

3-armed cubic filling

The 3-armed spirallohedron arises whenever n / j = 3. It is reminiscent of a cube, and by the methods below, only three translation vectors are found. The method is somewhat as follows: the spirallohedron face-centers are found. Using the first two edges of each face, the surface normals are found. From this list we find the positions of faces with normals opposite to each given face. A vector of translation is found by subtracting the first face-center from the second. Then this list of translation vectors is reduced by taking its union (using the vunion function), and further reduced by discarding collinear vectors.There are always but three directions of translation for 3-armed spirallohedra, after discarding opposites.

There are also always six “patches” of faces in three patch-counterpatch pairs. Below, we create close-packings of i*j*k 3-armed spirallohedra, and examine their translation vectors, and their patches and counterpatches. Each patch is a “rhomboid region” of n/k by n/k rhombs.

Close-packing

(*create and preview 3-armed spirallohedron 'a'*)
n=12;k=3;
Show[a=Spirallohedron[n,k],
Axes->True,
ViewPoint->{0,0,2}];

faces=parse[a];(*get faces*)
c1=Map[ Chop[cent[#] ]&, faces];(*get face centers*)
(*now get cross products for all faces*)
c2=Map[ Chop[ Cross[ (#[[2]]-#[[1]]), (#[[3]]-#[[2]]) ] ]&, faces];

(*find indices of pairs of faces with opposite normals*)
z1=Table[
{k,Flatten[ Position[c2, x_ /; Chop[ x+c2[[k]] ]=={0,0,0},{1}, 1 ] ][[1]]},
{k,Length[c2]}];

Print["Total number of faces is ", Length[faces]]
dex=Union[ Sort /@ z1 ];(*unique pairs of indices, half number of faces*)

Print["Number of opposite pairs is ", Length[dex]];

q=Map[c1[[ #[[1]] ]]-c1[[ #[[2]] ]]&, z1];(*subtract every face center from its opposite*)

vectrans=vunion[q];(*discard duplicate translation vectors*)
Print["Total number of trans vectors is ", Length[vectrans]]

t=Union[ vectrans, SameTest->(Chop[#1+#2]=={0,0,0}&)];(*discard collinears*)
Print["Discarding opposites, number of trans vectors is ", Length[t]]

qaz=Table[(*find pairs of patches connected to the same trans vector*)
Flatten[ Position[q, x_ /; (Chop[ x+vectrans[[k]] ]=={0,0,0}) ,{1}] ],
{k,Length[vectrans]}];

Print["The faces fall into ", Length[qaz], " sets, each of length ",
Union[Map[Length[#]&, qaz]] ]

We preview the space-filling by applying each of the three translation vectors to our 3-armed spirallohedron a.

(*Preview translations*)
b1=trans[a, t[[1]] ];
b2=trans[a, t[[2]] ];
b3=trans[a, t[[3]] ];
Show[{a,b1,b2,b3},
Boxed->False,
Background->RGBColor[.4,.1,.1],
ViewPoint->2.25*unit[{0,0,2}]];
Print["There are ",Length[faces]," polygons bounding the spirallohedron."]
Print["There are ",Length[p1=parse[{a,b1,b2,b3}]]," polygons altogether."]
Print["There are ", Length[dePair[p1] ]," unshared polygons."]

Having visually confirmed that, at least, the spirallohedra do not interpenetrate, we create a list of multiples of the three translation vectors. A parallelepipedal array of points is created, the spirallohedra translated accordingly, and then, partly as a check to see if the spirallohedra are meeting face-to-face, and partly to reduce computational overload, all coincident pairs of faces are discarded. Only the exterior faces are displayed. For instance, in the example here, there are 2592 rhombs altogether, but only 864 remain after removing the coincident (interior) pairs.

latticepoints=makearray[{3,3,3},t];(*make an array of {i*j*k} points for translation*)
pack=Map[trans[a,#]&, latticepoints];(*can explode by writing "1.1*latticepoints"*)
Print["There are ",Length[f1=parse[pack]]," faces altogether."]
Print["There are ",Length[ f2=dePair[f1] ]," exterior faces."]

(*display array*)
Show[zaq=Graphics3D[{EdgeForm[],Polygon /@ f2}],
Boxed->False,
Background->RGBColor[.4,.1,.1],
ViewPoint->2.25*unit[ {0,0,1} ]];

Patches, Counterpatches, and produced translation vectors

Show all patch-counterpatch pairs in sequence.

(*reduce to (k(k-1)/2) pairs of patch and counterpatch*)
patches=Union[Sort/@Map[ Sort/@Transpose[ z1[[ # ]] ]&, qaz]];
facepatches=Map[Polygon/@faces[[#]]&, patches, {2} ];(*patch/counterpatch*)
patchnumber=Length[facepatches];

Table[(*look at all patches and counterpatches*)
Show[zaq=Graphics3D[ facepatches[[ i ]]  ],(*Range[i] builds up patches*)
ViewPoint->2.25*unit[{0,0,1}]],
{i,patchnumber}];

Show all translation vectors, produced through the face-centers.

(*make lines through face centers "c1"*)
sca=0.75;(*scale line from face to counterface*)
segs=Table[
vc=c1[[ dex[[ k,1 ]] ]] - c1[[ dex[[ k,2 ]] ]];(*face to counterface*)
Line[ {c1[[ dex[[ k,1 ]] ]]+(sca*vc), c1[[ dex[[ k,2 ]] ]]-(sca*vc)} ],
{k,Length[ dex ]}
];

Show[{a,Graphics3D[{
AbsolutePointSize[4], Map[Point[#]&, c1],
GrayLevel[0],AbsoluteDashing[{4,4}],segs}]
},
Boxed->False,
LightSources -> 
   {
   {{-10,8,2}, GrayLevel[.5]},
   {{0,0,5}, GrayLevel[.75]}
   },
Background->RGBColor[1,1,1],
PlotRange->All,
ViewPoint->2.25*unit[{0,0,-2}]];

We can specially color a patch-counterpatch pair.

chosenpatch=patches[[1]];(*our chosen patch/counterpatch*)
notpatch=Complement[ Range[Length[faces]], chosenpatch];(*all the other faces*)


(*apply differing colors to (1) patch/counterpatch,
and (2) the rest of the faces*)
xaq=Graphics3D[
{SurfaceColor[ RGBColor[1,1,0]],
   Polygon /@ faces[[ Flatten[chosenpatch,1] ]],
SurfaceColor[RGBColor[.5,.5,.5]],
   Polygon /@ faces[[ notpatch ]]
}
];


Show[xaq,
LightSources -> 
   {
   {{-1,1,2}, RGBColor[1, 1, 1]},
   {{1,-2,0}, RGBColor[1, 1, 1]}
   },
Boxed->False,
Background->GrayLevel[0],
ViewPoint->2.25*unit[ {0,-.5,2} ] ];

Or, we can assign three different colors to the three pairs of patches and counterpatches.

(*apply differing colors to all three patch/counterpatch pairs*)
xaq=Graphics3D[
{SurfaceColor[ RGBColor[1,1,0]],
   Polygon /@ faces[[ Flatten[ patches[[1]],1 ] ]],
SurfaceColor[RGBColor[1,0,0]],
   Polygon /@ faces[[ Flatten[ patches[[2]],1 ] ]],
SurfaceColor[RGBColor[0,0,1]],
   Polygon /@ faces[[ Flatten[ patches[[3]],1 ] ]]
}
];


Show[xaq,
LightSources -> 
   {
   {{-1,1,2}, RGBColor[1, 1, 1]},
   {{1,-2,0}, RGBColor[1, 1, 1]}
   },
Boxed->False,
Background->GrayLevel[0],
ViewPoint->2.25*unit[ {0,-1,2} ] ];

4-armed dodecahedral filling

The 4-armed spirallohedron arises whenever k=4. It is reminiscent of a rhombic dodecahedron, and by the methods below, six translation vectors are found. These methods are the same as those discussed in the section on 3-armed cubic filling, above.

There are also always twelve “patches” of faces, in six patch-counterpatch pairs. Below, we create close-packings of i*j*k 4-armed spirallohedra, and examine their translation vectors, and their patches and counterpatches. Each patch is a “rhomboid region” of n/k by n/k rhombs.

Close-packing

(*create and preview 4-armed spirallohedron 'a'*)
n=12;
Show[a=Spirallohedron[n,4],
Axes->True,
ViewPoint->{0,0,2}];

faces=parse[a];(*get faces*)
c1=Map[ Chop[cent[#] ]&, faces];(*get face centers*)
(*now get cross products for all faces*)
c2=Map[ Chop[ Cross[ (#[[2]]-#[[1]]), (#[[3]]-#[[2]]) ] ]&, faces];

(*find indices of pairs of faces with opposite normals*)
z1=Table[
{k,Flatten[ Position[c2, x_ /; Chop[ x+c2[[k]] ]=={0,0,0},{1}, 1 ] ][[1]]},
{k,Length[c2]}];

Print["Total number of faces is ", Length[faces]]
dex=Union[ Sort /@ z1 ];(*unique pairs of indices, half number of faces*)

Print["Number of opposite pairs is ", Length[dex]];

q=Map[c1[[ #[[1]] ]]-c1[[ #[[2]] ]]&, z1];(*subtract every face center from its opposite*)

vectrans=vunion[q];(*discard duplicate translation vectors*)
Print["Total number of trans vectors is ", Length[vectrans]]

t=Union[ vectrans, SameTest->(Chop[#1+#2]=={0,0,0}&)];(*discard collinears*)
Print["Discarding opposites, number of trans vectors is ", Length[t]]

qaz=Table[(*find pairs of patches connected to the same trans vector*)
Flatten[ Position[q, x_ /; (Chop[ x+vectrans[[k]] ]=={0,0,0}) ,{1}] ],
{k,Length[vectrans]}];

Print["The faces fall into ", Length[qaz], " sets, each of length ",
Union[Map[Length[#]&, qaz]] ]

We preview the space-filling by applying each of the three translation vectors to our 3-armed spirallohedron a.

(*Preview translations*)
b1=trans[a, t[[1]] ];
b2=trans[a, t[[2]] ];
b3=trans[a, t[[3]] ];
Show[{a,b1,b2,b3},
Boxed->False,
Background->RGBColor[.4,.1,.1],
ViewPoint->2.25*unit[{0,0,2}]];
Print["There are ",Length[faces]," polygons bounding the spirallohedron."]
Print["There are ",Length[f1=parse[{a,b1,b2,b3}]]," polygons altogether."]
Print["There are ", Length[dePair[f1] ]," unshared polygons."]

Having visually confirmed that, at least, the spirallohedra do not interpenetrate, we create a list of multiples of the three translation vectors. A parallelepipedal array of points is created, the spirallohedra translated accordingly, and then, partly as a check to see if the spirallohedra are meeting face-to-face, and partly to reduce computational overload, all coincident pairs of faces are discarded. Only the exterior faces are displayed. For instance, in the example here, there are 2592 rhombs altogether, but only 864 remain after removing the coincident (interior) pairs.

latticepoints=makearray[{3,3,3},t];(*make an array of {i*j*k} points for translation*)
pack=Map[trans[a,#]&, latticepoints];(*can explode by writing "1.1*latticepoints"*)
Print["There are ",Length[f1=parse[pack]]," faces altogether."]
Print["There are ",Length[ f2=dePair[f1] ]," exterior faces."]

(*display array*)
Show[zaq=Graphics3D[{EdgeForm[],Polygon /@ f2}],
Boxed->False,
Background->RGBColor[.4,.1,.1],
ViewPoint->2.25*unit[ {0,0,1} ]];

Patches, Counterpatches, and produced translation vectors

Show all patch-counterpatch pairs in sequence.

(*reduce to (k(k-1)/2) pairs of patch and counterpatch*)
patches=Union[Sort/@Map[ Sort/@Transpose[ z1[[ # ]] ]&, qaz]];
facepatches=Map[Polygon/@faces[[#]]&, patches, {2} ];(*patch/counterpatch*)
patchnumber=Length[facepatches];

Table[(*look at all patches and counterpatches*)
Show[zaq=Graphics3D[ facepatches[[ i ]]  ],(*Range[i] builds up patches*)
ViewPoint->2.25*unit[{0,0,1}]],
{i,patchnumber}];

Show all translation vectors, produced through the face-centers.

(*make lines through face centers "c1"*)
sca=0.75;(*scale line from face to counterface*)
segs=Table[
vc=c1[[ dex[[ k,1 ]] ]] - c1[[ dex[[ k,2 ]] ]];(*face to counterface*)
Line[ {c1[[ dex[[ k,1 ]] ]]+(sca*vc), c1[[ dex[[ k,2 ]] ]]-(sca*vc)} ],
{k,Length[ dex ]}
];

Show[{a,Graphics3D[{
AbsolutePointSize[4], Map[Point[#]&, c1],
GrayLevel[0],AbsoluteDashing[{4,4}],segs}]
},
Boxed->False,
LightSources -> 
   {
   {{-10,8,2}, GrayLevel[.5]},
   {{0,0,5}, GrayLevel[.75]}
   },
Background->RGBColor[1,1,1],
PlotRange->All,
ViewPoint->2.25*unit[{0,0,-2}]];

We can specially color a patch-counterpatch pair.

chosenpatch=patches[[1]];(*our chosen patch/counterpatch*)
notpatch=Complement[ Range[Length[faces]], chosenpatch];(*all the other faces*)


(*apply differing colors to (1) patch/counterpatch,
and (2) the rest of the faces*)
xaq=Graphics3D[
{SurfaceColor[ RGBColor[1,1,0]],
   Polygon /@ faces[[ Flatten[chosenpatch,1] ]],
SurfaceColor[RGBColor[.5,.5,.5]],
   Polygon /@ faces[[ notpatch ]]
}
];


Show[xaq,
LightSources -> 
   {
   {{-1,1,2}, RGBColor[1, 1, 1]},
   {{1,-2,0}, RGBColor[1, 1, 1]}
   },
Boxed->False,
Background->GrayLevel[0],
ViewPoint->2.25*unit[ {0,-.5,2} ] ];

Or, we can assign six different colors to the six pairs of patches and counterpatches.

(*apply differing colors to all three patch/counterpatch pairs*)
xaq=Graphics3D[
{SurfaceColor[ RGBColor[1,1,0]],
   Polygon /@ faces[[ Flatten[ patches[[1]],1 ] ]],
SurfaceColor[RGBColor[1,0,0]],
   Polygon /@ faces[[ Flatten[ patches[[2]],1 ] ]],
SurfaceColor[RGBColor[0,0,1]],
   Polygon /@ faces[[ Flatten[ patches[[3]],1 ] ]],
SurfaceColor[ RGBColor[1,0,0]],
   Polygon /@ faces[[ Flatten[ patches[[4]],1 ] ]],
SurfaceColor[RGBColor[0,0,.5]],
   Polygon /@ faces[[ Flatten[ patches[[5]],1 ] ]],
SurfaceColor[RGBColor[0,.5,0]],
   Polygon /@ faces[[ Flatten[ patches[[6]],1 ] ]]
}
];


Show[xaq,
LightSources -> 
   {
   {{-1,1,2}, RGBColor[1, 1, 1]},
   {{1,-2,0}, RGBColor[1, 1, 1]}
   },
Boxed->False,
Background->GrayLevel[0],
ViewPoint->2.25*unit[ {0,0,2} ] ];

Combined 3- and 6-armed filling

This packing of 3- and 6-armed spirallohedra was discovered by Frans Snik, physics student in the Netherlands. It took me a long time to construct it. After some failures I even doubted it could exist, and believed the 3-armed spirallohedra must interpenetrate.

The 6-armed spirallohedra fit together along their “equatorial” face-patches, in parallel layers, each layer constructed in the manner of hexagons tiling the plane. The 3-armed spirallohedra fit into the “hollows” where any three 6-armed spirallohedra in one plane layer meet, with the added twist, that they are in two orientations, one spirallohedron being rotated 60 degrees (equivalent to 180 degrees) with respect to the other around its axis of 3-fold symmetry. The two types occupy alternate hollows around any one 6-armed spirallohedron. Thus there are alternating “plane” layers of 6-armed and 3-armed spirallohedra.

Interestingly, this packing arises for all n a multiple of 6, including the case when n and k equal 6, and the “degenerate” 6-armed spirallohedron arises, which is an ordinary convex polar zonohedron bounded by 30 faces.

We begin by creating the two 3-armed spirallohedra; by reversing signs on only {x,y} of {x,y,z}, we effect a rotation of 180 degrees, equivalent to a rotation of 60 degrees. We name the two a1 and a2. The 6-armed spirallohedron is named b. If you wish to create the n=6 polar zonohedron, set n=6 and substitute PolarZonohedron[n] for Spirallohedron[n,k].

(*set 'n'*)
n=12;

(*create and preview 3-armed spirallohedron 'a1'*)
k=3;
Show[a1=Spirallohedron[n,k],
Axes->True,
ViewPoint->2.25*unit[ {0,0,1} ]];

(*create and preview the rotated 3-armed spirallohedron 'a2'*)
Show[a2=Graphics3D[ Map[{-1,-1, 1}*#&, a1[[1]], {-2}] ],
Axes->True,
ViewPoint->2.25*unit[ {0,0,1} ]];

(*create and preview the 6-armed spirallohedron 'b'*)
k=6;
Show[b=Spirallohedron[n,k],
Axes->True,
ViewPoint->2.25*unit[ {0,0,1} ]];

Now we find the translation vectors for the 6-armed spirallohedron b.

faces=parse[b];(*get faces*)
c1=Map[ Chop[cent[#] ]&, faces];(*get face centers*)
(*now get cross products for all faces*)
c2=Map[ Chop[ Cross[ (#[[2]]-#[[1]]), (#[[3]]-#[[2]]) ] ]&, faces];

(*find indices of pairs of faces with opposite normals*)
z1=Table[
{k,Flatten[ Position[c2, x_ /; Chop[ x+c2[[k]] ]=={0,0,0},{1}, 1 ] ][[1]]},
{k,Length[c2]}];

Print["Total number of faces is ", Length[faces]]
dex=Union[ Sort /@ z1 ];(*unique pairs of indices, half number of faces*)

Print["Number of opposite pairs is ", Length[dex]];

q=Map[c1[[ #[[1]] ]]-c1[[ #[[2]] ]]&, z1];(*subtract every face center from its opposite*)

vectrans=vunion[q];(*discard duplicate translation vectors*)
Print["Total number of trans vectors is ", Length[vectrans]]

t=Union[ vectrans, SameTest->(Chop[#1+#2]=={0,0,0}&)];(*discard collinears*)
Print["Discarding opposites, number of trans vectors is ", Length[t]]

qaz=Table[(*find pairs of patches connected to the same trans vector*)
Flatten[ Position[q, x_ /; (Chop[ x+vectrans[[k]] ]=={0,0,0}) ,{1}] ],
{k,Length[vectrans]}];

Print["The faces fall into ", Length[qaz], " sets, each of length ",
Union[Map[Length[#]&, qaz]] ]

We select those translation vectors in the xy plane.

(*select those translation vectors with z=0*)
t2=Select[Chop[t], #[[3]]==0&]

Now, create a planar array of 6-armed spirallohedra.Only two of the three vectors t2 are used; the function make2Darray makes a 3X3 rhombic grid of points, which we use to translate the spirallohedron. We choose 3X3 so that one 6-armed spirallohedron falls on the center of the array, and we can use its “corners” or extreme equatorial vertices for the positions of the prototypical 3-armed spirallohedra. It’s just a lazy way to do it. Someday I'll fix it.

latticepoints=make2Darray[{3,3},t2];(*make an array of {i*j} points for translation*)
pack=Map[trans[b,#]&, latticepoints];(*can explode by writing "1.1*latticepoints"*)
Print["There are ",Length[f1=parse[pack]]," faces altogether."]
Print["There are ",Length[ f2=dePair[f1] ]," exterior faces."]
zaq=Graphics3D[Polygon /@ f2];

(*display array*)
Show[{zaq},
Boxed->True,
Axes->True,
ViewPoint->2.25*unit[ {0,-0,1} ]];

Next we find the six extreme equatorial points, and order them into the hexagon poly, so that we can use points {1,3,5} for one type of 3-armed spirallohedron, and points {2,4,6} for the other. To actually move the 3-armed spirallohedra around on their own translation vectors, we define hexvectors.

(*all points on equator of 6-armed "b"*)
spec=vunion[Select[Chop[Flatten[parse[b],1] ], #[[3]]==0&]]
Print["There are ", Length[spec], " points on the equator."]
hexdist=Max[Map[mag[#]&, spec]];(*greatest distance->points on proper hexagon*)
Print["The greatest distance is ", hexdist]

(*points on hexagon*)
newt=Select[spec, (Abs[mag[#]-hexdist]<.001)&];(*find points*)
poly=convexhull[newt];(*order points*)
(*now get the vectors of translation*)
hexvectors={poly[[3]]-poly[[1]], poly[[4]]-poly[[2]], poly[[5]]-poly[[3]]};

(*look at hexagon, ordered vertices*)
r=hexdist*1.3;(*for PlotRange*)
Show[Graphics3D[ Table[{Point[ poly[[i]] ],Text[ToString[i],poly[[i]]*1.1]} ,{i,Length[poly]}] ],
Boxed->True,
Axes->True,
PlotRange->{{-r,r},{-r,r},{-1,1}},
ViewPoint->2.25*unit[ {0,-0,1} ]];

Here we actually move the two types of 3-armed spirallohedra into their “initial” position. They form a ring of six.
I am not quite sure why the order {a1, a2} always works. Maybe it doesn't always work.
We also check the difference between the total number of faces, and the number remaining after all coincident faces have been removed, using dePair.

(*define the rise of one edge at default pitch=35.2643 degrees*)
h1=N[Sqrt[1/3]];
(*OK ... put them all in position first*)
(*this really ought to fail some of the time! orientation ...*)
obj1=trans[a1, poly[[1]]+{0,0,(n/2)*h1}];
obj2=trans[a2, poly[[2]]+{0,0,(n/2)*h1}];
obj3=trans[a1, poly[[3]]+{0,0,(n/2)*h1}];
obj4=trans[a2, poly[[4]]+{0,0,(n/2)*h1}];
obj5=trans[a1, poly[[5]]+{0,0,(n/2)*h1}];
obj6=trans[a2, poly[[6]]+{0,0,(n/2)*h1}];
Print["There are ",Length[ff1=parse[{obj1,obj2,obj3,obj4,obj5,obj6}]]," faces altogether."]
Print["There are ",Length[ ff2=dePair[ff1] ]," exterior faces."]

(*display the ring of 3-armed spirallohedra*)
Show[Graphics3D[Polygon /@ ff2],
Boxed->True,
Axes->True,
ViewPoint->2.25*unit[ {0,0,2} ]];

Test the packing by fitting the ring of 3-armed spirallohedra to two 6-armed spirallohedra, one above, one below.

(*translate one 6-armed to higher position*)
bhi=trans[b, {0,0,n*h1}];

(*display ring of 3-armed with one 6-armed above, one below*)
Show[{obj1,obj2,obj3,obj4,obj5,obj6,b,bhi},
Boxed->True,
Axes->True,
ViewPoint->2.25*unit[ {0,0,2} ]];

Test the packing by fitting just two 3-armed spirallohedra together, an a1 to an a2, and finding the number of coincident faces. Interestingly, a1 and a2 do not meet along a common “patch” of faces, but share two half-patches.

Show[{obj1,obj2},
Boxed->True,
Axes->True,
ViewPoint->2.25*unit[ {0,0,2} ]];
Print["There are ",Length[ff1=parse[{obj1,obj2}]]," faces altogether."]
Print["There are ",Length[ ff2=dePair[ff1] ]," exterior faces."]
Print["There are ", Length[ff1]-Length[ff2]," coincident faces."]

Test the packing by fitting just one 3-armed spirallohedron to one 6-armed, and finding the number of coincident faces. The patches of the 6-armed are (when n=12) made of four rhombs each, while those of the 3-armed are made of 16 rhombs each. We examine the projection closely, and conclude that, when n=12, two half-patches of a 3-armed spirallohedron meet four full patches of a 6-armed spirallohedron.

Show[{obj1,b},
Boxed->True,
Axes->True,
ViewPoint->2.25*unit[ {0,2,1} ]];
Print["There are ",Length[ff1=parse[{obj1,b}]]," faces altogether."]
Print["There are ",Length[ ff2=dePair[ff1] ]," exterior faces."]
Print["There are ", Length[ff1]-Length[ff2]," coincident faces."]

At last we go ahead and make a rhombic array of 3-armed spirallohedra. We do this by first translating an a1 from its “prototypical” position at one corner of the equatorial hexagon of the “central 6-armed spirallohedron,” along a line, in both directions. We make an linear array named arr1 of length tr+1. Then we translate this linear array by the second vector to make a rhombic array arr2.

The same process is followed with the rotated 3-armed spirallohedron a2, and we obtain the rhombic array arr4.

Then coincident faces are removed (Mathematica often crashes, if too many coincident polygons are passed to display), and this finall array, triarr, is shown.

tr=2;

arr1=Table[
trans[obj1, (k*hexvectors[[1]])],
{k,-tr,tr,1}];
arr2=Table[
trans[arr1, (k*hexvectors[[2]])],
{k,-tr,tr,1}];

arr3=Table[
trans[obj2, (k*hexvectors[[1]])],
{k,-tr,tr,1}];
arr4=Table[
trans[arr3, (k*hexvectors[[2]])],
{k,-tr,tr,1}];
ff3=dePair[ parse[{arr2,arr4}] ];
triarr=Graphics3D[Polygon /@ ff3];

Show[triarr,
Boxed->True,
Axes->True,
ViewPoint->2.25*unit[ {0,0,1} ]];

Finally, combine the rhombic array of 3-armed spirallohedra with the original array of 6-armed, and with one copy of that array, translated to the next level up.

(*translate the array of 6-armed to next level up*)
hizaq=trans[zaq, {0,0,n*h1}];

Show[{zaq,triarr,hizaq},
Boxed->True,
Axes->True,
ViewPoint->2.25*unit[ {0,0,1} ]];

5-armed icosahedral array

The 5-armed spirallohedron arises whenever k=5. It is reminiscent of an n=5 polar zonohedron, and by the methods below,ten translation vectors are found. These methods are the same as those discussed in the section on 3-armed cubic filling, above.

There are also always twenty “patches” of faces, in ten patch-counterpatch pairs. Below, we create arrays of i*j*k 5-armed spirallohedra, and examine their translation vectors, and their patches and counterpatches. Each patch is a “rhomboid region” of n/k by n/k rhombs.

Array

(*create and preview 3-armed spirallohedron 'a'*)
n=10;k=5;
Show[a=Spirallohedron[n,k],
Axes->True,
ViewPoint->{0,0,2}];

faces=parse[a];(*get faces*)
c1=Map[ Chop[cent[#] ]&, faces];(*get face centers*)
(*now get cross products for all faces*)
c2=Map[ Chop[ Cross[ (#[[2]]-#[[1]]), (#[[3]]-#[[2]]) ] ]&, faces];

(*find indices of pairs of faces with opposite normals*)
z1=Table[
{k,Flatten[ Position[c2, x_ /; Chop[ x+c2[[k]] ]=={0,0,0},{1}, 1 ] ][[1]]},
{k,Length[c2]}];

Print["Total number of faces is ", Length[faces]]
dex=Union[ Sort /@ z1 ];(*unique pairs of indices, half number of faces*)

Print["Number of opposite pairs is ", Length[dex]];

q=Map[c1[[ #[[1]] ]]-c1[[ #[[2]] ]]&, z1];(*subtract every face center from its opposite*)

vectrans=vunion[q];(*discard duplicate translation vectors*)
Print["Total number of trans vectors is ", Length[vectrans]]

t=Union[ vectrans, SameTest->(Chop[#1+#2]=={0,0,0}&)];(*discard collinears*)
Print["Discarding opposites, number of trans vectors is ", Length[t]]

qaz=Table[(*find pairs of patches connected to the same trans vector*)
Flatten[ Position[q, x_ /; (Chop[ x+vectrans[[k]] ]=={0,0,0}) ,{1}] ],
{k,Length[vectrans]}];

Print["The faces fall into ", Length[qaz], " sets, each of length ",
Union[Map[Length[#]&, qaz]] ]

We preview the space-filling by applying the first three translation vectors to our 5-armed spirallohedron a.

(*Preview translations*)
b1=trans[a, t[[1]] ];
b2=trans[a, t[[2]] ];
b3=trans[a, t[[3]] ];
Show[{a,b1,b2,b3},
Boxed->False,
Background->RGBColor[.4,.1,.1],
ViewPoint->2.25*unit[{0,0,2}]];
Print["There are ",Length[faces]," polygons bounding the spirallohedron."]
Print["There are ",Length[f1=parse[{a,b1,b2,b3}]]," polygons altogether."]
Print["There are ", Length[dePair[f1] ]," unshared polygons."]

Having visually confirmed that, at least, the spirallohedra do not interpenetrate, we create a list of multiples of the three translation vectors. A parallelepipedal array of points is created, the spirallohedra translated accordingly, and then, partly as a check to see if the spirallohedra are meeting face-to-face, and partly to reduce computational overload, all coincident pairs of faces are discarded. Only the exterior faces are displayed. For instance, in the example here, there are 2592 rhombs altogether, but only 864 remain after removing the coincident (interior) pairs.

latticepoints=makearray[{3,3,3},t];(*make an array of {i*j*k} points for translation*)
pack=Map[trans[a,#]&, latticepoints];(*can explode by writing "1.1*latticepoints"*)
Print["There are ",Length[f1=parse[pack]]," faces altogether."]
Print["There are ",Length[ f2=dePair[f1] ]," exterior faces."]

(*display array*)
Show[zaq=Graphics3D[{EdgeForm[],Polygon /@ f2}],
Boxed->False,
Background->RGBColor[.4,.1,.1],
ViewPoint->2.25*unit[ {0,0,1} ]];

Patches, Counterpatches, and produced translation vectors

Show all patch-counterpatch pairs in sequence.

(*reduce to (k(k-1)/2) pairs of patch and counterpatch*)
patches=Union[Sort/@Map[ Sort/@Transpose[ z1[[ # ]] ]&, qaz]];
facepatches=Map[Polygon/@faces[[#]]&, patches, {2} ];(*patch/counterpatch*)
patchnumber=Length[facepatches];

Table[(*look at all patches and counterpatches*)
Show[zaq=Graphics3D[ facepatches[[ i ]]  ],(*Range[i] builds up patches*)
ViewPoint->2.25*unit[{0,0,1}]],
{i,patchnumber}];

Show all translation vectors, produced through the face-centers.

(*make lines through face centers "c1"*)
sca=0.75;(*scale line from face to counterface*)
segs=Table[
vc=c1[[ dex[[ k,1 ]] ]] - c1[[ dex[[ k,2 ]] ]];(*face to counterface*)
Line[ {c1[[ dex[[ k,1 ]] ]]+(sca*vc), c1[[ dex[[ k,2 ]] ]]-(sca*vc)} ],
{k,Length[ dex ]}
];

Show[{a,Graphics3D[{
AbsolutePointSize[4], Map[Point[#]&, c1],
GrayLevel[0],AbsoluteDashing[{4,4}],segs}]
},
Boxed->False,
LightSources -> 
   {
   {{-10,8,2}, GrayLevel[.5]},
   {{0,0,5}, GrayLevel[.75]}
   },
Background->RGBColor[1,1,1],
PlotRange->All,
ViewPoint->2.25*unit[{0,0,-2}]];

We can specially color a patch-counterpatch pair.

chosenpatch=patches[[1]];(*our chosen patch/counterpatch*)
notpatch=Complement[ Range[Length[faces]], chosenpatch];(*all the other faces*)


(*apply differing colors to (1) patch/counterpatch,
and (2) the rest of the faces*)
xaq=Graphics3D[
{SurfaceColor[ RGBColor[1,1,0]],
   Polygon /@ faces[[ Flatten[chosenpatch,1] ]],
SurfaceColor[RGBColor[.5,.5,.5]],
   Polygon /@ faces[[ notpatch ]]
}
];


Show[xaq,
LightSources -> 
   {
   {{-1,1,2}, RGBColor[1, 1, 1]},
   {{1,-2,0}, RGBColor[1, 1, 1]}
   },
Boxed->False,
Background->GrayLevel[0],
ViewPoint->2.25*unit[ {0,-.5,2} ] ];

Or, we can assign ten different colors to the ten pairs of patches and counterpatches.

(*apply differing colors to all three patch/counterpatch pairs*)
xaq=Graphics3D[
{SurfaceColor[ RGBColor[1,1,0]],
   Polygon /@ faces[[ Flatten[ patches[[1]],1 ] ]],
SurfaceColor[RGBColor[1,0,0]],
   Polygon /@ faces[[ Flatten[ patches[[2]],1 ] ]],
SurfaceColor[RGBColor[0,0,1]],
   Polygon /@ faces[[ Flatten[ patches[[3]],1 ] ]],
SurfaceColor[ RGBColor[1,0,0]],
   Polygon /@ faces[[ Flatten[ patches[[4]],1 ] ]],
SurfaceColor[RGBColor[0,0,.5]],
   Polygon /@ faces[[ Flatten[ patches[[5]],1 ] ]],
SurfaceColor[RGBColor[0,.5,0]],
   Polygon /@ faces[[ Flatten[ patches[[6]],1 ] ]],
SurfaceColor[RGBColor[1,1,1]],
   Polygon /@ faces[[ Flatten[ patches[[7]],1 ] ]],
SurfaceColor[ RGBColor[.5,.5,.5]],
   Polygon /@ faces[[ Flatten[ patches[[8]],1 ] ]],
SurfaceColor[RGBColor[1,0,.5]],
   Polygon /@ faces[[ Flatten[ patches[[9]],1 ] ]],
SurfaceColor[RGBColor[0,0,.2]],
   Polygon /@ faces[[ Flatten[ patches[[10]],1 ] ]]   
}
];


Show[xaq,
LightSources -> 
   {
   {{-1,1,2}, RGBColor[1, 1, 1]},
   {{1,-2,0}, RGBColor[1, 1, 1]}
   },
Boxed->False,
Background->GrayLevel[0],
ViewPoint->2.25*unit[ {0,0,2} ] ];

Filling selected regions of space

When creating close-packings of polyhedra, it is often convenient to assemble the list of translation vectors first, then actually translate the polyhedron. A single spirallohedron can be a rather husky chunk of data, and an array of spirallohedra can place large demands upon the computer’s processor and RAM. Here we create a rather large list of translation points (10^3), but use only a subset to actually translate our spirallohedron.

(*create simplest unit 4-armed spirallohedron*)
Show[a=Spirallohedron[8,4],
Axes->True,
ViewPoint->{0,0,2}];

faces=parse[a];(*get faces*)
c1=Map[ Chop[cent[#] ]&, faces];(*get face centers*)
(*now get cross products for all faces*)
c2=Map[ Chop[ Cross[ (#[[2]]-#[[1]]), (#[[3]]-#[[2]]) ] ]&, faces];

(*find indices of pairs of faces with opposite normals*)
z1=Table[
{k,Flatten[ Position[c2, x_ /; Chop[ x+c2[[k]] ]=={0,0,0},{1}, 1 ] ][[1]]},
{k,Length[c2]}];

Print["Total number of faces is ", Length[faces]]
dex=Union[ Sort /@ z1 ];(*unique pairs of indices, half number of faces*)

Print["Number of opposite pairs is ", Length[dex]];

q=Map[c1[[ #[[1]] ]]-c1[[ #[[2]] ]]&, z1];(*subtract every face center from its opposite*)

vectrans=vunion[q];(*discard duplicate translation vectors*)
Print["Total number of trans vectors is ", Length[vectrans]]

t=Union[ vectrans, SameTest->(Chop[#1+#2]=={0,0,0}&)];(*discard collinears*)
Print["Discarding opposites, number of trans vectors is ", Length[t]]

qaz=Table[(*find pairs of patches connected to the same trans vector*)
Flatten[ Position[q, x_ /; (Chop[ x+vectrans[[k]] ]=={0,0,0}) ,{1}] ],
{k,Length[vectrans]}];

Print["The faces fall into ", Length[qaz], " sets, each of length ",
Union[Map[Length[#]&, qaz]] ]

latticepoints=makearray[{10,10,10},t[[{1,2,3}]] ];(*make an array of {i*j*k} points for translation*)
Print["There are ", Length[latticepoints]," points, with maximum distance ",
Max[Map[mag[#]&, latticepoints]]]

Make a Planar Array

(*select all points near the xy plane*)
region=Select[latticepoints, (Abs[ #[[3]] ]<1)&];
Length[region]

Or, a Spherical Array

(*make a hollow ball of spirallohedra with an opening*)
region=Select[latticepoints, (mag[#]>6 && mag[#]<10 && #[[3]]<5)&];
Length[region]

pack=Map[trans[a,#]&, region];(*can explode by writing "1.1*latticepoints"*)
Print["There are ",Length[f1=parse[pack]]," faces altogether."]
Print["There are ",Length[ f2=dePair[f1] ]," exterior faces."]

(*display array*)
Show[zaq=Graphics3D[{EdgeForm[],Polygon /@ f2}],
Boxed->False,
Background->RGBColor[.4,.1,.1],
ViewPoint->2.25*unit[ {0,-.25,1} ]];

Conchoids and Holy Conchoids

I call them “conchoids” because their curves remind me of shells.

If collinear translation vectors are allowed, a three-armed spirallohedron will always yield a set of six vectors. Translations along these give a packing of spirallohedra with open “cups” at the north and south poles. By filling these cups with two more spirallohedra, and deleting interior faces, we obtain a “conchoid” polyhedron. In the example below, the conchoid is bounded by 672 rhombs. Note that it mimics a rhombic hexahedron, the eight spirallohedra approximating the positions of the eight vertices of the hexahedron.

4-armed conchoids can also be created. When k>4, the spirallohedra interpenetrate, but interesting forms result, anyway.

(*create 3-armed spirallohedron*)
n=12;
Show[a=Spirallohedron[n,3],
Axes->True,
ViewPoint->{0,0,2}];

faces=parse[a];(*get faces*)
c1=Map[ Chop[cent[#] ]&, faces];(*get face centers*)
(*now get cross products for all faces*)
c2=Map[ Chop[ Cross[ (#[[2]]-#[[1]]), (#[[3]]-#[[2]]) ] ]&, faces];

(*find indices of pairs of faces with opposite normals*)
z1=Table[
{k,Flatten[ Position[c2, x_ /; Chop[ x+c2[[k]] ]=={0,0,0},{1}, 1 ] ][[1]]},
{k,Length[c2]}];

Print["Total number of faces is ", Length[faces]]
dex=Union[ Sort /@ z1 ];(*unique pairs of indices, half number of faces*)

Print["Number of opposite pairs is ", Length[dex]];

q=Map[c1[[ #[[1]] ]]-c1[[ #[[2]] ]]&, z1];(*subtract every face center from its opposite*)

vectrans=vunion[q];(*discard duplicate translation vectors*)
Print["Total number of trans vectors is ", Length[vectrans]]

t=Union[ vectrans, SameTest->(Chop[#1+#2]=={0,0,0}&)];(*discard collinears*)
Print["Discarding opposites, number of trans vectors is ", Length[t]]

qaz=Table[(*find pairs of patches connected to the same trans vector*)
Flatten[ Position[q, x_ /; (Chop[ x+vectrans[[k]] ]=={0,0,0}) ,{1}] ],
{k,Length[vectrans]}];

Print["The faces fall into ", Length[qaz], " sets, each of length ",
Union[Map[Length[#]&, qaz]] ]

cap=Chop[Plus@@vectors[n,pitch]];(*vector for polar positions*)
allconch=Join[ Map[trans[a,#]&, vectrans], {a}, {trans[a,cap]}, {trans[a,-cap]} ];
Print["There are ", Length[ f1=parse[allconch] ], " polygons altogether."]
Print["There are ", Length[ f2=dePair[f1] ], " polygons after deleting coincident pairs."]

(*display conchoid*)
Show[zaq=Graphics3D[{EdgeForm[],Polygon /@ f2}],
Boxed->False,
Background->RGBColor[.4,.1,.1],
ViewPoint->2.25*unit[{0,-.25,0}]];

Holy Conchoid

holyconch=Join[ Map[trans[a,#]&, vectrans], {a}];
Print["There are ", Length[ f1=parse[holyconch] ], " polygons altogether."]
Print["There are ", Length[ f2=dePair[f1] ], " polygons after deleting coincident pairs."]

(*display conchoid*)
Show[zaq=Graphics3D[{EdgeForm[],Polygon /@ f2}],
Boxed->False,
Background->RGBColor[.4,.1,.1],
ViewPoint->2.25*unit[{0,-.25,2}]];

Rhombic surfaces; degenerate spirallohedra

Sometimes accidents are fruitful. In this case, the “illegal” used of k=2 leads to a twisted patch of rhombs. These enclose no volume, and cannot fill space; but they can be joined to form an infinite undulating surface. The methods seem strange, even to me, but for some reason they work.

(*create surface patch where k=2*)
Show[a=Spirallohedron[12,2],
Axes->True,
ViewPoint->{2,0,0}];

(*find candidates for translation vectors*)
faces=parse[a];(*get faces*)

(*sort vertices by z coordinate*)
t=Sort[vunion[ Flatten[faces,1] ],#1[[3]]<#2[[3]]&];

t1=Last@t;(*largest z coordinate*)

(*now add vertex with largest z to all vertices with z=0*)
q1=(t1+#)& /@ Select[t,#[[3]]==0&];

(*test translation vectors*)
vpt=Cross[q1[[1]],t1];(*viewpoint orthogonal to the two vectors*)

b1=trans[a, 2*t1 ];(*straight up*)
b2=trans[a, q1[[1]] ];(*use first of the q1's*)

Show[{a,b1,b2},
Boxed->False,
Background->RGBColor[.4,.1,.1],
ViewPoint->2.25*unit[ vpt ]];

i=4;(*iterator: there will be ((2*i)+1)^2 patches*)
b4=Table[k*q1[[1]],{k,-i,i}];(*one line of patches*)

(*now translate the whole line, repeatedly, making a parallelogrammic region*)
b5=Flatten[Table[Map[k*2*t1+#&, b4],{k,-i,i}],1];

(*select a subregion of patches, with an optional hole, to display*)
t2=Select[b5, (mag[#]<12 (*&& mag[#]>10*))&];
Print["There will be ",Length[t2]," patches altogether."];

b6=Map[trans[a,#]&, t2];(*translate patches*)
zaq=Graphics3D[{EdgeForm[],Cases[b6,Polygon[x_],Infinity]}];

(*display selected patches*)
Show[zaq,
Boxed->False,
Background->RGBColor[.4,.1,.1],
ViewPoint->2.25*unit[ vpt ] ];

Coloring by area

An example, if one does not yet exist.

(*avoid 35.2643 pitch, it makes faces of different types have equal areas*)
Show[a=Spirallohedron[12,3,35],
Axes->True,
ViewPoint->{0,0,2}];

Sort polygons by surface area.

(*split the polygons into lists according to area*)
s=colorSort[a];

Various schemes to build a list of colors; choose one.

(*relate rgb colors to indices*)
z=Length[s];
colorlist=Table[ If[Mod[i,2]==0,
	RGBColor[ .2,.2,i/z ],
		RGBColor[ i/z,.2,.2 ]], {i,z}];

(*Or, use Hue.*)
z=Length[s];
colorlist=Table[ Hue[ i/z, 1, 1-(i/z)], {i,z}];

(*relate Hue to indices*)
z=Length[s];
colorlist=Table[ If[Mod[i,2]==0 || Mod[i,5]==0,
	Hue[ i/z, 1, 1-(i/(z))],
		Hue[ i/z, 1, (i/z)]], {i,z}];

(*if the number of colors needed is 10 or less*) colorlist = {RGBColor[1,  .25,  .25], RGBColor[ .25,  .25, 1], 
          RGBColor[1, 1, 0], RGBColor[ .5,  .5,  .5], RGBColor[1, 0, 1], 
          RGBColor[ .5, 0, 0], RGBColor[0, 0,  .5], RGBColor[ .5,  .5, 0], 
          RGBColor[ .5, 0,  .5], RGBColor[0,  .5, 0]};

(*Or, assemble a list of random colors of the proper length.*)
colorlist=Table[RGBColor[Random[Real,{0,1}],
		Random[Real,{0,1}],
		Random[Real,{0,1}]],{Length[s]}];

Display Colored Polygons

xaq=Graphics3D[
{EdgeForm[],
	Table[ {SurfaceColor[colorlist[[i]],GrayLevel[.5], 5 ],
		Polygon /@ s[[ i ]]}, {i,Length[s]}]
}];
	
Show[xaq,
LightSources -> 
   {
   {{-1,0,2}, RGBColor[1, 1, 1]},
   {{1,1,0}, RGBColor[1, 1, 1]}
   },
Boxed->False,
Background->GrayLevel[0],
ViewPoint->2.25*unit[ {-5,-2,0} ] ];

Ruled Surfaces

The curved, spiraling surfaces between the ridges of spirallohedra are approximations of ruled surfaces, which seem related to helicoids. Below, an illustration.

n=24;
pitch=35.2643;
k=4;
a=Show[
{RSruled[ spiverts[n, k, pitch] ],
ridges[n, k, pitch],
Spirallohedron[n, k, pitch]},
Boxed->False,
ViewPoint->{0,-2,0}
];

Arrays of k-gons; the Vertical Lines

3- and 4-armed pirallohedra close-pack; very well; but so also do their suites of k-gons fit together. Below, a suite of k-gons is found, the same translation vectors are found by the same methods, an array of lattice points constructed, and the k-gonal suite is translated into the lattice positions.

(*choose n, 3- or 4-armed spirallohedron*)
n=24;k=3;pitch=35.2643;

(*display spirallohedron*)
Show[a=Spirallohedron[n,k],
Axes->True,
ViewPoint->2.25*unit[ {0,0,1}]
];

(*get corners of k-gons*)
gonverts=RSgons[spiverts[n,k,pitch]];

(*Make all k-gons of ruled surface*)
allbases=Table[
Polygon[convexhull[
gonverts[[k]]
]],
{k,1,n-1}];
Show[ab=Graphics3D[ allbases ],
PlotRange->All,
ViewPoint->2.25*unit[ {0,0,1} ]
];

(*get translation vectors "t"*)

faces=parse[a];(*get faces*)
c1=Map[ Chop[cent[#] ]&, faces];(*get face centers*)
(*now get cross products for all faces*)
c2=Map[ Chop[ Cross[ (#[[2]]-#[[1]]), (#[[3]]-#[[2]]) ] ]&, faces];

(*find indices of pairs of faces with opposite normals*)
z1=Table[
{k,Flatten[ Position[c2, x_ /; Chop[ x+c2[[k]] ]=={0,0,0},{1}, 1 ] ][[1]]},
{k,Length[c2]}];

Print["Total number of faces is ", Length[faces]]
dex=Union[ Sort /@ z1 ];(*unique pairs of indices, half number of faces*)

Print["Number of opposite pairs is ", Length[dex]];

q=Map[c1[[ #[[1]] ]]-c1[[ #[[2]] ]]&, z1];(*subtract every face center from its opposite*)

vectrans=vunion[q];(*discard duplicate translation vectors*)
Print["Total number of trans vectors is ", Length[vectrans]]

t=Union[ vectrans, SameTest->(Chop[#1+#2]=={0,0,0}&)];(*discard collinears*)
Print["Discarding opposites, number of trans vectors is ", Length[t]]

Make Array, Display

(*make lattice, translate copies of k-gon suites*)
latticepoints=makearray[{5,5,5},t];(*make an array of {i*j*k} points for translation*)
pack=Map[trans[ab,#]&, latticepoints];(*can explode by writing "1.1*latticepoints"*)
Print["There are ",Length[f1=parse[pack]]," faces altogether."]

(*display array*)
Show[zaq=Graphics3D[{EdgeForm[RGBColor[1,1,0]],Polygon /@ f1}],
Boxed->False,
Background->RGBColor[0,0,.3],
ViewPoint->2.25*unit[ {0,-1,1} ]];

Truncate Array, Display

(*select all polygons with z coordinate less than 5.25*)
region=Select[f1, (Abs[ #[[1,3]] ]<5.25)&];
Length[region]

(*display array*)
Show[zaq=Graphics3D[{EdgeForm[RGBColor[1,1,0]],Polygon /@ region}],
Boxed->False,
Background->RGBColor[0,0,.3],
ViewPoint->2.25*unit[ {0,-1,1} ]];

Vertical Lines

The mysterious “vertical lines” may be seen for all n even, and all k, by setting n and k to desired values in the first cell below. These lines connect the vertices of the “equatorial” faces of the “parent” polar zonohedron. The transformation polar zonohedron->spirallohedra sends a “horizontal” cycle of n equatorial faces into k vertical lines of n/k of these same faces. I presume the same vertical lines exist for n odd, but since the equator of the parent polar zonohedron crosses two cycles of n equal faces, the vertical lines likewise will span two almost-vertical stacks of the corresponding faces in the spirallohedron.

(*choose n even*)
n=12;k=3;pitch=35;

(*get corners of k-gons*)
gonverts=RSgons[spiverts[n,k,pitch]];

(*n even*)
j=n/k;
goncount=((j/2)*(k-2));
 q2=Join[ {gonverts[[ goncount ]]},{ 
      RotateRight[ gonverts[[ (n-goncount) ]] ]} ]; 
q3=Graphics3D[{AbsoluteThickness[4],Line/@Transpose[q2]}];

(*display spirallohedron*)
Show[{a=Spirallohedron[n,k],q3},
Axes->True,
ViewPoint->2.25*unit[ {0,-2,1}]
];


(*Make all k-gons of ruled surface*)
allbases=Table[
Polygon[convexhull[
gonverts[[k]]
]],
{k,1,n-1}];
Show[{q3,ab=Graphics3D[ allbases ]},
PlotRange->All,
ViewPoint->2.25*unit[ {0,-2,1} ]
];

Plane Sections

Needs work!! I have my doubts about which of my several versions of the section function to use. The sections sometimes are not accurate. I believe that this is because my methods assume that each section of a spirallohedron is just one polygon, when in fact a section can be two polygons.

To see an example of such failure, set

a=Spirallohedron[12,3];
v1=N[unit[ {0,-1,0} ]];(*line perpendicular to plane of section*)
dist=N[1*Sqrt[3]];(*distance of sectioning plane from origin*)

Just One Spirallohedron

(*create spirallohedron 'a'*)
Show[a=Spirallohedron[12,3],
Axes->True,
ViewPoint->{0,0,2}];

v1=N[unit[ {0,0,1} ]];(*line perpendicular to plane of section*)
dist=N[1.5];(*distance of sectioning plane from origin*)

(*get the faces and truncated faces to one side of the section*)
f=Chop[Cases[a,Polygon[x_]:>x,Infinity]];
s1=Table[
e1=tosegs[ f[[k]] ];
q1=vunion[Select[Flatten[e1,1], Chop[(#.v1)-dist]>=0&]];(*all on one side of plane*)
a1=Select[Map[section[#, dist, v1 ]&, e1 ], VectorQ[#]&];
b1=vunion[ Join[q1,a1] ];(*it was necessary to discard duplicates*)
If[Length[b1]>2, Polygon[convexhull[ b1 ]], {}],
{k,Length[f]}];

(*the truncated spirallohedron*)
Show[spirallochunk=Graphics3D[ Cases[s1,Polygon[x_],Infinity] ],
PlotRange->All,
ViewPoint->2.25*unit[ -v1 ]
];

(*now get the often non-convex section itself*)
f=Chop[Cases[a,Polygon[x_]:>x,Infinity]];

segments=Cases[Table[
e1=tosegs[ f[[k]] ];
a1=vunion[ Chop[Select[Map[section[#, dist, v1 ]&, e1 ], VectorQ[#]&]] ],
{k,Length[f]}],
{{_,_,_},{_,_,_}}];

If[Length[segments]>2,

segverts=Chop[N[vunion[ Flatten[segments,1] ]]];

z1=Table[
Flatten[
	Join[
	Position[segverts, x:{_,_,_} /; Chop[x-segments[[k,1]]]=={0,0,0} ],
	Position[segverts, x:{_,_,_} /; Chop[x-segments[[k,2]]]=={0,0,0} ]
	],
1],
{k,Length[segments]}];

z2=Join[z1,Map[Reverse[#]&, z1] ];

beg=z2[[1]];
r=Table[
a1=Join[beg, Flatten[Select[z2,#[[1]]==beg[[2]] && #[[2]] != beg[[1]]&],1]];
beg=Take[a1,-2];
a1,
{Length[z1]}];

p=Join[ r[[1,{1,2}]], Map[Last[#]&, Drop[r,-2] ] ];
s2=Polygon[ Map[segverts[[#]]&, p] ] ];(*end If*)

(*split into triangles and remove all edges*)
spirallosect=Graphics3D[ {EdgeForm[], Polygon /@ toTriangles[ s2[[1]] ]}];
Show[spirallosect,
PlotRange->All,
Boxed->False,
ViewPoint->2.25*unit[ -v1 ]
];

Show[{spirallochunk,spirallosect},
PlotRange->All,
Boxed->False,
ViewPoint->2.25*unit[ v1 ]
];

An Array of  Spirallohedra

(*create spirallohedron 'a'*)
Show[a=Spirallohedron[12,3],
Axes->True,
ViewPoint->{0,0,2}];

(*find translation vectors*)

faces=Chop[a[[1]]] /.Polygon[x_]:>x;(*get faces*)
c1=Map[ Chop[cent[#] ]&, faces];(*get face centers*)
s1=Map[ tosegs[#]&, faces];(*get all edges*)

(*now find face normals*)
c2=Map[ Chop[ Cross[ (#[[1,2]]-#[[1,1]]), (#[[2,2]]-#[[2,1]]) ] ]&, s1];

(*find pairs of faces with opposite normals*)
z1=Table[
p3=Flatten[ Position[c2, x_ /; Chop[ x+c2[[k]] ]=={0,0,0},{1} ] ];
c1[[k]]-First@c1[[p3]],(*a vector of translation*)
{k,Length[c2]}];

(*discard coincident and collinear translation vectors*)
t=Union[ vunion[z1], SameTest->(Chop[Cross[#1, #2]]=={0,0,0}&)]
Length[t]

latticepoints=makearray[{2,2,2},t];(*make an array of {i*j*k} points for translation*)
pack=Map[trans[a,#]&, 1.1*latticepoints];(*can explode by writing "1.1*latticepoints"*)
Print["There are ",Length[f1=parse[pack]]," faces altogether."]
Print["There are ",Length[ f2=dePair[f1] ]," exterior faces."]

(*
Cross[ unit[ t[[1]] ], unit[ t[[2]] ] ](*use this to get section perpendicular to two translation vectors*)
*)

v1=N[ unit[Cross[ unit[ t[[1]] ], unit[ t[[3]] ] ] ] ];(*line perpendicular to plane of section*)
dist=3; (*distance of sectioning plane from origin*)

tot=Table[(*OK all this does is get all polyhedra and polygons to one side of section*)
f=Chop[Cases[pack[[j]],Polygon[x_]:>x,Infinity]];
Table[
e1=tosegs[ f[[k]] ];
q1=vunion[Select[Flatten[e1,1], Chop[(#.v1)-dist]>=0&]];(*all on one side of plane*)
a1=Select[Map[section[#, dist, v1 ]&, e1 ], VectorQ[#]&];
b1=vunion[ Join[q1,a1] ];(*it was necessary to discard duplicates*)
If[Length[b1]>2, Polygon[convexhull[ b1 ]], {}],
{k,Length[f]}],
{j,Length[pack]}];

Show[zaq=Graphics3D[ Cases[tot,Polygon[x_],Infinity] ],
PlotRange->All,
ViewPoint->2.25*unit[ -v1 ]
];

pat=Table[(*and now this intricate piece of work gets the (possibly non-convex) sections themselves*)
f=Chop[Cases[pack[[j]],Polygon[x_]:>x,Infinity]];

segments=Cases[Table[
e1=tosegs[ f[[k]] ];
a1=vunion[ Chop[Select[Map[section[#, dist, v1 ]&, e1 ], VectorQ[#]&]] ],
{k,Length[f]}],
{{_,_,_},{_,_,_}}];

If[Length[segments]>2,

segverts=Chop[N[vunion[ Flatten[segments,1] ]]];

z1=Table[
Flatten[
	Join[
	Position[segverts, x:{_,_,_} /; Chop[x-segments[[k,1]]]=={0,0,0} ],
	Position[segverts, x:{_,_,_} /; Chop[x-segments[[k,2]]]=={0,0,0} ]
	],
1],
{k,Length[segments]}];

z2=Join[z1,Map[Reverse[#]&, z1] ];

beg=z2[[1]];
r=Table[
a=Join[beg, Flatten[Select[z2,#[[1]]==beg[[2]] && #[[2]] != beg[[1]]&],1]];
beg=Take[a,-2];
a,
{Length[z1]}];

p=Join[ r[[1,{1,2}]], Map[Last[#]&, Drop[r,-2] ] ];
Polygon[ Map[segverts[[#]]&, p] ] ](*end If*),
{j,Length[pack]}];
Length[ sect=Cases[pat,Polygon[x_]] ]

(*so, split into triangles and remove all edges*)
poly=Polygon /@ Flatten[ Map[toTriangles[#]&, Cases[sect,Polygon[x_]:>x,Infinity] ], 1];

Show[Graphics3D[ {EdgeForm[], poly} ],
PlotRange->All,
Boxed->False,
ViewPoint->2.25*unit[ -v1 ]
];

Show[{Graphics3D[ {EdgeForm[], poly} ],zaq},
PlotRange->All,
Boxed->False,
ViewPoint->2.25*unit[ -v1 ]
];

(*everything below has to do with making prism-things and seems broken!*)
allsegs=Map[tosegs[#]&, Cases[sect,Polygon[x_]:>x,Infinity] ];

tops=Flatten[allsegs,1];

bots=Map[{{0,0,-10.25},{0,0,-10.25}}+#&, tops ];

lats=MapThread[Join[Reverse[#1],#2]&, {tops, bots} ];

Show[zaq=Graphics3D[
{Cases[tot, Polygon[x_],Infinity],
Polygon /@ lats,
{EdgeForm[],SurfaceColor[RGBColor[1,1,0]],
Polygon /@ Flatten[poly,1] },
{GrayLevel[0],AbsoluteThickness[1],Line /@ Flatten[allsegs,1]}
}
],
PlotRange->All,
Boxed->False,
ViewPoint->2.25*unit[ {-1,0,1.5} ]
];

Show[zaq,
Boxed->False,
Background->GrayLevel[0],
(*ViewVertical->unit[{0,-1,0}],*)
ViewPoint->2.25*unit[ {0,-2,1.5} ] ];

Volume I

Clear[n,k]

Briefly, let the volume of a k-armed spirallohedron be denoted (Subscript[V, k]), the volume of the circumscribing n-merous polar zonohedron, (Subscript[V, n]), the area of a regular k-gon, (Subscript[A, k]), and the area of the circumscribing regular n-gon, (Subscript[A, n]), such that k>2 and k is a factor of n; then the volume of the spirallohedron is

(Subscript[V, k]) = ((Subscript[A, k])/(Subscript[A, n]))·(Subscript[V, n])

Prologue

The volume of the general polar zonohedron of n directions which edges occur, edge length edge and angle-between-edges-and-a-plane-orthogonal-to-the-symmetry-axis pitch is

(1/2) (edge^3) (n^2) cot((π/n)) (cos^2)(
            pitch) sin(pitch)

What is the volume of a k-armed spirallohedron? Perhaps we should adopt the usage, an {n,k}-armed spirallohedron, since there are still n directions in which edges occur. Then it is true that an {n,k}-armed spirallohedron may always be inscribed within an {n}-polar zonohedron of the same edge length and pitch, so that certain k sets of edges coincide, and true also that for any given polar zonohedron, if one or more related k-armed spirallohedra exist, they have, at any rate, a lesser volume.

In fact, spirallohedra can be said to arise by subtracting rhombic hexahedra from the exteriors of polar zonohedra. The fewer the “arms,” the less the volume. For the special case when k equals n/2, a simple expression gives the volume of the spirallohedron:

(1/2) (edge^3) (n^2) cot((π/n)) cos((π/((n/2)))) (cos^2)(pitch) sin(pitch)

The new term is highlighted.

I began efforts to find an expression giving the volume of the general {n,k}-armed spirallohedron as follows. In the example below, a spirallohedron is constructed, its series of k-gons is also constructed, and plane sections are taken at each level in a cumbersome exercise, so that a sort of “antiprismoid” is cut out from the rest of the spirallohedron (or a right regular k-gonal pyramid, from each pole), and the appropriate k-gonal bases or base added. Then, treating these as convex polyhedra, although they are not, their volumes are found. Finally, the volume is compared to that of the related polar zonohedron, giving a ratio.

In the next subsection, results are tabulated. There were many unsolved ratios.

Example of method to measure volume

This shows the method by which I split the spirallohedra into chunks, and summed their volumes.

(*create spirallohedron 'a'*)
n=10;pitch=35.2643;j=2;k=5;
Show[a=Spirallohedron[n,k],
Axes->True,
ViewPoint->{0,0,2}];

(*get vertices*)
verts=vunion[Flatten[ parse[a], 1]];
z=Union[ Map[#[[3]]&, verts], SameTest->(Chop[#1-#2]==0&) ]
Length[z]

(*get corners of k-gons*)
gonverts=Reverse[RSgons[spiverts[n,k,pitch]]];
Length[gonverts]

(*Make all k-gons of ruled surface*)
allbases=Table[
Polygon[convexhull[
gonverts[[k]]
]],
{k,n,2,-1}];
Show[zaq=Graphics3D[ allbases ],
PlotRange->All,
ViewPoint->2.25*unit[ {0,-1,0} ]
];
Length[allbases]

v1=N[unit[ {0,0,1} ]];(*vector perpendicular to plane of section*)

layers=Table[

mydist=z[[ i ]];(*corresponds to allbases[[ stage-1 ]]*)

f=Chop[parse[a]];(*get all faces*)

t=Table[(*make a list of polygons to one side of section*)
e1=tosegs[ f[[k]] ];(*make polygon into list of segments*)

(*all on one side of plane*)
q1=vunion[Select[Flatten[e1,1], Chop[(#.v1)-mydist]>=0&]];

(*now get actual segment sections (points)*)
a1=Chop[ Select[Map[section[#, mydist, v1 ]&, e1 ], VectorQ[#]&] ];
(*Print[ Length[q1],Length[a1] ];what is going on?*)

b1=vunion[ Join[q1,a1] ];(*it was necessary to discard duplicates*)

(*avoid degenerate cases, make polygon*)
If[Length[b1]>2, Polygon[convexhull[ b1 ]], {}],
{k,Length[f]}];

poly1=(*now get non-empty, Polygon members of list*)
Cases[t,Polygon[x_],Infinity];


mydist=z[[ i+1 ]];
f=Chop[parse[poly1]];(*get all faces*)

t=Table[(*make a list of polygons to one side of section*)
e1=tosegs[ f[[k]] ];(*make polygon into list of segments*)

(*all on one side of plane*)
q1=vunion[Select[Flatten[e1,1], Chop[(#.v1)-mydist]<=0&]];

(*now get actual segment sections (points)*)
a1=Chop[ Select[Map[section[#, mydist, v1 ]&, e1 ], VectorQ[#]&] ];
(*Print[ Length[q1],Length[a1] ];what is going on?*)

b1=vunion[ Join[q1,a1] ];(*it was necessary to discard duplicates*)

(*avoid degenerate cases, make polygon*)
If[Length[b1]>2, Polygon[convexhull[ b1 ]], {}],
{k,Length[f]}];

poly2=(*now get non-empty, Polygon members of list*)
Cases[t,Polygon[x_],Infinity];

(*Show[Graphics3D[ {poly2} ],
PlotRange->All,
ViewPoint->2.25*unit[ {0,-2,0} ]
];*)
poly2,
{i,n,1,-1}];

h1=Join[layers[[n]], {allbases[[1]]} ];(*pyramid*)
h12=Join[layers[[1]], {allbases[[n-1]]} ];(*pyramid*)
h=Table[(*antiprismoid chunks*)
Join[layers[[k]], {allbases[[n-k]],allbases[[n-k+1]]} ],
{k,2,Length[layers]-1}];

polychunks=Join[{h1},h,{h12}];(*assemble all chunks*)

Show[Graphics3D[ polychunks ],
PlotRange->All,
ViewPoint->2.25*unit[ {0,0,-1} ]
];

(*show chunks exploded apart slightly*)
chunkcents=Chop[Map[cent[vunion[ Flatten[ parse[ # ],1]]]&, polychunks]];
explodechunks=MapThread[trans[#1,.2*#2]&, {polychunks, chunkcents}];
Show[Graphics3D[ explodechunks ],
PlotRange->All,
ViewPoint->2.25*unit[ {0,-1,0} ]
];

(*volume of the circumscribing polar zonohedron*)
Print["n=",n]
Print["j=",j]
Print["k=",k]
edge=1;(*of unit edge length*)

pzvol = N[(edge^3)*(n^2)*(.5)*(Cot[Pi/n])*(Cos[Degree pitch]^2)*(Sin[Degree pitch])];
Print["Volume of circumscribing polar zonohedron is ", pzvol];

vols=Table[(*I am doubtful that the volumes are right*)
hullvol[ polychunks[[k]] ],
{k,Length[polychunks]}];

spivol=Plus@@vols;
Print["Volume of spirallohedron is ", spivol];

rv1=spivol/pzvol;
rv2=pzvol/spivol;
voldiff=pzvol-spivol;

Print["Ratio of spirallohedron to polar zonohedron is ", rv1];
Print["Difference between spirallohedron and polar zonohedron is ", voldiff];
Print["Volume of each of the k subtracted regions of hexahedra is ", subhex=voldiff/k];
Print["Ratio of a hexahedral region to polar zonohedron is ", subhex/pzvol];

Table of Initial Results

These were the results obtained by splitting spirallohedra into chunks, and then summing the volumes of the chunks. The third column is the volume ratio, spirallohedron/polar zonohedron; the fourth column gives the expression found, if any, which gave that ratio.

The solution is found!

Having failed to discover expressions for some cases, I finally set that job aside and considered the broad outlines. The ratios I sought had nothing to do with pitch; I had verified this by finding the volumes of various spirallohedra for various pitches. Finally it struck me that, just as the vertices of a k-armed spirallohedron fall into planes orthogonal to the symmetry axis, and the k vertices where these planes intersect the k ridges form the corners of a regular k-gon, so also the vertices of the circumscribing polar zonohedron lie on those same planes, and form the corners of regular n-gons. And, just as the k-armed spirallohedron can be inscribed within its related polar zonohedron, so also each k-gon is inscribed within its coplanar n-gon.

What, then, is the ratio of the area of a k-gon to that of its circumscribing n-gon? What is (Subscript[A, k])/(Subscript[A, n])?

A regular n-gon, of unit circumradius, divides into 2n equal-area right triangles of hypotenuse 1, base Sin[Pi/n] (half the side of the n-gon) and alternate side Cos[Pi/n] (in-radius of the n-gon). Two triangles equal a rectangle of sides Sin[Pi/n}, Cos[Pi/n]. Hence the area of the n-gon is n*Sin[Pi/n]*Cos[Pi/n].

n=12;
k=4;
area1=N[ n*Sin[Pi/n]*Cos[Pi/n] ];
area2=N[ k*Sin[Pi/k]*Cos[Pi/k] ];
Print["Ratio of areas, ",k,"-gon to ",n,"-gon, is ", area2/area1]
Clear[n,k]

The above example gives the ratio of the volume of an n=12, 4-armed spirallohedron to the volume of its circumscribing n=12 polar zonohedron, which was one of my problematic cases. Testing many other cases, including those for which I had found no expression, I found that I obtained the identical ratios. Hence, without proof, I think the ratio (Subscript[A, k])/(Subscript[A, n]) is the key to the volume of the general {n,k}-armed spirallohedron.

It only remained to combine this ratio with my expression giving the volume of the general polar zonohedron, that is, to multiply the volume of a polar zonohedron by

((k sin((π/k)) cos((π/k)))/(n sin((π/n)) cos((π/n))))

Now, we can substitute k(Sin[2*Pi/k]/2, and similarly for n, in the expression above.
Letting Mathematica combine terms, sort out identities, and simplify, we obtain an expression for the volume of the general k-armed spirallohedron..

Clear[n,j,k,edge,pitch]

(*volume polar zonohedron determined by n vectors,
edge length 'edge', complement to angle between edges and symmetry axis 'pitch'*)

polarvolume=edge^3 * n^2 * 1/2 * Cot[Pi/n] * Cos[pitch]^2 * Sin[pitch]

(*ratio, area regular k-gon to area circumscribing n-gon*)

polygonratio=((k*Sin[2*Pi/k])/2) / ((n*Sin[2*Pi/n])/2)

(*volume k-armed spirallohedron inscribed with polar zonohedron*)

spiralvolume=FullSimplify[ polygonratio*polarvolume ]

Hence, volume of {n,k}-armed Rhombic Spirallohedron:

\(1\/4\ edge\^3\ k\ n\ \(\(cos\^2\)(
              pitch)\)\ \(\(csc\^2\)(π\/n)\)\ \(sin(\(2\ π\)\/k)\)\ \(sin(pitch)\)\)

Below, we test this same expression in InputForm. Note that the example—a 3-armed, n=9 spirallohedron—is one of my problem cases, as seen in the Table of Results subsection, above.

n=9;edge=1;pitch=35.2643;k=3;

spirvol=(edge^3*k*n*Cos[Degree pitch]^2*Csc[Pi/n]^2*Sin[(2*Pi)/k]*Sin[Degree pitch])/4

(*Volume of Polar Zonohedron*)
pozvol=(edge^3)*(n^2)*(.5)*(Cot[Pi/n])*(Cos[Degree pitch]^2)*(Sin[Degree pitch])

spirvol/pozvol

Finally, it is apparent that the ratios of volumes (i.e., V[spirallohedron] / V[polar zonohedron]) for the various k—{k=3, 4, 5 ,...}—converge upon certain values. Consider the ratio of a 3-armed spirallohedron to its circumscribing polar zonohedron, as n grows large:

n=2400;edge=1;pitch=35.2643;k=3;

(*Volume of Spirallohedron*)
spirvol=(edge^3*k*n*Cos[Degree pitch]^2*Csc[Pi/n]^2*Sin[(2*Pi)/k]*Sin[Degree pitch])/4

(*Volume of Polar Zonohedron*)
pozvol=(edge^3)*(n^2)*(.5)*(Cot[Pi/n])*(Cos[Degree pitch]^2)*(Sin[Degree pitch])

Print["The ratio is ", spirvol/pozvol]

As n->infinity, this ratio is simply

((k sin((π/k)) cos((π/k)))/π)

Tabulating these limiting volume ratios, we have

Table[
N[ k*Sin[Pi/k]*Cos[Pi/k]/Pi ],
{k,3,13}]

Volume II: Hexahedral Composition of Spirallohedra

This was originally a separate notebook, with its own Initialization section. It has been rolled into the main Spirallohedron notebook intact. Included are the creation of k-regions from hexahedra, the creation of polar zonohedra from hexahedra, the creation of 3- and 4-armed spirallohedra from hexahedra, and the discovery and construction of nested, hidden “quasi-spirallohedra” within spirallohedra.

Discussion

Warning: my PolarZonohedron and Spirallohedron functions accept the variable pitch in degrees, and set pitch by default to 35.2643+ degrees. In many examples below, however, pitch is given in radian measure, for the purposes of certain embedded functions. It is always set to 35.2643+ degrees, in order to allow for construction of relevant spirallohedra, etc. If one wants to change pitch in this section, one must use considerable caution. My recommendation: leave pitch alone. It does not bear greatly upon the essentially combinatorial subjects explored herein.

Rhombic spirallohedra can be constructed by subtracting rhombic hexahedra from polar zonohedra. This notebook examines such constructions, and issues bearing upon the volumes of polar zonohedra, rhombic spirallohedra, and rhombic hexahedra.


Polar zonohedra are determined by a set of n vectors, and may be dissected into the “tetrahedral” number of n-things-taken-three-at-a-time rhombic hexahedra. In Mathematica we can express this tetrahedral number as follows:

Table[Print["When n=",n," there are ",Binomial[n,3]," hexahedra."],
{n,3,12}];

When n=3 or n=4, there is only one “type” of hexahedron, but more typically there are several to many types. Just as the polar zonohedron is determined by n vectors, so also each individual hexahedron is determined by some three of those n vectors. We can use the “KSubsets” function from Steven Skiena’s Combinatorica package to make a list of the triples of vector indices. For example, when n=5, we obtain this list of indices:

KSubsets[Range[5],3]

Polar zonohedra are bounded by n(n-1) faces. Any two of the n vectors determines a face, and, in reverse order, the equal and opposite counter-face. There are n-1 “rows” or cycles” of n faces, with row 1 equivalent to row n-1, row 2 equivalent to row n-2, and so on. When n is even, there is an “equatorial” row with no equivalent. That is, when n=6, we have five rows: row 1 (equiv.row 5), row 2 (equiv.row 4), and row 3 (equatorial row).

In the list of indices above, we can find the “face indices” of a hexahedron by subtracting pairs of indices in turn. For instance, the first triple is {1,2,3}. 2-1=1, 3-2=1, 3-1=2, hence, this hexahedron, in face-index form, can be represented as [1,1,2]. For some purposes it is useful to ignore equivalence between rows, in other cases one might want to reveal equivalence. The triple {1,4,5} provides an example. Ignoring equivalence, we have face-indices [3,1,4], but applying equivalence (and reducing to lowest terms) we have [2,1,1]. By sorting these indices we get [1,1,2] and can conclude that the hexahedra determined by triples {1,2,3} and {1,4,5} are the same.

If we actually test the equivalence of hexahedra {1,2,3} and {1,4,5} by finding their volumes, we find that they are at least equal in volume. However, the problem of deciding exactly how many different types of hexahedra do exist, for some given n, is a little difficult. One might take the approach that, if two hexahedra have the same volume, then they are of the same type. That, unfortunately, is not always the case. It is almost always the case, but not always. There exist “hexahedral identities” when a hexahedron of one type has the same volume as a hexahedron of another type. They are bounded by different faces, yet have the same volume. The lowest n for which this occurs is 12.

Such problems can be explored in more detail in the sections below. Eventually, we do dissect a polar zonohedron into its entire suite of rhombic hexahedra, not just mathematically, but in the sense of constructing each particular hexahedron, with its six rhombic faces, and displaying all or part of them. We also construct “k-regions” of rhombic hexahedra, being that exact set of hexahedra which must be subtracted from a polar zonohedron to create the k “valleys” bounding a k-armed spirallohedron. Finally, we construct 3- and 4-armed spirallohedra from hexahedra, and discover certain hidden quasi-spirallohedron, nested within.

Volumes of Hexahedra, a First Look

The n vectors which determine a polar zonohedron are made using my “vectors" function, which was designed to take the angle-between-edges-and-a-plane-orthogonal-to-the-symmetry-axis ‘pitch’ in degrees; so there is a little extra confusion below, where pitch is given in radian measure. We set n, pitch, edge length, make the n vectors, find the triples of indices, and then find the volumes of each hexahedron by taking the determinant of vector triple, and display results in the form {triple, volume}.

n=6;
pitch=ArcTan[ Sqrt[1/2] ];(*at this pitch, orthogonal shadow of hypercube*)
edge=1;(*set to (3/n)^.5 to obtain a "normalized eutactic star"*)
v=vectors[n,1/Degree pitch,edge];(*correct for radian measure*)

triples=KSubsets[Range[n],3];(*the idea here is that the hexahedra must be sorted by index not by volume*)
indexedhexvols=Table[{triples[[i]],
Det[ v[[ triples[[i]] ]] ]},{i,Length[triples]}]

We see that there are three different volumes scattered within the list. We reorder by volume.

dex1=Split[
	Sort[indexedhexvols, (#1[[2]]<#2[[2]]&)],
(Chop[ #1[[2]]-#2[[2]] ]==0&)]

Now we reorder each volume set from lowest to highest triple..

dex2=Map[Sort[# ]&, dex1]

Finally, we get the “representative” index triple for each volume. It so happens that when n<11, the ordering of volumes is the same as the ordering of representative vector index triples. Go back to the top and set n=11 to see an out-of-order set of volumes.

dex3=Sort[Map[First[#]&, dex2]](*indices in lowest terms, volumes for that class*)
Print["When n=",n," the volumes are ordered? ",OrderedQ[newdiffvols=Map[#[[2]]&, dex3]]]

Distinct Hexahedra and Their Volumes

When we use the approach in the section above, we do not pick up the rare hexahedral identities. The first instance of such an identity is when n=12. By converting the vector-index triples to face-index triples in lowest terms, so that when a face of row n-k appears, we treat it as a face of row k, we find that although there are 12 different types of hexahedra when n=12, there are only 11 different volumes. Experience has shown that the ratios between hexahedral volumes are unaffected by pitch, or by edge length, so to experiment further, change only n.

Binomial[n,3]

KSubsets[Range[n],3]

(*set n, pitch, edge length, make vectors*)
n=11;
pitch=ArcTan[ Sqrt[1/2] ];(*at this pitch, orthogonal shadow of hypercube*)
edge=1;(*set to (3/n)^.5 to obtain a "normalized eutactic star"*)
v=vectors[n,1/Degree pitch,edge];(*correct for radian measure*)

hemi=If[EvenQ[n],n/2,((n-1)/2)];(*face-number of equatorial cycle(s)*)
triples=KSubsets[Range[n],3];(*each triple determines a hexahedron*)

facehex=Map[Sort[ cycleface[{#[[2]]-#[[1]], #[[3]]-#[[2]], #[[3]]-#[[1]]}] ]&, triples];(*preserves order of triples while converting to face-indices*)
Print["The lowest-term face-index triples are ", dexface=Union[facehex] ];
Print[];
Print["There are ", numhex=Length[dexface], " different types of hexahedra"];

q1=Table[(*find positions in general list of each type of face-indexed hexahedron*)
Flatten[Position[ facehex, dexface[[j]]]],
{j,Length[dexface]}];

Print["The ",Length[triples]," hexahedra fall into sets of ",hexsets=Map[Length[#]&, q1]];(*this many of each type, ordered from lowest face-index to highest*)
Print[];Print[];

rephex=Map[ First[#]&, q1];(*just choose the first of each sublist to represent that type*)
(*go ahead and get representative volumes*)
rephexvols=Map[ Det[ v[[ # ]] ]&, triples[[rephex]] ];

(*now go back to vector indices for each type*)
Print["The vector triples are ", triples[[rephex]]];
Print[];Print[];
(*recall the corresponding face-index triples*)
Print["The face triples are ", dexface];
Print[];Print[];

(*and find the volumes for each type, in face-index order*)
Print["The corresponding volumes are ", rephexvols];

Print[];
Print["There are ",
Length[ Union[rephexvols,SameTest->(Chop[ #1-#2 ]==0&)] ],
" different volumes."]

From the above we see that although there are 12 different types of hexahedra, there are only eleven different volumes. Volume 4 equals volume 6, so that (by face-index) hexahedron [1,4,5] has the same volume as hexahedron [2,2,4].

Another way of computing the volume of a hexahedron is, given the three vectors which determine it, v1, v2, v3, say, to take Abs[ v1 . Cross[ v2, v3] ].

Det[ v[[ {1,2,3} ]] ]==Abs[ v[[1]]. Cross[ v[[2]], v[[3]] ] ]

Clear[nn,i,j,k,v1,v2,v3,h]
v1={Cos[i*2*Pi/nn],Sin[i*2*Pi/nn],h};
v2={Cos[j*2*Pi/nn],Sin[j*2*Pi/nn],h};
v3={Cos[k*2*Pi/nn],Sin[k*2*Pi/nn],h};

Simplify[Abs[ v1.Cross[v2,v3] ]]

When we use Mathematica's symbolic capabilities on the expression Abs[ v1 . Cross[ v2, v3] ] we obtain an interesting result. The details are a little involved, but in essence, the n vectors as I define them are of the form {Cos[i*2*Pi/n], Sin[i*2*Pi/n], h}, and when these are used symbolically, we find that the “core symbolic terms” (neglecting edge length and pitch), for a hexahedron of face-indices [i,j,k], are ( 4* Sin[i*Pi/n] * Sin[j*Pi/n] * Sin[k*Pi/n] ). See above. For our purposes let’s forget the ‘4’ for now.

(*now we do not reduce face-indices to lowest terms*)
faceindices=Map[
{#[[2]]-#[[1]], #[[3]]-#[[2]], #[[3]]-#[[1]]}&, triples[[rephex]] ];
Print["The face indices, not reduced to lowest terms, are ", faceindices];
Print[];Print[];
Print["The core symbolic terms for hexahedral volumes are ", Map[sinproduct[#]&, faceindices]];

When we give Mathematica an expression such as Cos[Pi/6], however, it will if possible return some other equivalent. See below.

Cos[Pi/6]

Hence, to make a table of vector index triples, face index triples, core symbolic expressions, and actual numeric volumes, a special “sinproduct” function was cobbed together, to force the desired expressions to appear. Note that the face-index triples are not in lowest equivalent terms.

symbs=Map[{
#,
{#[[2]]-#[[1]], #[[3]]-#[[2]], #[[3]]-#[[1]]},
sinproduct[ {#[[2]]-#[[1]], #[[3]]-#[[2]], #[[3]]-#[[1]]} ],
N[symhexvol[#, n, pitch ]]
}&, triples[[rephex]] ];


TraditionalForm[
TableForm[symbs,
TableSpacing->{4,2},
TableHeadings->{None,
{"Vectors", "Faces", "Core Symbolic Expression","Numeric Volume"}
}]
]

Now, recalling that we have obtained the number of hexahedra for each different type, we can approach the total volume of a polar zonohedron as a sum of multiples of the volumes for each type. Note that here no attempt was made to force Mathematica to use the “core symbolic expressions.”

MapThread[{
#1, symhexvol[#2, n, pitch ]
}&, {hexsets,triples[[rephex]]} ]

Here we actually form each particular product, and display in TraditionalForm.

TraditionalForm[TableForm[
MapThread[{
(#1) * (symhexvol[#2, n, pitch ])
}&, {hexsets,triples[[rephex]]} ]
]]

We add the products to obtain the volume of the polar zonohedron, and then Simplify the resulting expression.

pv=Plus@@MapThread[{
(#1) * (symhexvol[#2, n, pitch ])
}&, {hexsets,triples[[rephex]]} ]
TraditionalForm[Simplify[pv]]

Here we try a different form, and then Simplify the resulting expression.

qq=TrigToExp[Plus@@MapThread[{
(#1) * (symhexvol[#2, n, pitch ])
}&, {hexsets,triples[[rephex]]} ]
]
TraditionalForm[Simplify[qq]]

We have obtained the volume of the polar zonohedron.

N[pv]

Volumes of Polar Zonohedra, Spirallohedra, and k-regions

We obtain symbolic representations of various expressions and use them to construct functions to find the volumes of polar zonohedra, rhombic spirallohedra, and k-regions.

Years of effort led gradually to the discovery of the following expression for the volume of the general polar zonohedron, where n is the number of vectors, pitch is the complement to the angle between edges and symmetry axis, and edge is the edge length.

Clear[n,j,k,edge,pitch,pzvol, spivol]

(*volume polar zonohedron determined by n vectors,
edge length 'edge', complement to angle between edges and symmetry axis 'pitch'*)

polarvolume=edge^3 * n^2 * 1/2 * Cot[Pi/n] * Cos[pitch]^2 * Sin[pitch]

(*We cast the expression as a function accpeting n,pitch and edge as inputs*)
(*note: pitch is assumed to be given in degrees*)

pzvol[n_,pitch_,edge_]:=(edge^3*n^2*Cos[Degree*pitch]^2*Cot[Pi/n]*Sin[Degree*pitch])/2

See the Volume section above for details about how the volume of the general k-armed spirallohedron was found. It involves the ratio between a regular k-gon and a (circumscribing) regular n-gon. A regular polygon of j sides can be divided into j equal isosceles triangles. Each such triangle is half of a rhomb with two interior angles of 360/j degrees. The area of such a rhomb is Sin[360/j]; hence the area of the j-gon is j*Sin[360/j]/2.

(*ratio, area regular k-gon to area circumscribing n-gon*)

polygonratio=((k*Sin[2*Pi/k])/2) / ((n*Sin[2*Pi/n])/2)

The product of the polygon ratio and the volume of the circumscribing polar zonohedron gives the volume of the inscribed k-armed spirallohedron.

(*volume k-armed spirallohedron inscribed with polar zonohedron*)

spiralvolume=FullSimplify[ polygonratio*polarvolume ]

(*volume of k-armed spirallohedron*)
(*note: pitch is assumed to be given in degrees*)

spivol[n_,k_,pitch_,edge_]:=(edge^3*k*n*Cos[Degree*pitch]^2*Csc[Pi/n]^2*Sin[(2*Pi)/k]*Sin[Degree*pitch])/4

Since k-armed spirallohedra can be constructed by subtracting k regions of hexahedra from the circumscribing polar zonohedron, we find an expression for the volume of one of the k regions.

(*volume of one of the k subtracted regions of hexahedra*)

kvol=FullSimplify[ (polarvolume-spiralvolume) / k ]

(*volume of one of the subtracted k regions of hexahedra*)
(*note: pitch is assumed to be given in degrees*)

kvolume[n_,k_,pitch_,edge_]:=(edge^3*n*Cos[Degree*pitch]^2*Csc[Pi/n]^2*(-(k*Sin[(2*Pi)/k]) + n*Sin[(2*Pi)/n])*Sin[Degree*pitch])/(4*k)

Here are some examples.

pzvol[12, 35.2643, 1]
spivol[12, 3, 35.2643, 1]
kvolume[12, 3, 35.2643, 1]

Make Hexahedral k-Regions

Here we construct the hexahedra forming one of the k regions subtracted from a polar zonohedron to make a spirallohedron.

Once again we adopt n=12. A hexahedron function was written expressly for this construction. The function accepts three vectors, vi, vj, and vk, and a reference point-of-beginning (the “top” of the hexahedron) ref. It creates the six faces and updates ref, defining the top of the next hexahedron in turn. In a close-packing of polyhedra every face is the join between two polyhedra; in this local, finite close-packing of hexahedra, which are zonohedra, every face is part of a “family” or set of equal and parallel “interfaces.” The first face created by the hexahedron function is propagated through a series of interfaces in nested loops. That is to say, consider the vector index triples, which are “canonically” ordered from least to greatest. The list of integer triples typically begins with a sequence such as {{1,2,3}, {1,2,4}, {1,2,5}, ...}. In the construction below, the common members of such a triple sequence, in this case {1,2}, form the continuous series of equal interfaces.

(*set n, pitch, edge length, make vectors, find triples*)
n=12;
pitch=ArcTan[ Sqrt[1/2] ];(*at this pitch, orthogonal shadow of hypercube*)
edge=1;(*set to (3/n)^.5 to obtain a "normalized eutactic star"*)
v=vectors[n,1/Degree pitch,edge];(*correct for radian measure*)
triples=KSubsets[Range[n],3];(*each triple determines a hexahedron*)

(*try out hexahedron function*)
Show[newhex[ v[[{1,2,3}]], {0,0,0}],
Axes->True,
ViewPoint->{0,-1,3}];

Now we create the vector indices for one k-region, in this case, k=3, but we call it karm. Nested loops are used, the innermost defining integer triples {i,j,k}. The nested output is flattened, the number of triples found, and the resulting region volume is first calculated as a sum of hexahedral volumes, then using the kvolume function. The two results are tested for equality.

karm=3;(*set number of arms of spirallohedron, karm a divisor of n*)
q=Table[(*set up list of vector triples for region hexahedra*)
Table[
	Table[{i,j,k},
	{k,j+1,n}],
{j,i+1,(n/karm)}],
{i,1,(n/karm)-1}]

regiontriples=Flatten[q,2];
Print["There are ", Length[regiontriples], " hexahedra in the region."];
Print["The volume of the hexahedra is ", vol1=Plus@@Map[ Det[ v[[ # ]] ]&, regiontriples ]];
Print["The volume of the region is ", vol2=N[kvolume[n,karm,pitch,edge]]];(*volume of one of k regions*)
Print["Computed volume equal to measured volume? ", vol1==vol2];

Now we actually make each hexahedron. In oder to force the k-region to align with a k-armed spirallohedron centered on the origin, the negative of the point “pt” is defined as the point-of-beginning. However, this will not suffice for all the hexahedra needed, and observation and experimentation led to the somewhat awkward construction for other points-of-beginning, -pt+(Plus@@v[[ Range[i-1] ]]). 

We also go ahead and make the associated k-armed spirallohedron, and display, first, the region itself, and then next, the region and the spirallohedron, with the region translated slightly to separate it.

pt=Chop[ (Plus@@v)/2 ];(*create region hexahedra*)
zaq=Table[
Table[If[i==1,mark=-pt,mark=-pt+(Plus@@v[[ Range[i-1] ]]) ];
	Table[ newhex[ v[[{i,j,k}]],mark],
	{k,j+1,n}],
{j,i+1,n/karm}],
{i,1,(n/karm)-1}];

a=Spirallohedron[n,karm];(*create k-armed spirallohedron*)


Show[zaq,(*show region hexahedra*)
Axes->True,
ViewPoint->{0,0,-3}];

Show[{a,trans[zaq, {-1,1,0}]},(*show spirallohedron and translated region*)
Axes->True,
ViewPoint->2.25*unit[{0,-.5, -2}]];

Here the hexahedra of the k-region are exploded apart slightly and displayed, alone and with the associated spirallohedron. The explosion is performed by first flattening the nested lists of hexahedra into one long list. A scale term sca is set, the center of each hexahedron is found, and that hexahedron is translated by the product of sca and that center.

Length[h1=Flatten[zaq,2]];
sca=0.25;
xaq=Map[ trans[#, sca*cent[ vunion[ Flatten[ parse[ # ], 1] ] ] ]&, h1];
Show[xaq,(*show hexahedra*)
Axes->True,
ViewPoint->2.25*unit[{0,0,3}]];

Show[{a,xaq},(*show hexahedra with spirallohedron*)
Axes->True,
ViewPoint->2.25*unit[{0,0,3}]];

Analyze Triples

We have constructed one k-region, in the example above, k=3. Here it is assumed that the section Make Hexahedral k-regions has been evaluated with k=3, and we form the vector indices for the other two regions by addition. After adding we sort so that the resulting triples go from lowest to highest terms. Then we find all those triples which are in the global list triples, comprising the entire polar zonohedron, but not in any of the regions. These are those triples which correspond to the hexahedra which compose the core spirallohedron, and we name this sublist core.

h2=(n/3)-1;
h3=(2*n/3)-1;
s1=regiontriples;
s2=Sort/@Map[ 1+Mod[{h2,h2,h2}+#, n]&, regiontriples];
s3=Sort/@Map[ 1+Mod[{h3,h3,h3}+#, n]&, regiontriples];
Print["The total number of hexahedra in all regions is ", r1=Length[s4=Union[s1,s2,s3]]];
Print[];
Print["The indices for hexahedra in the core spirallohedron are ", core=Complement[ triples, s4 ]];
Print[];
Print["The number of hexahedra in the core spirallohedron is ",r2=Length[core]];

Now we look at the average volumes of region hexahedra and core hexahedra.

Print["The volume of the regional hexahedra is ", vr=3*Plus@@Map[ Det[ v[[ # ]] ]&, regiontriples ]];

Print["Average volume of a regional hexahedron is ", vr/r1];

Print[];

Print["The volume of the core hexahedra is ", vc=Plus@@Map[ Det[ v[[ # ]] ]&, core ]];

Print["Average volume of a core hexahedron is ", vc/r2];

Print[];
Print["The ratio between the two volumes is ", vc/vr];
Print["The ratio between core volume and polar zonohedron is ", vc/(vr+vc)];

Make a Hexahedral Polar Zonohedron

In an unintended consequence of making a hexahedral k-region, by altering the lengths of the nested loops we can construct a polar zonohedron dissected into its tetrahedral number of hexahedra!

(*set n, pitch, edge length, make vectors, find triples*)
n=12;
pitch=ArcTan[ Sqrt[1/2] ];(*at this pitch, orthogonal shadow of hypercube*)
edge=1;(*set to (3/n)^.5 to obtain a "normalized eutactic star"*)
v=vectors[n,1/Degree pitch,edge];(*correct for radian measure*)
triples=KSubsets[Range[n],3];(*each triple determines a hexahedron*)

pt=Chop[ (Plus@@v)/2 ];(*create region hexahedra*)
zaq=Table[
Table[If[i==1,mark=-pt,mark=-pt+(Plus@@v[[ Range[i-1] ]]) ];
	Table[ newhex[ v[[{i,j,k}]],mark],
	{k,j+1,n}],
{j,i+1,n-1}],
{i,1,n-2}];

Print["There are by theory ", Binomial[n,3], " hexahedra, and we have constructed ", Length[Flatten[zaq,2]], " hexahedra."];

Show[zaq,(*show hexahedra*)
Axes->True,
ViewPoint->2.25*unit[{0,0,1}]];

Since the list “zaq” of hexahedra is nested, we can display each principal sublist.

Table[Show[zaq[[k]],(*a sublist of hexahedra*)
Axes->True,
ViewPoint->2.25*unit[{0,0,3}]],
{k,Length[zaq]}];

We can display each sublist within any principal sublist.

Table[Show[zaq[[1,k]],(*a sublist of hexahedra*)
Axes->True,
ViewPoint->2.25*unit[{0,0,3}]],
{k,Length[ zaq[[1]] ]}];

And we can explode the hexahedra apart.

h1=Flatten[zaq,2];
sca=0.5;
xaq=Map[ trans[#, sca*cent[ vunion[ Flatten[ parse[ # ], 1] ] ] ]&, h1];
Show[xaq,(*show hexahedra*)
Axes->True,
ViewPoint->2.25*unit[{0,0,3}]];

All Regions; Hexahedral Spirallohedra; Rinds and Layers

By further modifying the nested loops giving a “representative” k-region, we can construct all k of the k-regions. Then, subtracting the hexahedra from all k-regions from the hexahedra of the circumscribing polar zonohedron, we obtain the index triples for the “core” spiirallohedron. Analysis of these indices yielded cumbersome nested loops which give the hexahedral composition of 3- and 4-armed spirallohedra; no such construction is yet found for the “higher” k-armed spirallohedra. We sort hexahedra by volume and assign different colors to them. Finally, we discover hidden “quasi-spirallohedra” nested within the hexahedrally composed spirallohedra.

Further notes will be found in the particular subsections.

Construction of k-regions

By example we set n=15, k=3. We call k “karm” to avoid confusion with the nested loops, {i,j,k}.

(*set n, pitch, edge length, make vectors, find triples*)
n=15;
pitch=ArcTan[ Sqrt[1/2] ];(*at this pitch, orthogonal shadow of hypercube*)
edge=1;(*set to (3/n)^.5 to obtain a "normalized eutactic star"*)
v=vectors[n,1/Degree pitch,edge];(*correct for radian measure*)
triples=KSubsets[Range[n],3];(*each triple determines a hexahedron*)
karm=3;
kstep=n/karm;

We construct all k of the k-regions, and display them {1, 2, ..., k} at a time.

pt=Chop[ (Plus@@v)/2 ];(*sum of vectors, divided by 2*)

allregions=Table[

Table[
  Table[If[i==1,mark=-pt,mark=-pt+(Plus@@v[[ Range[sk+1,i-1] ]]) ];
	Table[ newhex[ v[[ 1+Mod[{i-1,j-1,k-1}, n] ]],mark],
	{k,j+1,n+sk}],
  {j, i+1, sk+(n/karm) }],
{i, sk+1, sk+((n/karm)-1) }],

{sk,0,n-kstep,kstep}];

Table[
Show[allregions[[ Range[i] ]],(*show two regions*)
Axes->True,
ViewPoint->{0,0,-3}],
{i, Length[allregions]}];

We construct the vector indices alone for all regions, eliminate these from the master list of indices for the circumscribing polar zonohedron, and thus obtain the index triples for the “core” k-armed spirallohedron. If you wish to see the list of core triples, remove the semi-colon from the end of the line core=Complement[triples, Sort/@regflat ];.

regdex=Table[

Table[
  Table[
	Table[ 1+Mod[{i-1,j-1,k-1}, n],
	{k,j+1,n+sk}],
  {j, i+1, sk+(n/karm) }],
{i, sk+1, sk+((n/karm)-1) }],

{sk,0,n-kstep,kstep}];

regflat=Flatten[Map[ Flatten[#,2]&, regdex],1];
core=Complement[triples, Sort/@regflat ];
Print["There are ", c=Length[core], " hexahedra in the core spirallohedron."];

(*we can compute the number of hexahedra*)
hexacore[n,karm]

Analyze Core Triples

We have created the list core in the section above, and also defined n, etc. Here we merely apply some of the same types of analysis as in the section “Distinct Hexahedra” above. There is nothing done here which is needed below, in subsequent sections.

hemi=If[EvenQ[n],n/2,((n-1)/2)];(*face-number of equatorial cycle(s)*)

spifacehex=Map[Sort[ cycleface[{#[[2]]-#[[1]], #[[3]]-#[[2]], #[[3]]-#[[1]]}] ]&, core];(*preserves order of triples while converting to face-indices*)
Print["The lowest-term face-index triples are ", spidexface=Union[spifacehex] ];
Print[];
Print["There are ", spinumhex=Length[spidexface], " different types of hexahedra"];

spiq1=Table[(*find positions in general list of each type of face-indexed hexahedron*)
Flatten[Position[ spifacehex, spidexface[[j]]]],
{j,Length[spidexface]}];

Print["The ",Length[core]," hexahedra fall into sets of ",spihexsets=Map[Length[#]&, spiq1]];(*this many of each type, ordered from lowest face-index to highest*)
Print[];Print[];

spirephex=Map[ First[#]&, spiq1];(*just choose the first of each sublist to represent that type*)
(*go ahead and get representative volumes*)
spirephexvols=Map[ Det[ v[[ # ]] ]&, core[[spirephex]] ];

(*now go back to vector indices for each type*)
Print["The vector triples are ", core[[spirephex]]];
Print[];Print[];
(*recall the corresponding face-index triples*)
Print["The face triples are ", spidexface];
Print[];Print[];

(*and find the volumes for each type, in face-index order*)
Print["The corresponding volumes are ", spirephexvols];

Print[];
Print["There are ",
Length[ Union[spirephexvols,SameTest->(Chop[ #1-#2 ]==0&)] ],
" different volumes."]

(*now we do not reduce face-indices to lowest terms*)
spifaceindices=Map[
{#[[2]]-#[[1]], #[[3]]-#[[2]], #[[3]]-#[[1]]}&, core[[spirephex]] ];
Print["The face indices, not reduced to lowest terms, are ", spifaceindices];
Print[];Print[];
Print["The core symbolic terms for hexahedral volumes are ", Map[sinproduct[#]&, spifaceindices]];

spisymbs=Map[{
#,
{#[[2]]-#[[1]], #[[3]]-#[[2]], #[[3]]-#[[1]]},
sinproduct[ {#[[2]]-#[[1]], #[[3]]-#[[2]], #[[3]]-#[[1]]} ],
N[symhexvol[#, n, pitch ]]
}&, core[[spirephex]] ];


TraditionalForm[
TableForm[spisymbs,
TableSpacing->{4,2},
TableHeadings->{None,
{"Vectors", "Faces", "Core Symbolic Expression","Numeric Volume"}
}]
]

Create Spirallohedron

So far, the nested loop structure below only works for 3- and 4-armed spirallohedra.

(*set n, pitch, edge length, make vectors, find triples*)
n=12;
pitch=ArcTan[ Sqrt[1/2] ];(*at this pitch, orthogonal shadow of hypercube*)
edge=1;
v=vectors[n,1/Degree pitch,edge];(*correct for radian measure*)
karm=3;
kstep=n/karm;

pt=Chop[ (Plus@@v)/2 ];
(*works for all 3- and 4-armed*)
spirhex=Table[
  Table[
  If[i<2,mark=-pt,mark=-pt+(Plus@@v[[ Range[1,i-1] ]]) ];
  If[j>(kstep)+1 && j<(2*kstep)+1,
    mark=mark+(Plus@@v[[ Range[(kstep)+1,j-1] ]]),
	mark=mark+(Plus@@v[[ Range[(2*kstep)+1,j-1] ]]) ];	
	Table[ newhex[ v[[ 1+Mod[{i-1,j-1,k-1}, n] ]], mark],
	{k, If[j<(2*kstep)+1, (2*kstep)+1, (3*kstep)+1], n}],
  {j, If[i<kstep+1, kstep+1, (2*kstep)+1], (n-kstep) }],
{i, 1, n-(2*kstep) }];

Show[spirhex,(*show spirallohedron composed by hexahedra*)
Axes->True,
ViewPoint->{0,0,3}];

We can, of course, explode the hexahedra apart.

h1=Flatten[spirhex,2];
sca=0.25;
xaq=Map[ trans[#, sca*cent[ vunion[ Flatten[ parse[ # ], 1] ] ] ]&, h1];
Show[xaq,(*show hexahedra*)
Axes->True,
ViewPoint->2.25*unit[{0,0,3}]];

Color Hexahedra by Volume

Having succeeded in creating 3- and 4-armed spirallohedra from rhombic hexahedra, I wondered what the distribution of volume—of the different types of hexahedra, roughly—might look like. I modified the nested loop to give the output {indices, hexahedron}, to facilitate finding the volume of any given hexahedron.

We use n=24, k=3, to get a reasonably complicated spirallohedron for a test subject.

(*set n, pitch, edge length, make vectors, find triples*)
n=24;
pitch=ArcTan[ Sqrt[1/2] ];(*at this pitch, orthogonal shadow of hypercube*)
edge=1;(*set to (3/n)^.5 to obtain a "normalized eutactic star"*)
v=vectors[n,1/Degree pitch,edge];(*correct for radian measure*)
karm=3;
kstep=n/karm;

(*creat list {triple, hexahedron}*)
pt=Chop[ (Plus@@v)/2 ];
(*works for all 3- and 4-armed*)
dexspirhex=Table[
  Table[
  If[i<2,mark=-pt,mark=-pt+(Plus@@v[[ Range[1,i-1] ]]) ];
  If[j>(kstep)+1 && j<(2*kstep)+1,
    mark=mark+(Plus@@v[[ Range[(kstep)+1,j-1] ]]),
	mark=mark+(Plus@@v[[ Range[(2*kstep)+1,j-1] ]]) ];	
	Table[ {1+Mod[{i-1,j-1,k-1}, n], newhex[ v[[ 1+Mod[{i-1,j-1,k-1}, n] ]], mark] },
	{k, If[j<(2*kstep)+1, (2*kstep)+1, (3*kstep)+1], n}],
  {j, If[i<kstep+1, kstep+1, (2*kstep)+1], (n-kstep) }],
{i, 1, n-(2*kstep) }];

Here we look at the hexahedra, exploded slightly, to see if they actually exist.

h1=Flatten[dexspirhex,2];
sca=0.25;
xaq=Map[ trans[#[[2]], sca*cent[ vunion[ Flatten[ parse[ #[[2]] ], 1] ] ] ]&, h1];
Show[xaq,(*show hexahedra*)
Axes->True,
ViewPoint->2.25*unit[{0,0,3}]];

Now we apply the volumeSort function to split the hexahedra into sublists according to volume.

h2=Flatten[dexspirhex,2];
Print["There are ",Length[h2], " hexahedra."];
Print[];
shex=volumeSort[ h2 ];
Print[];
Print["The hexahedra fall into sets of ", Map[Length[#]&, shex], " according to volume."];
Print[];
Print["There are ", colnum=Length[shex], " different volumes."];
colorlist=Table[GrayLevel[k/colnum],{k, Length[shex]}];

Make a list of shades of gray, and display.

colorlist=Table[GrayLevel[k/colnum],{k, Length[shex]}];

qaq=Graphics3D[
Table[
	Table[ {SurfaceColor[colorlist[[i]],GrayLevel[0], 2 ],
		First@shex[[i,k]]},{k, Length[ shex[[i]] ]}],
		{i,Length[shex]}]
];

Show[qaq,
Boxed->False,
Background->GrayLevel[0],
LightSources -> 
   {
   {{-1,1,0}, RGBColor[1, 1, 1]},
   {{1,0,1}, RGBColor[1, 1, 1]}
   },
ViewPoint->2.5*unit[{0, 0, 1}]];

Display the colored-by-volume hexahedra exploded slightly.

h1=Flatten[qaq[[1]],1];
Length[h1]
sca=0.25;
xaq=Graphics3D[{EdgeForm[],
Map[ trans[#, sca*cent[ vunion[ Flatten[ parse[ # ], 1] ] ] ]&, h1]}
];
Show[xaq,(*show hexahedra*)
Boxed->False,
Background->RGBColor[0, 0, .5],
LightSources -> 
   {
   {{-1,-1,0}, GrayLevel[1]},
   {{1,0,1}, GrayLevel[1]}
   },
ViewPoint->2.25*unit[{0,0,1}]];

Rinds, Layers, and Hidden Quasi-spirallohedra

I have been thinking about the arrangement of rhombic hexahedra within polar zonohedra for a long time. It is quite difficult to visualize. When n is greater than, say, 8, there are many different arrangements possible. One can perform other types of dissection upon polar zonohedra, too, for instance, into a mixture of rhombic dodecahedra and rhombic hexahedra, or into many different types of rhombic zonohedra.

As mentioned above, the ratios between hexahedra volumes are fixed, so let us imagine now that pitch is 35.2643+ degrees, so that any given type of hexahedron is at its maximum volume.

Whenever n is a multiple of 3, there are n/3 cubes within the mix of hexahedra. The positions of these cubes within the polar zonohedron are not tightly constrained. One is free, for instance, to locate the cubes along the symmetry axis. For that matter, whenever n is a multiple of 4, there are n/4 of Kepler’s Rhombic Dodecahedra within a polar zonohedron, each formed from four equal hexahedra. We are free to locate these dodecahedra along the symmetry axis of the polar zonohedron, and let all other hexahedra accommodate their positions to these dodecahedra.

3-armed spirallohedra clearly have cubes embedded at their north and south poles; 4-armed spirallohedra have Kepler’s Rhombic Dodecahedra embedded at their north and south poles; and the general k-armed spirallohedron have n=k polar zonohedra embedded at their north and south poles. Are the positions of these zonohedra, within spirallohedra, more tightly constrained than they are within the circumscribing polar zonohedra? This seems likely.

Here we discover that, whatever the total range of positions of hexahedra or higher rhombic zonohedra may be, we can find “quasi-spirallohedra” nested within 3- and 4-armed spirallohedra. For reasons which are not yet clear to me, the hidden quasi-spirallohedra are more symmetrical in 3-armed spirallohedra. When n is large, n an odd multiple of 3, we develop a series of nested 3-armed quasi-spirallohedra, with the two end members of the series being the outermost, “true” spirallohedron, while the innermost 3-armed quasi-spirallohedron is simply a single cube. Were one to animate this series in reverse one would see a cube morphing into a 3-armed spirallohedron.

However, the equivalent series, with n a multiple of 4, would not show Kepler’s Rhombic Dodecahedron morphing into a 4-armed spirallohedron. This will be apparent if one experiments.

The initial example is n=15, k=3. Try these other pairs of {n, k}: {16, 4}, {12, 3}, {12, 4}.

(*set n, pitch, edge length, make vectors, find triples*)
n=15;
pitch=ArcTan[ Sqrt[1/2] ];(*at this pitch, orthogonal shadow of hypercube*)
edge=1;(*set to (3/n)^.5 to obtain a "normalized eutactic star"*)
v=vectors[n,1/Degree pitch,edge];(*correct for radian measure*)
karm=3;(*number of arms of spirallohedron*)
kstep=n/karm;(*length of step from one arm to next*)

To facilitate selection of hexahedral layers, we form the expression {index triple, center of hexahedron, hexahedron}.

(*create list {triple, hexahedron center, hexahedron}*)
pt=Chop[ (Plus@@v)/2 ];
(*works for all 3- and 4-armed*)
dexspirhex=Table[
  Table[
  If[i<2,mark=-pt,mark=-pt+(Plus@@v[[ Range[1,i-1] ]]) ];
  If[j>(kstep)+1 && j<(2*kstep)+1,
    mark=mark+(Plus@@v[[ Range[(kstep)+1,j-1] ]]),
	mark=mark+(Plus@@v[[ Range[(2*kstep)+1,j-1] ]]) ];
		
	Table[
	{
	1+Mod[{i-1,j-1,k-1}, n],(*index triple*)
	mark+(Plus@@v[[ 1+Mod[{i-1,j-1,k-1}, n] ]]/2),(*hexahedron center*)
	newhex[ v[[ 1+Mod[{i-1,j-1,k-1}, n] ]], mark](*hexahedron*)
	},
	
	{k, If[j<(2*kstep)+1, (2*kstep)+1, (3*kstep)+1], n}],
  {j, If[i<kstep+1, kstep+1, (2*kstep)+1], (n-kstep) }],
{i, 1, n-(2*kstep) }];

(*flatten nested list of {triples, centers, hexahedra}*)
allhex=Flatten[dexspirhex,2];

We are going to find the outermost layer of hexahedra by selecting all hexahedra with centers closest to face-centers of the spirallohedron.

(*find face-centers of spirallohedron*)
a=Spirallohedron[n,karm];
fcents=Chop[Map[ cent[ vunion[#] ]&,  parse[a]]];

(*find those hexahedra closest to face centers; outermost layer*)
selhex=Table[
Select[allhex, mag[#[[2]]-fcents[[k]] ]<0.6&],
{k,Length[fcents]}];

Now we discard duplicate hexahedra, and then subtract the layer from the spirallohedron.

(*discard duplicate hexahedra and flatten list*)
Length[ layer1=Union[Flatten[selhex,1], SameTest->(#1[[1]]==#2[[1]]&)] ]

(*now find all hexahedra not in outermost layer*)
Length[ innerhex=Complement[ allhex, layer1, SameTest->(#1[[1]]==#2[[1]]&) ] ]

(*hexahedra of hidden quasi-spirallohedron*)
hidden=Map[#[[3]]&, innerhex];

Show[hidden,(*show hexahedra of quasi-spirallohedron*)
Axes->True,
ViewPoint->2.25*unit[{0,0,3}]];

Finally we discard all coincident faces, giving a quasi-spirallohedron as a list of its exterior faces..

Print["There are ",Length[f1=parse[hidden]]," polygons altogether."]
Print["There are ", Length[f2=dePair[f1] ]," unshared polygons."]
newsurface=Graphics3D[Polygon /@ f2];

Show[newsurface,(*show faces of hidden quasi-spirallohedron*)
Axes->True,
ViewPoint->2.25*unit[{0,0,3}]];

Now we can find the next hidden quasi-spirallohedron. For 3-armed spirallohedra, if n=15, as when n is odd, generally, the innermost quasi-spirallohedron is a single cube. It is as though n is descending in steps of six: 15, 9, 3. When n is even, the innermost quasi-spirallohedron is a quasi n=6. Below, we simply overwrite pre-existing variables so that the process of hidden quasi-spirallohedron extraction can be iterated; each time we finish by finding the exterior faces, and name this “newsurface.” Similarly, each time we name the hidden hexahedra “innerhex.”

(*find face-centers of spirallohedron*)
a=newsurface;
fcents=Chop[Map[ cent[ vunion[#] ]&,  parse[a]]];

(*find those hexahedra closest to face centers; outermost layer*)
selhex=Table[
Select[innerhex, mag[#[[2]]-fcents[[k]] ]<0.6&],
{k,Length[fcents]}];

(*discard duplicate hexahedra and flatten list*)
Length[ layer1=Union[Flatten[selhex,1], SameTest->(#1[[1]]==#2[[1]]&)] ]

(*now find all hexahedra not in outermost layer*)
Length[ innerhex=Complement[ innerhex, layer1, SameTest->(#1[[1]]==#2[[1]]&) ] ]

(*hexahedra of hidden quasi-spirallohedron*)
hidden=Map[#[[3]]&, innerhex];

Show[hidden,(*show hexahedra of quasi-spirallohedron*)
Axes->True,
ViewPoint->2.25*unit[{0,0,3}]];

Print["There are ",Length[f1=parse[hidden]]," polygons altogether."]
Print["There are ", Length[f2=dePair[f1] ]," unshared polygons."]
newsurface=Graphics3D[Polygon /@ f2];
Show[newsurface,(*show faces of hidden quasi-spirallohedron*)
Axes->True,
ViewPoint->2.25*unit[{0,0,3}]];

Now we will fully iterate the process, and create an animation. We must have already set n and karm and made the spirallohedron of hexahedra first. The discovery of the hidden quasi-spirallohedra, below, becomes a somewhat lengthy process, when n>15.

newsurface=Spirallohedron[n,karm];
innerhex=allhex;
jlim=Floor[(n-3)/6];(*number of steps*)

quasi=Join[{newsurface},
Table[

(*find face-centers of spirallohedron*)
fcents=Chop[Map[ cent[ vunion[#] ]&,  parse[newsurface]]];

(*find those hexahedra closest to face centers; outermost layer*)
selhex=Table[
Select[innerhex, mag[#[[2]]-fcents[[k]] ]<0.6&],
{k,Length[fcents]}];

(*discard duplicate hexahedra and flatten list*)
layer1=Union[Flatten[selhex,1], SameTest->(#1[[1]]==#2[[1]]&)];

(*now find all hexahedra not in outermost layer*)
innerhex=Complement[ innerhex, layer1, SameTest->(#1[[1]]==#2[[1]]&) ];

(*hexahedra of hidden quasi-spirallohedron*)
hidden=Map[#[[3]]&, innerhex];
f1=parse[hidden];
f2=dePair[f1];
newsurface=Graphics3D[Polygon /@ f2],
{jlim}]
];

Table[(*display series of polyhedra*)
Show[ quasi[[k]],
Boxed->False,
ViewPoint->{0,0,2.5}],
{k, Length[quasi]}];

Experimental

Here we make k-regions with the output of the nested loops in the form {index triple, center of hexahedron, hexahedron}, which allows us to, say, display only those hexahedra with centers below a certain level.

By example we set n=15, k=3. We call k “karm” to avoid confusion with the nested loops, {i,j,k}.

(*set n, pitch, edge length, make vectors, find triples*)
n=15;
pitch=ArcTan[ Sqrt[1/2] ];(*at this pitch, orthogonal shadow of hypercube*)
edge=1;
v=vectors[n,1/Degree pitch,edge];(*correct for radian measure*)
karm=3;
kstep=n/karm;

We construct all k of the k-regions.

pt=Chop[ (Plus@@v)/2 ];(*sum of vectors, divided by 2*)

allregions=Table[

Table[
  Table[If[i==1,mark=-pt,mark=-pt+(Plus@@v[[ Range[sk+1,i-1] ]]) ];
	Table[
	{
	1+Mod[{i-1,j-1,k-1}, n],(*index triple*)
	mark+(Plus@@v[[ 1+Mod[{i-1,j-1,k-1}, n] ]]/2),(*hexahedron center*)
	newhex[ v[[ 1+Mod[{i-1,j-1,k-1}, n] ]], mark](*hexahedron*)
	},
	{k,j+1,n+sk}],
  {j, i+1, sk+(n/karm) }],
{i, sk+1, sk+((n/karm)-1) }],

{sk,0,n-kstep,kstep}];

We flatten the nested loops and then, for no particular reason, remove coincident faces from the first region, and display the resulting “region-polyhedron,” identical in outward appearance to the region of hexahedra. It is interesting to use the colorSort function on such a region-polyhedron (this region-polyhedron has been named “regpoly”; go to the “Color by Area” section, above).

(*flatten nested list of {triples, centers, hexahedra}*)
hexreg=Map[Flatten[#,2]&, allregions];

f1=dePair[ parse[ hexreg[[1]] ] ];
Show[regpoly=Graphics3D[Polygon/@f1],
Axes->True,
ViewPoint->{0,0,-3}];

Make the k-armed spirallohedron, then select all regional hexahedra with z coordinate less than zero, and display spirallohedron with sublist of regional hexahedra.

a=Spirallohedron[n,karm];

h1=Table[
Map[#[[3]]&, Select[hexreg[[k]], #[[2,3]]<0&] ],
{k,Length[hexreg]}];

Show[xaq={a,h1},
Axes->True,
ViewPoint->2.25*unit[{0,0,1}]];