#ifndef VECTOR3D_H_26B75ED0
#define VECTOR3D_H_26B75ED0

#include <ostream>
#include <cmath>

class Vector3D {
  public:
    double x, y, z;

    Vector3D() : x(0.0), y(0.0), z(0.0) {}
    Vector3D(double x, double y, double z) : x(x), y(y), z(z) {}
    Vector3D(double c) : x(c), y(c), z(c) {}
    Vector3D(const Vector3D& v) : x(v.x), y(v.y), z(v.z) {}

    inline double& operator[] (const int& index)
    {
      // Yuck
      return ( &x )[ index ];
    }

    inline const double& operator[] (const int& index) const
    {
      return ( &x )[ index ];
    }

    inline bool operator==(const Vector3D& v) const
    {
      return v.x == x && v.y == y && v.z == z;
    }

    inline Vector3D operator-() const
    {
      return Vector3D( -x, -y, -z );
    }

    inline Vector3D operator+(const Vector3D& v) const
    {
      return Vector3D( x + v.x, y + v.y, z + v.z );
    }

    inline Vector3D operator-(const Vector3D& v) const
    {
      return Vector3D( x - v.x, y - v.y, z - v.z );
    }

    inline Vector3D operator*(const double& c) const
    {
      return Vector3D( x * c, y * c, z * c );
    }

    inline Vector3D operator/(const double& c) const
    {
      const double rc = 1.0/c;
      return Vector3D( rc * x, rc * y, rc * z );
    }

    inline void operator+=(const Vector3D& v)
    {
      x += v.x;
      y += v.y;
      z += v.z;
    }

    inline void operator-=(const Vector3D& v)
    {
      x -= v.x;
      y -= v.y;
      z -= v.z;
    }

    inline void operator*=(const double& c)
    {
      x *= c;
      y *= c;
      z *= c;
    }

    inline void operator/=(const double& c)
    {
      (*this) *= ( 1./c );
    }

    inline double norm() const
    {
      return sqrt( x*x + y*y + z*z );
    }

    inline double norm2() const
    {
      return x*x + y*y + z*z;
    }

    inline Vector3D unit() const
    {
      double rNorm = 1. / sqrt( x*x + y*y + z*z );
      return Vector3D( rNorm*x, rNorm*y, rNorm*z );
    }

    inline void normalize()
    {
      (*this) /= norm();
    }

};

// left scalar multiplication
inline Vector3D operator* (const double& c, const Vector3D& v)
{
  return Vector3D( c * v.x, c * v.y, c * v.z );
}

// dot product (a.k.a. inner or scalar product)
inline double dot(const Vector3D& u, const Vector3D& v)
{
  return u.x*v.x + u.y*v.y + u.z*v.z ;
}

// cross product
inline Vector3D cross(const Vector3D& u, const Vector3D& v)
{
  return Vector3D( u.y*v.z - u.z*v.y,
                   u.z*v.x - u.x*v.z,
                   u.x*v.y - u.y*v.x );
}

// prints components
std::ostream& operator<<(std::ostream& os, const Vector3D& v);

#endif /* end of include guard: VECTOR3D_H_26B75ED0 */