#ifndef __JVERSOR3D__
#define __JVERSOR3D__

#include <limits>
#include <cmath>

#include "JMath/JMath.hh"
#include "JMath/JConstants.hh"


/**
 * \author mdejong
 */

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

namespace JGEOMETRY3D {

  using JMATH::JMath;


  /**
   * Data structure for normalised vector in three dimensions.
   */
  class JVersor3D :
    public JMath<JVersor3D>
  {
  public:
    /**
     * Default constructor.
     * This constructor yields a unit z-vector.  
     */
    JVersor3D() :
      __dx(0.0),
      __dy(0.0),
      __dz(1.0)
    {}


    /**
     * Constructor.
     *
     * \param  dx            dx value
     * \param  dy            dy value
     * \param  dz            dz value
     */
    JVersor3D(const double dx,
              const double dy,
              const double dz) :
      __dx(dx),
      __dy(dy),
      __dz(dz)
    {
      normalise();
    }

    
    /**
     * Negate versor.
     *
     * \return               this versor
     */
    JVersor3D& negate()
    {
      __dx = -__dx;
      __dy = -__dy;
      __dz = -__dz;

      return *this;
    }


    /**
     * Check equality.
     *
     * \param  versor        versor
     * \param  precision     precision
     * \return               true if versors are equal; else false
     */
    bool equals(const JVersor3D& versor,
                const double     precision = std::numeric_limits<double>::min()) const
    {
      return (fabs(getDX() - versor.getDX()) <= precision &&
              fabs(getDY() - versor.getDY()) <= precision &&
              fabs(getDZ() - versor.getDZ()) <= precision);
    }


    /**
     * Get x direction.
     *
     * \return               x direction
     */
    double getDX() const
    {
      return __dx;
    }


    /**
     * Get y direction.
     *
     * \return               y direction
     */
    double getDY() const
    {
      return __dy;
    }


    /**
     * Get z direction.
     *
     * \return               z direction
     */
    double getDZ() const
    {
      return __dz;
    }


    /**
     * Get theta angle.
     *
     * \return               theta angle [rad]
     */
    double getTheta() const 
    {
      if      (__dz > +1.0)
	return 0.0;
      else if (__dz < -1.0)
	return JMATH::PI;
      else
	return acos(__dz);
    }


    /**
     * Get phi angle.
     *
     * \return               phi angle [rad]
     */
    double getPhi() const 
    { 
      return atan2(__dy, __dx);
    }


    /**
     * Get dot product.
     *
     * \param  versor        versor
     * \return               dot product
     */
    double getDot(const JVersor3D& versor) const
    {
      return
	getDX() * versor.getDX() + 
	getDY() * versor.getDY() +
	getDZ() * versor.getDZ();
    }


    /**
     * Get cross product.
     * Note that this versor should not overlap with the first or second versor,
     *
     * \param  first         first  versor
     * \param  second        second versor
     * \return               this   versor
     */
    JVersor3D& getCross(const JVersor3D& first,
                        const JVersor3D& second)
    {
      __dx = first .getDY() * second.getDZ()  -  second.getDY() * first .getDZ();
      __dy = second.getDX() * first .getDZ()  -  first .getDX() * second.getDZ();
      __dz = first .getDX() * second.getDY()  -  second.getDX() * first .getDY();

      normalise();

      return *this;
    }


    /**
     * Normalise versor.
     * This operation may set the result to the unit z-vector.
     *
     * \return               this versor
     */
    JVersor3D& normalise()
    {
      const double v = sqrt(getDX()*getDX() + getDY()*getDY() + getDZ()*getDZ());

      if (v != 0.0) {
	__dx /= v;
	__dy /= v;
	__dz /= v;
      }

      return *this;
    }

  protected:
    double __dx;
    double __dy;
    double __dz;
  };


  static const JVersor3D JVersor3X_t(1,0,0);   //!< unit x-vector
  static const JVersor3D JVersor3Y_t(0,1,0);   //!< unit y-vector
  static const JVersor3D JVersor3Z_t(0,0,1);   //!< unit z-vector
}

#endif