////////////////////////////////////////////////////////////////////////////////
//
//  File:        vect.h
//  Created by:  Jerzy Majewski - jmajewsk@meil.pw.edu.pl
//  Modified by: 
//    2005-07-28 JM: modification of definitions of friend template operators
//					 in order to get rid of warning generated by gcc 
//					 (tested with ver. 3.3.3 and MSVC++ .NET)
//
// Copyright notice:
//    Copyright (C) 2000-2003 Jerzy Majewski 
// 
//    This is free software. You can redistribute it and/or
//    modify it under the terms of the GNU Lesser General Public License
//    as published by the Free Software Foundation.
//
//    This software is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY. See the GNU Lesser General Public License 
//    for more details.
//
////////////////////////////////////////////////////////////////////////////////

#ifndef __VECT_H__
#define __VECT_H__


#include "mgdecl.h"


template <class ELEM_TYPE, Dimension DIM> 
class Vect;

template <class ELEM_TYPE, Dimension DIM> 
Vect<ELEM_TYPE,DIM> operator *( const ELEM_TYPE&,  const Vect<ELEM_TYPE,DIM>&);

template <class ELEM_TYPE, Dimension DIM> 
Vect<ELEM_TYPE,DIM> operator *( const Vect<ELEM_TYPE,DIM>&, const ELEM_TYPE&);

template <class ELEM_TYPE, Dimension DIM> 
Vect<ELEM_TYPE,DIM> operator /( const Vect<ELEM_TYPE,DIM>&, const ELEM_TYPE&);

template <class ELEM_TYPE, Dimension DIM> 
ELEM_TYPE	operator *( const Vect<ELEM_TYPE,DIM>&, const Vect<ELEM_TYPE,DIM>&);   // dot product

template <class ELEM_TYPE, Dimension DIM> 
Vect<ELEM_TYPE,DIM> operator +( const Vect<ELEM_TYPE,DIM>&, const Vect<ELEM_TYPE,DIM>&);

template <class ELEM_TYPE, Dimension DIM> 
Vect<ELEM_TYPE,DIM> operator -( const Vect<ELEM_TYPE,DIM>&, const Vect<ELEM_TYPE,DIM>&);

template <class ELEM_TYPE, Dimension DIM> 
Vect<ELEM_TYPE,DIM> operator %( const Vect<ELEM_TYPE,DIM>&, const Vect<ELEM_TYPE,DIM>&);   // vector product




//////////////////////////////////////////////////////////////////////
//	class Vect
//////////////////////////////////////////////////////////////////////

template <class ELEM_TYPE, Dimension DIM> 
class Vect
{
public:
    //////////////////////////////////////////////////////////////////
	// constructors
	Vect( const ELEM_TYPE& x);
	Vect( const ELEM_TYPE& x, const ELEM_TYPE& y);
	Vect( const ELEM_TYPE& x, const ELEM_TYPE& y, const ELEM_TYPE& z);
	Vect( const ELEM_TYPE& x, const ELEM_TYPE& y, const ELEM_TYPE& z, const ELEM_TYPE& w);
	Vect( const Vect<ELEM_TYPE, DIM>& vec);
	Vect();

    //////////////////////////////////////////////////////////////////
	// declaration of friend two argument operators
	friend Vect<ELEM_TYPE, DIM>	operator  *<>( const ELEM_TYPE&,  const Vect<ELEM_TYPE, DIM>&);
	friend Vect<ELEM_TYPE, DIM> operator  *<>( const Vect<ELEM_TYPE, DIM>&, const ELEM_TYPE&);
	friend Vect<ELEM_TYPE, DIM> operator  /<>( const Vect<ELEM_TYPE, DIM>&, const ELEM_TYPE&);
	friend ELEM_TYPE			operator  *<>( const Vect<ELEM_TYPE, DIM>&, const Vect<ELEM_TYPE, DIM>&);   //scalar mult.
	friend Vect<ELEM_TYPE, DIM>	operator  +<>( const Vect<ELEM_TYPE, DIM>&, const Vect<ELEM_TYPE, DIM>&);
	friend Vect<ELEM_TYPE, DIM>	operator  -<>( const Vect<ELEM_TYPE, DIM>&, const Vect<ELEM_TYPE, DIM>&);
	friend Vect<ELEM_TYPE, DIM>	operator  %<>( const Vect<ELEM_TYPE, DIM>&, const Vect<ELEM_TYPE, DIM>&);   //Vect mult.

    //////////////////////////////////////////////////////////////////
	// one argument operators
	Vect<ELEM_TYPE, DIM>&	operator  =( const Vect<ELEM_TYPE, DIM> &vec);
	Vect<ELEM_TYPE, DIM>&	operator +=( const Vect<ELEM_TYPE, DIM>&);
	Vect<ELEM_TYPE, DIM>&	operator -=( const Vect<ELEM_TYPE, DIM>&);
	Vect<ELEM_TYPE, DIM>&	operator *=( const ELEM_TYPE&);
	Vect<ELEM_TYPE, DIM>&	operator /=( const ELEM_TYPE&);

    //////////////////////////////////////////////////////////////////
	// misc functions
	ELEM_TYPE				module() const;
	Vect<ELEM_TYPE, DIM>	versor() const;
	
	
	const ELEM_TYPE&	X( const MGInt& i) const	{ return mtab[i];}
	ELEM_TYPE&			rX( const MGInt& i)			{ return mtab[i];}

	const ELEM_TYPE&	X() const	{ return mtab[0];}
	const ELEM_TYPE&	Y() const	{ ASSERT(DIM>1); return mtab[1];}
	const ELEM_TYPE&	Z() const	{ ASSERT(DIM>2); return mtab[2];}
	const ELEM_TYPE&	W() const	{ ASSERT(DIM>3); return mtab[3];}
	
	ELEM_TYPE&			rX()		{ return mtab[0];}
	ELEM_TYPE&			rY()		{ ASSERT(DIM>1); return mtab[1];}
	ELEM_TYPE&			rZ()		{ ASSERT(DIM>2); return mtab[2];}
	ELEM_TYPE&			rW()		{ ASSERT(DIM>3); return mtab[3];}

protected:
	ELEM_TYPE mtab[DIM];
};


//////////////////////////////////////////////////////////////////////
typedef Vect<MGFloat,DIM_1D> Vect1D;
typedef Vect<MGFloat,DIM_2D> Vect2D;
typedef Vect<MGFloat,DIM_3D> Vect3D;
typedef Vect<MGFloat,DIM_4D> Vect4D;




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

template <class ELEM_TYPE, Dimension DIM> 
inline Vect<ELEM_TYPE, DIM>::Vect( const Vect<ELEM_TYPE, DIM> &vec)
{
	ASSERT( DIM>0 && DIM<5);
	for ( MGInt i=0; i<DIM; ++i)
		mtab[i] = vec.mtab[i];
}

template <class ELEM_TYPE, Dimension DIM> 
inline Vect<ELEM_TYPE, DIM>::Vect( const ELEM_TYPE& x)
{
		ASSERT( DIM == DIM_1D);
		mtab[0] = x;
}

template <class ELEM_TYPE, Dimension DIM> 
inline Vect<ELEM_TYPE, DIM>::Vect( const ELEM_TYPE& x, const ELEM_TYPE& y)
{
		ASSERT( DIM == DIM_2D);
		mtab[0] = x;
		mtab[1] = y;
}

template <class ELEM_TYPE, Dimension DIM> 
inline Vect<ELEM_TYPE, DIM>::Vect( const ELEM_TYPE& x, const ELEM_TYPE& y, const ELEM_TYPE& z)
{
		ASSERT( DIM == DIM_3D);
		mtab[0] = x;
		mtab[1] = y;
		mtab[2] = z;
}

template <class ELEM_TYPE, Dimension DIM> 
inline Vect<ELEM_TYPE, DIM>::Vect( const ELEM_TYPE& x, const ELEM_TYPE& y, const ELEM_TYPE& z, const ELEM_TYPE& w)
{
		ASSERT( DIM == DIM_4D);
		mtab[0] = x;
		mtab[1] = y;
		mtab[2] = z;
		mtab[3] = w;
}


template <class ELEM_TYPE, Dimension DIM> 
inline Vect<ELEM_TYPE, DIM>::Vect()
{
	ASSERT( DIM > DIM_NONE && DIM <= DIM_4D);
	for ( MGInt i=0; i<DIM; mtab[i++]=(ELEM_TYPE)0.0);
}


template <class ELEM_TYPE, Dimension DIM> 
inline Vect<ELEM_TYPE, DIM>& Vect<ELEM_TYPE, DIM>::operator =( const Vect<ELEM_TYPE, DIM> &vec)
{
	for ( MGInt i=0; i<DIM; ++i)
		mtab[i] = vec.mtab[i];
	return *this;
}

template <class ELEM_TYPE, Dimension DIM> 
inline Vect<ELEM_TYPE, DIM>& Vect<ELEM_TYPE, DIM>::operator+=( const Vect<ELEM_TYPE, DIM>& vec)
{
	for ( MGInt i=0; i<DIM; ++i)
		mtab[i] += vec.mtab[i];
	return *this;
}


template <class ELEM_TYPE, Dimension DIM> 
inline Vect<ELEM_TYPE, DIM>& Vect<ELEM_TYPE, DIM>::operator-=( const Vect<ELEM_TYPE, DIM>& vec)
{
	for ( MGInt i=0; i<DIM; ++i)
		mtab[i] -= vec.mtab[i];
	return *this;
}

template <class ELEM_TYPE, Dimension DIM> 
inline Vect<ELEM_TYPE, DIM>& Vect<ELEM_TYPE, DIM>::operator*=( const ELEM_TYPE& doub)
{
	for ( MGInt i=0; i<DIM; ++i)
		mtab[i] *= doub;
	return *this;
}

template <class ELEM_TYPE, Dimension DIM> 
inline Vect<ELEM_TYPE, DIM>& Vect<ELEM_TYPE, DIM>::operator/=( const ELEM_TYPE& doub)
{
	for ( MGInt i=0; i<DIM; ++i)
		mtab[i] /= doub;
	return *this;
}

template <class ELEM_TYPE, Dimension DIM> 
inline ELEM_TYPE Vect<ELEM_TYPE, DIM>::module() const
{
	return sqrt( (*this)*(*this));
}

template <class ELEM_TYPE, Dimension DIM> 
inline Vect<ELEM_TYPE, DIM> Vect<ELEM_TYPE, DIM>::versor() const
{
	return (*this / module() );
}



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


inline bool operator == ( const Vect2D& c1, const Vect2D& c2 )
{
	if ( fabs( c1.X() - c2.X() ) < ZERO &&
		 fabs( c1.Y() - c2.Y() ) < ZERO )

		return true;
	else
		return false;
}


inline bool operator < ( const Vect2D& c1, const Vect2D& c2 )
{
	if ( fabs( c1.X() - c2.X() ) > ZERO )
	{
		if ( c1.X() + ZERO < c2.X() )
			return true;
		else 
			return false;
	}
	else
	{
		if ( fabs( c1.Y() - c2.Y() ) > ZERO )
		{
			if ( c1.Y() + ZERO < c2.Y() )
				return true;
			else 
				return false;
		}
		else
		{
			return false;
		}
	}
}




inline bool operator == ( const Vect3D& c1, const Vect3D& c2 )
{
	if ( fabs( c1.X() - c2.X() ) < ZERO &&
		 fabs( c1.Y() - c2.Y() ) < ZERO &&
		 fabs( c1.Z() - c2.Z() ) < ZERO )

		return true;
	else
		return false;
}


inline bool operator < ( const Vect3D& c1, const Vect3D& c2 )
{
	if ( fabs( c1.X() - c2.X() ) > ZERO )
	{
		if ( c1.X() + ZERO < c2.X() )
			return true;
		else 
			return false;
	}
	else
	{
		if ( fabs( c1.Y() - c2.Y() ) > ZERO )
		{
			if ( c1.Y() + ZERO < c2.Y() )
				return true;
			else 
				return false;
		}
		else
		{
			if ( fabs( c1.Z() - c2.Z() ) > ZERO )
			{
				if ( c1.Z() + ZERO < c2.Z() )
					return true;
				else 
					return false;
			}
			else
			{
				return false;
			}
		}
	}
}



//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////
// specialization of the Vect operators for DIM=1,2,3,4


//////////////////////////////////////////////////////////////////////
// Vect1D
//////////////////////////////////////////////////////////////////////
template <>
inline Vect1D operator*( const MGFloat& comp, const Vect1D& vec)
{
	return Vect1D( vec.mtab[0] * comp );
}


template <>
inline Vect1D operator*( const Vect1D& vec, const MGFloat& comp)
{
	return Vect1D( vec.mtab[0] * comp );
}


template <>
inline Vect1D operator/( const Vect1D& vec, const MGFloat& comp)
{
	return Vect1D( vec.mtab[0] / comp );
}


template <>
inline MGFloat operator*( const Vect1D& vec1, const Vect1D& vec2)
{
	return vec1.mtab[0] * vec2.mtab[0];
}


template <>
inline Vect1D operator+( const Vect1D& vec1, const Vect1D& vec2)
{
	return Vect1D( vec1.mtab[0] + vec2.mtab[0] );
}


template <>
inline Vect1D operator-( const Vect1D& vec1, const Vect1D& vec2)
{
	return Vect1D( vec1.mtab[0] - vec2.mtab[0] );
}




//////////////////////////////////////////////////////////////////////
// Vect2D
//////////////////////////////////////////////////////////////////////
template <>
inline Vect2D operator*( const MGFloat& comp, const Vect2D& vec)
{
	return Vect2D( vec.mtab[0] * comp, vec.mtab[1] * comp );
}


template <>
inline Vect2D operator*( const Vect2D& vec, const MGFloat& comp)
{
	return Vect2D( vec.mtab[0] * comp, vec.mtab[1] * comp );
}


template <>
inline Vect2D operator/( const Vect2D& vec, const MGFloat& comp)
{
	return Vect2D( vec.mtab[0] / comp, vec.mtab[1] / comp );
}


template <>
inline MGFloat operator*( const Vect2D& vec1, const Vect2D& vec2)
{
	return vec1.mtab[0] * vec2.mtab[0] + vec1.mtab[1] * vec2.mtab[1];
}


template <>
inline Vect2D operator+( const Vect2D& vec1, const Vect2D& vec2)
{
	return Vect2D( vec1.mtab[0] + vec2.mtab[0], vec1.mtab[1] + vec2.mtab[1] );
}


template <>
inline Vect2D operator-( const Vect2D& vec1, const Vect2D& vec2)
{
	return Vect2D( vec1.mtab[0] - vec2.mtab[0], vec1.mtab[1] - vec2.mtab[1] );
}


//////////////////////////////////////////////////////////////////////
// Vect3D
//////////////////////////////////////////////////////////////////////
template <>
inline Vect3D operator*( const MGFloat& comp, const Vect3D& vec)
{
	return Vect3D( vec.mtab[0] * comp, vec.mtab[1] * comp, vec.mtab[2] * comp );
}


template <>
inline Vect3D operator*( const Vect3D& vec, const MGFloat& comp)
{
	return Vect3D( vec.mtab[0] * comp, vec.mtab[1] * comp, vec.mtab[2] * comp );
}


template <>
inline Vect3D operator/( const Vect3D& vec, const MGFloat& comp)
{
	return Vect3D( vec.mtab[0] / comp, vec.mtab[1] / comp, vec.mtab[2] / comp );
}


template <>
inline MGFloat operator*( const Vect3D& vec1, const Vect3D& vec2)
{
	return vec1.mtab[0]*vec2.mtab[0] + vec1.mtab[1]*vec2.mtab[1] + vec1.mtab[2]*vec2.mtab[2];
}


template <>
inline Vect3D operator+( const Vect3D& vec1, const Vect3D& vec2)
{
	return Vect3D( vec1.mtab[0] + vec2.mtab[0], vec1.mtab[1] + vec2.mtab[1], vec1.mtab[2] + vec2.mtab[2] );
}


template <>
inline Vect3D operator-( const Vect3D& vec1, const Vect3D& vec2)
{
	return Vect3D( vec1.mtab[0] - vec2.mtab[0], vec1.mtab[1] - vec2.mtab[1], vec1.mtab[2] - vec2.mtab[2] );
}


template <>
inline Vect3D operator%( const Vect3D& vec1, const Vect3D& vec2)
{
	MGFloat x,y,z;
	x = vec1.mtab[1]*vec2.mtab[2] - vec1.mtab[2]*vec2.mtab[1];
	y = vec1.mtab[2]*vec2.mtab[0] - vec1.mtab[0]*vec2.mtab[2];
	z = vec1.mtab[0]*vec2.mtab[1] - vec1.mtab[1]*vec2.mtab[0];
	return Vect3D(x,y,z);
}



//////////////////////////////////////////////////////////////////////
// Vect4D
//////////////////////////////////////////////////////////////////////
template <>
inline Vect4D operator*( const MGFloat& comp, const Vect4D& vec)
{
	return Vect4D( vec.mtab[0] * comp, vec.mtab[1] * comp, vec.mtab[2] * comp, vec.mtab[3] * comp );
}


template <>
inline Vect4D operator*( const Vect4D& vec, const MGFloat& comp)
{
	return Vect4D( vec.mtab[0] * comp, vec.mtab[1] * comp, vec.mtab[2] * comp, vec.mtab[3] * comp );
}


template <>
inline Vect4D operator/( const Vect4D& vec, const MGFloat& comp)
{
	return Vect4D( vec.mtab[0] / comp, vec.mtab[1] / comp, vec.mtab[2] / comp, vec.mtab[3] / comp );
}


template <>
inline MGFloat operator*( const Vect4D& vec1, const Vect4D& vec2)
{
	return vec1.mtab[0]*vec2.mtab[0] + vec1.mtab[1]*vec2.mtab[1] + vec1.mtab[2]*vec2.mtab[2] + vec1.mtab[3]*vec2.mtab[3];
}


template <>
inline Vect4D operator+( const Vect4D& vec1, const Vect4D& vec2)
{
	return Vect4D( vec1.mtab[0] + vec2.mtab[0], vec1.mtab[1] + vec2.mtab[1], vec1.mtab[2] + vec2.mtab[2], vec1.mtab[3] + vec2.mtab[3] );
}


template <>
inline Vect4D operator-( const Vect4D& vec1, const Vect4D& vec2)
{
	return Vect4D( vec1.mtab[0] - vec2.mtab[0], vec1.mtab[1] - vec2.mtab[1], vec1.mtab[2] - vec2.mtab[2], vec1.mtab[3] - vec2.mtab[3] );
}



//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////
template <class ELEM_TYPE, Dimension DIM> 
inline Vect<ELEM_TYPE,DIM> operator-( const Vect<ELEM_TYPE,DIM>& vec1, const Vect<ELEM_TYPE,DIM>& vec2)
{
	Vect<ELEM_TYPE,DIM>	v;
	for ( MGInt i=0; i<DIM; ++i)
		v.mtab[i] = vec1.mtab[i] - vec2.mtab[i];
	return v;
}

template <class ELEM_TYPE, Dimension DIM> 
inline Vect<ELEM_TYPE,DIM> operator+( const Vect<ELEM_TYPE,DIM>& vec1, const Vect<ELEM_TYPE,DIM>& vec2)
{
	Vect<ELEM_TYPE,DIM>	v;
	for ( MGInt i=0; i<DIM; ++i)
		v.mtab[i] = vec1.mtab[i] + vec2.mtab[i];
	return v;
}

template <class ELEM_TYPE, Dimension DIM> 
inline ELEM_TYPE operator*( const Vect<ELEM_TYPE,DIM>& vec1, const Vect<ELEM_TYPE,DIM>& vec2)
{
	ELEM_TYPE	e(0);
	for ( MGInt i=0; i<DIM; ++i)
		e += vec1.mtab[i] * vec2.mtab[i];
	return e;
}

template <class ELEM_TYPE, Dimension DIM> 
inline Vect<ELEM_TYPE,DIM> operator*( const ELEM_TYPE& e, const Vect<ELEM_TYPE,DIM>& vec)
{
	Vect<ELEM_TYPE,DIM>	v;
	for ( MGInt i=0; i<DIM; ++i)
		v.mtab[i] = e * vec.mtab[i];
	return v;
}

template <class ELEM_TYPE, Dimension DIM> 
inline Vect<ELEM_TYPE,DIM> operator*( const Vect<ELEM_TYPE,DIM>& vec, const ELEM_TYPE& e)
{
	Vect<ELEM_TYPE,DIM>	v;
	for ( MGInt i=0; i<DIM; ++i)
		v.mtab[i] = e * vec.mtab[i];
	return v;
}

template <class ELEM_TYPE, Dimension DIM> 
inline Vect<ELEM_TYPE,DIM> operator/( const Vect<ELEM_TYPE,DIM>& vec, const ELEM_TYPE& e)
{
	Vect<ELEM_TYPE,DIM>	v;
	for ( MGInt i=0; i<DIM; ++i)
		v.mtab[i] = vec.mtab[i] / e;
	return v;
}



#endif // __VECT_H__