Merge branch 'integration' into devel

Jürgen Knödlseder
2 parents 42502f5c c3f23172
Showing 36 changed files with 195 additions and 115 deletions
include/GApplicationPars.hpp
include/GArf.hpp
include/GEnergies.hpp
include/GFits.hpp
include/GFitsHeader.hpp
include/GModel.hpp
include/GModelSpatial.hpp
include/GModelSpectral.hpp
include/GModelSpectralNodes.hpp
include/GModelTemporal.hpp
include/GModels.hpp
include/GNodeArray.hpp
include/GObservations.hpp
include/GOptimizerPars.hpp
include/GPha.hpp
include/GPhotons.hpp
include/GSkyRegions.hpp
include/GTestSuite.hpp
include/GTestSuites.hpp
include/GTimes.hpp
 /***************************************************************************
  *               GApplicationPars.hpp - Application parameters             *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2014 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2015 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -160,7 +160,7 @@ const GApplicationPar&amp; GApplicationPars::operator[](const int&amp; index) const
 inline
 int GApplicationPars::size(void) const
 {
-    return (m_pars.size());
+    return (int)m_pars.size();
 }
  
  
 /***************************************************************************
  *              GArf.hpp - XSPEC Auxiliary Response File class             *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2013-2014 by Juergen Knoedlseder                         *
+ *  copyright (C) 2013-2015 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -134,7 +134,7 @@ const double&amp; GArf::operator[](const int&amp; index) const
 inline
 int GArf::size(void) const
 {
-    return (m_specresp.size());
+    return (int)m_specresp.size();
 }
  
  
@@ -149,7 +149,7 @@ const GEnergy&amp; GEnergies::operator[](const int&amp; index) const
 inline
 int GEnergies::size(void) const
 {
-    return (m_energies.size());
+    return (int)m_energies.size();
 }
  
  
 /***************************************************************************
  *                       GFits.hpp - FITS file class                       *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2008-2014 by Juergen Knoedlseder                         *
+ *  copyright (C) 2008-2015 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -231,7 +231,7 @@ const GFitsHDU* GFits::operator[](const std::string&amp; extname) const
 inline
 int GFits::size(void) const
 {
-    return (m_hdu.size());
+    return (int)m_hdu.size();
 }
  
  
 /***************************************************************************
  *            GFitsHeader.hpp - FITS header cards container class          *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2008-2014 by Juergen Knoedlseder                         *
+ *  copyright (C) 2008-2015 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -175,7 +175,7 @@ const GFitsHeaderCard&amp; GFitsHeader::operator[](const std::string&amp; keyname) const
 inline
 int GFitsHeader::size(void) const
 {
-    return (m_cards.size());
+    return (int)m_cards.size();
 }
  
  
 /***************************************************************************
  *              GModel.hpp - Abstract virtual model base class             *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2009-2014 by Juergen Knoedlseder                         *
+ *  copyright (C) 2009-2015 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -212,7 +212,7 @@ const GModelPar&amp; GModel::operator[](const int&amp; index) const
 inline
 int GModel::size(void) const
 {
-    return (m_pars.size());
+    return (int)m_pars.size();
 }
  
  
@@ -157,7 +157,7 @@ const GModelPar&amp; GModelSpatial::operator[](const int&amp; index) const
 inline
 int GModelSpatial::size(void) const
 {
-    return (m_pars.size());
+    return (int)m_pars.size();
 }
  
 #endif /* GMODELSPATIAL_HPP */
 /***************************************************************************
  *         GModelSpectral.hpp - Abstract spectral model base class         *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2009-2014 by Juergen Knoedlseder                         *
+ *  copyright (C) 2009-2015 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -155,7 +155,7 @@ const GModelPar&amp; GModelSpectral::operator[](const int&amp; index) const
 inline
 int GModelSpectral::size(void) const
 {
-    return (m_pars.size());
+    return (int)m_pars.size();
 }
  
 #endif /* GMODELSPECTRAL_HPP */
 /***************************************************************************
  *          GModelSpectralNodes.hpp - Spectral nodes model class           *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012-2014 by Juergen Knoedlseder                         *
+ *  copyright (C) 2012-2015 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -201,7 +201,7 @@ std::string GModelSpectralNodes::type(void) const
 inline
 int GModelSpectralNodes::nodes(void) const
 {
-    return (m_energies.size());
+    return (int)m_energies.size();
 }
  
 #endif /* GMODELSPECTRALNODES_HPP */
@@ -134,7 +134,7 @@ const GModelPar&amp; GModelTemporal::operator[](const int&amp; index) const
 inline
 int GModelTemporal::size(void) const
 {
-    return (m_pars.size());
+    return (int)m_pars.size();
 }
  
 #endif /* GMODELTEMPORAL_HPP */
 /***************************************************************************
  *                     GModels.hpp - Model container class                 *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2009-2014 by Juergen Knoedlseder                         *
+ *  copyright (C) 2009-2015 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -248,7 +248,7 @@ const GModel* GModels::operator[](const int&amp; index) const
 inline
 int GModels::size(void) const
 {
-    return (m_models.size());
+    return (int)m_models.size();
 }
  
  
@@ -184,7 +184,7 @@ const double&amp; GNodeArray::operator[](const int&amp; index) const
 inline
 int GNodeArray::size(void) const
 {
-    return (m_node.size());
+    return (int)m_node.size();
 }
  
  
@@ -238,7 +238,7 @@ const GObservation* GObservations::operator[](const int&amp; index) const
 inline
 int GObservations::size(void) const
 {
-    return (m_obs.size());
+    return (int)m_obs.size();
 }
  
  
 /***************************************************************************
  *         GOptimizerPars.hpp - Optimizer parameter container class        *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2009-2014 by Juergen Knoedlseder                         *
+ *  copyright (C) 2009-2015 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -152,7 +152,7 @@ const GOptimizerPar* GOptimizerPars::operator[](const int&amp; index) const
 inline
 int GOptimizerPars::size(void) const
 {
-    return (m_pars.size());
+    return (int)m_pars.size();
 }
  
  
 /***************************************************************************
  *                GPha.hpp - XSPEC Pulse Height Analyzer class             *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2013-2014 by Juergen Knoedlseder                         *
+ *  copyright (C) 2013-2015 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -154,7 +154,7 @@ const double&amp; GPha::operator[](const int&amp; index) const
 inline
 int GPha::size(void) const
 {
-    return (m_counts.size());
+    return (int)m_counts.size();
 }
  
  
 /***************************************************************************
  *                   GPhotons.hpp - Photon container class                 *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012-2014 by Juergen Knoedlseder                         *
+ *  copyright (C) 2012-2015 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -99,7 +99,7 @@ std::string GPhotons::classname(void) const
 inline
 int GPhotons::size(void) const
 {
-    return m_photons.size();
+    return (int)m_photons.size();
 }
  
  
 /***************************************************************************
  *                  GSkyRegions.hpp - Sky region container class           *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2013-2014 by Pierrick Martin                             *
+ *  copyright (C) 2013-2015 by Pierrick Martin                             *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -154,7 +154,7 @@ const GSkyRegion* GSkyRegions::operator[](const int&amp; index) const
 inline
 int GSkyRegions::size(void) const
 {
-    return (m_regions.size());
+    return (int)m_regions.size();
 }
  
  
 /***************************************************************************
  *             GTestSuite.hpp - Abstract test suite base class             *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012-2014 by Jean-Baptiste Cayrou                        *
+ *  copyright (C) 2012-2015 by Jean-Baptiste Cayrou                        *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -129,7 +129,7 @@ protected:
 inline
 int GTestSuite::size(void) const
 {
-    return m_tests.size();
+    return (int)m_tests.size();
 }
  
  
 /***************************************************************************
  *              GTestSuites.hpp - Test suite container class               *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012-2014 by Jean-Baptiste Cayrou                        *
+ *  copyright (C) 2012-2015 by Jean-Baptiste Cayrou                        *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -155,7 +155,7 @@ const GTestSuite* GTestSuites::operator[](const int&amp; index) const
 inline
 int GTestSuites::size(void) const
 {
-    return m_testsuites.size();
+    return (int)m_testsuites.size();
 }
  
  
 /***************************************************************************
  *                     GTimes.hpp - Time container class                   *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012-2014 by Juergen Knoedlseder                         *
+ *  copyright (C) 2012-2015 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -99,7 +99,7 @@ std::string GTimes::classname(void) const
 inline
 int GTimes::size(void) const
 {
-    return m_times.size();
+    return (int)m_times.size();
 }
  
  
@@ -138,7 +138,7 @@ protected:
     std::string         m_hub_path;  //!< Hub path
     std::string         m_version;   //!< The version of the SAMP Standard Profile implemented by the hub
     std::string         m_hub_id;    //!< Hub identifier used by the hub when it sends message itself rather than forwarding from others
-    socklen_t           m_socket;    //!< Hub socket
+    int                 m_socket;    //!< Hub socket
     bool                m_shutdown;  //!< Shutdown request
     std::vector<client> m_clients;   //!< Clients
 };
 /***************************************************************************
  *                GXmlNode.hpp - Abstract XML node base class              *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2014 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2015 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -125,7 +125,7 @@ protected:
 inline
 int GXmlNode::size(void) const
 {
-    return m_nodes.size();
+    return (int)m_nodes.size();
 }
  
  
@@ -128,7 +128,7 @@ std::string GCTAEventList::classname(void) const
 inline
 int GCTAEventList::size(void) const
 {
-    return (m_events.size());
+    return (int)m_events.size();
 }
  
  
@@ -140,7 +140,7 @@ int GCTAEventList::size(void) const
 inline
 int GCTAEventList::number(void) const
 {
-    return (m_events.size());
+    return (int)m_events.size();
 }
  
  
 /***************************************************************************
  *          GCTAModelRadial.hpp - Radial model abstract base class         *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2011-2014 by Juergen Knoedlseder                         *
+ *  copyright (C) 2011-2015 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -73,7 +73,7 @@ public:
     virtual std::string      print(const GChatter& chatter = NORMAL) const = 0;
  
     // Methods
-    int size(void) const { return m_pars.size(); }
+    int size(void) const { return (int)m_pars.size(); }
  
 protected:
     // Protected methods
 /***************************************************************************
  *       GCTAModelRadialPolynom.hpp - Radial Polynom CTA model class       *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2011-2014 by Juergen Knoedlseder                         *
+ *  copyright (C) 2011-2015 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -139,7 +139,7 @@ std::string GCTAModelRadialPolynom::type(void) const
 inline
 int GCTAModelRadialPolynom::size(void) const
 {
-    return (m_coeffs.size());
+    return (int)m_coeffs.size();
 }
  
 #endif /* GCTAMODELRADIALPOLYNOM_HPP */
 /***************************************************************************
  *    GCTAOnOffObservations.hpp - CTA on-off observation container class   *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2013-2014 by Pierrick Martin                             *
+ *  copyright (C) 2013-2015 by Pierrick Martin                             *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -158,7 +158,7 @@ const GCTAOnOffObservation* GCTAOnOffObservations::operator[](const int&amp; index) 
 inline
 int GCTAOnOffObservations::size(void) const
 {
-    return (m_obs.size());
+    return (int)m_obs.size();
 }
  
  
@@ -864,11 +864,15 @@ double cta_irf_elliptical_kern_omega::eval(const double&amp; omega)
         double arg2            = m_semimajor_rad * sinus;
         double r_ellipse       = m_semiminor_rad * m_semimajor_rad /
                                  std::sqrt(arg1*arg1 + arg2*arg2);
-        std::cout << "*** WARNING: cta_irf_elliptical_kern_omega::eval";
+        std::cout << "*** WARNING: cta_irf_elliptical_kern_omega::eval:";
         std::cout << " zero model for (rho,omega)=(";
         std::cout << m_rho*gammalib::rad2deg << ",";
         std::cout << omega*gammalib::rad2deg << ")";
-        std::cout << " rho-r_ellipse=" << (m_rho-r_ellipse) << " radians";
+        std::cout << " semiminor=" << m_semiminor_rad;
+        std::cout << " semimajor=" << m_semimajor_rad;
+        std::cout << " posang=" << m_model->posangle() * gammalib::deg2rad;
+        std::cout << " rel_posang=" << diff_angle * gammalib::deg2rad;
+        std::cout << " rho/r_ellipse=" << (m_rho/r_ellipse);
         std::cout << std::endl;
     }
     #endif
@@ -278,7 +278,7 @@ int GLATEventCube::ebins(void) const
 inline
 int GLATEventCube::ndiffrsp(void) const
 {
-    return m_srcmap.size();
+    return (int)m_srcmap.size();
 }
  
 #endif /* GLATEVENTCUBE_HPP */
@@ -115,7 +115,7 @@ std::string GLATEventList::classname(void) const
 inline
 int GLATEventList::size(void) const
 {
-    return m_events.size();
+    return (int)m_events.size();
 }
  
  
@@ -127,7 +127,7 @@ int GLATEventList::size(void) const
 inline
 int GLATEventList::number(void) const
 {
-    return m_events.size();
+    return (int)m_events.size();
 }
  
  
@@ -164,7 +164,7 @@ bool GLATResponse::use_tdisp(void) const
 inline
 int GLATResponse::size(void) const
 {
-    return m_aeff.size();
+    return (int)m_aeff.size();
 }
  
  
 /***************************************************************************
  *          GMWLSpectrum.hpp  -  Multi-wavelength spectrum class           *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2014 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2015 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -119,7 +119,7 @@ std::string GMWLSpectrum::classname(void) const
 inline
 int GMWLSpectrum::size(void) const
 {
-    return m_data.size();
+    return (int)m_data.size();
 }
  
  
@@ -144,7 +144,7 @@ int GMWLSpectrum::dim(void) const
 inline
 int GMWLSpectrum::naxis(const int& axis) const
 {
-    return m_data.size();
+    return (int)m_data.size();
 }
  
  
@@ -138,8 +138,8 @@ public:
     GSkyMap       operator-(const GSkyMap& map) const;
     GSkyMap       operator*(const GSkyMap& map) const;
     GSkyMap       operator/(const GSkyMap& map) const;
-    double&       operator()(const int& index, const int& map = 0);
-    double&       operator()(const GSkyPixel& pixel, const int& map = 0);
+    double        operator()(const int& index, const int& map = 0);
+    double        operator()(const GSkyPixel& pixel, const int& map = 0);
     double        operator()(const GSkyDir& dir, const int& map = 0) const;
  
     // Methods
@@ -36,6 +36,11 @@
 #include "GModelSpatialRegistry.hpp"
  
 /* __ Constants __________________________________________________________ */
+namespace {
+    const double c_theta_max    = 3.0; //!< semiaxis multiplied for theta_max()
+    const double c_max_exponent = 0.5 * c_theta_max * c_theta_max;
+    const double c_fraction     = 1.0 - std::exp(-c_max_exponent);
+}
  
 /* __ Globals ____________________________________________________________ */
 const GModelSpatialEllipticalGauss g_elliptical_gauss_seed;
@@ -48,6 +53,7 @@ const GModelSpatialRegistry g_elliptical_gauss_registry(&amp;g_elliptical_gauss_seed
 /* __ Macros _____________________________________________________________ */
  
 /* __ Coding definitions _________________________________________________ */
+#define G_SMALL_ANGLE_APPROXIMATION      //!< Use small angle approximation
  
 /* __ Debug definitions __________________________________________________ */
  
@@ -73,14 +79,14 @@ GModelSpatialEllipticalGauss::GModelSpatialEllipticalGauss(void) :
  
  
 /***********************************************************************//**
- * @brief Elliptical Gauss constructor
+ * @brief Elliptical Gaussian constructor
  *
- * @param[in] dir Sky position of gauss centre.
+ * @param[in] dir Centre of elliptical Gaussian.
  * @param[in] semimajor Semi-major axis (degrees).
  * @param[in] semiminor Semi-minor axis (degrees).
  * @param[in] posangle Position angle of semi-major axis (degrees).
  *
- * Construct elliptical gauss model from sky position of the ellipse centre
+ * Construct elliptical Gaussian model from the ellipse centre direction
  * (@p dir), the @p semimajor and @p semiminor axes, and the position
  * angle (@p posangle).
  ***************************************************************************/
@@ -109,7 +115,7 @@ GModelSpatialEllipticalGauss::GModelSpatialEllipticalGauss(const GSkyDir&amp; dir,
  *
  * @param[in] xml XML element.
  *
- * Constructs elliptical gauss model by extracting information from an XML
+ * Constructs elliptical Gaussian model by extracting information from an XML
  * element. See the read() method for more information about the expected
  * structure of the XML element.
  ***************************************************************************/
@@ -130,7 +136,7 @@ GModelSpatialEllipticalGauss::GModelSpatialEllipticalGauss(const GXmlElement&amp; xm
 /***********************************************************************//**
  * @brief Copy constructor
  *
- * @param[in] model Elliptical gauss model.
+ * @param[in] model Elliptical Gaussian model.
  ***************************************************************************/
 GModelSpatialEllipticalGauss::GModelSpatialEllipticalGauss(const GModelSpatialEllipticalGauss& model) :
                               GModelSpatialElliptical(model)
@@ -169,7 +175,7 @@ GModelSpatialEllipticalGauss::~GModelSpatialEllipticalGauss(void)
  * @brief Assignment operator
  *
  * @param[in] model Elliptical gauss model.
- * @return Elliptical gauss model.
+ * @return Elliptical Gaussian model.
  ***************************************************************************/
 GModelSpatialEllipticalGauss& GModelSpatialEllipticalGauss::operator=(const GModelSpatialEllipticalGauss& model)
 {
@@ -202,7 +208,7 @@ GModelSpatialEllipticalGauss&amp; GModelSpatialEllipticalGauss::operator=(const GMod
  ==========================================================================*/
  
 /***********************************************************************//**
- * @brief Clear elliptical gauss model
+ * @brief Clear elliptical Gaussian model
  ***************************************************************************/
 void GModelSpatialEllipticalGauss::clear(void)
 {
@@ -222,13 +228,13 @@ void GModelSpatialEllipticalGauss::clear(void)
  
  
 /***********************************************************************//**
- * @brief Clone elliptical gauss model
+ * @brief Clone elliptical Gaussian model
  *
- * @return Pointer to deep copy of elliptical gauss model.
+ * @return Pointer to deep copy of elliptical Gaussian model.
  ***************************************************************************/
 GModelSpatialEllipticalGauss* GModelSpatialEllipticalGauss::clone(void) const
 {
-    // Clone elliptical gauss model
+    // Clone elliptical Gaussian model
     return new GModelSpatialEllipticalGauss(*this);
 }
  
@@ -242,44 +248,50 @@ GModelSpatialEllipticalGauss* GModelSpatialEllipticalGauss::clone(void) const
  * @param[in] time Photon arrival time.
  * @return Model value.
  *
- * Evaluates the spatial component for an elliptical gauss source model. The
- * source model is an elliptical gaussian function
+ * Evaluates the spatial component for an elliptical Gaussian source model.
+ *
+ * The elliptical Gaussian function is defined by
  * \f$S_{\rm p}(\theta, \phi | E, t)\f$, where
- * \f$\theta\f$ is the angular separation between elliptical gauss centre and
- * the actual location and \f$\phi\f$ the position angle with respect to the
- * model centre, counted counterclockwise from North.
+ * \f$\theta\f$ is the angular separation between elliptical Gaussian centre
+ * and the actual location and \f$\phi\f$ the position angle with respect to
+ * the model centre, counted counterclockwise from North.
  *
- * The function \f$f(\theta, \phi)\f$ is given by
+ * The function \f$S_{\rm p}(\theta, \phi | E, t)\f$ is given by
  *
  * \f[
- * S_{\rm p}(\theta, \phi | E, t) = \left \{
- *  \begin{array}{l l}
- *     {\tt m\_norm}
- *     & \mbox{if} \, \, \theta \le \theta_0 \\
- *     \\
- *     0 & \mbox{else}
- *  \end{array}
- *  \right .
+ * S_{\rm p}(\theta, \phi | E, t) = {\tt m\_norm} \times
+ * \exp \left( -\frac{\theta^2}{2 r_{rm eff}} \right)
  * \f]
  *
- * where \f$\theta_0\f$ is the effective radius of the ellipse on the sphere
- * given by
- *
- * \f[\theta_0\ =
- *    \frac{ab}{\sqrt{b^2 \cos^2(\phi-\phi_0) + a^2 \sin^2(\phi-\phi_0)}}\f]
+ * where the effective ellipse radius \f$r_{rm eff}\f$ towards a given
+ * position angle is given by
  *
- * and
+ * \f[
+ * r_{rm eff} = \frac{ab}
+ *                   {\sqrt{\left( a \sin (\phi - \phi_0) \right)^2 +
+ *                    \sqrt{\left( b \cos (\phi - \phi_0) \right)^2}
+ * \f]
+ * and  
  * \f$a\f$ is the semi-major axis of the ellipse,
  * \f$b\f$ is the semi-minor axis, and
  * \f$\phi_0\f$ is the position angle of the ellipse, counted
  * counterclockwise from North.
+ * \f${\tt m\_norm}\f$ is a normalization constant given by
  *
- * The normalisation constant \f${\tt m\_norm}\f$ which is the inverse of the
- * solid angle subtended by an ellipse is given by
+ * \f[
+ *    {\tt m\_norm} = \frac{1}{2 \pi \times a \times b}
+ * \f]
  *
- * @todo Quote formula for ellipse solid angle
+ * @warning
+ * The above formula for an elliptical Gaussian are accurate for small
+ * angles, with semimajor and semiminor axes below a few degrees. This
+ * should be fine for most use cases, but you should be aware that the
+ * ellipse will get distorted for larger angles.
  *
- * (see the update() method).
+ * @warning
+ * For numerical reasons the elliptical Gaussian will be truncated for
+ * \f$\theta\f$ angles that correspond to 3.0 times the effective ellipse
+ * radius.
  ***************************************************************************/
 double GModelSpatialEllipticalGauss::eval(const double&  theta,
                                           const double&  posangle,
@@ -295,6 +307,18 @@ double GModelSpatialEllipticalGauss::eval(const double&amp;  theta,
         // Update precomputation cache
         update();
  
+        // Perform computations to compute exponent
+        #if defined(G_SMALL_ANGLE_APPROXIMATION)
+        double rel_posangle  = posangle - this->posangle() * gammalib::deg2rad;
+        double cosinus       = std::cos(rel_posangle);
+        double sinus         = std::sin(rel_posangle);
+        double arg1          = m_minor_rad * cosinus;
+        double arg2          = m_major_rad * sinus;
+        double r_ellipse     = m_minor_rad * m_major_rad /
+                               std::sqrt(arg1*arg1 + arg2*arg2); //!< small angle
+        double r_relative    = theta/r_ellipse;
+        double exponent      = 0.5*r_relative*r_relative;
+        #else
         // Perform computations
         double sinphi = std::sin(posangle);
         double cosphi = std::cos(posangle);
@@ -305,10 +329,13 @@ double GModelSpatialEllipticalGauss::eval(const double&amp;  theta,
         double term3 = m_term3 * sinphi * cosphi;
  
         // Compute exponent
-        double exponent = theta * theta * (term1 + term2 + term3);
+        double exponent = theta * theta * (term1 + term2 + term3);        
+        #endif
  
-        // Set value
-        value = m_norm * std::exp(-exponent);
+        // Set value if the exponent is within the truncation region
+        if (exponent <= c_max_exponent) {
+            value = m_norm * std::exp(-exponent);
+        }
  
         // Compile option: Check for NaN/Inf
         #if defined(G_NAN_CHECK)
@@ -341,7 +368,7 @@ double GModelSpatialEllipticalGauss::eval(const double&amp;  theta,
  * @return Model value.
  *
  * Evaluates the function value. No gradient computation is implemented as
- * Elliptical models will be convolved with the instrument response and thus
+ * elliptical models will be convolved with the instrument response and thus
  * require the numerical computation of the derivatives.
  *
  * See the eval() method for more information.
@@ -364,7 +391,12 @@ double GModelSpatialEllipticalGauss::eval_gradients(const double&amp;  theta,
  * @param[in,out] ran Random number generator.
  * @return Sky direction.
  *
- * Draws an arbitrary sky position from the 2D gauss distribution.
+ * Draws an arbitrary sky direction from the elliptical Gaussian model.
+ *
+ * @warning
+ * For numerical reasons the elliptical Gaussian will be truncated for
+ * \f$\theta\f$ angles that correspond to 3.0 times the effective ellipse
+ * radius.
  ***************************************************************************/
 GSkyDir GModelSpatialEllipticalGauss::mc(const GEnergy& energy,
                                          const GTime&   time,
@@ -379,10 +411,18 @@ GSkyDir GModelSpatialEllipticalGauss::mc(const GEnergy&amp; energy,
     photon.time(time);
  
     // Draw gaussian offset from each axis
-    double theta1 = semimajor() * ran.normal();
-    double theta2 = semiminor() * ran.normal();
-
-    // Compute total offset from model centre
+    double ran_major;
+    double ran_minor;
+    do {
+        ran_major = ran.normal();
+    } while (ran_major > c_theta_max);
+    do {
+        ran_minor = ran.normal();
+    } while (ran_minor > c_theta_max);
+    double theta1 = semimajor() * ran_major;
+    double theta2 = semiminor() * ran_minor;
+
+    // Compute total offset from model centre in small angle approximation
     double theta = std::sqrt(theta1 * theta1 + theta2 * theta2);
  
     // Compute rotation angle, taking into account given position angle
@@ -401,7 +441,7 @@ GSkyDir GModelSpatialEllipticalGauss::mc(const GEnergy&amp; energy,
  
  
 /***********************************************************************//**
- * @brief Checks where model contains specified sky direction
+ * @brief Checks whether model contains specified sky direction
  *
  * @param[in] dir Sky direction.
  * @param[in] margin Margin to be added to sky direction (degrees)
@@ -428,16 +468,17 @@ bool GModelSpatialEllipticalGauss::contains(const GSkyDir&amp; dir,
  *
  * @return Returns maximum model radius.
  *
- * Returns \f$5 \sigma\f$ as approximate edge of the Gaussian. This limit
- * is of course arbitrary, but allows to limit the integration region for
- * response computation.
+ * Returns the maximum of 3.0 semimajor() and 3.0 semiminor() as
+ * approximate edge of the Gaussian. This limit is of course arbitrary, but
+ * allows to limit the integration region for response computation. The value
+ * of 3.0 has been determined by experiment (#1561).
  ***************************************************************************/
 double GModelSpatialEllipticalGauss::theta_max(void) const
 {
     // Set maximum model radius
     double theta_max = (semimajor() > semiminor())
-                       ? semimajor() * gammalib::deg2rad * 5.0
-                       : semiminor() * gammalib::deg2rad * 5.0;
+                       ? semimajor() * gammalib::deg2rad * c_theta_max
+                       : semiminor() * gammalib::deg2rad * c_theta_max;
  
     // Return value
     return theta_max;
@@ -684,7 +725,7 @@ void GModelSpatialEllipticalGauss::init_members(void)
 /***********************************************************************//**
  * @brief Copy class members
  *
- * @param[in] model Elliptical gauss model.
+ * @param[in] model Elliptical Gaussian model.
  ***************************************************************************/
 void GModelSpatialEllipticalGauss::copy_members(const GModelSpatialEllipticalGauss& model)
 {
@@ -724,10 +765,15 @@ void GModelSpatialEllipticalGauss::free_members(void)
  *
  * Computes the normalization
  * \f[
- *    {\tt m\_norm} = \frac{1}{2 \pi (1 - \cos a) (1 - \cos b)}
+ *    {\tt m\_norm} = \frac{1}{2 \pi \times a \times b}
  * \f]
+ * where 
+ * \f$a\f$ is the semi-major axis and
+ * \f$b\f$ is the semi-minor axis of the ellipse.
  *
- * @todo check this formula
+ * @warning
+ * The normalization of the elliptical Gaussian is only valid in the small
+ * angle approximation.
  ***************************************************************************/
 void GModelSpatialEllipticalGauss::update() const
 {
@@ -748,15 +794,20 @@ void GModelSpatialEllipticalGauss::update() const
         m_minor_rad = m_last_minor * gammalib::deg2rad;
         m_major_rad = m_last_major * gammalib::deg2rad;
  
-        // Perform precomputations
+        // Perform precomputations. Note that I verified the normalization
+        // by numerical integration of the resulting Gaussian. Note also
+        // also that the normalization is only correct in the small angle
+        // approximation.
         m_minor2     = m_minor_rad * m_minor_rad;
         m_major2     = m_major_rad * m_major_rad;
-        double denom = gammalib::twopi * m_minor_rad * m_major_rad;
+        double denom = gammalib::twopi * m_minor_rad * m_major_rad *
+                       c_fraction;
         m_norm       = (denom > 0.0) ? 1.0 / denom : 0.0;
  
     } // endif: update required
  
     // Update chache if position angle changed
+    #if !defined(G_SMALL_ANGLE_APPROXIMATION)
     if (m_last_posangle != posangle()) {
  
         // Signal parameter changes
@@ -788,6 +839,7 @@ void GModelSpatialEllipticalGauss::update() const
         m_term3 = 0.5 * (m_sin2pos / m_major2 - m_sin2pos / m_minor2);
  
     } // endif: something has changed
+    #endif
  
     // Return
     return;
@@ -420,12 +420,12 @@ double gammalib::erfinv(const double&amp; arg)
     double result;
  
     // Return NAN if out of range
-    if (std::abs(arg > 1.0)) {
+    if (std::abs(arg) > 1.0) {
         result = std::atof("NAN");
     }
  
     // Return maximum double if at range limit
-    else if (std::abs(arg == 1.0)) {
+    else if (std::abs(arg) == 1.0) {
         result = copysign(1.0, arg) * DBL_MAX;
     }
  
@@ -1252,7 +1252,24 @@ GSkyPixel GSkyMap::dir2pix(const GSkyDir&amp; dir) const
  * @exception GException::invalid_value
  *            No valid sky projection found.
  *
- * Returns the flux in the pixel with the specified @p index.
+ * Returns the flux in the pixel with the specified @p index. The flux is
+ * computed by integrating the intensity over the solid angle subtended by
+ * the pixel.
+ * Integration is done by dividing the pixel into wedges, computing the
+ * intensitites at the cornes of each wedge, averaging these intensities,
+ * and multiplying it with the solid angle of each wedge. This provides an
+ * approximation of the true pixel flux which is accurate to better than
+ * about 5%.
+ * 
+ * For a HealPix pixelisation, the pixel is divided into 12 wedges. For a
+ * WCS pixelisation, the pixel is divided into 8 wedges.
+ *
+ * @warning
+ * This method only returns correct values once the skymap is completely
+ * setup with values. Do not use this method during a setup operation as
+ * the method uses neighboring pixels for interpolation. If the map is not
+ * completely setup the neighboring pixels may still be empty, hence the
+ * flux interpolation will be wrong.
  ***************************************************************************/
 double GSkyMap::flux(const int& index, const int& map) const
 {
@@ -1290,6 +1307,13 @@ double GSkyMap::flux(const int&amp; index, const int&amp; map) const
  * 
  * For a HealPix pixelisation, the pixel is divided into 12 wedges. For a
  * WCS pixelisation, the pixel is divided into 8 wedges.
+ *
+ * @warning
+ * This method only returns correct values once the skymap is completely
+ * setup with values. Do not use this method during a setup operation as
+ * the method uses neighboring pixels for interpolation. If the map is not
+ * completely setup the neighboring pixels may still be empty, hence the
+ * flux interpolation will be wrong.
  ***************************************************************************/
 double GSkyMap::flux(const GSkyPixel& pixel, const int& map) const
 {
@@ -361,7 +361,7 @@ void GVOHub::start_hub(void)
     while (1) {
  
         // Accept connection from the client 
-    	socklen_t socket = accept(m_socket, (struct sockaddr *)&cli_addr, &clilen);
+    	int socket = accept(m_socket, (struct sockaddr *)&cli_addr, &clilen);
     	if (socket < 0) {
             std::string msg = "Client connection to socket not accepted.";
             throw GException::runtime_error(G_START_HUB, msg);