```/*
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)
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.

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.

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

*/

#include <math.h>

public:
double        Coords;
double        A ()  const { return Coords; }
double        B ()  const { return Coords; }
double        C ()  const { return Coords; }
double        D ()  const { return Coords; }
double        &A()        { return Coords; }
double        &B()        { return Coords; }
double        &C()        { return Coords; }
double        &D()        { return Coords; }

void          ZeroSumNormalize();

// dist and dot product are scaled so that basis vectors have unity length

double        DistFrom      (const Quadray &Other) const;
double        DistFromOrigin() const;

};

{
// 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;
}

{
// 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;
}

{
// 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

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

{
// 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

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

{
// 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

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

{
// 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;

Cross = 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 = 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 = 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 = 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(
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.;
}
};

{
// 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;

Rotated = QX.A();
Rotated = SumOfProducts(QX,0.,RC.F,RC.H,RC.G);
Rotated = SumOfProducts(QX,0.,RC.G,RC.F,RC.H);
Rotated = SumOfProducts(QX,0.,RC.H,RC.G,RC.F);

for (int I = 0; I < 4; I++) Coords[I] = Rotated[I];
}

```