Add margin to maximum light and phase curve values (#2149)

Jürgen Knödlseder
1 parent e13c97b5
Showing 5 changed files with 38 additions and 11 deletions
ChangeLog
NEWS
doc/source/admin/release_history/1.3.1.rst
src/model/GModelTemporalLightCurve.cpp
src/model/GModelTemporalPhaseCurve.cpp
@@ -3,6 +3,7 @@
         * Bug fix version 1.3.1 released
         ================================
  
+        Add margin to maximum light and phase curve values (#2149)
         Improve caching in GSkyDir (#2148)
         Make GCTAPsfKing save against invalid values (#2133)
  
@@ -63,7 +63,11 @@ None
  
 8. Model module
 ---------------
-None
+The normalisation values of light and phase curves shall not exceed values of 1
+but due to rounding problems when reading the FITS file it can happen that
+slightly larger value occur. To circumvent the issue a small margin of 1e-8 has
+been added to the maximum value of 1, well below the required precision for
+flux normalisations (#2149).
  
  
 9. Numerics module
@@ -10,6 +10,8 @@ GammaLib 1.3.1 is a bug fix for GammaLib release 1.3.
  
 The following bugs have been fixed:
  
+* [`2149 <https://cta-redmine.irap.omp.eu/issues/2149>`_] -
+  Add margin to maximum light and phase curve values
 * [`2133 <https://cta-redmine.irap.omp.eu/issues/2133>`_] -
   Make GCTAPsfKing save against invalid values
  
@@ -587,6 +587,9 @@ void GModelTemporalLightCurve::free_members(void)
  ***************************************************************************/
 void GModelTemporalLightCurve::load_nodes(const GFilename& filename)
 {
+    // Set maximum light curve normalization value, including a small margin
+    const double max_norm = 1.0 + 1.0e-8;
+
     // Clear nodes and values
     m_nodes.clear();
     m_values.clear();
@@ -630,7 +633,7 @@ void GModelTemporalLightCurve::load_nodes(const GFilename&amp; filename)
     int nodes = time_col->nrows();
  
     // Check that time values are in ascending order and that no node is
-    // larger than one
+    // larger than 1
     double last_time = -1.0;
     for (int i = 0; i < nodes; ++i) {
  
@@ -644,12 +647,12 @@ void GModelTemporalLightCurve::load_nodes(const GFilename&amp; filename)
             throw GException::invalid_value(G_LOAD_NODES, msg);
         }
  
-        // Check if value is <= 1
-        if (norm_col->real(i) > 1.0) {
+        // Check if value is smaller than maximum allowed normalisation
+        if (norm_col->real(i) > max_norm) {
             std::string msg = "Value "+gammalib::str(norm_col->real(i))+" at "
                               "time "+gammalib::str(time_col->real(i))+" is "
-                              "above 1. Please provide a light curve file with "
-                              "normalizations not exceeding 1.";
+                              "larger than 1. Please provide a light curve file "
+                              "with normalizations not exceeding 1.";
             throw GException::invalid_value(G_LOAD_NODES, msg);
         }
  
@@ -661,6 +664,13 @@ void GModelTemporalLightCurve::load_nodes(const GFilename&amp; filename)
         m_values.push_back(norm_col->real(i));
     }
  
+    // Make sure that no node exceeds 1
+    for (int i = 0; i < m_values.size(); ++i) {
+        if (m_values[i] > 1.0) {
+            m_values[i] = 1.0;
+        }
+    }
+
     // Set minimum and maximum times (assumes that times are ordered)
     m_tmin.set(time_col->real(0), m_timeref);
     m_tmax.set(time_col->real(time_col->nrows()-1), m_timeref);
@@ -688,6 +688,9 @@ void GModelTemporalPhaseCurve::free_members(void)
  ***************************************************************************/
 void GModelTemporalPhaseCurve::load_nodes(const GFilename& filename)
 {
+    // Set maximum phase curve normalization value, including a small margin
+    const double max_norm = 1.0 + 1.0e-8;
+
     // Clear nodes and values
     m_nodes.clear();
     m_values.clear();
@@ -728,7 +731,7 @@ void GModelTemporalPhaseCurve::load_nodes(const GFilename&amp; filename)
     int nodes = phase_col->nrows();
  
     // Check that phase values are in ascending order and comprised between
-    // 0 and 1, and that no node is larger than one
+    // 0 and 1, and that no node is larger than 1
     double last_phase = -1.0;
     for (int i = 0; i < nodes; ++i) {
  
@@ -751,12 +754,12 @@ void GModelTemporalPhaseCurve::load_nodes(const GFilename&amp; filename)
             throw GException::invalid_value(G_LOAD_NODES, msg);
         }
  
-        // Check if value is <= 1
-        if (norm_col->real(i) > 1.0) {
+        // Check if value is smaller than maximum allowed normalisation
+        if (norm_col->real(i) > max_norm) {
             std::string msg = "Value "+gammalib::str(norm_col->real(i))+" at "
                               "phase "+gammalib::str(phase_col->real(i))+" is "
-                              "above 1. Please provide a phase curve file with "
-                              "normalizations not exceeding 1.";
+                              "larger than 1. Please provide a phase curve file "
+                              "with normalizations not exceeding 1.";
             throw GException::invalid_value(G_LOAD_NODES, msg);
         }
  
@@ -784,6 +787,13 @@ void GModelTemporalPhaseCurve::load_nodes(const GFilename&amp; filename)
         m_values.push_back(norm_col->real(0));
     }
  
+    // Make sure that no node exceeds 1
+    for (int i = 0; i < m_values.size(); ++i) {
+        if (m_values[i] > 1.0) {
+            m_values[i] = 1.0;
+        }
+    }
+
     // Close FITS file
     fits.close();