#ifndef VECTOR_H
#define VECTOR_H VECTOR_H

#include <math.h>                   /* sqrt() */

#include <iostream>

class Vector
{
    public:
        Vector()
        {
            m_array[0] = 0.0;
            m_array[1] = 0.0;
            m_array[2] = 0.0;
        }

        Vector(double x, double y, double z)
        {
            m_array[0] = x;
            m_array[1] = y;
            m_array[2] = z;
        }

        Vector & normalize()
        {
            double length = mag();
            m_array[0] /= length;
            m_array[1] /= length;
            m_array[2] /= length;
            return *this;
        }

        double mag() const
        {
            return sqrt(m_array[0] * m_array[0]
                      + m_array[1] * m_array[1]
                      + m_array[2] * m_array[2]);
        }

        double mag2() const
        {
            return m_array[0] * m_array[0]
                 + m_array[1] * m_array[1]
                 + m_array[2] * m_array[2];
        }

        double dist_to(const Vector & other) const
        {
            return (other - *this).mag();
        }

        Vector proj(const Vector & target) const
        {
            Vector target_normalized = target;
            target_normalized.normalize();
            return target_normalized * ((*this) % target_normalized);
        }

        Vector reflect(const Vector & target) const
        {
            return (*this) - target * (2 * dot(target));
        }

        /*
         * from http://www.flipcode.com/archives/
         *   Reflections_and_Refraction_in_Raytracing.shtml
         * target: normal vector of surface
         * n1: refraction index of object we're coming from
         * n2: refraction index of object we're going into
         */
        Vector refract(const Vector & target, double n1, double n2) const
        {
            const double n = n1 / n2;
            const double cosI = dot(target);
            const double sinT2 = n * n * (1.0 - cosI * cosI);
            if (sinT2 > 1.0)
                return Vector(0.0, 0.0, 0.0);
            return (*this) * n - target * (n * cosI + sqrt(1.0 - sinT2));
        }

        Vector getPerpendicular() const
        {
            Vector t = *this;
            t.normalize();
            Vector p = t * Vector(0, 0, 1);
            if (p.mag() <= 0.1)
            {
                p = t * Vector(1, 0, 0);
            }
            return p;
        }

        Vector mult(const Vector & v2) const
        {
            return Vector(
                    m_array[0] * v2.m_array[0],
                    m_array[1] * v2.m_array[1],
                    m_array[2] * v2.m_array[2]);
        }

        Vector div(const Vector & v2) const
        {
            return Vector(
                    m_array[0] / v2.m_array[0],
                    m_array[1] / v2.m_array[1],
                    m_array[2] / v2.m_array[2]);
        }

        Vector operator-() const
        {
            return Vector(
                    -m_array[0],
                    -m_array[1],
                    -m_array[2]);
        }

        /* Compute the dot-product of two vectors */
        double dot(const Vector & v2) const
        {
            return m_array[0] * v2.m_array[0]
                + m_array[1] * v2.m_array[1]
                + m_array[2] * v2.m_array[2];
        }
        double operator%(const Vector & v2) const { return dot(v2); }

        /* Compute the cross-product of two vectors */
        Vector cross(const Vector & v2) const
        {
            return Vector(
                    m_array[1] * v2.m_array[2] - m_array[2] * v2.m_array[1],
                    m_array[2] * v2.m_array[0] - m_array[0] * v2.m_array[2],
                    m_array[0] * v2.m_array[1] - m_array[1] * v2.m_array[0]);
        }
        Vector operator*(const Vector & v2) const { return cross(v2); }

        Vector operator+(const Vector & v2) const
        {
            return Vector(
                    m_array[0] + v2.m_array[0],
                    m_array[1] + v2.m_array[1],
                    m_array[2] + v2.m_array[2]);
        }

        Vector operator-(const Vector & v2) const
        {
            return Vector(
                    m_array[0] - v2.m_array[0],
                    m_array[1] - v2.m_array[1],
                    m_array[2] - v2.m_array[2]);
        }

        Vector operator*(double scale) const
        {
            return Vector(
                    m_array[0] * scale,
                    m_array[1] * scale,
                    m_array[2] * scale);
        }

        Vector operator/(double scale) const
        {
            return Vector(
                    m_array[0] / scale,
                    m_array[0] / scale,
                    m_array[0] / scale);
        }

        Vector & operator+=(const Vector & v2)
        {
            m_array[0] += v2.m_array[0];
            m_array[1] += v2.m_array[1];
            m_array[2] += v2.m_array[2];
            return *this;
        }

        Vector & operator-=(const Vector & v2)
        {
            m_array[0] -= v2.m_array[0];
            m_array[1] -= v2.m_array[1];
            m_array[2] -= v2.m_array[2];
            return *this;
        }

        double & operator[](int idx) { return m_array[idx]; }
        double operator[](int idx) const { return m_array[idx]; }

        static Vector randomVector();

        static const Vector X;
        static const Vector Y;
        static const Vector Z;

    protected:
        double  m_array[3];
};

std::ostream & operator<<(std::ostream & out, const Vector & v);
static inline Vector operator*(double d, const Vector & v) { return v * d; }
static inline Vector operator/(double d, const Vector & v) { return v / d; }

#endif