#ifndef __JPOINT3D__
#define __JPOINT3D__

#include <cmath>

#include "JGeometry2D/JVector2D.hh"
#include "JGeometry3D/JVector3D.hh"
#include "JGeometry3D/JAngle3D.hh"
#include "JGeometry3D/JRotation3D.hh"


namespace JGEOMETRY3D {

  namespace {
    using JGEOMETRY2D::JVector2D;
  }
    

  /**
   * Data structure for point in three dimensions.
   */
  class JPoint3D :
    public JVector3D
  {
  public:
    /**
     * Default constructor.
     */
    JPoint3D() :
      JVector3D()
    {}

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

    
    /**
     * Constructor.
     *
     * \param  pos      2D position
     * \param  z        z value
     */
    JPoint3D(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
     */
    JPoint3D(const double x,
	     const double y,
	     const double z) :
      JVector3D(x,y,z)
    {}


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

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

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

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


    /**
     * Transform.
     *
     * \param  T        matrix
     */
    JPoint3D& transform(const JMatrix3D& T)
    {
      T.transform(__x, __y, __z);

      return *this;
    }

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

      return *this;
    }

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

      return *this;
    }

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

      return *this;
    }

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

      return *this;
    }

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

      return *this;
    }

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

      return *this;
    }

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

      return *this;
    }

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

      return *this;
    }


    /**
     * Transform position.
     *
     * The final position is obtained as follows:
     *   1/ rotation of the position according matrix R;
     *   2/ offset position with pos;
     *   3/ 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:
     *   1/ offset position with position pos;
     *   2/ 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 squared of distance to point.
     *
     * \param  point    JPoint3D
     * \return          distance
     */
    double getDistanceSquared(const JPoint3D& point) const
    {
      const double dx = getX() - point.getX();
      const double dy = getY() - point.getY();
      const double dz = getZ() - point.getZ();
      
      return dx*dx + dy*dy + dz*dz;
    }


    /**
     * Get distance to point.
     *
     * \param  point    JPoint3D
     * \return          distance
     */
    double getDistance(const JPoint3D& point) const
    {
      return sqrt(getDistanceSquared(point));
    }


    /**
     * Get dot product.
     *
     * \param  point     JPoint3D
     * \return           dot product
     */
    double getDot(const JPoint3D& point) const
    {
      return
	getX() * point.getX() + 
	getY() * point.getY() +
	getZ() * point.getZ();
    }
  };


  /**
   * Distance.
   *
   * \param  first     JPoint3D
   * \param  second    JPoint3D
   * \return           dot product
   */
  inline double distance(const JPoint3D& first,
			 const JPoint3D& second)
  {
    return first.getDistance(second);
  }


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

#endif