/***************************************************************************
 *   Copyright (C) 2007 by A.J. Tavakoli                                   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#ifndef VECTOR3_H
#define VECTOR3_H

#include <iostream>
#include <cmath>

#include "numeric.h"


namespace Math3D {


template<class Real>
class Vector3 {
public:
  // components of vector
  Real m[3];

  inline Vector3() { }
  inline Vector3(const Vector3 &v) { m[0] = v.m[0]; m[1] = v.m[1]; m[2] = v.m[2]; }
  inline Vector3(Real x, Real y, Real z) { m[0] = x; m[1] = y; m[2] = z; }
  inline Vector3(const Real *v) { m[0] = v[0]; m[1] = v[1]; m[2] = v[2]; }

   // overloaded operators  
  inline Real & operator [] (int i);
  inline const Real & operator [] (int i) const;
  inline Vector3 operator + (const Vector3 &v) const;
  inline Vector3 operator - (const Vector3 &v) const;
  inline const Vector3 & operator += (const Vector3 &v);
  inline const Vector3 & operator -= (const Vector3 &v);
  inline Real operator * (const Vector3 &v) const;
  inline Vector3 operator * (Real s) const;
  inline Vector3 operator / (Real s) const;
  inline const Vector3& operator *= (Real s);
  inline const Vector3& operator /= (Real s);

  // methods  
  inline Real mag() const;
  inline void normalize();
  inline Vector3 reflect(const Vector3 &n) const;
  inline void reflect(const Vector3 &n, Vector3 &result);
  inline void set(Real _x, Real _y, Real _z);
  inline Real sqrMag() const;

  static inline Vector3 cross(const Vector3 &u, const Vector3 &v);
};


template<class Real>
inline Real & Vector3<Real>::operator [] (int i) {
  // no error checking done for the sake of speed
  return m[i];
}


template<class Real>
inline const Real & Vector3<Real>::operator [] (int i) const {
  // no error checking done for the sake of speed
  return m[i];
}


// This overloaded + operator performs vector addition.  The Vector3 object is added
// to another Vector3 object, v, and the result is returned.
template<class Real>
inline Vector3<Real> Vector3<Real>::operator + (const Vector3<Real> &v) const {
   return Vector3(m[0] + v.m[0], m[1] + v.m[1], m[2] + v.m[2]);
}


// This overloaded - operator performs vector subtraction.  v is subtracted from the Vector3
// object, and the result is returned.
template<class Real>
inline Vector3<Real> Vector3<Real>::operator - (const Vector3<Real> &v) const {
  return Vector3<Real>(m[0] - v.m[0], m[1] - v.m[1], m[2] - v.m[2]);
}


// vector addition, adds v to the Vector3 object and stores the result in
// itself.
template<class Real>
inline const Vector3<Real>& Vector3<Real>::operator += (const Vector3<Real> &v) {
  m[0] += v.m[0];
  m[1] += v.m[1];
  m[2] += v.m[2];

  return (*this);
}


// vector subtraction, subtracts v from the Vector3 object and stores the result in
// itself.
template<class Real>
inline const Vector3<Real>& Vector3<Real>::operator -= (const Vector3<Real> &v) {
  m[0] -= v.m[0];
  m[1] -= v.m[1];
  m[2] -= v.m[2];

  return (*this);
}


// overloaded * operator for scalar multiplication, multiplies this Vector3 object by
// the scalar and returns the resulting Vector3.
template<class Real>
inline Vector3<Real> Vector3<Real>::operator * (Real s) const {
  return Vector3<Real>(m[0]*s, m[1]*s, m[2]*s);
}


// overloaded / operator for scalar division, divides this Vector3 object by
// the scalar and returns the resulting Vector3.
template<class Real>
inline Vector3<Real> Vector3<Real>::operator / (Real s) const {
  return Vector3<Real>(m[0]/s, m[1]/s, m[2]/s);
}


// another overloaded * for scalar multiplication, but with vector on right side
template<class Real>
inline Vector3<Real> operator * (Real s, const Vector3<Real> &v) {
  return Vector3<Real>(v.m[0]*s, v.m[1]*s, v.m[2]*s);
}


// multiplies Vector3 by a scalar, and stores the result in itself.
template<class Real>
inline const Vector3<Real>& Vector3<Real>::operator *= (Real s) {
  m[0] *= s;
  m[1] *= s;
  m[2] *= s;

  return (*this);
}


// divides Vector3 by a scalar, and stores the result in itself
template<class Real>
inline const Vector3<Real>& Vector3<Real>::operator /= (Real s) {
  m[0] /= s;
  m[1] /= s;
  m[2] /= s;

  return (*this);
}


// This overloaded * operator returns the dot product of the vector, and another vector v.
template<class Real>
inline Real Vector3<Real>::operator *(const Vector3<Real> &v) const {
  return m[0]*v.m[0] + m[1]*v.m[1] + m[2]*v.m[2];
}


// This function computes the cross product of two vectors and returns the result
template<class Real>
inline Vector3<Real> Vector3<Real>::cross(const Vector3<Real> &u, const Vector3<Real> &v) {
  Vector3 result;
  result.m[0] = u.m[1]*v.m[2] - u.m[2]*v.m[1];
  result.m[1] = u.m[2]*v.m[0] - u.m[0]*v.m[2];
  result.m[2] = u.m[0]*v.m[1] - u.m[1]*v.m[0];
  return result;
}


// returns the magnitude of the vector
template<class Real>
inline Real Vector3<Real>::mag() const {
  return rsqrt(m[0]*m[0] + m[1]*m[1] + m[2]*m[2]);
}


// normalizes the vector
template<class Real>
inline void Vector3<Real>::normalize() {
  Real len = mag();

  m[0] /= len;
  m[1] /= len;
  m[2] /= len;
}


// reflect vector about vector n and return result
template<class Real>
inline Vector3<Real> Vector3<Real>::reflect(const Vector3<Real> &n) const {
  return n*Real(2.0)*( n*(*this) ) - (*this);
}


// stores vector reflected about n in 'result'
template<class Real>
inline void Vector3<Real>::reflect(const Vector3<Real> &n, Vector3 &result) {
  result = n*Real(2.0)*( n*(*this) ) - (*this);
}


// returns the squared magnitude of the vector (this is faster than mag(), since it
// does not involve the computation of a square root).  useful in calculations where
// magnitude can be left squared.
template<class Real>
inline Real Vector3<Real>::sqrMag() const {
  return m[0]*m[0] + m[1]*m[1] + m[2]*m[2];
}


// assign the values of the vector's components to x, y, and z
template<class Real>
inline void Vector3<Real>::set(Real x, Real y, Real z) {
  m[0] = x;
  m[1] = y;
  m[2] = z;
}


// overloaded stream insertion operator for Vector3
template<class Real>
std::ostream& operator << (std::ostream &out, const Vector3<Real> &v) {
  out << "<" << v.m[0] << ", " << v.m[1] << ", " << v.m[2] << ">";
  return out;
}


} // namespace Math3D


#endif
