/*
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];
}