#ifndef __JCYLINDER3D__
#define __JCYLINDER3D__

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

#include "JIO/JSerialisable.hh"
#include "JLang/JManip.hh"
#include "JMath/JConstants.hh"
#include "JGeometry2D/JCircle2D.hh"
#include "JGeometry3D/JVector3D.hh"
#include "JGeometry3D/JAxis3D.hh"


/**
 * \author mdejong
 */

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

namespace JGEOMETRY3D {

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

  /**
   * Cylinder object.
   *
   * The cylinder consists of a 2D circle in the (x,y) plane and a range in z
   * (i.e.\ axis of cylinder is parallel to the z-axis).
   */
  class JCylinder3D :
    public JCircle2D
  {
  public:
    /**
     * Type definition of intersection.
     */
    typedef std::pair<double, double> intersection_type;


    /**
     * Default constructor.
     */
    JCylinder3D() :
      JCircle2D(),
      zmin(0.0),
      zmax(0.0)
    {}
    

    /**
     * Constructor.
     *
     * \param  circle      2D circle in (x,y)
     * \param  zmin        minimal z
     * \param  zmax        maximal z
     */
    JCylinder3D(const JCircle2D& circle,
		const double     zmin,
		const double     zmax) :
      JCircle2D(circle)
    {
      this->zmin = zmin;
      this->zmax = zmax;
    }
    

    /**
     * Constructor.
     *
     * Determines smallest enclosing cylinder for any number of points.
     *
     * \param  __begin     begin of data
     * \param  __end       end   of data
     * \param  precision   precision
     */
    template<class T>
    JCylinder3D(T __begin,
		T __end,
		const double precision = std::numeric_limits<double>::epsilon()) :
      JCircle2D(__begin, __end, precision),
      zmin(0.0),
      zmax(0.0)
    {
      if (__begin != __end) {

	zmin = std::numeric_limits<double>::max();
	zmax = std::numeric_limits<double>::lowest();
	
	for (T i = __begin; i != __end; ++i) {
	  if (i->getZ() < zmin) zmin = i->getZ();
	  if (i->getZ() > zmax) zmax = i->getZ();
	}
      }
    }


    /**
     * Get minimal z position.
     *
     * \return             minimal z position
     */
    double getZmin() const 
    { 
      return zmin;
    }


    /**
     * Get maximal z position.
     *
     * \return             maximal z position
     */
    double getZmax() const 
    { 
      return zmax;
    }


    /**
     * Set minimal z position.
     *
     * \param  zmin        minimal z position
     */
    void setZmin(const double zmin)
    { 
      this->zmin = zmin;
    }


    /**
     * Set maximal z position.
     *
     * \param  zmax        maximal z position
     */
    void setZmax(const double zmax)
    { 
      this->zmax = zmax;
    }


    /**
     * Add position.
     *
     * \param  pos         position
     * \return             this cylinder
     */
    JCylinder3D& add(const JVector3D& pos)
    {
      static_cast<JPosition2D&>(*this).add(JPosition2D(pos.getX(), pos.getY()));

      zmin += pos.getZ();
      zmax += pos.getZ();
					   
      return *this;
    }

    
    /**
     * Add (safety) margin.
     *
     * \param  D           margin
     */
    void addMargin(const double D)
    {
      __r  += D;
      zmin -= D;
      zmax += D;
    }

    
    /**
     * Get volume.
     *
     * \return             volume
     */
    inline double getVolume() const
    {
      return (getZmax() - getZmin()) * JMATH::PI * getRadius() * getRadius();
    }


    /**
     * Get centre.
     *
     * \return            centre
     */
    JPosition3D getCenter() const
    {
      return JPosition3D(getPosition(), (getZmax() - getZmin())/2.0);
    }

    
    /**
     * Check whether given point is inside cylinder.
     *
     * \param  pos         position
     * \return             true if inside; else false
     */
    inline bool is_inside(const JVector3D& pos) const
    {
      return (pos.getZ() >= getZmin() && 
	      pos.getZ() <= getZmax() &&
	      JCircle2D::is_inside(JVector2D(pos.getX(), pos.getY())));
    }


    /**
     * Get distance between cylinder wall and given position.
     *
     * \param  pos         position
     * \return             distance
     */
    inline double getDistance(const JVector3D& pos) const
    {
      JVector2D D(pos);

      D.sub(*this);

      double R = D.getLength();
    
      if (R > this->getRadius()) {

	R -= this->getRadius();

	double dz = 0.0;

	if      (pos.getZ() > this->getZmax())
	  dz = pos.getZ() - this->getZmax();
	else if (pos.getZ() < this->getZmin())
	  dz = this->getZmin() - pos.getZ();
	else 
	  return R;

	return sqrt(R*R + dz*dz);

      } else {

	if      (pos.getZ() > this->getZmax())
	  return  pos.getZ() - this->getZmax();
	else if (pos.getZ() < this->getZmin())
	  return this->getZmin() - pos.getZ();
	else 
	  return 0.0;
      }
    }


    /**
     * Get square of distance between cylinder wall and given position.
     *
     * \param  pos         position
     * \return             square of distance
     */
    inline double getDistanceSquared(const JVector3D& pos) const
    {
      const double d = getDistance(pos);

      return d*d;
    }


    /**
     * Get intersection points of axis with cylinder.
     *
     * \param  axis        axis
     * \return             up and down stream positions along axis
     */
    inline intersection_type getIntersection(const JAxis3D& axis) const
    {
      double path[] = { 0.0, 0.0 };

      if (fabs(axis.getDZ()) != 0.0) {

	// intersection with bottom or top surface

	const double Z[] = { 
	  axis.getDZ() > 0 ? this->getZmin() : this->getZmax(),
	  axis.getDZ() > 0 ? this->getZmax() : this->getZmin()
	};

	for (int i = 0; i != 2; ++i) {

	  const double u = (Z[i] - axis.getZ()) / axis.getDZ();
	  const double x =  axis.getX() + u * axis.getDX() - this->getX();
	  const double y =  axis.getY() + u * axis.getDY() - this->getY();
      
	  if (x*x + y*y <= this->getRadius() * this->getRadius()) {
	    path[i] = u;
	  }
	}
      }
  
      if (fabs(axis.getDZ()) != 1.0) {
    
	// intersection with cylinder wall

	const double x  = axis.getX() - this->getX();
	const double y  = axis.getY() - this->getY();
	const double dx = axis.getDX();
	const double dy = axis.getDY();
	const double R  = this->getRadius();

	const double a  =   (dx * dx +  dy * dy);
	const double b  = 2*(dx * x  +  dy * y);
	const double c  =   (x  * x  +  y  * y)  -  R * R;
    
	const double q  = b*b - 4*a*c;
    
	if (q >= 0) {
      
	  const double u[] = {
	    (-b - sqrt(q)) / (2*a),
	    (-b + sqrt(q)) / (2*a)
	  };
      
	  for (int i = 0; i != 2; ++i) {
	
	    const double z = axis.getZ() + u[i] * axis.getDZ();
	
	    if (z >= this->getZmin() && z <= this->getZmax()) {
	      path[i] = u[i];
	    }
	  }
	}
      }

      return std::make_pair(path[0], path[1]);
    }

    
    /**
     * Read cylinder from input stream.
     *
     * \param  in          input stream
     * \param  cylinder    cylinder
     * \return             input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JCylinder3D& cylinder)
    {
      in >> static_cast<JCircle2D&>(cylinder);
      in >> cylinder.zmin >> cylinder.zmax;
      
      return in;
    }


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

      out << static_cast<const JCircle2D&>(cylinder);
      out << ' ';
      out << format << cylinder.zmin;
      out << ' ';
      out << format << cylinder.zmax;

      return out;
    }


    /**
     * Read cylinder from input.
     *
     * \param  in          reader
     * \param  cylinder    cylinder
     * \return             reader
     */
    friend inline JReader& operator>>(JReader& in, JCylinder3D& cylinder)
    {
      in >> static_cast<JCircle2D&>(cylinder);
      in >> cylinder.zmin >> cylinder.zmax;

      return in;
    }


    /**
     * Write cylinder to output.
     *
     * \param  out         writer
     * \param  cylinder    cylinder
     * \return             writer
     */
    friend inline JWriter& operator<<(JWriter& out, const JCylinder3D& cylinder)
    {
      out << static_cast<const JCircle2D&>(cylinder);
      out << cylinder.zmin << cylinder.zmax;

      return out;
    }
    

  protected:
    double zmin;
    double zmax;
  };
}

#endif