/*
 * Copyright Staffan Gimåker 2006-2008.
 *
 * Distributed under the Boost Software License, Version 1.0.
 * (See accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 */

#ifndef __PEEKABOT_VECTOR_HH
#define __PEEKABOT_VECTOR_HH

/**
 * \file 
 *
 * \brief Three- and four-dimensional vector implementations.
 */

#include "Types.hh"
#include <cmath>
#include <ostream>


namespace peekabot
{


    template<class T = real_unit>
    class Vector3
    {
    public:
        typedef T UnitType;

        
        /**
         * \brief The vector's components.
         */
        T m_c[3];



        /**
         * \brief Default constructor.
         *
         * The coordinates of the constructed vector are undefined.
         */
        Vector3() throw() {}

        
        /**
         * \brief Copy constructor (deep copy).
         *
         * \param v The vector which to construct a copy of.
         */
        template<class U>
        Vector3(const Vector3<U> &v) throw()
        {
            // Consider memcpy...
            m_c[0] = v.m_c[0];
            m_c[1] = v.m_c[1]; 
            m_c[2] = v.m_c[2];
        }


        /**
         * \brief Convenience constructor for constructing a vector from 
         * three axis coordinates.
         */
        Vector3(T x, T y, T z) throw()
        {
            m_c[0] = x;
            m_c[1] = y;
            m_c[2] = z;
        }


        /**
         * \brief Vector component access method.
         *
         * Does not perform bounds checking! Be cautious, or expect a segfault.
         * Components are indexed from 0 and up.
         *
         * \param i The component to retrieve, must be in the range [0,2].
         * \return A (modifiable) reference to the <em>i</em>:th component.
         */
        T &operator()(int component) throw()
        {
            return m_c[component];
        }


        /**
         * \brief Vector component access method.
         *
         * Does not perform bounds checking! Be cautious, or expect a segfault.
         * Components are indexed from 0 and up.
         *
         * \param i The component to retrieve, must be in the range [0,2].
         * \return A copy of the <em>i</em>:th component.
         */
        const T &operator()(int component) const throw()
        {
            return m_c[component];
        }


        /**
         * \brief Calculate the dot product between two vectors.
         *
         * \param v The second vector in \f$<w,v>\f$.
         * \return The dot product between the vector and \c v.
         */
        template<class U>
        T operator*(const Vector3<U> &v) const throw()
        {
            return m_c[0] * v.m_c[0] + 
                   m_c[1] * v.m_c[1] + 
                   m_c[2] * v.m_c[2];
        }

        
        template<class U>
        Vector3<T> operator+(const Vector3<U> &v) const throw()
        {
            return Vector3<T>(m_c[0] + v.m_c[0], 
                              m_c[1] + v.m_c[1],
                              m_c[2] + v.m_c[2]);
        }


        template<class U>
        Vector3<T> operator-(const Vector3<U> &v) const throw()
        {
            return Vector3<T>(m_c[0] - v.m_c[0], 
                              m_c[1] - v.m_c[1],
                              m_c[2] - v.m_c[2]);
        }


        Vector3<T> operator*(T scalar) const throw()
        {
            return Vector3<T>(scalar*m_c[0], scalar*m_c[1], scalar*m_c[2]);
        }


        Vector3<T> operator/(T scalar) const throw()
        {
            return Vector3<T>(m_c[0]/scalar, m_c[1]/scalar, m_c[2]/scalar);
        }


        template<class U>
        Vector3<T> &operator+=(const Vector3<U> &v) throw()
        {
            m_c[0] += v.m_c[0];
            m_c[1] += v.m_c[1];
            m_c[2] += v.m_c[2];
            return *this;
        }


        template<class U>
        Vector3<T> &operator-=(const Vector3<U> &v) throw()
        {
            m_c[0] -= v.m_c[0];
            m_c[1] -= v.m_c[1];
            m_c[2] -= v.m_c[2];
            return *this;
        }


        Vector3<T> &operator*=(T scalar) throw()
        {
            m_c[0] *= scalar;
            m_c[1] *= scalar;
            m_c[2] *= scalar;
            return *this;
        }

        // TODO: fix the throw type, or make up some precondition..
        // whatever.
        Vector3<T> &operator/=(T scalar)
        {
            return *this *= 1/scalar;
        }


        Vector3<T> &negate() throw()
        {
            m_c[0] = -m_c[0];
            m_c[1] = -m_c[1];
            m_c[2] = -m_c[2];
            return *this;
        }

        
        Vector3<T> operator-() const throw()
        {
            return Vector3<T>(-m_c[0], -m_c[1], -m_c[2]);
        }


        T norm() const throw()
        {
            return static_cast<T>(sqrt(*this * *this));
        }


        T norm_squared() const throw()
        {
            return (*this) * (*this);
        }


        // if it's the 0 vector, well it won't be normalized...
        Vector3<T> &normalize() throw()
        {
            *this /= norm();
            return *this;
        }

        /**
         * \brief Calculate the square distance between two vectors.
         *
         * \param v The vector which to calculate the distance to.
         * \return The squared distance from \c *this to \c v.
         */
        template<class U>
        T dist_squared(const Vector3<U> &v) const throw()
        {
            return (*this - v) * (*this - v);
        }

        /**
         * \brief Calculate the distance between two vectors.
         *
         * \param v The vector which to calculate the distance to.
         * \return The distance from \c *this to \c v.
         */
        template<class U>
        T dist(const Vector3<U> &v) const throw()
        {
            return (*this - v).norm();
        }

        /**
         * \brief Calculate the cross product of the vector and the 
         * given parameter.
         *
         * \return <tt>*this x v</tt>.
         */
        Vector3<T> cross_product(const Vector3<T> &v) const throw()
        {
            return Vector3<T>(m_c[1]*v.m_c[2] - m_c[2]*v.m_c[1],
                              m_c[2]*v.m_c[0] - m_c[0]*v.m_c[2],
                              m_c[0]*v.m_c[1] - m_c[1]*v.m_c[0]);
        }


        template<class U>
        bool operator==(const Vector3<U> &v) const throw()
        {
            return (v.m_c[0] == m_c[0] &&
                    v.m_c[1] == m_c[1] &&
                    v.m_c[2] == m_c[2]);
        }


        template<class U, class S>
        bool equals(const Vector3<U> &v, S limit) const throw()
        {
            return ((S)fabs(v.m_c[0]-m_c[0]) <= limit &&
                    (S)fabs(v.m_c[1]-m_c[1]) <= limit &&
                    (S)fabs(v.m_c[2]-m_c[2]) <= limit);
        }
        
    };



    /**
     * \brief A four dimensional vector class.
     */
    template<class T = real_unit>
    class Vector4
    {
    public:
        typedef T UnitType;

        /**
         * \brief The vector's components.
         */
        T m_c[4];


        operator Vector3<T>() const throw()
        {
            return Vector3<T>(m_c[0], m_c[1], m_c[2]);
        }


        /**
         * \brief Default constructor.
         *
         * The coordinates of the constructed vector are undefined.
         */
        Vector4() throw() {}

        
        /**
         * \brief Copy constructor (deep copy).
         *
         * \param v The vector which to construct a copy of.
         */
        template<class U>
        Vector4(const Vector4<U> &v) throw()
        {
            // Consider memcpy...
            m_c[0] = v.m_c[0];
            m_c[1] = v.m_c[1]; 
            m_c[2] = v.m_c[2];
            m_c[3] = v.m_c[3];
        }


        /**
         * \brief Convenience constructor for constructing a vector from 
         * four axis coordinates.
         */
        Vector4(T x, T y, T z, T w) throw()
        {
            m_c[0] = x;
            m_c[1] = y; 
            m_c[2] = z;
            m_c[3] = w;
        }


        /**
         * \brief Convenience constructor, for constructing a vector from 
         * a 3D vector and an extra w-coordinate.
         */
        template<class U>
        Vector4(const Vector3<U> &v, T w) throw()
        {
            // Consider memcpy...
            m_c[0] = v.m_c[0];
            m_c[1] = v.m_c[1]; 
            m_c[2] = v.m_c[2];
            m_c[3] = w;
        }


        /**
         * \brief Vector component access method.
         *
         * Does not perform bounds checking! Be cautious, or expect a segfault.
         * Components are indexed from 0 and up.
         *
         * \param i The component to retrieve, must be in the range [0,3].
         * \return A (modifiable) reference to the <em>i</em>:th component.
         */
        T &operator()(int component) throw()
        {
            return m_c[component];
        }


        /**
         * \brief Vector component access method.
         *
         * Does not perform bounds checking! Be cautious, or expect a segfault.
         * Components are indexed from 0 and up.
         *
         * \param i The component to retrieve, must be in the range [0,3].
         * \return A copy of the <em>i</em>:th component.
         */
        const T &operator()(int component) const throw()
        {
            return m_c[component];
        }

        
        /**
         * \brief Calculate the dot product between two vectors.
         *
         * \param v The second vector in \f$<w,v>\f$.
         * \return The dot product between the vector and \c v.
         */
        template<class U>
        T operator*(const Vector4<U> &v) const throw()
        {
            return m_c[0] * v.m_c[0] + 
                   m_c[1] * v.m_c[1] + 
                   m_c[2] * v.m_c[2] + 
                   m_c[3] * v.m_c[3];
        }

        
        template<class U>
        Vector4<T> operator+(const Vector4<U> &v) const throw()
        {
            return Vector4<T>(m_c[0] + v.m_c[0], 
                              m_c[1] + v.m_c[1],
                              m_c[2] + v.m_c[2],
                              m_c[3] + v.m_c[3]);
        }


        template<class U>
        Vector4<T> operator-(const Vector4<U> &v) const throw()
        {
            return Vector4<T>(m_c[0] - v.m_c[0], 
                              m_c[1] - v.m_c[1],
                              m_c[2] - v.m_c[2],
                              m_c[3] - v.m_c[3]);
        }


        Vector4<T> operator*(T scalar) const throw()
        {
            return Vector4<T>(scalar*m_c[0], scalar*m_c[1], 
                              scalar*m_c[2], scalar*m_c[3]);
        }


        Vector4<T> operator/(T scalar) const throw()
        {
            return Vector4<T>(m_c[0]/scalar, m_c[1]/scalar, 
                              m_c[2]/scalar, m_c[3]/scalar);
        }


        template<class U>
        Vector4<T> &operator+=(const Vector4<U> &v) throw()
        {
            m_c[0] += v.m_c[0];
            m_c[1] += v.m_c[1];
            m_c[2] += v.m_c[2];
            m_c[3] += v.m_c[3];
            return *this;
        }


        template<class U>
        Vector4<T> &operator-=(const Vector4<U> &v) throw()
        {
            m_c[0] -= v.m_c[0];
            m_c[1] -= v.m_c[1];
            m_c[2] -= v.m_c[2];
            m_c[3] -= v.m_c[3];
            return *this;
        }


        Vector4<T> &operator*=(T scalar) throw()
        {
            m_c[0] *= scalar;
            m_c[1] *= scalar;
            m_c[2] *= scalar;
            m_c[3] *= scalar;
            return *this;
        }

        // TODO: fix the throw type, or make up some precondition..
        // whatever.
        Vector4<T> &operator/=(T scalar) throw()
        {
            return *this *= 1/scalar;
        }


        Vector4<T> &negate() throw()
        {
            m_c[0] = -m_c[0];
            m_c[1] = -m_c[1];
            m_c[2] = -m_c[2];
            m_c[3] = -m_c[3];
            return *this;
        }

        
        Vector4<T> operator-() const throw()
        {
            return Vector4<T>(-m_c[0], -m_c[1], -m_c[2], -m_c[3]);
        }   

        T norm() const throw()
        {
            return static_cast<T>(sqrt(*this * *this));
        }


        T norm_squared() const throw()
        {
            return (*this) * (*this);
        }


        // if it's the 0 vector, well it won't be normalized...
        Vector4<T> &normalize() throw()
        {
            *this /= norm();
            return *this;
        }

        /**
         * \brief Calculate the square distance between two vectors.
         *
         * \param v The vector which to calculate the distance to.
         * \return The squared distance from \c *this to \c v.
         */
        template<class U>
        T dist_squared(const Vector4<U> &v) const throw()
        {
            return (*this - v) * (*this - v);
        }

        /**
         * \brief Calculate the distance between two vectors.
         *
         * \param v The vector which to calculate the distance to.
         * \return The distance from \c *this to \c v.
         */
        template<class U>
        T dist(const Vector4<U> &v) const throw()
        {
            return (*this - v).norm();
        }


        template<class U>
        bool operator==(const Vector4<U> &v) const throw()
        {
            return (v.m_c[0] == m_c[0] &&
                    v.m_c[1] == m_c[1] &&
                    v.m_c[2] == m_c[2] &&
                    v.m_c[3] == m_c[3]);
        }
    };


    /**
     * \brief A two dimensional vector class.
     */
    template<class T = real_unit>
    class Vector2
    {
    public:
        typedef T UnitType;

        /**
         * \brief The vector's components.
         */
        T m_c[2];


        /**
         * \brief Default constructor.
         *
         * The coordinates of the constructed vector are undefined.
         */
        Vector2() throw() {}

        
        /**
         * \brief Copy constructor (deep copy).
         *
         * \param v The vector which to construct a copy of.
         */
        template<class U>
        Vector2(const Vector2<U> &v) throw()
        {
            m_c[0] = v.m_c[0];
            m_c[1] = v.m_c[1]; 
        }


        /**
         * \brief Convenience constructor for constructing a vector from 
         * two axis coordinates.
         */
        Vector2(T x, T y) throw()
        {
            m_c[0] = x;
            m_c[1] = y;
        }


        /**
         * \brief Vector component access method.
         *
         * Does not perform bounds checking! Be cautious, or expect a segfault.
         * Components are indexed from 0 and up.
         *
         * \param i The component to retrieve, must be in the range [0,1].
         * \return A (modifiable) reference to the <em>i</em>:th component.
         */
        T &operator()(int component) throw()
        {
            return m_c[component];
        }


        /**
         * \brief Vector component access method.
         *
         * Does not perform bounds checking! Be cautious, or expect a segfault.
         * Components are indexed from 0 and up.
         *
         * \param i The component to retrieve, must be in the range [0,1].
         * \return A copy of the <em>i</em>:th component.
         */
        const T &operator()(int component) const throw()
        {
            return m_c[component];
        }

        
        /**
         * \brief Calculate the dot product between two vectors.
         *
         * \param v The second vector in \f$<w,v>\f$.
         * \return The dot product between the vector and \c v.
         */
        template<class U>
        T operator*(const Vector2<U> &v) const throw()
        {
            return ( m_c[0] * v.m_c[0] + 
                     m_c[1] * v.m_c[1] );
        }

        
        template<class U>
        Vector2<T> operator+(const Vector2<U> &v) const throw()
        {
            return Vector2<T>(m_c[0] + v.m_c[0], 
                              m_c[1] + v.m_c[1]);
        }


        template<class U>
        Vector2<T> operator-(const Vector2<U> &v) const throw()
        {
            return Vector2<T>(m_c[0] - v.m_c[0], 
                              m_c[1] - v.m_c[1]);
        }


        Vector2<T> operator*(T scalar) const throw()
        {
            return Vector2<T>(scalar*m_c[0], scalar*m_c[1]);
        }


        Vector2<T> operator/(T scalar) const throw()
        {
            return Vector2<T>(m_c[0]/scalar, m_c[1]/scalar);
        }


        template<class U>
        Vector2<T> &operator+=(const Vector2<U> &v) throw()
        {
            m_c[0] += v.m_c[0];
            m_c[1] += v.m_c[1];
            return *this;
        }


        template<class U>
        Vector2<T> &operator-=(const Vector2<U> &v) throw()
        {
            m_c[0] -= v.m_c[0];
            m_c[1] -= v.m_c[1];
            return *this;
        }


        Vector2<T> &operator*=(T scalar) throw()
        {
            m_c[0] *= scalar;
            m_c[1] *= scalar;
            return *this;
        }

        // TODO: fix the throw type, or make up some precondition..
        // whatever.
        Vector2<T> &operator/=(T scalar) throw()
        {
            return *this *= 1/scalar;
        }


        Vector2<T> &negate() throw()
        {
            m_c[0] = -m_c[0];
            m_c[1] = -m_c[1];
            return *this;
        }

        
        Vector2<T> operator-() const throw()
        {
            return Vector2<T>(-m_c[0], -m_c[1]);
        }   

        T norm() const throw()
        {
            return static_cast<T>(sqrt(*this * *this));
        }


        T norm_squared() const throw()
        {
            return (*this) * (*this);
        }


        // if it's the 0 vector, well it won't be normalized...
        Vector2<T> &normalize() throw()
        {
            *this /= norm();
            return *this;
        }

        /**
         * \brief Calculate the square distance between two vectors.
         *
         * \param v The vector which to calculate the distance to.
         * \return The squared distance from \c *this to \c v.
         */
        template<class U>
        T dist_squared(const Vector2<U> &v) const throw()
        {
            return (*this - v) * (*this - v);
        }

        /**
         * \brief Calculate the distance between two vectors.
         *
         * \param v The vector which to calculate the distance to.
         * \return The distance from \c *this to \c v.
         */
        template<class U>
        T dist(const Vector2<U> &v) const throw()
        {
            return (*this - v).norm();
        }


        template<class U>
        bool operator==(const Vector2<U> &v) const throw()
        {
            return (v.m_c[0] == m_c[0] &&
                    v.m_c[1] == m_c[1]);
        }

        template<class U, class S>
        bool equals(const Vector2<U> &v, S limit) const throw()
        {
            return ((S)fabs(v.m_c[0]-m_c[0]) <= limit &&
                    (S)fabs(v.m_c[1]-m_c[1]) <= limit);
        }
    };


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

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

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



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


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

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




    typedef Vector3<real_unit> Vector3ru;
    typedef Vector3<real32> Vector3f;
    typedef Vector3<real64> Vector3d;
    /**
     * \brief The Vector type used internally in the render engine.
     */
    typedef Vector3<real_unit> Vector3r;

    typedef Vector4<real_unit> Vector4ru;
    typedef Vector4<real32> Vector4f;
    typedef Vector4<real64> Vector4d;
    /**
     * \brief The Vector type used internally in the render engine.
     */
    typedef Vector4<real_unit> Vector4r;

    typedef Vector2<float> Vector2f;


}


#endif // __PEEKABOT_VECTOR_HH
