#ifndef _SimpleVec_h__
#define _SimpleVec_h__


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


/** Vector class.
*/
class SimpleVec
{
public:

	///////////
	//-- members;

	float x; 
	float y;
	float z;


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

	void  SetX(float q) { x = q; }
	void  SetY(float q) { y = q; }
	void  SetZ(float q) { z = q; }

	float GetX() const { return x; }
	float GetY() const { return y; }
	float GetZ() const { return z; }


	//////////////
	//-- constructors

	SimpleVec() { }

	SimpleVec(float x, float y, float z): x(x), y(y), z(z)
	{}

	explicit SimpleVec(float v): x(v), y(v), z(v){}
	/** Copy constructor.
	*/
	SimpleVec(const SimpleVec &v): x(v.x), y(v.y), z(v.z)
	{}

	// Functions to get at the std::vector components
	float& operator[] (const int inx) 
	{
		return (&x)[inx];
	}

	operator const float*() const 
	{ 
		return (const float*) this; 
	}

	const float &operator[] (const int inx) const 
	{
		return *(&x + inx);
	}

	void ExtractVerts(float *px, float *py, int which) const;
  
	void SetValue(float a, float b, float c)
	{  
		x = a; y = b; z = c; 
	}

	void SetValue(float a) 
	{ 
		x = y = z = a; 
	}
 
	
	/** Returns the axis, where the std::vector has the largest value.
	*/
	int DrivingAxis(void) const;
	/** returns the axis, where the std::vector has the smallest value
	*/
	int TinyAxis(void) const;
	/** Returns largest component in this vector
	*/
	inline float MaxComponent(void) const 
	{
		return (x > y) ? ( (x > z) ? x : z) : ( (y > z) ? y : z);
	}
	/** Returns copy of this vector where all components are positiv.
	*/
	inline SimpleVec Abs(void) const 
	{
		return SimpleVec(fabs(x), fabs(y), fabs(z));
	}
	/** normalizes the std::vector of unit size corresponding to given std::vector
	*/
	inline void Normalize();
	/** Returns false if this std::vector has a nan component.
	*/
	bool CheckValidity() const;
	/**
	Find a right handed coordinate system with (*this) being
	the z-axis. For a right-handed system, U x V = (*this) holds.
	This implementation is here to avoid inconsistence and confusion
	when construction coordinate systems using ArbitraryNormal():
	In fact:
	V = ArbitraryNormal(N);
	U = CrossProd(V,N);
	constructs a right-handed coordinate system as well, BUT:
	1) bugs can be introduced if one mistakenly constructs a 
	left handed sytems e.g. by doing
	U = ArbitraryNormal(N);
	V = CrossProd(U,N);
	2) this implementation gives non-negative base vectors
	for (*this)==(0,0,1) |  (0,1,0) | (1,0,0), which is
	good for debugging and is not the case with the implementation
	using ArbitraryNormal()

	===> Using ArbitraryNormal() for constructing coord systems 
	is obsoleted by this method (<JK> 12/20/03).   
	*/
	void RightHandedBase(SimpleVec& U, SimpleVec& V) const;
	
	
	// Unary operators

	SimpleVec operator+ () const;
	SimpleVec operator- () const;

	// Assignment operators
	
	SimpleVec& operator+= (const SimpleVec &A);
	SimpleVec& operator-= (const SimpleVec &A);
	SimpleVec& operator*= (const SimpleVec &A);
	SimpleVec& operator*= (float A);
	SimpleVec& operator/= (float A);

	static inline SimpleVec UNIT_X() { return SimpleVec(1, 0, 0); }
	static inline SimpleVec UNIT_Y() { return SimpleVec(0, 1, 0); }
	static inline SimpleVec UNIT_Z() { return SimpleVec(0, 0, 1); }
	static inline SimpleVec ZERO() { return SimpleVec(0, 0, 0); }


	//////////
	//-- friends
	
	/**
	===> Using ArbitraryNormal() for constructing coord systems 
	===> is obsoleted by RightHandedBase() method (<JK> 12/20/03).   

	Return an arbitrary normal to `v'.
	In fact it tries v x (0,0,1) an if the result is too small,
	it definitely does v x (0,1,0). It will always work for
	non-degenareted std::vector and is much faster than to use
	TangentVectors.

	@param v(in) The std::vector we want to find normal for.
	@return The normal std::vector to v.
	*/
	friend inline SimpleVec ArbitraryNormal(const SimpleVec &v); 

	/// Transforms a std::vector to the global coordinate frame.
	/**
	Given a local coordinate frame (U,V,N) (i.e. U,V,N are 
	the x,y,z axes of the local coordinate system) and
	a std::vector 'loc' in the local coordiante system, this
	function returns a the coordinates of the same std::vector
	in global frame (i.e. frame (1,0,0), (0,1,0), (0,0,1).
	*/
	friend inline SimpleVec ToGlobalFrame(const SimpleVec& loc,
		                                const SimpleVec& U,
										const SimpleVec& V,
										const SimpleVec& N);

	/// Transforms a std::vector to a local coordinate frame.
	/**
	Given a local coordinate frame (U,V,N) (i.e. U,V,N are 
	the x,y,z axes of the local coordinate system) and
	a std::vector 'loc' in the global coordiante system, this
	function returns a the coordinates of the same std::vector
	in the local frame.
	*/
	friend inline SimpleVec ToLocalFrame(const SimpleVec& loc,
                                	   const SimpleVec& U,
									   const SimpleVec& V,
									   const SimpleVec& N);

	/// the magnitude=size of the std::vector
	friend inline float Magnitude(const SimpleVec &v);
	/// the squared magnitude of the std::vector .. for efficiency in some cases
	friend inline float SqrMagnitude(const SimpleVec &v);
	/// Magnitude(v1-v2)
	friend inline float Distance(const SimpleVec &v1, const SimpleVec &v2);
	/// SqrMagnitude(v1-v2)
	friend inline float SqrDistance(const SimpleVec &v1, const SimpleVec &v2);
	/// creates the std::vector of unit size corresponding to given std::vector
	friend inline SimpleVec Normalize(const SimpleVec &A);
	

	// construct view vectors .. DirAt is the main viewing direction
	// Viewer is the coordinates of viewer location, UpL is the std::vector.
	friend void ViewVectors(const SimpleVec &DirAt, 
							const SimpleVec &Viewer,
							const SimpleVec &UpL, 
							SimpleVec &ViewV,
							SimpleVec &ViewU, 
							SimpleVec &ViewN);

	// Given the intersection point `P', you have available normal `N'
	// of unit length. Let us suppose the incoming ray has direction `D'.
	// Then we can construct such two vectors `U' and `V' that
	// `U',`N', and `D' are coplanar, and `V' is perpendicular
	// to the vectors `N','D', and `V'. Then 'N', 'U', and 'V' create
	// the orthonormal base in space R3.
	friend void TangentVectors(SimpleVec &U, SimpleVec &V, // output
							   const SimpleVec &normal, // input
							   const SimpleVec &dirIncoming);


	// Binary operators
	
	friend inline SimpleVec operator+ (const SimpleVec &A, const SimpleVec &B);
	friend inline SimpleVec operator- (const SimpleVec &A, const SimpleVec &B);
	friend inline SimpleVec operator* (const SimpleVec &A, const SimpleVec &B);
	friend inline SimpleVec operator* (const SimpleVec &A, float B);
	friend inline SimpleVec operator* (float A, const SimpleVec &B);
	friend inline SimpleVec operator/ (const SimpleVec &A, const SimpleVec &B);

	friend inline int operator< (const SimpleVec &A, const SimpleVec &B);
	friend inline int operator<= (const SimpleVec &A, const SimpleVec &B);

	friend inline SimpleVec operator/ (const SimpleVec &A, float B);
	friend inline int operator== (const SimpleVec &A, const SimpleVec &B);
	friend inline float DotProd(const SimpleVec &A, const SimpleVec &B);
	friend inline SimpleVec CrossProd (const SimpleVec &A, const SimpleVec &B);

	friend std::ostream& operator<< (std::ostream &s, const SimpleVec &A);
	friend std::istream& operator>> (std::istream &s, SimpleVec &A);

	friend void Minimize(SimpleVec &min, const SimpleVec &candidate);
	friend void Maximize(SimpleVec &max, const SimpleVec &candidate);

	friend inline int EpsilonEqualV3(const SimpleVec &v1, const SimpleVec &v2, float thr);
	friend inline int EpsilonEqualV3(const SimpleVec &v1, const SimpleVec &v2);
};


inline SimpleVec ArbitraryNormal(const SimpleVec &N)
{
	float dist2 = N.x * N.x + N.y * N.y;

	if (dist2 > 0.0001f) 
	{
		float inv_size = 1.0f / sqrtf(dist2);
		return SimpleVec(N.y * inv_size, -N.x * inv_size, 0); // N x (0,0,1)
	}

	float inv_size = 1.0f / sqrtf(N.z * N.z + N.x * N.x);
	return SimpleVec(-N.z * inv_size, 0, N.x * inv_size); // N x (0,1,0)
}



inline float Magnitude(const SimpleVec &v) 
{
	return sqrtf(v.x * v.x + v.y * v.y + v.z * v.z);
}


inline float SqrMagnitude(const SimpleVec &v) 
{
	return v.x * v.x + v.y * v.y + v.z * v.z;
}


inline float Distance(const SimpleVec &v1, const SimpleVec &v2)
{
	//return sqrtf(sqrt(v1.x - v2.x) + sqrt(v1.y - v2.y) + sqrt(v1.z - v2.z));
	return Magnitude(v1 - v2);
}


inline float SqrDistance(const SimpleVec &v1, const SimpleVec &v2)
{
	//return sqrt(v1.x - v2.x) + sqrt(v1.y - v2.y) + sqrt(v1.z - v2.z);
	return SqrMagnitude(v1 - v2);
}


inline SimpleVec Normalize(const SimpleVec &A)
{
	return A * (1.0f / Magnitude(A));
}


inline float DotProd(const SimpleVec &A, const SimpleVec &B)
{
	return A.x * B.x + A.y * B.y + A.z * B.z;
}


// angle between two vectors with respect to a surface normal in the
// range [0 .. 2 * pi]
inline float Angle(const SimpleVec &A, const SimpleVec &B, const SimpleVec &norm)
{
	SimpleVec cross = CrossProd(A, B);

	float signedAngle;

	if (DotProd(cross, norm) > 0)
		signedAngle = atan2(-Magnitude(CrossProd(A, B)), DotProd(A, B));
	else
		signedAngle = atan2(Magnitude(CrossProd(A, B)), DotProd(A, B));

	if (signedAngle < 0)
		return 2 * 3.145f + signedAngle;

	return signedAngle;
}


inline SimpleVec SimpleVec::operator+() const
{
	return *this;
}


inline SimpleVec SimpleVec::operator-() const
{
  return SimpleVec(-x, -y, -z);
}


inline SimpleVec &SimpleVec::operator+=(const SimpleVec &A)
{
	x += A.x;  y += A.y;  z += A.z;
	return *this;
}


inline SimpleVec& SimpleVec::operator-=(const SimpleVec &A)
{
	x -= A.x;  y -= A.y;  z -= A.z;
	
	return *this;
}


inline SimpleVec& SimpleVec::operator*= (float A)
{
	x *= A;  y *= A;  z *= A;
	return *this;
}


inline SimpleVec& SimpleVec::operator/=(float A)
{
	float a = 1.0f / A;
	
	x *= a;  y *= a;  z *= a;
	
	return *this;
}


inline SimpleVec& SimpleVec::operator*= (const SimpleVec &A)
{
  x *= A.x;  y *= A.y;  z *= A.z;
  return *this;
}


inline SimpleVec operator+ (const SimpleVec &A, const SimpleVec &B)
{
	return SimpleVec(A.x + B.x, A.y + B.y, A.z + B.z);
}


inline SimpleVec operator- (const SimpleVec &A, const SimpleVec &B)
{
	return SimpleVec(A.x - B.x, A.y - B.y, A.z - B.z);
}


inline SimpleVec operator* (const SimpleVec &A, const SimpleVec &B)
{
	return SimpleVec(A.x * B.x, A.y * B.y, A.z * B.z);
}


inline SimpleVec operator* (const SimpleVec &A, float B)
{
  return SimpleVec(A.x * B, A.y * B, A.z * B);
}


inline SimpleVec operator* (float A, const SimpleVec &B)
{
	return SimpleVec(B.x * A, B.y * A, B.z * A);
}


inline SimpleVec operator/ (const SimpleVec &A, const SimpleVec &B)
{
	return SimpleVec(A.x / B.x, A.y / B.y, A.z / B.z);
}


inline SimpleVec operator/ (const SimpleVec &A, float B)
{
	float b = 1.0f / B;
	return SimpleVec(A.x * b, A.y * b, A.z * b);
}


inline int operator< (const SimpleVec &A, const SimpleVec &B)
{
	return A.x < B.x && A.y < B.y && A.z < B.z;
}


inline int operator<= (const SimpleVec &A, const SimpleVec &B)
{
	return A.x <= B.x && A.y <= B.y && A.z <= B.z;
}


// Might replace floating-point == with comparisons of
// magnitudes, if needed.
inline int operator== (const SimpleVec &A, const SimpleVec &B)
{ 
	return (A.x == B.x) && (A.y == B.y) && (A.z == B.z);
}


inline SimpleVec CrossProd (const SimpleVec &A, const SimpleVec &B)
{
	return SimpleVec(A.y * B.z - A.z * B.y,
		           A.z * B.x - A.x * B.z,
				   A.x * B.y - A.y * B.x);
}


inline void SimpleVec::Normalize()
{
	float sqrmag = x * x + y * y + z * z;
	
	if (sqrmag > 0.0f)
		(*this) *= 1.0f / sqrtf(sqrmag);
}


// Overload << operator for C++-style output
inline std::ostream& operator<< (std::ostream &s, const SimpleVec &A)
{
	return s << "(" << A.x << ", " << A.y << ", " << A.z << ")";
}


// Overload >> operator for C++-style input
inline std::istream& operator>> (std::istream &s, SimpleVec &A)
{
	char a;
	// read "(x, y, z)"
	return s >> a >> A.x >> a >> A.y >> a >> A.z >> a;
}


inline int EpsilonEqualV3(const SimpleVec &v1, const SimpleVec &v2, float thr)
{
	if (fabsf(v1.x-v2.x) > thr)
		return false;
	if (fabsf(v1.y-v2.y) > thr)
		return false;
	if (fabsf(v1.z-v2.z) > thr)
		return false;
	
	return true;
}


inline int EpsilonEqualV3(const SimpleVec &v1, const SimpleVec &v2)
{
	return EpsilonEqualV3(v1, v2, 1e-6f);
}


#endif