/*

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.
Quadrays have been independently devised by several people.
John H. Conway recognized quadrays as being an instance of
tetrahedral coordinates (the general n-dimensional term is 
"simplicial coordinates").

In August 1999, I derived native quadray formulas for a few common
geometrical operations.  In the process I (re-)discovered a 
correspondence between quadrays and Cartesian 4-space coordinates 
that allows quick derivation of many quadray formulas.  That
correspondence is described in another article.

A few quadray formulas are described below in C++, with comments about
the method used and how I derived it.  (I haven't included trivial
methods like translation and scaling.)  I chose program code because 
it's unambiguous, and because I'm too lazy to make nice GIFs for 
expressions with math symbols.  Note that I was also too lazy to 
code the cross product method, and just gave a formula in the comments.

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

   static double DistFromOrigin(const Quadray &QX);

   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 Quadray &QX)
{
   // returns distance from *this to (0,0,0,0)
   // method:  normalize, apply traditional Cartesian formula, and scale
   // how this was derived:  see http://www.minortriad.com/q4d.html

   Quadray   QXN;

   QXN.ZeroSumNormalize(QX);
   double Dist = 0.;
   for (int I = 0; I < 4; I++) Dist += QXN.Coords[I] * QXN.Coords[I];
   return sqrt(Dist * 4. / 3.);
}

double Quadray::DotProduct(const Quadray &QX,const Quadray &QY)
{
   // returns QX . QY  (traditional dot product)
   // 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  (traditional 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 (code to calculate it isn't written yet).
   // The top row contains 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() |
}

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.
   //
   // Adapting this code for methods that rotate about B, C, and D
   // is left to the reader as an exercise.
   //
   // How this was derived:  probably close to the most tedious way possible;
   //                        please find me an elegant derivation for this!
   //
   // Don't use to rotate in place; i.e. it won't work if (this == &QX).

   RotationCoeffs    RC(Theta);

   A() = QX.A();
   B() = SumOfProducts(QX,0.,RC.F,RC.H,RC.G);
   C() = SumOfProducts(QX,0.,RC.G,RC.F,RC.H);
   D() = SumOfProducts(QX,0.,RC.H,RC.G,RC.F);
}