#ifndef __JDIRECTION3D__
#define __JDIRECTION3D__

#include <cmath>

#include "JGeometry3D/JAngle3D.hh"
#include "JGeometry3D/JPoint3D.hh"
#include "JGeometry3D/JRotation3D.hh"


namespace JGEOMETRY3D {


  /**
   * Data structure for direction in three dimensions.
   */
  class JDirection3D :
    public JAngle3D 
  {
  public:
    /**
     * Default constructor.
     */
    JDirection3D() :
      JAngle3D()
    {}


    /**
     * Constructor.
     *
     * \param  theta    theta angle [rad]
     * \param  phi      phi angle [rad]
     */
    JDirection3D(const double theta,
		 const double phi) :
      JAngle3D(theta, phi)
    {}


    /**
     * Constructor.
     *
     * \param  vector   vector
     */
    JDirection3D(JVector3D vector) :
      JAngle3D()
    {
      __theta = acos (vector.getZ() / vector.getLength());
      __phi   = atan2(vector.getY(),vector.getX());
    }


    /**
     * Constructor.
     *
     * \param  x        x value
     * \param  y        y value
     * \param  z        z value
     */
    JDirection3D(const double x,
		 const double y,
		 const double z) :
      JAngle3D()
    {
      __theta = acos (z / sqrt(x*x + y*y + z*z));
      __phi   = atan2(y,x);
    }


    /**
     * Transform.
     *
     * \param  R        matrix
     */
    JDirection3D& transform(const JMatrix3D& R)
    {
      *this = JPoint3D(*this).transform(R);

      return *this;
    }


    /**
     * Rotate.
     *
     * \param  R        rotation matrix
     */
    JDirection3D& rotate(const JRotation3D& R)
    {
      *this = JPoint3D(*this).rotate(R);

      return *this;
    }


    /**
     * Rotate back.
     *
     * \param  R        rotation matrix
     */
    JDirection3D& rotate_back(const JRotation3D& R)
    {
      *this = JPoint3D(*this).rotate_back(R);

      return *this;
    }


    /**
     * Rotate around X-axis.
     *
     * \param  R        rotation matrix
     */
    JDirection3D& rotate(const JRotation3X& R)
    {
      *this = JPoint3D(*this).rotate(R);     

      return *this;
    }


    /**
     * Rotate back around X-axis.
     *
     * \param  R        rotation matrix
     */
    JDirection3D& rotate_back(const JRotation3X& R)
    {
      *this = JPoint3D(*this).rotate_back(R);

      return *this;
    }


    /**
     * Rotate around Y-axis.
     *
     * \param  R        rotation matrix
     */
    JDirection3D& rotate(const JRotation3Y& R)
    {
      *this = JPoint3D(*this).rotate(R);     

      return *this;
    }


    /**
     * Rotate back around Y-axis.
     *
     * \param  R        rotation matrix
     */
    JDirection3D& rotate_back(const JRotation3Y& R)
    {
      *this = JPoint3D(*this).rotate_back(R);

      return *this;
    }


    /**
     * Rotate around Z-axis.
     *
     * \param  R        rotation matrix
     */
    JDirection3D& rotate(const JRotation3Z& R)
    {
      *this = JPoint3D(*this).rotate(R);     

      return *this;
    }


    /**
     * Rotate back around Z-axis.
     *
     * \param  R        rotation matrix
     */
    JDirection3D& rotate_back(const JRotation3Z& R)
    {
      *this = JPoint3D(*this).rotate_back(R);

      return *this;
    }


    /**
     * Get dot product.
     *
     * \param  dir       JDirection3D
     * \return           dot product
     */
    double getDot(const JDirection3D& dir) const
    {
      return
	cos(getTheta()) * cos(dir.getTheta()) +
	sin(getTheta()) * sin(dir.getTheta()) *
	cos(getPhi() - dir.getPhi());
    }
  };


  /**
   * Dot product.
   *
   * \param  first     JDirection3D
   * \param  second    JDirection3D
   * \return           dot product
   */
  inline double dot(const JDirection3D& first,
                    const JDirection3D& second)
  {
    return first.getDot(second);
  }
}

#endif