Replace log flag by string in GEbounds class (#2881)

Jürgen Knödlseder
1 parent 2a198690
Showing 10 changed files with 357 additions and 145 deletions
ChangeLog
NEWS
doc/source/admin/release_history/1.7.rst
include/GEbounds.hpp
inst/cta/src/GCTAModelAeffBackground.cpp
pyext/GEbounds.i
src/obs/GEbounds.cpp
test/test_GObservation.cpp
test/test_GObservation.py
test/test_GXspec.cpp
@@ -27,6 +27,7 @@
         Fix GTime.utc() precision bug (#2917)
         Add GTime.datetime() method and add optional precision to GTime.utc() (#2909)
         Remove obsolete methods from GCTAResponseCache class (#2893)
+        Replace log flag by string in GEbounds class (#2881)
         Allow RoI centres different from pointing direction in unbinned analysis
           (#2695)
         Add GSkyDir::moon() method (#1855)
@@ -46,7 +46,8 @@ The following methods have been added:
 - GSkyMap::extract() - Extract section of sky map
 - GSkyMap::region_circle() - Get bounding region circle of sky map
 - GSkyRegionCircle::operator==() - Equality operator
-- GSkyRegionCircle::operator!=() - Non-equality operator 
+- GSkyRegionCircle::operator!=() - Non-equality operator
+- GEbounds::set() - Set energy boundaries for different methods 
 - GModelPar::has_factor_min() - Test if parameter has factor minimum
 - GModelPar::has_factor_max() - Test if parameter has factor maximum
 - GModelPar::remove_factor_min - Remove factor minimum
@@ -64,6 +65,8 @@ The following methods have been added:
 - GCOMDri::cone_content() - Get DRI content within event cone 
  
 The following methods have been removed:
+- GEbounds::set_lin() - method is now private
+- GEbounds::set_log() - method is now private
 - GModelSpatialDiffuseCube::maps()
 - GModelSpatialDiffuseCube::pixels()
  
@@ -78,6 +81,7 @@ The arguments for the following methods have been changed:
 - GTime::utc(int&) - Add optional precision value
 - GSkyDir::posang(GSkyDir&, std::string&) - Add optional coordinate system string
 - GSkyDir::posang_dir(GSkyDir&, std::string&) - Add optional coordinate system string
+- GEbounds::GEbounds(int&, GEnergy&, GEnergy&, std::string&, double&) - Replace bool by string and double
  
 The return value of the following methods has been changed:
 - none
@@ -178,6 +182,14 @@ Add GTime.datetime() method to the Python interface to convert times from and
 to datetime.datetime objects. Add an optional precision argument to the GTime.utc()
 to set the number of digits after the comma for the seconds (#2909).
  
+The interface of the GEbounds constructor that sets a given number of energy
+boundaries was changed. The bool argument that toggled between linear and
+logarithmic spacing was replaced by a string that can take the values "LIN",
+"LOG" or "POW". For the "POW" method an additional optional argument was added
+to specify the power-law index. The GEbounds::set() method was added with the
+same method arguments, and the former set_lin() and set_log() methods were
+made private (#2881).
+
 Add GTime.julian_epoch() methods that return the Julian epoch for a given time
 (#1854).
  
@@ -55,6 +55,8 @@ Improvements
 * Add ``GGti::overlap()`` method.
 * [`2909 <https://cta-redmine.irap.omp.eu/issues/2909>`_] -
   Add ``GTime.datetime()`` method and add optional precision to ``GTime.utc()`` method.
+* [`2881 <https://cta-redmine.irap.omp.eu/issues/2881>`_] -
+  Replace ``log`` flag by ``string`` in ``GEbounds`` class.
 * [`2695 <https://cta-redmine.irap.omp.eu/issues/2695>`_] -
   Allow RoI centres different from pointing direction in unbinned analysis.
 * [`1855 <https://cta-redmine.irap.omp.eu/issues/1855>`_] -
 /***************************************************************************
  *                GEbounds.hpp - Energy boundaries class                   *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2009-2017 by Juergen Knoedlseder                         *
+ *  copyright (C) 2009-2020 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -71,8 +71,11 @@ public:
     explicit GEbounds(const GEnergies& energies);
     GEbounds(const GEbounds& ebds);
     GEbounds(const GEnergy& emin, const GEnergy& emax);
-    GEbounds(const int& num, const GEnergy& emin, const GEnergy& emax,
-             const bool& log = true);
+    GEbounds(const int&         num,
+             const GEnergy&     emin,
+             const GEnergy&     emax,
+             const std::string& method = "LOG",
+             const double&      gamma  = 1.0);
     virtual ~GEbounds(void);
  
     // Operators
@@ -92,8 +95,11 @@ public:
     void           reserve(const int& num);
     void           extend(const GEbounds& ebds);
     void           set(const GEnergies& energies);
-    void           set_lin(const int& num, const GEnergy& emin, const GEnergy& emax);
-    void           set_log(const int& num, const GEnergy& emin, const GEnergy& emax);
+    void           set(const int&         num,
+                       const GEnergy&     emin,
+                       const GEnergy&     emax,
+                       const std::string& method = "LOG",
+                       const double&      gamma  = 1.0);
     void           load(const GFilename& filename);
     void           save(const GFilename& filename,
                         const bool& clobber = false,
@@ -118,7 +124,6 @@ public:
     bool           contains(const GEnergy& emin, const GEnergy& emax) const;
     std::string    print(const GChatter& chatter = NORMAL) const;
  
-
 protected:
     // Protected methods
     void init_members(void);
@@ -126,6 +131,11 @@ protected:
     void free_members(void);
     void set_attributes(void);
     void insert_eng(const int& index, const GEnergy& emin, const GEnergy& emax);
+    void set_lin(const int& num, const GEnergy& emin, const GEnergy& emax);
+    void set_log(const int& num, const GEnergy& emin, const GEnergy& emax);
+    void set_pow(const int& num, const GEnergy& emin, const GEnergy& emax,
+                 const double& gamma);
+
  
     // Protected data area
     int      m_num;         //!< Number of energy boundaries
@@ -563,8 +563,10 @@ GCTAEventList* GCTAModelAeffBackground::mc(const GObservation&amp; obs,
         // a fixed energy sampling for an instance of GModelSpectralNodes.
         // This is analogous to to the GCTAModelIrfBackground::mc method.
         // We make sure that only non-negative nodes get appended.
-        GEbounds spectral_ebounds =
-            GEbounds(m_n_mc_energies, ebounds.emin(), ebounds.emax(), true);
+        GEbounds spectral_ebounds(m_n_mc_energies,
+                                  ebounds.emin(),
+                                  ebounds.emax(),
+                                  "LOG");
         GModelSpectralNodes spectral;
         for (int i = 0; i < spectral_ebounds.size(); ++i) {
             GEnergy energy    = spectral_ebounds.elogmean(i);
 /***************************************************************************
  *                 GEbounds.i - Energy boundaries class                    *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2018 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2020 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -49,8 +49,11 @@ public:
     explicit GEbounds(const GEnergies& energies);
     GEbounds(const GEbounds& ebds);
     GEbounds(const GEnergy& emin, const GEnergy& emax);
-    GEbounds(const int& num, const GEnergy& emin, const GEnergy& emax,
-             const bool& log = true);
+    GEbounds(const int&         num,
+             const GEnergy&     emin,
+             const GEnergy&     emax,
+             const std::string& method = "LOG",
+             const double&      gamma  = 1.0);
     virtual ~GEbounds(void);
  
     // Methods
@@ -67,8 +70,11 @@ public:
     void           reserve(const int& num);
     void           extend(const GEbounds& ebds);
     void           set(const GEnergies& energies);
-    void           set_lin(const int& num, const GEnergy& emin, const GEnergy& emax);
-    void           set_log(const int& num, const GEnergy& emin, const GEnergy& emax);
+    void           set(const int&         num,
+                       const GEnergy&     emin,
+                       const GEnergy&     emax,
+                       const std::string& method = "LOG",
+                       const double&      gamma  = 1.0);
     void           load(const GFilename& filename);
     void           save(const GFilename& filename,
                         const bool& clobber = false,
 /***************************************************************************
  *                  GEbounds.cpp - Energy boundary class                   *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2009-2017 by Juergen Knoedlseder                         *
+ *  copyright (C) 2009-2020 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -41,8 +41,8 @@
 #include "GXmlElement.hpp"
  
 /* __ Method name definitions ____________________________________________ */
-#define G_SET_LIN               "GEbounds::set_lin(int&, GEnergy&, GEnergy&)"
-#define G_SET_LOG               "GEbounds::set_log(int&, GEnergy&, GEnergy&)"
+#define G_SET        "GEbounds::set(int&, GEnergy&, GEnergy&, std::string&, "\
+                                                                   "double&)"
 #define G_READ_XML                             "GEbounds::read(GXmlElement&)"
 #define G_WRITE_XML                           "GEbounds::write(GXmlElement&)"
 #define G_REMOVE                                     "GEbounds::remove(int&)"
@@ -54,6 +54,9 @@
 #define G_ELOGMEAN                                 "GEbounds::elogmean(int&)"
 #define G_EWIDTH                                     "GEbounds::ewidth(int&)"
 #define G_INSERT_ENG         "GEbounds::insert_eng(int&, GEnergy&, GEnergy&)"
+#define G_SET_LIN               "GEbounds::set_lin(int&, GEnergy&, GEnergy&)"
+#define G_SET_LOG               "GEbounds::set_log(int&, GEnergy&, GEnergy&)"
+#define G_SET_POW      "GEbounds::set_pow(int&, GEnergy&, GEnergy&, double&)"
  
 /* __ Macros _____________________________________________________________ */
  
@@ -186,25 +189,25 @@ GEbounds::GEbounds(const GEnergy&amp; emin, const GEnergy&amp; emax)
  * @param[in] num Number of energy intervals.
  * @param[in] emin Minimum energy of first interval.
  * @param[in] emax Maximum energy of last interval.
- * @param[in] log Use logarithmic spacing? (defaults to true).
+ * @param[in] method Energy spacing method (one of "LIN", "LOG" or "POW").
+ * @param[in] gamma Power law index for @p POW method.
  *
  * Constructs energy boundaries by defining @p num successive energy
- * intervals between @p emin and @p emax. The @p log parameter controls
- * whether the energy spacing is logarihmic (default) or linear.
+ * intervals between @p emin and @p emax. The @p method parameter controls
+ * the energy spacing of the energy boundaries. See the set() method for
+ * more information.
  ***************************************************************************/
-GEbounds::GEbounds(const int& num, const GEnergy& emin, const GEnergy& emax,
-                   const bool& log)
+GEbounds::GEbounds(const int&         num,
+                   const GEnergy&     emin,
+                   const GEnergy&     emax,
+                   const std::string& method,
+                   const double&      gamma)
 {
     // Initialise members
     init_members();
  
     // Set intervals
-    if (log) {
-        this->set_log(num, emin, emax);
-    }
-    else {
-        this->set_lin(num, emin, emax);
-    }
+    set(num, emin, emax, method, gamma);
  
     // Return
     return;
@@ -572,119 +575,44 @@ void GEbounds::set(const GEnergies&amp; energies)
  * @param[in] num Number of energy intervals.
  * @param[in] emin Minimum energy of first interval.
  * @param[in] emax Maximum energy of last interval.
+ * @param[in] method Energy spacing method (one of "LIN", "LOG" or "POW").
+ * @param[in] gamma Power law index for @p POW method.
  *
  * @exception GException::invalid_argument
  *            Minimum energy is larger than the maximum energy.
  *
- * Creates @p num linearly spaced energy boundaries running from @p emin to
- * @p emax. If @p num is not a positive integer, no energy
- * boundaries will be set.
- ***************************************************************************/
-void GEbounds::set_lin(const int& num, const GEnergy& emin, const GEnergy& emax)
-{
-    // Initialise members
-    clear();
-
-    // Continue only if the requested number of energy boundaries is positive
-    if (num > 0) {
-
-        // Throw an exception if the maximum energy is smaller than the
-        // minimum energy
-        if (emin > emax) {
-            std::string msg = "Minimum energy "+emin.print()+" is larger than "
-                              "maximum energy "+emax.print()+". Please provide "
-                              "a minimum energy that is not larger than the "
-                              "maximum energy.";
-            throw GException::invalid_argument(G_SET_LIN, msg);
-        }
-
-        // Compute bin width
-        GEnergy ebin = (emax - emin)/double(num);
-
-        // Append boundaries
-        GEnergy min = emin;
-        GEnergy max = emin + ebin;
-        for (int i = 0; i < num; ++i) {
-            append(min, max);
-            min += ebin;
-            max += ebin;
-        }
-
-        // Set attributes
-        set_attributes();
-
-    } // endif: number of energy boundaries was positive
-
-    // Return
-    return;
-}
-
-
-/***********************************************************************//**
- * @brief Set logarithmically spaced energy intervals
- *
- * @param[in] num Number of energy intervals.
- * @param[in] emin Minimum energy of first interval.
- * @param[in] emax Maximum energy of last interval.
+ * Sets energy boundaries by defining @p num successive energy intervals
+ * between @p emin and @p emax. The @p method parameter controls the energy
+ * spacing of the energy boundaries. See the set_lin(), set_log() and
+ * set_pow() methods for more information.
  *
- * @exception GException::invalid_argument
- *            Minimum or maximum energy are not positive or minimum energy
- *            is larger than the maximum energy.
- *
- * Creates @p num logarithmically spaced energy boundaries running from
- * @p emin to @p emax. If @p num is not a positive integer, no energy
- * boundaries will be set.
+ * If @p num is not a positive integer, no energy boundaries will be set.
  ***************************************************************************/
-void GEbounds::set_log(const int& num, const GEnergy& emin, const GEnergy& emax)
+void GEbounds::set(const int&         num,
+                   const GEnergy&     emin,
+                   const GEnergy&     emax,
+                   const std::string& method,
+                   const double&      gamma)
 {
-    // Initialise members
-    clear();
+    // Convert method to upper-case string
+    std::string umethod = gammalib::toupper(method);
  
-    // Continue only if the requested number of energy boundaries is positive
-    if (num > 0) {
-
-        // Throw an exception if the minimum or maximum energy is not positive
-        if (emin.MeV() <= 0.0) {
-            std::string msg = "Non-positive minimum energy "+emin.print()+
-                              " specified. Please provide a positive minimum "
-                              "energy value.";
-            throw GException::invalid_argument(G_SET_LOG, msg);
-        }
-        if (emax.MeV() <= 0.0) {
-            std::string msg = "Non-positive maximum energy "+emax.print()+
-                              " specified. Please provide a positive minimum "
-                              "energy value.";
-            throw GException::invalid_argument(G_SET_LOG, msg);
-        }
-
-        // Throw an exception if the maximum energy is smaller than the
-        // minimum energy
-        if (emin > emax) {
-            std::string msg = "Minimum energy "+emin.print()+" is larger than "
-                              "maximum energy "+emax.print()+". Please provide "
-                              "a minimum energy that is not larger than the "
-                              "maximum energy.";
-            throw GException::invalid_argument(G_SET_LOG, msg);
-        }
-
-        // Compute bin width
-        double elogmin = std::log10(emin.MeV());
-        double elogmax = std::log10(emax.MeV());
-        double elogbin = (elogmax - elogmin)/double(num);
-
-        // Append boundaries
-        GEnergy min;
-        GEnergy max;
-        for (int i = 0; i < num; ++i) {
-            min.MeV(std::pow(10.0, double(i)*elogbin   + elogmin));
-            max.MeV(std::pow(10.0, double(i+1)*elogbin + elogmin));
-            append(min, max);
-        }
-
-        // Set attributes
-        set_attributes();
-
-    } // endif: number of energy boundaries was positive
+    // Dispatch to corresponding set methos
+    if (umethod == "LIN") {
+        set_lin(num, emin, emax);
+    }
+    else if (umethod == "LOG") {
+        set_log(num, emin, emax);
+    }
+    else if (umethod == "POW") {
+        set_pow(num, emin, emax, gamma);
+    }
+    else {
+        std::string msg = "Invalid energy spacing method \""+umethod+"\" "
+                          "specified. Please provide one of \"LIN\", \"LOG\""
+                          ", or \"POW\".";
+        throw GException::invalid_argument(G_SET, msg);
+    }
  
     // Return
     return;
@@ -1458,6 +1386,245 @@ void GEbounds::insert_eng(const int&amp;     index,
 }
  
  
+/***********************************************************************//**
+ * @brief Set linearly spaced energy intervals
+ *
+ * @param[in] num Number of energy intervals.
+ * @param[in] emin Minimum energy of first interval.
+ * @param[in] emax Maximum energy of last interval.
+ *
+ * @exception GException::invalid_argument
+ *            Minimum energy is larger than the maximum energy.
+ *
+ * Creates \f$N\f$ linearly spaced energy boundaries running from
+ * \f$E_{\rm min}\f$ to \f$E_{\rm max}\f$ defined by
+ *
+ * \f[
+ *    E_{\rm min}^{(i)} = E_{\rm min} + \frac{E_{\rm max} - E_{\rm min}}{N}
+ *                        \times i
+ * \f]
+ *
+ * and
+ *
+ * \f[
+ *    E_{\rm max}^{(i)} = E_{\rm min} + \frac{E_{\rm max} - E_{\rm min}}{N}
+ *                        \times (i + 1)
+ * \f]
+ *
+ * If @p num is not positive the method will do nothing.
+ ***************************************************************************/
+void GEbounds::set_lin(const int& num, const GEnergy& emin, const GEnergy& emax)
+{
+    // Continue only if the requested number of energy boundaries is positive
+    if (num > 0) {
+
+        // Throw an exception if the maximum energy is smaller than the
+        // minimum energy
+        if (emin > emax) {
+            std::string msg = "Minimum energy "+emin.print()+" is larger than "
+                              "maximum energy "+emax.print()+". Please provide "
+                              "a minimum energy that is not larger than the "
+                              "maximum energy.";
+            throw GException::invalid_argument(G_SET_LIN, msg);
+        }
+
+        // Compute bin width
+        GEnergy ebin = (emax - emin)/double(num);
+
+        // Append boundaries
+        GEnergy min = emin;
+        GEnergy max = emin + ebin;
+        for (int i = 0; i < num; ++i) {
+            append(min, max);
+            min += ebin;
+            max += ebin;
+        }
+
+        // Set attributes
+        set_attributes();
+
+    } // endif: number of energy boundaries was positive
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Set logarithmically spaced energy intervals
+ *
+ * @param[in] num Number of energy intervals.
+ * @param[in] emin Minimum energy of first interval.
+ * @param[in] emax Maximum energy of last interval.
+ *
+ * @exception GException::invalid_argument
+ *            Minimum or maximum energy are not positive or minimum energy
+ *            is larger than the maximum energy.
+ *
+ * Creates \f$N\f$ logarithmically spaced energy boundaries running from
+ * \f$E_{\rm min}\f$ to \f$E_{\rm max}\f$ defined by
+ *
+ * \f[
+ *    E_{\rm min}^{(i)} = 10^{\log_{10} E_{\rm min} +
+ *                            \frac{\log_{10} E_{\rm max} -
+ *                                  \log_{10} E_{\rm min}}{N}
+ *                            \times i}
+ * \f]
+ *
+ * and
+ *
+ * \f[
+ *    E_{\rm max}^{(i)} = 10^{\log_{10} E_{\rm min} +
+ *                            \frac{\log_{10} E_{\rm max} -
+ *                                  \log_{10} E_{\rm min}}{N}
+ *                            \times (i + 1)}
+ * \f]
+ *
+ * If @p num is not positive the method will do nothing.
+ ***************************************************************************/
+void GEbounds::set_log(const int& num, const GEnergy& emin, const GEnergy& emax)
+{
+    // Initialise members
+    clear();
+
+    // Continue only if the requested number of energy boundaries is positive
+    if (num > 0) {
+
+        // Throw an exception if the minimum or maximum energy is not positive
+        if (emin.MeV() <= 0.0) {
+            std::string msg = "Non-positive minimum energy "+emin.print()+
+                              " specified. Please provide a positive minimum "
+                              "energy value.";
+            throw GException::invalid_argument(G_SET_LOG, msg);
+        }
+        if (emax.MeV() <= 0.0) {
+            std::string msg = "Non-positive maximum energy "+emax.print()+
+                              " specified. Please provide a positive minimum "
+                              "energy value.";
+            throw GException::invalid_argument(G_SET_LOG, msg);
+        }
+
+        // Throw an exception if the maximum energy is smaller than the
+        // minimum energy
+        if (emin > emax) {
+            std::string msg = "Minimum energy "+emin.print()+" is larger than "
+                              "maximum energy "+emax.print()+". Please provide "
+                              "a minimum energy that is not larger than the "
+                              "maximum energy.";
+            throw GException::invalid_argument(G_SET_LOG, msg);
+        }
+
+        // Compute bin width
+        double elogmin = std::log10(emin.MeV());
+        double elogmax = std::log10(emax.MeV());
+        double elogbin = (elogmax - elogmin)/double(num);
+
+        // Append boundaries
+        GEnergy min;
+        GEnergy max;
+        for (int i = 0; i < num; ++i) {
+            min.MeV(std::pow(10.0, double(i)*elogbin   + elogmin));
+            max.MeV(std::pow(10.0, double(i+1)*elogbin + elogmin));
+            append(min, max);
+        }
+
+        // Set attributes
+        set_attributes();
+
+    } // endif: number of energy boundaries was positive
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Set power-law spaced energy intervals
+ *
+ * @param[in] num Number of energy intervals.
+ * @param[in] emin Minimum energy of first interval.
+ * @param[in] emax Maximum energy of last interval.
+ * @param[in] gamma Power law index.
+ *
+ * @exception GException::invalid_argument
+ *            Minimum or maximum energy are not positive or minimum energy
+ *            is larger than the maximum energy.
+ *
+ * Creates \f$N\f$ power-law spaced energy boundaries running from
+ * \f$E_{\rm min}\f$ to \f$E_{\rm max}\f$.
+ *
+ * If @p num is not positive the method will do nothing.
+ ***************************************************************************/
+void GEbounds::set_pow(const int&     num,
+                       const GEnergy& emin,
+                       const GEnergy& emax,
+                       const double&  gamma)
+{
+    // Initialise members
+    clear();
+
+    // Continue only if the requested number of energy boundaries is positive
+    if (num > 0) {
+
+        // Throw an exception if the minimum or maximum energy is not positive
+        if (emin.MeV() <= 0.0) {
+            std::string msg = "Non-positive minimum energy "+emin.print()+
+                              " specified. Please provide a positive minimum "
+                              "energy value.";
+            throw GException::invalid_argument(G_SET_POW, msg);
+        }
+        if (emax.MeV() <= 0.0) {
+            std::string msg = "Non-positive maximum energy "+emax.print()+
+                              " specified. Please provide a positive minimum "
+                              "energy value.";
+            throw GException::invalid_argument(G_SET_POW, msg);
+        }
+
+        // Throw an exception if the maximum energy is smaller than the
+        // minimum energy
+        if (emin > emax) {
+            std::string msg = "Minimum energy "+emin.print()+" is larger than "
+                              "maximum energy "+emax.print()+". Please provide "
+                              "a minimum energy that is not larger than the "
+                              "maximum energy.";
+            throw GException::invalid_argument(G_SET_POW, msg);
+        }
+
+        // Precomputation
+        double a = 1.0 - gamma;
+        double c = (a == 0.0) ? 1.0 / (std::log(emax.MeV())   - std::log(emin.MeV()))
+                              : a   / (std::pow(emax.MeV(),a) - std::pow(emin.MeV(),a));
+        //double b = double(c*(num-1.0));
+        double b = double(c*num);
+
+        // Append boundaries
+        double e = emin.MeV();
+        for (int i = 0; i < num; ++i) {
+
+            // Compute upper bin boundary
+            double log_e_next = (a == 0.0)
+                                ? 1.0/b + std::log(e)
+                                : std::log(a/b + std::pow(e,a)) / a;
+            double e_next = std::exp(log_e_next);
+
+            // Append energy bin
+            append(GEnergy(e, "MeV"), GEnergy(e_next, "MeV"));
+
+            // Store upper boundary as next lower boundary
+            e = e_next;
+
+        }
+
+        // Set attributes
+        set_attributes();
+
+    } // endif: number of energy boundaries was positive
+
+    // Return
+    return;
+}
+
+
 /*==========================================================================
  =                                                                         =
  =                                Friends                                  =
@@ -201,7 +201,7 @@ void TestGObservation::test_ebounds(void)
  
     // Check linear boundaries
     ebds.clear();
-    ebds.set_lin(2, GEnergy(1.0, "MeV"), GEnergy(3.0, "MeV"));
+    ebds.set(2, GEnergy(1.0, "MeV"), GEnergy(3.0, "MeV"), "LIN");
     test_value(ebds.size(), 2, "GEbounds should have 2 elements.");
     test_assert(!ebds.is_empty(), "GEbounds should not be empty.");
     test_value(ebds.emin(0).MeV(), 1.0, 1.0e-10, "Bin 0 minimum energy should be 1.");
@@ -213,7 +213,19 @@ void TestGObservation::test_ebounds(void)
  
     // Check logarithmic boundaries
     ebds.clear();
-    ebds.set_log(2, GEnergy(1.0, "MeV"), GEnergy(100.0, "MeV"));
+    ebds.set(2, GEnergy(1.0, "MeV"), GEnergy(100.0, "MeV"), "LOG");
+    test_value(ebds.size(), 2, "GEbounds should have 2 elements.");
+    test_assert(!ebds.is_empty(), "GEbounds should not be empty.");
+    test_value(ebds.emin(0).MeV(), 1.0, 1.0e-10, "Bin 0 minimum energy should be 1.");
+    test_value(ebds.emin(1).MeV(), 10.0, 1.0e-10, "Bin 1 minimum energy should be 10.");
+    test_value(ebds.emax(0).MeV(), 10.0, 1.0e-10, "Bin 0 maximum energy should be 10.");
+    test_value(ebds.emax(1).MeV(), 100.0, 1.0e-10, "Bin 1 maximum energy should be 100.");
+    test_value(ebds.emin().MeV(), 1.0, 1.0e-10, "Minimum energy should be 1.");
+    test_value(ebds.emax().MeV(), 100.0, 1.0e-10, "Maximum energy should be 100.");
+
+    // Check power-law boundaries
+    ebds.clear();
+    ebds.set(2, GEnergy(1.0, "MeV"), GEnergy(100.0, "MeV"), "POW", 1.0);
     test_value(ebds.size(), 2, "GEbounds should have 2 elements.");
     test_assert(!ebds.is_empty(), "GEbounds should not be empty.");
     test_value(ebds.emin(0).MeV(), 1.0, 1.0e-10, "Bin 0 minimum energy should be 1.");
@@ -238,7 +250,7 @@ void TestGObservation::test_ebounds(void)
     test_value(ebds.emax().MeV(), 1000.0, 1.0e-10, "Maximum energy should be 1000.");
  
     // Check emean, elogmean and ewidth methods
-    ebds.set_log(1, GEnergy(1.0, "MeV"), GEnergy(10.0, "MeV"));
+    ebds.set(1, GEnergy(1.0, "MeV"), GEnergy(10.0, "MeV"), "LOG");
     test_value(ebds.emean(0).MeV(), 5.5, 1.0e-10, "Mean energy should be 5.5.");
     test_value(ebds.elogmean(0).MeV(), 3.16227766017, 1.0e-10, "Log mean energy should be 3.16227766017.");
     test_value(ebds.ewidth(0).MeV(), 9.0, 1.0e-10, "Energy width should be 9.0.");
@@ -259,7 +271,7 @@ void TestGObservation::test_ebounds(void)
     // Check constructing of linear invalid energy boundaries
     test_try("Test constructing of linear invalid energy boundaries");
     try {
-        GEbounds bad(10, GEnergy(100.0, "MeV"), GEnergy(10.0, "MeV"), false);
+        GEbounds bad(10, GEnergy(100.0, "MeV"), GEnergy(10.0, "MeV"), "LIN");
         test_try_failure("Constructing of linear invalid energy boundaries "
                          "shall throw an exception.");
     }
@@ -273,7 +285,7 @@ void TestGObservation::test_ebounds(void)
     // Check constructing of logarithmic invalid energy boundaries
     test_try("Test constructing of logarithmic invalid energy boundaries");
     try {
-        GEbounds bad(10, GEnergy(100.0, "MeV"), GEnergy(10.0, "MeV"), true);
+        GEbounds bad(10, GEnergy(100.0, "MeV"), GEnergy(10.0, "MeV"), "LOG");
         test_try_failure("Constructing of logarithmic invalid energy "
                          "boundaries shall throw an exception.");
     }
@@ -287,7 +299,7 @@ void TestGObservation::test_ebounds(void)
     // Check constructing of logarithmic invalid energy boundaries
     test_try("Test constructing of logarithmic invalid energy boundaries");
     try {
-        GEbounds bad(10, GEnergy(0.0, "MeV"), GEnergy(10.0, "MeV"), true);
+        GEbounds bad(10, GEnergy(0.0, "MeV"), GEnergy(10.0, "MeV"), "LOG");
         test_try_failure("Constructing of logarithmic invalid energy "
                          "boundaries shall throw an exception.");
     }
@@ -301,7 +313,7 @@ void TestGObservation::test_ebounds(void)
     // Check constructing of logarithmic invalid energy boundaries
     test_try("Test constructing of logarithmic invalid energy boundaries");
     try {
-        GEbounds bad(10, GEnergy(10.0, "MeV"), GEnergy(0.0, "MeV"), true);
+        GEbounds bad(10, GEnergy(10.0, "MeV"), GEnergy(0.0, "MeV"), "LOG");
         test_try_failure("Constructing of logarithmic invalid energy "
                          "boundaries shall throw an exception.");
     }
@@ -1370,7 +1382,7 @@ void TestGObservation::test_energies(void)
     test_value(energies[2].MeV(), 100.0, 1.0e-10, "Energy 2 should be 100 MeV.");
  
     // Check energy boundary set method
-    GEbounds ebds1(2, GEnergy(1.0, "MeV"), GEnergy(3.0, "MeV"), false);
+    GEbounds ebds1(2, GEnergy(1.0, "MeV"), GEnergy(3.0, "MeV"), "LIN");
     energies.set(ebds1);
     test_value(energies.size(), 3, "GEbounds constructor (3 elements)");
     test_value(energies[0].MeV(), 1.0, 1.0e-10, "GEbounds constructor (3 elements)");
 # ==========================================================================
 # This module performs unit tests for the GammaLib observation module.
 #
-# Copyright (C) 2012-2019 Juergen Knoedlseder
+# Copyright (C) 2012-2020 Juergen Knoedlseder
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -361,7 +361,7 @@ class Test(gammalib.GPythonTestSuite):
         """
         # Setup energies container
         ebounds = gammalib.GEbounds(10, gammalib.GEnergy(1.0,  'MeV'),
-                                        gammalib.GEnergy(11.0, 'MeV'), False)
+                                        gammalib.GEnergy(11.0, 'MeV'), 'LIN')
  
         # Test ebounds[start:end]
         self.test_value(len(ebounds[3:5]), 2)
 /***************************************************************************
  *                  test_GXspec.hpp - Test Xspec module                    *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2013-2017 by Juergen Knoedlseder                         *
+ *  copyright (C) 2013-2020 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -88,7 +88,7 @@ void TestGXspec::test_GPha(void)
     test_value(pha1.ancrfile(), "", "Check if empty PHA has empty ARF");
  
     // Test energy boundary constructor, filling and accessing
-    GEbounds ebds2(9, GEnergy(1.0, "TeV"), GEnergy(10.0, "TeV"), false);
+    GEbounds ebds2(9, GEnergy(1.0, "TeV"), GEnergy(10.0, "TeV"), "LIN");
     GPha     pha2(ebds2);
     pha2.exposure(1234.5);
     pha2.backfile("background");