JOmega3D.hh

Go to the documentation of this file.

#ifndef __JOMEGA3D__
#define __JOMEGA3D__
 
#include <cmath>
#include <utility>
#include <vector>
 
#include "JMath/JConstants.hh"
#include "JGeometry3D/JAngle3D.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JGeometry3D/JRotation3D.hh"
 
 
/**
 * \author mdejong
 */
 
namespace JGEOMETRY3D {}
namespace JPP { using namespace JGEOMETRY3D; }
 
namespace JGEOMETRY3D {
 
  using JMATH::PI;
 
 
  /**
   * Base class for direction set.
   */
  struct JOmega3D_t :
    public std::vector<JAngle3D> 
  {
    typedef std::pair<double, double>  range_type;
 
 
    /**
     * Get index of direction closest to given direction.
     *
     * \param  dir       direction
     * \return           index (-1 if error)
     */
    int find(const JAngle3D& 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) solid angle.\n
   * The set consists of approximately equally spaced directions.\n 
   * In this, the spacing in azimuth angle is zenith angle dependent. 
   */
  class JOmega3D :
    public JOmega3D_t
  {
  public:
    /**
     * Default constructor.
     */
    JOmega3D() :
      JOmega3D_t()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  dir           direction
     */
    JOmega3D(const JAngle3D& dir) :
      JOmega3D_t()
    {
      push_back(dir);
    }
 
 
    /**
     * Constructor.
     *
     * \param  grid          angular spacing   [rad]
     */
    JOmega3D(const double grid) :
      JOmega3D_t()
    {
      configure(JAngle3D(0.0,0.0), range_type(0.0,PI), grid);
    }
 
 
 
    /**
     * Constructor.
     *
     * \param  dir           principal direction
     * \param  theta         polar angle range [rad]
     * \param  grid          angular spacing   [rad]
     */
    JOmega3D(const JAngle3D&   dir,
             const range_type& theta,
             const double      grid) :
      JOmega3D_t()
    {
      configure(dir, theta, grid);
    }
 
 
    /**
     * Configure direction set.
     *
     * \param  dir           principal direction
     * \param  theta         polar angle range [rad]
     * \param  grid          angular spacing   [rad]
     */
    void configure(const JAngle3D&   dir,
                   const range_type& theta,
                   const double      grid)
    {
      clear();
 
      // sanity checks
 
      double thetaMin = theta.first;
      double thetaMax = theta.second;
 
      if (thetaMin < 0.0) { thetaMin = 0.0; }
      if (thetaMin > PI)  { thetaMin = PI; }
      if (thetaMax < 0.0) { thetaMax = 0.0; }
      if (thetaMax > PI)  { thetaMax = PI; }
 
      if (thetaMax > thetaMin) {
        
        const JRotation3D R(dir);
 
        const double rad  = thetaMax - thetaMin;
        const double bin  = rad / floor(rad/(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 = thetaMin; theta < thetaMax + 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;                                            // pole
 
          for (double phi = 0.0; phi < 2.0*PI - 0.5*step; phi += step) {
            push_back(JDirection3D(JAngle3D(theta,phi)).rotate_back(R));
          }
        }
      }
    }
  };
}
 
#endif