#ifndef __JCIRCLE2D__
#define __JCIRCLE2D__

#include <iterator>
#include <cmath>

#include "JIO/JSerialisable.hh"
#include "JGeometry2D/JPoint2D.hh"


namespace JGEOMETRY2D {

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


  /**
   * Data structure for circle in two dimensions.
   */
  class JCircle2D :
    public JPoint2D
  {
  public:
    /**
     * Default constructor.
     */
    JCircle2D() :
      JPoint2D(),
      __r(0.0)
    {}


    /**
     * Constructor.
     *
     * \param  point    JPoint2D
     * \param  r        radius
     */
    JCircle2D(const JPoint2D& point,
	      const double    r) :
      JPoint2D(point),
      __r(r)
    {}


    /**
     * Constructor.
     *
     * Determines smallest enclosing circle for 2 points.
     *
     * \param  p0       JPoint2D
     * \param  p1       JPoint2D
     */
    JCircle2D(const JPoint2D& p0,
	      const JPoint2D& p1) :
      JPoint2D(),
      __r(0.0)
    {
      __x = 0.5 * (p0.getX() + p1.getX());
      __y = 0.5 * (p0.getY() + p1.getY());

      const double dx = p0.getX() - this->getX();
      const double dy = p0.getY() - this->getY();
	
      __r = sqrt(dx*dx + dy*dy);
    }


    /**
     * Constructor.
     *
     * Determines smallest enclosing circle for 3 points.
     *
     * \param  p0       JPoint2D
     * \param  p1       JPoint2D
     * \param  p2       JPoint2D
     */
    JCircle2D(const JPoint2D& p0,
	      const JPoint2D& p1,
	      const JPoint2D& p2) :
      JPoint2D(),
      __r(0.0)
    {
      const double x0 = p2.getX() - p1.getX();
      const double x1 = p0.getX() - p2.getX();
      const double x2 = p1.getX() - p0.getX();

      const double y0 = p1.getY() - p2.getY();
      const double y1 = p2.getY() - p0.getY();
      const double y2 = p0.getY() - p1.getY();

      const double D  = 2.0 * (p0.getX()*y0 + p1.getX()*y1 + p2.getX()*y2);

      if (D != 0.0) {

	const double a = p0.getX()*p0.getX()  +  p0.getY()*p0.getY();
	const double b = p1.getX()*p1.getX()  +  p1.getY()*p1.getY();
	const double c = p2.getX()*p2.getX()  +  p2.getY()*p2.getY();

	__x = (a*y0 + b*y1 + c*y2) / D;
	__y = (a*x0 + b*x1 + c*x2) / D;
	
	const double dx = p0.getX() - this->getX();
	const double dy = p0.getY() - this->getY();
	
	__r = sqrt(dx*dx + dy*dy);

      } else {

	*this = JCircle2D(p0, p1);

	const JCircle2D c1(p1, p2);
	const JCircle2D c2(p0, p2);

	if (c1.getRadius() > this->getRadius()) *this = c1;
	if (c2.getRadius() > this->getRadius()) *this = c2;
      }
    }


    /**
     * Constructor.
     *
     * Determines smallest enclosing circle for any number of points.
     *
     * \param  __begin  begin of data
     * \param  __end    end   of data
     */
    template<class T>
    JCircle2D(T __begin,
	      T __end) :
      JPoint2D(),
      __r(0.0)
    {
      if (std::distance(__begin,__end) == 1) {

	__x = __begin->getX();
	__y = __begin->getY();
	__r = 0.0;

	return;
      }	

      if (std::distance(__begin,__end) > 1) {

	T i = __begin;
	T j = __begin;

	++j;

	*this = JCircle2D(*i, *j);

	while (++j != __end) {
	  if (getDistance(*j) > getRadius())
	    configure(__begin, j, *j);
	}

	return;
      }
    }


    /**
     * Read JCircle2D from input.
     *
     * \param  in         JReader
     * \param  circle     JCircle2D
     * \return            JReader
     */
    friend inline JReader& operator>>(JReader& in, JCircle2D& circle)
    {
      in >> static_cast<JPoint2D&>(circle);
      in >> circle.__r;

      return in;
    }


    /**
     * Write JCircle2D to output.
     *
     * \param  out        JWriter
     * \param  circle     JCircle2D
     * \return            JWriter
     */
    friend inline JWriter& operator<<(JWriter& out, const JCircle2D& circle)
    {
      out << static_cast<const JPoint2D&>(circle);
      out << circle.__r;

      return out;
    }


    /**
     * Get radius.
     *
     * \return          radius
     */
    double getRadius() const { return __r; }


  protected:
    double __r;


  private:
    /**
     * Determine smallest enclosing circle.
     *
     * \param  __begin  begin of data
     * \param  __end    end   of data
     * \param  p1       point on circle
     */
    template<class T>
    void configure(T __begin, 
		   T __end, 
		   const typename std::iterator_traits<T>::reference p1)
    {
      *this = JCircle2D(*__begin, p1);
      
      for (T i = __begin; ++i != __end; ) {
	if (this->getDistance(*i) > this->getRadius())
	  configure(__begin, i, *i, p1);
      }
    }


    /**
     * Determine smallest enclosing circle.
     *
     * \param  __begin  begin of data
     * \param  __end    end   of data
     * \param  p1       point on circle
     * \param  p2       point on circle
     */
    template<class T>
    void configure(T __begin, 
		   T __end, 
		   const typename std::iterator_traits<T>::reference p1, 
		   const typename std::iterator_traits<T>::reference p2)
    {
      *this = JCircle2D(*__begin, p1, p2);

      for (T i = __begin; ++i != __end; ) {
	if (this->getDistance(*i) > this->getRadius())
	  *this = JCircle2D(*i, p1, p2);
      }
    }
  };
}

#endif