#ifndef   geo_Vectors_H__
#define   geo_Vectors_H__

#include "../math/utility.h"
#include "../math/Vector.h"
#include "../math/Matrix.h"

#include <cmath>

namespace meow {

/*!
 * @brief 2D's vector
 *
 * @author cat_leopard
 */
template<class Scalar>
class Vector2D {
private:
  Scalar x_, y_;
public:
  //! @brief consturctor (0, 0)
  Vector2D(): x_(0), y_(0) {
  }

  //! @brief consturctor (from another Vector2D)
  Vector2D(Vector2D const& v): x_(v.x_), y_(v.y_) {
  }

  //! @brief constructor (s, s)
  Vector2D(Scalar const& s): x_(s), y_(s) {
  }

  //! @brief constructor (sx, sy)
  Vector2D(Scalar const& sx, Scalar const& sy): x_(sx), y_(sy) {
  }

  //! @brief constructor (from another Vector)
  Vector2D(Vector<Scalar> const& v): x_(v(0)), y_(v(1)) {
  }

  //! @brief constructor (from another Vector, i-th)
  Vector2D(Vector<Scalar> const& v, size_t i): x_(v(i)), y_(v(i + 1)) {
  }

  //! @brief destructor
  ~Vector2D() {
  }

  //! @brief copy
  Vector2D& copyFrom(Vector2D const& v) {
    return xy(v.x(), v.y());
  }

  size_t dimension() const { return 2; } //!< @brief return 2

  Scalar const& x() const { return x_; } //!< @brief access x
  Scalar const& y() const { return y_; } //!< @brief access y

  Scalar& xGet() { return x_; } //!< @brief access x with non constant reference
  Scalar& yGet() { return y_; } //!< @brief access y with non constant reference

  Scalar const& x(Scalar const& s) { x_ = s; return x(); } //!< @brief modify x
  Scalar const& y(Scalar const& s) { y_ = s; return y(); } //!< @brief modify y

  //! @brief modify x and y
  Vector2D& xy(Scalar const& sx, Scalar const& sy){
    x(sx);
    y(sy);
    return *this;
  }

  //! @brief access the \c i -th scalar (0 => x, 1 => y)
  Scalar const& scalar(size_t i) const {
    static Scalar tmp(0);
    return (i == 0 ? x() : (i == 1 ? y() : tmp));
  }

  //! @brief modivy the \c i -th scalar (0 => x, 1 => y)
  Scalar const& scalar(size_t i, Scalar const& s) {
    return (i == 0 ? x(s) : (i == 1 ? y(s) : s));
  }

  //! @brief return \a +(*this)
  Vector2D positive() const {
    return *this;
  }

  //! @brief return \a -(*this)
  Vector2D negative() const {
    return Vector2D(-x(), -y());
  }

  //! @brief return \a count-clockwise \a rotate \a 90 \a degree of itself
  Vector2D right()const{
    return Vector2D(-y(), x());
  }

  //! @brief return \a (*this)+v
  Vector2D add(Vector2D const& v) const {
    return Vector2D(x() + v.x(), y() + v.y());
  }

  //! @brief Let itself add v
  Vector2D& added(Vector2D const& v) {
    return xy(x() + v.x(), y() + v.y());
  }

  //! @brief return \a (*this)-v
  Vector2D sub(Vector2D const& v) const {
    return Vector2D(x() - v.x(), y() - v.y());
  }

  //! @brief Let itself substract v
  Vector2D& subed(Vector2D const& v) {
    return xy(x() - v.x(), y() - v.y());
  }

  //! @brief return \a (*this)*s , where s is a scalar
  Vector2D mul(Scalar const& s) const {
    return Vector2D(x() * s, y() * s);
  }

  //! @brief Let itself mulitple s
  Vector2D& muled(Scalar const& s) {
    return xy(x() * s, y() * s);
  }

  //! @brief return \a (*this)/s , where s is a scalar
  Vector2D div(Scalar const& s) const {
    return Vector2D(x() / s, y() / s);
  }

  //! @brief Let itself divide s
  Vector2D& dived(Scalar const& s) {
    return xy(x() / s, y() / s);
  }

  //! @brief same as dot(v)
  Scalar mul(Vector2D const& v) const {
    return dot(v);
  }

  //! @brief dot
  Scalar dot(Vector2D const& v) const {
    return x() * v.x() + y() * v.y();
  }

  //! @brief cross
  Scalar cross(Vector2D const& v) const {
    return x() * v.y() - y() * v.x();
  }

  //! @brief sqrt of length2
  Scalar length() const {
    return Scalar(sqrt(double(length2())));
  }

  //! @brief same as \a dot(*this)
  Scalar length2() const {
    return dot(*this);
  }

  //! @brief return normalize form of itself
  Vector2D normalize() const {
    return div(length());
  }

  //! @brief normalize itself
  Vector2D& normalized() {
    return dived(length());
  }

  //! @brief return rotate \a theta degree of itself
  Vector2D rotate(Scalar const& theta) const {
    Scalar cs(cos(-double(theta)));
    Scalar sn(sin(-double(theta)));
    Vector2D<Scalar> new_x(cs, sn);
    return Vector2D(new_x.dot(*this), new_x.cross(*this));
  }

  //! @brief Let itself rotate \a theta degree
  Vector2D& rotated(Scalar const& theta) {
    return copyFrom(rotate(theta));
  }

  //! @brief return reflect from given vector \a v
  Vector2D reflect(Vector2D const& v) const {
    return v.mul(v.dot(*this) * 2 / v.length2()).sub(*this);
  }

  //! @brief reflect itself given vector \a v
  Vector2D& reflected(Vector2D const& v) {
    return copyFrom(reflecte(v));
  }

  //! @brief return a 2x1 matrix form of itself
  Matrix<Scalar> matrix() const {
    static Matrix<Scalar> ret(2, 1, Scalar(0));
    ret(0, 0, x());
    ret(1, 0, y());
    return ret;
  }

  //! @brief return a 3x1 matrix form of itself
  Matrix<Scalar> matrix(Scalar const& homo) const {
    static Matrix<Scalar> ret(3, 1, Scalar(0));
    ret(0, 0, x());
    ret(1, 0, y());
    ret(2, 0, homo);
    return ret;
  }

  Scalar const& operator()(size_t n) const {
    return (n == 0 ? x() : y());
  }

  Vector2D& operator()(Scalar const& sx, Scalar const& sy) {
    return xy(sx, sy);
  }

  Vector2D operator+() const { return positive(); }
  Vector2D operator-() const { return negative(); }
  Vector2D operator~() const { return right   (); }

  Vector2D operator+(Vector2D const& v) const { return add(v); }
  Vector2D operator-(Vector2D const& v) const { return sub(v); }
  Vector2D operator*(Scalar   const& s) const { return mul(s); }
  Vector2D operator/(Scalar   const& s) const { return div(s); }
  Scalar   operator*(Vector2D const& v) const { return mul(v); }

  Vector2D& operator =(Vector2D const& v) { return copyFrom(v); }
  Vector2D& operator+=(Vector2D const& v) { return added(v); }
  Vector2D& operator-=(Vector2D const& v) { return subed(v); }
  Vector2D& operator*=(Scalar   const& s) { return muled(s); }
  Vector2D& operator/=(Scalar   const& s) { return dived(s); }
};

/*!
 * @brief 3D's vector
 *
 * @author cat_leopard
 */
template<class Scalar>
class Vector3D{
private:
  Scalar x_, y_, z_;
public:
  //! @brief consturctor (0, 0)
  Vector3D(): x_(0), y_(0), z_(0) {
  }

  //! @brief consturctor (from another Vector3D)
  Vector3D(Vector3D const& v): x_(v.x_), y_(v.y_), z_(v.z_) {
  }

  //! @brief constructor (s, s)
  Vector3D(Scalar const& s): x_(s), y_(s), z_(s) {
  }

  //! @brief constructor (sx, sy)
  Vector3D(Scalar const& sx,
           Scalar const& sy,
           Scalar const& sz): x_(sx), y_(sy), z_(sz) {
  }

  //! @brief constructor (from another Vector)
  Vector3D(Vector<Scalar> const& v): x_(v(0)), y_(v(1)), z_(v(2)) {
  }

  //! @brief constructor (from another Vector, i-th)
  Vector3D(Vector<Scalar> const& v, size_t i): x_(v(i)), y_(v(i+1)), z_(v(i+2)){
  }

  //! @brief destructor
  ~Vector3D() {
  }

  //! @brief copy
  Vector3D& copyFrom(Vector3D const& v) {
    return xyz(v.x(), v.y(), v.z());
  }

  size_t dimension() const { return 3; } //!< @brief return 2

  Scalar const& x() const { return x_; } //!< @brief access x
  Scalar const& y() const { return y_; } //!< @brief access y
  Scalar const& z() const { return z_; } //!< @brief access z

  Scalar& xGet() { return x_; } //!< @brief access x with non-constant reference
  Scalar& yGet() { return y_; } //!< @brief access y with non-constant reference
  Scalar& zGet() { return z_; } //!< @brief access z with non-constant reference

  Scalar const& x(Scalar const& s) { x_ = s; return x(); } //!< @brief modify x
  Scalar const& y(Scalar const& s) { y_ = s; return y(); } //!< @brief modify y
  Scalar const& z(Scalar const& s) { z_ = s; return z(); } //!< @brief modify z

  //! @brief modify x, y and z
  Vector3D& xyz(Scalar const& sx, Scalar const& sy, Scalar const& sz) {
    x(sx);
    y(sy);
    z(sz);
    return *this;
  }

  //! @brief access the \c i -th scalar (0 => x, 1 => y, 2 => z) {
  Scalar const& scalar(size_t i) const {
    static Scalar tmp(0);
    return (i == 0 ? x() : (i == 1 ? y() : (i == 2 ? z() : tmp)));
  }

  //! @brief modivy the \c i -th scalar (0 => x, 1 => y, 2 => z)
  Scalar const& scalar(size_t i, Scalar const& s) {
    return (i == 0 ? x(s) : (i == 1 ? y(s) : (i == 2 ? z(s) : s )));
  }

  //! @brief return \a +(*this)
  Vector3D positive() const {
    return *this;
  }

  //! @brief return \a -(*this)
  Vector3D negative() const {
    return Vector3D(-x(), -y(), -z());
  }

  //! @brief return \a (*this)+v
  Vector3D add(Vector3D const& v) const {
    return Vector3D(x() + v.x(), y() + v.y(), z() + v.z());
  }

  //! @brief Let itself add v
  Vector3D& added(Vector3D const& v) {
    return xyz(x() + v.x(), y() + v.y(), z() + v.z());
  }

  //! @brief return \a (*this)-v
  Vector3D sub(Vector3D const& v) const {
    return Vector3D(x() - v.x(), y() - v.y(), z() - v.z());
  }

  //! @brief Let itself substract v
  Vector3D& subed(Vector3D const& v) {
    return xyz(x() - v.x(), y() - v.y(), z() - v.z());
  }

  //! @brief return \a (*this)*s , where s is a scalar
  Vector3D mul(Scalar const& s) const {
    return Vector3D(x() * s, y() * s, z() * s);
  }

  //! @brief Let itself mulitple s
  Vector3D& muled(Scalar const& s) {
    return xyz(x() * s, y() * s, z() * s);
  }

  //! @brief return \a (*this)/s , where s is a scalar
  Vector3D div(Scalar const& s) const {
    return Vector3D(x() / s, y() / s, z() / s);
  }

  //! @brief Let itself divide s
  Vector3D& dived(Scalar const& s) {
    return xyz(x() / s, y() / s, z() / s);
  }

  //! @brief same as dot(v)
  Scalar mul(Vector3D const& v) const {
    return dot(v);
  }

  //! @brief dot
  Scalar dot(Vector3D const& v) const {
    return x() * v.x() + y() * v.y() + z() * v.z();
  }

  //! @brief cross
  Vector3D cross(Vector3D const& v) const {
    return Vector3D(y() * v.z() - z() * v.y(),
                    z() * v.x() - x() * v.z(),
                    x() * v.y() - y() * v.x());
  }

  //! @brief crossed
  Vector3D& crossed(Vector3D const& v) {
    return copyFrom(cross(v));
  }

  //! @brief sqrt of length2
  Scalar length() const {
    return Scalar(sqrt(double(length2())));
  }

  //! @brief same as \a dot(*this)
  Scalar length2() const {
    return dot(*this);
  }

  //! @brief return normalize form of itself
  Vector3D normalize() const {
    return div(length());
  }

  //! @brief normalize itself
  Vector3D& normalized() {
    return dived(length());
  }

  //! @brief return rotate \a theta degree by \a axis of itself
  Vector3D rotate(Vector3D const& axis, double theta) const {
    Vector3D a(axis.normalize());
    Vector3D xx(sub(a)        .mul(cos(theta)));
    Vector3D yy(a.cross(*this).mul(sin(theta)));
    return a.mul(a.dot(*this)).add(xx).add(yy);
  }

  //! @brief Let itself rotate \a theta degree
  Vector3D& rotated(Vector3D const& axis, double theta) {
    return copyFrom(rotate(axis, theta));
  }

  //! @brief return reflect from given vector \a v
  Vector3D reflect(Vector3D const& v) const {
    return v.mul(v.dot(*this) * 2 / v.length2()).sub(*this);
  }

  //! @brief reflect itself given vector \a v
  Vector3D& reflected(Vector3D const& v) {
    return copyFrom(reflecte(v));
  }

  //! @brief return a 3x1 matrix form of itself
  Matrix<Scalar> matrix() const {
    static Matrix<Scalar> ret(3, 1, Scalar(0));
    ret(0, 0, x());
    ret(1, 0, y());
    ret(2, 0, z());
    return ret;
  }

  //! @brief return a 4x1 matrix form of itself
  Matrix<Scalar> matrix(Scalar const& homo) const {
    static Matrix<Scalar> ret(4, 1, Scalar(0));
    ret(0, 0, x());
    ret(1, 0, y());
    ret(2, 0, z());
    ret(3, 0, homo);
    return ret;
  }

  Scalar const& operator()(size_t n) const {
    return (n == 0 ? x() : (n == 1 ? y() : z()));
  }

  Vector3D& operator()(Scalar const& sx, Scalar const& sy, Scalar const& sz) {
    return xyz(sx, sy, sz);
  }

  Vector3D operator+() const { return positive(); }
  Vector3D operator-() const { return negative(); }

  Vector3D operator+(Vector3D const& v) const { return add(v); }
  Vector3D operator-(Vector3D const& v) const { return sub(v); }
  Vector3D operator*(Scalar   const& s) const { return mul(s); }
  Vector3D operator/(Scalar   const& s) const { return div(s); }
  Scalar   operator*(Vector3D const& v) const { return mul(v); }

  Vector3D& operator =(Vector3D const& v) { return copyFrom(v); }
  Vector3D& operator+=(Vector3D const& v) { return added(v); }
  Vector3D& operator-=(Vector3D const& v) { return subed(v); }
  Vector3D& operator*=(Scalar   const& s) { return muled(s); }
  Vector3D& operator/=(Scalar   const& s) { return dived(s); }
};

} // meow

#endif // geo_Vectors_H__