#ifndef __JOMEGA3D__ #define __JOMEGA3D__ #include #include #include #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 { typedef std::pair 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