#ifndef __JPOSITION3D__
#define __JPOSITION3D__

#include <istream>
#include <ostream>
#include <cmath>

#include "JIO/JSerialisable.hh"
#include "JLang/JManip.hh"
#include "JGeometry2D/JVector2D.hh"
#include "JGeometry3D/JVector3D.hh"
#include "JGeometry3D/JAngle3D.hh"
#include "JGeometry3D/JVersor3D.hh"
#include "JGeometry3D/JVersor3Z.hh"
#include "JGeometry3D/JRotation3D.hh"
#include "JGeometry3D/JQuaternion3D.hh"


/**
 * \author mdejong
 */

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

namespace JGEOMETRY3D {

  using JIO::JReader;
  using JIO::JWriter;
  using JGEOMETRY2D::JVector2D;
    

  /**
   * Data structure for position in three dimensions.
   */
  class JPosition3D :
    public JVector3D
  {
  public:

    using JVector3D::getDot;
    using JVector3D::transform;
    

    /**
     * Default constructor.
     */
    JPosition3D() :
      JVector3D()
    {}

    
    /**
     * Constructor.
     *
     * \param  pos           position
     */
    JPosition3D(const JVector3D& pos) :
      JVector3D(pos)
    {}


    /**
     * Constructor.
     *
     * \param  angle         angle
     */
    JPosition3D(const JAngle3D& angle) :
      JVector3D(angle.getDX(),
		angle.getDY(),
		angle.getDZ())
    {}


    /**
     * Constructor.
     *
     * \param  dir           direction
     */
    JPosition3D(const JVersor3D& dir) :
      JVector3D(dir.getDX(),
		dir.getDY(),
		dir.getDZ())
    {}


    /**
     * Constructor.
     *
     * \param  dir           direction
     */
    JPosition3D(const JVersor3Z& dir) :
      JVector3D(dir.getDX(),
		dir.getDY(),
		dir.getDZ())
    {}

    
    /**
     * Constructor.
     *
     * \param  pos           2D position
     * \param  z             z value
     */
    JPosition3D(const JVector2D& pos,
		const double     z) :
      JVector3D(pos.getX(), pos.getY(), z)
    {}


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

    /**
     * Get position.
     *
     * \return               position
     */
    const JPosition3D& getPosition() const
    {
      return static_cast<const JPosition3D&>(*this);
    }
    

    /**
     * Get position.
     *
     * \return               position
     */
    JPosition3D& getPosition()
    {
      return static_cast<JPosition3D&>(*this);
    }
    

    /**
     * Set position.
     *
     * \param  pos           position
     */
    void setPosition(const JVector3D& pos)
    {
      static_cast<JVector3D&>(*this) = pos;
    }

    
    /**
     * Type conversion operator.
     *
     * \return               angle
     */
    operator JAngle3D() const
    {
      return JAngle3D(getX(), getY(), getZ());
    }


    /**
     * Type conversion operator.
     *
     * \return               direction
     */
    operator JVersor3D() const
    {
      return JVersor3D(getX(), getY(), getZ());
    }

    
    /**
     * Rotate.
     *
     * \param  R             rotation matrix
     * \return               this position
     */
    JPosition3D& rotate(const JRotation3D& R)
    {
      R.rotate(__x, __y, __z);

      return *this;
    }

    
    /**
     * Rotate back.
     *
     * \param  R             rotation matrix
     * \return               this position
     */
    JPosition3D& rotate_back(const JRotation3D& R)
    {
      R.rotate_back(__x, __y, __z);

      return *this;
    }

    
    /**
     * Rotate around X-axis.
     *
     * \param  R             rotation matrix
     * \return               this position
     */
    JPosition3D& rotate(const JRotation3X& R)
    {
      R.rotate(__y, __z);

      return *this;
    }

    
    /**
     * Rotate back around X-axis.
     *
     * \param  R             rotation matrix
     * \return               this position
     */
    JPosition3D& rotate_back(const JRotation3X& R)
    {
      R.rotate_back(__y, __z);

      return *this;
    }

    
    /**
     * Rotate around Y-axis.
     *
     * \param  R             rotation matrix
     * \return               this position
     */
    JPosition3D& rotate(const JRotation3Y& R)
    {
      R.rotate(__x, __z);

      return *this;
    }

    
    /**
     * Rotate back around Y-axis.
     *
     * \param  R             rotation matrix
     * \return               this position
     */
    JPosition3D& rotate_back(const JRotation3Y& R)
    {
      R.rotate_back(__x, __z);

      return *this;
    }

    
    /**
     * Rotate around Z-axis.
     *
     * \param  R             rotation matrix
     * \return               this position
     */
    JPosition3D& rotate(const JRotation3Z& R)
    {
      R.rotate(__x, __y);

      return *this;
    }

    
    /**
     * Rotate back around Z-axis.
     *
     * \param  R             rotation matrix
     * \return               this position
     */
    JPosition3D& rotate_back(const JRotation3Z& R)
    {
      R.rotate_back(__x, __y);

      return *this;
    }

    
    /**
     * Rotate.
     *
     * \param  Q             quaternion
     * \return               this position
     */
    JPosition3D& rotate(const JQuaternion3D& Q)
    {
      Q.rotate(__x, __y, __z);
      
      return *this;
    }


    /**
     * Rotate back.
     *
     * \param  Q             quaternion
     * \return               this position
     */
    JPosition3D& rotate_back(const JQuaternion3D& Q)
    {
      Q.rotate_back(__x, __y, __z);

      return *this;
    }


    /**
     * Transform position.
     *
     * The final position is obtained as follows:
     *   -# rotation of the position according matrix R;
     *   -# offset position with pos;
     *   -# rotation of position 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)
    {
      // rotate geometry to system with particle direction along z-axis

      rotate(R);

      // offset with respect to origin

      sub(pos);

      // rotate geometry to x-z plane

      __x = sqrt(__x*__x + __y*__y);
      __y = 0.0;
    }


    /**
     * Transform back position.
     *
     * The final position is obtained as follows:
     *   -# offset position with position pos;
     *   -# rotation of postion according matrix R;
     *
     * \param  R             rotation matrix
     * \param  pos           position of origin (before rotation)
     */
    void transform_back(const JRotation3D& R,
                        const JVector3D&   pos)
    {
      // offset with respect to origin

      add(pos);

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

      rotate_back(R);
    }


    /**
     * Get dot product.
     *
     * \param  angle         angle
     * \return               dot product
     */
    double getDot(const JAngle3D& angle) const
    {
      return
	getX() * angle.getDX() +
	getY() * angle.getDY() +
	getZ() * angle.getDZ();
    }

    
    /**
     * Get dot product.
     *
     * \param  dir           direction
     * \return               dot product
     */
    double getDot(const JVersor3D& dir) const
    {
      return
	getX() * dir.getDX() +
	getY() * dir.getDY() +
	getZ() * dir.getDZ();
    }

    
    /**
     * Get dot product.
     *
     * \param  dir           direction
     * \return               dot product
     */
    double getDot(const JVersor3Z& dir) const
    {
      return
	getX() * dir.getDX() +
	getY() * dir.getDY() +
	getZ() * dir.getDZ();
    }
    

    /**
     * Read position from input.
     *
     * \param  in            input stream
     * \param  position      position
     * \return               input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JPosition3D& position)
    {
      in >> position.__x  >> position.__y  >> position.__z;

      return in;
    }


    /**
     * Write position to output.
     *
     * \param  out           output stream
     * \param  position      position
     * \return               output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JPosition3D& position)
    {
      const JFormat format(out, getFormat<JPosition3D>(JFormat_t(9, 3, std::ios::fixed | std::ios::showpos)));

      out << format << position.getX() << ' '
	  << format << position.getY() << ' '
	  << format << position.getZ();

      return out;
    }


    /**
     * Read position from input.
     *
     * \param  in            reader
     * \param  position      position
     * \return               reader
     */
    friend inline JReader& operator>>(JReader& in, JPosition3D& position)
    {
      return in >> position.__x >> position.__y >> position.__z;
    }


    /**
     * Write position to output.
     *
     * \param  out           writer
     * \param  position      position
     * \return               writer
     */
    friend inline JWriter& operator<<(JWriter& out, const JPosition3D& position)
    {
      return out << position.__x << position.__y << position.__z;
    }
  };
}

#endif