From 717bee7d2614ddf1c1f04814e3eb37cc8ef87c9e Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Thu, 6 Jul 2017 14:51:03 +0200 Subject: [PATCH] Add margin to maximum light and phase curve values (#2149) --- ChangeLog | 1 + NEWS | 6 +++++- doc/source/admin/release_history/1.3.1.rst | 2 ++ src/model/GModelTemporalLightCurve.cpp | 20 +++++++++++++++----- src/model/GModelTemporalPhaseCurve.cpp | 20 +++++++++++++++----- 5 files changed, 38 insertions(+), 11 deletions(-) diff --git a/ChangeLog b/ChangeLog index 1168e20..a93d0c2 100644 --- a/ChangeLog +++ b/ChangeLog @@ -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) diff --git a/NEWS b/NEWS index fffaf36..9dd3635 100644 --- a/NEWS +++ b/NEWS @@ -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 diff --git a/doc/source/admin/release_history/1.3.1.rst b/doc/source/admin/release_history/1.3.1.rst index 4e20fc1..ac1bdfe 100644 --- a/doc/source/admin/release_history/1.3.1.rst +++ b/doc/source/admin/release_history/1.3.1.rst @@ -10,6 +10,8 @@ GammaLib 1.3.1 is a bug fix for GammaLib release 1.3. The following bugs have been fixed: +* [`2149 `_] - + Add margin to maximum light and phase curve values * [`2133 `_] - Make GCTAPsfKing save against invalid values diff --git a/src/model/GModelTemporalLightCurve.cpp b/src/model/GModelTemporalLightCurve.cpp index 71bcd00..59f016d 100644 --- a/src/model/GModelTemporalLightCurve.cpp +++ b/src/model/GModelTemporalLightCurve.cpp @@ -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& 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& 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& 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); diff --git a/src/model/GModelTemporalPhaseCurve.cpp b/src/model/GModelTemporalPhaseCurve.cpp index 18559bd..1097258 100644 --- a/src/model/GModelTemporalPhaseCurve.cpp +++ b/src/model/GModelTemporalPhaseCurve.cpp @@ -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& 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& 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& 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(); -- libgit2 0.21.4