#ifndef __JAXIS3D__
#define __JAXIS3D__

#include <istream>
#include <ostream>

#include "JGeometry3D/JPoint3D.hh"
#include "JGeometry3D/JVersor3D.hh"
#include "JGeometry3D/JRotation3D.hh"
#include "JIO/JSerialisable.hh"


namespace JGEOMETRY3D {

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


  /**
   * Axis object.
   *
   * An axis object consists of a position and a direction.
   */
  class JAxis3D : 
    public JPoint3D,
    public JVersor3D
  {
  public:
    /**
     * Default constructor.
     */
    JAxis3D() :
      JPoint3D (),
      JVersor3D()
    {}


    /**
     * Constructor.
     *
     * \param  pos          position
     */
    JAxis3D(const JVector3D& pos) :
      JPoint3D (pos), 
      JVersor3D()
    {}


    /**
     * Constructor.
     *
     * \param  pos          position
     * \param  dir          direction
     */
    JAxis3D(const JVector3D& pos,
	    const JVersor3D& dir) :
      JPoint3D (pos), 
      JVersor3D(dir)
    {}

    
    /**
     * Get axis.
     *
     * \return               axis
     */
    const JAxis3D& getAxis() const 
    { 
      return *this; 
    }

    
    /**
     * Set axis.
     *
     * \param  axis          axis
     */
    void setAxis(const JAxis3D& axis)
    { 
      *this = axis; 
    }
    

    /**
     * Rotate axis.
     *
     * \param  R             rotation matrix
     * \return               this axis
     */
    JAxis3D& rotate(const JRotation3D& R)
    {
      static_cast<JPoint3D& >(*this).rotate(R);
      static_cast<JVersor3D&>(*this).rotate(R);

      return *this;
    }


    /**
     * Rotate back axis.
     *
     * \param  R             rotation matrix
     * \return               this axis
     */
    JAxis3D& rotate_back(const JRotation3D& R)
    {
      static_cast<JPoint3D& >(*this).rotate_back(R);
      static_cast<JVersor3D&>(*this).rotate_back(R);

      return *this;
    }


    /**
     * Transform axis.
     *
     * The final position and direction are obtained as follows:
     *   1/ rotation of the position and the direction according matrix R;
     *   2/ offset position with pos;
     *   3/ rotation of position and direction around z-axis, such that final position lies in x-z plane;
     *
     * \param  R          rotation matrix
     * \param  pos        position of origin (after rotation)
     */
    void transform(const JRotation3D& R,
                   const JVector3D&   pos)
    {
      JPoint3D&  u = getPosition ();
      JVersor3D& v = getDirection();

      // rotate axis to system with particle direction along z-axis

      u.rotate(R);
      v.rotate(R);

      // offset with respect to origin

      u.sub(pos);

      // rotate axis to x-z plane

      const double phi = atan2(u.getY(), u.getX());

      const JRotation3Z r(-phi);

      u.rotate(r);
      v.rotate(r);
    }


    /**
     * Transform back axis.
     *
     * The final position and direction are obtained as follows:
     *   1/ offset position with position pos;
     *   2/ rotation of position and direction according matrix R;
     *
     * \param  R          rotation matrix
     * \param  pos        position of origin (before rotation)
     */
    void transform_back(const JRotation3D& R,
                        const JVector3D&   pos)
    {
      JPoint3D&  u = getPosition ();
      JVersor3D& v = getDirection();

      // offset with respect to origin

      u.add(pos);

      // rotate back axis to system with original particle direction

      u.rotate_back(R);
      v.rotate_back(R);
    }


    /**
     * Read JAxis from input.
     *
     * \param  in            input stream
     * \param  axis          axis
     * \return               input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JAxis3D& axis)
    {
      in >> axis.__x  >> axis.__y  >> axis.__z;
      in >> axis.__dx >> axis.__dy >> axis.__dz;

      return in;
    }


    /**
     * Write JAxis to output.
     *
     * \param  out           output stream
     * \param  axis          axis
     * \return               output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JAxis3D& axis)
    {
      out << axis.getX()  << ' ' << axis.getY()  << ' ' << axis.getZ();
      out << ' ';
      out << axis.getDX() << ' ' << axis.getDY() << ' ' << axis.getDZ();

      return out;
    }


    /**
     * Read JAxis from input.
     *
     * \param  in            JReader
     * \param  axis          axis
     * \return               JReader
     */
    friend inline JReader& operator>>(JReader& in, JAxis3D& axis)
    {
      in >> static_cast<JPoint3D&> (axis);
      in >> static_cast<JVersor3D&>(axis);

      return in;
    }


    /**
     * Write JAxis to output.
     *
     * \param  out           JWriter
     * \param  axis          axis
     * \return               JWriter
     */
    friend inline JWriter& operator<<(JWriter& out, const JAxis3D& axis)
    {
      out << static_cast<const JPoint3D&> (axis);
      out << static_cast<const JVersor3D&>(axis);

      return out;
    }
  };
}

#endif