// Copyright 2005 Keir Mierle.
// Licensed under the BSD license

#ifndef _VEC_H
#define _VEC_H

// Reduce code: We use this all over the place.
#define DimLoopReturn(op) \
	vec<T,dim> ret; \
	for (int i=0; i<dim; i++) \
		ret[i] = op; \
	return ret;

#define DimLoop(op) \
	for (int i=0; i<dim; i++) \
		op;

// Aid us in the gl calls - Makes it easy to use my vector class with OpenGL.
#define expand(v) (v).coords[0], (v).coords[1], (v).coords[2]
#define expandp(v) (v)->coords[0], (v)->coords[1], (v)->coords[2]

#include <iostream>
#include <math.h>

template <class T, int dim>
class vec
{
	public:
	T coords[dim];

	vec<T,dim>() { }
	vec<T,dim>(const vec<T,dim>& b) 
		{ DimLoop(coords[i] = b.coords[i]) }
	vec<T,dim>(const T b)
		{ DimLoop(coords[i] = b) }

	T &operator[]( const int ndx ) { return( coords[ndx] ); }

	const vec<T,dim>& operator= (const vec<T,dim>& param) 
	{
		DimLoop(coords[i] = param.coords[i])
		return *this;
	}

	const vec<T,dim>& operator= (const T a)
	{
		DimLoop(coords[i] = a)
		return *this;
	}

	// Arithmetic operators
	vec<T,dim> operator+ (const vec<T,dim>& param) const
		{ DimLoopReturn(coords[i] + param.coords[i]) }
	void operator+= (const vec<T,dim>& param) 
		{ DimLoop(coords[i] += param.coords[i]) }
	vec<T,dim> operator- (const vec<T,dim>& param) const
		{ DimLoopReturn(coords[i] - param.coords[i]) }
	void operator-= (const vec<T,dim>& param) 
		{ DimLoop(coords[i] -= param.coords[i]) }
	vec<T,dim> operator- ()  const
		{ DimLoopReturn(ret[i] = -coords[i]) }
	vec<T,dim> operator* (const vec<T,dim>& param)  const
		{ DimLoopReturn(coords[i] * param.coords[i]) }

	// Dot product (although please use dot(a,b) for clarity)
	T dot(const vec<T,dim>& param) 
	{
		T dp = 0;
		DimLoop( dp += coords[i] * param.coords[i] )
		return dp;
	}

	// Magnitude
	T mag() const
	{
		T accum = 0;
		DimLoop( accum += coords[i] * coords[i] )
		return sqrt(accum);
	}

	T magsq() const
	{
		T accum = 0;
		DimLoop( accum += coords[i] * coords[i] )
		return accum;
	}

	// Projection of this onto b, in the direction of b
	vec<T,dim> proj (const vec<T,dim>& b) 
		{ return (b/b.mag()) * (((*this)*b)/b.mag()); }
};

// Scalar operations: One glob so we can make this short.
#define VEC_SCALAR_OP(op) \
template <class T, int dim> \
static inline vec<T,dim> operator op (const vec<T,dim>& a, const T b) \
	{ DimLoopReturn(a.coords[i] op b) } \
template <class T, int dim> \
static inline vec<T,dim> operator op (const T b, const vec<T,dim>& a)  \
	{ DimLoopReturn(b op a.coords[i]) } \
template <class T, int dim> \
static inline void operator op##= (vec<T,dim>& a, const T b)  \
	{ DimLoop(a.coords[i] op##= b) } 

// Now implement those scalar <-> vector operations using the macro above.
VEC_SCALAR_OP(+)
VEC_SCALAR_OP(-)
VEC_SCALAR_OP(*)
VEC_SCALAR_OP(/)

// Now specialize down to the 3D case
template <class T>
class vec3 : public vec<T, 3>
{
	public:
	vec3<T>() { }
	vec3<T>(const vec<T,3>& b)
	{ vec<T,3>::operator=(b); }
	vec3<T>(T x, T y, T z)
	{
		this->coords[0] = x;
		this->coords[1] = y;
		this->coords[2] = z;
	}

	const vec3<T>& operator= (const vec3<T>& b)
		{ vec<T,3>::operator=(b); return *this; }
	const vec3<T>& operator= (const vec<T,3>& b)
		{ vec<T,3>::operator=(b); return *this; }
	const vec3<T>& operator= (const T b)
		{ vec<T,3>::operator=(b); return *this; }

	// Cross product
	vec3<T> cross (const vec3<T>& b) const
	{
		vec3<T> ret; 
		ret[0] = this->coords[1]*b.coords[2] - this->coords[2]*b.coords[1]; 
		ret[1] = this->coords[2]*b.coords[0] - this->coords[0]*b.coords[2];
		ret[2] = this->coords[0]*b.coords[1] - this->coords[1]*b.coords[0];
		return ret;
	}
	T x() const { return this->coords[0]; }
	T y() const { return this->coords[1]; }
	T z() const { return this->coords[2]; }
};

// Output routines
template <class T, int dim>
std::ostream& operator<< (std::ostream& os, vec<T,dim>& v)
{
	os << '(' ;
	for (int i=0;i<dim-1;i++)
	    os << v[i] << ", " ;
	return os << v[dim-1] << ')';
}

typedef vec3<float> vec3f;
typedef vec3<double> vec3d;

// It's nicer to write dot(a,b)
template <class T, int dim>
inline T dot(vec<T,dim> a,vec<T,dim> b)
{
	return a.dot(b);
}

// It's nicer to write cross(a,b)
template <class T>
inline vec3<T> cross(vec3<T> a,vec3<T> b)
{
	return a.cross(b);
}

#endif // _VEC_H