diff --git a/inst/cta/src/GCTAResponse_helpers.cpp b/inst/cta/src/GCTAResponse_helpers.cpp index d4cc4c8..6365f87 100644 --- a/inst/cta/src/GCTAResponse_helpers.cpp +++ b/inst/cta/src/GCTAResponse_helpers.cpp @@ -864,11 +864,15 @@ double cta_irf_elliptical_kern_omega::eval(const double& omega) double arg2 = m_semimajor_rad * sinus; double r_ellipse = m_semiminor_rad * m_semimajor_rad / std::sqrt(arg1*arg1 + arg2*arg2); - std::cout << "*** WARNING: cta_irf_elliptical_kern_omega::eval"; + std::cout << "*** WARNING: cta_irf_elliptical_kern_omega::eval:"; std::cout << " zero model for (rho,omega)=("; std::cout << m_rho*gammalib::rad2deg << ","; std::cout << omega*gammalib::rad2deg << ")"; - std::cout << " rho-r_ellipse=" << (m_rho-r_ellipse) << " radians"; + std::cout << " semiminor=" << m_semiminor_rad; + std::cout << " semimajor=" << m_semimajor_rad; + std::cout << " posang=" << m_model->posangle() * gammalib::deg2rad; + std::cout << " rel_posang=" << diff_angle * gammalib::deg2rad; + std::cout << " rho/r_ellipse=" << (m_rho/r_ellipse); std::cout << std::endl; } #endif diff --git a/src/model/GModelSpatialEllipticalGauss.cpp b/src/model/GModelSpatialEllipticalGauss.cpp index 41c85cc..73239fc 100644 --- a/src/model/GModelSpatialEllipticalGauss.cpp +++ b/src/model/GModelSpatialEllipticalGauss.cpp @@ -37,7 +37,7 @@ /* __ Constants __________________________________________________________ */ namespace { - const double c_theta_max = 2.0; //!< semiaxis multiplied for theta_max() + const double c_theta_max = 3.0; //!< semiaxis multiplied for theta_max() const double c_max_exponent = 0.5 * c_theta_max * c_theta_max; const double c_fraction = 1.0 - std::exp(-c_max_exponent); } @@ -53,6 +53,7 @@ const GModelSpatialRegistry g_elliptical_gauss_registry(&g_elliptical_gauss_seed /* __ Macros _____________________________________________________________ */ /* __ Coding definitions _________________________________________________ */ +#define G_SMALL_ANGLE_APPROXIMATION //!< Use small angle approximation /* __ Debug definitions __________________________________________________ */ @@ -78,14 +79,14 @@ GModelSpatialEllipticalGauss::GModelSpatialEllipticalGauss(void) : /***********************************************************************//** - * @brief Elliptical Gauss constructor + * @brief Elliptical Gaussian constructor * - * @param[in] dir Sky position of gauss centre. + * @param[in] dir Centre of elliptical Gaussian. * @param[in] semimajor Semi-major axis (degrees). * @param[in] semiminor Semi-minor axis (degrees). * @param[in] posangle Position angle of semi-major axis (degrees). * - * Construct elliptical gauss model from sky position of the ellipse centre + * Construct elliptical Gaussian model from the ellipse centre direction * (@p dir), the @p semimajor and @p semiminor axes, and the position * angle (@p posangle). ***************************************************************************/ @@ -114,7 +115,7 @@ GModelSpatialEllipticalGauss::GModelSpatialEllipticalGauss(const GSkyDir& dir, * * @param[in] xml XML element. * - * Constructs elliptical gauss model by extracting information from an XML + * Constructs elliptical Gaussian model by extracting information from an XML * element. See the read() method for more information about the expected * structure of the XML element. ***************************************************************************/ @@ -135,7 +136,7 @@ GModelSpatialEllipticalGauss::GModelSpatialEllipticalGauss(const GXmlElement& xm /***********************************************************************//** * @brief Copy constructor * - * @param[in] model Elliptical gauss model. + * @param[in] model Elliptical Gaussian model. ***************************************************************************/ GModelSpatialEllipticalGauss::GModelSpatialEllipticalGauss(const GModelSpatialEllipticalGauss& model) : GModelSpatialElliptical(model) @@ -174,7 +175,7 @@ GModelSpatialEllipticalGauss::~GModelSpatialEllipticalGauss(void) * @brief Assignment operator * * @param[in] model Elliptical gauss model. - * @return Elliptical gauss model. + * @return Elliptical Gaussian model. ***************************************************************************/ GModelSpatialEllipticalGauss& GModelSpatialEllipticalGauss::operator=(const GModelSpatialEllipticalGauss& model) { @@ -207,7 +208,7 @@ GModelSpatialEllipticalGauss& GModelSpatialEllipticalGauss::operator=(const GMod ==========================================================================*/ /***********************************************************************//** - * @brief Clear elliptical gauss model + * @brief Clear elliptical Gaussian model ***************************************************************************/ void GModelSpatialEllipticalGauss::clear(void) { @@ -227,13 +228,13 @@ void GModelSpatialEllipticalGauss::clear(void) /***********************************************************************//** - * @brief Clone elliptical gauss model + * @brief Clone elliptical Gaussian model * - * @return Pointer to deep copy of elliptical gauss model. + * @return Pointer to deep copy of elliptical Gaussian model. ***************************************************************************/ GModelSpatialEllipticalGauss* GModelSpatialEllipticalGauss::clone(void) const { - // Clone elliptical gauss model + // Clone elliptical Gaussian model return new GModelSpatialEllipticalGauss(*this); } @@ -249,46 +250,48 @@ GModelSpatialEllipticalGauss* GModelSpatialEllipticalGauss::clone(void) const * * Evaluates the spatial component for an elliptical Gaussian source model. * - * The source model is an elliptical Gaussian function + * The elliptical Gaussian function is defined by * \f$S_{\rm p}(\theta, \phi | E, t)\f$, where * \f$\theta\f$ is the angular separation between elliptical Gaussian centre * and the actual location and \f$\phi\f$ the position angle with respect to * the model centre, counted counterclockwise from North. * - * The function \f$f(\theta, \phi)\f$ is given by + * The function \f$S_{\rm p}(\theta, \phi | E, t)\f$ is given by * * \f[ - * S_{\rm p}(\theta, \phi | E, t) = \left \{ - * \begin{array}{l l} - * {\tt m\_norm} - * & \mbox{if} \, \, \theta \le \theta_0 \\ - * \\ - * 0 & \mbox{else} - * \end{array} - * \right . + * S_{\rm p}(\theta, \phi | E, t) = {\tt m\_norm} \times + * \exp \left( -\frac{\theta^2}{2 r_{rm eff}} \right) * \f] * - * where \f$\theta_0\f$ is the effective radius of the ellipse on the sphere - * given by - * - * \f[\theta_0\ = - * \frac{ab}{\sqrt{b^2 \cos^2(\phi-\phi_0) + a^2 \sin^2(\phi-\phi_0)}}\f] + * where the effective ellipse radius \f$r_{rm eff}\f$ towards a given + * position angle is given by * - * and + * \f[ + * r_{rm eff} = \frac{ab} + * {\sqrt{\left( a \sin (\phi - \phi_0) \right)^2 + + * \sqrt{\left( b \cos (\phi - \phi_0) \right)^2} + * \f] + * and * \f$a\f$ is the semi-major axis of the ellipse, * \f$b\f$ is the semi-minor axis, and * \f$\phi_0\f$ is the position angle of the ellipse, counted * counterclockwise from North. + * \f${\tt m\_norm}\f$ is a normalization constant given by * - * The normalisation constant \f${\tt m\_norm}\f$ which is the inverse of the - * solid angle subtended by an ellipse is given by * \f[ * {\tt m\_norm} = \frac{1}{2 \pi \times a \times b} * \f] * * @warning - * The normalization of the elliptical Gaussian is only valid in the small - * angle approximation. + * The above formula for an elliptical Gaussian are accurate for small + * angles, with semimajor and semiminor axes below a few degrees. This + * should be fine for most use cases, but you should be aware that the + * ellipse will get distorted for larger angles. + * + * @warning + * For numerical reasons the elliptical Gaussian will be truncated for + * \f$\theta\f$ angles that correspond to 3 times the effective ellipse + * radius. ***************************************************************************/ double GModelSpatialEllipticalGauss::eval(const double& theta, const double& posangle, @@ -304,6 +307,18 @@ double GModelSpatialEllipticalGauss::eval(const double& theta, // Update precomputation cache update(); + // Perform computations to compute exponent + #if defined(G_SMALL_ANGLE_APPROXIMATION) + double rel_posangle = posangle - this->posangle() * gammalib::deg2rad; + double cosinus = std::cos(rel_posangle); + double sinus = std::sin(rel_posangle); + double arg1 = m_minor_rad * cosinus; + double arg2 = m_major_rad * sinus; + double r_ellipse = m_minor_rad * m_major_rad / + std::sqrt(arg1*arg1 + arg2*arg2); + double r_relative = theta/r_ellipse; + double exponent = 0.5*r_relative*r_relative; + #else // Perform computations double sinphi = std::sin(posangle); double cosphi = std::cos(posangle); @@ -314,10 +329,13 @@ double GModelSpatialEllipticalGauss::eval(const double& theta, double term3 = m_term3 * sinphi * cosphi; // Compute exponent - double exponent = theta * theta * (term1 + term2 + term3); + double exponent = theta * theta * (term1 + term2 + term3); + #endif - // Set value - value = m_norm * std::exp(-exponent); + // Set value if the exponent is within the truncation region + if (exponent <= c_max_exponent) { + value = m_norm * std::exp(-exponent); + } // Compile option: Check for NaN/Inf #if defined(G_NAN_CHECK) @@ -350,7 +368,7 @@ double GModelSpatialEllipticalGauss::eval(const double& theta, * @return Model value. * * Evaluates the function value. No gradient computation is implemented as - * Elliptical models will be convolved with the instrument response and thus + * elliptical models will be convolved with the instrument response and thus * require the numerical computation of the derivatives. * * See the eval() method for more information. @@ -373,10 +391,12 @@ double GModelSpatialEllipticalGauss::eval_gradients(const double& theta, * @param[in,out] ran Random number generator. * @return Sky direction. * - * Draws an arbitrary sky direction from the 2D ellipital Gaussian model. - * The drawn sky direction will be truncated to \f$2 \sigma\f$ offsets from - * the ellipse centre, where \f$\sigma\f$ is the effective Gaussian size - * in the direction of a specific position angle. + * Draws an arbitrary sky direction from the elliptical Gaussian model. + * + * @warning + * For numerical reasons the elliptical Gaussian will be truncated for + * \f$\theta\f$ angles that correspond to 3 times the effective ellipse + * radius. ***************************************************************************/ GSkyDir GModelSpatialEllipticalGauss::mc(const GEnergy& energy, const GTime& time, @@ -448,12 +468,10 @@ bool GModelSpatialEllipticalGauss::contains(const GSkyDir& dir, * * @return Returns maximum model radius. * - * Returns the maximum of \f$2\f$ semimajor() and \f$2\f$ semiminor() as + * Returns the maximum of \f$3\f$ semimajor() and \f$3\f$ semiminor() as * approximate edge of the Gaussian. This limit is of course arbitrary, but * allows to limit the integration region for response computation. The value - * of 2 has been determined by experiment. Ideally, a large value should - * be used to properly take into account the tails of the distribution, but - * this apparently leads to fit convergence problems (#1561). + * of 3 has been determined by experiment (#1561). ***************************************************************************/ double GModelSpatialEllipticalGauss::theta_max(void) const { @@ -707,7 +725,7 @@ void GModelSpatialEllipticalGauss::init_members(void) /***********************************************************************//** * @brief Copy class members * - * @param[in] model Elliptical gauss model. + * @param[in] model Elliptical Gaussian model. ***************************************************************************/ void GModelSpatialEllipticalGauss::copy_members(const GModelSpatialEllipticalGauss& model) { @@ -756,9 +774,6 @@ void GModelSpatialEllipticalGauss::free_members(void) * @warning * The normalization of the elliptical Gaussian is only valid in the small * angle approximation. - * - * @todo Generalize normalization to the sphere (for the moment only accurate - * in the small angle approximation). ***************************************************************************/ void GModelSpatialEllipticalGauss::update() const { @@ -792,6 +807,7 @@ void GModelSpatialEllipticalGauss::update() const } // endif: update required // Update chache if position angle changed + #if !defined(G_SMALL_ANGLE_APPROXIMATION) if (m_last_posangle != posangle()) { // Signal parameter changes @@ -823,6 +839,7 @@ void GModelSpatialEllipticalGauss::update() const m_term3 = 0.5 * (m_sin2pos / m_major2 - m_sin2pos / m_minor2); } // endif: something has changed + #endif // Return return;