From 7a14eaae3c21796789d28a18980ea612e9528dde Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Mon, 8 Mar 2021 12:46:03 +0100 Subject: [PATCH] Add GSkyMap::flux(GRegion(s)) and GSkyMap::solidangle(GRegion(s)) methods (#1785) --- ChangeLog | 5 +++-- NEWS | 13 ++++++++++++- README.md | 2 +- doc/source/admin/release_history/2.0.rst | 4 +++- include/GSkyMap.hpp | 6 +++++- pyext/GSkyMap.i | 6 +++++- src/sky/GSkyMap.cpp | 130 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------- test/test_GSky.cpp | 20 +++++++++++++++++++- 8 files changed, 171 insertions(+), 15 deletions(-) diff --git a/ChangeLog b/ChangeLog index f0061fe..38d3880 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,4 +1,4 @@ -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) diff --git a/NEWS b/NEWS index 134e696..48c6f8a 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ 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 ------------------ diff --git a/README.md b/README.md index 27576a5..b0f3c49 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ 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/) diff --git a/doc/source/admin/release_history/2.0.rst b/doc/source/admin/release_history/2.0.rst index 555409a..e5cbe3d 100644 --- a/doc/source/admin/release_history/2.0.rst +++ b/doc/source/admin/release_history/2.0.rst @@ -1,6 +1,6 @@ .. _2.0: -GammaLib 2.0 release (3 March 2021) +GammaLib 2.0 release (8 March 2021) =================================== Overview @@ -56,6 +56,8 @@ Improvements * [`3203 `_] - Implement vectorised response computation and analytical gradients for the radial Gaussian model in a stacked analysis +* [`1785 `_] - + Add ``GSkyMap::flux(GRegion)``, ``GSkyMap::flux(GRegions)``, ``GSkyMap::solidangle(GRegion)`` and ``GSkyMap::solidangle(GRegions)`` methods * [`1300 `_] - Implement pre-computation cache in GModelSpatialGauss diff --git a/include/GSkyMap.hpp b/include/GSkyMap.hpp index 3f2ad1b..c5dae36 100644 --- a/include/GSkyMap.hpp +++ b/include/GSkyMap.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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; diff --git a/pyext/GSkyMap.i b/pyext/GSkyMap.i index 63957f4..2585389 100644 --- a/pyext/GSkyMap.i +++ b/pyext/GSkyMap.i @@ -1,7 +1,7 @@ /*************************************************************************** * 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; diff --git a/src/sky/GSkyMap.cpp b/src/sky/GSkyMap.cpp index 8e165ae..1344a46 100644 --- a/src/sky/GSkyMap.cpp +++ b/src/sky/GSkyMap.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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& GSkyMap::operator+=(const GSkyMap& map) } // endfor: looped over all layers } // endfor: looped over all pixels - + } // ... otherwise use interpolation method and add map @@ -507,7 +507,7 @@ GSkyMap& GSkyMap::operator-=(const GSkyMap& map) } // endfor: looped over all layers } // endfor: looped over all pixels - + } // ... otherwise use interpolation method and subtract map @@ -597,7 +597,7 @@ GSkyMap& GSkyMap::operator*=(const GSkyMap& map) } // endfor: looped over all layers } // endfor: looped over all pixels - + } // ... otherwise use interpolation method and multiply map @@ -701,7 +701,7 @@ GSkyMap& GSkyMap::operator/=(const GSkyMap& 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& pixel, const int& 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& pixel, const int& 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& 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; diff --git a/test/test_GSky.cpp b/test/test_GSky.cpp index fdde797..33b5e6d 100644 --- a/test/test_GSky.cpp +++ b/test/test_GSky.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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); -- libgit2 0.21.4