/*
Quadrays are a 4-coordinate system for mapping space, described in detail
on Kirby Urner's page at http://www.grunch.net/synergetics/quadintro.html.

Briefly:  quadrays use four basis vectors, configured in directions
from the center of a regular tetrahedron towards its four vertices.
A point at (a,b,c,d) is a linear combination of the four basis vectors.
These are also known as tetrahedral coordinates; the general 
n‑dimensional term is simplicial coordinates.

As the four quadray basis vectors sum to zero, multiples of (1,1,1,1)
may be added to quadray coordinates (a,b,c,d) without changing the 
point in space being referred to.  Various methods of normalizing
coordinates are possible.  Kirby Urner's quadintro.html page discusses
a form of normalization that minimizes the values of a, b, c, and d
while keeping them all nonnegative.  Choosing (a,b,c,d) such that
a+b+c+d=1 gives barycentric coordinates.  Choosing (a,b,c,d) such that
a+b+c+d=0 facilitates computation by exploiting an isomorphism
described on my page at http://minortriad.com/q4d.html.

A few quadray formulas are coded below in C++, with comments about
the method used and how I derived it.  (I haven't included trivial
methods like translation and scaling.)  This document, saved as text,
compiles as ANSI C++.  There is no warranty of fitness for purpose,
nor of anything else.

Kirby's page at http://www.grunch.net/synergetics/quadphil.html considers 
some philosophical and aesthetic implications of quadrays.  From my 
experience, one doesn't have to agree with all of Kirby's points to 
find quadrays interesting to play with.

Tom Ace

return to Tom's home page */ #include <math.h> class Quadray { public: double Coords[4]; double A () const { return Coords[0]; } double B () const { return Coords[1]; } double C () const { return Coords[2]; } double D () const { return Coords[3]; } double &A() { return Coords[0]; } double &B() { return Coords[1]; } double &C() { return Coords[2]; } double &D() { return Coords[3]; } void ZeroSumNormalize(); void ZeroSumNormalize(const Quadray &QX); // dist and dot product are scaled so that basis vectors have unity length double DistFrom (const Quadray &Other) const; double DistFromOrigin() const; static double DotProduct (const Quadray &QX,const Quadray &QY); void CrossProduct(const Quadray &QX,const Quadray &QY); void RotateAboutA(const Quadray &QX,double Theta); }; void Quadray::ZeroSumNormalize() { // normalizes *this to A+B+C+D==0 double Mean = (A() + B() + C() + D()) / 4.; for (int I = 0; I < 4; I++) Coords[I] -= Mean; } void Quadray::ZeroSumNormalize(const Quadray &QX) { // Sets *this to a quadray equivalent to QX, but with A+B+C+D==0 double Mean = (QX.A() + QX.B() + QX.C() + QX.D()) / 4.; for (int I = 0; I < 4; I++) Coords[I] = QX.Coords[I] - Mean; } double Quadray::DistFromOrigin() const { // returns distance from *this to (0,0,0,0) // method: normalize, use Pythagorean distance dist formula, and scale // how this was derived: see http://www.minortriad.com/q4d.html Quadray QN; QN.ZeroSumNormalize(*this); double DistSq = 0.; for (int I = 0; I < 4; I++) DistSq += QN.Coords[I] * QN.Coords[I]; return sqrt(DistSq * 4. / 3.); } double Quadray::DistFrom(const Quadray &Other) const { // returns distance from *this to Other // method: normalize, use Pythagorean distance formula, and scale // how this was derived: see http://www.minortriad.com/q4d.html Quadray QN; Quadray OtherN; QN .ZeroSumNormalize(*this); OtherN.ZeroSumNormalize(Other); double DistSq = 0.; for (int I = 0; I < 4; I++) { double Delta = QN.Coords[I] - OtherN.Coords[I]; DistSq += Delta * Delta; } return sqrt(DistSq * 4. / 3.); } double Quadray::DotProduct(const Quadray &QX,const Quadray &QY) { // returns QX . QY (dot product of two vectors) // method: normalize, apply traditional Cartesian formula, and scale // how this was derived: see http://www.minortriad.com/q4d.html Quadray QXN,QYN; QXN.ZeroSumNormalize(QX); QYN.ZeroSumNormalize(QY); double Dot = 0.; for (int I = 0; I < 4; I++) Dot += QXN.Coords[I] * QYN.Coords[I]; return Dot * 4. / 3.; } void Quadray::CrossProduct(const Quadray &QX,const Quadray &QY) { // sets *this to QX x QY (vector cross product) // method: calculate a determinant (by Laplace development on the top row) // (somewhat reminiscent of Cartesian cross product determinant) // how this was derived: by intuition; // verified (and k determined) empirically; // can be motivated by the 4-D correspondence // The determinant is as follows. The top row consists of the // four basis vectors; all other elements in the matrix are scalars. double k = sqrt(3.) / 3.; // | A B C D | // | | // | k k k k | // | | // | QX.A() QX.B() QX.C() QX.D() | // | | // | QY.A() QY.B() QY.C() QY.D() | double Cross[4]; Cross[0] = k * ( QX.C() * QY.D() - QX.D() * QY.C() + QX.D() * QY.B() - QX.B() * QY.D() + QX.B() * QY.C() - QX.C() * QY.B()); Cross[1] = k * ( QX.D() * QY.C() - QX.C() * QY.D() + QX.A() * QY.D() - QX.D() * QY.A() + QX.C() * QY.A() - QX.A() * QY.C()); Cross[2] = k * ( QX.A() * QY.B() - QX.B() * QY.A() + QX.D() * QY.A() - QX.A() * QY.D() + QX.B() * QY.D() - QX.D() * QY.B()); Cross[3] = k * ( QX.C() * QY.B() - QX.B() * QY.C() + QX.A() * QY.C() - QX.C() * QY.A() + QX.B() * QY.A() - QX.A() * QY.B()); for (int I = 0; I < 4; I++) Coords[I] = Cross[I]; } static inline double SumOfProducts( const Quadray &QX, double A,double B,double C,double D) { return QX.A() * A + QX.B() * B + QX.C() * C + QX.D() * D; } class RotationCoeffs { public: double F,G,H; RotationCoeffs(double Theta) // ctor { static const double RAD_PER_DEG = M_PI / 180.; F = (2. * cos((Theta ) * RAD_PER_DEG) + 1.) / 3.; G = (2. * cos((Theta - 120.) * RAD_PER_DEG) + 1.) / 3.; H = (2. * cos((Theta + 120.) * RAD_PER_DEG) + 1.) / 3.; } }; void Quadray::RotateAboutA(const Quadray &QX,double Theta) { // sets *this equal to QX, rotated Theta degrees about A // method: multiply the following matrix by QX: // // 1 0 0 0 where // 0 F H G F = (2 * cos(Theta ) + 1) / 3 // 0 G F H G = (2 * cos(Theta - 120) + 1) / 3 // 0 H G F H = (2 * cos(Theta + 120) + 1) / 3 // // This is an orthogonal matrix (its transpose is its inverse). // Note that F + G + H == 1 and F*F + G*G + H*H == 1. // // How this was derived: a bunch of algebra and trig. I later // developed a method for rotation about any // desired axis but I haven't coded that yet. RotationCoeffs RC(Theta); double Rotated[4]; Rotated[0] = QX.A(); Rotated[1] = SumOfProducts(QX,0.,RC.F,RC.H,RC.G); Rotated[2] = SumOfProducts(QX,0.,RC.G,RC.F,RC.H); Rotated[3] = SumOfProducts(QX,0.,RC.H,RC.G,RC.F); for (int I = 0; I < 4; I++) Coords[I] = Rotated[I]; }