#ifndef __JOMEGA2D__
#define __JOMEGA2D__

#include <cmath>
#include <vector>

#include "JMath/JConstants.hh"

#include "JGeometry2D/JAngle2D.hh"
#include "JGeometry2D/JRotation2D.hh"


/**
 * \author mdejong
 */

namespace JGEOMETRY2D {}
namespace JPP { using namespace JGEOMETRY2D; }

namespace JGEOMETRY2D {

  using JMATH::PI;


  /**
   * Base class for direction set.
   */
  struct JOmega2D_t :
    public std::vector<JAngle2D> 
  {
    /**
     * Get index of direction closest to given direction.
     *
     * \param  dir       direction
     * \return           index (-1 if error)
     */
    int find(const JAngle2D& dir) const
    {
      double dot   = -1.0;
      int    index = -1;

      for (const_iterator i = this->begin(); i != this->end(); ++i) {

        const double x = dir.getDot(*i);
        
        if (x > dot) {
          dot   = x;
          index = std::distance(this->begin(), i);
        }
      }

      return index;
    }
  };


  /**
   * Direction set covering (part of) circle.
   */
  class JOmega2D :
    public JOmega2D_t
  {
  public:
    /**
     * Default constructor.
     */
    JOmega2D() :
      JOmega2D_t()
    {}


    /**
     * Constructor.
     *
     * \param  grid          angular spacing   [rad]
     */
    JOmega2D(const double grid) :
      JOmega2D_t()
    {
      configure(JAngle2D(0.0), PI, grid);
    }



    /**
     * Constructor.
     *
     * \param  dir           principal direction
     * \param  phiMax        maximal angle     [rad]
     * \param  grid          angular spacing   [rad]
     */
    JOmega2D(const JAngle2D& dir,
	     const double    phiMax,
	     const double    grid) :
      JOmega2D_t()
    {
      configure(dir, phiMax, grid);
    }


    /**
     * Configure direction set.
     *
     * \param  dir           principal direction
     * \param  phiMax        maximal angle     [rad]
     * \param  grid          angular spacing   [rad]
     */
    void configure(const JAngle2D& dir,
		   const double    phiMax,
		   const double    grid)
    {
      clear();

      // sanity checks

      double phi_min = 0.0;
      double phi_max = phiMax;

      if (phi_max < 0.0) phi_max = 0.0;
      if (phi_max > PI)  phi_max = PI;

      if (phi_max > phi_min) {

	const double rad = phi_max - phi_min; 
	const double bin = rad / floor(rad/grid + 0.5);             // angle step size

	for (double phi = phi_min; phi < phi_max + 0.5*bin; phi += bin) {
	  
	  if        (phi      < 0.5*bin) {
	    push_back(dir);
	  } else if (PI - phi < 0.5*bin) {
	    push_back(dir + PI);
	  } else {
	    push_back(dir + phi);
	    push_back(dir - phi);
	  }
	}
      }
    }
  };
}

#endif