/*
	Rogue Trader.
	Copyright (c) 2007, Darragh Coy.
	All rights reserved.
	Original work(s), source code & author(s) details available at: http://www.darraghcoy.com

	Redistribution and use in source and binary forms, with or without modification, are permitted 
	provided that the following conditions are met:

		* 	Redistributions of source code must retain the above copyright notice, this list of 
			conditions and the following disclaimer.
		* 	Redistributions in binary form must reproduce the above copyright notice, this list of 
			conditions and the following disclaimer in the documentation.
		* 	Neither the name of the author nor the names of contributors may be used to endorse or 
			promote products derived from this software without specific prior written permission.
		* 	The use of any third party libraries, modules or data files provided with the source 
			form of the program must follow and respect the individual licensing terms of these 
			components. You do not have permission to modify or circumvent the licensing terms for 
			provided third party components unless the third party licenses specifically allow this. 
			If you do not agree to the terms of provided third party components then you do not have 
			permission to use them in your derived work(s).
		* 	The terms of this license do not apply to separate (but related) distributions if supplied 
			by the author(s), only to the contents of this source distribution, derived works, and the 
			binaries produced by either modified or unmodified versions of the source code. The terms 
			of this license do not also apply to any third party libraries or components provided with 
			the source distribution. Such third party elements are licensed under their own separate 
			agreements.

	THIS SOFTWARE IS PROVIDED BY DARRAGH COY ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, 
	BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
	ARE DISCLAIMED. IN NO EVENT SHALL DARRAGH COY BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 
	EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 
	OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY 
	OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING 
	IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

//##############################################################################
//##############################################################################
//##############################################################################
//##############################################################################
//
//
// Module Name:  MathLibrary
// Purpose:      C Module consisting of various maths helper functions.
//
//
//##############################################################################
//##############################################################################
//##############################################################################
//##############################################################################

#include "MathLibrary.h"

//------------------------------------------------------------------------------
// Seed for the random number generator:
//------------------------------------------------------------------------------

unsigned int seed;

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlDot2d()
Purpose:  Returns the dot product of two 2D vectors
Params:   two vectors
Returns:  dot product

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

float mlDot2d( const float * vec1 , const float * vec2 )
{
    return vec1[0] * vec2[0] + vec1[1] * vec2[1];
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlAngle2d()
Purpose:  Returns the angle a vector is pointing at (in radians) in range 0-2PI
Params:   vector
Returns:  angle or direction the vector points in

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

float mlAngle2d( const float * vec1 )
{
    //--------------------------------------------------------------------------
    // If both components are zero then return angle of zero by default
    //--------------------------------------------------------------------------

    if ( vec1[0] == 0 && vec1[1] == 0 ) return 0;

    //--------------------------------------------------------------------------
    // Produce a normalised version of the input vector
    //--------------------------------------------------------------------------

    float vecN[2] = { vec1[0] , vec1[1] }; mlNormalise2d(vecN);

    //--------------------------------------------------------------------------
    // Get the inverse sine of the y component of the vector
    //--------------------------------------------------------------------------

    float angle = asinf( vecN[1] );

    //--------------------------------------------------------------------------
    // If x coordinate is negative then angle is PI - angle
    //--------------------------------------------------------------------------

    if ( vec1[0] < 0 ) angle = PI - angle;

    //--------------------------------------------------------------------------
    // If angle is negative add PI*2 onto it to make positive
    //--------------------------------------------------------------------------

    if ( angle < 0 ) angle += PI2;

    //--------------------------------------------------------------------------
    // Finally return angle
    //--------------------------------------------------------------------------

    return angle;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlScale2d()
Purpose:  Scales input vector by given amount.
Params:   vector, amount to scale
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void mlScale2d( float * vec1 , const float scale  )
{
    vec1[0] *= scale;
    vec1[1] *= scale;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlRotate2D()
Purpose:  Rotates vector (anti clockwise) by given amount (in radians)
Params:   vector, amount to rotate by
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void mlRotate2d( float * vec1 , const float angle  )
{
    // Get the coordinates of the vector after rotation

    float rotated_x = vec1[0] * cosf(angle) + vec1[1] * -sinf(angle);
    float rotated_y = vec1[0] * sinf(angle) + vec1[1] *  cosf(angle);

    // Save the newly rotated coordinates

    vec1[0] = rotated_x;
    vec1[1] = rotated_y;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlRotate2DCW()
Purpose:  Rotates vector (clockwise) by given amount (in radians)
Params:   vector, amount to rotate by
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void mlRotate2dCW( float * vec1 , const float angle  )
{
    // Get the coordinates of the vector after rotation

    float rotated_x = vec1[0] *  cosf(angle) + vec1[1] * sinf(angle);
    float rotated_y = vec1[0] * -sinf(angle) + vec1[1] * cosf(angle);

    // Save the newly rotated coordinates

    vec1[0] = rotated_x;
    vec1[1] = rotated_y;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlTranslate2d()
Purpose:  Translates vec1 by vec2
Params:   two vectors
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void mlTranslate2d( float * vec1 , const float * vec2 )
{
    vec1[0] += vec1[0];
    vec1[1] += vec2[1];
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlNormalise2d()
Purpose:  Normalises the vector such that it is at unit lenght
Params:   vector to normalise
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void mlNormalise2d( float * vec1 )
{
    // If zero lenght then make vector point in direction of angle 0

    if ( vec1[0] == 0 && vec1[1] == 0 )
    {
        vec1[0] = 1;
        vec1[1] = 0;

        return;
    }

    // Otherwise get lenght of vector

    float lenght = sqrtf( vec1[0] * vec1[0] + vec1[1] * vec1[1] );

    // Normalise by dividing by lenght:

    vec1[0] /= lenght;
    vec1[1] /= lenght;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlLenght2D()
Purpose:  Returns lenght of a vector
Params:   vector
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

float mlLenght2d( const float * vec1 )
{
    return sqrtf( vec1[0] * vec1[0] + vec1[1] * vec1[1] );
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlMagnitude2d()
Purpose:  Returns magnitude of a vector
Params:   vector
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

float mlMagnitude2d( const float * vec1 )
{
    return vec1[0] * vec1[0] + vec1[1] * vec1[1];
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlPointInPolygon2d()
Purpose:  Tells if a point is inside a given polygon. Polygon must be convex for
          this to work correctly.
Params:   as commented
Returns:  if the point is inside the polygon:

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

bool mlPointInPolygon2d
(
    float       x ,         // x coordinate of point to check
    float       y ,         // y coordinate of point to check
    float *     coords ,    // pointer to array of x,y coordinates for the polygon
    unsigned    numPoints   // number of points in the polygon
)
{
    //--------------------------------------------------------------------------
    // Check if the point is on the inside of every edge:
    //--------------------------------------------------------------------------

    for ( unsigned i = 0 ; i < numPoints ; i++ )
    {
        // Get the edge vector for this edge:

        float edgeX = coords[ ( ( i + 1 ) % numPoints ) * 2 + 0 ] - coords[ i * 2 + 0 ];
        float edgeY = coords[ ( ( i + 1 ) % numPoints ) * 2 + 1 ] - coords[ i * 2 + 1 ];

        // Get the normal to this edge pointing inwards in the polygon

        float normalX = -edgeY;
        float normalY =  edgeX;

        // Get the test point relative to this point in polygon

        float relativeX = x - coords[ i * 2 + 0 ];
        float relativeY = y - coords[ i * 2 + 1 ];

        // Now do a dot product with the edge normal

        float dot;

        dot  = normalX * relativeX;
        dot += normalY * relativeY;

        // If the dot product is less than zero then return false: point is on outside of convex polygon

        if ( dot < 0 ) return false;
    }

    //--------------------------------------------------------------------------
    // If we have got to here then the point is inside: return true
    //--------------------------------------------------------------------------

    return true;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlDot3d()
Purpose:  Dot product between two vectors
Params:   two vectors
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

float mlDot3d( const float * vec1 , const float * vec2 )
{
    return vec1[0] * vec2[0] + vec1[1] * vec2[1] + vec1[2] * vec2[2];
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlScale3d()
Purpose:  Scales input vector by given amount.
Params:   vector, amount to scale
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void mlScale3d( float * vec1 , const float scale  )
{
    vec1[0] *= scale;
    vec1[1] *= scale;
    vec1[2] *= scale;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlTranslate3d()
Purpose:  Translates vec1 by vec2
Params:   two vectors
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void mlTranslate3d( float * vec1 , const float * vec2 )
{
    vec1[0] += vec2[0];
    vec1[1] += vec2[1];
    vec1[2] += vec2[2];
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlNormalise3d()
Purpose:  Normalises the vector such that it is at unit lenght
Params:   vector to normalise
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void mlNormalise3d( float * vec1 )
{
    // If zero lenght then make vector point in direction of angle 0

    if ( vec1[0] == 0 && vec1[1] == 0 && vec1[2] == 0 )
    {
        vec1[0] = 1;
        vec1[1] = 0;
        vec1[2] = 0;

        return;
    }

    // Otherwise get lenght of vector

    float lenght = sqrtf( vec1[0] * vec1[0] + vec1[1] * vec1[1] + vec1[2] * vec1[2] );

    // Normalise by dividing by lenght:

    vec1[0] /= lenght;
    vec1[1] /= lenght;
    vec1[2] /= lenght;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlMagnitude3d()
Purpose:  Returns magnitude of a vector
Params:   vector
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

float mlMagnitude3d( const float * vec1 )
{
    return vec1[0] * vec1[0] + vec1[1] * vec1[1] + vec1[2] * vec1[2];
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlLenght3D()
Purpose:  Returns lenght of a vector
Params:   vector
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

float mlLenght3d( const float * vec1 )
{
    return sqrt( vec1[0] * vec1[0] + vec1[1] * vec1[1] + vec1[2] * vec1[2] );
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlDistance3d()
Purpose:  Returns the distance between two 3D Points
Params:   two points
Returns:  distance between them

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

float mlDistance3d
(
    const float * vec1 ,
    const float * vec2
)
{
    // Get a vector going from vec2 to vec1

    float x = vec1[0] - vec2[0];
    float y = vec1[1] - vec2[1];
    float z = vec1[2] - vec2[2];

    // Now get the lenght of this vector

    float lenght;

    lenght  = x * x;
    lenght += y * y;
    lenght += z * z;

    lenght = sqrtf(lenght);

    // Return the lenght or distance betweent the two points

    return lenght;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlSquaredDistance3d()
Purpose:  Returns the squared distance between two 3D Points
Params:   two points
Returns:  distance between them

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

float mlSquaredDistance3d
(
    const float * vec1 ,
    const float * vec2
)
{
    // Get a vector going from vec2 to vec1

    float x = vec1[0] - vec2[0];
    float y = vec1[1] - vec2[1];
    float z = vec1[2] - vec2[2];

    // Now get the squared lenght of this vector

    float lenght;

    lenght  = x * x;
    lenght += y * y;
    lenght += z * z;

    // Return the squared lenght or squared distance betweent the two points

    return lenght;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlRotate3dX()
Purpose:  Rotates vector by angle degrees about x axis.
Params:   vector , angle to rotate by
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void mlRotate3dX( float * vec1 , const float angle )
{
    float yt = vec1[1] *   cosf(angle) + vec1[2] * - sinf(angle);
    float zt = vec1[1] * + sinf(angle) + vec1[2] * + cosf(angle);

    vec1[1] = yt;
    vec1[2] = zt;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlRotate3dY()
Purpose:  Rotates vector by angle degrees about y axis.
Params:   vector , angle to rotate by
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void mlRotate3dY( float * vec1 , const float angle )
{
    float xt = vec1[0] * cosf(angle) + vec1[2] * - sinf(angle);
    float zt = vec1[0] * sinf(angle) + vec1[2] * + cosf(angle);

    vec1[0] = xt;
    vec1[2] = zt;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlRotate3dZ()
Purpose:  Rotates vector by angle degrees about z axis.
Params:   vector , angle to rotate by
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void mlRotate3dZ( float * vec1 , const float angle )
{
    float xt = vec1[0] *   cosf(angle) + vec1[1] * + sinf(angle);
    float yt = vec1[0] * - sinf(angle) + vec1[1] * + cosf(angle);

    vec1[0] = xt;
    vec1[1] = yt;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlCross3d()
Purpose:  Does the cross product between two vectors and saves the answer to dest
Params:   two vectors, array to save answer to
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void  mlCross3d( const float * vec1 , const float * vec2 , float * dest )
{
    dest[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1];
    dest[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2];
    dest[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0];
}


/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlReflect3d()
Purpose:  Reflects vec1 through vec2
Params:   Vector to reflect. Vector to reflect through.
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void mlReflect3d
(
    float       * vec1,
    const float * vec2
)
{
    // Make a normalised copy of vec2

    float vec2N[3] =
    {
        vec2[0],
        vec2[1],
        vec2[2]
    };

    mlNormalise3d(vec2N);

    // Calculate vec1 dot vec2N by 2

    float dot = mlDot3d(vec1,vec2N) * 2;

    // Now the reflected vector is : ( 2 * dot ) Vec2N - vec1:

    vec1[0] = dot * vec2N[0] - vec1[0];
    vec1[1] = dot * vec2N[1] - vec1[1];
    vec1[2] = dot * vec2N[2] - vec1[2];
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlMatrixMult3x3()
Purpose:  Mutliplies matrix2 by matrix1 and saves result in matrix2
Params:   two matrices two multiply together
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void mlMatrixMult3x3( const float * matrix1, float * matrix2 )
{
    // Create a temporary matrix to save the result to

    float temp[9];

    // Do the multiplication

    temp[0] = matrix1[0] * matrix2[0] + matrix1[1] * matrix2[3] + matrix1[2] * matrix2[6];
    temp[1] = matrix1[0] * matrix2[1] + matrix1[1] * matrix2[4] + matrix1[2] * matrix2[7];
    temp[2] = matrix1[0] * matrix2[2] + matrix1[1] * matrix2[5] + matrix1[2] * matrix2[8];

    temp[3] = matrix1[3] * matrix2[0] + matrix1[4] * matrix2[3] + matrix1[5] * matrix2[6];
    temp[4] = matrix1[3] * matrix2[1] + matrix1[4] * matrix2[4] + matrix1[5] * matrix2[7];
    temp[5] = matrix1[3] * matrix2[2] + matrix1[4] * matrix2[5] + matrix1[5] * matrix2[8];

    temp[6] = matrix1[6] * matrix2[0] + matrix1[7] * matrix2[3] + matrix1[8] * matrix2[6];
    temp[7] = matrix1[6] * matrix2[1] + matrix1[7] * matrix2[4] + matrix1[8] * matrix2[7];
    temp[8] = matrix1[6] * matrix2[2] + matrix1[7] * matrix2[5] + matrix1[8] * matrix2[8];

    // Copy temp into matrix 2

    for ( unsigned i = 0 ; i < 9 ; i++ ) matrix2[i] = temp[i];
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlOrthogonalise3x3()
Purpose:  Forces all three column vectors in the given matrix to be perpendicular
          to each other (orthogonal).
Params:   matrix to orthogonalise
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void mlOrthogonalise3x3( float * matrix )
{
    // Get the original x and z axis vectors

    float xaxis[3] = { matrix[0] , matrix[3] , matrix[6] };
    float zaxis[3] = { matrix[2] , matrix[5] , matrix[8] };

    // Get the y axis vector from doing a cross product between the z and x axes

    float y[3]; mlCross3d( zaxis , xaxis , y );

    // Get the z axis vector from doing a cross product between the x and and new y axis

    float z[3]; mlCross3d( xaxis , y , z );

    // Save new y and z axis vectors

    matrix[1] = y[0];
    matrix[4] = y[1];
    matrix[7] = y[2];

    matrix[2] = z[0];
    matrix[5] = z[1];
    matrix[8] = z[2];
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlOrthogonalise3x3()
Purpose:  Normalises all three column vectors of the 3x3 matrix
Params:   matrix to normalise column vectors in
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void mlNormalise3x3( float * matrix )
{
    // Get the three column vectors of the matrix

    float u[3] = { matrix[0] , matrix[3] , matrix[6] };
    float v[3] = { matrix[1] , matrix[4] , matrix[7] };
    float n[3] = { matrix[2] , matrix[5] , matrix[8] };

    // Normalise all 3 vectors

    mlNormalise3d(u);
    mlNormalise3d(v);
    mlNormalise3d(n);

    // Save the vectors back to the matrix

    matrix[0] = u[0];
    matrix[3] = u[1];
    matrix[6] = u[2];

    matrix[1] = v[0];
    matrix[4] = v[1];
    matrix[7] = v[2];

    matrix[2] = n[0];
    matrix[5] = n[1];
    matrix[8] = n[2];
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlMatrixMult4x4()
Purpose:  Mutliplies matrix2 by matrix1 and saves result in matrix2
Params:   two matrices two multiply together
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void mlMatrixMult4x4( const float * matrix1, float * matrix2 )
{
    // Create a temporary matrix to save the result to

    float temp[16];

    // Do the multiplication

    temp[0 ] = matrix1[0 ] * matrix2[0 ] + matrix1[1 ] * matrix2[4 ] + matrix1[2 ] * matrix2[8 ] + matrix1[3 ] * matrix2[12];
    temp[1 ] = matrix1[0 ] * matrix2[1 ] + matrix1[1 ] * matrix2[5 ] + matrix1[2 ] * matrix2[9 ] + matrix1[3 ] * matrix2[13];
    temp[2 ] = matrix1[0 ] * matrix2[2 ] + matrix1[1 ] * matrix2[6 ] + matrix1[2 ] * matrix2[10] + matrix1[3 ] * matrix2[14];
    temp[3 ] = matrix1[0 ] * matrix2[3 ] + matrix1[1 ] * matrix2[7 ] + matrix1[2 ] * matrix2[11] + matrix1[3 ] * matrix2[15];

    temp[4 ] = matrix1[4 ] * matrix2[0 ] + matrix1[5 ] * matrix2[4 ] + matrix1[6 ] * matrix2[8 ] + matrix1[7 ] * matrix2[12];
    temp[5 ] = matrix1[4 ] * matrix2[1 ] + matrix1[5 ] * matrix2[5 ] + matrix1[6 ] * matrix2[9 ] + matrix1[7 ] * matrix2[13];
    temp[6 ] = matrix1[4 ] * matrix2[2 ] + matrix1[5 ] * matrix2[6 ] + matrix1[6 ] * matrix2[10] + matrix1[7 ] * matrix2[14];
    temp[7 ] = matrix1[4 ] * matrix2[3 ] + matrix1[5 ] * matrix2[7 ] + matrix1[6 ] * matrix2[11] + matrix1[7 ] * matrix2[15];

    temp[8 ] = matrix1[8 ] * matrix2[0 ] + matrix1[9 ] * matrix2[4 ] + matrix1[10] * matrix2[8 ] + matrix1[11] * matrix2[12];
    temp[9 ] = matrix1[8 ] * matrix2[1 ] + matrix1[9 ] * matrix2[5 ] + matrix1[10] * matrix2[9 ] + matrix1[11] * matrix2[13];
    temp[10] = matrix1[8 ] * matrix2[2 ] + matrix1[9 ] * matrix2[6 ] + matrix1[10] * matrix2[10] + matrix1[11] * matrix2[14];
    temp[11] = matrix1[8 ] * matrix2[3 ] + matrix1[9 ] * matrix2[7 ] + matrix1[10] * matrix2[11] + matrix1[11] * matrix2[15];

    temp[12] = matrix1[12] * matrix2[0 ] + matrix1[13] * matrix2[4 ] + matrix1[14] * matrix2[8 ] + matrix1[15] * matrix2[12];
    temp[13] = matrix1[12] * matrix2[1 ] + matrix1[13] * matrix2[5 ] + matrix1[14] * matrix2[9 ] + matrix1[15] * matrix2[13];
    temp[14] = matrix1[12] * matrix2[2 ] + matrix1[13] * matrix2[6 ] + matrix1[14] * matrix2[10] + matrix1[15] * matrix2[14];
    temp[15] = matrix1[12] * matrix2[3 ] + matrix1[13] * matrix2[7 ] + matrix1[14] * matrix2[11] + matrix1[15] * matrix2[15];

    // Copy temp into matrix 2

    for ( unsigned i = 0 ; i < 16 ; i++ ) matrix2[i] = temp[i];
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlTransform3d()
Purpose:  Transforms the 3 compenent vector by the 3x3 transfromation matrix.
Params:   transformation matrix, vector to transform.
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void mlTransform3d( float * vec , const float * matrix )
{
    // Temporary coordinates for the vector being transformed:

    float x;
    float y;
    float z;

    // Do the transforms

    x = vec[0] * matrix[0] + vec[1] * matrix[3] + vec[2] * matrix[6];
    y = vec[0] * matrix[1] + vec[1] * matrix[4] + vec[2] * matrix[7];
    z = vec[0] * matrix[2] + vec[1] * matrix[5] + vec[2] * matrix[8];

    // Save the results of the transforms

    vec[0] = x;
    vec[1] = y;
    vec[2] = z;
}


/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlTransform4d()
Purpose:  Transforms the homogenous vector by the homogenous 4x4 transfromation
          matrix.
Params:   transformation matrix, vector to transform.
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void mlTransform4d( float * vec , const float * matrix )
{
    // Temporary coordinates for the vector being transformed:

    float x;
    float y;
    float z;
    float w;

    // Do the transforms

    x = vec[0] * matrix[0] + vec[1] * matrix[4] + vec[2] * matrix[8 ] + vec[3] * matrix[12];
    y = vec[0] * matrix[1] + vec[1] * matrix[5] + vec[2] * matrix[9 ] + vec[3] * matrix[13];
    z = vec[0] * matrix[2] + vec[1] * matrix[6] + vec[2] * matrix[10] + vec[3] * matrix[14];
    w = vec[0] * matrix[3] + vec[1] * matrix[7] + vec[2] * matrix[11] + vec[3] * matrix[15];

    // Save the results of the transforms

    vec[0] = x;
    vec[1] = y;
    vec[2] = z;
    vec[3] = w;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlTranspose3x3()
Purpose:  Transposes the given 3x3 matrix
Params:   matrix to transpose
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void mlTranspose3x3( float * matrix )
{
    // Produce a temporary matrix with the transposed values

    float transpose[9];

    transpose[0] = matrix[0];
    transpose[1] = matrix[3];
    transpose[2] = matrix[6];

    transpose[3] = matrix[1];
    transpose[4] = matrix[4];
    transpose[5] = matrix[7];

    transpose[6] = matrix[2];
    transpose[7] = matrix[5];
    transpose[8] = matrix[8];

    // Copy the transposed matrix into the given matrix

    for ( unsigned i = 0 ; i < 9 ; i++ ) matrix[i] = transpose[i];
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlTranspose4x4()
Purpose:  Transposes the given 4x4 matrix
Params:   matrix to transpose
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void mlTranspose4x4( float * matrix )
{
    // Produce a temporary matrix with the transposed values

    float transpose[16];

    transpose[0 ] = matrix[0 ];
    transpose[1 ] = matrix[4 ];
    transpose[2 ] = matrix[8 ];
    transpose[3 ] = matrix[12];

    transpose[4 ] = matrix[1 ];
    transpose[5 ] = matrix[5 ];
    transpose[6 ] = matrix[9 ];
    transpose[7 ] = matrix[13];

    transpose[8 ] = matrix[2 ];
    transpose[9 ] = matrix[6 ];
    transpose[10] = matrix[10 ];
    transpose[11] = matrix[14];

    transpose[12] = matrix[3 ];
    transpose[13] = matrix[7 ];
    transpose[14] = matrix[11];
    transpose[15] = matrix[15];

    // Copy the transposed matrix into the given matrix

    for ( unsigned i = 0 ; i < 16 ; i++ ) matrix[i] = transpose[i];
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlPointDistanceToLine3d()
Purpose:  Finds the minimum distance from a point to a given line
Params:   start of line, line direction, point
Returns:  minimum distance between the point and the line

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

float mlPointDistanceToLine3d
(
    const float * point,
    const float * lineStart,
    const float * lineDirection
)
{
    //--------------------------------------------------------------------------
    // If the line direction is zero then return the distance between two points
    // instead...
    //--------------------------------------------------------------------------

    if ( mlMagnitude3d( lineDirection ) == 0 )
    {
        // Make a vector going from the line start to the point

        float vectorTo[3] =
        {
            point[0] - lineStart[0],
            point[1] - lineStart[1],
            point[2] - lineStart[2]
        };

        // Return the lenght of this line. This is the distance between between the two points:

        return mlLenght3d( vectorTo );
    }

    //--------------------------------------------------------------------------
    // Now for some maths theory:
    //
    // Let the point be 'p'
    // Let the line start point be 's'
    // Let the line direction be 'r'
    // Let the line itself be 'L'
    //
    // The equation of any point on the line L is:
    //
    // s + tr
    //
    // where t is a variable
    //
    // The closest point on the line L to p will be also found on a line which
    // is perpendicular to L and passing through p. Let p be a point on a plane
    // and r be the normal vector of the plane. Since r is perpendicular to the
    // mentioned line which pases thorugh p, it is also the normal vector of a
    // the plane which 'p' and the mentioned line are on.
    //
    // Any point q which is on this plane will satisfy the equation:
    //
    // r.q - r.p = 0
    //
    // The four component (xyzw) equation of this plane is hence:
    //
    // x = rx * qx
    // y = ry * qy
    // z = rz * qz
    // w = -r.p
    //
    // and x + y + z + w = 0
    //
    // Now the point we are looking for on line L is also on this plane. We
    // shall call this point 'a'. a also has the same line equation of L.
    //
    // Hence the four components of the plane become:
    //
    // x = rx * ( s + t(rx) )
    // y = ry * ( s + t(ry) )
    // y = rz * ( s + t(rz) )
    // w = -r.p
    //
    // simplifying and solving for t we find:
    //
    // r.p = r.( s + tr )
    // r.p = r.s + t.r.r
    // r.p - r.s = t.r.r
    //
    // t = ( r.p - r.s ) / ( r.r )
    //
    //--------------------------------------------------------------------------

    // Now given this theory find 't':

    float t;

    t  = mlDot3d( lineDirection , point         );
    t -= mlDot3d( lineDirection , lineStart     );
    t /= mlDot3d( lineDirection , lineDirection );

    // Now that we have 't' find the closest point on the line to the point given

    float closest[3] =
    {
        lineStart[0] + lineDirection[0] * t,
        lineStart[1] + lineDirection[1] * t,
        lineStart[2] + lineDirection[2] * t
    };

    // Describe this point relative to the point given, or get a vector going
    // from the point given to this newly figured out point:

    closest[0] -= point[0];
    closest[1] -= point[1];
    closest[2] -= point[2];

    // Now the shortest distance between the point and the line is the lenght
    // of this vector, return this distance:

    return mlLenght3d( closest );
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlPointSquaredDistanceToLine3d()
Purpose:  Finds the minimum squared distance from a point to a given line
Params:   start of line, line direction, point
Returns:  minimum distance between the point and the line

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

float mlPointSquaredDistanceToLine3d
(
    const float * point,
    const float * lineStart,
    const float * lineDirection
)
{
    //--------------------------------------------------------------------------
    // If the line direction is zero then return the magntiude between two points
    // instead...
    //--------------------------------------------------------------------------

    if ( mlMagnitude3d( lineDirection ) == 0 )
    {
        // Make a vector going from the line start to the point

        float vectorTo[3] =
        {
            point[0] - lineStart[0],
            point[1] - lineStart[1],
            point[2] - lineStart[2]
        };

        // Return the magintude of this line. This is the squared distance between between the two points:

        return mlMagnitude3d( vectorTo );
    }

    //--------------------------------------------------------------------------
    // Now for some maths theory:
    //
    // Let the point be 'p'
    // Let the line start point be 's'
    // Let the line direction be 'r'
    // Let the line itself be 'L'
    //
    // The equation of any point on the line L is:
    //
    // s + tr
    //
    // where t is a variable
    //
    // The closest point on the line L to p will be also found on a line which
    // is perpendicular to L and passing through p. Let p be a point on a plane
    // and r be the normal vector of the plane. Since r is perpendicular to the
    // mentioned line which pases thorugh p, it is also the normal vector of a
    // the plane which 'p' and the mentioned line are on.
    //
    // Any point q which is on this plane will satisfy the equation:
    //
    // r.q - r.p = 0
    //
    // The four component (xyzw) equation of this plane is hence:
    //
    // x = rx * qx
    // y = ry * qy
    // z = rz * qz
    // w = -r.p
    //
    // and x + y + z + w = 0
    //
    // Now the point we are looking for on line L is also on this plane. We
    // shall call this point 'a'. a also has the same line equation of L.
    //
    // Hence the four components of the plane become:
    //
    // x = rx * ( s + t(rx) )
    // y = ry * ( s + t(ry) )
    // y = rz * ( s + t(rz) )
    // w = -r.p
    //
    // simplifying and solving for t we find:
    //
    // r.p = r.( s + tr )
    // r.p = r.s + t.r.r
    // r.p - r.s = t.r.r
    //
    // t = ( r.p - r.s ) / ( r.r )
    //
    //--------------------------------------------------------------------------

    // Now given this theory find 't':

    float t;

    t  = mlDot3d( lineDirection , point         );
    t -= mlDot3d( lineDirection , lineStart     );
    t /= mlDot3d( lineDirection , lineDirection );

    // Now that we have 't' find the closest point on the line to the point given

    float closest[3] =
    {
        lineStart[0] + lineDirection[0] * t,
        lineStart[1] + lineDirection[1] * t,
        lineStart[2] + lineDirection[2] * t
    };

    // Describe this point relative to the point given, or get a vector going
    // from the point given to this newly figured out point:

    closest[0] -= point[0];
    closest[1] -= point[1];
    closest[2] -= point[2];

    // Now the shortest squared distance between the point and the line is the
    // magnitude of this vector, return this amount:

    return mlMagnitude3d( closest );
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlLineSphereCollision3d()
Purpose:  Does collision detection between a line segment and a sphere, moving
          the sphere till it does not collide with the line.
Params:   start of line, end of line, sphere position, sphere radius, a vector
          to save the point closest on the line to the center of the bounding sphere.
Returns:  if there was a collision or not.

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

bool mlSphereLineCollision3d
(
    const float *  lineStart,
    const float *  lineEnd,
    float *        spherePosition,
    float          sphereRadius
)
{
    //--------------------------------------------------------------------------
    // HOW TO FIND THE CLOSEST POINT ON A LINE FROM A GIVEN POINT:
    //
    // Let the point given be:          'p'
    // Let the line start point be:     's'
    // Let the line end point be:       'e'
    // Let the line direction be        'e' - 's'
    // Let the line itself be:          'L'
    //
    // The equation of any point on the line L is:
    //
    // s + tr
    //
    // where t is a variable
    //
    // The closest point on the line L to p will be also found on a line which
    // is perpendicular to L and passing through p. Let p be a point on a plane
    // through which this line skims across and r be the normal vector of the plane.
    // Since r is perpendicular to the mentioned perpendicular line which pases
    // thorugh p, it is also the normal vector of a the plane which 'p' and the
    // mentioned line are on.
    //
    // Any point q which is on this plane will satisfy the equation:
    //
    // r.q - r.p = 0
    //
    // where r is the plane normal, p a point on the plane and q any arbitrary point
    //
    // Now the point we are looking for on line L is also on this plane, in
    // fact it is where L intersects the plane. We shall call the point of
    // intersection 'a'. Using the equation of line L which can give the point 'a'
    // and subsituting into the equation of the previously mentioned plane, the
    // following is derived:
    //
    // r.( s + tr ) - r.p = 0
    //
    // simplifying and solving for t we find:
    //
    // r.p  = r.s + tr.r
    // tr.r = r.p - r.s
    // t    = ( r.p - r.s ) / ( r.r )
    //
    // Now using t the closest point on line L can be found to point P
    //
    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    // First find the line direction, r
    //--------------------------------------------------------------------------

    float r[3] =
    {
        lineEnd[0] - lineStart[0],
        lineEnd[1] - lineStart[1],
        lineEnd[2] - lineStart[2]
    };

    //--------------------------------------------------------------------------
    // Find the value for 't' which gives the point along the given line closest
    // to the bounding sphere position. Do this according to the maths theory
    // above.
    //--------------------------------------------------------------------------

    float r_dot_p =
    {
        r[0] * spherePosition[0]
        +
        r[1] * spherePosition[1]
        +
        r[2] * spherePosition[2]
    };

    float r_dot_s =
    {
        r[0] * lineStart[0]
        +
        r[1] * lineStart[1]
        +
        r[2] * lineStart[2]
    };

    float r_dot_r = r[0] * r[0] + r[1] * r[1] + r[2] * r[2];

    float t = ( r_dot_p - r_dot_s ) / ( r_dot_r );

    //--------------------------------------------------------------------------
    // Now that we have t and can find the closest point on the line to the
    // sphere, clamp t to between 0 and 1 so to make sure this point does not
    // over extend the line:
    //--------------------------------------------------------------------------

    t = t < 0 ? 0 : t;
    t = t > 1 ? 1 : t;

    //--------------------------------------------------------------------------
    // Now get the point on the line closest to the bounding sphere center
    //--------------------------------------------------------------------------

    float closestLinePoint[3];

    closestLinePoint[0] = lineStart[0] + r[0] * t;
    closestLinePoint[1] = lineStart[1] + r[1] * t;
    closestLinePoint[2] = lineStart[2] + r[2] * t;

    //--------------------------------------------------------------------------
    // Get the distance between the bounding sphere position and the closest point:
    //--------------------------------------------------------------------------

    float distance = mlDistance3d( closestLinePoint , spherePosition );

    // Save if there was a collision here:

    bool collision = false;

    //--------------------------------------------------------------------------
    // If the distance is less than the radius of the sphere then move the sphere
    // away from this point:
    //--------------------------------------------------------------------------

    if ( distance < sphereRadius )
    {
        //----------------------------------------------------------------------
        // Get a normalised vector going from the point to the sphere position
        //----------------------------------------------------------------------

        float moveVector[3] =
        {
            spherePosition[0] - closestLinePoint[0],
            spherePosition[1] - closestLinePoint[1],
            spherePosition[2] - closestLinePoint[2]
        };

        mlNormalise3d(moveVector);

        //----------------------------------------------------------------------
        // Scale the vector by how far the sphere needs to be moved to be an
        // adequate distance away from the line:
        //----------------------------------------------------------------------

        moveVector[0] *= sphereRadius - distance;
        moveVector[1] *= sphereRadius - distance;
        moveVector[2] *= sphereRadius - distance;

        //----------------------------------------------------------------------
        // Now move the sphere away from the point by this vector
        //----------------------------------------------------------------------

        spherePosition[0] += moveVector[0];
        spherePosition[1] += moveVector[1];
        spherePosition[2] += moveVector[2];

        //----------------------------------------------------------------------
        // Set the collision flag as true
        //----------------------------------------------------------------------

        collision = true;
    }

    // return if there was a collision

    return collision;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlLinePlaneIntersect3d()
Purpose:  Calculates where the given line intersects the plane given and saves
          to given intersection point.
Params:   Names more or less explain them. all paramters are 3 component vectors.
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void mlLinePlaneIntersect3d
(
    const float * lineStart,
    const float * lineDirection,
    const float * planeNormal,
    const float * planePoint,
    float       * intersectPoint
)
{
    //--------------------------------------------------------------------------
    // Math theory:
    //
    // the equation of a point which lies on a plane is given by:
    //
    // n.p = n.p0
    //
    // where n is the plane normal and p a point on the plane, and p0 the point
    // in question which lies on the plane.
    //
    // Now the equation of any point on a ray (or line) is:
    //
    // s + tv
    //
    // where s is the start point and v is the direction of the line and v is
    // a variable.
    //
    // Plugging this equation into the plane equation we get:
    //
    // n.p = n.( s + tv )
    //
    // Rearranging to solve for t we get:
    //
    // n.p = n.s + tn.v
    // tn.v = n.p - n.s
    // t = ( n.p - n.s ) / ( n.v )
    //
    // Using t we can then find the point along the line where it hits the plane
    //--------------------------------------------------------------------------

    // Find n.p , n.s , and n.v

    float n_dot_p = mlDot3d( planeNormal , planePoint     );
    float n_dot_s = mlDot3d( planeNormal , lineStart      );
    float n_dot_v = mlDot3d( planeNormal , lineDirection  );

    // Now if n.v is zero then make it nearly zero such that the interesection
    // point is at near infinity along the line given:

    if ( n_dot_v == 0 ) n_dot_v = 0.0000001f;

    // Now solve for t:

    float t = ( n_dot_p - n_dot_s ) / n_dot_v;

    // Use t to find the point of intersection on the plane:

    intersectPoint[0] = lineStart[0] + lineDirection[0] * t;
    intersectPoint[1] = lineStart[1] + lineDirection[1] * t;
    intersectPoint[2] = lineStart[2] + lineDirection[2] * t;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlPointToPlaneDistance3d()
Purpose:  Gives the distance from a point to a plane. may be negative depending
          on which side of the plane the point lies.
Params:   point to get distance for, plane normal, point on the plane
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

float mlPointToPlaneDistance3d
(
    const float * point,
    const float * planeNormal,
    const float * planePoint
)
{
    // Get the original point relative to the point on the plane

    float x = point[0] - planePoint[0];
    float y = point[1] - planePoint[1];
    float z = point[2] - planePoint[2];

    // Now get the distance to the plane using the dot product and plane normal

    float dist;

    dist  = x * planeNormal[0];
    dist += y * planeNormal[1];
    dist += z * planeNormal[2];

    // Return the the distance

    return dist;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlPointToPlaneSquaredDistance3d()
Purpose:  Gives the squared distance from a point to a plane
Params:   point to get distance for, plane normal, point on the plane
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

float mlPointToPlaneSquaredDistance3d
(
    const float * point,
    const float * planeNormal,
    const float * planePoint
)
{
    // Get the original point relative to the point on the plane

    float x = point[0] - planePoint[0];
    float y = point[1] - planePoint[1];
    float z = point[2] - planePoint[2];

    // Now get the distance to the plane using the dot product and plane normal

    float dist;

    dist  = x * planeNormal[0];
    dist += y * planeNormal[1];
    dist += z * planeNormal[2];

    // Square the distance

    dist *= dist;

    // Return the squared distance

    return dist;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlPointInsidePlane3d()
Purpose:  Tells if the point is on the 'inside' of a given plane
Params:   point, plane normal, point on the plane
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

bool mlPointInsidePlane3d
(
    const float * point,
    const float * planeNormal,
    const float * planePoint
)
{
    // Get the original point relative to the point on the plane

    float x = point[0] - planePoint[0];
    float y = point[1] - planePoint[1];
    float z = point[2] - planePoint[2];

    // Now get the distance to the plane using the dot product and plane normal

    float dist;

    dist  = x * planeNormal[0];
    dist += y * planeNormal[1];
    dist += z * planeNormal[2];

    // If the distance is positive or zero then point is on inside of plane

    return dist >= 0;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlPointInSphere3d()
Purpose:  Tells if a given point is within a sphere
Params:   point, sphere position, sphere radius
Returns:  if the point is in the sphere or not

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

bool mlPointInSphere3d
(
    const float * point,
    const float * spherePosition,
    float         sphereRadius
)
{
    // Square the sphere radius

    sphereRadius *= sphereRadius;

    // Get the vector from the point to sphere position

    float vector[3] =
    {
        point[0] - spherePosition[0],
        point[1] - spherePosition[1],
        point[2] - spherePosition[2]
    };

    // Get the square lenght of this vector

    float lenght =
    {
        vector[0] * vector[0]
        +
        vector[1] * vector[1]
        +
        vector[2] * vector[2]
    };

    // If the square lenght is less than the squared radius then point is inside

    return lenght < sphereRadius;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlSphereInHalfSpace3d()
Purpose:  Tells if a bounding sphere lies within or at least interescts a half
          space region formed by the given plane
Params:   sphere position, sphere radius, plane normal, plane point
Returns:  if the sphere interesects or lies inside the given half space

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

bool mlSphereInsidePlane3d
(
    const float * spherePosition,
    const float   sphereRadius,
    const float * planeNormal,
    const float * planePoint
)
{
    //--------------------------------------------------------------------------
    // Get the sphere position relative to the plane first
    //--------------------------------------------------------------------------

    float pos[3] =
    {
        spherePosition[0] - planePoint[0],
        spherePosition[1] - planePoint[1],
        spherePosition[2] - planePoint[2]
    };

    //--------------------------------------------------------------------------
    // Do a dot product with the plane normal
    //--------------------------------------------------------------------------

    float dot =
    {
        planeNormal[0] * pos[0]
        +
        planeNormal[1] * pos[1]
        +
        planeNormal[2] * pos[2]
    };

    //--------------------------------------------------------------------------
    // Now if the dot product is greater than minus the sphere radius then the
    // sphere lies in or intersects the half space formed by the plane:
    //--------------------------------------------------------------------------

    return dot > -sphereRadius;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlRandomFloat()
Purpose:  Gives back a random float between 0 and 1
Params:   none
Returns:  the randomnly picked float

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

float mlRandomFloat()
{
    // Return a random float

    return float( mlRandomInt() & 0x00ffffff ) / float(0x00ffffff);
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlRandomInt()
Purpose:  Gives back a random integer. Integer may be any value.
Params:   none
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

unsigned mlRandomInt()
{
    // Random genration constants: all of these are nice prime numbers which are
    // very good for giving the appearance of randomness:

    #define ADD 275604541.11f
    #define MUL 32452843.11f

    // Mess around with the seed with integers first:

    seed = unsigned(seed * MUL + ADD);

    // Do a little farting abouting using floats and the infinite receprical, PI

    float d = seed * MUL + ADD + seed * PI_INV;

    // Reinterpret the float as an int: dont' ask why, there are no hard or fast
    // rules to how you get 'random' numbers:

    seed = *((int*)(&d));

    // Return the random number:

    return seed;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlRandomSeed()
Purpose:  Sets the seed for the random number generator.
Params:   seed to set to.
Returns:  none

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

void mlRandomSeed( unsigned randomSeed )
{
    // set the seed

    seed = randomSeed;
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlRound()
Purpose:  Rounds the given float
Params:   float to round
Returns:  rounded version of the float

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

float mlRound( float value )
{
    // See if the float is positive or negative:

    if ( value < 0 )
    {
        // Negative:

        int intValue = (int)(value);

        if ( value - intValue <= -0.5f ) intValue--; return (float) intValue;
    }
    else
    {
        // Positive:

        int intValue = (int)(value);

        if ( value - intValue >= 0.5f ) intValue++; return (float) intValue;
    }
}

/*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

Name:     mlDecimalRound()
Purpose:  Rounds the given float to the given number of decimal places
Params:   float to round
Returns:  rounded version of the float

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

float mlDecimalRound( float value , unsigned decimalPlaces )
{
    float mul = powf( 10.0f , (float) decimalPlaces );

    value *= mul; mlRound(value); value /= mul;

    return value;
}
