Remove small angle approximation (always compute in spherical system).

Jürgen Knödlseder
1 parent 86b32056
Showing 5 changed files with 42 additions and 135 deletions
ChangeLog
NEWS
include/GModelSpatialRadialShell.hpp
pyext/GModelSpatialRadialShell.i
src/model/GModelSpatialRadialShell.cpp
-2015-07-03  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
+2015-08-30  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
  
         * Version 0.11.0 released
         =========================
  
+        Remove small angle approximation from GModelSpatialRadialShell
         Rename GSkymap to GSkyMap
         Add binary operators to GSkymap class
         GSkymap::operator/=(GSkymap) operator now ignores division by zero
@@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 0.11.0
  
 Author(s) : Juergen Knoedlseder
  
-3 July 2015
+30 August 2015
  
  
 1. Introduction
@@ -20,6 +20,8 @@ optimizer.
 The GSkymap class has been renamed to GSkyMap to be consistent with the other
 classes in the sky module.
  
+The small angle approximation has been removed from GModelSpatialRadialShell.
+
  
 3. Configuration
 -----------------
@@ -62,6 +64,9 @@ The algorithm used to draw random numbers in GModelSpatialDiffuseConst::mc()
 was wrong (sky directions were not distributed isotropically on the sky). This
 has now been fixed (#1480).
  
+The small angle approximation has been removed from GModelSpatialRadialShell
+as this complexifies the class and the approximation is not really needed.
+
  
 9. Numerics module
 ------------------
@@ -61,8 +61,7 @@ public:
     GModelSpatialRadialShell(void);
     GModelSpatialRadialShell(const GSkyDir& dir,
                              const double&  radius,
-                             const double&  width,
-                             const bool&    small_angle = true);
+                             const double&  width);
     explicit GModelSpatialRadialShell(const GXmlElement& xml);
     GModelSpatialRadialShell(const GModelSpatialRadialShell& model);
     virtual ~GModelSpatialRadialShell(void);
@@ -92,12 +91,10 @@ public:
     virtual std::string               print(const GChatter& chatter = NORMAL) const;
  
     // Other methods
-    double      radius(void) const;
-    double      width(void) const;
-    const bool& small_angle(void) const;
-    void        radius(const double& radius);
-    void        width(const double& width);
-    void        small_angle(const bool& small_angle);
+    double radius(void) const;
+    double width(void) const;
+    void   radius(const double& radius);
+    void   width(const double& width);
  
 protected:
     // Protected methods
@@ -111,15 +108,14 @@ protected:
     // Protected members
     GModelPar       m_radius;        //!< Inner shell radius (deg)
     GModelPar       m_width;         //!< Shell thickness (deg)
-    bool            m_small_angle;   //!< Use small angle approximate formulae
  
     // Cached members used for pre-computations
     mutable double  m_last_radius;   //!< Last shell radius (deg)
     mutable double  m_last_width;    //!< Last shell width (deg)
     mutable double  m_theta_in;      //!< Inner shell radius (rad)
-    mutable double  m_x_in;          //!< m_theta_in^2 or sin(m_theta_in)^2
+    mutable double  m_x_in;          //!< sin(m_theta_in)^2
     mutable double  m_theta_out;     //!< Outer shell radius (rad)
-    mutable double  m_x_out;         //!< m_theta_out^2 or sin(m_theta_out)^2
+    mutable double  m_x_out;         //!< sin(m_theta_out)^2
     mutable double  m_norm;          //!< Shell normalization
 };
  
@@ -207,35 +203,4 @@ void GModelSpatialRadialShell::width(const double&amp; width)
     return;
 }
  
-
-/***********************************************************************//**
- * @brief Return shell approximation flag
- *
- * @return True if small angle approximation is taken.
- *
- * Signals if the small angle approximation is taken.
- ***************************************************************************/
-inline
-const bool& GModelSpatialRadialShell::small_angle(void) const
-{
-    return (m_small_angle);
-}
-
-
-/***********************************************************************//**
- * @brief Set shell approximation flag
- *
- * @param[in] small_angle Shell approximation flag.
- *
- * Sets the shell approximation flag. True specifies that the small angle
- * approximation shall be taken, false specifies that the correct spherical
- * computation shall be used.
- ***************************************************************************/
-inline
-void GModelSpatialRadialShell::small_angle(const bool& small_angle)
-{
-    m_small_angle = small_angle;
-    return;
-}
-
 #endif /* GMODELSPATIALRADIALSHELL_HPP */
@@ -40,8 +40,7 @@ public:
     GModelSpatialRadialShell(void);
     GModelSpatialRadialShell(const GSkyDir& dir,
                              const double&  radius,
-                             const double&  width,
-                             const bool&    small_angle = true);
+                             const double&  width);
     explicit GModelSpatialRadialShell(const GXmlElement& xml);
     GModelSpatialRadialShell(const GModelSpatialRadialShell& model);
     virtual ~GModelSpatialRadialShell(void);
@@ -69,10 +68,8 @@ public:
     // Other methods
     double  radius(void) const;
     double  width(void) const;
-    bool    small_angle(void) const;
     void    radius(const double& radius);
     void    width(const double& width);
-    void    small_angle(const bool& small_angle);
 };
  
  
@@ -77,21 +77,13 @@ GModelSpatialRadialShell::GModelSpatialRadialShell(void) : GModelSpatialRadial()
  * @param[in] dir Sky position of shell centre.
  * @param[in] radius Inner shell radius (degrees).
  * @param[in] width Shell width (degrees).
- * @param[in] small_angle Use small angle approximation (default to true).
  *
  * Constructs the shell model from the shell centre (@p dir), the inner
- * shell @p radius, and the shell @p width. The optional parameters
- * @p small_angle specifies whether the small angle approximation
- * should be used or not. The small angle approximation is (a little) faster,
- * but incorrect for shells larger than a few degrees.
- *
- * @todo Measure speed difference between small angle approximation and
- *       full computation.
+ * shell @p radius, and the shell @p width.
  ***************************************************************************/
 GModelSpatialRadialShell::GModelSpatialRadialShell(const GSkyDir& dir,
                                                    const double&  radius,
-                                                   const double&  width,
-                                                   const bool&    small_angle) :
+                                                   const double&  width) :
                           GModelSpatialRadial()
 {
     // Initialise members
@@ -101,7 +93,6 @@ GModelSpatialRadialShell::GModelSpatialRadialShell(const GSkyDir&amp; dir,
     this->dir(dir);
     this->radius(radius);
     this->width(width);
-    this->small_angle(small_angle);
  
     // Return
     return;
@@ -249,28 +240,7 @@ GModelSpatialRadialShell* GModelSpatialRadialShell::clone(void) const
  * model is a radial function \f$f(\theta)\f$, where \f$\theta\f$ is the
  * angular separation between shell centre and the actual location.
  *
- * In the small angle approximation,
- * \f[
- * S_{\rm p}(\vec{p} | E, t) = {\tt m\_norm} \left \{
- *  \begin{array}{l l}
- *     \displaystyle
- *     \sqrt{ \theta_{\rm out}^2 - \theta^2 } -
- *     \sqrt{ \theta_{\rm in}^2  - \theta^2 }
- *     & \mbox{if $\theta \le \theta_{\rm in}$} \\
- *     \\
- *    \displaystyle
- *     \sqrt{ \theta_{\rm out}^2 - \theta^2 }
- *     & \mbox{if $\theta_{\rm in} < \theta \le \theta_{\rm out}$} \\
- *     \\
- *    \displaystyle
- *    0 & \mbox{if $\theta > \theta_{\rm out}$}
- *  \end{array}
- *  \right .
- * \f]
- * where \f${\tt m\_norm}\f$ is a normalization constant (see the update()
- * method).
- *
- * In the general form that is also valid for large angles
+ * The function is given by
  * \f[
  * S_{\rm p}(\vec{p} | E, t) = {\tt m\_norm} \left \{
  *  \begin{array}{l l}
@@ -299,28 +269,20 @@ double GModelSpatialRadialShell::eval(const double&amp;  theta,
     // Update precomputation cache
     update();
  
-    // Set x appropriately for the small angle approximation or not
-    double x;
-    if (m_small_angle) {
-        x = theta * theta;
-    }
-    else {
-        x  = std::sin(theta);
-        x *= x;
-    }
+    // Set x appropriately
+    double x = std::sin(theta);
+    x *= x;
  
-    // Compute value
-    double value = 0.0;
+    // Compute result
+    double result = 0.0;
     if (x < m_x_out) {
-        value = std::sqrt(m_x_out - x);
+        result = std::sqrt(m_x_out - x);
         if (x < m_x_in) {
-            value -= std::sqrt(m_x_in - x);
+            result -= std::sqrt(m_x_in - x);
         }
+        result *= m_norm;
     }
  
-    // Normalise value
-    double result = m_norm * value;
-
     // Compile option: Check for NaN/Inf
     #if defined(G_NAN_CHECK)
     if (gammalib::is_notanumber(result) || gammalib::is_infinite(result)) {
@@ -328,11 +290,9 @@ double GModelSpatialRadialShell::eval(const double&amp;  theta,
         std::cout << "(theta=" << theta << "): NaN/Inf encountered";
         std::cout << " (result=" << result;
         std::cout << ", m_norm=" << m_norm;
-        std::cout << ", value=" << value;
         std::cout << ", x=" << x;
         std::cout << ", m_x_out=" << m_x_out;
         std::cout << ", m_x_in=" << m_x_in;
-        std::cout << ", m_small_angle=" << m_small_angle;
         std::cout << ")" << std::endl;
     }
     #endif
@@ -381,17 +341,15 @@ GSkyDir GModelSpatialRadialShell::mc(const GEnergy&amp; energy,
  
     // Simulate offset from photon arrival direction
     #if defined(G_DEBUG_MC)
-    int    n_samples = 0;
+    int    n_samples     = 0;
     #endif
-    double theta_max     = this->theta_max();
-    double sin_theta_max = std::sin(theta_max);
-    double u_max         = (m_small_angle) ? m_norm * sin_theta_max * theta_max
-                                           : m_norm * sin_theta_max * sin_theta_max;
+    double sin_theta_max = std::sin(m_theta_out);
+    double u_max         = m_norm * sin_theta_max * sin_theta_max;
     double value         = 0.0;
     double u             = 1.0;
     double theta         = 0.0;
     do {
-        theta = ran.uniform() * theta_max;
+        theta = ran.uniform() * m_theta_out;
         value = eval(theta, energy, time) * std::sin(theta);
         u     = ran.uniform() * u_max;
         #if defined(G_DEBUG_MC)
@@ -684,9 +642,6 @@ void GModelSpatialRadialShell::init_members(void)
     m_pars.push_back(&m_radius);
     m_pars.push_back(&m_width);
  
-    // Initialise other members
-    m_small_angle = false;
-
     // Initialise precomputation cache. Note that zero values flag
     // uninitialised as a zero radius and width shell is not meaningful
     m_last_radius = 0.0;
@@ -716,7 +671,6 @@ void GModelSpatialRadialShell::copy_members(const GModelSpatialRadialShell&amp; mode
     // Copy members
     m_radius      = model.m_radius;
     m_width       = model.m_width;
-    m_small_angle = model.m_small_angle;
  
     // Copy precomputation cache
     m_last_radius = model.m_last_radius;
@@ -753,13 +707,6 @@ void GModelSpatialRadialShell::free_members(void)
  * m_theta_in contains the inner shell radius in radians, while m_theta_out
  * contains the outer shell radius in radians.
  *
- * In the small angle approximation, 
- * \f${\tt m\_x\_in} = {\tt m\_theta\_in}^2\f$,
- * \f${\tt m\_x\_out} = {\tt m\_theta\_out}^2\f$, and
- * \f[{\tt m\_norm} = \frac{3}{2 \pi}
- *    \frac{1}{\left( {\tt m\_theta\_out}^3 - {\tt m\_theta\_in}^3 \right)}\f].
- *
- * In the general case that is also valid for large angles,
  * \f${\tt m\_x\_in} = \sin^2 {\tt m\_theta\_in}\f$,
  * \f${\tt m\_x\_out} = \sin^2 {\tt m\_theta\_out}\f$, and
  * \f[{\tt m\_norm} = \frac{1}{2 \pi}
@@ -786,25 +733,17 @@ void GModelSpatialRadialShell::update() const
         m_last_width  = width();
  
         // Perform precomputations
-        m_theta_in   = radius() * gammalib::deg2rad;
-        m_theta_out  = (radius() + width()) * gammalib::deg2rad;
-        if (m_small_angle) {
-            m_x_in       = m_theta_in  * m_theta_in;
-            m_x_out      = m_theta_out * m_theta_out;
-            double denom = c1 * (m_x_out*m_theta_out - m_x_in*m_theta_in);
-            m_norm       = (denom > 0.0) ? 1.0 / denom : 0.0;
-        } 
-        else {
-            double sin_theta_in  = std::sin(m_theta_in);
-            double sin_theta_out = std::sin(m_theta_out);
-            double term1         = (f1(m_theta_out) - f1(m_theta_in)) * c2;
-            double term2         = f2(m_theta_out);
-            double term3         = f2(m_theta_in);
-            double denom         = gammalib::twopi * (term1 + term2 - term3);
-            m_norm               = (denom > 0.0) ? 1.0 / denom : 0.0;
-            m_x_in               = sin_theta_in*sin_theta_in;
-            m_x_out              = sin_theta_out * sin_theta_out;
-        }
+        m_theta_in           = radius() * gammalib::deg2rad;
+        m_theta_out          = (radius() + width()) * gammalib::deg2rad;
+        double sin_theta_in  = std::sin(m_theta_in);
+        double sin_theta_out = std::sin(m_theta_out);
+        double term1         = (f1(m_theta_out) - f1(m_theta_in)) * c2;
+        double term2         = f2(m_theta_out);
+        double term3         = f2(m_theta_in);
+        double denom         = gammalib::twopi * (term1 + term2 - term3);
+        m_norm               = (denom > 0.0) ? 1.0 / denom : 0.0;
+        m_x_in               = sin_theta_in  * sin_theta_in;
+        m_x_out              = sin_theta_out * sin_theta_out;
  
     } // endif: update required