#ifndef __JOMEGA3D__
#define __JOMEGA3D__

#include <cmath>
#include <vector>

#include "JTools/JConstants.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JGeometry3D/JRotation3D.hh"


namespace JGEOMETRY3D {

  namespace {
    using JTOOLS::PI;
  }


  /**
   * Type definition of direction set.
   */
  typedef std::vector<JDirection3D>  JOmega3D_t;


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


    /**
     * Constructor.
     *
     * \param  grid          angular step size [rad]
     */
    JOmega3D(const double grid) :
      JOmega3D_t()
    {
      configure(JDirection3D(0,0,1), PI, grid);
    }



    /**
     * Constructor.
     *
     * \param  dir           principal direction
     * \param  thetaMax      maximal polar angle [rad]
     * \param  grid          angular step size   [rad]
     */
    JOmega3D(const JDirection3D& dir,
	     const double        thetaMax,
	     const double        grid) :
      JOmega3D_t()
    {
      configure(dir, thetaMax, grid);
    }


    /**
     * Configure direction set.
     *.
     * \param  dir           principal direction
     * \param  thetaMax      maximal polar angle [rad]
     * \param  grid          angular step size   [rad]
     */
    void configure(const JDirection3D& dir,
		   const double         thetaMax,
		   const double         grid)
    {
      clear();

      // sanity checks

      double theta_max = thetaMax; 

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

      const JRotation3D R(dir);

      const double bin  = PI / floor(PI/(1.4*grid) + 0.5);          // polar angle step size
      const double unit = 2.0 * sqrt(1.0 - cos(grid));              // azimuth angle unit step size

      for (double theta = 0; theta < theta_max + 0.5*bin; theta += bin) {

	double step;

        if (theta > 0.5*bin && PI - theta > 0.5*bin)
          step = PI / floor(PI*sin(theta)/unit + 0.5);              // polar angle dependent step size
        else
          step = 2.0*PI;

        for (double phi = 0; phi < 2.0*PI - 0.5*step; phi += step)
          push_back(JDirection3D(theta,phi).rotate_back(R));
      }
    }
  };
}

#endif