Add GSkyMap::flux(GRegion(s)) and GSkyMap::solidangle(GRegion(s)) methods (#1785)

Jürgen Knödlseder
1 parent 2659914b
Showing 8 changed files with 171 additions and 15 deletions
ChangeLog
NEWS
README.md
doc/source/admin/release_history/2.0.rst
include/GSkyMap.hpp
pyext/GSkyMap.i
src/sky/GSkyMap.cpp
test/test_GSky.cpp
-2021-03-04
+2021-03-08
         * Version 2.0.0 released
         ========================
@@ -40,7 +40,8 @@
         Implement analytical gradients for the radial Gaussian model in a
           stacked analysis (#3203)
         Add GSkyModel flux_error() and eflux_error() methods (#3074)
-        Remove COMPTEL internal response caching (#2973).
+        Remove COMPTEL internal response caching (#2973)
+        Add GSkyMap::flux(GRegion(s)) and GSkyMap::solidangle(GRegion(s)) methods (#1785)
         Implement pre-computation cache in GModelSpatialGauss (#1300)
 New Features and Important Changes in GammaLib 2.0.0
-4 March 2021
+8 March 2021
 1. Introduction
@@ -60,6 +60,10 @@ The following methods have been added:
 - GModelSpectralFunc::intensity()
 - GModelSpectralFunc::save()
 - GSkyMap::correlate()
+- GSkyMap::flux(GSkyRegion&, int&)
+- GSkyMap::flux(GSkyRegions&, int&)
+- GSkyMap::solidangle(GSkyRegion&)
+- GSkyMap::solidangle(GSkyRegions&)
 - GOptimizerLM::npars()
 - GOptimizerLM::nfree()
 - GXmlElement(std::string&, int&)
@@ -223,6 +227,13 @@ method now returns an instance of GSkyDirs (#3324).
 Added GSkyRegionRect class to implement a rectangular region on the sky (#3248).
+Added GSkyMap::flux(GSkyRegion&, int&) and GSkyMap::flux(GSkyRegions&, int&) methods
+to compute the total flux within a given sky region or a list of sky regions. The flux
+returns the sum of the intensity times the solidangle for all pixels that are contained
+within the sky region(s). Containment is tested using the GSkyMap::contains() methods.
+Equivalent methods GSkyMap::solidangle(GSkyRegion&) and GSkyMap::solidangle(GSkyRegions&)
+were added to compute the total solid angle within specified sky region(s) (#1785).
+
 13. Support module
 ------------------
 GammaLib information
 ====================
-* Version: 2.0.0.dev (4 March 2021)
+* Version: 2.0.0.dev (8 March 2021)
 [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=gammalib-integrate-os)](https://cta-jenkins.irap.omp.eu/job/gammalib-integrate-os/)
 .. _2.0:
-GammaLib 2.0 release (3 March 2021)
+GammaLib 2.0 release (8 March 2021)
 ===================================
 Overview
@@ -56,6 +56,8 @@ Improvements
 * [`3203 <https://cta-redmine.irap.omp.eu/issues/3203>`_] -
   Implement vectorised response computation and analytical gradients for
   the radial Gaussian model in a stacked analysis
+* [`1785 <https://cta-redmine.irap.omp.eu/issues/1785>`_] -
+  Add ``GSkyMap::flux(GRegion)``, ``GSkyMap::flux(GRegions)``, ``GSkyMap::solidangle(GRegion)`` and ``GSkyMap::solidangle(GRegions)`` methods
 * [`1300 <https://cta-redmine.irap.omp.eu/issues/1300>`_] -
   Implement pre-computation cache in GModelSpatialGauss
 /***************************************************************************
  *                       GSkyMap.hpp - Sky map class                       *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2020 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2021 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -162,9 +162,13 @@ public:
     GNdarray                counts(void) const;
     double                  flux(const int& index, const int& map = 0) const;
     double                  flux(const GSkyPixel& pixel, const int& map = 0) const;
+    double                  flux(const GSkyRegion& region, const int& map = 0) const;
+    double                  flux(const GSkyRegions& regions, const int& map = 0) const;
     GNdarray                flux(void) const;
     double                  solidangle(const int& index) const;
     double                  solidangle(const GSkyPixel& pixel) const;
+    double                  solidangle(const GSkyRegion& region) const;
+    double                  solidangle(const GSkyRegions& regions) const;
     bool                    contains(const GSkyDir& dir) const;
     bool                    contains(const GSkyPixel& pixel) const;
     bool                    overlaps(const GSkyRegion& region) const;
 /***************************************************************************
  *                         GSkyMap.i - Sky map class                       *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2020 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2021 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -97,9 +97,13 @@ public:
     GNdarray                counts(void);
     double                  flux(const int& index, const int& map = 0) const;
     double                  flux(const GSkyPixel& pixel, const int& map = 0) const;
+    double                  flux(const GSkyRegion& region, const int& map = 0) const;
+    double                  flux(const GSkyRegions& regions, const int& map = 0) const;
     GNdarray                flux(void);
     double                  solidangle(const int& index) const;
     double                  solidangle(const GSkyPixel& pixel) const;
+    double                  solidangle(const GSkyRegion& region) const;
+    double                  solidangle(const GSkyRegions& regions) const;
     bool                    contains(const GSkyDir& dir) const;
     bool                    contains(const GSkyPixel& pixel) const;
     bool                    overlaps(const GSkyRegion& region) const;
 /***************************************************************************
  *                       GSkyMap.cpp - Sky map class                       *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2020 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2021 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -417,7 +417,7 @@ GSkyMap&amp; GSkyMap::operator+=(const GSkyMap&amp; map)
             } // endfor: looped over all layers
         } // endfor: looped over all pixels
-    
+
     }
     // ... otherwise use interpolation method and add map
@@ -507,7 +507,7 @@ GSkyMap&amp; GSkyMap::operator-=(const GSkyMap&amp; map)
             } // endfor: looped over all layers
         } // endfor: looped over all pixels
-    
+
     }
     // ... otherwise use interpolation method and subtract map
@@ -597,7 +597,7 @@ GSkyMap&amp; GSkyMap::operator*=(const GSkyMap&amp; map)
             } // endfor: looped over all layers
         } // endfor: looped over all pixels
-    
+
     }
     // ... otherwise use interpolation method and multiply map
@@ -701,7 +701,7 @@ GSkyMap&amp; GSkyMap::operator/=(const GSkyMap&amp; map)
             } // endfor: looped over all layers
         } // endfor: looped over all pixels
-    
+
     }
     // ... otherwise use interpolation method and divide map
@@ -1724,6 +1724,66 @@ double GSkyMap::flux(const GSkyPixel&amp; pixel, const int&amp; map) const
 /***********************************************************************//**
+ * @brief Returns flux within sky region
+ *
+ * @param[in] region Sky region.
+ * @param[in] map Map index [0,...,nmaps()-1].
+ * @return Flux in sky region.
+ *
+ * The method returns the total flux in all sky map pixels that fall within
+ * the specified region. Containment is tested using the
+ * GSkyRegion::contains() method. Sky map values are multiplied by the solid
+ * angle of the pixel.
+ ***************************************************************************/
+double GSkyMap::flux(const GSkyRegion& region, const int& map) const
+{
+    // Initialise flux
+    double flux = 0.0;
+
+    // Loop over all sky pixels and sum flux of all pixels that fall within
+    // the sky region
+    for (int i = 0; i < m_num_pixels; ++i) {
+        if (region.contains(inx2dir(i))) {
+            flux += (*this)(i,map) * solidangle(i);
+        }
+    }
+
+    // Return flux
+    return flux;
+}
+
+
+/***********************************************************************//**
+ * @brief Returns flux within sky regions
+ *
+ * @param[in] regions Sky regions.
+ * @param[in] map Map index [0,...,nmaps()-1].
+ * @return Flux in sky region.
+ *
+ * The method returns the total flux in all sky map pixels that fall within
+ * the specified regions. Containment is tested using the
+ * GSkyRegions::contains() method. Sky map values are multiplied by the solid
+ * angle of the pixel.
+ ***************************************************************************/
+double GSkyMap::flux(const GSkyRegions& regions, const int& map) const
+{
+    // Initialise flux
+    double flux = 0.0;
+
+    // Loop over all sky pixels and sum flux of all pixels that fall within
+    // the sky region
+    for (int i = 0; i < m_num_pixels; ++i) {
+        if (regions.contains(inx2dir(i))) {
+            flux += (*this)(i,map) * solidangle(i);
+        }
+    }
+
+    // Return flux
+    return flux;
+}
+
+
+/***********************************************************************//**
  * @brief Returns array with total flux for sky maps
  *
  * @return Array of pixel flux sums.
@@ -1734,7 +1794,7 @@ double GSkyMap::flux(const GSkyPixel&amp; pixel, const int&amp; map) const
 GNdarray GSkyMap::flux(void) const
 {
     // Initialise output array
-    GNdarray  flux_array = GNdarray(m_num_maps);
+    GNdarray flux_array = GNdarray(m_num_maps);
     // Loop over maps and store sums in output array
     for (int i = 0; i < m_num_maps; ++i) {
@@ -1836,6 +1896,62 @@ double GSkyMap::solidangle(const GSkyPixel&amp; pixel) const
 /***********************************************************************//**
+ * @brief Returns solid angle within sky region
+ *
+ * @param[in] region Sky region.
+ * @return Solid angle within sky region.
+ *
+ * The method returns the total solid angle of all sky map pixels that fall
+ * within the specified region. Containment is tested using the
+ * GSkyRegion::contains() method.
+ ***************************************************************************/
+double GSkyMap::solidangle(const GSkyRegion& region) const
+{
+    // Initialise flux
+    double solidangle = 0.0;
+
+    // Loop over all sky pixels and sum solid angle of all pixels that fall
+    // within the sky region
+    for (int i = 0; i < m_num_pixels; ++i) {
+        if (region.contains(inx2dir(i))) {
+            solidangle += this->solidangle(i);
+        }
+    }
+
+    // Return solidangle
+    return solidangle;
+}
+
+
+/***********************************************************************//**
+ * @brief Returns flux within sky regions
+ *
+ * @param[in] region Sky region.
+ * @return Solid angle within sky region.
+ *
+ * The method returns the total solid angle of all sky map pixels that fall
+ * within the specified region. Containment is tested using the
+ * GSkyRegions::contains() method.
+ ***************************************************************************/
+double GSkyMap::solidangle(const GSkyRegions& regions) const
+{
+    // Initialise flux
+    double solidangle = 0.0;
+
+    // Loop over all sky pixels and sum solid angle of all pixels that fall
+    // within the sky regions
+    for (int i = 0; i < m_num_pixels; ++i) {
+        if (regions.contains(inx2dir(i))) {
+            solidangle += this->solidangle(i);
+        }
+    }
+
+    // Return solidangle
+    return solidangle;
+}
+
+
+/***********************************************************************//**
  * @brief Set sky projection
  *
  * @param[in] proj Sky projection.
@@ -2266,7 +2382,7 @@ GSkyRegionCircle GSkyMap::region_circle(void) const
         // Consider edges of pixel with maximum distance
         if (max_index != -1) {
-        
+
             // Set upper pixel edge just beyond the upper boundary since the
             // upper boundary is excluded
             const double up = 0.4999999;
 /***************************************************************************
  *                  test_GSky.cpp - Test sky module                        *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2020 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2021 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -1284,6 +1284,24 @@ void TestGSky::test_GSkyMap(void)
     GNdarray flux_spectrum = map_src.flux();
 	test_value(flux_spectrum.size(), map_src.nmaps(), "Test flux() method");
+    // Test flux in region methods
+    GSkyRegionCircle    gc_circle(266.40499, -28.93617, 3.0);
+    GSkyRegionRectangle gc_box(266.40499, -28.93617, 4.0, 2.0, -58.59867);
+    GSkyRegionMap       gc_map(&gc_box);
+    GSkyRegions         gc_regions;
+    gc_regions.append(gc_circle);
+    gc_regions.append(gc_box);
+    test_value(map_src.flux(gc_circle),  0.49208055, "Test flux(GSkyRegionCircle) method");
+    test_value(map_src.flux(gc_box),     0.12306231, "Test flux(GSkyRegionRectangle) method");
+    test_value(map_src.flux(gc_map),     0.12306231, "Test flux(GSkyRegionMap) method");
+    test_value(map_src.flux(gc_regions), 0.49208055, "Test flux(GSkyRegions) method");
+
+    // Test solidangle in region methods
+    test_value(map_src.solidangle(gc_circle),  0.0097441693, "Test solidangle(GSkyRegionCircle) method");
+    test_value(map_src.solidangle(gc_box),     0.0024368775, "Test solidangle(GSkyRegionRectangle) method");
+    test_value(map_src.solidangle(gc_map),     0.0024368775, "Test solidangle(GSkyRegionMap) method");
+    test_value(map_src.solidangle(gc_regions), 0.0097441693, "Test solidangle(GSkyRegions) method");
+
     // Test map number changing
     GSkyMap map_more = map_src;
     map_more.nmaps(4);