#ifndef __SL_VECTOR_3D_H
#define __SL_VECTOR_3D_H

#include "slMathTraits.h"

template<class T>
class vector3D {
 private :
  T data[3];
 public :
  static vector3D<T> ZERO;
  static vector3D<T> X, NEG_X;
  static vector3D<T> Y, NEG_Y;
  static vector3D<T> Z, NEG_Z;
 public :
  T& operator [](int index);
  const T& operator [](int index)const;
 public :
  vector3D<T>& operator  =(const vector3D<T>& v);
  vector3D<T>& operator +=(const vector3D<T>& v);
  vector3D<T>& operator -=(const vector3D<T>& v);
  vector3D<T>& operator *=(T scalar);
  vector3D<T>& operator /=(T scalar);
 public :
  vector3D<T> operator +(const vector3D<T>& v)const;
  vector3D<T> operator -(const vector3D<T>& v)const;
  vector3D<T> operator *(const vector3D<T>& v)const;
  vector3D<T> operator /(const vector3D<T>& v)const;
  vector3D<T> operator *(T scalar)const;
  vector3D<T> operator /(T scalar)const;
 public :
  vector3D<T> operator +(void)const;
  vector3D<T> operator -(void)const;
 public :
  bool operator ==(const vector3D<T>& v)const;
  bool operator !=(const vector3D<T>& v)const;
  bool operator !(void)const;
 public :
  const T* c_ptr(void)const;
  T length(void)const;
  T length_squared(void)const;
  void normalize(void);
 public :
  T dot(const vector3D<T>& v);
  vector3D<T> cross(const vector3D<T>& v);
  T scalar_triple(const vector3D<T>& v1, const vector3D<T>& v2);
  vector3D<T> vector_triple(const vector3D<T>& v1, const vector3D<T>& v2);
 public :
  friend T abs(const vector3D<T>& v);
  friend T dot_product(const vector3D<T>& v1, const vector3D<T>& v2);
  friend vector3D<T> cross_product(const vector3D<T>& v1, const vector3D<T>& v2);
  friend T scalar_triple_product(const vector3D<T>& v1, const vector3D<T>& v2, const vector3D<T>& v3);
  friend vector3D<T> vector_triple_product(const vector3D<T>& v1, const vector3D<T>& v2, const vector3D<T>& v3);
  friend vector3D<T> min(const vector3D<T>& v1, const vector3D<T>& v2);
  friend vector3D<T> max(const vector3D<T>& v1, const vector3D<T>& v2);
 public :
  vector3D();
  vector3D(T x, T y, T z);
  explicit vector3D(T fill);
  explicit vector3D(const T* v);
  vector3D(const vector3D& v);
 ~vector3D();
};

template<class T> vector3D<T> vector3D<T>::ZERO = vector3D<T>(math_traits<T>::zero(), math_traits<T>::zero(), math_traits<T>::zero());
template<class T> vector3D<T> vector3D<T>::X = vector3D<T>(math_traits<T>::one(), math_traits<T>::zero(), math_traits<T>::zero());
template<class T> vector3D<T> vector3D<T>::Y = vector3D<T>(math_traits<T>::zero(), math_traits<T>::one(), math_traits<T>::zero());
template<class T> vector3D<T> vector3D<T>::Z = vector3D<T>(math_traits<T>::zero(), math_traits<T>::zero(), math_traits<T>::one());
template<class T> vector3D<T> vector3D<T>::NEG_X = vector3D<T>(math_traits<T>::neg_one(), math_traits<T>::zero(), math_traits<T>::zero());
template<class T> vector3D<T> vector3D<T>::NEG_Y = vector3D<T>(math_traits<T>::zero(), math_traits<T>::neg_one(), math_traits<T>::zero());
template<class T> vector3D<T> vector3D<T>::NEG_Z = vector3D<T>(math_traits<T>::zero(), math_traits<T>::zero(), math_traits<T>::neg_one());

template<class T>
inline vector3D<T>::vector3D()
{
}

template<class T>
inline vector3D<T>::vector3D(T x, T y, T z)
{
 data[0] = x;
 data[1] = y;
 data[2] = z;
}

template<class T>
inline vector3D<T>::vector3D(T fill)
{
 data[0] = fill;
 data[1] = fill;
 data[2] = fill;
}

template<class T>
inline vector3D<T>::vector3D(const T* v)
{
 data[0] = v[0];
 data[1] = v[1];
 data[2] = v[2];
}

template<class T>
inline vector3D<T>::vector3D(const vector3D<T>& v)
{
 data[0] = v.data[0];
 data[1] = v.data[1];
 data[2] = v.data[2];
}

template<class T>
inline vector3D<T>::~vector3D()
{
}

template<class T>
inline vector3D<T>& vector3D<T>::operator =(const vector3D<T>& v)
{
 if(this == &v) return *this;
 data[0] = v.data[0];
 data[1] = v.data[1];
 data[2] = v.data[2];
 return *this;
}

template<class T>
inline vector3D<T>& vector3D<T>::operator +=(const vector3D<T>& v)
{
 data[0] += v.data[0];
 data[1] += v.data[1];
 data[2] += v.data[2];
 return *this;
}

template<class T>
inline vector3D<T>& vector3D<T>::operator -=(const vector3D<T>& v)
{
 data[0] -= v.data[0];
 data[1] -= v.data[1];
 data[2] -= v.data[2];
 return *this;
}

template<class T>
inline vector3D<T>& vector3D<T>::operator *=(T scalar)
{
 data[0] *= scalar;
 data[1] *= scalar;
 data[2] *= scalar;
 return *this;
}

template<class T>
inline vector3D<T>& vector3D<T>::operator /=(T scalar)
{
 scalar = math_traits<T>::one()/scalar;
 data[0] *= scalar;
 data[1] *= scalar;
 data[2] *= scalar;
 return *this;
}

template<class T>
inline vector3D<T> vector3D<T>::operator +(const vector3D<T>& v)const
{
 return vector3D<T>(data[0] + v.data[0], data[1] + v.data[1], data[2] + v.data[2]);
}

template<class T>
inline vector3D<T> vector3D<T>::operator -(const vector3D<T>& v)const
{
 return vector3D<T>(data[0] - v.data[0], data[1] - v.data[1], data[2] - v.data[2]);
}

template<class T>
inline vector3D<T> vector3D<T>::operator *(const vector3D<T>& v)const
{
 return vector3D<T>(data[0] * v.data[0], data[1] * v.data[1], data[2] * v.data[2]);
}

template<class T>
inline vector3D<T> vector3D<T>::operator /(const vector3D<T>& v)const
{
 T scalar[3] = {
   math_traits<T>::one()/v.data[0],
   math_traits<T>::one()/v.data[1],
   math_traits<T>::one()/v.data[2]
 }; 
 return vector3D<T>(data[0] * scalar[0], data[1] * scalar[1], data[2] * scalar[2]);
}

template<class T>
inline vector3D<T> vector3D<T>::operator *(T scalar)const
{
 return vector3D<T>(data[0]*scalar, data[1]*scalar, data[2]*scalar);
}

template<class T>
inline vector3D<T> vector3D<T>::operator /(T scalar)const
{
 scalar = math_traits<T>::one()/scalar;
 return vector3D<T>(data[0]*scalar, data[1]*scalar, data[2]*scalar);
}

template<class T>
inline vector3D<T> vector3D<T>::operator +(void)const
{
 return *this;
}

template<class T>
inline vector3D<T> vector3D<T>::operator -(void)const
{
 return vector3D<T>(-data[0], -data[1], -data[2]);
}

template<class T>
inline bool vector3D<T>::operator ==(const vector3D<T>& v)const
{
 return ((data[0] == v.data[0]) &&
         (data[1] == v.data[1]) &&
         (data[2] == v.data[2]));
}

template<class T>
inline bool vector3D<T>::operator !=(const vector3D<T>& v)const
{
 return !(*this == v);
}

template<class T>
inline bool vector3D<T>::operator !(void)const
{
 return ((data[0] == math_traits<T>::zero()) &&
         (data[1] == math_traits<T>::zero()) &&
         (data[2] == math_traits<T>::zero()));
}

template<class T>
inline T& vector3D<T>::operator [](int index)
{
 return data[index];
}

template<class T>
inline const T& vector3D<T>::operator [](int index)const
{
 return data[index];
}

template<class T>
inline const T* vector3D<T>::c_ptr(void)const
{
 return &data[0];
}

template<class T>
inline T vector3D<T>::length(void)const
{
 return std::sqrt(data[0]*data[0] + data[1]*data[1] + data[2]*data[2]);
}

template<class T>
inline T vector3D<T>::length_squared(void)const
{
 return data[0]*data[0] + data[1]*data[1] + data[2]*data[2];
}

template<class T>
inline void vector3D<T>::normalize(void)
{
 T denom = math_traits<T>::one()/length();
 data[0] *= denom;
 data[1] *= denom;
 data[2] *= denom;
}

template<class T>
inline T vector3D<T>::dot(const vector3D<T>& v)
{
 return data[0]*v.data[0] + data[1]*v.data[1] + data[2]*v.data[2];
}

template<class T>
inline vector3D<T> vector3D<T>::cross(const vector3D<T>& v)
{
 vector3D<T> r;
 r.data[0] = data[1]*v.data[2] - data[2]*v.data[1];
 r.data[1] = data[2]*v.data[0] - data[0]*v.data[2];
 r.data[2] = data[0]*v.data[1] - data[1]*v.data[0];
 return r;
}

template<class T>
inline T vector3D<T>::scalar_triple(const vector3D<T>& v1, const vector3D<T>& v2)
{
 T temp[3] = {
  v1.data[1]*v2.data[2] - v1.data[2]*v2.data[1],
  v1.data[2]*v2.data[0] - v1.data[0]*v2.data[2],
  v1.data[0]*v2.data[1] - v1.data[1]*v2.data[0]
 };
 return (data[0]*temp[0] + data[1]*temp[1] + data[2]*temp[2]);
}

template<class T>
inline vector3D<T> vector3D<T>::vector_triple(const vector3D<T>& v1, const vector3D<T>& v2)
{
 T temp[3] = {
  v1.data[1]*v2.data[2] - v1.data[2]*v2.data[1],
  v1.data[2]*v2.data[0] - v1.data[0]*v2.data[2],
  v1.data[0]*v2.data[1] - v1.data[1]*v2.data[0]
 };
 vector3D<T> r;
 r.data[0] = data[1]*temp[2] - data[2]*temp[1];
 r.data[1] = data[2]*temp[0] - data[0]*temp[2];
 r.data[2] = data[0]*temp[1] - data[1]*temp[0];
 return r;
}

template<class T>
inline T abs(const vector3D<T>& v)
{
 return v.length();
}

template<class T>
inline T dot_product(const vector3D<T>& v1, const vector3D<T>& v2)
{
 return v1.dot(v2);
}

template<class T>
inline vector3D<T> cross_product(const vector3D<T>& v1, const vector3D<T>& v2)
{
 return v1.cross(v2);
}

template<class T>
inline T scalar_triple_product(const vector3D<T>& v1, const vector3D<T>& v2, const vector3D<T>& v3)
{
 T temp[3] = {
  v2.data[1]*v3.data[2] - v2.data[2]*v3.data[1],
  v2.data[2]*v3.data[0] - v2.data[0]*v3.data[2],
  v2.data[0]*v3.data[1] - v2.data[1]*v3.data[0]
 };
 return (v1.data[0]*temp[0] + v1.data[1]*temp[1] + v1.data[2]*temp[2]);
}

template<class T>
inline vector3D<T> vector_triple_product(const vector3D<T>& v1, const vector3D<T>& v2, const vector3D<T>& v3)
{
 T temp[3] = {
  v2.data[1]*v3.data[2] - v2.data[2]*v3.data[1],
  v2.data[2]*v3.data[0] - v2.data[0]*v3.data[2],
  v2.data[0]*v3.data[1] - v2.data[1]*v3.data[0]
 };
 vector3D<T> r;
 r.data[0] = v1.data[1]*temp[2] - v1.data[2]*temp[1];
 r.data[1] = v1.data[2]*temp[0] - v1.data[0]*temp[2];
 r.data[2] = v1.data[0]*temp[1] - v1.data[1]*temp[0];
 return r;
}

template<class T>
inline vector3D<T> min(const vector3D<T>& v1, const vector3D<T>& v2)
{
 return vector3D<T>(
  (v1.data[0] < v2.data[0]) ? v1.data[0] : v2.data[0], 
  (v1.data[1] < v2.data[1]) ? v1.data[1] : v2.data[1],
  (v1.data[2] < v2.data[2]) ? v1.data[2] : v2.data[2]);
}

template<class T>
inline vector3D<T> max(const vector3D<T>& v1, const vector3D<T>& v2)
{
 return vector3D<T>(
  (v2.data[0] < v1.data[0]) ? v1.data[0] : v2.data[0], 
  (v2.data[1] < v1.data[1]) ? v1.data[1] : v2.data[1],
  (v2.data[2] < v1.data[2]) ? v1.data[2] : v2.data[2]);
}

template<class T>
inline vector3D<T> operator *(T scalar, const vector3D<T>& v)
{
 return v*scalar;
}

template<class T>
inline std::ostream& operator <<(std::ostream& os, const vector3D<T>& v)
{
 os << "<" << v.data[0] << "," << v.data[1] << "," << v.data[2] << ">";
 return os;
}

#endif
