#ifndef __JVECTOR3D__
#define __JVECTOR3D__

#include <limits>
#include <utility>
#include <cmath>

#include "JMath/JMath.hh"
#include "JGeometry2D/JVector2D.hh"
#include "JGeometry3D/JMatrix3D.hh"


/**
 * \author mdejong
 */

namespace JGEOMETRY3D {}
namespace JPP { using namespace JGEOMETRY3D; }

namespace JGEOMETRY3D {

  using JMATH::JMath;
  using JGEOMETRY2D::JVector2D;
  using JGEOMETRY2D::JRangeX;                 //!< Type definition of range along x-axis.
  using JGEOMETRY2D::JRangeY;                 //!< Type definition of range along y-axis.
  typedef std::pair<double, double>  JRangeZ; //!< Type definition of range along z-axis.
  

  /**
   * Data structure for vector in three dimensions.
   *
   * This class implements the JMATH::JMath interface.
   */
  class JVector3D :
    public JMath<JVector3D>
  {
  public:
    /**
     * Default constructor.
     */
    JVector3D() :
      __x(0.0),
      __y(0.0),
      __z(0.0)
    {}


    /**
     * Constructor.
     *
     * \param  x             x value
     * \param  y             y value
     * \param  z             z value
     */
    JVector3D(const double x,
	      const double y,
	      const double z) :
      __x(x),
      __y(y),
      __z(z)
    {}


    /**
     * Constructor.
     *
     * \param  vector        (x,y) values
     * \param  z              z    value
     */
    JVector3D(const JVector2D& vector,
	      const double     z) :
      __x(vector.getX()),
      __y(vector.getY()),
      __z(z)
    {}

    
    /**
     * Type conversion operator.
     *
     * \return               JVector2D
     */
    operator JVector2D() const
    {
      return JVector2D(this->getX(), this->getY());
    }


    /**
     * Get x position.
     *
     * \return               x position
     */
    double getX() const 
    { 
      return __x;
    }

    /**
     * Get y position.
     *
     * \return               y position
     */
    double getY() const 
    { 
      return __y;
    }


    /**
     * Get z position.
     *
     * \return               z position
     */
    double getZ() const 
    { 
      return __z;
    }


    /** 
     * Negate vector.
     *
     * \return               this vector
     */
    JVector3D& negate()
    {
      __x = -__x;
      __y = -__y;
      __z = -__z;

      return *this;
    }


    /** 
     * Add vector.
     *
     * \param  vector        vector
     * \return               this vector
     */
    JVector3D& add(const JVector3D& vector)
    {
      __x += vector.getX();
      __y += vector.getY();
      __z += vector.getZ();

      return *this;
    }


    /** 
     * Subtract vector.
     *
     * \param  vector        vector
     * \return               this vector
     */
    JVector3D& sub(const JVector3D& vector)
    {
      __x -= vector.getX();
      __y -= vector.getY();
      __z -= vector.getZ();

      return *this;
    }


    /** 
     * Scale vector.
     *
     * \param  factor        multiplication factor
     * \return               this vector
     */
    JVector3D& mul(const double factor)
    {
      __x *= factor;
      __y *= factor;
      __z *= factor;

      return *this;
    }


    /** 
     * Scale vector.
     *
     * \param  factor        division factor
     * \return               this vector
     */
    JVector3D& div(const double factor)
    {
      __x /= factor;
      __y /= factor;
      __z /= factor;

      return *this;
    }


    /**
     * Transform.
     *
     * \param  T             matrix
     * \return               this vector
     */
    JVector3D& transform(const JMatrix3D& T)
    {
      T.transform(__x, __y, __z);

      return *this;
    }
    

    /** 
     * Check equality.
     *
     * \param  vector        vector
     * \param  precision     precision
     * \return               true if vectors are equal; else false
     */
    bool equals(const JVector3D& vector,
                const double     precision = std::numeric_limits<double>::min()) const
    {
      return (fabs(getX() - vector.getX()) <= precision &&
              fabs(getY() - vector.getY()) <= precision &&
              fabs(getZ() - vector.getZ()) <= precision);
    }


    /**
     * Get length squared.
     *
     * \return               square of length
     */
    double getLengthSquared() const 
    { 
      return getX()*getX() + getY()*getY() + getZ()*getZ();
    }


    /**
     * Get length.
     *
     * \return               length
     */
    double getLength() const 
    { 
      return sqrt(getLengthSquared());
    }


    /**
     * Get squared of distance to point.
     *
     * \param  pos           position
     * \return               square of distance
     */
    double getDistanceSquared(const JVector3D& pos) const
    {
      return JVector3D(pos).sub(*this).getLengthSquared();
    }


    /**
     * Get distance to point.
     *
     * \param  pos           position
     * \return               distance
     */
    double getDistance(const JVector3D& pos) const
    {
      return sqrt(getDistanceSquared(pos));
    }


    /**
     * Get dot product.
     *
     * \param  vector        vector
     * \return               dot product
     */
    double getDot(const JVector3D& vector) const
    {
      return
	getX() * vector.getX() +
	getY() * vector.getY() +
	getZ() * vector.getZ();
    }


    /**
     * Get cross product.
     * Note that this vector should not overlap with the first or second vector,
     *
     * \param  first         first  vector
     * \param  second        second vector
     * \return               this   vector
     */
    JVector3D& getCross(const JVector3D& first,
			const JVector3D& second)
    {
      __x = first .getY() * second.getZ()  -  second.getY() * first .getZ();
      __y = second.getX() * first .getZ()  -  first .getX() * second.getZ();
      __z = first .getX() * second.getY()  -  second.getX() * first .getY();

      return *this;
    }

  protected:
    double __x;
    double __y;
    double __z;
  };


  static const JVector3D JVector3X_t(1,0,0);   //!< unit x-vector
  static const JVector3D JVector3Y_t(0,1,0);   //!< unit y-vector
  static const JVector3D JVector3Z_t(0,0,1);   //!< unit z-vector
}

#endif