Merge branch 'devel' into hess

Jürgen Knödlseder
2 parents 8624151e 34fb46eb
Showing 24 changed files with 912 additions and 299 deletions
configure.ac
include/GEbounds.hpp
inst/cta/include/GCTAEventBin.hpp
inst/cta/include/GCTAException.hpp
inst/cta/src/GCTAEventBin.cpp
inst/cta/src/GCTAException.cpp
inst/lat/include/GLATEventBin.hpp
inst/lat/include/GLATException.hpp
inst/lat/src/GLATAeff.cpp
inst/lat/src/GLATEventBin.cpp
inst/lat/src/GLATEventCube.cpp
inst/lat/src/GLATException.cpp
inst/lat/src/GLATResponseTable.cpp
pyext/GEbounds.i
src/fits/GFits.cpp
src/fits/GFitsHDU.cpp
src/fits/GFitsImage.cpp
src/model/GModels.cpp
src/obs/GEbounds.cpp
src/obs/GEnergy.cpp
@@ -659,8 +659,8 @@ if test &quot;x$openmp_ok&quot; = &quot;xyes&quot; ; then
   # dynamically
   if test "$OPENMP_CXXFLAGS" = "-fopenmp" -a "$ac_sys_system" = "Linux"; then
     AC_MSG_CHECKING(whether we can dynamically load libgomp)
-    SAVED_LDFLAGS=$LDFLAGS
-    LDFLAGS="$LDFLAGS -ldl"
+    SAVED_LIBS=$LIBS
+    LIBS="-ldl"
     AC_TRY_RUN(AC_LANG_PROGRAM([#include <dlfcn.h>],
                                [if (dlopen("libgomp.so.1", RTLD_NOW) == 0) return 1;]),
                [AC_MSG_RESULT(yes)],
@@ -668,7 +668,7 @@ if test &quot;x$openmp_ok&quot; = &quot;xyes&quot; ; then
                 AC_MSG_RESULT(no)
                 AC_MSG_WARN([OpenMP detected but shared libgomp.so.1 library cannot be dlopen()ed.])
                 AC_MSG_WARN([If you use a gcc 4.2.x compiler, this is a known feature. Try to use gcc >= 4.3.])])
-    LDFLAGS=$SAVED_LDFLAGS
+    LIBS=$SAVED_LIBS
   fi
  
   # Special Darwin test to see whether OpenMP compiles without adding the
@@ -76,7 +76,7 @@ public:
     GEnergy     emax(int inx) const;
     GEnergy     emean(int inx) const;
     GEnergy     elogmean(int inx) const;
-    bool        isin(const GEnergy& eng) const;
+    bool        contains(const GEnergy& eng) const;
     std::string print(void) const;
  
  
 /***************************************************************************
- *                 GCTAEventBin.hpp  -  CTA event bin class                *
+ *                  GCTAEventBin.hpp - CTA event bin class                 *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010 by Jurgen Knodlseder                                *
+ *  copyright (C) 2010-2012 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -20,8 +20,8 @@
  ***************************************************************************/
 /**
  * @file GCTAEventBin.hpp
- * @brief GCTAEventBin class interface definition.
- * @author J. Knodlseder
+ * @brief CTA event bin class interface definition
+ * @author Juergen Knoedlseder
  */
  
 #ifndef GCTAEVENTBIN_HPP
@@ -67,18 +67,18 @@ public:
     virtual void               clear(void);
     virtual GCTAEventBin*      clone(void) const;
     virtual double             size(void) const;
-    virtual const GCTAInstDir& dir(void) const { return *m_dir; }
-    virtual const GEnergy&     energy(void) const { return *m_energy; }
-    virtual const GTime&       time(void) const { return *m_time; }
-    virtual double             counts(void) const { return *m_counts; }
+    virtual const GCTAInstDir& dir(void) const;
+    virtual const GEnergy&     energy(void) const;
+    virtual const GTime&       time(void) const;
+    virtual double             counts(void) const;
     virtual double             error(void) const;
-    virtual void               counts(const double& counts) { *m_counts=counts; }
+    virtual void               counts(const double& counts);
     virtual std::string        print(void) const;
  
     // Other methods
-    const double&  omega(void) const { return *m_omega; }
-    const GEnergy& ewidth(void) const { return *m_ewidth; }
-    const double&  ontime(void) const { return *m_ontime; }
+    const double&  omega(void) const;
+    const GEnergy& ewidth(void) const;
+    const double&  ontime(void) const;
  
 protected:
     // Protected methods
 /***************************************************************************
  *                GCTAException.hpp  - CTA exception handler               *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2011 by Jurgen Knodlseder                           *
+ *  copyright (C) 2010-2012 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -21,7 +21,7 @@
 /**
  * @file GCTAException.hpp
  * @brief CTA exception handler interface definition.
- * @author J. Knodlseder
+ * @author Juergen Knoedlseder
  */
  
 #ifndef GCTAEXCEPTION_HPP
@@ -65,6 +65,12 @@ public:
                         std::string message = "");
     };
  
+    // Event bin exceptions
+    class no_member : public GCTAExceptionHandler {
+    public:
+        no_member(const std::string& origin, const std::string& message);
+    };
+
     // Observation exceptions
     class no_response : public GCTAExceptionHandler {
     public:
 /***************************************************************************
- *                GCTAEventBin.cpp  -  CTA event bin class                 *
+ *                 GCTAEventBin.cpp - CTA event bin class                  *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010 by Jurgen Knodlseder                                *
+ *  copyright (C) 2010-2012 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -20,8 +20,8 @@
  ***************************************************************************/
 /**
  * @file GCTAEventBin.cpp
- * @brief GCTAEventBin class implementation.
- * @author J. Knodlseder
+ * @brief CTA event bin class implementation
+ * @author Juergen Knoedlseder
  */
  
 /* __ Includes ___________________________________________________________ */
@@ -31,9 +31,18 @@
 #include <string>
 #include <cmath>
 #include "GCTAEventBin.hpp"
+#include "GCTAException.hpp"
 #include "GTools.hpp"
  
 /* __ Method name definitions ____________________________________________ */
+#define G_DIR                                           "GCTAEventBin::dir()"
+#define G_ENERGY                                     "GCTAEventBin::energy()"
+#define G_TIME                                         "GCTAEventBin::time()"
+#define G_COUNTS_GET                                 "GCTAEventBin::counts()"
+#define G_COUNTS_SET                          "GCTAEventBin::counts(double&)"
+#define G_OMEGA                                       "GCTAEventBin::omega()"
+#define G_EWIDTH                                     "GCTAEventBin::ewidth()"
+#define G_ONTIME                                     "GCTAEventBin::ontime()"
  
 /* __ Macros _____________________________________________________________ */
  
@@ -185,8 +194,128 @@ double GCTAEventBin::size(void) const
  
  
 /***********************************************************************//**
+ * @brief Return instrument direction of event bin
+ *
+ * @return Instrument direction of event bin
+ *
+ * @exception GCTAException::no_member
+ *            Invalid instrument direction pointer.
+ *
+ * Returns reference to the instrument direction of the event bin.
+ ***************************************************************************/
+const GCTAInstDir& GCTAEventBin::dir(void) const
+{
+    // Throw an exception if instrument direction pointer is not valid
+    if (m_dir == NULL) {
+        throw GCTAException::no_member(G_DIR,
+                                       "Invalid instrument direction pointer.");
+    }
+
+    // Return instrument direction
+    return *m_dir;
+}
+
+
+/***********************************************************************//**
+ * @brief Return energy of event bin
+ *
+ * @return Energy of event bin
+ *
+ * @exception GCTAException::no_member
+ *            Invalid energy pointer.
+ *
+ * Returns reference to the energy of the event bin.
+ ***************************************************************************/
+const GEnergy& GCTAEventBin::energy(void) const
+{
+    // Throw an exception if energy pointer is not valid
+    if (m_energy == NULL) {
+        throw GCTAException::no_member(G_ENERGY,
+                                       "Invalid energy pointer.");
+    }
+
+    // Return energy
+    return *m_energy;
+}
+
+
+/***********************************************************************//**
+ * @brief Return time of event bin
+ *
+ * @return Time of event bin
+ *
+ * @exception GCTAException::no_member
+ *            Invalid time pointer.
+ *
+ * Returns reference to the time of the event bin.
+ ***************************************************************************/
+const GTime& GCTAEventBin::time(void) const
+{
+    // Throw an exception if time pointer is not valid
+    if (m_time == NULL) {
+        throw GCTAException::no_member(G_TIME,
+                                       "Invalid time pointer.");
+    }
+
+    // Return time
+    return *m_time;
+}
+
+
+/***********************************************************************//**
+ * @brief Return number of counts in event bin
+ *
+ * @return Number of counts in event bin
+ *
+ * @exception GCTAException::no_member
+ *            Invalid counts pointer.
+ *
+ * Returns reference to the number of counts in the event bin.
+ ***************************************************************************/
+double GCTAEventBin::counts(void) const
+{
+    // Throw an exception if counts pointer is not valid
+    if (m_counts == NULL) {
+        throw GCTAException::no_member(G_COUNTS_GET,
+                                       "Invalid counts pointer.");
+    }
+
+    // Return counts
+    return *m_counts;
+}
+
+
+/***********************************************************************//**
+ * @brief Set number of counts in event bin
+ *
+ * @param[in] counts Number of counts.
+ *
+ * @exception GCTAException::no_member
+ *            Invalid counts pointer.
+ *
+ * Set the number of counts in the event bin.
+ ***************************************************************************/
+void GCTAEventBin::counts(const double& counts)
+{
+    // Throw an exception if counts pointer is not valid
+    if (m_counts == NULL) {
+        throw GCTAException::no_member(G_COUNTS_SET,
+                                       "Invalid counts pointer.");
+    }
+
+    // Set number of counts in event bin
+    *m_counts = counts;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Return error in number of counts
  *
+ * @return Error in number of counts in event bin
+ *
  * Returns \f$\sqrt(counts+delta)\f$ as the uncertainty in the number of
  * counts in the bin. Adding delta avoids uncertainties of 0 which will
  * lead in the optimisation step to the exlusion of the corresponding bin.
@@ -206,6 +335,75 @@ double GCTAEventBin::error(void) const
  
  
 /***********************************************************************//**
+ * @brief Return solid angle of event bin
+ *
+ * @return Solid angle of event bin
+ *
+ * @exception GCTAException::no_member
+ *            Invalid solid angle pointer.
+ *
+ * Returns reference to the solid angle of the event bin.
+ ***************************************************************************/
+const double& GCTAEventBin::omega(void) const
+{
+    // Throw an exception if solid angle pointer is not valid
+    if (m_omega == NULL) {
+        throw GCTAException::no_member(G_OMEGA,
+                                       "Invalid solid angle pointer.");
+    }
+
+    // Return solid angle
+    return *m_omega;
+}
+
+
+/***********************************************************************//**
+ * @brief Return energy width of event bin
+ *
+ * @return Energy width of event bin
+ *
+ * @exception GCTAException::no_member
+ *            Invalid energy width pointer.
+ *
+ * Returns reference to the energy width of the event bin.
+ ***************************************************************************/
+const GEnergy& GCTAEventBin::ewidth(void) const
+{
+    // Throw an exception if energy width pointer is not valid
+    if (m_ewidth == NULL) {
+        throw GCTAException::no_member(G_EWIDTH,
+                                       "Invalid energy width pointer.");
+    }
+
+    // Return energy width
+    return *m_ewidth;
+}
+
+
+/***********************************************************************//**
+ * @brief Return ontime of event bin
+ *
+ * @return Ontime of event bin
+ *
+ * @exception GCTAException::no_member
+ *            Invalid ontime pointer.
+ *
+ * Returns reference to the ontime of the event bin.
+ ***************************************************************************/
+const double& GCTAEventBin::ontime(void) const
+{
+    // Throw an exception if ontime pointer is not valid
+    if (m_ontime == NULL) {
+        throw GCTAException::no_member(G_ONTIME,
+                                       "Invalid ontime pointer.");
+    }
+
+    // Return ontime
+    return *m_ontime;
+}
+
+
+/***********************************************************************//**
  * @brief Print event information
  ***************************************************************************/
 std::string GCTAEventBin::print(void) const
@@ -61,6 +61,26 @@ GCTAException::file_open_error::file_open_error(std::string origin,
  
  
 /***********************************************************************//**
+ * @brief Member not set.
+ *
+ * @param[in] origin Name of method that has thrown the exception.
+ * @param[in] message Error message.
+ ***************************************************************************/
+GCTAException::no_member::no_member(const std::string& origin,
+                                    const std::string& message)
+{
+    // Set origin and message
+    m_origin  = origin;
+    m_message = message;
+
+    // Return
+    return;
+}
+
+
+
+
+/***********************************************************************//**
  * @brief Pointing not set.
  *
  * @param[in] origin Name of method that has thrown the exception.
 /***************************************************************************
- *                 GLATEventBin.hpp  -  LAT event bin class                *
+ *                 GLATEventBin.hpp - LAT event bin class                  *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2009-2011 by Jurgen Knodlseder                           *
+ *  copyright (C) 2009-2012 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -21,7 +21,7 @@
 /**
  * @file GLATEventBin.hpp
  * @brief Fermi-LAT event bin class interface definition
- * @author J. Knodlseder
+ * @author Juergen Knoedlseder
  */
  
 #ifndef GLATEVENTBIN_HPP
@@ -62,18 +62,18 @@ public:
     virtual void               clear(void);
     virtual GLATEventBin*      clone(void) const;
     virtual double             size(void) const;
-    virtual const GLATInstDir& dir(void) const { return *m_dir; }
-    virtual const GEnergy&     energy(void) const { return *m_energy; }
-    virtual const GTime&       time(void) const { return *m_time; }
-    virtual double             counts(void) const { return *m_counts; }
+    virtual const GLATInstDir& dir(void) const;
+    virtual const GEnergy&     energy(void) const;
+    virtual const GTime&       time(void) const;
+    virtual double             counts(void) const;
     virtual double             error(void) const;
-    virtual void               counts(const double& counts) { *m_counts=counts; }
+    virtual void               counts(const double& counts);
     virtual std::string        print(void) const;
  
     // Other methods
-    const double&  omega(void) const { return *m_omega; }
-    const GEnergy& ewidth(void) const { return *m_ewidth; }
-    const double&  ontime(void) const { return *m_ontime; }
+    const double&  omega(void) const;
+    const GEnergy& ewidth(void) const;
+    const double&  ontime(void) const;
     const int&     index(void) const { return m_index; }
     const int&     ipix(void) const { return m_ipix; }
     const int&     ieng(void) const { return m_ieng; }
 /***************************************************************************
  *                GLATException.hpp  - LAT exception handler               *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2011 by Jurgen Knodlseder                           *
+ *  copyright (C) 2010-2012 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -21,7 +21,7 @@
 /**
  * @file GLATException.hpp
  * @brief LAT exception handler interface definition.
- * @author J. Knodlseder
+ * @author Juergen Knoedlseder
  */
  
 #ifndef GLATEXCEPTION_HPP
@@ -64,6 +64,10 @@ public:
         file_open_error(std::string origin, std::string filename,
                         std::string message = "");
     };
+    class no_member : public GLATExceptionHandler {
+    public:
+        no_member(const std::string& origin, const std::string& message);
+    };
     class no_sky : public GLATExceptionHandler {
     public:
         no_sky(std::string origin, std::string message = "");
@@ -539,8 +539,8 @@ void GLATAeff::copy_members(const GLATAeff&amp; aeff)
     m_back       = aeff.m_back;
  
     // Clone functors
-    m_eff_func1 = aeff.m_eff_func1->clone();
-    m_eff_func2 = aeff.m_eff_func2->clone();
+    m_eff_func1 = (aeff.m_eff_func1 != NULL) ? aeff.m_eff_func1->clone() : NULL;
+    m_eff_func2 = (aeff.m_eff_func2 != NULL) ? aeff.m_eff_func2->clone() : NULL;
  
     // Return
     return;
 /***************************************************************************
- *                GLATEventBin.cpp  -  LAT event bin class                 *
+ *                 GLATEventBin.cpp - LAT event bin class                  *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2009-2010 by Jurgen Knodlseder                           *
+ *  copyright (C) 2009-2012 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -20,8 +20,8 @@
  ***************************************************************************/
 /**
  * @file GLATEventBin.cpp
- * @brief GLATEventBin class implementation.
- * @author J. Knodlseder
+ * @brief Fermi-LAT event bin class implementation
+ * @author Juergen Knodlseder
  */
  
 /* __ Includes ___________________________________________________________ */
@@ -31,9 +31,18 @@
 #include <string>
 #include <cmath>
 #include "GLATEventBin.hpp"
+#include "GLATException.hpp"
 #include "GTools.hpp"
  
 /* __ Method name definitions ____________________________________________ */
+#define G_DIR                                           "GLATEventBin::dir()"
+#define G_ENERGY                                     "GLATEventBin::energy()"
+#define G_TIME                                         "GLATEventBin::time()"
+#define G_COUNTS_GET                                 "GLATEventBin::counts()"
+#define G_COUNTS_SET                          "GLATEventBin::counts(double&)"
+#define G_OMEGA                                       "GLATEventBin::omega()"
+#define G_EWIDTH                                     "GLATEventBin::ewidth()"
+#define G_ONTIME                                     "GLATEventBin::ontime()"
  
 /* __ Macros _____________________________________________________________ */
  
@@ -101,7 +110,8 @@ GLATEventBin::~GLATEventBin(void)
 /***********************************************************************//**
  * @brief Assignment operator
  *
- * @param[in] bin Event bin which should be assigned.
+ * @param[in] bin LAT event bin.
+ * @return Lat event bin.
  ***************************************************************************/
 GLATEventBin& GLATEventBin::operator= (const GLATEventBin& bin)
 {
@@ -157,7 +167,9 @@ void GLATEventBin::clear(void)
  
 /***********************************************************************//**
  * @brief Clone instance
-***************************************************************************/
+ *
+ * @return Deep copy of LAT event bin.
+ ***************************************************************************/
 GLATEventBin* GLATEventBin::clone(void) const
 {
     return new GLATEventBin(*this);
@@ -167,6 +179,8 @@ GLATEventBin* GLATEventBin::clone(void) const
 /***********************************************************************//**
  * @brief Return size of event bin
  *
+ * @return Size of event bin in units of sr MeV s.
+ *
  * The size of the event bin (units sr MeV s) is given by
  * \f[size = \Omega \times \Delta E \times \Delta T\f]
  * where
@@ -185,8 +199,128 @@ double GLATEventBin::size(void) const
  
  
 /***********************************************************************//**
+ * @brief Return instrument direction of event bin
+ *
+ * @return Instrument direction of event bin.
+ *
+ * @exception GLATException::no_member
+ *            Invalid instrument direction pointer.
+ *
+ * Returns reference to the instrument direction of the event bin.
+ ***************************************************************************/
+const GLATInstDir& GLATEventBin::dir(void) const
+{
+    // Throw an exception if instrument direction pointer is not valid
+    if (m_dir == NULL) {
+        throw GLATException::no_member(G_DIR,
+                                       "Invalid instrument direction pointer.");
+    }
+
+    // Return instrument direction
+    return *m_dir;
+}
+
+
+/***********************************************************************//**
+ * @brief Return energy of event bin
+ *
+ * @return Energy of event bin.
+ *
+ * @exception GLATException::no_member
+ *            Invalid energy pointer.
+ *
+ * Returns reference to the energy of the event bin.
+ ***************************************************************************/
+const GEnergy& GLATEventBin::energy(void) const
+{
+    // Throw an exception if energy pointer is not valid
+    if (m_energy == NULL) {
+        throw GLATException::no_member(G_ENERGY,
+                                       "Invalid energy pointer.");
+    }
+
+    // Return energy
+    return *m_energy;
+}
+
+
+/***********************************************************************//**
+ * @brief Return time of event bin
+ *
+ * @return Time of event bin.
+ *
+ * @exception GLATException::no_member
+ *            Invalid time pointer.
+ *
+ * Returns reference to the time of the event bin.
+ ***************************************************************************/
+const GTime& GLATEventBin::time(void) const
+{
+    // Throw an exception if time pointer is not valid
+    if (m_time == NULL) {
+        throw GLATException::no_member(G_TIME,
+                                       "Invalid time pointer.");
+    }
+
+    // Return time
+    return *m_time;
+}
+
+
+/***********************************************************************//**
+ * @brief Return number of counts in event bin
+ *
+ * @return Number of counts in event bin.
+ *
+ * @exception GLATException::no_member
+ *            Invalid counts pointer.
+ *
+ * Returns reference to the number of counts in the event bin.
+ ***************************************************************************/
+double GLATEventBin::counts(void) const
+{
+    // Throw an exception if counts pointer is not valid
+    if (m_counts == NULL) {
+        throw GLATException::no_member(G_COUNTS_GET,
+                                       "Invalid counts pointer.");
+    }
+
+    // Return counts
+    return *m_counts;
+}
+
+
+/***********************************************************************//**
+ * @brief Set number of counts in event bin
+ *
+ * @param[in] counts Number of counts.
+ *
+ * @exception GLATException::no_member
+ *            Invalid counts pointer.
+ *
+ * Set the number of counts in the event bin.
+ ***************************************************************************/
+void GLATEventBin::counts(const double& counts)
+{
+    // Throw an exception if counts pointer is not valid
+    if (m_counts == NULL) {
+        throw GLATException::no_member(G_COUNTS_SET,
+                                       "Invalid counts pointer.");
+    }
+
+    // Set number of counts in event bin
+    *m_counts = counts;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Return error in number of counts
  *
+ * @return Error in number of counts in event bin.
+ *
  * Returns \f$\sqrt(counts+delta)\f$ as the uncertainty in the number of
  * counts in the bin. Adding delta avoids uncertainties of 0 which will
  * lead in the optimisation step to the exlusion of the corresponding bin.
@@ -206,7 +340,78 @@ double GLATEventBin::error(void) const
  
  
 /***********************************************************************//**
+ * @brief Return solid angle of event bin
+ *
+ * @return Solid angle of event bin.
+ *
+ * @exception GLATException::no_member
+ *            Invalid solid angle pointer.
+ *
+ * Returns reference to the solid angle of the event bin.
+ ***************************************************************************/
+const double& GLATEventBin::omega(void) const
+{
+    // Throw an exception if solid angle pointer is not valid
+    if (m_omega == NULL) {
+        throw GLATException::no_member(G_OMEGA,
+                                       "Invalid solid angle pointer.");
+    }
+
+    // Return solid angle
+    return *m_omega;
+}
+
+
+/***********************************************************************//**
+ * @brief Return energy width of event bin
+ *
+ * @return Energy width of event bin.
+ *
+ * @exception GLATException::no_member
+ *            Invalid energy width pointer.
+ *
+ * Returns reference to the energy width of the event bin.
+ ***************************************************************************/
+const GEnergy& GLATEventBin::ewidth(void) const
+{
+    // Throw an exception if energy width pointer is not valid
+    if (m_ewidth == NULL) {
+        throw GLATException::no_member(G_EWIDTH,
+                                       "Invalid energy width pointer.");
+    }
+
+    // Return energy width
+    return *m_ewidth;
+}
+
+
+/***********************************************************************//**
+ * @brief Return ontime of event bin
+ *
+ * @return Ontime of event bin.
+ *
+ * @exception GLATException::no_member
+ *            Invalid ontime pointer.
+ *
+ * Returns reference to the ontime of the event bin.
+ ***************************************************************************/
+const double& GLATEventBin::ontime(void) const
+{
+    // Throw an exception if ontime pointer is not valid
+    if (m_ontime == NULL) {
+        throw GLATException::no_member(G_ONTIME,
+                                       "Invalid ontime pointer.");
+    }
+
+    // Return ontime
+    return *m_ontime;
+}
+
+
+/***********************************************************************//**
  * @brief Print event information
+ *
+ * @return String containing number of counts in event bin.
  ***************************************************************************/
 std::string GLATEventBin::print(void) const
 {
@@ -283,10 +488,3 @@ void GLATEventBin::free_members(void)
     // Return
     return;
 }
-
-
-/*==========================================================================
- =                                                                         =
- =                                Friends                                  =
- =                                                                         =
- ==========================================================================*/
 /***************************************************************************
- *                GLATEventCube.cpp  -  LAT event cube class               *
+ *                 GLATEventCube.cpp - LAT event cube class                *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2009-2012 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
@@ -21,7 +21,7 @@
 /**
  * @file GLATEventCube.cpp
  * @brief LAT event cube class implementation
- * @author J. Knoedlseder
+ * @author Juergen Knoedlseder
  */
  
 /* __ Includes ___________________________________________________________ */
@@ -111,6 +111,7 @@ GLATEventCube::~GLATEventCube(void)
  * @brief Assignment operator
  *
  * @param[in] cube LAT event cube.
+ * @return LAT event cube.
  ***************************************************************************/
 GLATEventCube& GLATEventCube::operator= (const GLATEventCube& cube)
 {
@@ -140,6 +141,7 @@ GLATEventCube&amp; GLATEventCube::operator= (const GLATEventCube&amp; cube)
  * @brief Event bin access operator
  *
  * @param[in] index Event index [0,...,size()-1].
+ * @return Pointer to event bin.
  *
  * Returns pointer to an event bin.
  ***************************************************************************/
@@ -157,6 +159,7 @@ GLATEventBin* GLATEventCube::operator[](const int&amp; index)
  * @brief Event bin access operator (const version)
  *
  * @param[in] index Event index [0,...,size()-1].
+ * @return Pointer to event bin.
  *
  * Returns pointer to an event bin.
  ***************************************************************************/
@@ -200,6 +203,8 @@ void GLATEventCube::clear(void)
  
 /***********************************************************************//**
  * @brief Clone instance
+ *
+ * @return Deep copy of LAT event cube.
  ***************************************************************************/
 GLATEventCube* GLATEventCube::clone(void) const
 {
@@ -209,6 +214,8 @@ GLATEventCube* GLATEventCube::clone(void) const
  
 /***********************************************************************//**
  * @brief Return number of bins in event cube
+ *
+ * @return Number of bins in event cube.
  ***************************************************************************/
 int GLATEventCube::size(void) const
 {
@@ -222,6 +229,8 @@ int GLATEventCube::size(void) const
  
 /***********************************************************************//**
  * @brief Return dimension of event cube
+ *
+ * @return Number of dimensions in event cube (either 2 or 3).
  ***************************************************************************/
 int GLATEventCube::dim(void) const
 {
@@ -236,7 +245,8 @@ int GLATEventCube::dim(void) const
 /***********************************************************************//**
  * @brief Return number of bins in axis
  *
- * @param[in] axis Axis.
+ * @param[in] axis Axis [0,...,dim()-1]
+ * @return Number of bins in specified axis
  *
  * @exception GException::out_of_range
  *            Axis is out of range.
@@ -247,8 +257,9 @@ int GLATEventCube::naxis(int axis) const
 {
     // Optionally check if the axis is valid
     #if defined(G_RANGE_CHECK)
-    if (axis < 0 || axis >= dim())
+    if (axis < 0 || axis >= dim()) {
         throw GException::out_of_range(G_NAXIS, axis, 0, dim()-1);
+    }
     #endif
  
     // Set result
@@ -368,6 +379,8 @@ void GLATEventCube::write(GFits&amp; file) const
  
 /***********************************************************************//**
  * @brief Return number of events in cube
+ *
+ * @return Total number of events in event cube.
  ***************************************************************************/
 int GLATEventCube::number(void) const
 {
@@ -379,8 +392,9 @@ int GLATEventCube::number(void) const
  
     // Sum event cube
     if (size() > 0 && pixels != NULL) {
-        for (int i = 0; i < size(); ++i)
+        for (int i = 0; i < size(); ++i) {
             number += pixels[i];
+        }
     }
  
     // Return
@@ -406,6 +420,8 @@ void GLATEventCube::map(const GSkymap&amp; map)
  
 /***********************************************************************//**
  * @brief Print event cube information
+ *
+ * @return String containing event cube information.
  ***************************************************************************/
 std::string GLATEventCube::print(void) const
 {
@@ -419,17 +435,29 @@ std::string GLATEventCube::print(void) const
     result.append(str(m_map.nx())+" x "+str(m_map.ny()));
     result.append("\n"+parformat("Number of energy bins")+str(ebins()));
     result.append("\n"+parformat("Number of events")+str(number()));
+
+    // Append time interval
     result.append("\n"+parformat("Time interval"));
-    if (gti().size() > 0)
+    if (gti().size() > 0) {
         result.append(str(tstart().met())+" - "+str(tstop().met()));
-    else
+    }
+    else {
         result.append("not defined");
+    }
+
+    // Append energy range
     result.append("\n"+parformat("Energy range"));
-    if (ebounds().size() > 0)
+    if (ebounds().size() > 0) {
         result.append(emin().print()+" - "+emax().print());
-    else
+    }
+    else {
         result.append("not defined");
-    result.append("\n"+m_map.wcs()->print());
+    }
+
+    // Append WCS
+    if (m_map.wcs() != NULL) {
+        result.append("\n"+m_map.wcs()->print());
+    }
  
     // Append source maps
     result.append("\n"+parformat("Number of source maps")+str(m_srcmap.size()));
@@ -450,7 +478,8 @@ std::string GLATEventCube::print(void) const
 /***********************************************************************//**
  * @brief Return name of diffuse model
  *
- * @param[in] index Diffuse model index (starting from 0).
+ * @param[in] index Diffuse model index [0,...,ndiffrsp()-1].
+ * @return Name of diffuse model.
  *
  * @exception GException::out_of_range
  *            Model index out of valid range.
@@ -461,8 +490,9 @@ std::string GLATEventCube::diffname(const int&amp; index) const
 {
     // Optionally check if the index is valid
     #if defined(G_RANGE_CHECK)
-    if (index < 0 || index >= ndiffrsp())
+    if (index < 0 || index >= ndiffrsp()) {
         throw GException::out_of_range(G_DIFFNAME, index, 0, ndiffrsp()-1);
+    }
     #endif
  
     // Return
@@ -473,7 +503,8 @@ std::string GLATEventCube::diffname(const int&amp; index) const
 /***********************************************************************//**
  * @brief Return diffuse response map
  *
- * @param[in] index Diffuse model index (starting from 0).
+ * @param[in] index Diffuse model index [0,...,ndiffrsp()-1].
+ * @return Pointer to diffuse response map.
  *
  * @exception GException::out_of_range
  *            Model index out of valid range.
@@ -484,8 +515,9 @@ GSkymap* GLATEventCube::diffrsp(const int&amp; index) const
 {
     // Optionally check if the index is valid
     #if defined(G_RANGE_CHECK)
-    if (index < 0 || index >= ndiffrsp())
+    if (index < 0 || index >= ndiffrsp()) {
         throw GException::out_of_range(G_DIFFRSP, index, 0, ndiffrsp()-1);
+    }
     #endif
  
     // Return
@@ -498,6 +530,7 @@ GSkymap* GLATEventCube::diffrsp(const int&amp; index) const
  *        direction that fits spatially into the event cube
  *
  * @param[in] srcDir Source direction.
+ * @return Maximum radius in degrees that fully fits into event cube.
  *
  * By computing the sky directions of the event cube boundaries, the maximum
  * radius is computed that fits fully within the event cube. This method is
@@ -513,8 +546,9 @@ double GLATEventCube::maxrad(const GSkyDir&amp; srcDir) const
     for (int ix = 0; ix < nx(); ++ix) {
         GSkyPixel pixel    = GSkyPixel(double(ix), double(iy));
         double    distance = m_map.xy2dir(pixel).dist_deg(srcDir);
-        if (distance < radius)
+        if (distance < radius) {
             radius = distance;
+        }
     }
  
     // Move along lower edge in longitude
@@ -522,8 +556,9 @@ double GLATEventCube::maxrad(const GSkyDir&amp; srcDir) const
     for (int ix = 0; ix < nx(); ++ix) {
         GSkyPixel pixel    = GSkyPixel(double(ix), double(iy));
         double    distance = m_map.xy2dir(pixel).dist_deg(srcDir);
-        if (distance < radius)
+        if (distance < radius) {
             radius = distance;
+        }
     }
  
     // Move along left edge in latitude
@@ -531,8 +566,9 @@ double GLATEventCube::maxrad(const GSkyDir&amp; srcDir) const
     for (int iy = 0; iy < ny(); ++iy) {
         GSkyPixel pixel    = GSkyPixel(double(ix), double(iy));
         double    distance = m_map.xy2dir(pixel).dist_deg(srcDir);
-        if (distance < radius)
+        if (distance < radius) {
             radius = distance;
+        }
     }
  
     // Move along right edge in latitude
@@ -540,8 +576,9 @@ double GLATEventCube::maxrad(const GSkyDir&amp; srcDir) const
     for (int iy = 0; iy < ny(); ++iy) {
         GSkyPixel pixel    = GSkyPixel(double(ix), double(iy));
         double    distance = m_map.xy2dir(pixel).dist_deg(srcDir);
-        if (distance < radius)
+        if (distance < radius) {
             radius = distance;
+        }
     }
  
     // Return radius
@@ -664,18 +701,21 @@ void GLATEventCube::read_srcmap(const GFitsImage* hdu)
         map->read(hdu);
  
         // Check that source map WCS is consistent with counts map WCS
-        if (*(m_map.wcs()) != *(map->wcs()))
+        if (*(m_map.wcs()) != *(map->wcs())) {
             throw GLATException::wcs_incompatible(G_READ_SRCMAP, hdu->extname());
+        }
  
         // Check that source map dimension is consistent with counts map
         // dimension
         if (m_map.nx() != map->nx() ||
-            m_map.ny() != map->ny())
+            m_map.ny() != map->ny()) {
             throw GLATException::wcs_incompatible(G_READ_SRCMAP, hdu->extname());
+        }
  
         // Check that source map has required number of energy bins
-        if (m_map.nmaps()+1 != map->nmaps())
+        if (m_map.nmaps()+1 != map->nmaps()) {
             throw GLATException::wcs_incompatible(G_READ_SRCMAP, hdu->extname());
+        }
  
         // Append source map to list of maps
         m_srcmap.push_back(map);
@@ -837,9 +877,10 @@ void GLATEventCube::set_energies(void)
 void GLATEventCube::set_times(void)
 {
     // Throw an error if GTI is empty
-    if (m_gti.size() < 1)
+    if (m_gti.size() < 1) {
         throw GLATException::no_gti(G_SET_TIMES, "Every LAT event cube needs"
                   " associated GTIs to allow the computation of the ontime.");
+    }
  
     // Compute mean time
     m_time = 0.5 * (gti().tstart() + gti().tstop());
@@ -874,17 +915,20 @@ void GLATEventCube::set_bin(const int&amp; index)
 {
     // Optionally check if the index is valid
     #if defined(G_RANGE_CHECK)
-    if (index < 0 || index >= size())
+    if (index < 0 || index >= size()) {
         throw GException::out_of_range(G_SET_BIN, index, 0, size()-1);
+    }
     #endif
  
     // Check for the existence of energies and energy widths
-    if (m_energies.size() != ebins() || m_ewidth.size() != ebins())
+    if (m_energies.size() != ebins() || m_ewidth.size() != ebins()) {
         throw GLATException::no_energies(G_SET_BIN);
+    }
  
     // Check for the existence of sky directions and solid angles
-    if (m_dirs.size() != npix() || m_omega.size() != npix())
+    if (m_dirs.size() != npix() || m_omega.size() != npix()) {
         throw GLATException::no_dirs(G_SET_BIN);
+    }
  
     // Get pixel and energy bin indices.
     m_bin.m_index = index;
@@ -904,10 +948,3 @@ void GLATEventCube::set_bin(const int&amp; index)
     // Return
     return;
 }
-
-
-/*==========================================================================
- =                                                                         =
- =                                Friends                                  =
- =                                                                         =
- ==========================================================================*/
@@ -59,6 +59,24 @@ GLATException::file_open_error::file_open_error(std::string origin,
  
  
 /***********************************************************************//**
+ * @brief Member not set.
+ *
+ * @param[in] origin Name of method that has thrown the exception.
+ * @param[in] message Error message.
+ ***************************************************************************/
+GLATException::no_member::no_member(const std::string& origin,
+                                    const std::string& message)
+{
+    // Set origin and message
+    m_origin  = origin;
+    m_message = message;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief No sky pixels found.
  *
  * @param[in] origin Name of method that has thrown the exception.
@@ -21,7 +21,7 @@
 /**
  * @file GLATResponseTable.cpp
  * @brief Fermi/LAT response table class implementation
- * @author J. Knoedlseder
+ * @author Juergen Knoedlseder
  */
  
 /* __ Includes ___________________________________________________________ */
@@ -159,7 +159,7 @@ void GLATResponseTable::clear(void)
  
 /***********************************************************************//**
  * @brief Clone instance
-***************************************************************************/
+ ***************************************************************************/
 GLATResponseTable* GLATResponseTable::clone(void) const
 {
     return new GLATResponseTable(*this);
@@ -290,9 +290,10 @@ int GLATResponseTable::index(const int&amp; ie, const int&amp; ic) const
     // Optionally check if the index is valid
     #if defined(G_RANGE_CHECK)
     if (ie < 0 || ie >= m_energy_num ||
-        ic < 0 || ic >= m_ctheta_num)
+        ic < 0 || ic >= m_ctheta_num) {
         throw GException::out_of_range(G_INDEX, ie, ic, 
                                        m_energy_num-1, m_ctheta_num-1);
+    }
     #endif
  
     // Compute index
@@ -323,8 +324,9 @@ double GLATResponseTable::energy(const int&amp; ie) const
 {
     // Optionally check if the index is valid
     #if defined(G_RANGE_CHECK)
-    if (ie < 0 || ie >= m_energy_num)
+    if (ie < 0 || ie >= m_energy_num) {
         throw GException::out_of_range(G_ENERGY, ie, 0, m_energy_num-1);
+    }
     #endif
  
     // Determine mean energy of bin
@@ -462,8 +464,9 @@ double GLATResponseTable::energy_lo(const int&amp; inx) const
 {
     // Optionally check if the index is valid
     #if defined(G_RANGE_CHECK)
-    if (inx < 0 || inx >= m_energy_num)
+    if (inx < 0 || inx >= m_energy_num) {
         throw GException::out_of_range(G_ENERGY_LO, inx, 0, m_energy_num-1);
+    }
     #endif
  
     // Return mean energy
@@ -483,8 +486,9 @@ double GLATResponseTable::energy_hi(const int&amp; inx) const
 {
     // Optionally check if the index is valid
     #if defined(G_RANGE_CHECK)
-    if (inx < 0 || inx >= m_energy_num)
+    if (inx < 0 || inx >= m_energy_num) {
         throw GException::out_of_range(G_ENERGY_HI, inx, 0, m_energy_num-1);
+    }
     #endif
  
     // Return mean energy
@@ -504,8 +508,9 @@ double GLATResponseTable::costheta_lo(const int&amp; inx) const
 {
     // Optionally check if the index is valid
     #if defined(G_RANGE_CHECK)
-    if (inx < 0 || inx >= m_ctheta_num)
+    if (inx < 0 || inx >= m_ctheta_num) {
         throw GException::out_of_range(G_COSTHETA_LO, inx, 0, m_ctheta_num-1);
+    }
     #endif
  
     // Return mean energy
@@ -525,8 +530,9 @@ double GLATResponseTable::costheta_hi(const int&amp; inx) const
 {
     // Optionally check if the index is valid
     #if defined(G_RANGE_CHECK)
-    if (inx < 0 || inx >= m_ctheta_num)
+    if (inx < 0 || inx >= m_ctheta_num) {
         throw GException::out_of_range(G_COSTHETA_HI, inx, 0, m_ctheta_num-1);
+    }
     #endif
  
     // Return mean energy
@@ -562,10 +568,12 @@ std::vector&lt;double&gt; GLATResponseTable::energies(void) const
     std::vector<double> energies;
  
     // Push energies on vector
-    energies.push_back(m_energy[m_logE.inx_left()]);
-    energies.push_back(m_energy[m_logE.inx_left()]);
-    energies.push_back(m_energy[m_logE.inx_right()]);
-    energies.push_back(m_energy[m_logE.inx_right()]);
+    if (m_energy_num > 0) {
+        energies.push_back(m_energy[m_logE.inx_left()]);
+        energies.push_back(m_energy[m_logE.inx_left()]);
+        energies.push_back(m_energy[m_logE.inx_right()]);
+        energies.push_back(m_energy[m_logE.inx_right()]);
+    }
  
     // Return vector
     return energies;
@@ -66,7 +66,7 @@ public:
     GEnergy   emax(int inx) const;
     GEnergy   emean(int inx) const;
     GEnergy   elogmean(int inx) const;
-    bool      isin(const GEnergy& eng) const;
+    bool      contains(const GEnergy& eng) const;
 };
  
  
 /***************************************************************************
- *                    GFits.cpp  - FITS file access class                  *
+ *                    GFits.cpp - FITS file access class                   *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2008-2012 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
@@ -348,10 +348,13 @@ void GFits::open(const std::string&amp; filename, bool create)
  * @exception GException::fits_file_exist
  *            Attemting to overwrite an existing file without having specified
  *            clobber=true.
+ * @exception GException::fits_file_not_open
+ *            FITS file needs to be opened before saving.
+ *
+ * Saves all HDUs of an open FITS file to disk. After saving, the FITS file
+ * remains open. Invoke the close() method if explicit closing is needed.
+ * Note that de-allocation of the GFits object also closes the FITS file.
  *
- * Saves all HDUs to the FITS file by looping over the HDU's save() method.
- * Saving does not mean that the FITS file is closed. Call the close()
- * method if explicit closing is needed (de-allocation also closes the file).
  * In the special case that no first HDU exists an empty primary image is
  * created.
  ***************************************************************************/
@@ -368,6 +371,12 @@ void GFits::save(bool clobber)
         throw GException::fits_file_exist(G_SAVE, m_filename);
     }
  
+    // If no FITS file has been opened then throw an error
+    if (m_fitsfile == NULL) {
+        throw GException::fits_file_not_open(G_SAVE, 
+              "Either open FITS file before saving or use saveto() method.");
+    }
+
     // If no HDUs exist then save an empty primary image.
     if (size() == 0) {
         int status = 0;
@@ -537,28 +546,33 @@ void GFits::append(const GFitsHDU&amp; hdu)
     // Clone FITS HDU
     GFitsHDU* ptr = hdu.clone();
  
-    // If FITS file exists then connect cloned HDU to FITS file
-    if (m_fitsfile != NULL) {
-        __fitsfile fptr  = *(FPTR(m_fitsfile));
-        fptr.HDUposition = n_hdu;
-        ptr->connect(&fptr);
-    }
-    else {
-        ptr->extno(n_hdu);
-    }
+    // Append HDU if it's valid
+    if (ptr != NULL) {
  
-    // Debug information
-    #if DEBUG
-    std::cout << "Append HDU (extno=" << n_hdu << ")." << std::endl;
-    #endif
+        // If FITS file exists then connect cloned HDU to FITS file
+        if (m_fitsfile != NULL) {
+            __fitsfile fptr  = *(FPTR(m_fitsfile));
+            fptr.HDUposition = n_hdu;
+            ptr->connect(&fptr);
+        }
+        else {
+            ptr->extno(n_hdu);
+        }
  
-    // Push back HDU
-    m_hdu.push_back(ptr);
+        // Debug information
+        #if DEBUG
+        std::cout << "Append HDU (extno=" << n_hdu << ")." << std::endl;
+        #endif
  
-    // Debug trailer
-    #if DEBUG
-    std::cout << "<-- GFits::append" << std::endl;
-    #endif
+        // Push back HDU
+        m_hdu.push_back(ptr);
+
+        // Debug trailer
+        #if DEBUG
+        std::cout << "<-- GFits::append" << std::endl;
+        #endif
+        
+    } // endif: HDU was valid
  
     // Return
     return;
@@ -358,6 +358,8 @@ void GFitsHDU::connect(void* vptr)
 /***********************************************************************//**
  * @brief Move FITS file pointer to HDU
  *
+ * @exception GException::fits_file_not_open
+ *            No FITS file has been opened.
  * @exception GException::fits_hdu_not_found
  *            Requested HDU not found.
  *
@@ -366,6 +368,12 @@ void GFitsHDU::connect(void* vptr)
  ***************************************************************************/
 void GFitsHDU::move_to_hdu(void)
 {
+    // Throw an exception if FITS file is not open
+    if (FPTR(m_fitsfile)->Fptr == NULL) {
+        throw GException::fits_file_not_open(G_MOVE_TO_HDU, 
+              "Open file before moving file pointer to HDU.");
+    }
+
     // Move to HDU
     int status = 0;
     status     = __ffmahd(FPTR(m_fitsfile), m_hdunum+1, NULL, &status);
@@ -382,13 +390,23 @@ void GFitsHDU::move_to_hdu(void)
 /***********************************************************************//**
  * @brief Get HDU type from FITS file
  *
+ * @exception GException::fits_file_not_open
+ *            No FITS file has been opened.
  * @exception GException::fits_hdu_not_found
  *            Requested HDU not found.
  ***************************************************************************/
 GFitsHDU::HDUType GFitsHDU::get_hdu_type(void) const
 {
+    // Initialise HDU type
+    int type = 0;
+
+    // Throw an exception if FITS file is not open
+    if (FPTR(m_fitsfile)->Fptr == NULL) {
+        throw GException::fits_file_not_open(G_GET_HDU_TYPE, 
+              "Open file before requesting HDU type.");
+    }
+
     // Get HDU type
-    int type   = 0;
     int status = 0;
     status     = __ffghdt(FPTR(m_fitsfile), &type, &status);
     if (status != 0) {
 /***************************************************************************
- *                  GFitsImage.cpp  - FITS image class                     *
+ *                   GFitsImage.cpp - FITS image class                     *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2008-2012 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
@@ -371,8 +371,9 @@ int GFitsImage::naxes(int axis) const
 {
     // Check if axis is within the range
     #if defined(G_RANGE_CHECK)
-    if (axis < 0 || axis >= m_naxis)
+    if (axis < 0 || axis >= m_naxis) {
         throw GException::out_of_range(G_NAXES, axis, 0, m_naxis-1);
+    }
     #endif
  
     // Get axis dimension
@@ -448,9 +449,10 @@ std::string GFitsImage::print(void) const
     result.append(parformat("Image type")+typecode(type())+"\n");
     result.append(parformat("Number of dimensions")+str(naxis())+"\n");
     result.append(parformat("Number of image pixels")+str(size()));
-    for (int i = 0; i < naxis(); ++i)
+    for (int i = 0; i < naxis(); ++i) {
         result.append("\n"+parformat("Number of bins in "+str(i)) +
                       str(naxes(i)));
+    }
  
     // Append header information
     result.append(+"\n"+m_header.print());
@@ -577,6 +579,9 @@ void GFitsImage::init_image_header(void)
  *
  * @param[in] vptr FITS file void pointer.
  *
+ * @exception GException::fits_error
+ *            FITS error.
+ *
  * Open FITS image in FITS file. Opening means connecting the FITS file
  * pointer to the image and reading the image and axes dimensions.
  ***************************************************************************/
@@ -585,8 +590,9 @@ void GFitsImage::open_image(void* vptr)
     // Move to HDU
     int status = 0;
     status     = __ffmahd(FPTR(vptr), (FPTR(vptr)->HDUposition)+1, NULL, &status);
-    if (status != 0)
+    if (status != 0) {
         throw GException::fits_error(G_OPEN_IMAGE, status);
+    }
  
     // Save the FITS file pointer and the HDU number
     FPTR_COPY(m_fitsfile, vptr);
@@ -594,8 +600,9 @@ void GFitsImage::open_image(void* vptr)
  
     // Get the image dimensions
     status = __ffgidm(FPTR(m_fitsfile), &m_naxis, &status);
-    if (status != 0)
+    if (status != 0) {
         throw GException::fits_error(G_OPEN_IMAGE, status);
+    }
  
     // Reset number of image pixels
     m_num_pixels = 0;
@@ -609,13 +616,15 @@ void GFitsImage::open_image(void* vptr)
  
         // Get the axes dimensions
         status = __ffgisz(FPTR(m_fitsfile), m_naxis, m_naxes, &status);
-        if (status != 0)
+        if (status != 0) {
             throw GException::fits_error(G_OPEN_IMAGE, status);
+        }
  
         // Calculate number of image pixels
         m_num_pixels = 1;
-        for (int i = 0; i < m_naxis; ++i)
+        for (int i = 0; i < m_naxis; ++i) {
             m_num_pixels *= m_naxes[i];
+        }
  
     } // endif: there is an image
  
@@ -632,6 +641,9 @@ void GFitsImage::open_image(void* vptr)
  * @param[in] nulval Pointer to pixel nul value.
  * @param[out] anynul Number of nul values encountered during loading.
  *
+ * @exception GException::fits_error
+ *            FITS error.
+ *
  * Load image pixels from FITS file.
  ***************************************************************************/
 void GFitsImage::load_image(int datatype, const void* pixels,
@@ -656,8 +668,9 @@ void GFitsImage::load_image(int datatype, const void* pixels,
         delete [] fpixel;
         delete [] lpixel;
         delete [] inc;
-        if (status != 0)
+        if (status != 0) {
             throw GException::fits_error(G_LOAD_IMAGE, status);
+        }
     }
  
     // Return
@@ -671,12 +684,23 @@ void GFitsImage::load_image(int datatype, const void* pixels,
  * @param[in] datatype Datatype of pixels to be saved
  * @param[in] pixels Pixel array to be saved
  *
+ * @exception GException::fits_file_not_open
+ *            No FITS file has been opened.
+ * @exception GException::fits_error
+ *            FITS error.
+ *
  * Save image pixels into FITS file. In case that the HDU does not exist it
  * is created. In case that the pixel array is empty no data are saved; all
  * image pixels will be empty in this case.
  ***************************************************************************/
 void GFitsImage::save_image(int datatype, const void* pixels)
 {
+    // Throw an exception if FITS file is not open
+    if (FPTR(m_fitsfile)->Fptr == NULL) {
+        throw GException::fits_file_not_open(G_SAVE_IMAGE, 
+              "Open file before saving the image.");
+    }
+
     // Move to HDU. We use here an explicit cfitsio moveto function since we
     // want to recover the error code ...
     int status = 0;
@@ -686,23 +710,27 @@ void GFitsImage::save_image(int datatype, const void* pixels)
     if (status == 107) {
         status = 0;
         status = __ffcrim(FPTR(m_fitsfile), m_bitpix, m_naxis, m_naxes, &status);
-        if (status != 0)
+        if (status != 0) {
             throw GException::fits_error(G_SAVE_IMAGE, status);
+        }
     }
-    else if (status != 0)
+    else if (status != 0) {
         throw GException::fits_error(G_SAVE_IMAGE, status);
+    }
  
     // If HDU seems to be empty then create it now. This is only needed for the
     // primary HDU, since __ffmahd gives no error if the primary HDU is empty.
     // By checking the number of keywords in the HDU we detect an empty HDU ...
     int num = 0;
     status  = __ffghsp(FPTR(m_fitsfile), &num, NULL, &status);
-    if (status != 0)
+    if (status != 0) {
         throw GException::fits_error(G_SAVE_IMAGE, status);
+    }
     if (num == 0) {
         status = __ffcrim(FPTR(m_fitsfile), m_bitpix, m_naxis, m_naxes, &status);
-        if (status != 0)
+        if (status != 0) {
             throw GException::fits_error(G_SAVE_IMAGE, status);
+        }
     }
  
     // Save the image pixels (if there are some ...)
@@ -717,8 +745,9 @@ void GFitsImage::save_image(int datatype, const void* pixels)
                          (void*)pixels, &status);
         delete [] fpixel;
         delete [] lpixel;
-        if (status != 0)
+        if (status != 0) {
             throw GException::fits_error(G_SAVE_IMAGE, status);
+        }
     }
  
     // Return
@@ -749,8 +778,9 @@ void GFitsImage::fetch_data(void)
         init_data();
  
         // If a FITS file is attached then load pixels from FITS file.
-        if (FPTR(m_fitsfile)->Fptr != NULL)
+        if (FPTR(m_fitsfile)->Fptr != NULL) {
             load_image(type(), ptr_data(), ptr_nulval(), &m_anynul);
+        }
  
     } // endif: there were pixels available
  
@@ -771,8 +801,9 @@ int GFitsImage::offset(const int&amp; ix) const
 {
     // Check if axis is within the range
     #if defined(G_RANGE_CHECK)
-    if (ix < 0 || ix >= m_num_pixels)
+    if (ix < 0 || ix >= m_num_pixels) {
         throw GException::out_of_range(G_OFFSET_1D, ix, 0, m_num_pixels-1);
+    }
     #endif
  
     // Return index
@@ -797,15 +828,18 @@ int GFitsImage::offset(const int&amp; ix) const
 int GFitsImage::offset(const int& ix, const int& iy) const
 {
     // Operator is only valid for 2D images
-    if (m_naxis < 2)
+    if (m_naxis < 2) {
         throw GException::fits_wrong_image_operator(G_OFFSET_2D, m_naxis, 2);
+    }
  
     // Check if axis is within the range
     #if defined(G_RANGE_CHECK)
-    if (ix < 0 || ix >= m_naxes[0])
+    if (ix < 0 || ix >= m_naxes[0]) {
         throw GException::out_of_range(G_OFFSET_2D, ix, 0, m_naxes[0]-1);
-    if (iy < 0 || iy >= m_naxes[1])
+    }
+    if (iy < 0 || iy >= m_naxes[1]) {
         throw GException::out_of_range(G_OFFSET_2D, iy, 0, m_naxes[1]-1);
+    }
     #endif
  
     // Return offset
@@ -831,17 +865,21 @@ int GFitsImage::offset(const int&amp; ix, const int&amp; iy) const
 int GFitsImage::offset(const int& ix, const int& iy, const int& iz) const
 {
     // Operator is only valid for 3D images
-    if (m_naxis < 3)
+    if (m_naxis < 3) {
         throw GException::fits_wrong_image_operator(G_OFFSET_3D, m_naxis, 3);
+    }
  
     // Check if axis is within the range
     #if defined(G_RANGE_CHECK)
-    if (ix < 0 || ix >= m_naxes[0])
+    if (ix < 0 || ix >= m_naxes[0]) {
         throw GException::out_of_range(G_OFFSET_3D, ix, 0, m_naxes[0]-1);
-    if (iy < 0 || iy >= m_naxes[1])
+    }
+    if (iy < 0 || iy >= m_naxes[1]) {
         throw GException::out_of_range(G_OFFSET_3D, iy, 0, m_naxes[1]-1);
-    if (iz < 0 || iz >= m_naxes[2])
+    }
+    if (iz < 0 || iz >= m_naxes[2]) {
         throw GException::out_of_range(G_OFFSET_3D, iz, 0, m_naxes[2]-1);
+    }
     #endif
  
     // Return offset
@@ -869,19 +907,24 @@ int GFitsImage::offset(const int&amp; ix, const int&amp; iy, const int&amp; iz,
                        const int& it) const
 {
     // Operator is only valid for 4D images
-    if (m_naxis < 4)
+    if (m_naxis < 4) {
         throw GException::fits_wrong_image_operator(G_OFFSET_4D, m_naxis, 4);
+    }
  
     // Check if axis is within the range
     #if defined(G_RANGE_CHECK)
-    if (ix < 0 || ix >= m_naxes[0])
+    if (ix < 0 || ix >= m_naxes[0]) {
         throw GException::out_of_range(G_OFFSET_4D, ix, 0, m_naxes[0]-1);
-    if (iy < 0 || iy >= m_naxes[1])
+    }
+    if (iy < 0 || iy >= m_naxes[1]) {
         throw GException::out_of_range(G_OFFSET_4D, iy, 0, m_naxes[1]-1);
-    if (iz < 0 || iz >= m_naxes[2])
+    }
+    if (iz < 0 || iz >= m_naxes[2]) {
         throw GException::out_of_range(G_OFFSET_4D, iz, 0, m_naxes[2]-1);
-    if (it < 0 || it >= m_naxes[3])
+    }
+    if (it < 0 || it >= m_naxes[3]) {
         throw GException::out_of_range(G_OFFSET_4D, it, 0, m_naxes[3]-1);
+    }
     #endif
  
     // Return offset
@@ -929,8 +972,9 @@ void GFitsImage::copy_members(const GFitsImage&amp; image)
     m_naxes = NULL;
     if (image.m_naxes != NULL && m_naxis > 0) {
         m_naxes = new long[m_naxis];
-        for (int i = 0; i < m_naxis; ++i)
+        for (int i = 0; i < m_naxis; ++i) {
             m_naxes[i] = image.m_naxes[i];
+        }
     }
  
     // Return
@@ -312,8 +312,13 @@ GModels* GModels::clone(void) const
  ***************************************************************************/
 void GModels::append(const GModel& model)
 {
-    // Append model
-    m_models.push_back(model.clone());
+    // Make deep copy of model
+    GModel* model_copy = model.clone();
+
+    // If copy is valid, append it to the container
+    if (model_copy != NULL) {
+        m_models.push_back(model_copy);
+    }
  
     // Set parameter pointers
     set_pointers();
@@ -578,7 +578,7 @@ GEnergy GEbounds::elogmean(int inx) const
  *
  * @param[in] eng Energy to be tested.
  ***************************************************************************/
-bool GEbounds::isin(const GEnergy& eng) const
+bool GEbounds::contains(const GEnergy& eng) const
 {
     // Initialise test
     bool found = false;
@@ -228,7 +228,7 @@ double GEnergy::TeV(void) const
 double GEnergy::log10keV(void) const
 {
     // Return log10 energy
-    return (log10keV()+3.0); 
+    return (log10MeV()+3.0); 
 }
  
  
 /***************************************************************************
- *            GOptimizerLM.cpp  -  Levenberg Marquardt optimizer           *
+ *             GOptimizerLM.cpp - Levenberg Marquardt optimizer            *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2009-2011 by Jurgen Knodlseder                           *
+ *  copyright (C) 2009-2012 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -21,7 +21,7 @@
 /**
  * @file GOptimizerLM.cpp
  * @brief Levenberg-Marquardt optimizer class implementation
- * @author J. Knodlseder
+ * @author Juergen Knoedlseder
  */
  
 /* __ Includes ___________________________________________________________ */
@@ -948,108 +948,108 @@ double GOptimizerLM::step_size(GVector* grad, GOptimizerPars* pars)
  ***************************************************************************/
 void GOptimizerLM::errors(GOptimizerFunction* fct, GOptimizerPars* pars)
 {
-    // Get number of parameters
-    int npars = pars->npars();
+    // Continue only if pointers are valid
+    if (fct != NULL && pars != NULL) {
  
-    // Perform final parameter evaluation
-    fct->eval(*pars);
+        // Get number of parameters
+        int npars = pars->npars();
  
-    // Fetch sparse matrix pointer. We have to do this after the eval()
-    // method since eval() will allocate new memory for the covariance
-    // matrix!
-    GSparseMatrix* covar = fct->covar();
+        // Perform final parameter evaluation
+        fct->eval(*pars);
  
-    // Save best fitting value
-    m_value = *(fct->value());
+        // Fetch sparse matrix pointer. We have to do this after the eval()
+        // method since eval() will allocate new memory for the covariance
+        // matrix!
+        GSparseMatrix* covar = fct->covar();
  
-    // Save covariance matrix
-    GSparseMatrix save_covar = GSparseMatrix(*covar);
+        // Save best fitting value
+        m_value = *(fct->value());
  
-    // Signal no diagonal element loading
-    bool diag_loaded = false;
+        // Save covariance matrix
+        GSparseMatrix save_covar = GSparseMatrix(*covar);
  
-    // Loop over error computation (maximum 2 turns)
-    for (int i = 0; i < 2; ++i) {
+        // Signal no diagonal element loading
+        bool diag_loaded = false;
  
-        // Solve: covar * X = unit
-        try {
-            covar->cholesky_decompose(1);
-            GVector unit(npars);
-            for (int ipar = 0; ipar < npars; ++ipar) {
-                unit[ipar] = 1.0;
-                GVector x  = covar->cholesky_solver(unit,1);
-                if (x[ipar] >= 0.0)
-                    pars->par(ipar).error(sqrt(x[ipar]));
-                else {
-                    pars->par(ipar).error(0.0);
-                    m_status = G_LM_BAD_ERRORS;
+        // Loop over error computation (maximum 2 turns)
+        for (int i = 0; i < 2; ++i) {
+
+            // Solve: covar * X = unit
+            try {
+                covar->cholesky_decompose(1);
+                GVector unit(npars);
+                for (int ipar = 0; ipar < npars; ++ipar) {
+                    unit[ipar] = 1.0;
+                    GVector x  = covar->cholesky_solver(unit,1);
+                    if (x[ipar] >= 0.0) {
+                        pars->par(ipar).error(sqrt(x[ipar]));
+                    }
+                    else {
+                        pars->par(ipar).error(0.0);
+                        m_status = G_LM_BAD_ERRORS;
+                    }
+                    unit[ipar] = 0.0;
                 }
-                unit[ipar] = 0.0;
-            }
-        }
-        catch (GException::matrix_zero &e) {
-            m_status = G_LM_SINGULAR;
-            if (m_logger != NULL) {
-                *m_logger << "GOptimizerLM::terminate: "
-                          << "All curvature matrix elements are zero."
-                          << std::endl;
             }
-            break;
-        }
-        catch (GException::matrix_not_pos_definite &e) {
-
-            // Load diagonal if this has not yet been tried
-            if (!diag_loaded) {
-
-                // Flag errors as inaccurate
-                m_status = G_LM_BAD_ERRORS;
+            catch (GException::matrix_zero &e) {
+                m_status = G_LM_SINGULAR;
                 if (m_logger != NULL) {
-                    *m_logger << "Non-Positive definite curvature matrix encountered."
-                              << std::endl;
-                    *m_logger << "Load diagonal elements with 1e-10."
-                              << " Fit errors may be inaccurate."
+                    *m_logger << "GOptimizerLM::terminate: "
+                              << "All curvature matrix elements are zero."
                               << std::endl;
                 }
+                break;
+            }
+            catch (GException::matrix_not_pos_definite &e) {
  
-                // Try now with diagonal loaded matrix
-                *covar = save_covar;
-                for (int ipar = 0; ipar < npars; ++ipar)
-                    (*covar)(ipar,ipar) += 1.0e-10;
+                // Load diagonal if this has not yet been tried
+                if (!diag_loaded) {
  
-                // Signal loading
-                diag_loaded = true;
+                    // Flag errors as inaccurate
+                    m_status = G_LM_BAD_ERRORS;
+                    if (m_logger != NULL) {
+                        *m_logger << "Non-Positive definite curvature matrix encountered."
+                                  << std::endl;
+                        *m_logger << "Load diagonal elements with 1e-10."
+                                  << " Fit errors may be inaccurate."
+                                  << std::endl;
+                    }
  
-                // Try again
-                continue;
+                    // Try now with diagonal loaded matrix
+                    *covar = save_covar;
+                    for (int ipar = 0; ipar < npars; ++ipar) {
+                        (*covar)(ipar,ipar) += 1.0e-10;
+                    }
  
-            } // endif: diagonal has not yet been loaded
+                    // Signal loading
+                    diag_loaded = true;
  
-            // ... otherwise signal an error
-            else {
-                m_status = G_LM_NOT_POSTIVE_DEFINITE;
-                if (m_logger != NULL) {
-                    *m_logger << "Non-Positive definite curvature matrix encountered,"
-                              << " even after diagonal loading." << std::endl;
+                    // Try again
+                    continue;
+
+                } // endif: diagonal has not yet been loaded
+
+                // ... otherwise signal an error
+                else {
+                    m_status = G_LM_NOT_POSTIVE_DEFINITE;
+                    if (m_logger != NULL) {
+                        *m_logger << "Non-Positive definite curvature matrix encountered,"
+                                  << " even after diagonal loading." << std::endl;
+                    }
+                    break;
                 }
-                break;
             }
-        }
-        catch (std::exception &e) {
-            throw;
-        }
+            catch (std::exception &e) {
+                throw;
+            }
  
-        // If no error occured then break now
-        break;
+            // If no error occured then break now
+            break;
+
+        } // endfor: looped over error computation
  
-    } // endfor: looped over error computation
+    } // endif: pointers were valid
  
     // Return
     return;
 }
-
-
-/*==========================================================================
- =                                                                         =
- =                                Friends                                  =
- =                                                                         =
- ==========================================================================*/
@@ -255,7 +255,7 @@ GException::wcs_hpx_bad_ordering::wcs_hpx_bad_ordering(std::string origin,
  * @brief Singular matrix encountered
  *
  * @param[in] origin Method that throws the error.
- * @param[in] naxis Number of matrix axes.
+ * @param[in] naxis Matrix size (naxis*naxis).
  * @param[in] mat Matrix (stored in vector).
  ***************************************************************************/
 GException::wcs_singular_matrix::wcs_singular_matrix(std::string                origin,
@@ -266,13 +266,20 @@ GException::wcs_singular_matrix::wcs_singular_matrix(std::string                
     m_origin  = origin;
     m_message = "Singular matrix encountered: mat=[";
     for (int i = 0, ij = 0; i < naxis; ++i) {
-        if (i > 0)
+        if (i > 0) {
             m_message += ",";
+        }
         m_message += "[";
         for (int j = 0; j < naxis; ++j, ++ij) {
-            if (j > 0)
+            if (j > 0) {
                 m_message += ",";
-            m_message += str(mat[ij]);
+            }
+            if (ij < mat.size()) {
+                m_message += str(mat[ij]);
+            }
+            else {
+                m_message += " ... ";
+            }
         }
         m_message += "]";
     }
 /***************************************************************************
- *                 GNodeArray.cpp  -  Array of nodes class                 *
+ *                  GNodeArray.cpp - Array of nodes class                  *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2008-2012 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
@@ -35,7 +35,9 @@
  
 /* __ Method name definitions ____________________________________________ */
 #define G_ACCESS                                "GNodeArray::operator[](int)"
-#define G_INTERPOLATE "GNodeArray::interpolate(std::vector<double>&, double&"
+#define G_INTERPOLATE "GNodeArray::interpolate(double&,std::vector<double>&)"
+#define G_SET_VALUE                          "GNodeArray::set_value(double&)"
+#define G_SETUP                                         "GNodeArray::setup()"
  
 /* __ Macros _____________________________________________________________ */
  
@@ -166,7 +168,8 @@ GNodeArray::~GNodeArray(void)
 /***********************************************************************//**
  * @brief Assignment operator
  *
- * @param[in] array Instance that should be assigned.
+ * @param[in] array Node array.
+ * @return Node array.
  ***************************************************************************/
 GNodeArray& GNodeArray::operator= (const GNodeArray& array)
 {
@@ -193,6 +196,7 @@ GNodeArray&amp; GNodeArray::operator= (const GNodeArray&amp; array)
  * @brief Node access operator
  *
  * @param[in] index Node index [0,...,size()-1].
+ * @return Node value.
  *
  * @exception GException::out_of_range
  *            Node index is out of range.
@@ -201,8 +205,9 @@ double&amp; GNodeArray::operator[](const int&amp; index)
 {
     // Compile option: raise an exception if index is out of range
     #if defined(G_RANGE_CHECK)
-    if (index < 0 || index >= size())
+    if (index < 0 || index >= size()) {
         throw GException::out_of_range(G_ACCESS, index, size()-1);
+    }
     #endif
  
     // Return node
@@ -214,6 +219,7 @@ double&amp; GNodeArray::operator[](const int&amp; index)
  * @brief Node access operator (const version)
  *
  * @param[in] index Node index [0,...,size()-1].
+ * @return Node value.
  *
  * @exception GException::out_of_range
  *            Node index is out of range.
@@ -222,8 +228,9 @@ const double&amp; GNodeArray::operator[](const int&amp; index) const
 {
     // Compile option: raise an exception if index is out of range
     #if defined(G_RANGE_CHECK)
-    if (index < 0 || index >= size())
+    if (index < 0 || index >= size()) {
         throw GException::out_of_range(G_ACCESS, index, size()-1);
+    }
     #endif
  
     // Return node
@@ -239,7 +246,7 @@ const double&amp; GNodeArray::operator[](const int&amp; index) const
  
 /***********************************************************************//**
  * @brief Clear instance
-***************************************************************************/
+ ***************************************************************************/
 void GNodeArray::clear(void)
 {
     // Free class members
@@ -255,7 +262,9 @@ void GNodeArray::clear(void)
  
 /***********************************************************************//**
  * @brief Clone instance
-***************************************************************************/
+ *
+ * @return Deep copy of node array.
+ ***************************************************************************/
 GNodeArray* GNodeArray::clone(void) const
 {
     return new GNodeArray(*this);
@@ -267,6 +276,8 @@ GNodeArray* GNodeArray::clone(void) const
  *
  * @param[in] num Number of nodes
  * @param[in] array Node values \f$x_i\f$.
+ *
+ * Setup node array from an array of double precision values.
  ***************************************************************************/
 void GNodeArray::nodes(const int& num, const double* array)
 {
@@ -277,8 +288,9 @@ void GNodeArray::nodes(const int&amp; num, const double* array)
     init_members();
  
     // Set node values
-    for (int i = 0; i < num; ++i)
+    for (int i = 0; i < num; ++i) {
         m_node.push_back(array[i]);
+    }
  
     // Setup node distances and linear array handling
     setup();
@@ -292,6 +304,8 @@ void GNodeArray::nodes(const int&amp; num, const double* array)
  * @brief Set node array from vector
  *
  * @param[in] vector Vector from which node array will be built.
+ *
+ * Setup node array from a vector of values.
  ***************************************************************************/
 void GNodeArray::nodes(const GVector& vector)
 {
@@ -302,8 +316,9 @@ void GNodeArray::nodes(const GVector&amp; vector)
     init_members();
  
     // Set node values
-    for (int i = 0; i < vector.size(); ++i)
+    for (int i = 0; i < vector.size(); ++i) {
         m_node.push_back(vector[i]);
+    }
  
     // Setup node distances and linear array handling
     setup();
@@ -317,6 +332,8 @@ void GNodeArray::nodes(const GVector&amp; vector)
  * @brief Set node array from vector
  *
  * @param[in] vector Vector from which node array will be built.
+ *
+ * Setup node array from a vector of double precision values.
  ***************************************************************************/
 void GNodeArray::nodes(const std::vector<double>& vector)
 {
@@ -373,13 +390,15 @@ double GNodeArray::interpolate(const double&amp; value,
                                const std::vector<double>& vector) const
 {
     // Throw exception if there are not enough nodes
-    if (m_node.size() < 2)
+    if (m_node.size() < 2) {
         throw GException::not_enough_nodes(G_INTERPOLATE, m_node.size());
+    }
  
     // Throw exception if vectors have not the same size
-    if (m_node.size() != vector.size())
+    if (m_node.size() != vector.size()) {
         throw GException::vector_mismatch(G_INTERPOLATE, m_node.size(),
                                           vector.size());
+    }
  
     // Set interpolation value (circumvent const correctness)
     const_cast<GNodeArray*>(this)->set_value(value);
@@ -398,6 +417,9 @@ double GNodeArray::interpolate(const double&amp; value,
  *
  * @param[in] value Value for which the interpolation should be done.
  *
+ * @exception GException::not_enough_nodes
+ *            At least two nodes are required for setting up the factors
+ *
  * Set the indices that bound the specified value and the corresponding
  * weighting factors for linear interpolation. If the array has a linear
  * form (i.e. the nodes are equidistant), an analytic formula is used to
@@ -408,55 +430,60 @@ void GNodeArray::set_value(const double&amp; value)
 {
     // Get number of nodes
     int nodes = m_node.size();
+
+    // Throw an exception if less than 2 nodes are available
+    if (nodes < 2) {
+        throw GException::not_enough_nodes(G_SET_VALUE, nodes);
+    }
  
-    // Continue only if we have at least 2 nodes
-    if (nodes > 1) {
-    
-        // If array is linear then get left index from analytic formula
-        if (m_is_linear) {
+    // If array is linear then get left index from analytic formula
+    if (m_is_linear) {
  
-            // Set left index
-            m_inx_left = int(m_linear_slope * value + m_linear_offset);
+        // Set left index
+        m_inx_left = int(m_linear_slope * value + m_linear_offset);
+        
+        // Keep index in valid range
+        if (m_inx_left < 0)             m_inx_left = 0;
+        else if (m_inx_left >= nodes-1) m_inx_left = nodes - 2;
  
-            // Keep index in valid range
-            if (m_inx_left < 0)             m_inx_left = 0;
-            else if (m_inx_left >= nodes-1) m_inx_left = nodes - 2;
+    } // endif: array is linear
  
-        } // endif: array is linear
+    // ... otherwise search the relevant indices by bisection
+    else {
+    
+        // Set left index if value is before first node
+        if (value < m_node[0]) {
+            m_inx_left = 0;
+        }
  
-        // ... otherwise search the relevant indices by bisection
+        // Set left index if value is after last node
+        else if (value >  m_node[nodes-1]) {
+            m_inx_left = nodes - 2;
+        }
+
+        // Set left index by bisection
         else {
-            // Set left index if value is before first node
-            if (value < m_node[0])
-                m_inx_left = 0;
-
-            // Set left index if value is after last node
-            else if (value >  m_node[nodes-1])
-                m_inx_left = nodes - 2;
-
-            // Set left index by bisection
-            else {
-                int low  = 0;
-                int high = nodes - 1;
-                while ((high - low) > 1) {
-                    int mid = (low+high) / 2;
-                    if (m_node[mid] > value)
-                        high = mid;
-                    else if (m_node[mid] <= value)
-                        low = mid;
+            int low  = 0;
+            int high = nodes - 1;
+            while ((high - low) > 1) {
+                int mid = (low+high) / 2;
+                if (m_node[mid] > value) {
+                    high = mid;
                 }
-                m_inx_left = low;
-            } // endelse: did bisection
-        }
+                else if (m_node[mid] <= value) {
+                    low = mid;
+                }
+            }
+            m_inx_left = low;
+        } // endelse: did bisection
+    }
  
-        // Set right index
-        m_inx_right = m_inx_left + 1;
+    // Set right index
+    m_inx_right = m_inx_left + 1;
  
-        // Set weighting factors
-        m_wgt_right = (value - m_node[m_inx_left]) / m_step[m_inx_left];
-        m_wgt_left  = 1.0 - m_wgt_right;
-    
-    } // endif: there were at least 2 nodes
+    // Set weighting factors
+    m_wgt_right = (value - m_node[m_inx_left]) / m_step[m_inx_left];
+    m_wgt_left  = 1.0 - m_wgt_right;
  
     // Return
     return;
@@ -493,7 +520,9 @@ std::string GNodeArray::print(void) const
     for (int i = 0; i < size(); ++i) {
         result.append("\n"+parformat("Node "+str(i)));
         result.append(str(m_node[i]));
-        result.append(" (delta="+str(m_step[i])+")");
+        if (i < m_step.size()) {
+            result.append(" (delta="+str(m_step[i])+")");
+        }
     }
  
     // Return result
@@ -515,7 +544,7 @@ void GNodeArray::init_members(void)
     // Initialise members
     m_node.clear();
     m_step.clear();
-    m_is_linear     = 0;
+    m_is_linear     = false;
     m_linear_slope  = 0.0;
     m_linear_offset = 0.0;
     m_inx_left      = 0;
@@ -563,6 +592,11 @@ void GNodeArray::free_members(void)
  
 /***********************************************************************//**
  * @brief Compute distance array and linear slope/offset
+ *
+ * Precomputes values for fast interpolation. The precomputation requires
+ * at least 2 nodes to be present in the node array. If less then two
+ * nodes are present, the distance vector m_step will be empty and no
+ * computation is done.
  ***************************************************************************/
 void GNodeArray::setup(void)
 {
@@ -571,13 +605,14 @@ void GNodeArray::setup(void)
  
     // Get number of nodes
     int nodes = m_node.size();
-    
-    // Continue only if we have nodes at least 2 nodes
+
+    // Continue only if there are at least 2 nodes
     if (nodes > 1) {
-    
+
         // Setup distance array between subsequent nodes
-        for (int i = 0; i < nodes-1; ++i)
+        for (int i = 0; i < nodes-1; ++i) {
             m_step.push_back(m_node[i+1] - m_node[i]);
+        }
  
         // Evaluate linear slope and offset
         m_linear_slope  = double(nodes-1) / (m_node[nodes-1] - m_node[0]);
@@ -592,9 +627,10 @@ void GNodeArray::setup(void)
                 break;
             }
         }
-        
-    } // endif: there were at least 2 nodes
+
+    } // endif: there were at least two nodes
  
     // Return
     return;
 }
+
@@ -258,8 +258,8 @@ void TestGSky::test_GWcslib(void)
         try {
  
             // Allocate projection CEL and GAL
-            GWcslib *cel = (GWcslib*)registry.alloc(registry.code(i));
-            GWcslib *gal = (GWcslib*)registry.alloc(registry.code(i));
+            GWcslib *cel = dynamic_cast<GWcslib*>(registry.alloc(registry.code(i)));
+            GWcslib *gal = dynamic_cast<GWcslib*>(registry.alloc(registry.code(i)));
  
             // Throw an error if allocation failed
             if (cel == NULL || gal == NULL) {