diff --git a/include/GApplicationPars.hpp b/include/GApplicationPars.hpp index 54795a2..51673fd 100644 --- a/include/GApplicationPars.hpp +++ b/include/GApplicationPars.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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& GApplicationPars::operator[](const int& index) const inline int GApplicationPars::size(void) const { - return (m_pars.size()); + return (int)m_pars.size(); } diff --git a/include/GArf.hpp b/include/GArf.hpp index d934a12..6809e15 100644 --- a/include/GArf.hpp +++ b/include/GArf.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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& GArf::operator[](const int& index) const inline int GArf::size(void) const { - return (m_specresp.size()); + return (int)m_specresp.size(); } diff --git a/include/GEnergies.hpp b/include/GEnergies.hpp index c0f5365..3b97cad 100644 --- a/include/GEnergies.hpp +++ b/include/GEnergies.hpp @@ -149,7 +149,7 @@ const GEnergy& GEnergies::operator[](const int& index) const inline int GEnergies::size(void) const { - return (m_energies.size()); + return (int)m_energies.size(); } diff --git a/include/GFits.hpp b/include/GFits.hpp index 72beb3c..75d7152 100644 --- a/include/GFits.hpp +++ b/include/GFits.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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& extname) const inline int GFits::size(void) const { - return (m_hdu.size()); + return (int)m_hdu.size(); } diff --git a/include/GFitsHeader.hpp b/include/GFitsHeader.hpp index d469fd6..7e6b233 100644 --- a/include/GFitsHeader.hpp +++ b/include/GFitsHeader.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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& GFitsHeader::operator[](const std::string& keyname) const inline int GFitsHeader::size(void) const { - return (m_cards.size()); + return (int)m_cards.size(); } diff --git a/include/GModel.hpp b/include/GModel.hpp index a21c6fd..8eefb99 100644 --- a/include/GModel.hpp +++ b/include/GModel.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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& GModel::operator[](const int& index) const inline int GModel::size(void) const { - return (m_pars.size()); + return (int)m_pars.size(); } diff --git a/include/GModelSpatial.hpp b/include/GModelSpatial.hpp index afcc720..73bd91e 100644 --- a/include/GModelSpatial.hpp +++ b/include/GModelSpatial.hpp @@ -157,7 +157,7 @@ const GModelPar& GModelSpatial::operator[](const int& index) const inline int GModelSpatial::size(void) const { - return (m_pars.size()); + return (int)m_pars.size(); } #endif /* GMODELSPATIAL_HPP */ diff --git a/include/GModelSpectral.hpp b/include/GModelSpectral.hpp index a44b720..1f6ec3e 100644 --- a/include/GModelSpectral.hpp +++ b/include/GModelSpectral.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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& GModelSpectral::operator[](const int& index) const inline int GModelSpectral::size(void) const { - return (m_pars.size()); + return (int)m_pars.size(); } #endif /* GMODELSPECTRAL_HPP */ diff --git a/include/GModelSpectralNodes.hpp b/include/GModelSpectralNodes.hpp index 4d84f6c..8a18a7f 100644 --- a/include/GModelSpectralNodes.hpp +++ b/include/GModelSpectralNodes.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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 */ diff --git a/include/GModelTemporal.hpp b/include/GModelTemporal.hpp index 9fb11f7..efc5dba 100644 --- a/include/GModelTemporal.hpp +++ b/include/GModelTemporal.hpp @@ -134,7 +134,7 @@ const GModelPar& GModelTemporal::operator[](const int& index) const inline int GModelTemporal::size(void) const { - return (m_pars.size()); + return (int)m_pars.size(); } #endif /* GMODELTEMPORAL_HPP */ diff --git a/include/GModels.hpp b/include/GModels.hpp index dc248c7..a1015da 100644 --- a/include/GModels.hpp +++ b/include/GModels.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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& index) const inline int GModels::size(void) const { - return (m_models.size()); + return (int)m_models.size(); } diff --git a/include/GNodeArray.hpp b/include/GNodeArray.hpp index 9375eaa..74d0d87 100644 --- a/include/GNodeArray.hpp +++ b/include/GNodeArray.hpp @@ -184,7 +184,7 @@ const double& GNodeArray::operator[](const int& index) const inline int GNodeArray::size(void) const { - return (m_node.size()); + return (int)m_node.size(); } diff --git a/include/GObservations.hpp b/include/GObservations.hpp index 6cf8c8c..f7f3809 100644 --- a/include/GObservations.hpp +++ b/include/GObservations.hpp @@ -238,7 +238,7 @@ const GObservation* GObservations::operator[](const int& index) const inline int GObservations::size(void) const { - return (m_obs.size()); + return (int)m_obs.size(); } diff --git a/include/GOptimizerPars.hpp b/include/GOptimizerPars.hpp index 16f066a..a7336c3 100644 --- a/include/GOptimizerPars.hpp +++ b/include/GOptimizerPars.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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& index) const inline int GOptimizerPars::size(void) const { - return (m_pars.size()); + return (int)m_pars.size(); } diff --git a/include/GPha.hpp b/include/GPha.hpp index 7c7657e..d3d5d53 100644 --- a/include/GPha.hpp +++ b/include/GPha.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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& GPha::operator[](const int& index) const inline int GPha::size(void) const { - return (m_counts.size()); + return (int)m_counts.size(); } diff --git a/include/GPhotons.hpp b/include/GPhotons.hpp index 115e6dc..78f06c0 100644 --- a/include/GPhotons.hpp +++ b/include/GPhotons.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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(); } diff --git a/include/GSkyRegions.hpp b/include/GSkyRegions.hpp index d99ae2c..9086dc8 100644 --- a/include/GSkyRegions.hpp +++ b/include/GSkyRegions.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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& index) const inline int GSkyRegions::size(void) const { - return (m_regions.size()); + return (int)m_regions.size(); } diff --git a/include/GTestSuite.hpp b/include/GTestSuite.hpp index fb2f45a..5b5e05c 100644 --- a/include/GTestSuite.hpp +++ b/include/GTestSuite.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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(); } diff --git a/include/GTestSuites.hpp b/include/GTestSuites.hpp index cdf06bd..6d7f40b 100644 --- a/include/GTestSuites.hpp +++ b/include/GTestSuites.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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& index) const inline int GTestSuites::size(void) const { - return m_testsuites.size(); + return (int)m_testsuites.size(); } diff --git a/include/GTimes.hpp b/include/GTimes.hpp index c3ce558..7bd8360 100644 --- a/include/GTimes.hpp +++ b/include/GTimes.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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(); } diff --git a/include/GVOHub.hpp b/include/GVOHub.hpp index 4d33037..3ab9670 100644 --- a/include/GVOHub.hpp +++ b/include/GVOHub.hpp @@ -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 m_clients; //!< Clients }; diff --git a/include/GXmlNode.hpp b/include/GXmlNode.hpp index 3eca44c..a7d0d1d 100644 --- a/include/GXmlNode.hpp +++ b/include/GXmlNode.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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(); } diff --git a/inst/cta/include/GCTAEventList.hpp b/inst/cta/include/GCTAEventList.hpp index ff3cdc1..53ea830 100644 --- a/inst/cta/include/GCTAEventList.hpp +++ b/inst/cta/include/GCTAEventList.hpp @@ -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(); } diff --git a/inst/cta/include/GCTAModelRadial.hpp b/inst/cta/include/GCTAModelRadial.hpp index f5b5afd..cf29afe 100644 --- a/inst/cta/include/GCTAModelRadial.hpp +++ b/inst/cta/include/GCTAModelRadial.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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 diff --git a/inst/cta/include/GCTAModelRadialPolynom.hpp b/inst/cta/include/GCTAModelRadialPolynom.hpp index eab3775..ebfcf8b 100644 --- a/inst/cta/include/GCTAModelRadialPolynom.hpp +++ b/inst/cta/include/GCTAModelRadialPolynom.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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 */ diff --git a/inst/cta/include/GCTAOnOffObservations.hpp b/inst/cta/include/GCTAOnOffObservations.hpp index 4de9d99..40917f5 100644 --- a/inst/cta/include/GCTAOnOffObservations.hpp +++ b/inst/cta/include/GCTAOnOffObservations.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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& index) inline int GCTAOnOffObservations::size(void) const { - return (m_obs.size()); + return (int)m_obs.size(); } diff --git a/inst/cta/src/GCTAResponse_helpers.cpp b/inst/cta/src/GCTAResponse_helpers.cpp index d4cc4c8..6365f87 100644 --- a/inst/cta/src/GCTAResponse_helpers.cpp +++ b/inst/cta/src/GCTAResponse_helpers.cpp @@ -864,11 +864,15 @@ double cta_irf_elliptical_kern_omega::eval(const double& 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 diff --git a/inst/lat/include/GLATEventCube.hpp b/inst/lat/include/GLATEventCube.hpp index 166912d..ffe3389 100644 --- a/inst/lat/include/GLATEventCube.hpp +++ b/inst/lat/include/GLATEventCube.hpp @@ -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 */ diff --git a/inst/lat/include/GLATEventList.hpp b/inst/lat/include/GLATEventList.hpp index 865550a..9e2611a 100644 --- a/inst/lat/include/GLATEventList.hpp +++ b/inst/lat/include/GLATEventList.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(); } diff --git a/inst/lat/include/GLATResponse.hpp b/inst/lat/include/GLATResponse.hpp index 66b4b05..5d0ad53 100644 --- a/inst/lat/include/GLATResponse.hpp +++ b/inst/lat/include/GLATResponse.hpp @@ -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(); } diff --git a/inst/mwl/include/GMWLSpectrum.hpp b/inst/mwl/include/GMWLSpectrum.hpp index 5a818d9..f54d3a4 100644 --- a/inst/mwl/include/GMWLSpectrum.hpp +++ b/inst/mwl/include/GMWLSpectrum.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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(); } diff --git a/pyext/GSkyMap.i b/pyext/GSkyMap.i index eb3ea8f..50f3dce 100644 --- a/pyext/GSkyMap.i +++ b/pyext/GSkyMap.i @@ -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 diff --git a/src/model/GModelSpatialEllipticalGauss.cpp b/src/model/GModelSpatialEllipticalGauss.cpp index 1728024..bf63455 100644 --- a/src/model/GModelSpatialEllipticalGauss.cpp +++ b/src/model/GModelSpatialEllipticalGauss.cpp @@ -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(&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& 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& 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& 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& 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& 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& 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& 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& 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& 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& 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; diff --git a/src/numerics/GMath.cpp b/src/numerics/GMath.cpp index 40fd5e2..16e3854 100644 --- a/src/numerics/GMath.cpp +++ b/src/numerics/GMath.cpp @@ -420,12 +420,12 @@ double gammalib::erfinv(const double& 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; } diff --git a/src/sky/GSkyMap.cpp b/src/sky/GSkyMap.cpp index 95287c9..edd9ed9 100644 --- a/src/sky/GSkyMap.cpp +++ b/src/sky/GSkyMap.cpp @@ -1252,7 +1252,24 @@ GSkyPixel GSkyMap::dir2pix(const GSkyDir& 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& index, const int& 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 { diff --git a/src/vo/GVOHub.cpp b/src/vo/GVOHub.cpp index 08c35c6..e04be5f 100644 --- a/src/vo/GVOHub.cpp +++ b/src/vo/GVOHub.cpp @@ -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);