#ifndef __JANGLE3D__
#define __JANGLE3D__

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

#include "JIO/JSerialisable.hh"
#include "JLang/JManip.hh"
#include "JMath/JMath.hh"
#include "JMath/JConstants.hh"


/**
 * \author mdejong
 */

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

namespace JGEOMETRY3D {

  using JIO::JReader;
  using JIO::JWriter;
  using JMATH::JMath;


  /**
   * Data structure for angles in three dimensions.
   * This class serves as input to the rotation matrix JRotation3D.
   */
  class JAngle3D :
    public JMath<JAngle3D>
  {
  public:
    /**
     * Default constructor.
     */
    JAngle3D() :
      __theta(0.0),
      __phi  (0.0)
    {}


    /**
     * Constructor.
     *
     * \param  theta         theta angle [rad]
     * \param  phi           phi   angle [rad]
     */
    JAngle3D(const double theta,
	     const double phi) :
      __theta(theta),
      __phi  (phi)
    {}


    /**
     * Constructor.
     *
     * \param  x             x value
     * \param  y             y value
     * \param  z             z value
     */
    JAngle3D(const double x,
	     const double y,
	     const double z) :
      __theta(0.0),
      __phi  (0.0)
    {
      const double v = x*x + y*y + z*z;

      if (v != 0.0) {      
	__theta = acos (z / sqrt(v));
	__phi   = atan2(y, x);
      }
    }


    /**
     * Get theta angle.
     *
     * \return               theta angle
     */
    double getTheta() const 
    { 
      return __theta;
    }


    /**
     * Get phi angle.
     *
     * \return               phi angle
     */
    double getPhi() const 
    { 
      return __phi;
    }


    /**
     * Get x direction.
     *
     * \return               x direction
     */
    double getDX() const 
    { 
      return sin(__theta) * cos(__phi);
    }


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


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

    
    /**
     * Negate angle.
     *
     * \return               this angle
     */
    JAngle3D& negate()
    {
      __theta = -__theta;
      __phi   = -__phi;

      return *this;
    }


    /**
     * Add angle.
     *
     * \param  angle         angle
     * \return               this angle
     */
    JAngle3D& add(const JAngle3D& angle)
    {
      __theta += angle.getTheta();
      __phi   += angle.getPhi();

      return *this;
    }


    /**
     * Subtract angle.
     *
     * \param  angle         angle
     * \return               this angle
     */
    JAngle3D& sub(const JAngle3D& angle)
    {
      __theta -= angle.getTheta();
      __phi   -= angle.getPhi();

      return *this;
    }

    /**
     * Scale angle.
     *
     * \param  factor        multiplication factor
     * \return               this angle
     */
    JAngle3D& mul(const double factor)
    {
      __theta *= factor;
      __phi   *= factor;

      return *this;
    }


    /**
     * Scale angle.
     *
     * \param  factor        division factor
     * \return               this angle
     */
    JAngle3D& div(const double factor)
    {
      __theta /= factor;
      __phi   /= factor;

      return *this;
    }

    
    /**
     * Check equality.
     *
     * \param  angle         angle
     * \param  precision     precision
     * \return               true if angles are equal; else false
     */
    bool equals(const JAngle3D& angle,
                const double    precision = std::numeric_limits<double>::min()) const
    {
      return (fabs(getTheta() - angle.getTheta()) <= precision &&
              fabs(getPhi()   - angle.getPhi())   <= precision);
    }


    /**
     * Get dot product.
     *
     * \param  angle         angle
     * \return               dot product
     */
    double getDot(const JAngle3D& angle) const
    {
      return
	cos(getTheta()) * cos(angle.getTheta()) +
	sin(getTheta()) * sin(angle.getTheta()) *
	cos(getPhi() - angle.getPhi());
    }


    /**
     * Normalise angles.
     *
     * - theta angle will be between 0 and &pi;
     * - phi angle will be between -&pi; and +&pi;
     *
     * \return               this angle
     */
    JAngle3D& normalise()
    {
      using JMATH::PI;

      if (__theta > PI)  { do { __theta -= PI; } while (__theta > PI); }
      if (__theta < 0.0) { do { __theta += PI; } while (__theta < 0.0); }

      if (__phi > +PI) { do { __phi -= 2*PI; } while (__phi > +PI); }
      if (__phi < -PI) { do { __phi += 2*PI; } while (__phi < -PI); }

      return *this;
    }

    
    /**
     * Write angle from input.
     *
     * \param  in            input stream
     * \param  angle         angle
     * \return               input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JAngle3D& angle)
    {
      in >> angle.__theta >> angle.__phi;

      return in;
    }


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

      out << format << angle.getTheta() << ' '
	  << format << angle.getPhi();

      return out;
    }


    /**
     * Read angle from input.
     *
     * \param  in            reader
     * \param  angle         angle
     * \return               reader
     */
    friend inline JReader& operator>>(JReader& in, JAngle3D& angle)
    {
      in >> angle.__theta;
      in >> angle.__phi;

      return in;
    }


    /**
     * Write angle to output.
     *
     * \param  out           writer
     * \param  angle         angle
     * \return               writer
     */
    friend inline JWriter& operator<<(JWriter& out, const JAngle3D& angle)
    {
      out << angle.__theta;
      out << angle.__phi;

      return out;
    }

  protected:
    double __theta;
    double __phi;
  };
}

#endif