JOmega3D.hh

Go to the documentation of this file.

#ifndef __JOMEGA3D__
#define __JOMEGA3D__

#include <cmath>
#include <vector>

#include "JTools/JConstants.hh"
#include "JTools/JAngleRange.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 JTOOLS::PI;
  using JTOOLS::JAngleRange;


  /**
   * Base class for direction set.
   */
  struct JOmega3D_t :
    public std::vector<JAngle3D> 
  {
    /**
     * 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.
   */
  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 step size [rad]
     */
    JOmega3D(const double grid) :
      JOmega3D_t()
    {
      configure(JAngle3D(0.0,0.0), JAngleRange(0.0,PI), grid);
    }



    /**
     * Constructor.
     *
     * \param  dir           principal direction
     * \param  theta         polar angle range [rad]
     * \param  grid          angular step size [rad]
     */
    JOmega3D(const JAngle3D&    dir,
             const JAngleRange& 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 step size [rad]
     */
    void configure(const JAngle3D&    dir,
                   const JAngleRange& theta,
                   const double       grid)
    {
      clear();

      // sanity checks

      double thetaMin = theta.getLowerLimit();
      double thetaMax = theta.getUpperLimit();

      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