Consistently use small angle approximation in GModelSpatialEllipticalGauss and t…

…runcate ellipse at 3 sigma (#1561).

Consistently use small angle approximation in GModelSpatialEllipticalGauss and t…
…runcate ellipse at 3 sigma (#1561).
Jürgen Knödlseder
1 parent 2a49e210
Showing 2 changed files with 70 additions and 49 deletions
inst/cta/src/GCTAResponse_helpers.cpp
src/model/GModelSpatialEllipticalGauss.cpp
@@ -864,11 +864,15 @@ double cta_irf_elliptical_kern_omega::eval(const double&amp; 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
@@ -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(&amp;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&amp; 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&amp; 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&amp; 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&amp;  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&amp;  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&amp;  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&amp;  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&amp; 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;