#ifndef __JOMEGA2D__
#define __JOMEGA2D__

#include <cmath>
#include <vector>

#include "JTools/JConstants.hh"
#include "JGeometry2D/JDirection2D.hh"
#include "JGeometry2D/JRotation2D.hh"


namespace JGEOMETRY2D {

  namespace {
    using JTOOLS::PI;
  }


  /**
   * Type definition of direction set.
   */
  typedef std::vector<JDirection2D>  JOmega2D_t;


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


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



    /**
     * Constructor.
     *
     * \param  dir           principal direction
     * \param  phiMax        maximal angle     [rad]
     * \param  grid          angular step size [rad]
     */
    JOmega2D(const JDirection2D& 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 step size [rad]
     */
    void configure(const JDirection2D& dir,
		   const double        phiMax,
		   const double        grid)
    {
      clear();

      // sanity checks

      double phi_max = phiMax;

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

      const JRotation2D R(dir);

      const double bin  = PI / floor(PI/grid + 0.5);             // angle step size

      for (double phi = 0.0; phi < phi_max + 0.5*bin; phi += bin) {

	if (phi > 0.5*bin && PI - phi > 0.5*bin) {
	  push_back(JDirection2D(+phi).rotate_back(R));
	  push_back(JDirection2D(-phi).rotate_back(R));
	} else
	  push_back(JDirection2D(phi).rotate_back(R));
      }
    }
  };
}

#endif