#include <string.h>
#include <math.h>

typedef double DBL;

typedef DBL VECTOR [3];

typedef struct Ray_Struct RAY;

struct Ray_Struct
{
  VECTOR Initial;
  VECTOR Direction;
  int Index;
  unsigned long Optimisiation_Flags;
};

enum
{
	X = 0,
	Y = 1,
	Z = 2,
	T = 3
};

typedef struct byte_xyz BYTE_XYZ;

struct byte_xyz {
  unsigned char x, y, z;
};

#define Assign_Vector(d,s)  memcpy((d),(s),sizeof(VECTOR))

#define VCross(a,b,c) {(a)[X]=(b)[Y]*(c)[Z]-(b)[Z]*(c)[Y]; \
                       (a)[Y]=(b)[Z]*(c)[X]-(b)[X]*(c)[Z]; \
                       (a)[Z]=(b)[X]*(c)[Y]-(b)[Y]*(c)[X];}

#define VDot(a, b, c) {a=(b)[X]*(c)[X]+(b)[Y]*(c)[Y]+(b)[Z]*(c)[Z];}

#define VNormalizeEq(a) {DBL _tmp=1.0/sqrt((a)[X]*(a)[X]+(a)[Y]*(a)[Y]+(a)[Z]*(a)[Z]);(a)[X]*=_tmp;(a)[Y]*=_tmp;(a)[Z]*=_tmp;}

#define VAddScaledEq(V, k, V2)  \
  { (V)[X] += (k) * (V2)[X];    \
    (V)[Y] += (k) * (V2)[Y];    \
    (V)[Z] += (k) * (V2)[Z]; }

static void VUnpack(VECTOR dest_vec, BYTE_XYZ * pack_vec)
{
  dest_vec[X] = ((double)pack_vec->x * (1./ 255.))*2.-1.;
  dest_vec[Y] = ((double)pack_vec->y * (1./ 255.))*2.-1.;
  dest_vec[Z] = ((double)pack_vec->z * (1./ 255.));

  VNormalizeEq(dest_vec);   /* already good to about 1%, but we can do better */
}

extern BYTE_XYZ rad_samples[];

void ChooseRay(RAY *NewRay, VECTOR Normal, RAY *Ray, VECTOR Raw_Normal, int WhichRay)
{
  VECTOR random_vec, up, n2, n3;
  int i;
  DBL /*n,*/ NRay_Direction;

  Assign_Vector(NewRay->Direction, Normal);

  if ( fabs(fabs(NewRay->Direction[Z])- 1.) < .1 ) {
    /* too close to vertical for comfort, so use cross product with horizon */
    up[X] = 0.; up[Y] = 1.; up[Z] = 0.;
  }
  else
  {
    up[X] = 0.; up[Y] = 0.; up[Z] = 1.;
  }

  VCross(n2, NewRay->Direction, up);  VNormalizeEq(n2);
  VCross(n3, NewRay->Direction, n2);  VNormalizeEq(n3);

  /*i = (int)(FRAND()*1600);*/
  i = WhichRay;
  WhichRay = (WhichRay + 1) % 1600;

  VUnpack(random_vec, &rad_samples[i]);

  if ( fabs(NewRay->Direction[Z] - 1.) < .001 )         /* pretty well straight Z, folks */
  {
    /* we are within 1/20 degree of pointing in the Z axis. */
    /* use all vectors as is--they're precomputed this way */
    Assign_Vector(NewRay->Direction, random_vec);
  }
  else
  {
    NewRay->Direction[X] = n2[X]*random_vec[X] + n3[X]*random_vec[Y] + NewRay->Direction[X]*random_vec[Z];
    NewRay->Direction[Y] = n2[Y]*random_vec[X] + n3[Y]*random_vec[Y] + NewRay->Direction[Y]*random_vec[Z];
    NewRay->Direction[Z] = n2[Z]*random_vec[X] + n3[Z]*random_vec[Y] + NewRay->Direction[Z]*random_vec[Z];
  }

  /* if our new ray goes through, flip it back across raw_normal */

  VDot(NRay_Direction, NewRay->Direction, Raw_Normal);
  if (NRay_Direction < 0.0)
  {
    /* subtract 2*(projection of NRay.Direction onto Raw_Normal)
       from NRay.Direction */
    DBL Proj;
    Proj = NRay_Direction * -2;
    VAddScaledEq(NewRay->Direction, Proj, Raw_Normal);
  }

  VNormalizeEq(NewRay->Direction);
}