// (c) Copyright 2010 Terry Golubiewski, all rights reserved.

#ifndef CARTESIAN_VECTOR_H
#define CARTESIAN_VECTOR_H
#pragma once

#include <boost/array.hpp>
#include <stdexcept>
#include <boost/throw_exception.hpp>
#include <std0x/common_type.h>
#include <cmath>
#include <algorithm>

namespace cartesian {

struct uninitialized_t { };
static const uninitialized_t uninitialized = { };

template<int Dim=3, class Rep = double>
class vector {
public:
  static const int dim = Dim;
  typedef Rep rep;

  boost::array<rep, dim> coord;

  vector(uninitialized_t) { }

  vector() { coord.fill(rep()); }

  vector(const vector& v) : coord(v.coord) { }
  vector& operator=(const vector& v) { coord = v.coord; return *this; }

  template<class U>
  vector(const vector<dim, U>& v) { coord = v.coord; }

  template<class U>
  vector& operator=(const vector<dim, U>& v) { coord = v.coord; return *this; }

  template<class InIter>
  vector(InIter first_, InIter last_) {
    using namespace std;
    if (distance(first_, last_) != dim) {
      std::out_of_range e("cartesian::vector<> invalid input iterators");
      boost::throw_exception(e);
    }
    copy(first_, last_, coord.begin());
  } // iterator initialization

  template<int N> const rep& get() const { return coord[N]; }
  template<int N>       rep& get()       { return coord[N]; }

  template<int N> void set(rep v) { coord[N] = v; }
  
  rep& operator[](int i) { return coord[i]; }
  const rep& operator[](int i) const { return coord[i]; }

  template<class U>
  vector& operator+=(const vector<dim, U>& rhs) {
    for (int i=0; i!=dim; ++i)
      coord[i] += rhs[i];
    return *this;
  }

  template<class U>
  vector& operator-=(const vector<dim, U>& rhs) {
    for (int i=0; i!=dim; ++i)
      coord[i] -= rhs[i];
    return *this;
  }

  vector& operator*=(const rep& rhs) {
    for (int i=0; i!=dim; ++i)
      coord[i] *= rhs;
    return *this;
  }

  vector& operator/=(const rep& rhs) {
    for (int i=0; i!=dim; ++i)
      coord[i] /= rhs;
    return *this;
  }

  friend vector operator*(vector lhs, const rep& rhs)
    { return lhs *= rhs; }

  friend vector operator/(vector lhs, const rep& rhs)
    { return lhs /= rhs; }

  friend rep abs_squared(const vector& v) {
    rep sum = v[0] * v[0];
    for (int i=1; i!=dim; ++i) {
      const rep& x = v[i];
      sum += x * x;
    }
    return sum;
  } // abs_squared

  friend rep abs(const vector& v) {
    using namespace std;
    return sqrt(abs_squared(v));
  }

}; // vector

template<int Dim, class T1, class T2> inline
vector<Dim, typename std0x::common_type<T1, T2>::type>
operator+(const vector<Dim, T1>&lhs, const vector<Dim, T2>& rhs) {
  vector<Dim, typename std0x::common_type<T1, T2>::type> rval(lhs);
  return rval += rhs;
}

template<int Dim, class T1, class T2> inline
vector<Dim, typename std0x::common_type<T1, T2>::type>
operator-(const vector<Dim, T1>&lhs, const vector<Dim, T2>& rhs) {
  vector<Dim, typename std0x::common_type<T1, T2>::type> rval(lhs);
  return rval -= rhs;
}

template<int Dim, class T1, class T2> inline
typename std0x::common_type<T1, T2>::type
operator*(const vector<Dim, T1>& lhs, const vector<Dim, T2>& rhs) {
  typedef typename std0x::common_type<T1, T2>::type return_type;
  return_type sum = lhs[0] * rhs[0];
  for (int i=1; i!=dim; ++i)
    sum += lhs[i] * rhs[i];
  return sum;
} // dot product

template<int Dim, class T1, class T2> inline
bool operator==(const vector<Dim, T1>& lhs, const vector<Dim, T2>& rhs) {
  for (int i=0; i!=dim; ++i) {
    if (lhs[i] != rhs[i])
      return false;
  }
  return true;
}

template<int Dim, class T1, class T2> inline
bool operator!=(const vector<Dim, T1>& lhs, const vector<Dim, T2>& rhs)
  { return !(lhs == rhs); }

// Is it a good idea to order vectors based on their magnitude?

template<int Dim, class T1, class T2> inline
bool operator<(const vector<Dim, T1>& lhs, const vector<Dim, T2>& rhs)
  { return abs_squared(lhs) < abs_squared(rhs); }

template<int Dim, class T1, class T2> inline
bool operator>(const vector<Dim, T1>& lhs, const vector<Dim, T2>& rhs)
  { return rhs < lhs; }

template<int Dim, class T1, class T2> inline
bool operator>=(const vector<Dim, T1>& lhs, const vector<Dim, T2>& rhs)
  { return !(lhs < rhs); }

template<int Dim, class T1, class T2> inline
bool operator<=(const vector<Dim, T1>& lhs, const vector<Dim, T2>& rhs)
  { return !(rhs < lhs); }

template<class Rep> inline
vector<1, Rep> make_vector(Rep x) {
  vector<1, Rep> v(uninitialized);
  v.set<0>(x);
  return v;
} // make_vector

template<class T1, class T2> inline
vector<2, typename std0x::common_type<T1, T2>::type> make_vector(T1 x, T2 y) {
  vector<2, typename std0x::common_type<T1, T2>::type> v(uninitialized);
  v.set<0>(x);
  v.set<1>(y);
  return v;
} // make_vector

template<class T1, class T2, class T3>
struct common_type3 {
  typedef typename std0x::common_type<T1,
                                    typename std0x::common_type<T2, T3>::type
                                    >::type type;
}; // common_type3

template<class T1, class T2, class T3> inline
vector<3, typename common_type3<T1, T2, T3>::type>
make_vector(T1 x, T2 y, T3 z) {
  vector<3, typename common_type3<T1, T2, T3>::type> v(uninitialized);
  v.set<0>(x);
  v.set<1>(y);
  v.set<2>(z);
  return v;
} // make_vector

template<class T1, class T2> inline
vector<3, typename std0x::common_type<T1, T2>::type>
cross_product(const vector<3, T1>& a, const vector<3, T2>& b) {
  return make_vector(a[1] * b[2] - a[2] * b[1],
                     a[2] * b[0] - a[0] * b[2],
                     a[0] * b[1] - a[1] * b[0]);
} // cross_product

class vector2d: public vector<2, double> {
public:
  vector2d(double dx, double dy)
    : vector<2, double>(make_vector(dx, dy)) { }

  template<class U>
  vector2d(vector<2, U> v) : vector<2, double>(v) { }

  double& x() { return get<0>(); }
  double& y() { return get<1>(); }

  const double& x() const { return get<0>(); }
  const double& y() const { return get<1>(); }
}; // vector2d

class vector3d: public vector<3, double> {
public:
  vector3d(double dx, double dy, double dz=0)
    : vector<3, double>(make_vector(dx, dy, dz)) { }

  template<class U>
  vector3d(vector<3, U> v) : vector<3, double>(v) { }

  double& x() { return get<0>(); }
  double& y() { return get<1>(); }
  double& z() { return get<2>(); }

  const double& x() const { return get<0>(); }
  const double& y() const { return get<1>(); }
  const double& z() const { return get<2>(); }
}; // vector3d

} // cartesian

#endif