Renamed PowerLaw2 into PowerLawPhotonFlux and added GModelSpectralPlawEnergyFlux…

… including docs and unit tests

Renamed PowerLaw2 into PowerLawPhotonFlux and added GModelSpectralPlawEnergyFlux…
… including docs and unit tests
Michael Mayer
1 parent 548f06e8
Showing 17 changed files with 1772 additions and 237 deletions
doc/dev/maths/gammalib_maths.tex
doc/source/develop/coding/design.rst
doc/source/users/user_manual/modules/model.rst
include/GModelSpectralPlaw2.hpp → include/GModelSpectralPlawEnergyFlux.hpp
include/GModelSpectralPlawPhotonFlux.hpp
include/GammaLib.hpp
include/Makefile.am
pyext/GModelSpectralPlawEnergyFlux.i
pyext/GModelSpectralPlaw2.i → pyext/GModelSpectralPlawPhotonFlux.i
pyext/gammalib/model.i
src/model/GModelSpectralPlawEnergyFlux.cpp
src/model/GModelSpectralPlaw2.cpp → src/model/GModelSpectralPlawPhotonFlux.cpp
src/model/Makefile.am
test/data/model_point_plaw_eflux.xml
test/data/model_point_plaw2.xml → test/data/model_point_plaw_phflux.xml
test/test_GModel.cpp
test/test_GModel.hpp
@@ -347,9 +347,9 @@ The parameter value gradients for the power law are given by
  
  
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection{PowerLaw2}
+\subsection{PowerLawPhotonFlux}
  
-This flavour of the power law spectral model uses the integral flux $f$ within the
+This flavour of the power law spectral model uses the integral photon flux $f$ within the
 energy range $E_{\rm min}$ and $E_{\rm max}$ as free parameter instead of
 the normalization $k$.
 The integral flux $f$ is given by
@@ -380,7 +380,7 @@ is obtained from Eq.~(\ref{eq:fluxint}).
 Each of the 4 parameters is factorised into a scaling factor and a value, e.g.
 $f=f_s f_v$ and $\gamma=\gamma_s \gamma_v$.
 It is assumed that $E_{\rm min}$ and $E_{\rm max}$ are fixed parameters,
-and {\tt GModelSpectralPlaw2::eval\_gradients} returns valid gradients only for 
+and {\tt GModelSpectralPlawPhotonFlux::eval\_gradients} returns valid gradients only for 
 $f_v$ and $\gamma_v$.
  
 The flux value gradient $f_v$ is given by
@@ -429,4 +429,88 @@ Note that for $\gamma \ne -1$
           {\left( E_{\rm max}^{\gamma+1}-E_{\rm min}^{\gamma+1} \right)}
 \end{equation}
  
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{PowerLawEnergyFlux}
+
+This flavour of the power law spectral model uses the integral energy flux $F$ within the
+energy range $E_{\rm min}$ and $E_{\rm max}$ as free parameter instead of
+the normalization $k$.
+The integral energy flux $F$ is given by
+\begin{equation}
+F = \int_{E_{\rm min}}^{E_{\rm max}} I(E) {\rm d}E
+\label{eq:fluxint}
+\end{equation}
+The power law model is then defined by
+\begin{equation}
+I(E) = \tilde{k} \, E^{\gamma}
+\end{equation}
+where
+\begin{equation}
+   \tilde{k} = \left \{
+   \begin{array}{l l}
+      \displaystyle
+      \frac{F}{\ln E_{\rm max} - \ln E_{\rm min}} 
+        & \mbox{if $\gamma = -2$} \\
+     \\
+     \displaystyle
+      \frac{F (2+\gamma)}{E_{\rm max}^{\gamma+2}-E_{\rm min}^{\gamma+2}} 
+        & \mbox{else}
+   \end{array}
+   \right .
+\end{equation}
+is obtained from Eq.~(\ref{eq:fluxenergy}).
+
+Each of the 4 parameters is factorised into a scaling factor and a value, e.g.
+$F=F_s F_v$ and $\gamma=\gamma_s \gamma_v$.
+It is assumed that $E_{\rm min}$ and $E_{\rm max}$ are fixed parameters,
+and {\tt GModelSpectralPlawEnergyFlux::eval\_gradients} returns valid gradients only for 
+$F_v$ and $\gamma_v$.
+
+The flux value gradient $f_v$ is given by
+\begin{equation}
+\frac{\delta I}{\delta F_v} = 
+  \frac{\delta I}{\delta F} \frac{\delta F}{\delta F_v} = 
+  \frac{\delta I}{\delta F} F_s = 
+  \frac{\delta \tilde{k}}{\delta F} \, E^{\gamma} F_s = 
+  \frac{\tilde{k}}{F} \, E^{\gamma} F_s = 
+  \frac{I(E)}{F} F_s =
+  \frac{I(E)}{F_v}
+\end{equation}
+
+The index value gradient $\gamma_v$ is given by
+\begin{equation}
+\frac{\delta I}{\delta \gamma_v} =
+  \frac{\delta I}{\delta \gamma} \frac{\delta \gamma}{\delta \gamma_v} = 
+  \frac{\delta I}{\delta \gamma} \gamma_s = 
+  \left( \frac{\delta \tilde{k}}{\delta \gamma} \, E^{\gamma} + \tilde{k} \, E^{\gamma} \ln E \right)  \gamma_s =
+  \left( \frac{1}{\tilde{k}} \frac{\delta \tilde{k}}{\delta \gamma} + \ln E \right) \tilde{k} E^{\gamma} \gamma_s =
+  \left( \frac{1}{\tilde{k}} \frac{\delta \tilde{k}}{\delta \gamma} + \ln E \right) I(E) \gamma_s
+\end{equation}
+where
+\begin{equation}
+   \frac{\delta \tilde{k}}{\delta \gamma} = \left \{
+   \begin{array}{l l}
+     \displaystyle
+     0 & \mbox{if $\gamma = -2$} \\
+     \\
+     \\
+     \displaystyle
+      \frac{F \left( E_{\rm max}^{\gamma+2}-E_{\rm min}^{\gamma+2} \right) -
+              F (1+\gamma) \left( E_{\rm max}^{\gamma+2} \ln E_{\rm max} - 
+                                      E_{\rm min}^{\gamma+2} \ln E_{\rm min} \right)}
+              {\left( E_{\rm max}^{\gamma+2}-E_{\rm min}^{\gamma+2} \right)^2}
+         & \mbox{else}
+   \end{array}
+   \right .
+\end{equation}
+Note that for $\gamma \ne -2$
+\begin{equation}
+\frac{1}{\tilde{k}} \frac{\delta \tilde{k}}{\delta \gamma} =
+  \frac{1}{2+\gamma} -
+  \frac{\left( E_{\rm max}^{\gamma+2} \ln E_{\rm max} - 
+                    E_{\rm min}^{\gamma+2} \ln E_{\rm min} \right)}
+          {\left( E_{\rm max}^{\gamma+2}-E_{\rm min}^{\gamma+2} \right)}
+\end{equation}
+
 \end{document} 
@@ -397,11 +397,11 @@ non-accessible class members need then to be declared ``mutable`` to avoid
 compiler errors. Alternatively, the ``const_cast`` declaration can be used
 to allow member modifications within a ``const`` method.
  
-As example we show here part of the definition of :doxy:`GModelSpectralPlaw2`:
+As example we show here part of the definition of :doxy:`GModelSpectralPlawPhotonFlux`:
  
 .. code-block:: cpp
  
-    class GModelSpectralPlaw2 : public GModelSpectral {
+    class GModelSpectralPlawPhotonFlux : public GModelSpectral {
     public:
         virtual double eval(const GEnergy& srcEng) const;
         virtual void   read(const GXmlElement& xml);
@@ -430,11 +430,11 @@ As example we show here part of the definition of :doxy:`GModelSpectralPlaw2`:
         mutable double  m_last_g_index;    //!< Last spectral index gradient
  
 This class has an internal cache for precomputation, which is
-potentially updated when :doxy:`GModelSpectralPlaw2::eval` is called. Here the corresponding code:
+potentially updated when :doxy:`GModelSpectralPlawPhotonFlux::eval` is called. Here the corresponding code:
  
 .. code-block:: cpp
  
-    double GModelSpectralPlaw2::eval(const GEnergy& srcEng) const
+    double GModelSpectralPlawPhotonFlux::eval(const GEnergy& srcEng) const
     {
         // Update precomputed values
         update(srcEng);
@@ -445,17 +445,17 @@ format:
     <parameter name="Scale"     scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="0"/>
    </spectrum>
  
-An alternative power law function is defined by the :doxy:`GModelSpectralPlaw2`
+An alternative power law function is defined by the :doxy:`GModelSpectralPlawPhotonFlux`
 class that uses the integral photon flux as parameter rather than the
 Prefactor:
  
 .. math::
-    \frac{dN}{dE} = \frac{N(\gamma+1)E^{\gamma}}
+    \frac{dN}{dE} = \frac{F_{\rm Ph}(\gamma+1)E^{\gamma}}
                          {E_{\rm max}^{\gamma+1} - E_{\rm min}^{\gamma+1}}
  
 where the parameters in the XML definition have the following mappings:
  
-* :math:`N` = ``PhotonFlux``
+* :math:`F_{\rm Ph}` = ``PhotonFlux``
 * :math:`\gamma` = ``Index``
 * :math:`E_{\rm min}` = ``LowerLimit``
 * :math:`E_{\rm max}` = ``UpperLimit``
@@ -487,12 +487,47 @@ format:
 .. note::
  
    The UpperLimit and LowerLimit parameters are always treated as fixed and,
-   as should be apparent from this definition, the flux given by the Integral
+   as should be apparent from this definition, the flux given by the PhotonFlux
    parameter is over the range (LowerLimit, UpperLimit). Use of this model
-   allows the errors on the integrated flux to be evaluated directly by
+   allows the errors on the integrated photon flux to be evaluated directly by
    likelihood, obviating the need to propagate the errors if one is using
    the PowerLaw form.
  
+Another alternative power law function is defined by the :doxy:`GModelSpectralPlawEnergyFlux`
+class that uses the integral energy flux as parameter rather than the
+Prefactor:
+
+.. math::
+    \frac{dN}{dE} = \frac{F_{\rm E}(\gamma+2)E^{\gamma}}
+                         {E_{\rm max}^{\gamma+2} - E_{\rm min}^{\gamma+2}}
+
+where the parameters in the XML definition have the following mappings:
+
+* :math:`F_{\rm E}` = ``EnergyFlux``
+* :math:`\gamma` = ``Index``
+* :math:`E_{\rm min}` = ``LowerLimit``
+* :math:`E_{\rm max}` = ``UpperLimit``
+
+The XML format for specifying a power law defined by the integral energy flux
+is:
+
+.. code-block:: xml
+
+   <spectrum type="PowerLaw">
+    <parameter scale="1e-07" name="EnergyFlux" min="1e-07" max="1000.0"    value="1.0" free="1"/>
+    <parameter scale="1.0"   name="Index"      min="-5.0"  max="+5.0"      value="-2.0" free="1"/>
+    <parameter scale="1.0"   name="LowerLimit" min="10.0"  max="1000000.0" value="100.0" free="0"/>
+    <parameter scale="1.0"   name="UpperLimit" min="10.0"  max="1000000.0" value="500000.0" free="0"/>
+   </spectrum>
+
+.. note::
+
+   The UpperLimit and LowerLimit parameters are always treated as fixed and,
+   as should be apparent from this definition, the flux given by the EnergyFlux
+   parameter is over the range (LowerLimit, UpperLimit). Use of this model
+   allows the errors on the integrated energy flux to be evaluated directly by
+   likelihood, obviating the need to propagate the errors if one is using
+   the PowerLaw form.
  
 Exponentially cut-off power law
 ===============================
 /***************************************************************************
- *        GModelSpectralPlaw2.hpp - Spectral power law model class         *
+ *    GModelSpectralPlawEnergyFlux.hpp - Spectral power law model class    *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2011-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2016 by Michael Mayer                                    *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -19,13 +19,13 @@
  *                                                                         *
  ***************************************************************************/
 /**
- * @file GModelSpectralPlaw2.hpp
- * @brief Flux normalized power law spectral model class interface definition
- * @author Juergen Knoedlseder
+ * @file GModelSpectralPlawEnergyFlux.hpp
+ * @brief Energy flux normalized power law spectral model class interface definition
+ * @author Michael Mayer
  */
  
-#ifndef GMODELSPECTRALPLAW2_HPP
-#define GMODELSPECTRALPLAW2_HPP
+#ifndef GMODELSPECTRALPLAWENERGYFLUX_HPP
+#define GMODELSPECTRALPLAWENERGYFLUX_HPP
  
 /* __ Includes ___________________________________________________________ */
 #include <string>
@@ -40,60 +40,60 @@ class GXmlElement;
  
  
 /***********************************************************************//**
- * @class GModelSpectralPlaw2
+ * @class GModelSpectralPlawEnergyFlux
  *
- * @brief Flux normalized power law spectral model class
+ * @brief Energy flux normalized power law spectral model class
  *
  * This class implements a power law spectrum. The model is defined by
  *
  * \f[
- *    S_{\rm E}(E | t) = {\tt m\_integral}
+ *    S_{\rm E}(E | t) = {\tt m\_energy\_flux}
  *    \frac{{\tt m\_index}+1}
- *         {{\tt e\_max}^{{\tt m\_index}+1} -
- *          {\tt e\_min}^{{\tt m\_index}+1}}
+ *         {{\tt e\_max}^{{\tt m\_index}+2} -
+ *          {\tt e\_min}^{{\tt m\_index}+2}}
  *    E^{\tt m\_index}
  * \f]
  *
- * for \f${\tt m\_index} \ne -1\f$ and
+ * for \f${\tt m\_index} \ne -2\f$ and
  *
  * \f[
  *    S_{\rm E}(E | t) = 
- *    \frac{{\tt m\_integral}}
+ *    \frac{{\tt m\_energy\_flux}}
  *         {\log {\tt e\_max} - \log {\tt e\_min}}
  *    E^{\tt m\_index}
  * \f]
  *
- * for \f${\tt m\_index} = -1\f$, where
+ * for \f${\tt m\_index} = -2\f$, where
  * - \f${\tt e\_min}\f$ is the minimum energy of an interval,
  * - \f${\tt e\_max}\f$ is the maximum energy of an interval,
- * - \f${\tt m\_integral}\f$ is the integral flux between 
+ * - \f${\tt m\_energy\_flux}\f$ is the energy flux between
  *   \f${\tt e\_min}\f$ and \f${\tt e\_max}\f$, and
  * - \f${\tt m\_index}\f$ is the spectral index.
  ***************************************************************************/
-class GModelSpectralPlaw2 : public GModelSpectral {
+class GModelSpectralPlawEnergyFlux : public GModelSpectral {
  
 public:
     // Constructors and destructors
-    GModelSpectralPlaw2(void);
-    GModelSpectralPlaw2(const std::string& type,
-                        const std::string& integral,
+    GModelSpectralPlawEnergyFlux(void);
+    GModelSpectralPlawEnergyFlux(const std::string& type,
+                        const std::string& energy_flux,
                         const std::string& index,
                         const std::string& emin,
                         const std::string& emax);
-    GModelSpectralPlaw2(const double&  integral,
-                        const double&  index,
-                        const GEnergy& emin,
-                        const GEnergy& emax);
-    explicit GModelSpectralPlaw2(const GXmlElement& xml);
-    GModelSpectralPlaw2(const GModelSpectralPlaw2& model);
-    virtual ~GModelSpectralPlaw2(void);
+    explicit GModelSpectralPlawEnergyFlux(const double&  eflux,
+                                 const double&  index,
+                                 const GEnergy& emin,
+                                 const GEnergy& emax);
+    explicit GModelSpectralPlawEnergyFlux(const GXmlElement& xml);
+    GModelSpectralPlawEnergyFlux(const GModelSpectralPlawEnergyFlux& model);
+    virtual ~GModelSpectralPlawEnergyFlux(void);
  
     // Operators
-    virtual GModelSpectralPlaw2& operator=(const GModelSpectralPlaw2& model);
+    virtual GModelSpectralPlawEnergyFlux& operator=(const GModelSpectralPlawEnergyFlux& model);
  
     // Implemented pure virtual base class methods
     virtual void                 clear(void);
-    virtual GModelSpectralPlaw2* clone(void) const;
+    virtual GModelSpectralPlawEnergyFlux* clone(void) const;
     virtual std::string          classname(void) const;
     virtual std::string          type(void) const;
     virtual double               eval(const GEnergy& srcEng,
@@ -113,8 +113,8 @@ public:
     virtual std::string          print(const GChatter& chatter = NORMAL) const;
  
     // Other methods
-    double  integral(void) const;
-    void    integral(const double& integral);
+    double  energy_flux(void) const;
+    void    energy_flux(const double& energy_flux);
     double  index(void) const;
     void    index(const double& index);
     GEnergy emin(void) const;
@@ -125,45 +125,45 @@ public:
 protected:
     // Protected methods
     void init_members(void);
-    void copy_members(const GModelSpectralPlaw2& model);
+    void copy_members(const GModelSpectralPlawEnergyFlux& model);
     void free_members(void);
     void update(const GEnergy& srcEng) const;
  
     // Protected members
-    std::string     m_type;            //!< Model type
-    GModelPar       m_integral;        //!< Integral flux
-    GModelPar       m_index;           //!< Spectral index
-    GModelPar       m_emin;            //!< Lower energy limit (MeV)
-    GModelPar       m_emax;            //!< Upper energy limit (MeV)
-
-    // Cached members used for pre-computations
-    mutable double  m_log_emin;        //!< Log(emin)
-    mutable double  m_log_emax;        //!< Log(emax)
-    mutable double  m_pow_emin;        //!< emin^(index+1)
-    mutable double  m_pow_emax;        //!< emax^(index+1)
-    mutable double  m_norm;            //!< Power-law normalization (for pivot energy 1 MeV)
-    mutable double  m_g_norm;          //!< Power-law normalization gradient
-    mutable double  m_power;           //!< Power-law factor
-    mutable double  m_last_integral;   //!< Last integral flux
-    mutable double  m_last_index;      //!< Last spectral index (MeV)
-    mutable GEnergy m_last_emin;       //!< Last lower energy limit
-    mutable GEnergy m_last_emax;       //!< Last upper energy limit
-    mutable GEnergy m_last_energy;     //!< Last source energy
-    mutable double  m_last_value;      //!< Last function value
-    mutable double  m_last_g_integral; //!< Last integral flux gradient
-    mutable double  m_last_g_index;    //!< Last spectral index gradient
+	std::string     m_type;            //!< Model type
+	GModelPar       m_energy_flux;     //!< Energy flux
+	GModelPar       m_index;           //!< Spectral index
+	GModelPar       m_emin;            //!< Lower energy limit (MeV)
+	GModelPar       m_emax;            //!< Upper energy limit (MeV)
+
+	// Cached members used for pre-computations
+	mutable double  m_log_emin;        //!< Log(emin)
+	mutable double  m_log_emax;        //!< Log(emax)
+	mutable double  m_pow_emin;        //!< emin^(index+1)
+	mutable double  m_pow_emax;        //!< emax^(index+1)
+	mutable double  m_norm;            //!< Power-law normalization (for pivot energy 1 MeV)
+	mutable double  m_g_norm;          //!< Power-law normalization gradient
+	mutable double  m_power;           //!< Power-law factor
+	mutable double  m_last_phflux;     //!< Last energy flux
+	mutable double  m_last_index;      //!< Last spectral index (MeV)
+	mutable GEnergy m_last_emin;       //!< Last lower energy limit
+	mutable GEnergy m_last_emax;       //!< Last upper energy limit
+	mutable GEnergy m_last_energy;     //!< Last source energy
+	mutable double  m_last_value;      //!< Last function value
+	mutable double  m_last_g_phflux;   //!< Last energy flux gradient
+	mutable double  m_last_g_index;    //!< Last spectral index gradient
 };
  
  
 /***********************************************************************//**
  * @brief Return class name
  *
- * @return String containing the class name ("GModelSpectralPlaw2").
+ * @return String containing the class name ("GModelSpectralPlawEnergyFlux").
  ***************************************************************************/
 inline
-std::string GModelSpectralPlaw2::classname(void) const
+std::string GModelSpectralPlawEnergyFlux::classname(void) const
 {
-    return ("GModelSpectralPlaw2");
+    return ("GModelSpectralPlawEnergyFlux");
 }
  
  
@@ -175,43 +175,41 @@ std::string GModelSpectralPlaw2::classname(void) const
  * Returns the type of the spectral power law model.
  ***************************************************************************/
 inline
-std::string GModelSpectralPlaw2::type(void) const
+std::string GModelSpectralPlawEnergyFlux::type(void) const
 {
     return (m_type);
 }
  
  
 /***********************************************************************//**
- * @brief Return integral flux
+ * @brief Return energy flux
  *
- * @return Integral flux (ph/cm2/s).
+ * @return Energy flux (erg/cm2/s).
  *
- * Returns the integral flux.
+ * Returns the energy flux.
  ***************************************************************************/
 inline
-double GModelSpectralPlaw2::integral(void) const
+double GModelSpectralPlawEnergyFlux::energy_flux(void) const
 {
-    return (m_integral.value());
+    return (m_energy_flux.value());
 }
  
  
 /***********************************************************************//**
- * @brief Set integral flux
+ * @brief Set energy flux
  *
- * @param[in] integral Integral flux (ph/cm2/s).
+ * @param[in] energy_flux Energy flux (erg/cm2/s).
  *
- * Sets the integral flux.
+ * Sets the energy flux.
  ***************************************************************************/
 inline
-void GModelSpectralPlaw2::integral(const double& integral)
+void GModelSpectralPlawEnergyFlux::energy_flux(const double& energy_flux)
 {
-    m_integral.value(integral);
+    m_energy_flux.value(energy_flux);
     return;
 }
  
  
-
-
 /***********************************************************************//**
  * @brief Return power law index
  *
@@ -220,7 +218,7 @@ void GModelSpectralPlaw2::integral(const double&amp; integral)
  * Returns the power law index.
  ***************************************************************************/
 inline
-double GModelSpectralPlaw2::index(void) const
+double GModelSpectralPlawEnergyFlux::index(void) const
 {
     return (m_index.value());
 }
@@ -234,7 +232,7 @@ double GModelSpectralPlaw2::index(void) const
  * Sets the power law index.
  ***************************************************************************/
 inline
-void GModelSpectralPlaw2::index(const double& index)
+void GModelSpectralPlawEnergyFlux::index(const double& index)
 {
     m_index.value(index);
     return;
@@ -249,7 +247,7 @@ void GModelSpectralPlaw2::index(const double&amp; index)
  * Returns the minimum energy.
  ***************************************************************************/
 inline
-GEnergy GModelSpectralPlaw2::emin(void) const
+GEnergy GModelSpectralPlawEnergyFlux::emin(void) const
 {
     GEnergy energy;
     energy.MeV(m_emin.value());
@@ -265,7 +263,7 @@ GEnergy GModelSpectralPlaw2::emin(void) const
  * Sets the minimum energy.
  ***************************************************************************/
 inline
-void GModelSpectralPlaw2::emin(const GEnergy& emin)
+void GModelSpectralPlawEnergyFlux::emin(const GEnergy& emin)
 {
     m_emin.value(emin.MeV());
     return;
@@ -280,7 +278,7 @@ void GModelSpectralPlaw2::emin(const GEnergy&amp; emin)
  * Returns the maximum energy.
  ***************************************************************************/
 inline
-GEnergy GModelSpectralPlaw2::emax(void) const
+GEnergy GModelSpectralPlawEnergyFlux::emax(void) const
 {
     GEnergy energy;
     energy.MeV(m_emax.value());
@@ -296,10 +294,10 @@ GEnergy GModelSpectralPlaw2::emax(void) const
  * Sets the maximum energy.
  ***************************************************************************/
 inline
-void GModelSpectralPlaw2::emax(const GEnergy& emax)
+void GModelSpectralPlawEnergyFlux::emax(const GEnergy& emax)
 {
     m_emax.value(emax.MeV());
     return;
 }
  
-#endif /* GMODELSPECTRALPLAW2_HPP */
+#endif /* GMODELSPECTRALPLAWENERGYFLUX_HPP */
+/***************************************************************************
+ *   GModelSpectralPlawPhotonFlux.hpp - Spectral power law model class     *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2011-2016 by 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   *
+ *  the Free Software Foundation, either version 3 of the License, or      *
+ *  (at your option) any later version.                                    *
+ *                                                                         *
+ *  This program is distributed in the hope that it will be useful,        *
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
+ *  GNU General Public License for more details.                           *
+ *                                                                         *
+ *  You should have received a copy of the GNU General Public License      *
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file GModelSpectralPlawPhotonFlux.hpp
+ * @brief Flux normalized power law spectral model class interface definition
+ * @author Juergen Knoedlseder
+ */
+
+#ifndef GMODELSPECTRALPLAWPHOTONFLUX_HPP
+#define GMODELSPECTRALPLAWPHOTONFLUX_HPP
+
+/* __ Includes ___________________________________________________________ */
+#include <string>
+#include "GModelSpectral.hpp"
+#include "GModelPar.hpp"
+#include "GEnergy.hpp"
+
+/* __ Forward declarations _______________________________________________ */
+class GRan;
+class GTime;
+class GXmlElement;
+
+
+/***********************************************************************//**
+ * @class GModelSpectralPlawPhotonFlux
+ *
+ * @brief Photon flux normalized power law spectral model class
+ *
+ * This class implements a power law spectrum. The model is defined by
+ *
+ * \f[
+ *    S_{\rm E}(E | t) = {\tt m\_photon\_flux}
+ *    \frac{{\tt m\_index}+1}
+ *         {{\tt e\_max}^{{\tt m\_index}+1} -
+ *          {\tt e\_min}^{{\tt m\_index}+1}}
+ *    E^{\tt m\_index}
+ * \f]
+ *
+ * for \f${\tt m\_index} \ne -1\f$ and
+ *
+ * \f[
+ *    S_{\rm E}(E | t) = 
+ *    \frac{{\tt m\_photon\_flux}}
+ *         {\log {\tt e\_max} - \log {\tt e\_min}}
+ *    E^{\tt m\_index}
+ * \f]
+ *
+ * for \f${\tt m\_index} = -1\f$, where
+ * - \f${\tt e\_min}\f$ is the minimum energy of an interval,
+ * - \f${\tt e\_max}\f$ is the maximum energy of an interval,
+ * - \f${\tt m\_photon\_flux}\f$ is the photon flux between
+ *   \f${\tt e\_min}\f$ and \f${\tt e\_max}\f$, and
+ * - \f${\tt m\_index}\f$ is the spectral index.
+ ***************************************************************************/
+class GModelSpectralPlawPhotonFlux : public GModelSpectral {
+
+public:
+    // Constructors and destructors
+    GModelSpectralPlawPhotonFlux(void);
+    GModelSpectralPlawPhotonFlux(const std::string& type,
+                        const std::string& integral,
+                        const std::string& index,
+                        const std::string& emin,
+                        const std::string& emax);
+    GModelSpectralPlawPhotonFlux(const double&  integral,
+                        const double&  index,
+                        const GEnergy& emin,
+                        const GEnergy& emax);
+    explicit GModelSpectralPlawPhotonFlux(const GXmlElement& xml);
+    GModelSpectralPlawPhotonFlux(const GModelSpectralPlawPhotonFlux& model);
+    virtual ~GModelSpectralPlawPhotonFlux(void);
+
+    // Operators
+    virtual GModelSpectralPlawPhotonFlux& operator=(const GModelSpectralPlawPhotonFlux& model);
+
+    // Implemented pure virtual base class methods
+    virtual void                 clear(void);
+    virtual GModelSpectralPlawPhotonFlux* clone(void) const;
+    virtual std::string          classname(void) const;
+    virtual std::string          type(void) const;
+    virtual double               eval(const GEnergy& srcEng,
+                                      const GTime&   srcTime = GTime()) const;
+    virtual double               eval_gradients(const GEnergy& srcEng,
+                                                const GTime&   srcTime = GTime());
+    virtual double               flux(const GEnergy& emin,
+                                      const GEnergy& emax) const;
+    virtual double               eflux(const GEnergy& emin,
+                                       const GEnergy& emax) const;
+    virtual GEnergy              mc(const GEnergy& emin,
+                                    const GEnergy& emax,
+                                    const GTime&   time,
+                                    GRan&          ran) const;
+    virtual void                 read(const GXmlElement& xml);
+    virtual void                 write(GXmlElement& xml) const;
+    virtual std::string          print(const GChatter& chatter = NORMAL) const;
+
+    // Other methods
+    double  photon_flux(void) const;
+    void    photon_flux(const double& photon_flux);
+    double  index(void) const;
+    void    index(const double& index);
+    GEnergy emin(void) const;
+    void    emin(const GEnergy& emin);
+    GEnergy emax(void) const;
+    void    emax(const GEnergy& emax);
+
+protected:
+    // Protected methods
+    void init_members(void);
+    void copy_members(const GModelSpectralPlawPhotonFlux& model);
+    void free_members(void);
+    void update(const GEnergy& srcEng) const;
+
+    // Protected members
+    std::string     m_type;            //!< Model type
+    GModelPar       m_photon_flux;     //!< Photon flux
+    GModelPar       m_index;           //!< Spectral index
+    GModelPar       m_emin;            //!< Lower energy limit (MeV)
+    GModelPar       m_emax;            //!< Upper energy limit (MeV)
+
+    // Cached members used for pre-computations
+    mutable double  m_log_emin;        //!< Log(emin)
+    mutable double  m_log_emax;        //!< Log(emax)
+    mutable double  m_pow_emin;        //!< emin^(index+1)
+    mutable double  m_pow_emax;        //!< emax^(index+1)
+    mutable double  m_norm;            //!< Power-law normalization (for pivot energy 1 MeV)
+    mutable double  m_g_norm;          //!< Power-law normalization gradient
+    mutable double  m_power;           //!< Power-law factor
+    mutable double  m_last_phflux;     //!< Last photon flux
+    mutable double  m_last_index;      //!< Last spectral index (MeV)
+    mutable GEnergy m_last_emin;       //!< Last lower energy limit
+    mutable GEnergy m_last_emax;       //!< Last upper energy limit
+    mutable GEnergy m_last_energy;     //!< Last source energy
+    mutable double  m_last_value;      //!< Last function value
+    mutable double  m_last_g_phflux;   //!< Last photon flux gradient
+    mutable double  m_last_g_index;    //!< Last spectral index gradient
+};
+
+
+/***********************************************************************//**
+ * @brief Return class name
+ *
+ * @return String containing the class name ("GModelSpectralPlawPhotonFlux").
+ ***************************************************************************/
+inline
+std::string GModelSpectralPlawPhotonFlux::classname(void) const
+{
+    return ("GModelSpectralPlawPhotonFlux");
+}
+
+
+/***********************************************************************//**
+ * @brief Return model type
+ *
+ * @return Model type.
+ *
+ * Returns the type of the spectral power law model.
+ ***************************************************************************/
+inline
+std::string GModelSpectralPlawPhotonFlux::type(void) const
+{
+    return (m_type);
+}
+
+
+/***********************************************************************//**
+ * @brief Return photon flux
+ *
+ * @return Photon flux (ph/cm2/s).
+ *
+ * Returns the photon flux.
+ ***************************************************************************/
+inline
+double GModelSpectralPlawPhotonFlux::photon_flux(void) const
+{
+    return (m_photon_flux.value());
+}
+
+
+/***********************************************************************//**
+ * @brief Set photon flux
+ *
+ * @param[in] photon_flux Photon flux (ph/cm2/s).
+ *
+ * Sets the photon flux.
+ ***************************************************************************/
+inline
+void GModelSpectralPlawPhotonFlux::photon_flux(const double& photon_flux)
+{
+    m_photon_flux.value(photon_flux);
+    return;
+}
+
+
+
+
+/***********************************************************************//**
+ * @brief Return power law index
+ *
+ * @return Power law index.
+ *
+ * Returns the power law index.
+ ***************************************************************************/
+inline
+double GModelSpectralPlawPhotonFlux::index(void) const
+{
+    return (m_index.value());
+}
+
+
+/***********************************************************************//**
+ * @brief Set power law index 
+ *
+ * @param[in] index Power law index.
+ *
+ * Sets the power law index.
+ ***************************************************************************/
+inline
+void GModelSpectralPlawPhotonFlux::index(const double& index)
+{
+    m_index.value(index);
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Return minimum energy
+ *
+ * @return Minimum energy.
+ *
+ * Returns the minimum energy.
+ ***************************************************************************/
+inline
+GEnergy GModelSpectralPlawPhotonFlux::emin(void) const
+{
+    GEnergy energy;
+    energy.MeV(m_emin.value());
+    return energy;
+}
+
+
+/***********************************************************************//**
+ * @brief Set minimum energy
+ *
+ * @param[in] emin Minimum energy.
+ *
+ * Sets the minimum energy.
+ ***************************************************************************/
+inline
+void GModelSpectralPlawPhotonFlux::emin(const GEnergy& emin)
+{
+    m_emin.value(emin.MeV());
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Return maximum energy
+ *
+ * @return Maximum energy.
+ *
+ * Returns the maximum energy.
+ ***************************************************************************/
+inline
+GEnergy GModelSpectralPlawPhotonFlux::emax(void) const
+{
+    GEnergy energy;
+    energy.MeV(m_emax.value());
+    return energy;
+}
+
+
+/***********************************************************************//**
+ * @brief Set maximum energy
+ *
+ * @param[in] emax Maximum energy.
+ *
+ * Sets the maximum energy.
+ ***************************************************************************/
+inline
+void GModelSpectralPlawPhotonFlux::emax(const GEnergy& emax)
+{
+    m_emax.value(emax.MeV());
+    return;
+}
+
+#endif /* GMODELSPECTRALPLAWPHOTONFLUX_HPP */
@@ -215,7 +215,8 @@
 #include "GModelSpectralLogParabola.hpp"
 #include "GModelSpectralNodes.hpp"
 #include "GModelSpectralPlaw.hpp"
-#include "GModelSpectralPlaw2.hpp"
+#include "GModelSpectralPlawPhotonFlux.hpp"
+#include "GModelSpectralPlawenergyFlux.hpp"
 #include "GModelTemporal.hpp"
 #include "GModelTemporalRegistry.hpp"
 #include "GModelTemporalConst.hpp"
@@ -173,7 +173,8 @@ pkginclude_HEADERS = GammaLib.hpp \
                      GModelSpectralLogParabola.hpp \
                      GModelSpectralNodes.hpp \
                      GModelSpectralPlaw.hpp \
-                     GModelSpectralPlaw2.hpp \
+                     GModelSpectralPlawPhotonFlux.hpp \
+                     GModelSpectralPlawEnergyFlux.hpp \
                      GModelTemporal.hpp \
                      GModelTemporalRegistry.hpp \
                      GModelTemporalConst.hpp
+/***************************************************************************
+ *    GModelSpectralPlawEnergyFlux.i - Spectral power law model class      *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2016 by Michael Mayer                                    *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  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   *
+ *  the Free Software Foundation, either version 3 of the License, or      *
+ *  (at your option) any later version.                                    *
+ *                                                                         *
+ *  This program is distributed in the hope that it will be useful,        *
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
+ *  GNU General Public License for more details.                           *
+ *                                                                         *
+ *  You should have received a copy of the GNU General Public License      *
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file GModelSpectralPlawEnergyFlux.i
+ * @brief Energy flux normalized power law spectral model class Python interface definition
+ * @author Michael Mayer
+ */
+%{
+/* Put headers and other declarations here that are needed for compilation */
+#include "GModelSpectralPlawEnergyFlux.hpp"
+%}
+
+
+/***********************************************************************//**
+ * @class GModelSpectralPlawEnergyFlux
+ *
+ * @brief Flux normalized power law spectral model class
+ ***************************************************************************/
+class GModelSpectralPlawEnergyFlux : public GModelSpectral {
+public:
+    // Constructors and destructors
+    GModelSpectralPlawEnergyFlux(void);
+    GModelSpectralPlawEnergyFlux(const std::string& type,
+                                 const std::string& energy_flux,
+                                 const std::string& index,
+                                 const std::string& emin,
+                                 const std::string& emax);
+    explicit GModelSpectralPlawEnergyFlux(const double&  energy_flux,
+                                 const double&  index,
+                                 const GEnergy& emin,
+                                 const GEnergy& emax);
+    explicit GModelSpectralPlawEnergyFlux(const GXmlElement& xml);
+    GModelSpectralPlawEnergyFlux(const GModelSpectralPlawEnergyFlux& model);
+    virtual ~GModelSpectralPlawEnergyFlux(void);
+
+    // Implemented pure virtual methods
+    virtual void                 clear(void);
+    virtual GModelSpectralPlawEnergyFlux* clone(void) const;
+    virtual std::string          classname(void) const;
+    virtual std::string          type(void) const;
+    virtual double               eval(const GEnergy& srcEng,
+                                      const GTime&   srcTime = GTime()) const;
+    virtual double               eval_gradients(const GEnergy& srcEng,
+                                                const GTime&   srcTime = GTime());
+    virtual double               flux(const GEnergy& emin,
+                                      const GEnergy& emax) const;
+    virtual double               eflux(const GEnergy& emin,
+                                       const GEnergy& emax) const;
+    virtual GEnergy              mc(const GEnergy& emin,
+                                    const GEnergy& emax,
+                                    const GTime&   time,
+                                    GRan&          ran) const;
+    virtual void                 read(const GXmlElement& xml);
+    virtual void                 write(GXmlElement& xml) const;
+
+    // Other methods
+    double  energy_flux(void) const;
+    void    energy_flux(const double& energy_flux);
+    double  index(void) const;
+    void    index(const double& index);
+    GEnergy emin(void) const;
+    void    emin(const GEnergy& emin);
+    GEnergy emax(void) const;
+    void    emax(const GEnergy& emax);
+};
+
+
+/***********************************************************************//**
+ * @brief GModelSpectralPlawEnergyFlux class extension
+ ***************************************************************************/
+%extend GModelSpectralPlawEnergyFlux {
+    GModelSpectralPlawEnergyFlux copy() {
+        return (*self);
+    }
+};
 /***************************************************************************
- *         GModelSpectralPlaw2.i - Spectral power law model class          *
+ *    GModelSpectralPlawPhotonFlux.i - Spectral power law model class      *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2011-2016 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
@@ -19,41 +19,41 @@
  *                                                                         *
  ***************************************************************************/
 /**
- * @file GModelSpectralPlaw2.i
+ * @file GModelSpectralPlawPhotonFlux.i
  * @brief Flux normalized power law spectral model class Python interface definition
  * @author Juergen Knoedlseder
  */
 %{
 /* Put headers and other declarations here that are needed for compilation */
-#include "GModelSpectralPlaw2.hpp"
+#include "GModelSpectralPlawPhotonFlux.hpp"
 %}
  
  
 /***********************************************************************//**
- * @class GModelSpectralPlaw2
+ * @class GModelSpectralPlawPhotonFlux
  *
  * @brief Flux normalized power law spectral model class
  ***************************************************************************/
-class GModelSpectralPlaw2 : public GModelSpectral {
+class GModelSpectralPlawPhotonFlux : public GModelSpectral {
 public:
     // Constructors and destructors
-    GModelSpectralPlaw2(void);
-    GModelSpectralPlaw2(const std::string& type,
-                        const std::string& integral,
+    GModelSpectralPlawPhotonFlux(void);
+    GModelSpectralPlawPhotonFlux(const std::string& type,
+                        const std::string& photon_flux,
                         const std::string& index,
                         const std::string& emin,
                         const std::string& emax);
-    GModelSpectralPlaw2(const double&  integral,
+    GModelSpectralPlawPhotonFlux(const double&  photon_flux,
                         const double&  index,
                         const GEnergy& emin,
                         const GEnergy& emax);
-    explicit GModelSpectralPlaw2(const GXmlElement& xml);
-    GModelSpectralPlaw2(const GModelSpectralPlaw2& model);
-    virtual ~GModelSpectralPlaw2(void);
+    explicit GModelSpectralPlawPhotonFlux(const GXmlElement& xml);
+    GModelSpectralPlawPhotonFlux(const GModelSpectralPlawPhotonFlux& model);
+    virtual ~GModelSpectralPlawPhotonFlux(void);
  
     // Implemented pure virtual methods
     virtual void                 clear(void);
-    virtual GModelSpectralPlaw2* clone(void) const;
+    virtual GModelSpectralPlawPhotonFlux* clone(void) const;
     virtual std::string          classname(void) const;
     virtual std::string          type(void) const;
     virtual double               eval(const GEnergy& srcEng,
@@ -72,8 +72,8 @@ public:
     virtual void                 write(GXmlElement& xml) const;
  
     // Other methods
-    double  integral(void) const;
-    void    integral(const double& integral);
+    double  photon_flux(void) const;
+    void    photon_flux(const double& photon_flux);
     double  index(void) const;
     void    index(const double& index);
     GEnergy emin(void) const;
@@ -84,10 +84,10 @@ public:
  
  
 /***********************************************************************//**
- * @brief GModelSpectralPlaw2 class extension
+ * @brief GModelSpectralPlawPhotonFlux class extension
  ***************************************************************************/
-%extend GModelSpectralPlaw2 {
-    GModelSpectralPlaw2 copy() {
+%extend GModelSpectralPlawPhotonFlux {
+    GModelSpectralPlawPhotonFlux copy() {
         return (*self);
     }
 };
@@ -152,7 +152,8 @@
 %include "GModelSpectralLogParabola.i"
 %include "GModelSpectralNodes.i"
 %include "GModelSpectralPlaw.i"
-%include "GModelSpectralPlaw2.i"
+%include "GModelSpectralPlawPhotonFlux.i"
+%include "GModelSpectralPlawEnergyFlux.i"
 %include "GModelTemporal.i"
 %include "GModelTemporalRegistry.i"
 %include "GModelTemporalConst.i"
+/***************************************************************************
+ *   GModelSpectralPlawEnergyFlux.cpp - Spectral power law model class     *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2016 by Michael Mayer                                    *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  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   *
+ *  the Free Software Foundation, either version 3 of the License, or      *
+ *  (at your option) any later version.                                    *
+ *                                                                         *
+ *  This program is distributed in the hope that it will be useful,        *
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
+ *  GNU General Public License for more details.                           *
+ *                                                                         *
+ *  You should have received a copy of the GNU General Public License      *
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file GModelSpectralPlawEnergyFlux.cpp
+ * @brief Energy flux normalized power law spectral model class implementation
+ * @author Michael Mayer
+ */
+
+/* __ Includes ___________________________________________________________ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include <cmath>
+#include "GException.hpp"
+#include "GTools.hpp"
+#include "GMath.hpp"
+#include "GRan.hpp"
+#include "GModelSpectralPlawEnergyFlux.hpp"
+#include "GModelSpectralRegistry.hpp"
+
+/* __ Constants __________________________________________________________ */
+
+/* __ Globals ____________________________________________________________ */
+const GModelSpectralPlawEnergyFlux    g_spectral_plaw_eflux_seed("PowerLaw",
+                                                    "EnergyFlux",
+                                                    "Index",
+                                                    "LowerLimit",
+                                                    "UpperLimit");
+const GModelSpectralRegistry g_spectral_plaw_eflux_registry(&g_spectral_plaw_eflux_seed);
+
+/* __ Method name definitions ____________________________________________ */
+#define G_FLUX                "GModelSpectralPlawEnergyFlux::flux(GEnergy&, GEnergy&)"
+#define G_EFLUX              "GModelSpectralPlawEnergyFlux::eflux(GEnergy&, GEnergy&)"
+#define G_MC     "GModelSpectralPlawEnergyFlux::mc(GEnergy&, GEnergy&, GTime&, GRan&)"
+#define G_READ                      "GModelSpectralPlawEnergyFlux::read(GXmlElement&)"
+#define G_WRITE                    "GModelSpectralPlawEnergyFlux::write(GXmlElement&)"
+
+/* __ Macros _____________________________________________________________ */
+
+/* __ Coding definitions _________________________________________________ */
+
+/* __ Debug definitions __________________________________________________ */
+
+
+/*==========================================================================
+ =                                                                         =
+ =                        Constructors/destructors                         =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Void constructor
+ ***************************************************************************/
+GModelSpectralPlawEnergyFlux::GModelSpectralPlawEnergyFlux(void) : GModelSpectral()
+{
+    // Initialise members
+    init_members();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Model type and parameter name constructor
+ *
+ * @param[in] type Model type.
+ * @param[in] energy_flux Name of energy flux parameter.
+ * @param[in] index Name of index parameter.
+ * @param[in] emin Name of emin parameter.
+ * @param[in] emax Name of emax parameter.
+ ***************************************************************************/
+GModelSpectralPlawEnergyFlux::GModelSpectralPlawEnergyFlux(const std::string& type,
+                                         const std::string& energy_flux,
+                                         const std::string& index,
+                                         const std::string& emin,
+                                         const std::string& emax) :
+                     GModelSpectral()
+{
+    // Initialise members
+    init_members();
+
+    // Set model type
+    m_type = type;
+
+    // Set parameter names
+    m_energy_flux.name(energy_flux);
+    m_index.name(index);
+    m_emin.name(emin);
+    m_emax.name(emax);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Constructor
+ *
+ * @param[in] energy_flux Energy flux (erg/cm2/s).
+ * @param[in] index Power law index.
+ * @param[in] emin Minimum energy.
+ * @param[in] emax Maximum energy.
+ *
+ * Construct a spectral power law from the
+ * - energy flux (in erg/cm2/s),
+ * - spectral index,
+ * - minimum energy and
+ * - maximum energy.
+ ***************************************************************************/
+GModelSpectralPlawEnergyFlux::GModelSpectralPlawEnergyFlux(const double&  energy_flux,
+                                         const double&  index,
+                                         const GEnergy& emin,
+                                         const GEnergy& emax) :
+                     GModelSpectral()
+{
+    // Initialise members
+    init_members();
+
+    // Set parameters
+    m_energy_flux.value(energy_flux);
+    m_index.value(index);
+    m_emin.value(emin.MeV());
+    m_emax.value(emax.MeV());
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief XML constructor
+ *
+ * @param[in] xml XML element.
+ *
+ * Constructs energy flux normalized power law spectral model by extracting
+ * information from an XML element. See the read() method for more
+ * information about the expected structure of the XML element.
+ ***************************************************************************/
+GModelSpectralPlawEnergyFlux::GModelSpectralPlawEnergyFlux(const GXmlElement& xml) :
+                     GModelSpectral()
+{
+    // Initialise members
+    init_members();
+
+    // Read information from XML element
+    read(xml);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy constructor
+ *
+ * @param[in] model Spectral power law model.
+ ***************************************************************************/
+GModelSpectralPlawEnergyFlux::GModelSpectralPlawEnergyFlux(const GModelSpectralPlawEnergyFlux& model) :
+                     GModelSpectral(model)
+{
+    // Initialise members
+    init_members();
+
+    // Copy members
+    copy_members(model);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Destructor
+ ***************************************************************************/
+GModelSpectralPlawEnergyFlux::~GModelSpectralPlawEnergyFlux(void)
+{
+    // Free members
+    free_members();
+
+    // Return
+    return;
+}
+
+
+/*==========================================================================
+ =                                                                         =
+ =                                Operators                                =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Assignment operator
+ *
+ * @param[in] model Spectral power law model.
+ * @return Spectral power law model.
+ ***************************************************************************/
+GModelSpectralPlawEnergyFlux& GModelSpectralPlawEnergyFlux::operator=(const GModelSpectralPlawEnergyFlux& model)
+{
+    // Execute only if object is not identical
+    if (this != &model) {
+
+        // Copy base class members
+        this->GModelSpectral::operator=(model);
+
+        // Free members
+        free_members();
+
+        // Initialise members
+        init_members();
+
+        // Copy members
+        copy_members(model);
+
+    } // endif: object was not identical
+
+    // Return
+    return *this;
+}
+
+
+/*==========================================================================
+ =                                                                         =
+ =                              Public methods                             =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Clear power law model
+ ***************************************************************************/
+void GModelSpectralPlawEnergyFlux::clear(void)
+{
+    // Free class members (base and derived classes, derived class first)
+    free_members();
+    this->GModelSpectral::free_members();
+
+    // Initialise members
+    this->GModelSpectral::init_members();
+    init_members();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Clone power law model
+ *
+ * @return Pointer to deep copy of power law model
+ ***************************************************************************/
+GModelSpectralPlawEnergyFlux* GModelSpectralPlawEnergyFlux::clone(void) const
+{
+    // Clone power law model
+    return new GModelSpectralPlawEnergyFlux(*this);
+}
+
+
+/***********************************************************************//**
+ * @brief Evaluate function
+ *
+ * @param[in] srcEng True photon energy.
+ * @param[in] srcTime True photon arrival time.
+ * @return Model value (ph/cm2/s/MeV).
+ *
+ * Computes
+ *
+ * \f[
+ *    S_{\rm E}(E | t) = {\tt m\_energy\_flux}
+ *    \frac{{\tt m\_index}+2}
+ *         {{\tt e\_max}^{{\tt m\_index}+2} -
+ *          {\tt e\_min}^{{\tt m\_index}+2}}
+ *    E^{\tt m\_index}
+ * \f]
+ *
+ * for \f${\tt m\_index} \ne -2\f$ and
+ *
+ * \f[
+ *    S_{\rm E}(E | t) = 
+ *    \frac{{\tt m\_energy\_flux}}
+ *         {\log {\tt e\_max} - \log {\tt e\_min}}
+ *    E^{\tt m\_index}
+ * \f]
+ *
+ * for \f${\tt m\_index} = -2\f$, where
+ * - \f${\tt e\_min}\f$ is the minimum energy of an interval,
+ * - \f${\tt e\_max}\f$ is the maximum energy of an interval,
+ * - \f${\tt m\_energy\_flux}\f$ is the energy flux between
+ *   \f${\tt e\_min}\f$ and \f${\tt e\_max}\f$, and
+ * - \f${\tt m\_index}\f$ is the spectral index.
+ ***************************************************************************/
+double GModelSpectralPlawEnergyFlux::eval(const GEnergy& srcEng,
+                                 const GTime&   srcTime) const
+{
+    // Update precomputed values
+    update(srcEng);
+
+    // Compute function value
+    double value = m_energy_flux.value() * gammalib::erg2MeV * m_norm * m_power;
+
+    // Compile option: Check for NaN/Inf
+    #if defined(G_NAN_CHECK)
+    if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
+        std::cout << "*** ERROR: GModelSpectralPlawEnergyFlux::eval";
+        std::cout << "(srcEng=" << srcEng;
+        std::cout << ", srcTime=" << srcTime << "):";
+        std::cout << " NaN/Inf encountered";
+        std::cout << " (value=" << value;
+        std::cout << ", energy flux=" << energy_flux();
+        std::cout << ", m_norm=" << m_norm;
+        std::cout << ", m_power=" << m_power;
+        std::cout << ")" << std::endl;
+    }
+    #endif
+
+    // Return
+    return value;
+}
+
+
+/***********************************************************************//**
+ * @brief Evaluate function and gradients
+ *
+ * @param[in] srcEng True photon energy.
+ * @param[in] srcTime True photon arrival time.
+ * @return Model value (ph/cm2/s/MeV).
+ *
+ * Computes
+ *
+ * \f[
+ *    S_{\rm E}(E | t) = {\tt m\_energy\_flux}
+ *    \frac{{\tt m\_index}+2}
+ *         {{\tt e\_max}^{{\tt m\_index}+2} -
+ *          {\tt e\_min}^{{\tt m\_index}+2}}
+ *    E^{\tt m\_index}
+ * \f]
+ *
+ * for \f${\tt m\_index} \ne -2\f$ and
+ *
+ * \f[
+ *    S_{\rm E}(E | t) = 
+ *    \frac{{\tt m\_energy\_flux}}
+ *         {\log {\tt e\_max} - \log {\tt e\_min}}
+ *    E^{\tt m\_index}
+ * \f]
+ *
+ * for \f${\tt m\_index} = -2\f$, where
+ * - \f${\tt e\_min}\f$ is the minimum energy of an interval,
+ * - \f${\tt e\_max}\f$ is the maximum energy of an interval,
+ * - \f${\tt m\_energy\_flux}\f$ is the energy flux between
+ *   \f${\tt e\_min}\f$ and \f${\tt e\_max}\f$, and
+ * - \f${\tt m\_index}\f$ is the spectral index.
+ *
+ * The method also evaluates the partial derivatives of the model with
+ * respect to the parameters using
+ *
+ * \f[
+ *    \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_energy\_flux}} =
+ *      \frac{S_{\rm E}(E | t)}{{\tt m\_energy\_flux}}
+ * \f]
+ *
+ * \f[
+ *    \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_index}} =
+ *      S_{\rm E}(E | t) \,
+ *      \left( \frac{1}{{\tt m\_index}+2} -
+ *             \frac{\log({\tt e\_max}) {\tt e\_max}^{{\tt m\_index}+2} -
+ *                   \log({\tt e\_min}) {\tt e\_min}^{{\tt m\_index}+2}}
+ *                        {{\tt e\_max}^{{\tt m\_index}+2} -
+ *                         {\tt e\_min}^{{\tt m\_index}+2}}
+ *                 + \ln(E) \right)
+ * \f]
+ *
+ * for \f${\tt m\_index} \ne -2\f$ and
+ *
+ * \f[
+ *    \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_index}} =
+ *      S_{\rm E}(E | t) \, \ln(E)
+ * \f]
+ *
+ * for \f${\tt m\_index} = -2\f$.
+ *
+ * No partial derivatives are supported for the energy boundaries.
+ ***************************************************************************/
+double GModelSpectralPlawEnergyFlux::eval_gradients(const GEnergy& srcEng,
+                                           const GTime&   srcTime)
+{
+    // Initialise gradients
+    double g_energy_flux = 0.0;
+    double g_index    = 0.0;
+    
+    // Update precomputed values
+    update(srcEng);
+
+    // Compute function value
+    double value = m_energy_flux.value() * gammalib::erg2MeV * m_norm * m_power;
+
+    // Integral flux gradient
+    if (m_energy_flux.is_free() && m_energy_flux.factor_value() > 0.0) {
+         g_energy_flux = value / m_energy_flux.factor_value();
+    }
+
+    // Index gradient
+    if (m_index.is_free()) {
+        g_index = value * (m_g_norm + gammalib::ln10 * srcEng.log10MeV()) *
+                  m_index.scale();
+    }
+
+    // Set gradients
+    m_energy_flux.factor_gradient(g_energy_flux);
+    m_index.factor_gradient(g_index);
+
+    // Compile option: Check for NaN/Inf
+    #if defined(G_NAN_CHECK)
+    if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
+        std::cout << "*** ERROR: GModelSpectralPlawEnergyFlux::eval_gradients";
+        std::cout << "(srcEng=" << srcEng;
+        std::cout << ", srcTime=" << srcTime << "):";
+        std::cout << " NaN/Inf encountered";
+        std::cout << " (value=" << value;
+        std::cout << ", energy flux=" << energy_flux();
+        std::cout << ", m_norm=" << m_norm;
+        std::cout << ", m_power=" << m_power;
+        std::cout << ", g_energy_flux=" << g_energy_flux;
+        std::cout << ", g_index=" << g_index;
+        std::cout << ")" << std::endl;
+    }
+    #endif
+
+    // Return
+    return value;
+}
+
+
+/***********************************************************************//**
+ * @brief Returns model photon flux between [emin, emax] (units: ph/cm2/s)
+ *
+ * @param[in] emin Minimum photon energy.
+ * @param[in] emax Minimum photon energy.
+ * @return Photon flux (ph/cm2/s).
+ *
+ * Computes
+ *
+ * \f[
+ *    \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) dE
+ * \f]
+ *
+ * where
+ * - [@p emin, @p emax] is an energy interval, and
+ * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
+ ***************************************************************************/
+double GModelSpectralPlawEnergyFlux::flux(const GEnergy& emin,
+                                 const GEnergy& emax) const
+{
+    // Initialise flux
+    double flux = 0.0;
+
+    // Compute only if integration range is valid
+    if (emin < emax) {
+
+    	// Compute power law normalization
+    	double norm = 0.0;
+    	if (index() != -2.0) {
+			double gamma        = m_index.value() + 2.0;
+			double pow_ref_emin = std::pow(this->emin().MeV(), gamma);
+			double pow_ref_emax = std::pow(this->emax().MeV(), gamma);
+			norm = m_energy_flux.value() * gammalib::erg2MeV * gamma / (pow_ref_emax - pow_ref_emin);
+		}
+		else {
+			double log_ref_emin = std::log(this->emin().MeV());
+			double log_ref_emax = std::log(this->emax().MeV());
+			norm = m_energy_flux.value() * gammalib::erg2MeV / (log_ref_emax - log_ref_emin);
+		}
+
+    	// Compute photon flux
+		if (index() != -1.0) {
+			double gamma    = m_index.value() + 1.0;
+			double pow_emin = std::pow(emin.MeV(), gamma);
+			double pow_emax = std::pow(emax.MeV(), gamma);
+			flux = norm / gamma * (pow_emax - pow_emin);
+		}
+
+		// Case B: Index is -1
+		else {
+			double log_emin = std::log(emin.MeV());
+			double log_emax = std::log(emax.MeV());
+			flux = norm * (log_emax - log_emin);
+		}
+
+    } // endif: integration range was valid
+
+    // Return flux
+    return flux;
+}
+
+
+/***********************************************************************//**
+ * @brief Returns model energy flux between [emin, emax] (units: erg/cm2/s)
+ *
+ * @param[in] emin Minimum photon energy.
+ * @param[in] emax Minimum photon energy.
+ * @return Energy flux (erg/cm2/s).
+ *
+ * Computes
+ *
+ * \f[
+ *    \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) E \, dE
+ * \f]
+ *
+ * where
+ * - [@p emin, @p emax] is an energy interval, and
+ * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
+ ***************************************************************************/
+double GModelSpectralPlawEnergyFlux::eflux(const GEnergy& emin,
+                                  const GEnergy& emax) const
+{
+    // Initialise flux
+    double eflux = 0.0;
+
+    // Compute only if integration range is valid
+    if (emin < emax) {
+
+        // Compute power law normalization
+        double norm;
+        if (index() != -2.0) {
+            double gamma        = m_index.value() + 2.0;
+            double pow_ref_emin = std::pow(this->emin().MeV(), gamma);
+            double pow_ref_emax = std::pow(this->emax().MeV(), gamma);
+            double pow_emin     = std::pow(emin.MeV(), gamma);
+            double pow_emax     = std::pow(emax.MeV(), gamma);
+            double factor       = (pow_emax - pow_emin) /
+                                              (pow_ref_emax - pow_ref_emin);
+            eflux               = m_energy_flux.value() * factor;
+        }
+        else {
+        	double log_emin     = std::log(emin.MeV());
+			double log_emax     = std::log(emax.MeV());
+			double log_ref_emin = std::log(this->emin().MeV());
+			double log_ref_emax = std::log(this->emax().MeV());
+			double factor       = (log_emax - log_emin) /
+								  (log_ref_emax - log_ref_emin);
+			eflux               = m_energy_flux.value() * factor;
+        }
+
+    } // endif: integration range was valid
+
+    // Return flux
+    return eflux;
+}
+
+
+/***********************************************************************//**
+ * @brief Returns MC energy between [emin, emax]
+ *
+ * @param[in] emin Minimum photon energy.
+ * @param[in] emax Maximum photon energy.
+ * @param[in] time True photon arrival time.
+ * @param[in,out] ran Random number generator.
+ * @return Energy.
+ *
+ * @exception GException::erange_invalid
+ *            Energy range is invalid (emin < emax required).
+ *
+ * Returns Monte Carlo energy by randomly drawing from a power law.
+ ***************************************************************************/
+GEnergy GModelSpectralPlawEnergyFlux::mc(const GEnergy& emin,
+                                const GEnergy& emax,
+                                const GTime&   time,
+                                GRan&          ran) const
+{
+    // Throw an exception if energy range is invalid
+    if (emin >= emax) {
+        throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(),
+              "Minimum energy < maximum energy required.");
+    }
+
+    // Allocate energy
+    GEnergy energy;
+
+    // Case A: Index is not -1
+    if (index() != -1.0) {
+        double exponent = index() + 1.0;
+        double e_max    = std::pow(emax.MeV(), exponent);
+        double e_min    = std::pow(emin.MeV(), exponent);
+        double u        = ran.uniform();
+        double eng      = (u > 0.0) 
+                          ? std::exp(std::log(u * (e_max - e_min) + e_min) / exponent)
+                          : 0.0;
+        energy.MeV(eng);
+    }
+
+    // Case B: Index is -1
+    else {
+        double e_max = std::log(emax.MeV());
+        double e_min = std::log(emin.MeV());
+        double u     = ran.uniform();
+        double eng   = std::exp(u * (e_max - e_min) + e_min);
+        energy.MeV(eng);
+    }
+
+    // Return energy
+    return energy;
+}
+
+
+/***********************************************************************//**
+ * @brief Read model from XML element
+ *
+ * @param[in] xml XML element.
+ *
+ * Reads the spectral information from an XML element.
+ ***************************************************************************/
+void GModelSpectralPlawEnergyFlux::read(const GXmlElement& xml)
+{
+    // Get parameter pointers
+    const GXmlElement* eflux = gammalib::xml_get_par(G_READ, xml, m_energy_flux.name());
+    const GXmlElement* index = gammalib::xml_get_par(G_READ, xml, m_index.name());
+    const GXmlElement* emin  = gammalib::xml_get_par(G_READ, xml, m_emin.name());
+    const GXmlElement* emax  = gammalib::xml_get_par(G_READ, xml, m_emax.name());
+
+    // Read parameters
+    m_energy_flux.read(*eflux);
+    m_index.read(*index);
+    m_emin.read(*emin);
+    m_emax.read(*emax);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Write model into XML element
+ *
+ * @param[in] xml XML element.
+ *
+ * @exception GException::model_invalid_spectral
+ *            Existing XML element is not of type "PowerLaw"
+ *
+ * Writes the spectral information into an XML element.
+ ***************************************************************************/
+void GModelSpectralPlawEnergyFlux::write(GXmlElement& xml) const
+{
+    // Set model type
+    if (xml.attribute("type") == "") {
+        xml.attribute("type", type());
+    }
+
+    // Verify model type
+    if (xml.attribute("type") != type()) {
+        throw GException::model_invalid_spectral(G_WRITE, xml.attribute("type"),
+              "Spectral model is not of type \""+type()+"\".");
+    }
+
+    // Get XML parameters
+    GXmlElement* eflux = gammalib::xml_need_par(G_WRITE, xml, m_energy_flux.name());
+    GXmlElement* index = gammalib::xml_need_par(G_WRITE, xml, m_index.name());
+    GXmlElement* emin  = gammalib::xml_need_par(G_WRITE, xml, m_emin.name());
+    GXmlElement* emax  = gammalib::xml_need_par(G_WRITE, xml, m_emax.name());
+
+    // Write parameters
+    m_energy_flux.write(*eflux);
+    m_index.write(*index);
+    m_emin.write(*emin);
+    m_emax.write(*emax);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Print power law information
+ *
+ * @param[in] chatter Chattiness (defaults to NORMAL).
+ * @return String containing power law information.
+ ***************************************************************************/
+std::string GModelSpectralPlawEnergyFlux::print(const GChatter& chatter) const
+{
+    // Initialise result string
+    std::string result;
+
+    // Continue only if chatter is not silent
+    if (chatter != SILENT) {
+
+        // Append header
+        result.append("=== GModelSpectralPlawEnergyFlux ===");
+
+        // Append information
+        result.append("\n"+gammalib::parformat("Number of parameters"));
+        result.append(gammalib::str(size()));
+        for (int i = 0; i < size(); ++i) {
+            result.append("\n"+m_pars[i]->print(chatter));
+        }
+
+    } // endif: chatter was not silent
+
+    // Return result
+    return result;
+}
+
+
+/*==========================================================================
+ =                                                                         =
+ =                             Private methods                             =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Initialise class members
+ ***************************************************************************/
+void GModelSpectralPlawEnergyFlux::init_members(void)
+{
+	// Initialise model type
+	m_type = "PowerLaw";
+
+    // Initialise energy flux
+    m_energy_flux.clear();
+    m_energy_flux.name("EnergyFlux");
+    m_energy_flux.unit("erg/cm2/s");
+    m_energy_flux.scale(1.0);
+    m_energy_flux.value(1.0);       // default: 1.0
+    m_energy_flux.range(0.0, 10.0); // range:   [0,10]
+    m_energy_flux.free();
+    m_energy_flux.gradient(0.0);
+    m_energy_flux.has_grad(true);
+
+    // Initialise powerlaw index
+    m_index.clear();
+    m_index.name("Index");
+    m_index.scale(1.0);
+    m_index.value(-2.0);        // default: -2.0
+    m_index.range(-10.0,+10.0); // range:   [-10,+10]
+    m_index.free();
+    m_index.gradient(0.0);
+    m_index.has_grad(true);
+
+    // Initialise lower limit
+    m_emin.clear();
+    m_emin.name("LowerLimit");
+    m_emin.unit("MeV");
+    m_emin.scale(1.0);
+    m_emin.value(100.0);         // default: 100
+    m_emin.range(0.001, 1.0e15); // range:   [0.001, 1e15]
+    m_emin.fix();
+    m_emin.gradient(0.0);
+    m_emin.has_grad(false);
+
+    // Initialise upper limit
+    m_emax.clear();
+    m_emax.name("UpperLimit");
+    m_emax.unit("MeV");
+    m_emax.scale(1.0);
+    m_emax.value(500000.0);      // default: 500000
+    m_emin.range(0.001, 1.0e15); // range:   [0.001, 1e15]
+    m_emax.fix();
+    m_emax.gradient(0.0);
+    m_emax.has_grad(false);
+
+    // Set parameter pointer(s)
+    m_pars.clear();
+    m_pars.push_back(&m_energy_flux);
+    m_pars.push_back(&m_index);
+    m_pars.push_back(&m_emin);
+    m_pars.push_back(&m_emax);
+
+    // Initialise last parameters (for fast computation)
+    m_log_emin     = 0.0;
+    m_log_emax     = 0.0;
+    m_pow_emin     = 0.0;
+    m_pow_emax     = 0.0;
+    m_norm         = 0.0;
+    m_power        = 0.0;
+    m_last_index   = 1000.0;
+    m_last_emin.MeV(0.0);
+    m_last_emax.MeV(0.0);
+    m_last_energy.MeV(0.0);
+    m_last_value   = 0.0;
+    m_last_g_index = 0.0;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy class members
+ *
+ * @param[in] model Spectral power law model.
+ ***************************************************************************/
+void GModelSpectralPlawEnergyFlux::copy_members(const GModelSpectralPlawEnergyFlux& model)
+{
+    // Copy members
+	m_type        = model.m_type;
+    m_energy_flux = model.m_energy_flux;
+    m_index       = model.m_index;
+    m_emin        = model.m_emin;
+    m_emax        = model.m_emax;
+
+    // Set parameter pointer(s)
+    m_pars.clear();
+    m_pars.push_back(&m_energy_flux);
+    m_pars.push_back(&m_index);
+    m_pars.push_back(&m_emin);
+    m_pars.push_back(&m_emax);
+
+    // Copy bookkeeping information
+    m_log_emin     = model.m_log_emin;
+    m_log_emax     = model.m_log_emax;
+    m_pow_emin     = model.m_pow_emin;
+    m_pow_emax     = model.m_pow_emax;
+    m_norm         = model.m_norm;
+    m_power        = model.m_power;
+    m_last_index   = model.m_last_index;
+    m_last_emin    = model.m_last_emin;
+    m_last_emax    = model.m_last_emax;
+    m_last_energy  = model.m_last_energy;
+    m_last_value   = model.m_last_value;
+    m_last_g_index = model.m_last_g_index;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Delete class members
+ ***************************************************************************/
+void GModelSpectralPlawEnergyFlux::free_members(void)
+{
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Update precomputed values
+ *
+ * @param[in] srcEng Source energy
+ ***************************************************************************/
+void GModelSpectralPlawEnergyFlux::update(const GEnergy& srcEng) const
+{
+    // Compute index+1
+    double gamma = index() + 2.0;
+
+    // Change in spectral index?
+    if (index() != m_last_index) {
+
+        // Save actual spectral index
+        m_last_index = index();
+
+        // Change in energy boundaries?
+        if (emin() != m_last_emin || emax() != m_last_emax) {
+            m_log_emin  = std::log(emin().MeV());
+            m_last_emin = emin();
+            m_log_emax  = std::log(emax().MeV());
+            m_last_emax = emax();
+        }
+
+        // Compute normalization factors
+        if (gamma != 0.0) {
+            m_pow_emin  = std::pow(emin().MeV(), gamma);
+            m_pow_emax  = std::pow(emax().MeV(), gamma);
+            double d    = m_pow_emax - m_pow_emin;
+            m_norm      = gamma / d;
+            m_g_norm    = 1.0/gamma -
+                          (m_pow_emax*m_log_emax - m_pow_emin*m_log_emin)/d;
+        }
+        else {
+            m_norm   = 1.0 / (m_log_emax - m_log_emin);
+            m_g_norm = 0.0;
+        }
+
+        // Update power law factor
+        m_power       = std::pow(srcEng.MeV(), index());
+        m_last_energy = srcEng;
+
+    } // endif: change in spectral index
+
+    // ... no change in spectral index
+    else {
+
+        // Change in energy boundaries?
+        if (emin() != m_last_emin || emax() != m_last_emax) {
+
+            // Update energy boundaries
+            m_log_emin  = std::log(emin().MeV());
+            m_last_emin = emin();
+            m_log_emax  = std::log(emax().MeV());
+            m_last_emax = emax();
+        
+            // Compute power law normalization
+            if (gamma != 0.0) {
+                m_pow_emin  = std::pow(emin().MeV(), gamma);
+                m_pow_emax  = std::pow(emax().MeV(), gamma);
+                double d    = m_pow_emax - m_pow_emin;
+                m_norm      = gamma / d;
+                m_g_norm    = 1.0/gamma -
+                              (m_pow_emax*m_log_emax - m_pow_emin*m_log_emin)/d;
+            }
+            else {
+                m_norm = 1.0 / (m_log_emax - m_log_emin);
+                m_g_norm = 0.0;
+            }
+
+        } // endif: change in energy boundaries
+
+        // Change in energy?
+        if (srcEng != m_last_energy) {
+            m_power       = std::pow(srcEng.MeV(), index());
+            m_last_energy = srcEng;
+        }
+
+    } // endelse: no change in spectral index
+
+    // Return
+    return;
+}
 /***************************************************************************
- *         GModelSpectralPlaw2.cpp - Spectral power law model class        *
+ *   GModelSpectralPlawPhotonFlux.cpp - Spectral power law model class     *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2010-2016 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
@@ -19,7 +19,7 @@
  *                                                                         *
  ***************************************************************************/
 /**
- * @file GModelSpectralPlaw2.cpp
+ * @file GModelSpectralPlawPhotonFlux.cpp
  * @brief Flux normalized power law spectral model class implementation
  * @author Juergen Knoedlseder
  */
@@ -33,33 +33,33 @@
 #include "GTools.hpp"
 #include "GMath.hpp"
 #include "GRan.hpp"
-#include "GModelSpectralPlaw2.hpp"
+#include "GModelSpectralPlawPhotonFlux.hpp"
 #include "GModelSpectralRegistry.hpp"
  
 /* __ Constants __________________________________________________________ */
  
 /* __ Globals ____________________________________________________________ */
-const GModelSpectralPlaw2    g_spectral_plaw2_seed1("PowerLaw",
+const GModelSpectralPlawPhotonFlux    g_spectral_plaw_phflux_seed1("PowerLaw",
                                                     "PhotonFlux",
                                                     "Index",
                                                     "LowerLimit",
                                                     "UpperLimit");
-const GModelSpectralRegistry g_spectral_plaw2_registry1(&g_spectral_plaw2_seed1);
+const GModelSpectralRegistry g_spectral_plaw_phflux_registry1(&g_spectral_plaw_phflux_seed1);
 #if defined(G_LEGACY_XML_FORMAT)
-const GModelSpectralPlaw2    g_spectral_plaw2_seed2("PowerLaw2",
+const GModelSpectralPlawPhotonFlux    g_spectral_plaw_phflux_seed2("PowerLaw2",
                                                     "Integral",
                                                     "Index",
                                                     "LowerLimit",
                                                     "UpperLimit");
-const GModelSpectralRegistry g_spectral_plaw2_registry2(&g_spectral_plaw2_seed2);
+const GModelSpectralRegistry g_spectral_plaw_phflux_registry2(&g_spectral_plaw_phflux_seed2);
 #endif
  
 /* __ Method name definitions ____________________________________________ */
-#define G_FLUX                "GModelSpectralPlaw2::flux(GEnergy&, GEnergy&)"
-#define G_EFLUX              "GModelSpectralPlaw2::eflux(GEnergy&, GEnergy&)"
-#define G_MC     "GModelSpectralPlaw2::mc(GEnergy&, GEnergy&, GTime&, GRan&)"
-#define G_READ                      "GModelSpectralPlaw2::read(GXmlElement&)"
-#define G_WRITE                    "GModelSpectralPlaw2::write(GXmlElement&)"
+#define G_FLUX                "GModelSpectralPlawPhotonFlux::flux(GEnergy&, GEnergy&)"
+#define G_EFLUX              "GModelSpectralPlawPhotonFlux::eflux(GEnergy&, GEnergy&)"
+#define G_MC     "GModelSpectralPlawPhotonFlux::mc(GEnergy&, GEnergy&, GTime&, GRan&)"
+#define G_READ                      "GModelSpectralPlawPhotonFlux::read(GXmlElement&)"
+#define G_WRITE                    "GModelSpectralPlawPhotonFlux::write(GXmlElement&)"
  
 /* __ Macros _____________________________________________________________ */
  
@@ -79,7 +79,7 @@ const GModelSpectralRegistry g_spectral_plaw2_registry2(&amp;g_spectral_plaw2_seed2)
  *
  * Constructs empty power law photon flux model.
  ***************************************************************************/
-GModelSpectralPlaw2::GModelSpectralPlaw2(void) : GModelSpectral()
+GModelSpectralPlawPhotonFlux::GModelSpectralPlawPhotonFlux(void) : GModelSpectral()
 {
     // Initialise members
     init_members();
@@ -93,13 +93,13 @@ GModelSpectralPlaw2::GModelSpectralPlaw2(void) : GModelSpectral()
  * @brief Model type and parameter name constructor
  *
  * @param[in] type Model type.
- * @param[in] integral Name of integral parameter.
+ * @param[in] photon_flux Name of photon flux parameter.
  * @param[in] index Name of index parameter.
  * @param[in] emin Name of emin parameter.
  * @param[in] emax Name of emax parameter.
  ***************************************************************************/
-GModelSpectralPlaw2::GModelSpectralPlaw2(const std::string& type,
-                                         const std::string& integral,
+GModelSpectralPlawPhotonFlux::GModelSpectralPlawPhotonFlux(const std::string& type,
+                                         const std::string& photon_flux,
                                          const std::string& index,
                                          const std::string& emin,
                                          const std::string& emax) :
@@ -112,7 +112,7 @@ GModelSpectralPlaw2::GModelSpectralPlaw2(const std::string&amp; type,
     m_type = type;
  
     // Set parameter names
-    m_integral.name(integral);
+    m_photon_flux.name(photon_flux);
     m_index.name(index);
     m_emin.name(emin);
     m_emax.name(emax);
@@ -125,18 +125,18 @@ GModelSpectralPlaw2::GModelSpectralPlaw2(const std::string&amp; type,
 /***********************************************************************//**
  * @brief Constructor
  *
- * @param[in] integral Integral flux (ph/cm2/s).
+ * @param[in] photon_flux Photon flux (ph/cm2/s).
  * @param[in] index Power law index.
  * @param[in] emin Minimum energy.
  * @param[in] emax Maximum energy.
  *
  * Construct a spectral power law from the
- * - integral flux (in ph/cm2/s),
+ * - integral photon flux (in ph/cm2/s),
  * - spectral index,
  * - minimum energy and
  * - maximum energy.
  ***************************************************************************/
-GModelSpectralPlaw2::GModelSpectralPlaw2(const double&  integral,
+GModelSpectralPlawPhotonFlux::GModelSpectralPlawPhotonFlux(const double& photon_flux,
                                          const double&  index,
                                          const GEnergy& emin,
                                          const GEnergy& emax) :
@@ -146,7 +146,7 @@ GModelSpectralPlaw2::GModelSpectralPlaw2(const double&amp;  integral,
     init_members();
  
     // Set parameters
-    m_integral.value(integral);
+    m_photon_flux.value(photon_flux);
     m_index.value(index);
     m_emin.value(emin.MeV());
     m_emax.value(emax.MeV());
@@ -165,7 +165,7 @@ GModelSpectralPlaw2::GModelSpectralPlaw2(const double&amp;  integral,
  * information from an XML element. See the read() method for more
  * information about the expected structure of the XML element.
  ***************************************************************************/
-GModelSpectralPlaw2::GModelSpectralPlaw2(const GXmlElement& xml) :
+GModelSpectralPlawPhotonFlux::GModelSpectralPlawPhotonFlux(const GXmlElement& xml) :
                      GModelSpectral()
 {
     // Initialise members
@@ -184,7 +184,7 @@ GModelSpectralPlaw2::GModelSpectralPlaw2(const GXmlElement&amp; xml) :
  *
  * @param[in] model Spectral power law model.
  ***************************************************************************/
-GModelSpectralPlaw2::GModelSpectralPlaw2(const GModelSpectralPlaw2& model) :
+GModelSpectralPlawPhotonFlux::GModelSpectralPlawPhotonFlux(const GModelSpectralPlawPhotonFlux& model) :
                      GModelSpectral(model)
 {
     // Initialise members
@@ -201,7 +201,7 @@ GModelSpectralPlaw2::GModelSpectralPlaw2(const GModelSpectralPlaw2&amp; model) :
 /***********************************************************************//**
  * @brief Destructor
  ***************************************************************************/
-GModelSpectralPlaw2::~GModelSpectralPlaw2(void)
+GModelSpectralPlawPhotonFlux::~GModelSpectralPlawPhotonFlux(void)
 {
     // Free members
     free_members();
@@ -223,7 +223,7 @@ GModelSpectralPlaw2::~GModelSpectralPlaw2(void)
  * @param[in] model Spectral power law model.
  * @return Spectral power law model.
  ***************************************************************************/
-GModelSpectralPlaw2& GModelSpectralPlaw2::operator=(const GModelSpectralPlaw2& model)
+GModelSpectralPlawPhotonFlux& GModelSpectralPlawPhotonFlux::operator=(const GModelSpectralPlawPhotonFlux& model)
 {
     // Execute only if object is not identical
     if (this != &model) {
@@ -256,7 +256,7 @@ GModelSpectralPlaw2&amp; GModelSpectralPlaw2::operator=(const GModelSpectralPlaw2&amp; m
 /***********************************************************************//**
  * @brief Clear power law model
  ***************************************************************************/
-void GModelSpectralPlaw2::clear(void)
+void GModelSpectralPlawPhotonFlux::clear(void)
 {
     // Free class members (base and derived classes, derived class first)
     free_members();
@@ -276,10 +276,10 @@ void GModelSpectralPlaw2::clear(void)
  *
  * @return Pointer to deep copy of power law model
  ***************************************************************************/
-GModelSpectralPlaw2* GModelSpectralPlaw2::clone(void) const
+GModelSpectralPlawPhotonFlux* GModelSpectralPlawPhotonFlux::clone(void) const
 {
     // Clone power law model
-    return new GModelSpectralPlaw2(*this);
+    return new GModelSpectralPlawPhotonFlux(*this);
 }
  
  
@@ -293,7 +293,7 @@ GModelSpectralPlaw2* GModelSpectralPlaw2::clone(void) const
  * Computes
  *
  * \f[
- *    S_{\rm E}(E | t) = {\tt m\_integral}
+ *    S_{\rm E}(E | t) = {\tt m\_photon\_flux}
  *    \frac{{\tt m\_index}+1}
  *         {{\tt e\_max}^{{\tt m\_index}+1} -
  *          {\tt e\_min}^{{\tt m\_index}+1}}
@@ -304,7 +304,7 @@ GModelSpectralPlaw2* GModelSpectralPlaw2::clone(void) const
  *
  * \f[
  *    S_{\rm E}(E | t) = 
- *    \frac{{\tt m\_integral}}
+ *    \frac{{\tt m\_photon\_flux}}
  *         {\log {\tt e\_max} - \log {\tt e\_min}}
  *    E^{\tt m\_index}
  * \f]
@@ -312,28 +312,28 @@ GModelSpectralPlaw2* GModelSpectralPlaw2::clone(void) const
  * for \f${\tt m\_index} = -1\f$, where
  * - \f${\tt e\_min}\f$ is the minimum energy of an interval,
  * - \f${\tt e\_max}\f$ is the maximum energy of an interval,
- * - \f${\tt m\_integral}\f$ is the integral flux between 
+ * - \f${\tt m\_photon\_flux}\f$ is the photon flux flux between
  *   \f${\tt e\_min}\f$ and \f${\tt e\_max}\f$, and
  * - \f${\tt m\_index}\f$ is the spectral index.
  ***************************************************************************/
-double GModelSpectralPlaw2::eval(const GEnergy& srcEng,
+double GModelSpectralPlawPhotonFlux::eval(const GEnergy& srcEng,
                                  const GTime&   srcTime) const
 {
     // Update precomputed values
     update(srcEng);
  
     // Compute function value
-    double value = m_integral.value() * m_norm * m_power;
+    double value = m_photon_flux.value() * m_norm * m_power;
  
     // Compile option: Check for NaN/Inf
     #if defined(G_NAN_CHECK)
     if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
-        std::cout << "*** ERROR: GModelSpectralPlaw2::eval";
+        std::cout << "*** ERROR: GModelSpectralPlawPhotonFlux::eval";
         std::cout << "(srcEng=" << srcEng;
         std::cout << ", srcTime=" << srcTime << "):";
         std::cout << " NaN/Inf encountered";
         std::cout << " (value=" << value;
-        std::cout << ", integral=" << integral();
+        std::cout << ", photon_flux=" << photon_flux();
         std::cout << ", m_norm=" << m_norm;
         std::cout << ", m_power=" << m_power;
         std::cout << ")" << std::endl;
@@ -355,7 +355,7 @@ double GModelSpectralPlaw2::eval(const GEnergy&amp; srcEng,
  * Computes
  *
  * \f[
- *    S_{\rm E}(E | t) = {\tt m\_integral}
+ *    S_{\rm E}(E | t) = {\tt m\_photon\_flux}
  *    \frac{{\tt m\_index}+1}
  *         {{\tt e\_max}^{{\tt m\_index}+1} -
  *          {\tt e\_min}^{{\tt m\_index}+1}}
@@ -366,7 +366,7 @@ double GModelSpectralPlaw2::eval(const GEnergy&amp; srcEng,
  *
  * \f[
  *    S_{\rm E}(E | t) = 
- *    \frac{{\tt m\_integral}}
+ *    \frac{{\ttm\_photon\_flux}}
  *         {\log {\tt e\_max} - \log {\tt e\_min}}
  *    E^{\tt m\_index}
  * \f]
@@ -374,7 +374,7 @@ double GModelSpectralPlaw2::eval(const GEnergy&amp; srcEng,
  * for \f${\tt m\_index} = -1\f$, where
  * - \f${\tt e\_min}\f$ is the minimum energy of an interval,
  * - \f${\tt e\_max}\f$ is the maximum energy of an interval,
- * - \f${\tt m\_integral}\f$ is the integral flux between 
+ * - \f${\tt m\_photon\_flux}\f$ is the photon flux between
  *   \f${\tt e\_min}\f$ and \f${\tt e\_max}\f$, and
  * - \f${\tt m\_index}\f$ is the spectral index.
  *
@@ -382,8 +382,8 @@ double GModelSpectralPlaw2::eval(const GEnergy&amp; srcEng,
  * respect to the parameters using
  *
  * \f[
- *    \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_integral}} =
- *      \frac{S_{\rm E}(E | t)}{{\tt m\_integral}}
+ *    \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_photon\_flux}} =
+ *      \frac{S_{\rm E}(E | t)}{{\tt m\_photon\_flux}}
  * \f]
  *
  * \f[
@@ -408,22 +408,22 @@ double GModelSpectralPlaw2::eval(const GEnergy&amp; srcEng,
  *
  * No partial derivatives are supported for the energy boundaries.
  ***************************************************************************/
-double GModelSpectralPlaw2::eval_gradients(const GEnergy& srcEng,
+double GModelSpectralPlawPhotonFlux::eval_gradients(const GEnergy& srcEng,
                                            const GTime&   srcTime)
 {
     // Initialise gradients
-    double g_integral = 0.0;
+    double g_photon_flux = 0.0;
     double g_index    = 0.0;
  
     // Update precomputed values
     update(srcEng);
  
     // Compute function value
-    double value = m_integral.value() * m_norm * m_power;
+    double value = m_photon_flux.value() * m_norm * m_power;
  
-    // Integral flux gradient
-    if (m_integral.is_free() && m_integral.factor_value() > 0.0) {
-         g_integral = value / m_integral.factor_value();
+    // Photon flux gradient
+    if (m_photon_flux.is_free() && m_photon_flux.factor_value() > 0.0) {
+         g_photon_flux = value / m_photon_flux.factor_value();
     }
  
     // Index gradient
@@ -433,21 +433,21 @@ double GModelSpectralPlaw2::eval_gradients(const GEnergy&amp; srcEng,
     }
  
     // Set gradients
-    m_integral.factor_gradient(g_integral);
+    m_photon_flux.factor_gradient(g_photon_flux);
     m_index.factor_gradient(g_index);
  
     // Compile option: Check for NaN/Inf
     #if defined(G_NAN_CHECK)
     if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
-        std::cout << "*** ERROR: GModelSpectralPlaw2::eval_gradients";
+        std::cout << "*** ERROR: GModelSpectralPlawPhotonFlux::eval_gradients";
         std::cout << "(srcEng=" << srcEng;
         std::cout << ", srcTime=" << srcTime << "):";
         std::cout << " NaN/Inf encountered";
         std::cout << " (value=" << value;
-        std::cout << ", integral=" << integral();
+        std::cout << ", photon_flux=" << photon_flux();
         std::cout << ", m_norm=" << m_norm;
         std::cout << ", m_power=" << m_power;
-        std::cout << ", g_integral=" << g_integral;
+        std::cout << ", g_photon_flux=" << g_photon_flux;
         std::cout << ", g_index=" << g_index;
         std::cout << ")" << std::endl;
     }
@@ -475,7 +475,7 @@ double GModelSpectralPlaw2::eval_gradients(const GEnergy&amp; srcEng,
  * - [@p emin, @p emax] is an energy interval, and
  * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
  ***************************************************************************/
-double GModelSpectralPlaw2::flux(const GEnergy& emin,
+double GModelSpectralPlawPhotonFlux::flux(const GEnergy& emin,
                                  const GEnergy& emax) const
 {
     // Initialise flux
@@ -493,7 +493,7 @@ double GModelSpectralPlaw2::flux(const GEnergy&amp; emin,
             double pow_ref_emax = std::pow(this->emax().MeV(), gamma);
             double factor       = (pow_emax - pow_emin) /
                                   (pow_ref_emax - pow_ref_emin);
-            flux                = m_integral.value() * factor;
+            flux                = m_photon_flux.value() * factor;
         }
  
         // Case B: Index is -1
@@ -504,7 +504,7 @@ double GModelSpectralPlaw2::flux(const GEnergy&amp; emin,
             double log_ref_emax = std::log(this->emax().MeV());
             double factor       = (log_emax - log_emin) /
                                   (log_ref_emax - log_ref_emin);
-            flux                = m_integral.value() * factor;
+            flux                = m_photon_flux.value() * factor;
         }
  
     } // endif: integration range was valid
@@ -531,7 +531,7 @@ double GModelSpectralPlaw2::flux(const GEnergy&amp; emin,
  * - [@p emin, @p emax] is an energy interval, and
  * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
  ***************************************************************************/
-double GModelSpectralPlaw2::eflux(const GEnergy& emin,
+double GModelSpectralPlawPhotonFlux::eflux(const GEnergy& emin,
                                   const GEnergy& emax) const
 {
     // Initialise flux
@@ -546,12 +546,12 @@ double GModelSpectralPlaw2::eflux(const GEnergy&amp; emin,
             double gamma        = m_index.value() + 1.0;
             double pow_ref_emin = std::pow(this->emin().MeV(), gamma);
             double pow_ref_emax = std::pow(this->emax().MeV(), gamma);
-            norm = m_integral.value() * gamma / (pow_ref_emax - pow_ref_emin);
+            norm = m_photon_flux.value() * gamma / (pow_ref_emax - pow_ref_emin);
         }
         else {
             double log_ref_emin = std::log(this->emin().MeV());
             double log_ref_emax = std::log(this->emax().MeV());
-            norm = m_integral.value() / (log_ref_emax - log_ref_emin);
+            norm = m_photon_flux.value() / (log_ref_emax - log_ref_emin);
         }
  
         // Compute energy flux
@@ -593,7 +593,7 @@ double GModelSpectralPlaw2::eflux(const GEnergy&amp; emin,
  *
  * Returns Monte Carlo energy by randomly drawing from a power law.
  ***************************************************************************/
-GEnergy GModelSpectralPlaw2::mc(const GEnergy& emin,
+GEnergy GModelSpectralPlawPhotonFlux::mc(const GEnergy& emin,
                                 const GEnergy& emax,
                                 const GTime&   time,
                                 GRan&          ran) const
@@ -640,16 +640,16 @@ GEnergy GModelSpectralPlaw2::mc(const GEnergy&amp; emin,
  *
  * Reads the spectral information from an XML element.
  ***************************************************************************/
-void GModelSpectralPlaw2::read(const GXmlElement& xml)
+void GModelSpectralPlawPhotonFlux::read(const GXmlElement& xml)
 {
     // Get parameter pointers
-    const GXmlElement* integral = gammalib::xml_get_par(G_READ, xml, m_integral.name());
-    const GXmlElement* index    = gammalib::xml_get_par(G_READ, xml, m_index.name());
-    const GXmlElement* emin     = gammalib::xml_get_par(G_READ, xml, m_emin.name());
-    const GXmlElement* emax     = gammalib::xml_get_par(G_READ, xml, m_emax.name());
+    const GXmlElement* phflux = gammalib::xml_get_par(G_READ, xml, m_photon_flux.name());
+    const GXmlElement* index  = gammalib::xml_get_par(G_READ, xml, m_index.name());
+    const GXmlElement* emin   = gammalib::xml_get_par(G_READ, xml, m_emin.name());
+    const GXmlElement* emax   = gammalib::xml_get_par(G_READ, xml, m_emax.name());
  
     // Read parameters
-    m_integral.read(*integral);
+    m_photon_flux.read(*phflux);
     m_index.read(*index);
     m_emin.read(*emin);
     m_emax.read(*emax);
@@ -665,11 +665,11 @@ void GModelSpectralPlaw2::read(const GXmlElement&amp; xml)
  * @param[in] xml XML element.
  *
  * @exception GException::model_invalid_spectral
- *            Existing XML element is not of type "PowerLaw2"
+ *            Existing XML element is not of type "PowerLaw"
  *
  * Writes the spectral information into an XML element.
  ***************************************************************************/
-void GModelSpectralPlaw2::write(GXmlElement& xml) const
+void GModelSpectralPlawPhotonFlux::write(GXmlElement& xml) const
 {
     // Set model type
     if (xml.attribute("type") == "") {
@@ -683,13 +683,13 @@ void GModelSpectralPlaw2::write(GXmlElement&amp; xml) const
     }
  
     // Get XML parameters
-    GXmlElement* integral = gammalib::xml_need_par(G_WRITE, xml, m_integral.name());
-    GXmlElement* index    = gammalib::xml_need_par(G_WRITE, xml, m_index.name());
-    GXmlElement* emin     = gammalib::xml_need_par(G_WRITE, xml, m_emin.name());
-    GXmlElement* emax     = gammalib::xml_need_par(G_WRITE, xml, m_emax.name());
+    GXmlElement* phflux = gammalib::xml_need_par(G_WRITE, xml, m_photon_flux.name());
+    GXmlElement* index  = gammalib::xml_need_par(G_WRITE, xml, m_index.name());
+    GXmlElement* emin   = gammalib::xml_need_par(G_WRITE, xml, m_emin.name());
+    GXmlElement* emax   = gammalib::xml_need_par(G_WRITE, xml, m_emax.name());
  
     // Write parameters
-    m_integral.write(*integral);
+    m_photon_flux.write(*phflux);
     m_index.write(*index);
     m_emin.write(*emin);
     m_emax.write(*emax);
@@ -705,7 +705,7 @@ void GModelSpectralPlaw2::write(GXmlElement&amp; xml) const
  * @param[in] chatter Chattiness (defaults to NORMAL).
  * @return String containing power law information.
  ***************************************************************************/
-std::string GModelSpectralPlaw2::print(const GChatter& chatter) const
+std::string GModelSpectralPlawPhotonFlux::print(const GChatter& chatter) const
 {
     // Initialise result string
     std::string result;
@@ -714,7 +714,7 @@ std::string GModelSpectralPlaw2::print(const GChatter&amp; chatter) const
     if (chatter != SILENT) {
  
         // Append header
-        result.append("=== GModelSpectralPlaw2 ===");
+        result.append("=== GModelSpectralPlawPhotonFlux ===");
  
         // Append information
         result.append("\n"+gammalib::parformat("Number of parameters"));
@@ -739,21 +739,21 @@ std::string GModelSpectralPlaw2::print(const GChatter&amp; chatter) const
 /***********************************************************************//**
  * @brief Initialise class members
  ***************************************************************************/
-void GModelSpectralPlaw2::init_members(void)
+void GModelSpectralPlawPhotonFlux::init_members(void)
 {
     // Initialise model type
     m_type = "PowerLaw";
  
-    // Initialise integral flux
-    m_integral.clear();
-    m_integral.name("PhotonFlux");
-    m_integral.unit("ph/cm2/s");
-    m_integral.scale(1.0);
-    m_integral.value(1.0);       // default: 1.0
-    m_integral.range(0.0, 10.0); // range:   [0,10]
-    m_integral.free();
-    m_integral.gradient(0.0);
-    m_integral.has_grad(true);
+    // Initialise photon flux
+    m_photon_flux.clear();
+    m_photon_flux.name("PhotonFlux");
+    m_photon_flux.unit("ph/cm2/s");
+    m_photon_flux.scale(1.0);
+    m_photon_flux.value(1.0);       // default: 1.0
+    m_photon_flux.range(0.0, 10.0); // range:   [0,10]
+    m_photon_flux.free();
+    m_photon_flux.gradient(0.0);
+    m_photon_flux.has_grad(true);
  
     // Initialise powerlaw index
     m_index.clear();
@@ -789,7 +789,7 @@ void GModelSpectralPlaw2::init_members(void)
  
     // Set parameter pointer(s)
     m_pars.clear();
-    m_pars.push_back(&m_integral);
+    m_pars.push_back(&m_photon_flux);
     m_pars.push_back(&m_index);
     m_pars.push_back(&m_emin);
     m_pars.push_back(&m_emax);
@@ -801,13 +801,13 @@ void GModelSpectralPlaw2::init_members(void)
     m_pow_emax        = 0.0;
     m_norm            = 0.0;
     m_power           = 0.0;
-    m_last_integral   = 0.0;
+    m_last_phflux     = 0.0;
     m_last_index      = 1000.0;
     m_last_emin.MeV(0.0);
     m_last_emax.MeV(0.0);
     m_last_energy.MeV(0.0);
     m_last_value      = 0.0;
-    m_last_g_integral = 0.0;
+    m_last_g_phflux   = 0.0;
     m_last_g_index    = 0.0;
  
     // Return
@@ -820,18 +820,18 @@ void GModelSpectralPlaw2::init_members(void)
  *
  * @param[in] model Spectral power law model.
  ***************************************************************************/
-void GModelSpectralPlaw2::copy_members(const GModelSpectralPlaw2& model)
+void GModelSpectralPlawPhotonFlux::copy_members(const GModelSpectralPlawPhotonFlux& model)
 {
     // Copy members
-    m_type     = model.m_type;
-    m_integral = model.m_integral;
-    m_index    = model.m_index;
-    m_emin     = model.m_emin;
-    m_emax     = model.m_emax;
+    m_type        = model.m_type;
+    m_photon_flux = model.m_photon_flux;
+    m_index       = model.m_index;
+    m_emin        = model.m_emin;
+    m_emax        = model.m_emax;
  
     // Set parameter pointer(s)
     m_pars.clear();
-    m_pars.push_back(&m_integral);
+    m_pars.push_back(&m_photon_flux);
     m_pars.push_back(&m_index);
     m_pars.push_back(&m_emin);
     m_pars.push_back(&m_emax);
@@ -843,13 +843,13 @@ void GModelSpectralPlaw2::copy_members(const GModelSpectralPlaw2&amp; model)
     m_pow_emax        = model.m_pow_emax;
     m_norm            = model.m_norm;
     m_power           = model.m_power;
-    m_last_integral   = model.m_last_integral;
+    m_last_phflux   = model.m_last_phflux;
     m_last_index      = model.m_last_index;
     m_last_emin       = model.m_last_emin;
     m_last_emax       = model.m_last_emax;
     m_last_energy     = model.m_last_energy;
     m_last_value      = model.m_last_value;
-    m_last_g_integral = model.m_last_g_integral;
+    m_last_g_phflux   = model.m_last_g_phflux;
     m_last_g_index    = model.m_last_g_index;
  
     // Return
@@ -860,7 +860,7 @@ void GModelSpectralPlaw2::copy_members(const GModelSpectralPlaw2&amp; model)
 /***********************************************************************//**
  * @brief Delete class members
  ***************************************************************************/
-void GModelSpectralPlaw2::free_members(void)
+void GModelSpectralPlawPhotonFlux::free_members(void)
 {
     // Return
     return;
@@ -872,7 +872,7 @@ void GModelSpectralPlaw2::free_members(void)
  *
  * @param[in] srcEng Source energy
  ***************************************************************************/
-void GModelSpectralPlaw2::update(const GEnergy& srcEng) const
+void GModelSpectralPlawPhotonFlux::update(const GEnergy& srcEng) const
 {
     // Compute index+1
     double gamma = index() + 1.0;
@@ -56,7 +56,8 @@ sources = GModelPar.cpp \
           GModelSpectralLogParabola.cpp \
           GModelSpectralNodes.cpp \
           GModelSpectralPlaw.cpp \
-          GModelSpectralPlaw2.cpp \
+          GModelSpectralPlawPhotonFlux.cpp \
+          GModelSpectralPlawEnergyFlux.cpp \
           GModelTemporal.cpp \
           GModelTemporalRegistry.cpp \
           GModelTemporalConst.cpp \
+<?xml version="1.0" standalone="no"?>
+<source_library title="source library">
+  <source name="Crab" type="PointSource">
+    <spectrum type="PowerLaw">
+      <parameter scale="1e-07" name="EnergyFlux" min="1e-07" max="1000.0"    value="1.0" free="1"/>
+      <parameter scale="1.0"   name="Index"      min="-5.0"  max="+5.0"      value="-2.0" free="1"/>
+      <parameter scale="1.0"   name="LowerLimit" min="10.0"  max="1000000.0" value="100.0" free="0"/>
+      <parameter scale="1.0"   name="UpperLimit" min="10.0"  max="1000000.0" value="500000.0" free="0"/>
+    </spectrum>
+    <spatialModel type="PointSource">
+      <parameter free="0" max="360" min="-360" name="RA" scale="1" value="83.6331" />
+      <parameter free="0" max="90" min="-90" name="DEC" scale="1" value="22.0145" />
+    </spatialModel>
+  </source>
+</source_library>
@@ -58,7 +58,8 @@ void TestGModel::set(void)
     m_xml_model_point_const        = datadir + "/model_point_const.xml";
     m_xml_model_point_gauss        = datadir + "/model_point_gauss.xml";
     m_xml_model_point_plaw         = datadir + "/model_point_plaw.xml";
-    m_xml_model_point_plaw2        = datadir + "/model_point_plaw2.xml";
+    m_xml_model_point_plaw_phflux  = datadir + "/model_point_plaw_phflux.xml";
+    m_xml_model_point_plaw_eflux   = datadir + "/model_point_plaw_eflux.xml";
     m_xml_model_point_eplaw        = datadir + "/model_point_eplaw.xml";
     m_xml_model_point_bplaw        = datadir + "/model_point_bplaw.xml";
     m_xml_model_point_supeplaw     = datadir + "/model_point_supeplaw.xml";
@@ -124,8 +125,10 @@ void TestGModel::set(void)
            "Test GModelSpectralGauss");
     append(static_cast<pfunction>(&TestGModel::test_plaw),
            "Test GModelSpectralPlaw");
-    append(static_cast<pfunction>(&TestGModel::test_plaw2),
-           "Test GModelSpectralPlaw2");
+    append(static_cast<pfunction>(&TestGModel::test_plaw_phflux),
+           "Test GModelSpectralPlawPhotonFkux");
+    append(static_cast<pfunction>(&TestGModel::test_plaw_eflux),
+               "Test GModelSpectralPlawEnergyFlux");
     append(static_cast<pfunction>(&TestGModel::test_eplaw),
            "Test GModelSpectralExpPlaw");
     append(static_cast<pfunction>(&TestGModel::test_supeplaw),
@@ -1551,35 +1554,35 @@ void TestGModel::test_plaw(void)
  
  
 /***********************************************************************//**
- * @brief Test GModelSpectralPlaw2 class
+ * @brief Test GModelSpectralPlawPhotonFlux class
  ***************************************************************************/
-void TestGModel::test_plaw2(void)
+void TestGModel::test_plaw_phflux(void)
 {
     // Test void constructor
-    GModelSpectralPlaw2 model1;
+    GModelSpectralPlawPhotonFlux model1;
     test_value(model1.type(), "PowerLaw", "Check type of void model");
  
     // Test value constructor
-    GModelSpectralPlaw2 model2(2.0, -2.1, GEnergy(10.0, "MeV"), GEnergy(100.0, "MeV"));
-    test_value(model2.integral(), 2.0);
+    GModelSpectralPlawPhotonFlux model2(2.0, -2.1, GEnergy(10.0, "MeV"), GEnergy(100.0, "MeV"));
+    test_value(model2.photon_flux(), 2.0);
     test_value(model2.index(), -2.1);
     test_value(model2.emin().MeV(), 10.0);
     test_value(model2.emax().MeV(), 100.0);
-    
+
     // Test XML constructor and value
-    GXml                xml(m_xml_model_point_plaw2);
+    GXml                xml(m_xml_model_point_plaw_phflux);
     GXmlElement*        element = xml.element(0)->element(0)->element("spectrum", 0);
-    GModelSpectralPlaw2 model3(*element);
+    GModelSpectralPlawPhotonFlux model3(*element);
     test_value(model3.size(), 4);
     test_value(model3.type(), "PowerLaw", "Check model type");
-    test_value(model3.integral(), 1.0e-7);
+    test_value(model3.photon_flux(), 1.0e-7);
     test_value(model3.index(), -2.0);
     test_value(model3.emin().MeV(), 100.0);
     test_value(model3.emax().MeV(), 500000.0);
  
     // Test integral method
-    model3.integral(2.1e-7);
-    test_value(model3.integral(), 2.1e-7);
+    model3.photon_flux(2.1e-7);
+    test_value(model3.photon_flux(), 2.1e-7);
  
     // Test index method
     model3.index(-2.3);
@@ -1612,6 +1615,67 @@ void TestGModel::test_plaw2(void)
  
  
 /***********************************************************************//**
+ * @brief Test GModelSpectralPlawEnergyFlux class
+ ***************************************************************************/
+void TestGModel::test_plaw_eflux(void)
+{
+    // Test void constructor
+    GModelSpectralPlawEnergyFlux model1;
+    test_value(model1.type(), "PowerLaw", "Check type of void model");
+
+    // Test value constructor
+    GModelSpectralPlawEnergyFlux model2(2.0, -2.1, GEnergy(10.0, "MeV"), GEnergy(100.0, "MeV"));
+    test_value(model2.energy_flux(), 2.0);
+    test_value(model2.index(), -2.1);
+    test_value(model2.emin().MeV(), 10.0);
+    test_value(model2.emax().MeV(), 100.0);
+    
+    // Test XML constructor and value
+    GXml                xml(m_xml_model_point_plaw_eflux);
+    GXmlElement*        element = xml.element(0)->element(0)->element("spectrum", 0);
+    GModelSpectralPlawEnergyFlux model3(*element);
+    test_value(model3.size(), 4);
+    test_value(model3.type(), "PowerLaw", "Check model type");
+    test_value(model3.energy_flux(), 1.0e-7);
+    test_value(model3.index(), -2.0);
+    test_value(model3.emin().MeV(), 100.0);
+    test_value(model3.emax().MeV(), 500000.0);
+
+    // Test integral method
+    model3.energy_flux(2.1e-7);
+    test_value(model3.energy_flux(), 2.1e-7);
+
+    // Test index method
+    model3.index(-2.3);
+    test_value(model3.index(), -2.3);
+
+    // Test emin method
+    model3.emin(GEnergy(10.0, "MeV"));
+    test_value(model3.emin().MeV(), 10.0);
+
+    // Test emax method
+    model3.emax(GEnergy(10.0, "MeV"));
+    test_value(model3.emax().MeV(), 10.0);
+
+    // Test operator access
+    const char* strarray[] = {"EnergyFlux", "Index", "LowerLimit", "UpperLimit"};
+    for (int i = 0; i < 4; ++i) {
+        std::string keyname(strarray[i]);
+        model3[keyname].remove_range(); // To allow setting of any value
+        model3[keyname].value(2.1);
+        model3[keyname].error(1.9);
+        model3[keyname].gradient(0.8);
+        test_value(model3[keyname].value(), 2.1);
+        test_value(model3[keyname].error(), 1.9);
+        test_value(model3[keyname].gradient(), 0.8);
+    }
+
+    // Exit test
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Test GModelSpectralExpPlaw class
  ***************************************************************************/
 void TestGModel::test_eplaw(void)
@@ -2094,15 +2158,16 @@ void TestGModel::test_spatial_model(void)
 void TestGModel::test_spectral_model(void)
 {
     // Test spectral models XML interface
-    test_xml_model("GModelSpectralConst",        m_xml_model_point_const);
-    test_xml_model("GModelSpectralPlaw",         m_xml_model_point_plaw);
-    test_xml_model("GModelSpectralPlaw2",        m_xml_model_point_plaw2);
-    test_xml_model("GModelSpectralExpPaw",       m_xml_model_point_eplaw);
-    test_xml_model("GModelSpectralBrokenPlaw",   m_xml_model_point_bplaw);
-    test_xml_model("GModelSpectralSuperExpPlaw", m_xml_model_point_supeplaw);
-    test_xml_model("GModelSpectralLogParabola",  m_xml_model_point_logparabola);
-    test_xml_model("GModelSpectralNodes",        m_xml_model_point_nodes);
-    test_xml_model("GModelSpectralFunc",         m_xml_model_point_filefct);
+    test_xml_model("GModelSpectralConst",          m_xml_model_point_const);
+    test_xml_model("GModelSpectralPlaw",           m_xml_model_point_plaw);
+    test_xml_model("GModelSpectralPlawPhotonFlux", m_xml_model_point_plaw_phflux);
+    test_xml_model("GModelSpectralPlawEnergyFlux", m_xml_model_point_plaw_eflux);
+    test_xml_model("GModelSpectralExpPaw",         m_xml_model_point_eplaw);
+    test_xml_model("GModelSpectralBrokenPlaw",     m_xml_model_point_bplaw);
+    test_xml_model("GModelSpectralSuperExpPlaw",   m_xml_model_point_supeplaw);
+    test_xml_model("GModelSpectralLogParabola",    m_xml_model_point_logparabola);
+    test_xml_model("GModelSpectralNodes",          m_xml_model_point_nodes);
+    test_xml_model("GModelSpectralFunc",           m_xml_model_point_filefct);
  
     // Return
     return;
@@ -2823,13 +2888,13 @@ void TestGModel::test_legacy_model_point_plaw2(void)
     GModels models(m_xml_legacy_point_plaw2);
  
     // Extract spectral component
-    GModelSpectralPlaw2* spectral = static_cast<GModelSpectralPlaw2*>(
+    GModelSpectralPlawPhotonFlux* spectral = static_cast<GModelSpectralPlawPhotonFlux*>(
                                       static_cast<GModelSky*>(
                                         models[0])->spectral());
  
     // Test model values
     test_value(spectral->type(), "PowerLaw2", "Check model type");
-    test_value(spectral->integral(), 1.0e-7, 1.0e-7, "Check integral photon flux");
+    test_value(spectral->photon_flux(), 1.0e-7, 1.0e-7, "Check integral photon flux");
     test_value(spectral->index(), -2.0, 1.0e-7, "Check index");
     test_value(spectral->emin().MeV(), 100.0, 1.0e-7, "Check minimum energy");
     test_value(spectral->emax().MeV(), 500000.0, 1.0e-7, "Check maximum energy");
@@ -2840,12 +2905,12 @@ void TestGModel::test_legacy_model_point_plaw2(void)
     models.load("test_xml_legacy_point_plaw2.xml");
  
     // Extract spectral component
-    spectral = static_cast<GModelSpectralPlaw2*>(
+    spectral = static_cast<GModelSpectralPlawPhotonFlux*>(
                  static_cast<GModelSky*>(models[0])->spectral());
  
     // Re-test model values
     test_value(spectral->type(), "PowerLaw2", "Check model type");
-    test_value(spectral->integral(), 1.0e-7, 1.0e-7, "Check integral photon flux");
+    test_value(spectral->photon_flux(), 1.0e-7, 1.0e-7, "Check integral photon flux");
     test_value(spectral->index(), -2.0, 1.0e-7, "Check index");
     test_value(spectral->emin().MeV(), 100.0, 1.0e-7, "Check minimum energy");
     test_value(spectral->emax().MeV(), 500000.0, 1.0e-7, "Check maximum energy");
@@ -62,7 +62,8 @@ public:
     void                test_const(void);
     void                test_gauss(void);
     void                test_plaw(void);
-    void                test_plaw2(void);
+    void                test_plaw_phflux(void);
+    void                test_plaw_eflux(void);
     void                test_eplaw(void);
     void                test_bplaw(void);
     void                test_supeplaw(void);
@@ -102,7 +103,8 @@ private:
     std::string m_xml_model_point_const;
     std::string m_xml_model_point_gauss;
     std::string m_xml_model_point_plaw;
-    std::string m_xml_model_point_plaw2;
+    std::string m_xml_model_point_plaw_phflux;
+    std::string m_xml_model_point_plaw_eflux;
     std::string m_xml_model_point_eplaw;
     std::string m_xml_model_point_bplaw;
     std::string m_xml_model_point_supeplaw;