void VECTOR3D_Scale(float k, VECTOR3D_PTR va)
{
// this function scales a vector by the constant k,
// and modifies the original

// multiply each component by scaling factor
va->x*=k;
va->y*=k;
va->z*=k;

} // end VECTOR3D_Scale

/////////////////////////////////////////////////////////////

void VECTOR3D_Scale(float k, VECTOR3D_PTR va, VECTOR3D_PTR vscaled)
{
// this function scales a vector by the constant k,
// leaves the original unchanged, and returns the result
// in vres as well as on the stack

// multiply each component by scaling factor
vscaled->x = k*va->x;
vscaled->y = k*va->y;
vscaled->z = k*va->z;

} // end VECTOR3D_Scale

//////////////////////////////////////////////////////////////

float VECTOR3D_Dot(VECTOR3D_PTR va, VECTOR3D_PTR vb)
{
// computes the dot product between va and vb
return( (va->x * vb->x) + (va->y * vb->y) + (va->z * vb->z) );
} // end VECTOR3D_DOT
/////////////////////////////////////////////////////////////

void VECTOR3D_Cross(VECTOR3D_PTR va, VECTOR3D_PTR vb, VECTOR3D_PTR vn)
{
// this function computes the cross product between va and vb
// and returns the vector that is perpendicular to each in vn

vn->x =  ( (va->y * vb->z) - (va->z * vb->y) );
vn->y = -( (va->x * vb->z) - (va->z * vb->x) );
vn->z =  ( (va->x * vb->y) - (va->y * vb->x) ); 

} // end VECTOR3D_Cross

/////////////////////////////////////////////////////////////

VECTOR3D VECTOR3D_Cross(VECTOR3D_PTR va, VECTOR3D_PTR vb)
{
// this function computes the cross product between va and vb
// and returns the vector that is perpendicular to each

VECTOR3D vn;

vn.x =  ( (va->y * vb->z) - (va->z * vb->y) );
vn.y = -( (va->x * vb->z) - (va->z * vb->x) );
vn.z =  ( (va->x * vb->y) - (va->y * vb->x) ); 

// return result
return(vn);

} // end VECTOR3D_Cross

//////////////////////////////////////////////////////////////

float VECTOR3D_Length(VECTOR3D_PTR va)
{
// computes the magnitude of a vector, slow

return( (float)sqrtf(va->x*va->x + va->y*va->y + va->z*va->z) );

} // end VECTOR3D_Length

///////////////////////////////////////////////////////////////

float VECTOR3D_Length_Fast(VECTOR3D_PTR va)
{
// computes the magnitude of a vector using an approximation
// very fast
return( Fast_Distance_3D(va->x, va->y, va->z) );

} // end VECTOR3D_Length_Fast

///////////////////////////////////////////////////////////////

void VECTOR3D_Normalize(VECTOR3D_PTR va)
{
// normalizes the sent vector in placew

// compute length
float length = sqrtf(va->x*va->x + va->y*va->y + va->z*va->z);

// test for zero length vector
// if found return zero vector
if (length < EPSILON_E5) 
   return;

float length_inv = 1/length;

// compute normalized version of vector
va->x*=length_inv;
va->y*=length_inv;
va->z*=length_inv;

} // end VECTOR3D_Normalize

///////////////////////////////////////////////////////////////

void VECTOR3D_Normalize(VECTOR3D_PTR va, VECTOR3D_PTR vn)
{
// normalizes the sent vector and returns the result in vn

VECTOR3D_ZERO(vn);

// compute length
float length = VECTOR3D_Length(va);

// test for zero length vector
// if found return zero vector
if (length < EPSILON_E5) 
   return;

float length_inv = 1.0/length;

// compute normalized version of vector
vn->x = va->x*length_inv;
vn->y = va->y*length_inv;
vn->z = va->z*length_inv;

} // end VECTOR3D_Normalize

///////////////////////////////////////////////////////////////

void VECTOR3D_Build(VECTOR3D_PTR init,VECTOR3D_PTR term,VECTOR3D_PTR result)
{
// this function creates a vector from two vectors (or points)
//  in 3D space
result->x = term->x - init->x;
result->y = term->y - init->y;
result->z = term->z - init->z;

} // end VECTOR3D_Build

/////////////////////////////////////////////////////////////

float VECTOR3D_CosTh(VECTOR3D_PTR va, VECTOR3D_PTR vb)
{
// this function returns the cosine of the angle between
// two vectors. Note, we could compute the actual angle,
// many many times, in further calcs we will want ultimately
// compute cos of the angle, so why not just leave it!
return(VECTOR3D_Dot(va,vb)/(VECTOR3D_Length(va)*VECTOR3D_Length(vb)));

} // end VECTOR3D_CosTh

///////////////////////////////////////////////////////////////

void VECTOR3D_Print(VECTOR3D_PTR va, char *name="v")
{
// this function prints out a VECTOR3D

Write_Error("\n%s=[",name);
for (int index=0; index<3; index++)
    Write_Error("%f, ",va->M[index]);
Write_Error("]");

} // end VECTOR3D_Print

////////////////////////////////////////////////////////////////

// these are the 4D version of the vector functions, they
// assume that the vectors are 3D with a w, so w is left
// out of all the operations

void VECTOR4D_Build(VECTOR4D_PTR init, VECTOR4D_PTR term, VECTOR4D_PTR result)
{
// build a 4d vector
result->x = term->x - init->x;
result->y = term->y - init->y;
result->z = term->z - init->z;
result->w = 1;

} // end VECTOR4D_Build

////////////////////////////////////////////////////////////////

void VECTOR4D_Add(VECTOR4D_PTR va, VECTOR4D_PTR vb, VECTOR4D_PTR vsum)
{
// this function adds va+vb and return it in vsum
vsum->x = va->x + vb->x;
vsum->y = va->y + vb->y;
vsum->z = va->z + vb->z;
vsum->w = 1;

} // end VECTOR4D_Add

////////////////////////////////////////////////////////////

VECTOR4D VECTOR4D_Add(VECTOR4D_PTR va, VECTOR4D_PTR vb)
{
// this function adds va+vb and returns the result on 
// the stack
VECTOR4D vsum;

vsum.x = va->x + vb->x;
vsum.y = va->y + vb->y;
vsum.z = va->z + vb->z;
vsum.w = 1;

// return result
return(vsum);

} // end VECTOR4D_Add

////////////////////////////////////////////////////////////

void VECTOR4D_Sub(VECTOR4D_PTR va, VECTOR4D_PTR vb, VECTOR4D_PTR vdiff)
{
// this function subtracts va-vb and return it in vdiff
// the stack
vdiff->x = va->x - vb->x;
vdiff->y = va->y - vb->y;
vdiff->z = va->z - vb->z;
vdiff->w = 1;

} // end VECTOR4D_Sub

////////////////////////////////////////////////////////////

VECTOR4D VECTOR4D_Sub(VECTOR4D_PTR va, VECTOR4D_PTR vb)
{
// this function subtracts va-vb and returns the result on 
// the stack
VECTOR4D vdiff;

vdiff.x = va->x - vb->x;
vdiff.y = va->y - vb->y;
vdiff.z = va->z - vb->z;
vdiff.w = 1;

// return result
return(vdiff);                      

} // end VECTOR4D_Sub

////////////////////////////////////////////////////////////

void VECTOR4D_Scale(float k, VECTOR4D_PTR va)
{
// this function scales a vector by the constant k,
// in place , note w is left unchanged

// multiply each component by scaling factor
va->x*=k;
va->y*=k;
va->z*=k;
va->w = 1;

} // end VECTOR4D_Scale

/////////////////////////////////////////////////////////////

void VECTOR4D_Scale(float k, VECTOR4D_PTR va, VECTOR4D_PTR vscaled)
{
// this function scales a vector by the constant k,
// leaves the original unchanged, and returns the result
// in vres as well as on the stack

// multiply each component by scaling factor
vscaled->x = k*va->x;
vscaled->y = k*va->y;
vscaled->z = k*va->z;
vscaled->w = 1;

} // end VECTOR4D_Scale

//////////////////////////////////////////////////////////////

float VECTOR4D_Dot(VECTOR4D_PTR va, VECTOR4D_PTR vb)
{
// computes the dot product between va and vb
return( (va->x * vb->x) + (va->y * vb->y) + (va->z * vb->z) );
} // end VECTOR4D_DOT

/////////////////////////////////////////////////////////////

void VECTOR4D_Cross(VECTOR4D_PTR va, 
                    VECTOR4D_PTR vb,
                    VECTOR4D_PTR vn)
{
// this function computes the cross product between va and vb
// and returns the vector that is perpendicular to each in vn

vn->x =  ( (va->y * vb->z) - (va->z * vb->y) );
vn->y = -( (va->x * vb->z) - (va->z * vb->x) );
vn->z =  ( (va->x * vb->y) - (va->y * vb->x) ); 
vn->w = 1;

} // end VECTOR4D_Cross

/////////////////////////////////////////////////////////////

VECTOR4D VECTOR4D_Cross(VECTOR4D_PTR va, VECTOR4D_PTR vb)
{
// this function computes the cross product between va and vb
// and returns the vector that is perpendicular to each

VECTOR4D vn;

vn.x =  ( (va->y * vb->z) - (va->z * vb->y) );
vn.y = -( (va->x * vb->z) - (va->z * vb->x) );
vn.z =  ( (va->x * vb->y) - (va->y * vb->x) ); 
vn.w = 1;

// return result
return(vn);

} // end VECTOR4D_Cross

//////////////////////////////////////////////////////////////

float VECTOR4D_Length(VECTOR4D_PTR va)
{
// computes the magnitude of a vector, slow

return(sqrtf(va->x*va->x + va->y*va->y + va->z*va->z) );

} // end VECTOR4D_Length

///////////////////////////////////////////////////////////////

float VECTOR4D_Length_Fast(VECTOR4D_PTR va)
{
// computes the magnitude of a vector using an approximation
// very fast
return( Fast_Distance_3D(va->x, va->y, va->z) );

} // end VECTOR4D_Length_Fast

///////////////////////////////////////////////////////////////

void VECTOR4D_Normalize(VECTOR4D_PTR va)
{
// normalizes the sent vector and returns the result

// compute length
float length = sqrtf(va->x*va->x + va->y*va->y + va->z*va->z);

// test for zero length vector
// if found return zero vector
if (length < EPSILON_E5) 
   return;

float length_inv = 1.0/length;

// compute normalized version of vector
va->x*=length_inv;
va->y*=length_inv;
va->z*=length_inv;
va->w = 1;

} // end VECTOR4D_Normalize

///////////////////////////////////////////////////////////////

void VECTOR4D_Normalize(VECTOR4D_PTR va, VECTOR4D_PTR vn)
{
// normalizes the sent vector and returns the result in vn

VECTOR4D_ZERO(vn);

// compute length
float length = sqrt(va->x*va->x + va->y*va->y + va->z*va->z);

// test for zero length vector
// if found return zero vector
if (length < EPSILON_E5) 
   return;

float length_inv = 1.0/length;

// compute normalized version of vector
vn->x = va->x*length_inv;
vn->y = va->y*length_inv;
vn->z = va->z*length_inv;
vn->w = 1;

} // end VECTOR4D_Normalize

///////////////////////////////////////////////////////////////

float VECTOR4D_CosTh(VECTOR4D_PTR va, VECTOR4D_PTR vb)
{
// this function returns the cosine of the angle between
// two vectors. Note, we could compute the actual angle,
// many many times, in further calcs we will want ultimately
// compute cos of the angle, so why not just leave it!
return(VECTOR4D_Dot(va,vb)/(VECTOR4D_Length(va)*VECTOR4D_Length(vb)));

} // end VECTOR4D_CosTh

////////////////////////////////////////////////////////////

void VECTOR4D_Print(VECTOR4D_PTR va, char *name="v")
{
// this function prints out a VECTOR4D

Write_Error("\n%s[",name);
for (int index=0; index<4; index++)
    Write_Error("%f, ",va->M[index]);
Write_Error("]");

} // end VECTOR4D_Print

////////////////////////////////////////////////////////////////

void Mat_Init_2X2(MATRIX2X2_PTR ma, 
                  float m00, float m01,
                  float m10, float m11)

{
// this function fills a 2x2 matrix with the sent data in 
// row major form
ma->M00 = m00; ma->M01 = m01; 
ma->M10 = m10; ma->M11 = m11; 

} // end Mat_Init_2X2

/////////////////////////////////////////////////////////////////

void Mat_Add_2X2(MATRIX2X2_PTR ma, MATRIX2X2_PTR mb, MATRIX2X2_PTR msum)
{
// this function adds two 2x2 matrices together and stores
// the result in msum
msum->M00 = ma->M00+mb->M00;
msum->M01 = ma->M01+mb->M01;
msum->M10 = ma->M10+mb->M10;
msum->M11 = ma->M11+mb->M11;

} // end Mat_Add_2X2

/////////////////////////////////////////////////////////////////

void Mat_Mul_2X2(MATRIX2X2_PTR ma, MATRIX2X2_PTR mb, MATRIX2X2_PTR mprod)