diff --git a/include/GModelSpectralConst.hpp b/include/GModelSpectralConst.hpp
index 219b595..3dc0a3b 100644
--- a/include/GModelSpectralConst.hpp
+++ b/include/GModelSpectralConst.hpp
@@ -12,7 +12,7 @@
  ***************************************************************************/
 /**
  * @file GModelSpectralConst.hpp
- * @brief GModelSpectralConst class interface definition.
+ * @brief Constant spectral model class interface definition
  * @author J. Knodlseder
  */
 
@@ -30,15 +30,15 @@
 /***********************************************************************//**
  * @class GModelSpectralConst
  *
- * @brief Spectral constant interface definition
+ * @brief Constant spectral model class
  *
  * This class implements a constant as the spectral component of the
- * gamma-ray data model. The function is defined as
+ * gamma-ray sky model. The function is defined as
  * \f[I(E)=norm\f]
  * where
  * \f$norm\f$ is the normalization constant.
  ***************************************************************************/
-class GModelSpectralConst  : public GModelSpectral {
+class GModelSpectralConst : public GModelSpectral {
 
 public:
     // Constructors and destructors
@@ -48,20 +48,19 @@ public:
     virtual ~GModelSpectralConst(void);
 
     // Operators
-    GModelSpectralConst& operator= (const GModelSpectralConst& model);
+    virtual GModelSpectralConst& operator=(const GModelSpectralConst& model);
 
     // Implemented pure virtual methods
-    void                 clear(void);
-    GModelSpectralConst* clone(void) const;
-    int                  size(void) const { return m_npars; }
-    std::string          type(void) const { return "ConstantValue"; }
-    double               eval(const GEnergy& srcEng);
-    double               eval_gradients(const GEnergy& srcEng);
-    double               flux(const GEnergy& emin, const GEnergy& emax) const;
-    GEnergy              mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
-    void                 read(const GXmlElement& xml);
-    void                 write(GXmlElement& xml) const;
-    std::string          print(void) const;
+    virtual void                 clear(void);
+    virtual GModelSpectralConst* clone(void) const;
+    virtual std::string          type(void) const { return "ConstantValue"; }
+    virtual double               eval(const GEnergy& srcEng) const;
+    virtual double               eval_gradients(const GEnergy& srcEng) const;
+    virtual double               flux(const GEnergy& emin, const GEnergy& emax) const;
+    virtual GEnergy              mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
+    virtual void                 read(const GXmlElement& xml);
+    virtual void                 write(GXmlElement& xml) const;
+    virtual std::string          print(void) const;
 
     // Other methods
     double norm(void) const { return m_norm.real_value(); }
@@ -72,13 +71,8 @@ protected:
     void copy_members(const GModelSpectralConst& model);
     void free_members(void);
 
-    // Implemented pure virtual methods
-    GModelPar** par(void) { return m_par; }
-
     // Protected members
-    int        m_npars;     //!< Number of parameters
-    GModelPar* m_par[1];    //!< Pointers to parameters
-    GModelPar  m_norm;      //!< Normalization factor
+    GModelPar m_norm;      //!< Normalization factor
 };
 
 #endif /* GMODELSPECTRALCONST_HPP */
diff --git a/include/GModelSpectralExpPlaw.hpp b/include/GModelSpectralExpPlaw.hpp
index 17dbb76..32d187f 100644
--- a/include/GModelSpectralExpPlaw.hpp
+++ b/include/GModelSpectralExpPlaw.hpp
@@ -12,7 +12,7 @@
  ***************************************************************************/
 /**
  * @file GModelSpectralExpPlaw.hpp
- * @brief GModelSpectralExpPlaw class interface definition.
+ * @brief Exponential cut off power law spectral class interface definition
  * @author J. Knodlseder
  */
 
@@ -20,6 +20,7 @@
 #define GMODELSPECTRALEXPPLAW_HPP
 
 /* __ Includes ___________________________________________________________ */
+#include <string>
 #include "GModelPar.hpp"
 #include "GModelSpectral.hpp"
 #include "GEnergy.hpp"
@@ -29,10 +30,10 @@
 /***********************************************************************//**
  * @class GModelSpectralExpPlaw
  *
- * @brief Spectral power law interface definition.
+ * @brief Exponential cut off power law spectral class
  *
  * This class implements a power law as the spectral component of the
- * gamma-ray data model. The power law is defined as
+ * gamma-ray sky model. The power law is defined as
  * \f[I(E)=norm (E/pivot)^{index} \exp(-E/ecut)\f]
  * where
  * \f$norm\f$ is the normalization or prefactor,
@@ -40,7 +41,7 @@
  * \f$index\f$ is the spectral index, and
  * \f$ecut\f$ is the cut off energy.
  ***************************************************************************/
-class GModelSpectralExpPlaw  : public GModelSpectral {
+class GModelSpectralExpPlaw : public GModelSpectral {
 
 public:
     // Constructors and destructors
@@ -51,20 +52,19 @@ public:
     virtual ~GModelSpectralExpPlaw(void);
 
     // Operators
-    GModelSpectralExpPlaw& operator= (const GModelSpectralExpPlaw& model);
+    virtual GModelSpectralExpPlaw& operator= (const GModelSpectralExpPlaw& model);
 
     // Implemented pure virtual methods
-    void                   clear(void);
-    GModelSpectralExpPlaw* clone(void) const;
-    int                    size(void) const { return m_npars; }
-    std::string            type(void) const { return "ExpCutoff"; }
-    double                 eval(const GEnergy& srcEng);
-    double                 eval_gradients(const GEnergy& srcEng);
-    double                 flux(const GEnergy& emin, const GEnergy& emax) const;
-    GEnergy                mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
-    void                   read(const GXmlElement& xml);
-    void                   write(GXmlElement& xml) const;
-    std::string            print(void) const;
+    virtual void                   clear(void);
+    virtual GModelSpectralExpPlaw* clone(void) const;
+    virtual std::string            type(void) const { return "ExpCutoff"; }
+    virtual double                 eval(const GEnergy& srcEng) const;
+    virtual double                 eval_gradients(const GEnergy& srcEng) const;
+    virtual double                 flux(const GEnergy& emin, const GEnergy& emax) const;
+    virtual GEnergy                mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
+    virtual void                   read(const GXmlElement& xml);
+    virtual void                   write(GXmlElement& xml) const;
+    virtual std::string            print(void) const;
 
     // Other methods
     void   autoscale(void);
@@ -79,16 +79,11 @@ protected:
     void copy_members(const GModelSpectralExpPlaw& model);
     void free_members(void);
 
-    // Implemented pure virtual methods
-    GModelPar** par(void) { return m_par; }
-
     // Protected members
-    int        m_npars;           //!< Number of parameters
-    GModelPar* m_par[4];          //!< Pointers to parameters
-    GModelPar  m_norm;            //!< Normalization factor
-    GModelPar  m_index;           //!< Spectral index
-    GModelPar  m_ecut;            //!< Exponential cut off energy
-    GModelPar  m_pivot;           //!< Pivot energy
+    GModelPar m_norm;            //!< Normalization factor
+    GModelPar m_index;           //!< Spectral index
+    GModelPar m_ecut;            //!< Exponential cut off energy
+    GModelPar m_pivot;           //!< Pivot energy
 };
 
 #endif /* GMODELSPECTRALEXPPLAW_HPP */
diff --git a/include/GModelSpectralFunc.hpp b/include/GModelSpectralFunc.hpp
index 019dd35..b71f7f4 100644
--- a/include/GModelSpectralFunc.hpp
+++ b/include/GModelSpectralFunc.hpp
@@ -32,7 +32,7 @@
 /***********************************************************************//**
  * @class GModelSpectralFunc
  *
- * @brief Spectral function interface definition
+ * @brief Spectral function model class
  *
  * This class implements an arbitrary function as the spectral component of
  * the model. The function is defined as
@@ -40,7 +40,7 @@
  * where
  * \f$norm\f$ is the normalization of the function.
  ***************************************************************************/
-class GModelSpectralFunc  : public GModelSpectral {
+class GModelSpectralFunc : public GModelSpectral {
 
 public:
     // Constructors and destructors
@@ -51,20 +51,19 @@ public:
     virtual ~GModelSpectralFunc(void);
 
     // Operators
-    GModelSpectralFunc& operator= (const GModelSpectralFunc& model);
+    virtual GModelSpectralFunc& operator=(const GModelSpectralFunc& model);
 
     // Implemented pure virtual methods
-    void                clear(void);
-    GModelSpectralFunc* clone(void) const;
-    int                 size(void) const { return m_npars; }
-    std::string         type(void) const { return "FileFunction"; }
-    double              eval(const GEnergy& srcEng);
-    double              eval_gradients(const GEnergy& srcEng);
-    double              flux(const GEnergy& emin, const GEnergy& emax) const;
-    GEnergy             mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
-    void                read(const GXmlElement& xml);
-    void                write(GXmlElement& xml) const;
-    std::string         print(void) const;
+    virtual void                clear(void);
+    virtual GModelSpectralFunc* clone(void) const;
+    virtual std::string         type(void) const { return "FileFunction"; }
+    virtual double              eval(const GEnergy& srcEng) const;
+    virtual double              eval_gradients(const GEnergy& srcEng) const;
+    virtual double              flux(const GEnergy& emin, const GEnergy& emax) const;
+    virtual GEnergy             mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
+    virtual void                read(const GXmlElement& xml);
+    virtual void                write(GXmlElement& xml) const;
+    virtual std::string         print(void) const;
 
     // Other methods
     double norm(void) const { return m_norm.real_value(); }
@@ -76,12 +75,7 @@ protected:
     void free_members(void);
     void load_nodes(const std::string& filename);
 
-    // Implemented pure virtual methods
-    GModelPar** par(void) { return m_par; }
-
     // Protected members
-    int                 m_npars;     //!< Number of parameters
-    GModelPar*          m_par[1];    //!< Pointers to parameters
     GModelPar           m_norm;      //!< Normalization factor
     GNodeArray          m_nodes;     //!< Nodes of function
     std::vector<double> m_values;    //!< Function values at nodes
diff --git a/include/GModelSpectralPlaw.hpp b/include/GModelSpectralPlaw.hpp
index 974d76a..ed60bd8 100644
--- a/include/GModelSpectralPlaw.hpp
+++ b/include/GModelSpectralPlaw.hpp
@@ -12,7 +12,7 @@
  ***************************************************************************/
 /**
  * @file GModelSpectralPlaw.hpp
- * @brief GModelSpectralPlaw class interface definition.
+ * @brief Power law spectral model class interface definition
  * @author J. Knodlseder
  */
 
@@ -20,6 +20,7 @@
 #define GMODELSPECTRALPLAW_HPP
 
 /* __ Includes ___________________________________________________________ */
+#include <string>
 #include "GModelPar.hpp"
 #include "GModelSpectral.hpp"
 #include "GEnergy.hpp"
@@ -29,17 +30,17 @@
 /***********************************************************************//**
  * @class GModelSpectralPlaw
  *
- * @brief Spectral power law interface definition.
+ * @brief Power law spectral model class
  *
  * This class implements a power law as the spectral component of the
- * gamma-ray data model. The power law is defined as
+ * gamma-ray sky model. The power law is defined as
  * \f[I(E)=norm (E/pivot)^{index}\f]
  * where
  * \f$norm\f$ is the normalization or prefactor,
  * \f$pivot\f$ is the pivot energy, and
  * \f$index\f$ is the spectral index.
  ***************************************************************************/
-class GModelSpectralPlaw  : public GModelSpectral {
+class GModelSpectralPlaw : public GModelSpectral {
 
 public:
     // Constructors and destructors
@@ -50,20 +51,19 @@ public:
     virtual ~GModelSpectralPlaw(void);
 
     // Operators
-    GModelSpectralPlaw& operator= (const GModelSpectralPlaw& model);
+    virtual GModelSpectralPlaw& operator=(const GModelSpectralPlaw& model);
 
     // Implemented pure virtual methods
-    void                clear(void);
-    GModelSpectralPlaw* clone(void) const;
-    int                 size(void) const { return m_npars; }
-    std::string         type(void) const { return "PowerLaw"; }
-    double              eval(const GEnergy& srcEng);
-    double              eval_gradients(const GEnergy& srcEng);
-    double              flux(const GEnergy& emin, const GEnergy& emax) const;
-    GEnergy             mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
-    void                read(const GXmlElement& xml);
-    void                write(GXmlElement& xml) const;
-    std::string         print(void) const;
+    virtual void                clear(void);
+    virtual GModelSpectralPlaw* clone(void) const;
+    virtual std::string         type(void) const { return "PowerLaw"; }
+    virtual double              eval(const GEnergy& srcEng) const;
+    virtual double              eval_gradients(const GEnergy& srcEng) const;
+    virtual double              flux(const GEnergy& emin, const GEnergy& emax) const;
+    virtual GEnergy             mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
+    virtual void                read(const GXmlElement& xml);
+    virtual void                write(GXmlElement& xml) const;
+    virtual std::string         print(void) const;
 
     // Other methods
     void   autoscale(void);
@@ -77,15 +77,10 @@ protected:
     void copy_members(const GModelSpectralPlaw& model);
     void free_members(void);
 
-    // Implemented pure virtual methods
-    GModelPar** par(void) { return m_par; }
-
     // Protected members
-    int        m_npars;           //!< Number of parameters
-    GModelPar* m_par[3];          //!< Pointers to parameters
-    GModelPar  m_norm;            //!< Normalization factor
-    GModelPar  m_index;           //!< Spectral index
-    GModelPar  m_pivot;           //!< Pivot energy
+    GModelPar m_norm;            //!< Normalization factor
+    GModelPar m_index;           //!< Spectral index
+    GModelPar m_pivot;           //!< Pivot energy
 };
 
 #endif /* GMODELSPECTRALPLAW_HPP */
diff --git a/include/GModelSpectralPlaw2.hpp b/include/GModelSpectralPlaw2.hpp
index e4aae5a..f8367c9 100644
--- a/include/GModelSpectralPlaw2.hpp
+++ b/include/GModelSpectralPlaw2.hpp
@@ -12,7 +12,7 @@
  ***************************************************************************/
 /**
  * @file GModelSpectralPlaw2.hpp
- * @brief GModelSpectralPlaw2 class interface definition.
+ * @brief Flux normalized power law spectral model class interface definition
  * @author J. Knodlseder
  */
 
@@ -20,6 +20,7 @@
 #define GMODELSPECTRALPLAW2_HPP
 
 /* __ Includes ___________________________________________________________ */
+#include <string>
 #include "GModelPar.hpp"
 #include "GModelSpectral.hpp"
 #include "GEnergy.hpp"
@@ -29,10 +30,10 @@
 /***********************************************************************//**
  * @class GModelSpectralPlaw2
  *
- * @brief Spectral power law interface definition.
+ * @brief Flux normalized power law spectral model class
  *
  * This class implements a power law as the spectral component of the
- * gamma-ray data model.
+ * gamma-ray sky model.
  * The power law function is defined as
  * \f[I(E)=integral (index+1)/(emax^{index+1}-emin^{index+1}) E^{index}\f]
  * for \f$index \ne -1\f$ and
@@ -43,7 +44,7 @@
  * \f$emin\f$ is the lower energy limit, and
  * \f$emax\f$ is the upper energy limit.
  ***************************************************************************/
-class GModelSpectralPlaw2  : public GModelSpectral {
+class GModelSpectralPlaw2 : public GModelSpectral {
 
 public:
     // Constructors and destructors
@@ -54,20 +55,19 @@ public:
     virtual ~GModelSpectralPlaw2(void);
 
     // Operators
-    GModelSpectralPlaw2& operator= (const GModelSpectralPlaw2& model);
+    virtual GModelSpectralPlaw2& operator= (const GModelSpectralPlaw2& model);
 
     // Implemented pure virtual methods
-    void                 clear(void);
-    GModelSpectralPlaw2* clone(void) const;
-    int                  size(void) const { return m_npars; }
-    std::string          type(void) const { return "PowerLaw2"; }
-    double               eval(const GEnergy& srcEng);
-    double               eval_gradients(const GEnergy& srcEng);
-    double               flux(const GEnergy& emin, const GEnergy& emax) const;
-    GEnergy              mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
-    void                 read(const GXmlElement& xml);
-    void                 write(GXmlElement& xml) const;
-    std::string          print(void) const;
+    virtual void                 clear(void);
+    virtual GModelSpectralPlaw2* clone(void) const;
+    virtual std::string          type(void) const { return "PowerLaw2"; }
+    virtual double               eval(const GEnergy& srcEng) const;
+    virtual double               eval_gradients(const GEnergy& srcEng) const;
+    virtual double               flux(const GEnergy& emin, const GEnergy& emax) const;
+    virtual GEnergy              mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
+    virtual void                 read(const GXmlElement& xml);
+    virtual void                 write(GXmlElement& xml) const;
+    virtual std::string          print(void) const;
 
     // Other methods
     double integral(void) const { return m_integral.real_value(); }
@@ -82,31 +82,26 @@ protected:
     void free_members(void);
     void update(const GEnergy& srcEng);
 
-    // Implemented pure virtual methods
-    GModelPar** par(void) { return m_par; }
-
     // Protected members
-    int        m_npars;           //!< Number of parameters
-    GModelPar* m_par[4];          //!< Pointers to parameters
-    GModelPar  m_integral;        //!< Integral flux
-    GModelPar  m_index;           //!< Spectral index
-    GModelPar  m_emin;            //!< Lower energy limit (MeV)
-    GModelPar  m_emax;            //!< Upper energy limit (MeV)
-    double     m_log_emin;        //!< Log(emin)
-    double     m_log_emax;        //!< Log(emax)
-    double     m_pow_emin;        //!< emin^(index+1)
-    double     m_pow_emax;        //!< emax^(index+1)
-    double     m_norm;            //!< Power-law normalization (for pivot energy 1 MeV)
-    double     m_g_norm;          //!< Power-law normalization gradient
-    double     m_power;           //!< Power-law factor
-    double     m_last_integral;   //!< Last integral flux
-    double     m_last_index;      //!< Last spectral index (MeV)
-    double     m_last_emin;       //!< Last lower energy limit (MeV)
-    double     m_last_emax;       //!< Last upper energy limit (MeV)
-    GEnergy    m_last_energy;     //!< Last source energy
-    double     m_last_value;      //!< Last function value
-    double     m_last_g_integral; //!< Last integral flux gradient
-    double     m_last_g_index;    //!< Last spectral index gradient
+    GModelPar m_integral;        //!< Integral flux
+    GModelPar m_index;           //!< Spectral index
+    GModelPar m_emin;            //!< Lower energy limit (MeV)
+    GModelPar m_emax;            //!< Upper energy limit (MeV)
+    double    m_log_emin;        //!< Log(emin)
+    double    m_log_emax;        //!< Log(emax)
+    double    m_pow_emin;        //!< emin^(index+1)
+    double    m_pow_emax;        //!< emax^(index+1)
+    double    m_norm;            //!< Power-law normalization (for pivot energy 1 MeV)
+    double    m_g_norm;          //!< Power-law normalization gradient
+    double    m_power;           //!< Power-law factor
+    double    m_last_integral;   //!< Last integral flux
+    double    m_last_index;      //!< Last spectral index (MeV)
+    double    m_last_emin;       //!< Last lower energy limit (MeV)
+    double    m_last_emax;       //!< Last upper energy limit (MeV)
+    GEnergy   m_last_energy;     //!< Last source energy
+    double    m_last_value;      //!< Last function value
+    double    m_last_g_integral; //!< Last integral flux gradient
+    double    m_last_g_index;    //!< Last spectral index gradient
 };
 
 #endif /* GMODELSPECTRALPLAW2_HPP */
diff --git a/include/GModelTemporalConst.hpp b/include/GModelTemporalConst.hpp
index d8788bb..c15753d 100644
--- a/include/GModelTemporalConst.hpp
+++ b/include/GModelTemporalConst.hpp
@@ -12,7 +12,7 @@
  ***************************************************************************/
 /**
  * @file GModelTemporalConst.hpp
- * @brief GModelTemporalConst class interface definition.
+ * @brief Constant temporal model class interface definition
  * @author J. Knodlseder
  */
 
@@ -20,22 +20,24 @@
 #define GMODELTEMPORALCONST_HPP
 
 /* __ Includes ___________________________________________________________ */
+#include <string>
 #include "GModelPar.hpp"
 #include "GModelTemporal.hpp"
 #include "GTime.hpp"
+#include "GXmlElement.hpp"
 
 
 /***********************************************************************//**
  * @class GModelTemporalConst
  *
- * @brief Interface definition for the constant model class
+ * @brief Constant temporal model class
  *
  * This class implements the temporal component of the factorised model for
  * a model that is constant in time. It has a single parameter norm that is
  * used for scaling, and that normally is set to 1. The parameter is not
  * supposed to be fitted.
  ***************************************************************************/
-class GModelTemporalConst  : public GModelTemporal {
+class GModelTemporalConst : public GModelTemporal {
 
 public:
     // Constructors and destructors
@@ -44,21 +46,19 @@ public:
     virtual ~GModelTemporalConst(void);
 
     // Operators
-    GModelTemporalConst& operator= (const GModelTemporalConst& model);
+    virtual GModelTemporalConst& operator= (const GModelTemporalConst& model);
 
     // Implemented virtual methods
-    void                 clear(void);
-    GModelTemporalConst* clone(void) const;
-    int                  size(void) const { return m_npars; }
-    std::string          type(void) const { return "Constant"; }
-    double               eval(const GTime& srcTime);
-    double               eval_gradients(const GTime& srcTime);
-    GTimes               mc(const double& rate,
-                            const GTime& tmin, const GTime& tmax,
-                            GRan& ran);
-    void                 read(const GXmlElement& xml);
-    void                 write(GXmlElement& xml) const;
-    std::string          print(void) const;
+    virtual void                 clear(void);
+    virtual GModelTemporalConst* clone(void) const;
+    virtual std::string          type(void) const { return "Constant"; }
+    virtual double               eval(const GTime& srcTime) const;
+    virtual double               eval_gradients(const GTime& srcTime) const;
+    virtual GTimes               mc(const double& rate, const GTime& tmin,
+                                    const GTime& tmax, GRan& ran) const;
+    virtual void                 read(const GXmlElement& xml);
+    virtual void                 write(GXmlElement& xml) const;
+    virtual std::string          print(void) const;
 
     // Other methods
     double norm(void) const { return m_norm.real_value(); }
@@ -69,13 +69,8 @@ protected:
     void copy_members(const GModelTemporalConst& model);
     void free_members(void);
 
-    // Implemented pure virtual methods
-    GModelPar** par(void) { return m_par; }
-
     // Protected members
-    int        m_npars;   //!< Number of parameters
-    GModelPar* m_par[1];  //!< Pointers to parameters
-    GModelPar  m_norm;    //!< Constant
+    GModelPar m_norm;    //!< Constant
 };
 
 #endif /* GMODELTEMPORALCONST_HPP */
diff --git a/pyext/GModelSpectralExpPlaw.i b/pyext/GModelSpectralExpPlaw.i
index d4f2dec..89fa01f 100644
--- a/pyext/GModelSpectralExpPlaw.i
+++ b/pyext/GModelSpectralExpPlaw.i
@@ -12,7 +12,7 @@
  ***************************************************************************/
 /**
  * @file GModelSpectralExpPlaw.i
- * @brief GModelSpectralExpPlaw class python interface
+ * @brief Exponential cut off power law spectral class Python interface definition
  * @author J. Knodlseder
  */
 %{
@@ -25,9 +25,9 @@
 /***********************************************************************//**
  * @class GModelSpectralExpPlaw
  *
- * @brief Implements an exponentially cut off power law.
+ * @brief Exponential cut off power law spectral class
  ***************************************************************************/
-class GModelSpectralExpPlaw  : public GModelSpectral {
+class GModelSpectralExpPlaw : public GModelSpectral {
 public:
     // Constructors and destructors
     GModelSpectralExpPlaw(void);
@@ -38,23 +38,22 @@ public:
     virtual ~GModelSpectralExpPlaw(void);
 
     // Implemented pure virtual methods
-    void                   clear(void);
-    GModelSpectralExpPlaw* clone(void) const;
-    int                    size(void) const { return m_npars; }
-    std::string            type(void) const { return "ExpCutoff"; }
-    double                 eval(const GEnergy& srcEng);
-    double                 eval_gradients(const GEnergy& srcEng);
-    double                 flux(const GEnergy& emin, const GEnergy& emax) const;
-    GEnergy                mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
-    void                   read(const GXmlElement& xml);
-    void                   write(GXmlElement& xml) const;
+    virtual void                   clear(void);
+    virtual GModelSpectralExpPlaw* clone(void) const;
+    virtual std::string            type(void) const;
+    virtual double                 eval(const GEnergy& srcEng) const;
+    virtual double                 eval_gradients(const GEnergy& srcEng) const;
+    virtual double                 flux(const GEnergy& emin, const GEnergy& emax) const;
+    virtual GEnergy                mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
+    virtual void                   read(const GXmlElement& xml);
+    virtual void                   write(GXmlElement& xml) const;
 
     // Other methods
     void   autoscale(void);
-    double norm(void) const { return m_norm.real_value(); }
-    double index(void) const { return m_index.real_value(); }
-    double ecut(void) const { return m_ecut.real_value(); }
-    double pivot(void) const { return m_pivot.real_value(); }
+    double norm(void) const;
+    double index(void) const;
+    double ecut(void) const;
+    double pivot(void) const;
 };
 
 
diff --git a/pyext/GModelSpectralPlaw.i b/pyext/GModelSpectralPlaw.i
index c0870f6..348e6ba 100644
--- a/pyext/GModelSpectralPlaw.i
+++ b/pyext/GModelSpectralPlaw.i
@@ -1,5 +1,5 @@
 /***************************************************************************
- *   GModelSpectralPlaw.i  -  Spectral power law model class python I/F    *
+ *         GModelSpectralPlaw.i  -  Spectral power law model class         *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2009-2011 by Jurgen Knodlseder                           *
  * ----------------------------------------------------------------------- *
@@ -12,7 +12,7 @@
  ***************************************************************************/
 /**
  * @file GModelSpectralPlaw.i
- * @brief GModelSpectralPlaw class python interface
+ * @brief Power law spectral model class Python interface definition
  * @author J. Knodlseder
  */
 %{
@@ -25,9 +25,9 @@
 /***********************************************************************//**
  * @class GModelSpectralPlaw
  *
- * @brief Powerlaw python interface definition
+ * @brief Power law spectral model class
  ***************************************************************************/
-class GModelSpectralPlaw  : public GModelSpectral {
+class GModelSpectralPlaw : public GModelSpectral {
 public:
     // Constructors and destructors
     GModelSpectralPlaw(void);
@@ -37,22 +37,21 @@ public:
     virtual ~GModelSpectralPlaw(void);
 
     // Implemented pure virtual methods
-    void                clear(void);
-    GModelSpectralPlaw* clone(void) const;
-    int                 size(void) const { return m_npars; }
-    std::string         type(void) const { return "PowerLaw"; }
-    double              eval(const GEnergy& srcEng);
-    double              eval_gradients(const GEnergy& srcEng);
-    double              flux(const GEnergy& emin, const GEnergy& emax) const;
-    GEnergy             mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
-    void                read(const GXmlElement& xml);
-    void                write(GXmlElement& xml) const;
+    virtual void                clear(void);
+    virtual GModelSpectralPlaw* clone(void) const;
+    virtual std::string         type(void) const;
+    virtual double              eval(const GEnergy& srcEng) const;
+    virtual double              eval_gradients(const GEnergy& srcEng) const;
+    virtual double              flux(const GEnergy& emin, const GEnergy& emax) const;
+    virtual GEnergy             mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
+    virtual void                read(const GXmlElement& xml);
+    virtual void                write(GXmlElement& xml) const;
 
     // Other methods
     void   autoscale(void);
-    double norm(void) const { return m_norm.real_value(); }
-    double index(void) const { return m_index.real_value(); }
-    double pivot(void) const { return m_pivot.real_value(); }
+    double norm(void) const;
+    double index(void) const;
+    double pivot(void) const;
 };
 
 
diff --git a/pyext/GModelSpectralPlaw2.i b/pyext/GModelSpectralPlaw2.i
index b4e2854..f21b594 100644
--- a/pyext/GModelSpectralPlaw2.i
+++ b/pyext/GModelSpectralPlaw2.i
@@ -1,5 +1,5 @@
 /***************************************************************************
- *   GModelSpectralPlaw2.i  -  Spectral power law model class python I/F   *
+ *        GModelSpectralPlaw2.i  -  Spectral power law model class         *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2011 by Jurgen Knodlseder                                *
  * ----------------------------------------------------------------------- *
@@ -12,7 +12,7 @@
  ***************************************************************************/
 /**
  * @file GModelSpectralPlaw2.i
- * @brief GModelSpectralPlaw2 class python interface
+ * @brief Flux normalized power law spectral model class Python interface definition
  * @author J. Knodlseder
  */
 %{
@@ -25,9 +25,9 @@
 /***********************************************************************//**
  * @class GModelSpectralPlaw2
  *
- * @brief Powerlaw python interface definition
+ * @brief Flux normalized power law spectral model class
  ***************************************************************************/
-class GModelSpectralPlaw2  : public GModelSpectral {
+class GModelSpectralPlaw2 : public GModelSpectral {
 public:
     // Constructors and destructors
     GModelSpectralPlaw2(void);
@@ -37,16 +37,15 @@ public:
     virtual ~GModelSpectralPlaw2(void);
 
     // Implemented pure virtual methods
-    void                 clear(void);
-    GModelSpectralPlaw2* clone(void) const;
-    int                  size(void) const;
-    std::string          type(void) const;
-    double               eval(const GEnergy& srcEng);
-    double               eval_gradients(const GEnergy& srcEng);
-    double               flux(const GEnergy& emin, const GEnergy& emax) const;
-    GEnergy              mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
-    void                 read(const GXmlElement& xml);
-    void                 write(GXmlElement& xml) const;
+    virtual void                 clear(void);
+    virtual GModelSpectralPlaw2* clone(void) const;
+    virtual std::string          type(void) const;
+    virtual double               eval(const GEnergy& srcEng) const;
+    virtual double               eval_gradients(const GEnergy& srcEng) const;
+    virtual double               flux(const GEnergy& emin, const GEnergy& emax) const;
+    virtual GEnergy              mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
+    virtual void                 read(const GXmlElement& xml);
+    virtual void                 write(GXmlElement& xml) const;
 
     // Other methods
     double integral(void) const;
diff --git a/pyext/GModelTemporalConst.i b/pyext/GModelTemporalConst.i
index c43e955..b23d934 100644
--- a/pyext/GModelTemporalConst.i
+++ b/pyext/GModelTemporalConst.i
@@ -1,5 +1,5 @@
 /***************************************************************************
- *   GModelTemporalConst.i  -  Temporal constant model class python I/F    *
+ *         GModelTemporalConst.i  -  Temporal constant model class         *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2009-2011 by Jurgen Knodlseder                           *
  * ----------------------------------------------------------------------- *
@@ -12,7 +12,7 @@
  ***************************************************************************/
 /**
  * @file GModelTemporalConst.i
- * @brief GModelTemporalConst class python interface
+ * @brief Constant temporal model class Python interface definition
  * @author J. Knodlseder
  */
 %{
@@ -25,7 +25,7 @@
 /***********************************************************************//**
  * @class GModelTemporalConst
  *
- * @brief Python interface definition for the constant model class.
+ * @brief Constant temporal model class
  ***************************************************************************/
 class GModelTemporalConst  : public GModelTemporal {
 public:
@@ -35,20 +35,18 @@ public:
     virtual ~GModelTemporalConst(void);
 
     // Implemented virtual methods
-    void                 clear(void);
-    GModelTemporalConst* clone(void) const;
-    int                  size(void) const { return m_npars; }
-    std::string          type(void) const { return "Constant"; }
-    double               eval(const GTime& srcTime);
-    double               eval_gradients(const GTime& srcTime);
-    GTimes               mc(const double& rate,
-                            const GTime& tmin, const GTime& tmax,
-                            GRan& ran);
-    void                 read(const GXmlElement& xml);
-    void                 write(GXmlElement& xml) const;
+    virtual void                 clear(void);
+    virtual GModelTemporalConst* clone(void) const;
+    virtual std::string          type(void) const;
+    virtual double               eval(const GTime& srcTime) const;
+    virtual double               eval_gradients(const GTime& srcTime) const;
+    virtual GTimes               mc(const double& rate, const GTime& tmin,
+                                    const GTime& tmax, GRan& ran) const;
+    virtual void                 read(const GXmlElement& xml);
+    virtual void                 write(GXmlElement& xml) const;
 
     // Other methods
-    double norm(void) const { return m_norm.real_value(); }
+    double norm(void) const;
 };
 
 
diff --git a/src/model/GModelSpatialCube.cpp b/src/model/GModelSpatialCube.cpp
index b03a17f..054c5b7 100644
--- a/src/model/GModelSpatialCube.cpp
+++ b/src/model/GModelSpatialCube.cpp
@@ -435,9 +435,6 @@ void GModelSpatialCube::free_members(void)
  *
  * @param[in] filename File name.
  *
- * @exception GException::feature_not_implemented
- *            Map cube loading not yet implemented.
- *
  * @todo Implement method. 
  ***************************************************************************/
 void GModelSpatialCube::load_cube(const std::string& filename)
@@ -445,9 +442,6 @@ void GModelSpatialCube::load_cube(const std::string& filename)
     // Set filename
     m_filename = filename;
 
-    // Dump warning that method is not yet implemented
-    throw GException::feature_not_implemented(G_LOAD_CUBE);
-
     // Return
     return;
 }
diff --git a/src/model/GModelSpectralConst.cpp b/src/model/GModelSpectralConst.cpp
index e1184d5..27b76bf 100644
--- a/src/model/GModelSpectralConst.cpp
+++ b/src/model/GModelSpectralConst.cpp
@@ -12,7 +12,7 @@
  ***************************************************************************/
 /**
  * @file GModelSpectralConst.cpp
- * @brief GModelSpectralConst class implementation.
+ * @brief Constant spectral model class implementation
  * @author J. Knodlseder
  */
 
@@ -54,7 +54,7 @@ const GModelSpectralRegistry g_spectral_const_registry(&g_spectral_const_seed);
  ***************************************************************************/
 GModelSpectralConst::GModelSpectralConst(void) : GModelSpectral()
 {
-    // Initialise private members for clean destruction
+    // Initialise members
     init_members();
 
     // Return
@@ -63,16 +63,18 @@ GModelSpectralConst::GModelSpectralConst(void) : GModelSpectral()
 
 
 /***********************************************************************//**
- * @brief Constructor
+ * @brief XML constructor
  *
  * @param[in] xml XML element.
  *
- * Construct an arbitrary spectral function from a XML element.
+ * Creates instance of a constant spectral model by extracting information
+ * from an XML element. See GModelSpectralConst::read() for more information
+ * about the expected structure of the XML element.
  ***************************************************************************/
 GModelSpectralConst::GModelSpectralConst(const GXmlElement& xml) :
                      GModelSpectral()
 {
-    // Initialise private members for clean destruction
+    // Initialise members
     init_members();
 
     // Read information from XML element
@@ -86,12 +88,12 @@ GModelSpectralConst::GModelSpectralConst(const GXmlElement& xml) :
 /***********************************************************************//**
  * @brief Copy constructor
  *
- * @param[in] model Spectral component.
+ * @param[in] model Spectral constant model.
  ***************************************************************************/
 GModelSpectralConst::GModelSpectralConst(const GModelSpectralConst& model) :
                      GModelSpectral(model)
 {
-    // Initialise private members for clean destruction
+    // Initialise members
     init_members();
 
     // Copy members
@@ -124,7 +126,7 @@ GModelSpectralConst::~GModelSpectralConst(void)
 /***********************************************************************//**
  * @brief Assignment operator
  *
- * @param[in] model Spectral component.
+ * @param[in] model Spectral constant model.
  ***************************************************************************/
 GModelSpectralConst& GModelSpectralConst::operator= (const GModelSpectralConst& model)
 {
@@ -137,7 +139,7 @@ GModelSpectralConst& GModelSpectralConst::operator= (const GModelSpectralConst& 
         // Free members
         free_members();
 
-        // Initialise private members for clean destruction
+        // Initialise members
         init_members();
 
         // Copy members
@@ -193,7 +195,7 @@ GModelSpectralConst* GModelSpectralConst::clone(void) const
  * where
  * \f$norm\f$ is the normalization of the function.
  ***************************************************************************/
-double GModelSpectralConst::eval(const GEnergy& srcEng)
+double GModelSpectralConst::eval(const GEnergy& srcEng) const
 {
     // Compute function value
     double value = norm();
@@ -219,16 +221,16 @@ double GModelSpectralConst::eval(const GEnergy& srcEng)
  * The partial derivative of the normalization value is given by
  * \f[dI/dn_v=n_s\f]
  ***************************************************************************/
-double GModelSpectralConst::eval_gradients(const GEnergy& srcEng)
+double GModelSpectralConst::eval_gradients(const GEnergy& srcEng) const
 {
     // Compute function value
     double value = norm();
 
     // Compute partial derivatives of the parameter values
-    double g_norm = (m_norm.isfree())  ? m_norm.scale() : 0.0;
+    double g_norm = (m_norm.isfree()) ? m_norm.scale() : 0.0;
 
-    // Set gradients
-    m_norm.gradient(g_norm);
+    // Set gradients (circumvent const correctness)
+    ((GModelSpectralConst*)this)->m_norm.gradient(g_norm);
 
     // Return
     return value;
@@ -387,7 +389,7 @@ std::string GModelSpectralConst::print(void) const
     result.append("=== GModelSpectralConst ===\n");
     result.append(parformat("Number of parameters")+str(size()));
     for (int i = 0; i < size(); ++i)
-        result.append("\n"+m_par[i]->print());
+        result.append("\n"+m_pars[i]->print());
 
     // Return result
     return result;
@@ -405,17 +407,19 @@ std::string GModelSpectralConst::print(void) const
  ***************************************************************************/
 void GModelSpectralConst::init_members(void)
 {
-    // Initialise parameters
-    m_npars  = 1;
-    m_par[0] = &m_norm;
-
     // Initialise powerlaw normalisation
-    m_norm = GModelPar();
+    m_norm.clear();
     m_norm.name("Value");
     m_norm.value(1.0);
     m_norm.scale(1.0);
     m_norm.range(0.0, 1000.0);
     m_norm.free();
+    m_norm.gradient(0.0);
+    m_norm.hasgrad(true);
+
+    // Set parameter pointer(s)
+    m_pars.clear();
+    m_pars.push_back(&m_norm);
 
     // Return
     return;
@@ -425,13 +429,17 @@ void GModelSpectralConst::init_members(void)
 /***********************************************************************//**
  * @brief Copy class members
  *
- * @param[in] model Spectral model component.
+ * @param[in] model Spectral constant model.
  ***************************************************************************/
 void GModelSpectralConst::copy_members(const GModelSpectralConst& model)
 {
-    // Copy model parameters
+    // Copy members
     m_norm = model.m_norm;
 
+    // Set parameter pointer(s)
+    m_pars.clear();
+    m_pars.push_back(&m_norm);
+
     // Return
     return;
 }
diff --git a/src/model/GModelSpectralExpPlaw.cpp b/src/model/GModelSpectralExpPlaw.cpp
index 33fa1eb..6157b07 100644
--- a/src/model/GModelSpectralExpPlaw.cpp
+++ b/src/model/GModelSpectralExpPlaw.cpp
@@ -12,7 +12,7 @@
  ***************************************************************************/
 /**
  * @file GModelSpectralExpPlaw.cpp
- * @brief GModelSpectralExpPlaw class implementation.
+ * @brief Exponential cut off power law spectral class implementation
  * @author J. Knodlseder
  */
 
@@ -56,7 +56,7 @@ const GModelSpectralRegistry g_spectral_expplaw_registry(&g_spectral_expplaw_see
  ***************************************************************************/
 GModelSpectralExpPlaw::GModelSpectralExpPlaw(void) : GModelSpectral()
 {
-    // Initialise private members for clean destruction
+    // Initialise members
     init_members();
 
     // Return
@@ -69,7 +69,7 @@ GModelSpectralExpPlaw::GModelSpectralExpPlaw(void) : GModelSpectral()
  *
  * @param[in] norm Power law normalization.
  * @param[in] index Power law index.
- * @param[in] ecut Cut off energy.
+ * @param[in] ecut Cut off energy (in MeV).
  *
  * Construct an exponential cut off power law from a normalization value,
  * a spectral index and a cut off energy.
@@ -78,7 +78,7 @@ GModelSpectralExpPlaw::GModelSpectralExpPlaw(const double& norm,
                                              const double& index,
                                              const double& ecut) : GModelSpectral()
 {
-    // Initialise private members for clean destruction
+    // Initialise members
     init_members();
 
     // Set parameters
@@ -92,15 +92,18 @@ GModelSpectralExpPlaw::GModelSpectralExpPlaw(const double& norm,
 
 
 /***********************************************************************//**
- * @brief Constructor
+ * @brief XML constructor
  *
- * @param[in] xml XML element containing position information.
+ * @param[in] xml XML element.
  *
- * Construct a spectral power law from a XML element.
+ * Creates instance of an exponential cut off power law spectral model by
+ * extracting information from an XML element.
+ * See GModelSpectralExpPlaw::read() for more information about the expected
+ * structure of the XML element.
  ***************************************************************************/
 GModelSpectralExpPlaw::GModelSpectralExpPlaw(const GXmlElement& xml) : GModelSpectral()
 {
-    // Initialise private members for clean destruction
+    // Initialise members
     init_members();
 
     // Read information from XML element
@@ -114,12 +117,12 @@ GModelSpectralExpPlaw::GModelSpectralExpPlaw(const GXmlElement& xml) : GModelSpe
 /***********************************************************************//**
  * @brief Copy constructor
  *
- * @param[in] model Model from which the instance should be built.
+ * @param[in] model Exponential cut off power law model.
  ***************************************************************************/
 GModelSpectralExpPlaw::GModelSpectralExpPlaw(const GModelSpectralExpPlaw& model) :
                                                         GModelSpectral(model)
 {
-    // Initialise private members for clean destruction
+    // Initialise members
     init_members();
 
     // Copy members
@@ -152,7 +155,7 @@ GModelSpectralExpPlaw::~GModelSpectralExpPlaw(void)
 /***********************************************************************//**
  * @brief Assignment operator
  *
- * @param[in] model Model which should be assigned.
+ * @param[in] model Exponential cut off power law model.
  ***************************************************************************/
 GModelSpectralExpPlaw& GModelSpectralExpPlaw::operator= (const GModelSpectralExpPlaw& model)
 {
@@ -165,7 +168,7 @@ GModelSpectralExpPlaw& GModelSpectralExpPlaw::operator= (const GModelSpectralExp
         // Free members
         free_members();
 
-        // Initialise private members for clean destruction
+        // Initialise members
         init_members();
 
         // Copy members
@@ -229,7 +232,7 @@ GModelSpectralExpPlaw* GModelSpectralExpPlaw::clone(void) const
  * in the futur. Furthermore, the method expects that pivot!=0 and ecut!=0.
  * Otherwise Inf or NaN may result.
  ***************************************************************************/
-double GModelSpectralExpPlaw::eval(const GEnergy& srcEng)
+double GModelSpectralExpPlaw::eval(const GEnergy& srcEng) const
 {
     // Compute function value
     double energy = srcEng.MeV();
@@ -266,11 +269,11 @@ double GModelSpectralExpPlaw::eval(const GEnergy& srcEng)
  * \f[dI/dp_v=norm (E/pivot)^{index} \exp(-E/ecut) (-index) / p_v\f]
  *
  * @todo For the moment the pivot and the cut off energies are fixed to
- * units of MeV. This may not be ideal and should eventually be improved
- * in the futur. Furthermore, the method expects that pivot!=0 and ecut!=0.
- * Otherwise Inf or NaN may result.
+ *       units of MeV. This may not be ideal and should eventually be
+ *       improved in the futur. Furthermore, the method expects that
+ *       pivot!=0 and ecut!=0. Otherwise Inf or NaN may result.
  ***************************************************************************/
-double GModelSpectralExpPlaw::eval_gradients(const GEnergy& srcEng)
+double GModelSpectralExpPlaw::eval_gradients(const GEnergy& srcEng) const
 {
     // Compute function value
     double energy = srcEng.MeV();
@@ -285,11 +288,11 @@ double GModelSpectralExpPlaw::eval_gradients(const GEnergy& srcEng)
     double g_ecut  = (m_ecut.isfree())  ? value * e_cut/m_ecut.value() : 0.0;
     double g_pivot = (m_pivot.isfree()) ? -value * index() / m_pivot.value() : 0.0;
 
-    // Set gradients
-    m_norm.gradient(g_norm);
-    m_index.gradient(g_index);
-    m_ecut.gradient(g_ecut);
-    m_pivot.gradient(g_pivot);
+    // Set gradients (circumvent const correctness)
+    ((GModelSpectralExpPlaw*)this)->m_norm.gradient(g_norm);
+    ((GModelSpectralExpPlaw*)this)->m_index.gradient(g_index);
+    ((GModelSpectralExpPlaw*)this)->m_ecut.gradient(g_ecut);
+    ((GModelSpectralExpPlaw*)this)->m_pivot.gradient(g_pivot);
 
     // Return
     return value;
@@ -302,6 +305,9 @@ double GModelSpectralExpPlaw::eval_gradients(const GEnergy& srcEng)
  * @param[in] emin Minimum photon energy.
  * @param[in] emax Maximum photon energy.
  *
+ * @exception GException::feature_not_implemented
+ *            Method not yet implemented.
+ *
  * Computes
  * \f[\int_{E_{\rm min}}^{E_{\rm max}} I(E) dE\f]
  * where
@@ -309,7 +315,7 @@ double GModelSpectralExpPlaw::eval_gradients(const GEnergy& srcEng)
  * energy, respectively, and
  * \f$I(E)\f$ is the spectral model (units: ph/cm2/s/MeV).
  *
- * @todo Implement method
+ * @todo Implement method.
  ***************************************************************************/
 double GModelSpectralExpPlaw::flux(const GEnergy& emin, const GEnergy& emax) const
 {
@@ -328,7 +334,10 @@ double GModelSpectralExpPlaw::flux(const GEnergy& emin, const GEnergy& emax) con
  * @param[in] emax Maximum photon energy.
  * @param[in] ran Random number generator.
  *
- * @todo To be implemented
+ * @exception GException::feature_not_implemented
+ *            Method not yet implemented.
+ *
+ * @todo Implement method.
  ***************************************************************************/
 GEnergy GModelSpectralExpPlaw::mc(const GEnergy& emin, const GEnergy& emax,
                                   GRan& ran) const
@@ -347,8 +356,9 @@ GEnergy GModelSpectralExpPlaw::mc(const GEnergy& emin, const GEnergy& emax,
 /***********************************************************************//**
  * @brief Autoscale normalization
  *
- * Set the internal scaling of the normalization parameter so that the
- * parameter value equals 1.
+ * Based on the actual value of the m_norm parameter, set the scale of m_norm
+ * so that the value will be 1. If minimum and/or maximum value boundaries
+ * exist, the boundaries are also modified accordingly.
  ***************************************************************************/
 void GModelSpectralExpPlaw::autoscale(void)
 {
@@ -537,7 +547,7 @@ std::string GModelSpectralExpPlaw::print(void) const
     result.append("=== GModelSpectralExpPlaw ===\n");
     result.append(parformat("Number of parameters")+str(size()));
     for (int i = 0; i < size(); ++i)
-        result.append("\n"+m_par[i]->print());
+        result.append("\n"+m_pars[i]->print());
 
     // Return result
     return result;
@@ -555,43 +565,51 @@ std::string GModelSpectralExpPlaw::print(void) const
  ***************************************************************************/
 void GModelSpectralExpPlaw::init_members(void)
 {
-    // Initialise parameters
-    m_npars  = 4;
-    m_par[0] = &m_norm;
-    m_par[1] = &m_index;
-    m_par[2] = &m_ecut;
-    m_par[3] = &m_pivot;
-
     // Initialise powerlaw normalisation
-    m_norm = GModelPar();
+    m_norm.clear();
     m_norm.name("Prefactor");
     m_norm.unit("ph/cm2/s/MeV");
     m_norm.value(1.0);
     m_norm.scale(1.0);
     m_norm.free();
+    m_norm.gradient(0.0);
+    m_norm.hasgrad(true);
 
     // Initialise powerlaw index
-    m_index = GModelPar();
+    m_index.clear();
     m_index.name("Index");
     m_index.value(-2.0);
     m_index.scale(1.0);
     m_index.free();
+    m_index.gradient(0.0);
+    m_index.hasgrad(true);
 
     // Initialise cut off energy
-    m_ecut = GModelPar();
+    m_ecut.clear();
     m_ecut.name("Cutoff");
     m_ecut.unit("MeV");
     m_ecut.value(1000.0);
     m_ecut.scale(1.0);
     m_ecut.free();
+    m_ecut.gradient(0.0);
+    m_ecut.hasgrad(true);
 
     // Initialise pivot energy
-    m_pivot = GModelPar();
+    m_pivot.clear();
     m_pivot.name("PivotEnergy");
     m_pivot.unit("MeV");
     m_pivot.value(100.0);
     m_pivot.scale(1.0);
     m_pivot.fix();
+    m_pivot.gradient(0.0);
+    m_pivot.hasgrad(true);
+
+    // Set parameter pointer(s)
+    m_pars.clear();
+    m_pars.push_back(&m_norm);
+    m_pars.push_back(&m_index);
+    m_pars.push_back(&m_ecut);
+    m_pars.push_back(&m_pivot);
 
     // Return
     return;
@@ -601,17 +619,23 @@ void GModelSpectralExpPlaw::init_members(void)
 /***********************************************************************//**
  * @brief Copy class members
  *
- * @param[in] model GModelSpectralExpPlaw members which should be copied.
+ * @param[in] model Exponential cut off power law model.
  ***************************************************************************/
 void GModelSpectralExpPlaw::copy_members(const GModelSpectralExpPlaw& model)
 {
-    // Copy model parameters (we do not need to copy the rest since it is
-    // static)
+    // Copy members
     m_norm  = model.m_norm;
     m_index = model.m_index;
     m_ecut  = model.m_ecut;
     m_pivot = model.m_pivot;
 
+    // Set parameter pointer(s)
+    m_pars.clear();
+    m_pars.push_back(&m_norm);
+    m_pars.push_back(&m_index);
+    m_pars.push_back(&m_ecut);
+    m_pars.push_back(&m_pivot);
+
     // Return
     return;
 }
diff --git a/src/model/GModelSpectralFunc.cpp b/src/model/GModelSpectralFunc.cpp
index 2a14c34..5fc51c1 100644
--- a/src/model/GModelSpectralFunc.cpp
+++ b/src/model/GModelSpectralFunc.cpp
@@ -58,7 +58,7 @@ const GModelSpectralRegistry g_spectral_func_registry(&g_spectral_func_seed);
  ***************************************************************************/
 GModelSpectralFunc::GModelSpectralFunc(void) : GModelSpectral()
 {
-    // Initialise private members for clean destruction
+    // Initialise members
     init_members();
 
     // Return
@@ -67,17 +67,18 @@ GModelSpectralFunc::GModelSpectralFunc(void) : GModelSpectral()
 
 
 /***********************************************************************//**
- * @brief Constructor
+ * @brief File constructor
  *
  * @param[in] filename File name of nodes.
  *
- * Construct an arbitrary spectral function from a list of nodes that is
- * found in the specified file.
+ * Creates instance of a spectral function model from a list of nodes that is
+ * found in the specified file. See GModelSpectralFunc::load_nodes() for more
+ * information about the expected structure of the file.
  ***************************************************************************/
 GModelSpectralFunc::GModelSpectralFunc(const std::string& filename) :
                     GModelSpectral()
 {
-    // Initialise private members for clean destruction
+    // Initialise members
     init_members();
 
     // Load nodes
@@ -89,16 +90,18 @@ GModelSpectralFunc::GModelSpectralFunc(const std::string& filename) :
 
 
 /***********************************************************************//**
- * @brief Constructor
+ * @brief XML constructor
  *
  * @param[in] xml XML element.
  *
- * Construct an arbitrary spectral function from a XML element.
+ * Creates instance of a spectral function model by extracting information
+ * from an XML element. See GModelSpectralFunc::read() for more information
+ * about the expected structure of the XML element.
  ***************************************************************************/
 GModelSpectralFunc::GModelSpectralFunc(const GXmlElement& xml) :
                     GModelSpectral()
 {
-    // Initialise private members for clean destruction
+    // Initialise members
     init_members();
 
     // Read information from XML element
@@ -112,12 +115,12 @@ GModelSpectralFunc::GModelSpectralFunc(const GXmlElement& xml) :
 /***********************************************************************//**
  * @brief Copy constructor
  *
- * @param[in] model Spectral component.
+ * @param[in] model Spectral function model.
  ***************************************************************************/
 GModelSpectralFunc::GModelSpectralFunc(const GModelSpectralFunc& model) :
                     GModelSpectral(model)
 {
-    // Initialise private members for clean destruction
+    // Initialise members
     init_members();
 
     // Copy members
@@ -150,7 +153,7 @@ GModelSpectralFunc::~GModelSpectralFunc(void)
 /***********************************************************************//**
  * @brief Assignment operator
  *
- * @param[in] model Spectral component.
+ * @param[in] model Spectral function model.
  ***************************************************************************/
 GModelSpectralFunc& GModelSpectralFunc::operator= (const GModelSpectralFunc& model)
 {
@@ -163,7 +166,7 @@ GModelSpectralFunc& GModelSpectralFunc::operator= (const GModelSpectralFunc& mod
         // Free members
         free_members();
 
-        // Initialise private members for clean destruction
+        // Initialise members
         init_members();
 
         // Copy members
@@ -221,7 +224,7 @@ GModelSpectralFunc* GModelSpectralFunc::clone(void) const
  * Note that the node energies are stored as log10 of energy in units of
  * MeV.
  ***************************************************************************/
-double GModelSpectralFunc::eval(const GEnergy& srcEng)
+double GModelSpectralFunc::eval(const GEnergy& srcEng) const
 {
     // Interpolate function
     double func = m_nodes.interpolate(srcEng.log10MeV(), m_values);
@@ -250,7 +253,7 @@ double GModelSpectralFunc::eval(const GEnergy& srcEng)
  * The partial derivative of the normalization value is given by
  * \f[dI/dn_v=n_s f(E)\f]
  ***************************************************************************/
-double GModelSpectralFunc::eval_gradients(const GEnergy& srcEng)
+double GModelSpectralFunc::eval_gradients(const GEnergy& srcEng) const
 {
     // Interpolate function
     double func = m_nodes.interpolate(srcEng.log10MeV(), m_values);
@@ -261,8 +264,8 @@ double GModelSpectralFunc::eval_gradients(const GEnergy& srcEng)
     // Compute partial derivatives of the parameter values
     double g_norm  = (m_norm.isfree())  ? m_norm.scale() * func : 0.0;
 
-    // Set gradients
-    m_norm.gradient(g_norm);
+    // Set gradients (circumvent const correctness)
+    ((GModelSpectralFunc*)this)->m_norm.gradient(g_norm);
 
     // Return
     return value;
@@ -275,6 +278,9 @@ double GModelSpectralFunc::eval_gradients(const GEnergy& srcEng)
  * @param[in] emin Minimum photon energy.
  * @param[in] emax Maximum photon energy.
  *
+ * @exception GException::feature_not_implemented
+ *            Method not yet implemented
+ *
  * Computes
  * \f[\int_{E_{\rm min}}^{E_{\rm max}} I(E) dE\f]
  * where
@@ -301,6 +307,9 @@ double GModelSpectralFunc::flux(const GEnergy& emin, const GEnergy& emax) const
  * @param[in] emax Maximum photon energy.
  * @param[in] ran Random number generator.
  *
+ * @exception GException::feature_not_implemented
+ *            Method not yet implemented
+ *
  * @todo To be implemented
  ***************************************************************************/
 GEnergy GModelSpectralFunc::mc(const GEnergy& emin, const GEnergy& emax,
@@ -426,7 +435,7 @@ std::string GModelSpectralFunc::print(void) const
     result.append(parformat("Number of nodes")+str(m_nodes.size()));
     result.append(parformat("Number of parameters")+str(size()));
     for (int i = 0; i < size(); ++i)
-        result.append("\n"+m_par[i]->print());
+        result.append("\n"+m_pars[i]->print());
 
     // Return result
     return result;
@@ -444,18 +453,20 @@ std::string GModelSpectralFunc::print(void) const
  ***************************************************************************/
 void GModelSpectralFunc::init_members(void)
 {
-    // Initialise parameters
-    m_npars  = 1;
-    m_par[0] = &m_norm;
-
     // Initialise powerlaw normalisation
-    m_norm = GModelPar();
+    m_norm.clear();
     m_norm.name("Normalization");
     m_norm.min(0.0);
     m_norm.max(1000.0);
     m_norm.value(1.0);
     m_norm.scale(1.0);
     m_norm.free();
+    m_norm.gradient(0.0);
+    m_norm.hasgrad(true);
+
+    // Set parameter pointer(s)
+    m_pars.clear();
+    m_pars.push_back(&m_norm);
 
     // Initialise other members
     m_nodes.clear();
@@ -470,17 +481,20 @@ void GModelSpectralFunc::init_members(void)
 /***********************************************************************//**
  * @brief Copy class members
  *
- * @param[in] model Spectral model component.
+ * @param[in] model Spectral function model.
  ***************************************************************************/
 void GModelSpectralFunc::copy_members(const GModelSpectralFunc& model)
 {
-    // Copy model parameters (we do not need to copy the rest since it is
-    // static)
+    // Copy members
     m_norm     = model.m_norm;
     m_nodes    = model.m_nodes;
     m_values   = model.m_values;
     m_filename = model.m_filename;
 
+    // Set parameter pointer(s)
+    m_pars.clear();
+    m_pars.push_back(&m_norm);
+
     // Return
     return;
 }
diff --git a/src/model/GModelSpectralPlaw.cpp b/src/model/GModelSpectralPlaw.cpp
index 9855693..eaf59f9 100644
--- a/src/model/GModelSpectralPlaw.cpp
+++ b/src/model/GModelSpectralPlaw.cpp
@@ -12,7 +12,7 @@
  ***************************************************************************/
 /**
  * @file GModelSpectralPlaw.cpp
- * @brief GModelSpectralPlaw class implementation.
+ * @brief Power law spectral model class implementation
  * @author J. Knodlseder
  */
 
@@ -74,7 +74,7 @@ GModelSpectralPlaw::GModelSpectralPlaw(void) : GModelSpectral()
 GModelSpectralPlaw::GModelSpectralPlaw(const double& norm,
                                        const double& index) : GModelSpectral()
 {
-    // Initialise private members for clean destruction
+    // Initialise members
     init_members();
 
     // Set parameters
@@ -87,15 +87,17 @@ GModelSpectralPlaw::GModelSpectralPlaw(const double& norm,
 
 
 /***********************************************************************//**
- * @brief Constructor
+ * @brief XML constructor
  *
  * @param[in] xml XML element containing position information.
  *
- * Construct a spectral power law from a XML element.
+ * Creates instance of a power law spectral model by extracting information
+ * from an XML element. See GModelSpectralPlaw::read() for more information
+ * about the expected structure of the XML element.
  ***************************************************************************/
 GModelSpectralPlaw::GModelSpectralPlaw(const GXmlElement& xml) : GModelSpectral()
 {
-    // Initialise private members for clean destruction
+    // Initialise members
     init_members();
 
     // Read information from XML element
@@ -109,12 +111,12 @@ GModelSpectralPlaw::GModelSpectralPlaw(const GXmlElement& xml) : GModelSpectral(
 /***********************************************************************//**
  * @brief Copy constructor
  *
- * @param[in] model Model from which the instance should be built.
+ * @param[in] model Spectral power law model.
  ***************************************************************************/
 GModelSpectralPlaw::GModelSpectralPlaw(const GModelSpectralPlaw& model) :
                                                         GModelSpectral(model)
 {
-    // Initialise private members for clean destruction
+    // Initialise members
     init_members();
 
     // Copy members
@@ -147,9 +149,9 @@ GModelSpectralPlaw::~GModelSpectralPlaw(void)
 /***********************************************************************//**
  * @brief Assignment operator
  *
- * @param[in] model Model which should be assigned.
+ * @param[in] model Spectral power law model.
  ***************************************************************************/
-GModelSpectralPlaw& GModelSpectralPlaw::operator= (const GModelSpectralPlaw& model)
+GModelSpectralPlaw& GModelSpectralPlaw::operator=(const GModelSpectralPlaw& model)
 {
     // Execute only if object is not identical
     if (this != &model) {
@@ -160,7 +162,7 @@ GModelSpectralPlaw& GModelSpectralPlaw::operator= (const GModelSpectralPlaw& mod
         // Free members
         free_members();
 
-        // Initialise private members for clean destruction
+        // Initialise members
         init_members();
 
         // Copy members
@@ -223,7 +225,7 @@ GModelSpectralPlaw* GModelSpectralPlaw::clone(void) const
  * Furthermore, the method expects that energy!=0. Otherwise Inf or NaN
  * may result.
  ***************************************************************************/
-double GModelSpectralPlaw::eval(const GEnergy& srcEng)
+double GModelSpectralPlaw::eval(const GEnergy& srcEng) const
 {
     // Compute function value
     double energy = srcEng.MeV() / pivot();
@@ -260,7 +262,7 @@ double GModelSpectralPlaw::eval(const GEnergy& srcEng)
  * Furthermore, the method expects that energy!=0. Otherwise Inf or NaN
  * may result.
  ***************************************************************************/
-double GModelSpectralPlaw::eval_gradients(const GEnergy& srcEng)
+double GModelSpectralPlaw::eval_gradients(const GEnergy& srcEng) const
 {
     // Compute function value
     double energy = srcEng.MeV() / pivot();
@@ -272,10 +274,10 @@ double GModelSpectralPlaw::eval_gradients(const GEnergy& srcEng)
     double g_index = (m_index.isfree()) ? value * m_index.scale() * std::log(energy) : 0.0;
     double g_pivot = (m_pivot.isfree()) ? -value * index() / m_pivot.value() : 0.0;
 
-    // Set gradients
-    m_norm.gradient(g_norm);
-    m_index.gradient(g_index);
-    m_pivot.gradient(g_pivot);
+    // Set gradients (circumvent const correctness)
+    ((GModelSpectralPlaw*)this)->m_norm.gradient(g_norm);
+    ((GModelSpectralPlaw*)this)->m_index.gradient(g_index);
+    ((GModelSpectralPlaw*)this)->m_pivot.gradient(g_pivot);
 
     // Return
     return value;
@@ -301,7 +303,8 @@ double GModelSpectralPlaw::flux(const GEnergy& emin, const GEnergy& emax) const
     double flux = norm() * pow(pivot(), -index());
     if (index() != -1.0) {
         double exponent = index() + 1.0;
-        flux *= (std::pow(emax.MeV(), exponent) - std::pow(emin.MeV(), exponent)) / exponent;
+        flux *= (std::pow(emax.MeV(), exponent)-std::pow(emin.MeV(), exponent)) /
+                exponent;
     }
     else
         flux *= (std::log(emax.MeV()) - std::log(emin.MeV()));
@@ -355,8 +358,9 @@ GEnergy GModelSpectralPlaw::mc(const GEnergy& emin, const GEnergy& emax,
 /***********************************************************************//**
  * @brief Autoscale normalization
  *
- * Set the internal scaling of the normalization parameter so that the
- * parameter value equals 1.
+ * Based on the actual value of the m_norm parameter, set the scale of m_norm
+ * so that the value will be 1. If minimum and/or maximum value boundaries
+ * exist, the boundaries are also modified accordingly.
  ***************************************************************************/
 void GModelSpectralPlaw::autoscale(void)
 {
@@ -530,7 +534,7 @@ std::string GModelSpectralPlaw::print(void) const
     result.append("=== GModelSpectralPlaw ===\n");
     result.append(parformat("Number of parameters")+str(size()));
     for (int i = 0; i < size(); ++i)
-        result.append("\n"+m_par[i]->print());
+        result.append("\n"+m_pars[i]->print());
 
     // Return result
     return result;
@@ -548,34 +552,40 @@ std::string GModelSpectralPlaw::print(void) const
  ***************************************************************************/
 void GModelSpectralPlaw::init_members(void)
 {
-    // Initialise parameters
-    m_npars  = 3;
-    m_par[0] = &m_norm;
-    m_par[1] = &m_index;
-    m_par[2] = &m_pivot;
-
     // Initialise powerlaw normalisation
-    m_norm = GModelPar();
+    m_norm.clear();
     m_norm.name("Prefactor");
     m_norm.unit("ph/cm2/s/MeV");
     m_norm.value(1.0);
     m_norm.scale(1.0);
     m_norm.free();
+    m_norm.gradient(0.0);
+    m_norm.hasgrad(true);
 
     // Initialise powerlaw index
-    m_index = GModelPar();
+    m_index.clear();
     m_index.name("Index");
     m_index.value(-2.0);
     m_index.scale(1.0);
     m_index.free();
+    m_index.gradient(0.0);
+    m_index.hasgrad(true);
 
     // Initialise pivot energy
-    m_pivot = GModelPar();
+    m_pivot.clear();
     m_pivot.name("PivotEnergy");
     m_pivot.unit("MeV");
     m_pivot.value(100.0);
     m_pivot.scale(1.0);
     m_pivot.fix();
+    m_pivot.gradient(0.0);
+    m_pivot.hasgrad(true);
+
+    // Set parameter pointer(s)
+    m_pars.clear();
+    m_pars.push_back(&m_norm);
+    m_pars.push_back(&m_index);
+    m_pars.push_back(&m_pivot);
 
     // Return
     return;
@@ -589,12 +599,17 @@ void GModelSpectralPlaw::init_members(void)
  ***************************************************************************/
 void GModelSpectralPlaw::copy_members(const GModelSpectralPlaw& model)
 {
-    // Copy model parameters (we do not need to copy the rest since it is
-    // static)
+    // Copy members
     m_norm  = model.m_norm;
     m_index = model.m_index;
     m_pivot = model.m_pivot;
 
+    // Set parameter pointer(s)
+    m_pars.clear();
+    m_pars.push_back(&m_norm);
+    m_pars.push_back(&m_index);
+    m_pars.push_back(&m_pivot);
+
     // Return
     return;
 }
diff --git a/src/model/GModelSpectralPlaw2.cpp b/src/model/GModelSpectralPlaw2.cpp
index 4ec23d2..695f0ee 100644
--- a/src/model/GModelSpectralPlaw2.cpp
+++ b/src/model/GModelSpectralPlaw2.cpp
@@ -12,7 +12,7 @@
  ***************************************************************************/
 /**
  * @file GModelSpectralPlaw2.cpp
- * @brief GModelSpectralPlaw2 class implementation.
+ * @brief Flux normalized power law spectral model class implementation
  * @author J. Knodlseder
  */
 
@@ -55,7 +55,7 @@ const GModelSpectralRegistry g_spectral_plaw2_registry(&g_spectral_plaw2_seed);
  ***************************************************************************/
 GModelSpectralPlaw2::GModelSpectralPlaw2(void) : GModelSpectral()
 {
-    // Initialise private members for clean destruction
+    // Initialise members
     init_members();
 
     // Return
@@ -75,7 +75,7 @@ GModelSpectralPlaw2::GModelSpectralPlaw2(const double& integral,
                                          const double& index) :
                                          GModelSpectral()
 {
-    // Initialise private members for clean destruction
+    // Initialise members
     init_members();
 
     // Set parameters
@@ -88,16 +88,18 @@ GModelSpectralPlaw2::GModelSpectralPlaw2(const double& integral,
 
 
 /***********************************************************************//**
- * @brief Constructor
+ * @brief XML constructor
  *
  * @param[in] xml XML element.
  *
- * Construct a spectral power law from a XML element.
+ * Creates instance of a flux normalized power law spectral model by
+ * extracting information from an XML element. See GModelSpectralPlaw2::read()
+ * for more information about the expected structure of the XML element.
  ***************************************************************************/
 GModelSpectralPlaw2::GModelSpectralPlaw2(const GXmlElement& xml) :
                      GModelSpectral()
 {
-    // Initialise private members for clean destruction
+    // Initialise members
     init_members();
 
     // Read information from XML element
@@ -116,7 +118,7 @@ GModelSpectralPlaw2::GModelSpectralPlaw2(const GXmlElement& xml) :
 GModelSpectralPlaw2::GModelSpectralPlaw2(const GModelSpectralPlaw2& model) :
                      GModelSpectral(model)
 {
-    // Initialise private members for clean destruction
+    // Initialise members
     init_members();
 
     // Copy members
@@ -162,7 +164,7 @@ GModelSpectralPlaw2& GModelSpectralPlaw2::operator= (const GModelSpectralPlaw2& 
         // Free members
         free_members();
 
-        // Initialise private members for clean destruction
+        // Initialise members
         init_members();
 
         // Copy members
@@ -223,10 +225,10 @@ GModelSpectralPlaw2* GModelSpectralPlaw2::clone(void) const
  * \f$emin\f$ is the lower energy limit, and
  * \f$emax\f$ is the upper energy limit.
  ***************************************************************************/
-double GModelSpectralPlaw2::eval(const GEnergy& srcEng)
+double GModelSpectralPlaw2::eval(const GEnergy& srcEng) const
 {
     // Update precomputed values
-    update(srcEng);
+    ((GModelSpectralPlaw2*)this)->update(srcEng);
 
     // Compute function value
     double value = integral() * m_norm * m_power;
@@ -265,17 +267,17 @@ double GModelSpectralPlaw2::eval(const GEnergy& srcEng)
  * \f[dI/df_v=I(E) / f_v\f]
  * \f[dI/di_v=I(E) \log(E) i_s\f]
  * for \f$index = -1\f$.
- * The partial derivatives for the energy boundaries are fixed to 0 as it is
- * not expected the fit these parameters.
+ *
+ * No partial derivatives are supported for the energy boundaries.
  ***************************************************************************/
-double GModelSpectralPlaw2::eval_gradients(const GEnergy& srcEng)
+double GModelSpectralPlaw2::eval_gradients(const GEnergy& srcEng) const
 {
     // Initialise gradients
     double g_integral = 0.0;
     double g_index    = 0.0;
     
     // Update precomputed values
-    update(srcEng);
+    ((GModelSpectralPlaw2*)this)->update(srcEng);
 
     // Compute function value
     double value = integral() * m_norm * m_power;
@@ -289,10 +291,8 @@ double GModelSpectralPlaw2::eval_gradients(const GEnergy& srcEng)
         g_index = value * (m_g_norm + ln10*srcEng.log10MeV()) * m_index.scale();
 
     // Set gradients
-    m_integral.gradient(g_integral);
-    m_index.gradient(g_index);
-    m_emin.gradient(0.0);
-    m_emax.gradient(0.0);
+    ((GModelSpectralPlaw2*)this)->m_integral.gradient(g_integral);
+    ((GModelSpectralPlaw2*)this)->m_index.gradient(g_index);
 
     // Return
     return value;
@@ -526,7 +526,7 @@ std::string GModelSpectralPlaw2::print(void) const
     result.append("=== GModelSpectralPlaw2 ===\n");
     result.append(parformat("Number of parameters")+str(size()));
     for (int i = 0; i < size(); ++i)
-        result.append("\n"+m_par[i]->print());
+        result.append("\n"+m_pars[i]->print());
 
     // Return result
     return result;
@@ -544,47 +544,55 @@ std::string GModelSpectralPlaw2::print(void) const
  ***************************************************************************/
 void GModelSpectralPlaw2::init_members(void)
 {
-    // Initialise parameters
-    m_npars  = 4;
-    m_par[0] = &m_integral;
-    m_par[1] = &m_index;
-    m_par[2] = &m_emin;
-    m_par[3] = &m_emax;
-
     // Initialise integral flux
-    m_integral = GModelPar();
+    m_integral.clear();
     m_integral.name("Integral");
     m_integral.unit("ph/cm2/s");
     m_integral.value(1.0);
     m_integral.scale(1.0);
     m_integral.range(0.0, 10.0);
     m_integral.free();
+    m_integral.gradient(0.0);
+    m_integral.hasgrad(true);
 
     // Initialise powerlaw index
-    m_index = GModelPar();
+    m_index.clear();
     m_index.name("Index");
     m_index.value(-2.0);
     m_index.scale(1.0);
     m_index.range(-5.0, -0.1);
     m_index.free();
+    m_index.gradient(0.0);
+    m_index.hasgrad(true);
 
     // Initialise lower limit
-    m_emin = GModelPar();
+    m_emin.clear();
     m_emin.name("LowerLimit");
     m_emin.unit("MeV");
     m_emin.value(100.0);
     m_emin.scale(1.0);
     m_emin.range(0.001, 1.0e15);
     m_emin.fix();
+    m_emin.gradient(0.0);
+    m_emin.hasgrad(false);
 
     // Initialise upper limit
-    m_emax = GModelPar();
+    m_emax.clear();
     m_emax.name("UpperLimit");
     m_emax.unit("MeV");
     m_emax.value(500000.0);
     m_emax.scale(1.0);
     m_emax.range(0.001, 1.0e15);
     m_emax.fix();
+    m_emax.gradient(0.0);
+    m_emax.hasgrad(false);
+
+    // Set parameter pointer(s)
+    m_pars.clear();
+    m_pars.push_back(&m_integral);
+    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;
@@ -614,12 +622,19 @@ void GModelSpectralPlaw2::init_members(void)
  ***************************************************************************/
 void GModelSpectralPlaw2::copy_members(const GModelSpectralPlaw2& model)
 {
-    // Copy model parameters (no need to copy others since they are static)
+    // Copy members
     m_integral = model.m_integral;
     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_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;
diff --git a/src/model/GModelTemporalConst.cpp b/src/model/GModelTemporalConst.cpp
index b2ebf85..7937683 100644
--- a/src/model/GModelTemporalConst.cpp
+++ b/src/model/GModelTemporalConst.cpp
@@ -12,7 +12,7 @@
  ***************************************************************************/
 /**
  * @file GModelTemporalConst.cpp
- * @brief GModelTemporalConst class implementation.
+ * @brief Constant temporal model class interface implementation
  * @author J. Knodlseder
  */
 
@@ -51,7 +51,7 @@ const GModelTemporalRegistry g_temporal_const_registry(&g_temporal_const_seed);
  ***************************************************************************/
 GModelTemporalConst::GModelTemporalConst(void) : GModelTemporal()
 {
-    // Initialise private members for clean destruction
+    // Initialise members
     init_members();
 
     // Return
@@ -62,12 +62,12 @@ GModelTemporalConst::GModelTemporalConst(void) : GModelTemporal()
 /***********************************************************************//**
  * @brief Copy constructor
  *
- * @param[in] model Model from which the instance should be built.
+ * @param[in] model Constant temporal model
  ***************************************************************************/
 GModelTemporalConst::GModelTemporalConst(const GModelTemporalConst& model) :
     GModelTemporal(model)
 {
-    // Initialise private members for clean destruction
+    // Initialise members
     init_members();
 
     // Copy members
@@ -100,7 +100,7 @@ GModelTemporalConst::~GModelTemporalConst(void)
 /***********************************************************************//**
  * @brief Assignment operator
  *
- * @param[in] model Model which should be assigned.
+ * @param[in] model Constant temporal model
  ***************************************************************************/
 GModelTemporalConst& GModelTemporalConst::operator= (const GModelTemporalConst& model)
 {
@@ -113,7 +113,7 @@ GModelTemporalConst& GModelTemporalConst::operator= (const GModelTemporalConst& 
         // Free members
         free_members();
 
-        // Initialise private members for clean destruction
+        // Initialise members
         init_members();
 
         // Copy members
@@ -167,10 +167,13 @@ GModelTemporalConst* GModelTemporalConst::clone(void) const
  * This method implements the temporal component of a constant model. It
  * simply returns the normalization constant (which typically is set to 1).
  ***************************************************************************/
-double GModelTemporalConst::eval(const GTime& srcTime)
+double GModelTemporalConst::eval(const GTime& srcTime) const
 {
+    // Compute function value
+    double value = norm();
+
     // Return
-    return (norm());
+    return value;
 }
 
 
@@ -179,17 +182,30 @@ double GModelTemporalConst::eval(const GTime& srcTime)
  *
  * @param[in] srcTime True photon arrival time (not used).
  *
- * This method implements the temporal component of a constant model. It
- * simply returns the normalization constant (which typically is set to 1)
- * and a gradient of 0 (as it is not expected to fit this parameter).
+ * The temporal model is defined as
+ * \f[I(E)=norm\f]
+ * where
+ * \f$norm=n_s n_v\f$ is the normalization of the function.
+ * Note that the normalization is factorised into a scaling factor and a
+ * value and that the method is expected to return the gradient with respect
+ * to the parameter value \f$n_v\f$.
+ *
+ * The partial derivative of the normalization value is given by
+ * \f[dI/dn_v=n_s\f]
  ***************************************************************************/
-double GModelTemporalConst::eval_gradients(const GTime& srcTime)
+double GModelTemporalConst::eval_gradients(const GTime& srcTime) const
 {
-    // Set gradient to 0
-    m_norm.gradient(0.0);
+    // Compute function value
+    double value = norm();
+
+    // Compute partial derivatives of the parameter values
+    double g_norm = (m_norm.isfree()) ? m_norm.scale() : 0.0;
+
+    // Set gradients (circumvent const correctness)
+    ((GModelTemporalConst*)this)->m_norm.gradient(g_norm);
 
     // Return
-    return (norm());
+    return value;
 }
 
 
@@ -205,7 +221,7 @@ double GModelTemporalConst::eval_gradients(const GTime& srcTime)
  * event rate that is specified by the rate parameter.
  ***************************************************************************/
 GTimes GModelTemporalConst::mc(const double& rate, const GTime&  tmin,
-                               const GTime&  tmax, GRan& ran)
+                               const GTime&  tmax, GRan& ran) const
 {
     // Allocates empty vector of times
     GTimes times;
@@ -240,9 +256,9 @@ GTimes GModelTemporalConst::mc(const double& rate, const GTime&  tmin,
 /***********************************************************************//**
  * @brief Read model from XML element
  *
- * @param[in] xml XML element containing power law model information.
+ * @param[in] xml XML element.
  *
- * @todo Not yet implemented.
+ * @todo To be implemented.
  ***************************************************************************/
 void GModelTemporalConst::read(const GXmlElement& xml)
 {
@@ -254,9 +270,9 @@ void GModelTemporalConst::read(const GXmlElement& xml)
 /***********************************************************************//**
  * @brief Write model into XML element
  *
- * @param[in] xml XML element into which model information is written.
+ * @param[in] xml XML element.
  *
- * @todo Not yet implemented.
+ * @todo To be implemented.
  ***************************************************************************/
 void GModelTemporalConst::write(GXmlElement& xml) const
 {
@@ -277,7 +293,7 @@ std::string GModelTemporalConst::print(void) const
     result.append("=== GModelTemporalConst ===\n");
     result.append(parformat("Number of parameters")+str(size()));
     for (int i = 0; i < size(); ++i)
-        result.append("\n"+m_par[i]->print());
+        result.append("\n"+m_pars[i]->print());
 
     // Return result
     return result;
@@ -295,16 +311,18 @@ std::string GModelTemporalConst::print(void) const
  ***************************************************************************/
 void GModelTemporalConst::init_members(void)
 {
-    // Initialise parameters
-    m_npars  = 1;
-    m_par[0] = &m_norm;
-
     // Initialise normalisation parameter
-    m_norm = GModelPar();
+    m_norm.clear();
     m_norm.name("Constant");
     m_norm.unit("(relative value)");
     m_norm.value(1.0);
     m_norm.fix();
+    m_norm.gradient(0.0);
+    m_norm.hasgrad(true);
+
+    // Set parameter pointer(s)
+    m_pars.clear();
+    m_pars.push_back(&m_norm);
 
     // Return
     return;
@@ -314,13 +332,16 @@ void GModelTemporalConst::init_members(void)
 /***********************************************************************//**
  * @brief Copy class members
  *
- * @param[in] model GModelTemporalConst members which should be copied.
+ * @param[in] model Constant temporal model
  ***************************************************************************/
 void GModelTemporalConst::copy_members(const GModelTemporalConst& model)
 {
-    // Copy model parameters (we do not need to copy the rest since it is
-    // static)
-    m_norm  = model.m_norm;
+    // Copy members
+    m_norm = model.m_norm;
+
+    // Set parameter pointer(s)
+    m_pars.clear();
+    m_pars.push_back(&m_norm);
 
     // Return
     return;