#ifndef lint
static char *sccsid = "@(#)vector.c	1.1 (Steve Hill) 3/9/90";
#endif

/* vector.c
 *
 * This module provides a selection of vector operations.
 *
 * Operations over two vectors generally take the form
 *
 * op(v1, v2, v3)
 *
 * and perform v1 := v2 op v3.  They also return v1.
 */

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

#include "basetype.h"
#include "cartesian.h"
#include "error.h"
#include "malloc.h"
#include "vector.h"


/* GetVector
 *
 * Private function to allocate a vector structure.
 */

static
vector_t *
GetVector()
{
	vector_t	*vector;

	vector = (vector_t *) malloc(sizeof(vector_t));
	if (vector == VectorNull)
		FatalError(__FILE__, __LINE__, "GetVector: Out of memory");

	return(vector);
}


/* Vector
 *
 * Allocate and initialise a vector to (x, y, z).
 */

vector_t *
Vector(x, y, z)
real_t	x, y, z;
{
	vector_t	*vector;

	vector = GetVector();

	vector->x = x;
	vector->y = y;
	vector->z = z;

	return(vector);
}


/* VectorCopy
 *
 * Return pointer to a duplicate vector to argument.
 */

vector_t *
VectorCopy(vector)
vector_t	*vector;
{
	if (vector == VectorNull)
		return(VectorNull);

	return(Vector(vector->x, vector->y, vector->z));
}


/* VectorUnit
 *
 * Allocate and return pointer to unit vector along specified axis.
 */

vector_t *
VectorUnit(axis)
axis_t	axis;
{
	switch (axis)
	{
	case X_AXIS:
		return(Vector(REAL_ONE, REAL_ZERO, REAL_ZERO));
	case Y_AXIS:
		return(Vector(REAL_ZERO, REAL_ONE, REAL_ZERO));
	case Z_AXIS:
		return(Vector(REAL_ZERO, REAL_ZERO, REAL_ONE));
	default:
		FatalError(__FILE__, __LINE__, "VectorUnit: Bad axis");
	}
	/*NOTREACHED*/
}


/* VectorLet
 *
 * Assign vector2 to vector1.
 */

vector_t *
VectorLet(vector1, vector2)
vector_t	*vector1, *vector2;
{
	vector1->x = vector2->x;
	vector1->y = vector2->y;
	vector1->z = vector2->z;

	return(vector1);
}


/* VectorZero
 *
 * Allocate and return pointer to zero vector.
 */

vector_t *
VectorZero()
{
	return(Vector(REAL_ZERO, REAL_ZERO, REAL_ZERO));
}


/* VectorNeg
 *
 * Negate a vector.  Returns pointer to result vector.
 */

vector_t *
VectorNeg(vector1, vector2)
vector_t	*vector1, *vector2;
{
	if (vector1 == VectorNull)
		vector1 = GetVector();

	vector1->x = -(vector2->x);
	vector1->y = -(vector2->y);
	vector1->z = -(vector2->z);

	return(vector1);
}


/* VectorAdd
 *
 * Add vector vector2 and vector3.  Return result in vector1 (allocate
 * if necessary).
 */

vector_t *
VectorAdd(vector1, vector2, vector3)
vector_t	*vector1, *vector2, *vector3;
{
	if (vector1 == VectorNull)
		vector1 = GetVector();

	vector1->x = vector2->x + vector3->x;
	vector1->y = vector2->y + vector3->y;
	vector1->z = vector2->z + vector3->z;

	return(vector1);
}


/* VectorSub
 *
 * Subtract vector vector3 from vector2.  Return result in vector1 (allocate
 * if necessary).
 */

vector_t *
VectorSub(vector1, vector2, vector3)
vector_t	*vector1, *vector2, *vector3;
{
	if (vector1 == VectorNull)
		vector1 = GetVector();

	vector1->x = vector2->x - vector3->x;
	vector1->y = vector2->y - vector3->y;
	vector1->z = vector2->z - vector3->z;

	return(vector1);
}


/* VectorMul
 *
 * Multiply vector3 and vector2.
 * Return result in v1 (allocate if necessary).
 */

vector_t *
VectorMul(vector1, vector2, vector3)
vector_t	*vector1, *vector2, *vector3;
{
	if (vector1 == VectorNull)
		vector1 = GetVector();

	vector1->x = vector2->x * vector3->x;
	vector1->y = vector2->y * vector3->y;
	vector1->z = vector2->z * vector3->z;

	return(vector1);
}


/* VectorDot
 *
 * Return dot product of two vectors.
 */

real_t
VectorDot(vector1, vector2)
vector_t	*vector1, *vector2;
{
	return((real_t)
	       (vector1->x * vector2->x +
	        vector1->y * vector2->y +
	        vector1->z * vector2->z));
}


/* VectorCross
 *
 * Form cross product of vector2 and vector3 - return in vector1 (allocate
 * if necessary).
 */

vector_t *
VectorCross(vector1, vector2, vector3)
vector_t	*vector1, *vector2, *vector3;
{
	if (vector1 == VectorNull)
		vector1 = GetVector();

	vector1->x = (vector2->y * vector3->z) - (vector2->z * vector3->y);
	vector1->y = (vector2->z * vector3->x) - (vector2->x * vector3->z);
	vector1->z = (vector2->x * vector3->y) - (vector2->y * vector3->x);

	return(vector1);
}


/* VectorScale
 *
 * Multply vector2 by a scalar.  Return in vector1 (allocate if necessary).
 */

vector_t *
VectorScale(vector1, k, vector2)
vector_t	*vector1, *vector2;
real_t		k;
{
	if (vector1 == VectorNull)
		vector1 = GetVector();

	vector1->x = k * vector2->x;
	vector1->y = k * vector2->y;
	vector1->z = k * vector2->z;

	return(vector1);
}


/* VectorNormalise
 *
 * Produce a normalised vector.  A normalised vector has modulus of
 * one.  Return normalised vector in vector1 (allocate if necessary).
 */

vector_t *
VectorNormalise(vector1, vector2)
vector_t	*vector1, *vector2;
{
	real_t	d;

	if (vector1 == VectorNull)
		vector1 = GetVector();

	d = (real_t) sqrt((double) VectorDot(vector2, vector2));

	if (d != REAL_ZERO)
	{
		vector1->x = vector2->x / d;
		vector1->y = vector2->y / d;
		vector1->z = vector2->z / d;
	}

	return(vector1);
}


/* VectorFree
 *
 * Free a vector.
 */

void
VectorFree(vector)
vector_t	*vector;
{
	if (vector == VectorNull)
		return;

	(void) free((char *) vector);
}


/* VectorPrint
 *
 * print a vector to a file.
 */

void
VectorPrint(file, vector)
FILE		*file;
vector_t	*vector;
{
	if (vector == VectorNull)
	{
		fprintf(file, "Null Vector\n");
		return;
	}

	fprintf(file, "(%f, %f, %f)", vector->x, vector->y, vector->z);
}