#ifndef __JVERSOR3D__
#define __JVERSOR3D__

#include <cmath>

#include "JIO/JSerialisable.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JGeometry3D/JRotation3D.hh"


namespace JGEOMETRY3D {

  namespace {
    using JIO::JReader;
    using JIO::JWriter;
  }


  /**
   * Data structure for normalised vector in three dimensions.
   */
  class JVersor3D
  {
  public:
    /**
     * Default constructor.
     */
    JVersor3D() :
      __dx(0.0),
      __dy(0.0),
      __dz(1.0)
    {}


    /**
     * Constructor.
     *
     * \param  angle    angle
     */
    JVersor3D(const JAngle3D& angle) :
      __dx(angle.getDX()),
      __dy(angle.getDY()),
      __dz(angle.getDZ())
    {}


    /**
     * 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();
    }


    /**
     * Get direction.
     *
     * \return               direction
     */
    const JVersor3D& getDirection() const 
    { 
      return static_cast<const JVersor3D&>(*this); 
    }

    
    /**
     * Get direction.
     *
     * \return               direction
     */
    JVersor3D& getDirection()
    { 
      return static_cast<JVersor3D&>(*this); 
    }

    
    /**
     * Set direction.
     *
     * \param  dir           direction
     */
    void setDirection(const JVersor3D& dir)
    { 
      static_cast<JVersor3D&>(*this) = dir; 
    }


    /**
     * Read JVersor3D from input.
     *
     * \param  in       JReader
     * \param  versor   JVersor3D
     * \return          JReader
     */
    friend inline JReader& operator>>(JReader& in, JVersor3D& versor)
    {
      in >> versor.__dx;
      in >> versor.__dy;
      in >> versor.__dz;

      return in;
    }


    /**
     * Write JVersor3D to output.
     *
     * \param  out      JWriter
     * \param  versor   JVersor3D
     * \return          JWriter
     */
    friend inline JWriter& operator<<(JWriter& out, const JVersor3D& versor)
    {
      out << versor.__dx;
      out << versor.__dy;
      out << versor.__dz;

      return out;
    }


    /**
     * Type conversion operator.
     *
     * \return          JAngle3D
     */
    operator JDirection3D() const
    {
      return JDirection3D(getDX(), getDY(), getDZ());
    }


    /**
     * 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 
    { 
      return acos(__dz);
    }


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

  
    /**
     * Transform.
     *
     * \param  T        matrix
     */
    JVersor3D& transform(const JMatrix3D& T)
    {
      T.transform(__dx, __dy, __dz);

      return *this;
    }

    
    /**
     * Rotate.
     *
     * \param  R        rotation matrix
     */
    JVersor3D& rotate(const JRotation3D& R)
    {
      R.rotate(__dx, __dy, __dz);

      return *this;
    }

    
    /**
     * Rotate back.
     *
     * \param  R        rotation matrix
     */
    JVersor3D& rotate_back(const JRotation3D& R)
    {
      R.rotate_back(__dx, __dy, __dz);

      return *this;
    }


    /**
     * Rotate around X-axis.
     *
     * \param  R        rotation matrix
     */
    JVersor3D& rotate(const JRotation3X& R)
    {
      R.rotate(__dy, __dz);

      return *this;
    }

    
    /**
     * Rotate back around X-axis.
     *
     * \param  R        rotation matrix
     */
    JVersor3D& rotate_back(const JRotation3X& R)
    {
      R.rotate_back(__dy, __dz);

      return *this;
    }

    
    /**
     * Rotate around Y-axis.
     *
     * \param  R        rotation matrix
     */
    JVersor3D& rotate(const JRotation3Y& R)
    {
      R.rotate(__dx, __dz);

      return *this;
    }

    
    /**
     * Rotate back around Y-axis.
     *
     * \param  R        rotation matrix
     */
    JVersor3D& rotate_back(const JRotation3Y& R)
    {
      R.rotate_back(__dx, __dz);

      return *this;
    }

    
    /**
     * Rotate around Z-axis.
     *
     * \param  R        rotation matrix
     */
    JVersor3D& rotate(const JRotation3Z& R)
    {
      R.rotate(__dx, __dy);

      return *this;
    }

    
    /**
     * Rotate back around Z-axis.
     *
     * \param  R        rotation matrix
     */
    JVersor3D& rotate_back(const JRotation3Z& R)
    {
      R.rotate_back(__dx, __dy);

      return *this;
    }


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


    /**
     * Normalise versor.
     */
    void normalise()
    {
      const double v = sqrt(getDX()*getDX() + getDY()*getDY() + getDZ()*getDZ());

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


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


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

#endif