/**************************************************************************
 *{@C
 *      Copyright:      1988-2020 Paul Obermeier (obermeier@poSoft.de)
 *
 *      Module:         Utilities
 *      Filename:       UT_VecBasic.c
 *
 *      Author:         Paul Obermeier
 *
 *      Description:	Basic operations on points in 2 and 3 dimensional
 *                      space with Float32 and Float64 precision.
 *			Note:
 *                      The functions working on Float64 values are called
 *                      UT_DVec*.
 *
 *      Additional documentation:
 *                      None.
 *
 *	Exported functions:
 *			UT_VecCopy2D
 *			UT_VecCopy3D
 *			UT_VecSwap2D
 *			UT_VecSwap3D
 *			UT_VecAdd2D
 *			UT_VecAdd3D
 *			UT_VecSub2D
 *			UT_VecSub3D
 *			UT_VecScale2D
 *			UT_VecScale3D
 *			UT_VecCrossProd2D
 *			UT_VecCrossProd3D
 *			UT_VecDotProd2D
 *			UT_VecDotProd3D
 *			UT_VecLength2D
 *			UT_VecLength3D
 *			UT_VecAddScaled2D
 *			UT_VecAddScaled3D
 *			UT_VecLinComb2D
 *			UT_VecLinComb3D
 *			UT_VecReflect2D
 *			UT_VecReflect3D
 *			UT_TfmCopy2D
 *			UT_DTfmCopy2D
 *			UT_TfmCopy3D
 *			UT_DTfmCopy3D
 *			UT_TfmApply3D
 *			UT_TfmDApply3D
 *			UT_DTfmApply3D
 *			UT_DTfmDApply3D
 *			UT_VecDist2D
 *			UT_VecDist3D
 *			UT_VecUnit2D
 *			UT_VecUnit3D
 *			UT_VecRefract2D
 *			UT_VecRefract3D
 *			UT_VecDecompose2D
 *			UT_VecDecompose3D
 *
 **************************************************************************/

#include <stdio.h>
#include <math.h>

#include "UT_Compat.h"
#include "UT_Const.h"
#include "UT_Macros.h"

#include "UT_Vector.h"

/* Two and three dimensional null vectors */

const UT_Vec2D     UT_VecNull2D  =   {0.0, 0.0};
const UT_DVec2D    UT_DVecNull2D =   {0.0, 0.0};

const UT_Vec3D     UT_VecNull3D  =   {0.0, 0.0, 0.0};
const UT_DVec3D    UT_DVecNull3D =   {0.0, 0.0, 0.0};


/* Two and three dimensional vectors containing all ones */

const UT_Vec2D     UT_VecOne2D  =    {1.0, 1.0};
const UT_DVec2D    UT_DVecOne2D =    {1.0, 1.0};

const UT_Vec3D     UT_VecOne3D  =    {1.0, 1.0, 1.0};
const UT_DVec3D    UT_DVecOne3D =    {1.0, 1.0, 1.0};

/***************************************************************************
 *[@e
 *	Name:		UT_VecCopy2D
 *
 *	Usage:		Copy a 2D vector.
 *
 *	Synopsis:	void UT_VecCopy2D (const UT_Vec2D a, UT_Vec2D b)
 *
 *	Description:	Copy a 2D vector: b = a
 *
 *			There is an identical function "UT_DVecCopy2D" for
 *			dealing with vectors in double precision. 
 *
 *	Return value:	None.
 *
 *	See also:
 *
 ***************************************************************************/

void UT_VecCopy2D (const UT_Vec2D a, UT_Vec2D b)
{
    b[0] = a[0];
    b[1] = a[1];
    return;
}

void UT_DVecCopy2D (const UT_DVec2D a, UT_DVec2D b)
{
    b[0] = a[0];
    b[1] = a[1];
    return;
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecCopy3D
 *
 *	Usage:		Copy a 3D vector.
 *
 *	Synopsis:	void UT_VecCopy3D (const UT_Vec3D a, UT_Vec3D b)
 *
 *	Description:	Copy a 3D vector: b = a
 *
 *			There is an identical function "UT_DVecCopy3D" for
 *			dealing with vectors in double precision. 
 *
 *	Return value:	None.
 *
 *	See also:
 *
 ***************************************************************************/

void UT_VecCopy3D (const UT_Vec3D a, UT_Vec3D b)
{
    b[0] = a[0];
    b[1] = a[1];
    b[2] = a[2];
    return;
}

void UT_DVecCopy3D (const UT_DVec3D a, UT_DVec3D b)
{
    b[0] = a[0];
    b[1] = a[1];
    b[2] = a[2];
    return;
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecSwap2D
 *
 *	Usage:		Exchange two 2D vectors.
 *
 *	Synopsis:	void UT_VecSwap2D (UT_Vec2D a, UT_Vec2D b)
 *
 *	Description:	Exchange two 2D vectors: (a, b) = (b, a)
 *
 *			There is an identical function "UT_DVecSwap2D" for
 *			dealing with vectors in double precision. 
 *
 *	Return value:	None.
 *
 *	See also:
 *
 ***************************************************************************/

void UT_VecSwap2D (UT_Vec2D a, UT_Vec2D b)
{
    Float32 tmp;

    tmp  = a[0];
    a[0] = b[0];
    b[0] = tmp;
    tmp  = a[1];
    a[1] = b[1];
    b[1] = tmp;
    return;
}

void UT_DVecSwap2D (UT_DVec2D a, UT_DVec2D b)
{
    Float64 tmp;

    tmp  = a[0];
    a[0] = b[0];
    b[0] = tmp;
    tmp  = a[1];
    a[1] = b[1];
    b[1] = tmp;
    return;
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecSwap3D
 *
 *	Usage:		Exchange two 3D vectors.
 *
 *	Synopsis:	void UT_VecSwap3D (UT_Vec3D a, UT_Vec3D b)
 *
 *	Description:	Exchange two 3D vectors: (a, b) = (b, a)
 *
 *			There is an identical function "UT_DVecSwap3D" for
 *			dealing with vectors in double precision.
 *
 *	Return value:	None.
 *
 *	See also:
 *
 ***************************************************************************/

void UT_VecSwap3D (UT_Vec3D a, UT_Vec3D b)
{
    Float32 tmp;

    tmp  = a[0];
    a[0] = b[0];
    b[0] = tmp;
    tmp  = a[1];
    a[1] = b[1];
    b[1] = tmp;
    tmp  = a[2];
    a[2] = b[2];
    b[2] = tmp;
    return;
}

void UT_DVecSwap3D (UT_DVec3D a, UT_DVec3D b)
{
    Float64 tmp;

    tmp  = a[0];
    a[0] = b[0];
    b[0] = tmp;
    tmp  = a[1];
    a[1] = b[1];
    b[1] = tmp;
    tmp  = a[2];
    a[2] = b[2];
    b[2] = tmp;
    return;
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecAdd2D
 *
 *	Usage:		Add two 2D vectors.
 *
 *	Synopsis:	void UT_VecAdd2D (const UT_Vec2D a, const UT_Vec2D b, UT_Vec2D r)
 *
 *	Description:	Add two 2D vectors: r = a + b
 *
 *			A, b and r may refer to the same memory location.
 *
 *			There is an identical function "UT_DVecAdd2D" for
 *			dealing with vectors in double precision.
 *
 *	Return value:	None.
 *
 *	See also:
 *
 ***************************************************************************/

void UT_VecAdd2D (const UT_Vec2D a, const UT_Vec2D b, UT_Vec2D r)
{
    r[0] = a[0] + b[0];
    r[1] = a[1] + b[1];
    return;
}

void UT_DVecAdd2D (const UT_DVec2D a, const UT_DVec2D b, UT_DVec2D r)
{
    r[0] = a[0] + b[0];
    r[1] = a[1] + b[1];
    return;
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecAdd3D
 *
 *	Usage:		Add two 3D vectors.
 *
 *	Synopsis:	void UT_VecAdd3D (const UT_Vec3D a, const UT_Vec3D b, UT_Vec3D r)
 *
 *	Description:	Add two 3D vectors: r = a + b
 *
 *			A, b and r may refer to the same memory location.
 *
 *			There is an identical function "UT_DVecAdd3D" for
 *			dealing with vectors in double precision.
 *
 *	Return value:	None.
 *
 *	See also:
 *
 ***************************************************************************/

void UT_VecAdd3D (const UT_Vec3D a, const UT_Vec3D b, UT_Vec3D r)
{
    r[0] = a[0] + b[0];
    r[1] = a[1] + b[1];
    r[2] = a[2] + b[2];
    return;
}

void UT_DVecAdd3D (const UT_DVec3D a, const UT_DVec3D b, UT_DVec3D r)
{
    r[0] = a[0] + b[0];
    r[1] = a[1] + b[1];
    r[2] = a[2] + b[2];
    return;
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecSub2D
 *
 *	Usage:		Subtract two 2D vectors.
 *
 *	Synopsis:	void UT_VecSub2D (const UT_Vec2D a, const UT_Vec2D b, UT_Vec2D r)
 *
 *	Description:	Subtract two 2D vectors: r = a - b
 *
 *			A, b and r may refer to the same memory location.
 *
 *			There is an identical function "UT_DVecSub2D" for
 *			dealing with vectors in double precision.
 *
 *	Return value:	None.
 *
 *	See also:
 *
 ***************************************************************************/

void UT_VecSub2D (const UT_Vec2D a, const UT_Vec2D b, UT_Vec2D r)
{
    r[0] = a[0] - b[0];
    r[1] = a[1] - b[1];
    return;
}

void UT_DVecSub2D (const UT_DVec2D a, const UT_DVec2D b, UT_DVec2D r)
{
    r[0] = a[0] - b[0];
    r[1] = a[1] - b[1];
    return;
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecSub3D
 *
 *	Usage:		Subtract two 3D vectors.
 *
 *	Synopsis:	void UT_VecSub3D (const UT_Vec3D a, const UT_Vec3D b, UT_Vec3D r)
 *
 *	Description:	Subtract two 3D vectors: r = a - b
 *
 *			A, b and r may refer to the same memory location.
 *
 *			There is an identical function "UT_DVecSub3D" for
 *			dealing with vectors in double precision.
 *
 *	Return value:	None.
 *
 *	See also:
 *
 ***************************************************************************/

void UT_VecSub3D (const UT_Vec3D a, const UT_Vec3D b, UT_Vec3D r)
{
    r[0] = a[0] - b[0];
    r[1] = a[1] - b[1];
    r[2] = a[2] - b[2];
    return;
}

void UT_DVecSub3D (const UT_DVec3D a, const UT_DVec3D b, UT_DVec3D r)
{
    r[0] = a[0] - b[0];
    r[1] = a[1] - b[1];
    r[2] = a[2] - b[2];
    return;
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecScale2D
 *
 *	Usage:		Multiply a 2D vector and a real number.
 *
 *	Synopsis:	void UT_VecScale2D (Float32 m, const UT_Vec2D a, UT_Vec2D r)
 *
 *	Description:	Multiply a 2D vector and a real number: r = a*m
 *
 *			There is an identical function "UT_DVecScale2D" for
 *			dealing with vectors in double precision.
 *
 *	Return value:	None.
 *
 *	See also:
 *
 ***************************************************************************/

void UT_VecScale2D (Float32 m, const UT_Vec2D a, UT_Vec2D r)
{
    r[0] = a[0] * m;
    r[1] = a[1] * m;
    return;
}

void UT_DVecScale2D (Float64 m, const UT_DVec2D a, UT_DVec2D r)
{
    r[0] = a[0] * m;
    r[1] = a[1] * m;
    return;
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecScale3D
 *
 *	Usage:		Multiply a 3D vector and a real number.
 *
 *	Synopsis:	void UT_VecScale3D (Float32 m, const UT_Vec3D a, UT_Vec3D r)
 *
 *	Description:	Multiply a 3D vector and a real number: r = a*m
 *
 *			There is an identical function "UT_DVecScale3D" for
 *			dealing with vectors in double precision.
 *
 *	Return value:	None.
 *
 *	See also:
 *
 ***************************************************************************/

void UT_VecScale3D (Float32 m, const UT_Vec3D a, UT_Vec3D r)
{
    r[0] = a[0] * m;
    r[1] = a[1] * m;
    r[2] = a[2] * m;
    return;
}

void UT_DVecScale3D (Float64 m, const UT_DVec3D a, UT_DVec3D r)
{
    r[0] = a[0] * m;
    r[1] = a[1] * m;
    r[2] = a[2] * m;
    return;
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecCrossProd2D
 *
 *	Usage:		Calculate the vector product of two 2D vectors.
 *
 *	Synopsis:	Float32 UT_VecCrossProd2D
 *				(const UT_Vec2D a, const UT_Vec2D b)
 *
 *	Description:	Calculate the vector product of two 2D vectors:
 *
 *			{0, 0, r[2]} = {a[0], a[1], 0} x {b[0], b[1], 0}.
 *
 *			There is an identical function "UT_DVecCrossProd2D" for
 *			dealing with vectors in double precision.
 *
 *	Return value:	r[2].
 *
 *	See also:
 *
 ***************************************************************************/

Float32 UT_VecCrossProd2D (const UT_Vec2D a, const UT_Vec2D b)
{
     return a[0] * b[1] - b[0] * a[1];
}

Float64 UT_DVecCrossProd2D (const UT_DVec2D a, const UT_DVec2D b)
{
     return a[0] * b[1] - b[0] * a[1];
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecCrossProd3D
 *
 *	Usage:		Calculate the vector product of two 3D vectors.
 *
 *	Synopsis: 	void UT_VecCrossProd3D
 *				(const UT_Vec3D a, const UT_Vec3D b, UT_Vec3D r)
 *
 *	Description:	Calculate the vector product of two 3D vectors:
 *
 *			R must not refer to the same memory location
 *			as a or b.
 *
 *			There is an identical function "UT_DVecCrossProd3D" for
 *			dealing with vectors in double precision.
 *
 *	Return value:	None.
 *
 *	See also:
 *
 ***************************************************************************/

void UT_VecCrossProd3D (const UT_Vec3D a, const UT_Vec3D b, UT_Vec3D r)
{
    r[0] = a[1] * b[2] - b[1] * a[2];
    r[1] = a[2] * b[0] - b[2] * a[0];
    r[2] = a[0] * b[1] - b[0] * a[1];
    return;
}

void UT_DVecCrossProd3D (const UT_DVec3D a, const UT_DVec3D b, UT_DVec3D r)
{
    r[0] = a[1] * b[2] - b[1] * a[2];
    r[1] = a[2] * b[0] - b[2] * a[0];
    r[2] = a[0] * b[1] - b[0] * a[1];
    return;
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecDotProd2D
 *
 *	Usage:		Calculate the dot product of two 2D vectors.
 *
 *	Synopsis:	Float32 UT_VecDotProd2D (const UT_Vec2D a, const UT_Vec2D b)
 *
 *	Description:	Calculate the dot product of two 2D vectors.
 *
 *			There is an identical function "UT_DVecDotProd2D" for
 *			dealing with vectors in double precision.
 *
 *	Return value:	a . b
 *
 *	See also:
 *
 ***************************************************************************/

Float32 UT_VecDotProd2D (const UT_Vec2D a, const UT_Vec2D b)
{
     return a[0] * b[0] + a[1] * b[1];
}

Float64 UT_DVecDotProd2D (const UT_DVec2D a, const UT_DVec2D b)
{
     return a[0] * b[0] + a[1] * b[1];
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecDotProd3D
 *
 *	Usage:		Calculate the dot product of two 3D vectors.
 *
 *	Synopsis:	Float32 UT_VecDotProd3D (const UT_Vec3D a, const UT_Vec3D b)
 *
 *	Description:	Calculate the dot product of two 3D vectors.
 *
 *			There is an identical function "UT_DVecDotProd3D" for
 *			dealing with vectors in double precision.
 *
 *	Return value:	a . b
 *
 *	See also:
 *
 ***************************************************************************/

Float32 UT_VecDotProd3D (const UT_Vec3D a, const UT_Vec3D b)
{
     return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}

Float64 UT_DVecDotProd3D (const UT_DVec3D a, const UT_DVec3D b)
{
     return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecLength2D
 *
 *	Usage:		Calculate the length of a 2D vector.
 *
 *	Synopsis:	Float32 UT_VecLength2D (const UT_Vec2D a)
 *
 *	Description:	Calculate the length of a 2D vector.
 *
 *			There is an identical function "UT_DVecLength2D" for
 *			dealing with vectors in double precision. 
 *
 *	Return value:	length (a)
 *
 *	See also:
 *
 ***************************************************************************/

Float32 UT_VecLength2D (const UT_Vec2D a)
{
     return sqrt (a[0] * a[0] + a[1] * a[1]);
}

Float64 UT_DVecLength2D (const UT_DVec2D a)
{
     return sqrt (a[0] * a[0] + a[1] * a[1]);
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecLength3D
 *
 *	Usage:		Calculate the length of a 3D vector.
 *
 *	Synopsis:	Float32 UT_VecLength3D (const UT_Vec3D a)
 *
 *	Description:	Calculate the length of a 3D vector.
 *
 *			There is an identical function "UT_DVecLength3D" for
 *			dealing with vectors in double precision. 
 *
 *	Return value:	length (a)
 *
 *	See also:
 *
 ***************************************************************************/

Float32 UT_VecLength3D (const UT_Vec3D a)
{
     return sqrt (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]);
}

Float64 UT_DVecLength3D (const UT_DVec3D a)
{
     return sqrt (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]);
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecAddScaled2D
 *
 *	Usage:		Add two 2D vectors after scaling one of them.
 *
 *	Synopsis:	void UT_VecAddScaled2D
 *				(const UT_Vec2D v1, const UT_Vec2D v2, Float32 a,
 *				 UT_Vec2D r)
 *
 *	Description:	Add two 2D vectors after scaling one of them:
 *			r = a*v1 + v2
 *
 *			There is an identical function "UT_DVecAddScaled2D" for
 *			dealing with vectors in double precision. 
 *
 *	Return value:	None.
 *
 *	See also:
 *
 ***************************************************************************/

void UT_VecAddScaled2D (const UT_Vec2D v1, const UT_Vec2D v2, Float32 a, UT_Vec2D r)
{
    r[0] = a * v1[0] + v2[0];
    r[1] = a * v1[1] + v2[1];
    return;
}

void UT_DVecAddScaled2D (const UT_DVec2D v1, const UT_DVec2D v2, Float64 a, UT_DVec2D r)
{
    r[0] = a * v1[0] + v2[0];
    r[1] = a * v1[1] + v2[1];
    return;
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecAddScaled3D
 *
 *	Usage:		Add two 3D vectors after scaling one of them.
 *
 *	Synopsis:	void UT_VecAddScaled3D
 *				(const UT_Vec3D v1, const UT_Vec3D v2, Float32 a,
 *				 UT_Vec3D r)
 *
 *	Description:	Add two 3D vectors after scaling one of them:
 *			r = a*v1 + v2
 *
 *			There is an identical function "UT_DVecAddScaled3D" for
 *			dealing with vectors in double precision. 
 *
 *	Return value:	None.
 *
 *	See also:
 *
 ***************************************************************************/

void UT_VecAddScaled3D (const UT_Vec3D v1, const UT_Vec3D v2, Float32 a, UT_Vec3D r)
{
    r[0] = a * v1[0] + v2[0];
    r[1] = a * v1[1] + v2[1];
    r[2] = a * v1[2] + v2[2];
    return;
}

void UT_DVecAddScaled3D (const UT_DVec3D v1, const UT_DVec3D v2, Float64 a, UT_DVec3D r)
{
    r[0] = a * v1[0] + v2[0];
    r[1] = a * v1[1] + v2[1];
    r[2] = a * v1[2] + v2[2];
    return;
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecLinComb2D
 *
 *	Usage:		Compute a linear combination of two 2D vectors.
 *
 *	Synopsis:	void UT_VecLinComb2D
 *				(const UT_Vec2D v1, const UT_Vec2D v2,
 *				Float32 a, Float32 b, UT_Vec2D r)
 *
 *	Description:	Compute a linear combination of two 2D vectors:
 *			r = a*v1 + b*v2
 *
 *			There is an identical function "UT_DVecLinComb2D" for
 *			dealing with vectors in double precision. 
 *
 *	Return value:	None.
 *
 *	See also:
 *
 ***************************************************************************/

void UT_VecLinComb2D (const UT_Vec2D v1, const UT_Vec2D v2, Float32 a, Float32 b, UT_Vec2D r)
{
    r[0] = a * v1[0] + b * v2[0];
    r[1] = a * v1[1] + b * v2[1];
    return;
}

void UT_DVecLinComb2D
	(const UT_DVec2D v1, const UT_DVec2D v2,
	 Float64 a, Float64 b, UT_DVec2D r)
{
    r[0] = a * v1[0] + b * v2[0];
    r[1] = a * v1[1] + b * v2[1];
    return;
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecLinComb3D
 *
 *	Usage:		Compute a linear combination of two 3D vectors.
 *
 *	Synopsis: 	void UT_VecLinComb3D
 *				(const UT_Vec3D v1, const UT_Vec3D v2,
 *				 Float32 a, Float32 b, UT_Vec3D r)
 *
 *	Description:	Compute a linear combination of two 3D vectors:
 *			r = a*v1 + b*v2
 *
 *			There is an identical function "UT_DVecLinComb3D" for
 *			dealing with vectors in double precision. 
 *
 *	Return value:	None.
 *
 *	See also:
 *
 ***************************************************************************/

void UT_VecLinComb3D (const UT_Vec3D v1, const UT_Vec3D v2, Float32 a, Float32 b, UT_Vec3D r)
{
    r[0] = a * v1[0] + b * v2[0];
    r[1] = a * v1[1] + b * v2[1];
    r[2] = a * v1[2] + b * v2[2];
    return;
}

void UT_DVecLinComb3D
	(const UT_DVec3D v1, const UT_DVec3D v2,
	 Float64 a, Float64 b, UT_DVec3D r)
{
    r[0] = a * v1[0] + b * v2[0];
    r[1] = a * v1[1] + b * v2[1];
    r[2] = a * v1[2] + b * v2[2];
    return;
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecReflect2D
 *
 *	Usage:		Calculate a reflected 2D vector.
 *
 *	Synopsis:	void UT_VecReflect2D
 *				(const UT_Vec2D v, const UT_Vec2D n, UT_Vec2D r)
 *
 *	Description:	Calculates the direction r of a light ray v that
 *			is reflected on a plane (rsp. straight line) with
 *			vector n as its normal (n must have length 1.0).
 *
 *			There is an identical function "UT_DVecReflect2D" for
 *			dealing with vectors in double precision. 
 *
 *	Return value:	None.
 *
 *	See also:
 *
 ***************************************************************************/

void UT_VecReflect2D (const UT_Vec2D v, const UT_Vec2D n, UT_Vec2D r)
{
    Float32 t;

    t = 2.0 * (n[0] * v[0] + n[1] * v[1]);
    r[0] = v[0] - t * n[0];
    r[1] = v[1] - t * n[1];
    return;
}

void UT_DVecReflect2D (const UT_DVec2D v, const UT_DVec2D n, UT_DVec2D r)
{
    Float64 t;

    t = 2.0 * (n[0] * v[0] + n[1] * v[1]);
    r[0] = v[0] - t * n[0];
    r[1] = v[1] - t * n[1];
    return;
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecReflect3D
 *
 *	Usage:		Calculate a reflected 3D vector.
 *
 *	Synopsis:	void UT_VecReflect3D
 *				(const UT_Vec3D v, const UT_Vec3D n, UT_Vec3D r)
 *
 *	Description:	Calculates the direction r of a light ray v that
 *			is reflected on a plane (rsp. straight line) with
 *			vector n as its normal (n must have length 1.0).
 *
 *			There is an identical function "UT_DVecReflect3D" for
 *			dealing with vectors in double precision. 
 *
 *	Return value:	None.
 *
 *	See also:
 *
 ***************************************************************************/

void UT_VecReflect3D (const UT_Vec3D v, const UT_Vec3D n, UT_Vec3D r)
{
    Float32 t;

    t = 2.0 * (n[0] * v[0] + n[1] * v[1] + n[2] * v[2]);
    r[0] = v[0] - t * n[0];
    r[1] = v[1] - t * n[1];
    r[2] = v[2] - t * n[2];
    return;
}

void UT_DVecReflect3D (const UT_DVec3D v, const UT_DVec3D n, UT_DVec3D r)
{
    Float64 t;

    t = 2.0 * (n[0] * v[0] + n[1] * v[1] + n[2] * v[2]);
    r[0] = v[0] - t * n[0];
    r[1] = v[1] - t * n[1];
    r[2] = v[2] - t * n[2];
    return;
}

/***************************************************************************
 *[@e
 *      Name:           UT_TfmCopy2D
 *
 *      Usage:          Copy a 2D transformation.
 *
 *      Synopsis:       void UT_TfmCopy2D (const UT_Tfm2D S, UT_Tfm2D T)
 *
 *      Description:    Copy a 2D transformation: T = S
 *
 *      Return value:   None.
 *
 *      See also:
 *
 ***************************************************************************/

void UT_TfmCopy2D (const UT_Tfm2D S, UT_Tfm2D T)
{
    T[0][0] = S[0][0];
    T[0][1] = S[0][1];
    T[1][0] = S[1][0];
    T[1][1] = S[1][1];
    T[2][0] = S[2][0];
    T[2][1] = S[2][1];
    return;
}

void UT_DTfmCopy2D (const UT_DTfm2D S, UT_DTfm2D T)
{
    T[0][0] = S[0][0];
    T[0][1] = S[0][1];
    T[1][0] = S[1][0];
    T[1][1] = S[1][1];
    T[2][0] = S[2][0];
    T[2][1] = S[2][1];
    return;
}

/***************************************************************************
 *[@e
 *      Name:           UT_TfmCopy3D
 *
 *      Usage:          Copy a 3D transformation.
 *
 *      Synopsis:       void UT_TfmCopy3D (const UT_Tfm3D S, UT_Tfm3D T)
 *
 *      Description:    Copy a 3D transformation: T = S
 *
 *      Return value:   None.
 *
 *      See also:
 *
 ***************************************************************************/

void UT_TfmCopy3D (const UT_Tfm3D S, UT_Tfm3D T)
{
    T[0][0] = S[0][0];
    T[0][1] = S[0][1];
    T[0][2] = S[0][2];
    T[1][0] = S[1][0];
    T[1][1] = S[1][1];
    T[1][2] = S[1][2];
    T[2][0] = S[2][0];
    T[2][1] = S[2][1];
    T[2][2] = S[2][2];
    T[3][0] = S[3][0];
    T[3][1] = S[3][1];
    T[3][2] = S[3][2];
    return;
}

/***************************************************************************
 *[@e
 *      Name:           UT_DTfmCopy3D
 *
 *      Usage:          Copy a 3D transformation.
 *
 *      Synopsis:       void UT_DTfmCopy3D (const UT_DTfm3D S, UT_DTfm3D T)
 *
 *      Description:    Copy a 3D transformation: T = S
 *
 *      Return value:   None.
 *
 *      See also:
 *
 ***************************************************************************/

void UT_DTfmCopy3D (const UT_DTfm3D S, UT_DTfm3D T)
{
    T[0][0] = S[0][0];
    T[0][1] = S[0][1];
    T[0][2] = S[0][2];
    T[1][0] = S[1][0];
    T[1][1] = S[1][1];
    T[1][2] = S[1][2];
    T[2][0] = S[2][0];
    T[2][1] = S[2][1];
    T[2][2] = S[2][2];
    T[3][0] = S[3][0];
    T[3][1] = S[3][1];
    T[3][2] = S[3][2];
    return;
}

/***************************************************************************
 *[@e
 *      Name:           UT_TfmApply3D
 *
 *      Usage:          Apply a 3D transformation matrix to a vector.
 *
 *      Synopsis:       void UT_TfmApply3D
 *                              (const UT_Vec3D v, const UT_Tfm3D T, UT_Vec3D r)
 *
 *      Description:    The point represented by v is moved in space
 *                      according to transformation T (r = vT).
 *                      R must not refer to the same memory location as v.
 *
 *                      There is an identical function "UT_DTfmApply3D" for
 *                      dealing with vectors in double precision.
 *
 *      Return value:   None.
 *
 *      See also:
 *
 ***************************************************************************/

void UT_TfmApply3D (const UT_Vec3D v, const UT_Tfm3D T, UT_Vec3D r)
{
    r[0] = v[0] * T[0][0] + v[1] * T[1][0] + v[2] * T[2][0] + T[3][0];
    r[1] = v[0] * T[0][1] + v[1] * T[1][1] + v[2] * T[2][1] + T[3][1];
    r[2] = v[0] * T[0][2] + v[1] * T[1][2] + v[2] * T[2][2] + T[3][2];
    return;
}

void UT_DTfmApply3D (const UT_Vec3D v, const UT_DTfm3D T, UT_Vec3D r)
{
    r[0] = v[0] * T[0][0] + v[1] * T[1][0] + v[2] * T[2][0] + T[3][0];
    r[1] = v[0] * T[0][1] + v[1] * T[1][1] + v[2] * T[2][1] + T[3][1];
    r[2] = v[0] * T[0][2] + v[1] * T[1][2] + v[2] * T[2][2] + T[3][2];
    return;
}

void UT_TfmDApply3D (const UT_DVec3D v, const UT_Tfm3D T, UT_DVec3D r)
{
    r[0] = v[0] * T[0][0] + v[1] * T[1][0] + v[2] * T[2][0] + T[3][0];
    r[1] = v[0] * T[0][1] + v[1] * T[1][1] + v[2] * T[2][1] + T[3][1];
    r[2] = v[0] * T[0][2] + v[1] * T[1][2] + v[2] * T[2][2] + T[3][2];
    return;
}

void UT_DTfmDApply3D (const UT_DVec3D v, const UT_DTfm3D T, UT_DVec3D r)
{
    r[0] = v[0] * T[0][0] + v[1] * T[1][0] + v[2] * T[2][0] + T[3][0];
    r[1] = v[0] * T[0][1] + v[1] * T[1][1] + v[2] * T[2][1] + T[3][1];
    r[2] = v[0] * T[0][2] + v[1] * T[1][2] + v[2] * T[2][2] + T[3][2];
    return;
}
/***************************************************************************
 *[@e
 *	Name:		UT_VecDist2D
 *
 *	Usage:		Calculate the distance between two 2D points.
 *
 *	Synopsis:	Float32 UT_VecDist2D (const UT_Vec2D a, const UT_Vec2D b)
 *
 *	Description:	Calculate the distance between two 2D points.
 *
 *			There is an identical function "UT_DVecDist2D" for
 *			dealing with vectors in double precision. 
 *
 *	Return value:	length (a - b).
 *
 *	See also:
 *
 ***************************************************************************/

Float32 UT_VecDist2D (const UT_Vec2D a, const UT_Vec2D b)
{
    Float32 t0, t1;

    t0 = a[0] - b[0];
    t1 = a[1] - b[1];
    return sqrt (t0 * t0 + t1 * t1);
}

Float64 UT_DVecDist2D (const UT_DVec2D a, const UT_DVec2D b)
{
    Float64 t0, t1;

    t0 = a[0] - b[0];
    t1 = a[1] - b[1];
    return sqrt (t0 * t0 + t1 * t1);
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecDist3D
 *
 *	Usage:		Calculate the distance between two 3D points.
 *
 *	Synopsis:	Float32 UT_VecDist3D (const UT_Vec3D a, const UT_Vec3D b)
 *
 *	Description:	Calculate the distance between two 3D points.
 *
 *			There is an identical function "UT_DVecDist3D" for
 *			dealing with vectors in double precision.
 *
 *	Return value:	length (a - b).
 *
 *	See also:
 *
 ***************************************************************************/

Float32 UT_VecDist3D (const UT_Vec3D a, const UT_Vec3D b)
{
    Float32 t0, t1, t2;

    t0 = a[0] - b[0];
    t1 = a[1] - b[1];
    t2 = a[2] - b[2];
    return sqrt (t0 * t0 + t1 * t1 + t2 * t2);
}

Float64 UT_DVecDist3D (const UT_DVec3D a, const UT_DVec3D b)
{
    Float64 t0, t1, t2;

    t0 = a[0] - b[0];
    t1 = a[1] - b[1];
    t2 = a[2] - b[2];
    return sqrt (t0 * t0 + t1 * t1 + t2 * t2);
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecUnit2D
 *
 *	Usage:		Scale a 2D vector so that its length becomes 1.0.
 *
 *	Synopsis:	UT_Bool UT_VecUnit2D (UT_Vec2D v)
 *
 *	Description:	Scale a 2D vector so that its length becomes 1.0.
 *			v = 1 / max (MinNrmFloat32, length(v)) * v
 *
 *			There is an identical function "UT_DVecUnit2D" for
 *			dealing with vectors in double precision. 
 *
 *	Return value:	UT_True if length(v) was greater than MinNrmFloat32
 *			before v was scaled, otherwise UT_False.
 *
 *	See also:
 *
 ***************************************************************************/

UT_Bool UT_VecUnit2D (UT_Vec2D v)
{
    Float32 t;

    /* To avoid numerical difficulties, do not replace the divisions
       by multiplications! */

    if ((t = sqrt (v[0] * v[0] + v[1] * v[1])) < MinNrmFloat32) {
	v[0] /= MinNrmFloat32;
	v[1] /= MinNrmFloat32;
	return UT_False;
    }
    v[0] /= t;
    v[1] /= t;
    return UT_True;
}

UT_Bool UT_DVecUnit2D (UT_DVec2D v)
{
    Float64 t;

    /* To avoid numerical difficulties, do not replace the divisions
       by multiplications! */

    if ((t = sqrt (v[0] * v[0] + v[1] * v[1])) < MinNrmFloat64) {
	v[0] /= MinNrmFloat64;
	v[1] /= MinNrmFloat64;
	return UT_False;
    }
    v[0] /= t;
    v[1] /= t;
    return UT_True;
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecUnit3D
 *
 *	Usage:		Scale a 3D vector so that its length becomes 1.0.
 *
 *	Synopsis:	UT_Bool UT_VecUnit3D (UT_Vec3D v)
 *
 *	Description:	Scale a 3D vector so that its length becomes 1.0.
 *			v = 1 / max (MinNrmFloat32, length(v)) * v
 *
 *			There is an identical function "UT_DVecUnit3D" for
 *			dealing with vectors in double precision. 
 *
 *	Return value:	UT_True if length(v) was greater than MinNrmFloat32
 *			before v was scaled, otherwise UT_False.
 *
 *	See also:
 *
 ***************************************************************************/

UT_Bool UT_VecUnit3D (UT_Vec3D v)
{
    Float32 t;

    /* To avoid numerical difficulties, do not replace the divisions
       by multiplications! */

    if ((t = sqrt (v[0] * v[0] + v[1] * v[1] + v[2] * v[2])) < MinNrmFloat32) {
	v[0] /= MinNrmFloat32;
	v[1] /= MinNrmFloat32;
	v[2] /= MinNrmFloat32;
	return UT_False;
    }
    v[0] /= t;
    v[1] /= t;
    v[2] /= t;
    return UT_True;
}

UT_Bool UT_DVecUnit3D (UT_DVec3D v)
{
    Float64 t;

    /* To avoid numerical difficulties, do not replace the divisions
       by multiplications! */

    if ((t = sqrt (v[0] * v[0] + v[1] * v[1] + v[2] * v[2])) < MinNrmFloat64) {
	v[0] /= MinNrmFloat64;
	v[1] /= MinNrmFloat64;
	v[2] /= MinNrmFloat64;
	return UT_False;
    }
    v[0] /= t;
    v[1] /= t;
    v[2] /= t;
    return UT_True;
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecRefract2D
 *
 *	Usage:		Calculate a refracted 2D vector.
 *
 *	Synopsis:	UT_Bool UT_VecRefract2D
 *				(const UT_Vec2D v, const UT_Vec2D n,
 *				 Float32 i1, Float32 i2, UT_Vec2D r)
 *
 *	Description:	Computes the direction r of a light ray v after
 *			passing through a plane (rsp. straight line) with
 *			normal vector n (n must have length 1.0).
 *			The index of refraction on that side of the plane
 *			to which n points is i1, the index of refraction
 *			on the other side is i2.
 *
 *			There is an identical function "UT_DVecRefract2D" for
 *			dealing with vectors in double precision. 
 *
 *	Return value:	If total reflection of v occurs at the plane,
 *			UT_VecRefract returns UT_False without attempting
 *			to find r.
 *			If total reflection of v at the plane does not occur,
 *			UT_VecRefract computes r, scales it to length 1.0 and
 *			returns UT_True.
 *
 *	See also:
 *
 ***************************************************************************/

UT_Bool UT_VecRefract2D (const UT_Vec2D v, const UT_Vec2D n, Float32 i1, Float32 i2, UT_Vec2D r)
{
    Float32 vn, ss, t, f;
    UT_Vec2D s;

    vn = v[0] * n[0] + v[1] * n[1];
    if (vn > 0.0) {
	f = i2 / i1;
	s[0] = (v[0] - n[0] * vn) * f;
	s[1] = (v[1] - n[1] * vn) * f;
	if (1.0 <= (ss = s[0] * s[0] + s[1] * s[1])) {
	    return UT_False;
        }
	t = sqrt (1.0 - ss);
	r[0] = s[0] + n[0] * t;
	r[1] = s[1] + n[1] * t;
    } else {
	f = i1 / i2;
	s[0] = (v[0] - n[0] * vn) * f;
	s[1] = (v[1] - n[1] * vn) * f;
	if (1.0 <= (ss = s[0] * s[0] + s[1] * s[1])) {
	    return UT_False;
        }
	t = sqrt (1.0 - ss);
	r[0] = s[0] - n[0] * t;
	r[1] = s[1] - n[1] * t;
    }
    return UT_True;
}

UT_Bool UT_DVecRefract2D
	(const UT_DVec2D v, const UT_DVec2D n,
	 Float64 i1, Float64 i2, UT_DVec2D r)
{
    Float64 vn, ss, t, f;
    UT_DVec2D s;

    vn = v[0] * n[0] + v[1] * n[1];
    if (vn > 0.0) {
	f = i2 / i1;
	s[0] = (v[0] - n[0] * vn) * f;
	s[1] = (v[1] - n[1] * vn) * f;
	if (1.0 <= (ss = s[0] * s[0] + s[1] * s[1])) {
	    return UT_False;
        }
	t = sqrt (1.0 - ss);
	r[0] = s[0] + n[0] * t;
	r[1] = s[1] + n[1] * t;
    } else {
	f = i1 / i2;
	s[0] = (v[0] - n[0] * vn) * f;
	s[1] = (v[1] - n[1] * vn) * f;
	if (1.0 <= (ss = s[0] * s[0] + s[1] * s[1])) {
	    return UT_False;
        }
	t = sqrt (1.0 - ss);
	r[0] = s[0] - n[0] * t;
	r[1] = s[1] - n[1] * t;
    }
    return UT_True;
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecRefract3D
 *
 *	Usage:		Calculate a refracted 3D vector.
 *
 *	Synopsis:	UT_Bool UT_VecRefract3D
 *				(const UT_Vec3D v, const UT_Vec3D n,
 *				 Float32 i1, Float32 i2, UT_Vec3D r)
 *
 *	Description:	Computes the direction r of a light ray v after
 *			passing through a plane (rsp. straight line) with
 *			normal vector n (n must have length 1.0).
 *			The index of refraction on that side of the plane
 *			to which n points is i1, the index of refraction
 *			on the other side is i2.
 *
 *			There is an identical function "UT_DVecRefract3D" for
 *			dealing with vectors in double precision. 
 *
 *	Return value:	If total reflection of v occurs at the plane,
 *			UT_VecRefract returns UT_False without attempting
 *			to find r.
 *			If total reflection of v at the plane does not occur,
 *			UT_VecRefract computes r, scales it to length 1.0 and
 *			returns UT_True.
 *
 *	See also:
 *
 ***************************************************************************/

UT_Bool UT_VecRefract3D (const UT_Vec3D v, const UT_Vec3D n, Float32 i1, Float32 i2, UT_Vec3D r)
{
    Float32 vn, ss, t, f;
    UT_Vec3D s;

    vn = v[0] * n[0] + v[1] * n[1] + v[2] * n[2];
    if (vn > 0.0) {
	f = i2 / i1;
	s[0] = (v[0] - n[0] * vn) * f;
	s[1] = (v[1] - n[1] * vn) * f;
	s[2] = (v[2] - n[2] * vn) * f;
	if (1.0 <= (ss = s[0] * s[0] + s[1] * s[1] + s[2] * s[2])) {
	    return UT_False;
        }
	t = sqrt (1.0 - ss);
	r[0] = s[0] + n[0] * t;
	r[1] = s[1] + n[1] * t;
	r[2] = s[2] + n[2] * t;
    } else {
	f = i1 / i2;
	s[0] = (v[0] - n[0] * vn) * f;
	s[1] = (v[1] - n[1] * vn) * f;
	s[2] = (v[2] - n[2] * vn) * f;
	if (1.0 <= (ss = s[0] * s[0] + s[1] * s[1] + s[2] * s[2])) {
	    return UT_False;
        }
	t = sqrt (1.0 - ss);
	r[0] = s[0] - n[0] * t;
	r[1] = s[1] - n[1] * t;
	r[2] = s[2] - n[2] * t;
    }
    return UT_True;
}

UT_Bool UT_DVecRefract3D
	(const UT_DVec3D v, const UT_DVec3D n,
	 Float64 i1, Float64 i2, UT_DVec3D r)
{
    Float64 vn, ss, t, f;
    UT_DVec3D s;

    vn = v[0] * n[0] + v[1] * n[1] + v[2] * n[2];
    if (vn > 0.0) {
	f = i2 / i1;
	s[0] = (v[0] - n[0] * vn) * f;
	s[1] = (v[1] - n[1] * vn) * f;
	s[2] = (v[2] - n[2] * vn) * f;
	if (1.0 <= (ss = s[0] * s[0] + s[1] * s[1] + s[2] * s[2])) {
	    return UT_False;
        }
	t = sqrt (1.0 - ss);
	r[0] = s[0] + n[0] * t;
	r[1] = s[1] + n[1] * t;
	r[2] = s[2] + n[2] * t;
    } else {
	f = i1 / i2;
	s[0] = (v[0] - n[0] * vn) * f;
	s[1] = (v[1] - n[1] * vn) * f;
	s[2] = (v[2] - n[2] * vn) * f;
	if (1.0 <= (ss = s[0] * s[0] + s[1] * s[1] + s[2] * s[2])) {
	    return UT_False;
        }
	t = sqrt (1.0 - ss);
	r[0] = s[0] - n[0] * t;
	r[1] = s[1] - n[1] * t;
	r[2] = s[2] - n[2] * t;
    }
    return UT_True;
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecDecompose2D
 *
 *	Usage:		Decompose a 2D vector.
 *
 *	Synopsis:	void UT_VecDecompose2D
 *				(const UT_Vec2D v, const UT_Vec2D n,
 *				 UT_Vec2D m, UT_Vec2D p)
 *
 *	Description:	Given a 2D vector n of length 1.0, decompose a vector
 *			v into two components m and p, so that v = m+p and
 *			m is a multiple of n and p is perpendicular to n.
 *
 *			If you are not interested in one of the vectors m or
 *			p, specify NULL.
 *
 *			There is an identical function "UT_DVecDecompose2D" for
 *			dealing with vectors in double precision. 
 *
 *	Return value:	None.
 *
 *	See also:
 *
 ***************************************************************************/

void UT_VecDecompose2D (const UT_Vec2D v, const UT_Vec2D n, UT_Vec2D m, UT_Vec2D p)
{
    Float32 t;

    t = n[0] * v[0] + n[1] * v[1];
    if (m) {
	m[0] = n[0] * t;
	m[1] = n[1] * t;
	if (p) {
	    p[0] = v[0] - m[0];
	    p[1] = v[1] - m[1];
	}
	return;
    }
    if (p) {
	p[0] = v[0] - n[0] * t;
	p[1] = v[1] - n[1] * t;
    }
    return;
}

void UT_DVecDecompose2D (const UT_DVec2D v, const UT_DVec2D n, UT_DVec2D m, UT_DVec2D p)
{
    Float64 t;

    t = n[0] * v[0] + n[1] * v[1];
    if (m) {
	m[0] = n[0] * t;
	m[1] = n[1] * t;
	if (p) {
	    p[0] = v[0] - m[0];
	    p[1] = v[1] - m[1];
	}
	return;
    }
    if (p) {
	p[0] = v[0] - n[0] * t;
	p[1] = v[1] - n[1] * t;
    }
    return;
}

/***************************************************************************
 *[@e
 *	Name:		UT_VecDecompose3D
 *
 *	Usage:		Decompose a 3D vector.
 *
 *	Synopsis:	void UT_VecDecompose3D
 *				(const UT_Vec3D v, const UT_Vec3D n,
 *				 UT_Vec3D m, UT_Vec3D p)
 *
 *	Description:	Given a 3D vector n of length 1.0, decompose a vector
 *			v into two components m and p, so that v = m+p and
 *			m is a multiple of n and p is perpendicular to n.
 *
 *			If you are not interested in one of the vectors m or
 *			p, specify NULL.
 *
 *			There is an identical function "UT_DVecDecompose3D" for
 *			dealing with vectors in double precision. 
 *
 *	Return value:	None.
 *
 *	See also:
 *
 ***************************************************************************/

void UT_VecDecompose3D (const UT_Vec3D v, const UT_Vec3D n, UT_Vec3D m, UT_Vec3D p)
{
    Float32 t;

    t = n[0] * v[0] + n[1] * v[1] + n[2] * v[2];
    if (m) {
	m[0] = n[0] * t;
	m[1] = n[1] * t;
	m[2] = n[2] * t;
	if (p) {
	    p[0] = v[0] - m[0];
	    p[1] = v[1] - m[1];
	    p[2] = v[2] - m[2];
	}
	return;
    }
    if (p) {
	p[0] = v[0] - n[0] * t;
	p[1] = v[1] - n[1] * t;
	p[2] = v[2] - n[2] * t;
    }
    return;
}

void UT_DVecDecompose3D (const UT_DVec3D v, const UT_DVec3D n, UT_DVec3D m, UT_DVec3D p)
{
    Float64 t;

    t = n[0] * v[0] + n[1] * v[1] + n[2] * v[2];
    if (m) {
	m[0] = n[0] * t;
	m[1] = n[1] * t;
	m[2] = n[2] * t;
	if (p) {
	    p[0] = v[0] - m[0];
	    p[1] = v[1] - m[1];
	    p[2] = v[2] - m[2];
	}
	return;
    }
    if (p) {
	p[0] = v[0] - n[0] * t;
	p[1] = v[1] - n[1] * t;
	p[2] = v[2] - n[2] * t;
    }
    return;
}
