Clearly vector addition and negation are being handled according to some different (non-XYZ) rules. Specifically, to obtain the normal form of a quadray, we subtract its lowest coordinate from the other 4, an operation with no net effect because (d,d,d,d)= (0,0,0,0). Here's an algorithm:
def norm(self,plist):
"""Normalize such that 4-tuple all non-negative members."""
return tuple(map(sub,plist,[min(plist)]*4))
In this method, plist is a passed 4-tuple, e.g. (-1,-2,0,4), with min(plist) returning the lowest coordinate e.g. -2. The expression [min(plist)]*4, returns a 4-element list, e.g. [-2,-2,-2,-2]. The "map" operator then pairs each member of plist with a corresponding minimum, and does subtraction e.g. (-1-(-2), -2-(-2), 0-(-2),4-(-2)) = (1,0,2,6). The resulting tuple is normalized form of the quadray. This method is invoked whenever a quadray is initialized, perhaps as a result of vector addition or subtraction.
NOMENCLATUREWhether we use the term "quadray" is not especially critical and depends on the namespace of the author. David Chako originally introduced his 4-tuple vector algebra on Synergetics-L in December, 1996, without using this term, whereas D. Lloyd Jarmusch mentioned using "quadrays" as early as 1981 for essentially the same apparatus. Josef Hasslberger made a similar beginning under the heading of "tetra space coordinates."
At one point I played with calling the six Cartesian vectors (i,j,k,-i,-j,-k) "e-rays", because they go from (0,0,0) at the center of a tetrahedron through its six mid-edges (hence "e" for edge).
By extension, quadrays could be called "v-rays", because they go through the tetrahedron's four vertices, or, equivalently, "f-rays" through the face-centers (Hasslberger's depiction) -- "equivalently" because the tetrahedron is self-dual, i.e. "through the faces" is "through the vertices" of the dual, also a regular tetrahedron.
More important than what they're called is how they work, and here again we have some divergence in their implementation.
QUADRAYS, TETRAYS & SIMPLICIAL COORDINATESThe most detailed elaboration of quadrays of which I am personally aware occured on Synergetics-L, subsequent to Chako's getting the ball rolling, and in this context the research forked into quadrays and tetrays as distinct topics, the difference being that tetrays use the identity (1,1,1,1) to tick off an additional, fifth "temporal direction".
Possibly we'll end up classifying quadrays as a kind of "simplicial coordinate system", where "simplicial" derives from "simplex", another name for the tetrahedron. The differences in implementation may reduce to conventions around how we choose to normalize our 4-tuples.
ALTERNATIVE NORMALIZATIONSQuadrays as originally defined by Chako used only non-negative integers, meaning {a,b,c,0} did not sum to any specific constant k = a+b+c+d. On the other hand, because (d,d,d,d) may be added to (a,b,c,d) with no net effect -- (d,d,d,d)= (0,0,0,0) -- we have the option to normalize to any integer we like. The same rules apply if we allow (a,b,c,d) to be real or floating point numbers.
Tom Ace used the k=0 normalization to derive a dot product, and a more streamlined distance method, proving the importance of remaining flexible about how we normalize.
def norm0(self):
"""Normalize such that sum of 4-tuple members = 0"""
q = self.quadray()
return tuple(map(sub,q,[reduce(add,q)/4.0]*4))
def dot(self,v1):
"""Return the dot product of self with another vector.
return a scalar"""
return 0.5 * reduce(add,map(mul,v1.xyz,self.xyz))
def length(self):
"""Return this vector's length"""
return self.dot(self) ** 0.5
MORE QUADRAY METHODS
Ace also came up with a cross product method and 4x4 rotation matrices, with (1,0,0,0)... (0,0,0,1) as the axes of rotation, none of which depend on normalization to zero. Here's a cross product method:
def cross(self,v1):
"""Return the cross product of self with another vector.
return a Qvector"""
A=Qvector((1,0,0,0))
B=Qvector((0,1,0,0))
C=Qvector((0,0,1,0))
D=Qvector((0,0,0,1))
a1,b1,c1,d1 = v1.quadray()
a2,b2,c2,d2 = self.quadray()
k= (2.0**0.5)/4.0
sum = (A*c1*d2 - A*d1*c2 - A*b1*d2 + A*b1*c2
+ A*b2*d1 - A*b2*c1 - B*c1*d2 + B*d1*c2
+ b1*C*d2 - b1*D*c2 - b2*C*d1 + b2*D*c1
+ a1*B*d2 - a1*B*c2 - a1*C*d2 + a1*D*c2
+ a1*b2*C - a1*b2*D - a2*B*d1 + a2*B*c1
+ a2*C*d1 - a2*D*c1 - a2*b1*C + a2*b1*D)
return k*sum
As Tom Ace explains at his website, the above is simply an expansion of the following determinant:
A=(1,0,0,0) B=(0,1,0,0) C=(0,0,1,0) D=(0,0,0,1) | A B C D | | k k k k | | a1 b1 c1 d1 | | a2 b2 c2 d2 |VOLUME METHODS
Ace also developed a volume method consistent with the usual equivalence in XYZ of a scalar triple product with a determinant i.e.
| Ax Ay Ax |
A dot (B cross C) = | Bx By Bz | (A,B,C are vectors)
| Cx Cy Cz |
The scalar triple product likewise gives the volume of a parallelepiped with edges ABC from a common vertex, or of a corresponding tetrahedron, if we multiply the above scalar by k=(1/6) -- plus take the absolute value if we want to eliminate a possible negative sign.
The volume of tetrahedron ABCD using quadray coordinates is likewise |k*det[M]| where AB, AC, AD = quadrays (a1,b1,c1,d1),(a2,b2,c2,d2),(a3,b3,c3,d3) and M is the matrix:
| 1 1 1 1 |
| a1 b1 c1 d1 |
| a2 b2 c2 d2 |
| a3 b3 c3 d3 |
If we set k=(1/4), then the volume of a regular tetrahedron with all edges 1 would be 1 (see below). An regular icosahedron consisting of 20 irregular tetrahedra, with radials of ~0.951 and outer edges = 1, would have a volume of about 18.51.
>>> ico = Icosahedron() # create Icosahedron object >>> ico.vertices["O1"].length() # a center-to-vertex vector 0.95105651629515353 >>> ico.volume # compute volume 18.512295868219159
The Icosahedron's constructor method sets the starting volume
using the expression:
def vol1(qvectors): """Return the volume of a tetrahedron, given 3 quadrays thanks to Tom Ace""" return (1./4.)*abs(det(qvectors))
Vectors define direction. Any 3 non-coplanar edges of a tetrahedron, expressed as vectors, will provide sufficient information to determine the other three edges. These vectors may either share a common vertex, or define an "open triangle" zig-zag. Either way, they provide sufficient input for various vector-based volume methods.
In contrast, simple edge lengths, being scalar quantities, do not define direction. To specify three non-coplanar lengths does not uniquely determine the remaining three edges. If we label a tetrahedron's six edges a,b,c,d,e,f -- a,b,c connecting a common vertex to an opposite triangle with edges d,e,f -- then we will need all six values in order to compute a specific volume.
Leonhard Euler long ago solved the problem of how to express a tetrahedron's volume solely in terms of its edge lengths. Java programmer Gerald de Jong rederived the method and expressed his result in the form of a Java applet.
The computation involves using "closed triangles" (edges from the same face), "open triangles" (3-edge zig-zags), and "opposite pairs" (edges with no endpoints in common). Gerald's expression is also specifically designed to return a volume of 1 when the six edges of the tetrahedron are 1.
Here's a Python version of the algorithm:
def vol2(edges): """Return the volume of a tetrahedron, 6 edge lengths Thanks to Leonhard Euler and Gerald de Jong """ a2,b2,c2,d2,e2,f2 = map(mul,edges,edges) open = ( f2 * a2 * b2 + d2 * a2 * c2 + a2 * b2 * e2 + c2 * b2 * d2 + e2 * c2 * a2 + f2 * c2 * b2 + e2 * d2 * a2 + b2 * d2 * f2 + b2 * e2 * f2 + d2 * e2 * c2 + a2 * f2 * e2 + d2 * f2 * c2 ) closed = ( a2 * b2 * d2 + d2 * e2 * f2 + b2 * c2 * e2 + a2 * c2 * f2 ) oppos = ( a2 * e2 * (a2 + e2) + b2 * f2 * (b2 + f2) + c2 * d2 * (c2 + d2)) return (abs(open - closed - oppos)/2.0)**0.5
This method has nothing specifically to do with quadrays, since
length is an attribute of the generic Vector class. I mention it
more because of its historical role in the quadray thread on
David Chako conjectured that all tetrahedra with vertices at the centers of fcc spheres might have whole number volumes, provided the initial tetrahedron -- defined by four intertangent fcc spheres -- were defined as volumetric unity.
Until Robert Gray developed an algebraic proof for this conjecture, de Jong's Java applet and a MathCad worksheet featuring Euler's original formula, provided an initial source of empirical evidence for its validity (i.e. failed to turn up any counter-examples). Other volume methods would have served, but this was how events actually unfolded.
In terms of quadrays, fcc sphere centers may be expressed as sums of vectors with coordinates {2,1,1,0}. For example, if s were an fcc sphere center, (2,1,1,0) + (1,1,2,0) + (1,2,0,1) = (4,4,3,1) = (3,3,2,0) would be another fcc sphere center relative to s.
I will return to quadrays in a sphere packing context in the
notes on crystallography below.
