diff --git a/include/GModelSpatialRadialProfile.hpp b/include/GModelSpatialRadialProfile.hpp index 957478b..6fcc978 100644 --- a/include/GModelSpatialRadialProfile.hpp +++ b/include/GModelSpatialRadialProfile.hpp @@ -31,6 +31,7 @@ #include #include "GModelSpatialRadial.hpp" #include "GNodeArray.hpp" +#include "GSkyRegionCircle.hpp" /* __ Forward declaration ________________________________________________ */ class GXmlElement; @@ -67,17 +68,16 @@ public: virtual std::string print(const GChatter& chatter = NORMAL) const = 0; // Implemented pure virtual base class methods - virtual double eval(const double& theta, - const GEnergy& energy, - const GTime& time) const; - virtual double eval_gradients(const double& theta, - const GEnergy& energy, - const GTime& time) const; - virtual GSkyDir mc(const GEnergy& energy, - const GTime& time, - GRan& ran) const; - virtual bool contains(const GSkyDir& dir, - const double& margin = 0.0) const; + virtual double eval(const double& theta, + const GEnergy& energy, + const GTime& time, + const bool& gradients = false) const; + virtual GSkyDir mc(const GEnergy& energy, + const GTime& time, + GRan& ran) const; + virtual bool contains(const GSkyDir& dir, + const double& margin = 0.0) const; + virtual GSkyRegion* region(void) const; // Implement other methods int num_nodes(void) const; @@ -90,9 +90,11 @@ protected: void free_members(void); int cache_index(void) const; virtual double profile_value(const double& theta) const = 0; + void set_region(void) const; // Protected members - int m_num_nodes; //!< Number of profile nodes + int m_num_nodes; //!< Number of profile nodes + mutable GSkyRegionCircle m_region; //!< Bounding circle // Pre-computed radial profile struct profile { @@ -136,4 +138,19 @@ void GModelSpatialRadialProfile::num_nodes(const int& number) return; } + +/***********************************************************************//** + * @brief Return boundary sky region + * + * @return Boundary sky region. + * + * Returns a sky region that fully encloses the spatial model component. + ***************************************************************************/ +inline +GSkyRegion* GModelSpatialRadialProfile::region(void) const +{ + set_region(); + return (&m_region); +} + #endif /* GMODELSPATIALRADIALPROFILE_HPP */ diff --git a/include/GModelSpatialRadialProfileDMBurkert.hpp b/include/GModelSpatialRadialProfileDMBurkert.hpp index 815a7ed..4242e82 100644 --- a/include/GModelSpatialRadialProfileDMBurkert.hpp +++ b/include/GModelSpatialRadialProfileDMBurkert.hpp @@ -1,7 +1,7 @@ /*************************************************************************** - * GModelSpatialRadialProfileDMBurkert.hpp - DM Einasto profile class * + * GModelSpatialRadialProfileDMBurkert.hpp - DM Burkert profile class * * ----------------------------------------------------------------------- * - * copyright (C) 2016 by Juergen Knoedlseder * + * copyright (C) 2016 by Nathan Kelley-Hoskins * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -20,7 +20,7 @@ ***************************************************************************/ /** * @file GModelSpatialRadialProfileDMBurkert.hpp - * @brief Dark Matter Einasto profile model class interface definition + * @brief Dark Matter Burkert profile model class interface definition * @author Nathan Kelley-Hoskins */ @@ -41,10 +41,10 @@ class GXmlElement; /************************************************************************** * @class GModelSpatialRadialProfileDMBurkert * - * @brief Radial Dark Matter Einasto profile source model class + * @brief Radial Dark Matter Burkert profile source model class * * This class implements the spatial component of the factorised source - * model for a radial Dark Matter profile, using an Einasto density halo. + * model for a radial Dark Matter profile, using an Burkert density halo. ***************************************************************************/ class GModelSpatialRadialProfileDMBurkert : public GModelSpatialRadialProfile { @@ -59,19 +59,19 @@ public: virtual GModelSpatialRadialProfileDMBurkert& operator=(const GModelSpatialRadialProfileDMBurkert& model); // Implemented pure virtual base class methods - virtual void clear(void); + virtual void clear(void); virtual GModelSpatialRadialProfileDMBurkert* clone(void) const; - virtual std::string classname(void) const; - virtual std::string type(void) const; - virtual double theta_min(void) const; - virtual double theta_max(void) const; - virtual void read(const GXmlElement& xml); - virtual void write(GXmlElement& xml) const; - virtual std::string print(const GChatter& chatter = NORMAL) const; + virtual std::string classname(void) const; + virtual std::string type(void) const; + virtual double theta_min(void) const; + virtual double theta_max(void) const; + virtual void read(const GXmlElement& xml); + virtual void write(GXmlElement& xml) const; + virtual std::string print(const GChatter& chatter = NORMAL) const; // Other methods - double scale_radius(void) const ; - void scale_radius(const double& scale_radius ) ; + double scale_radius(void) const; + void scale_radius(const double& scale_radius); protected: // Protected methods @@ -83,29 +83,29 @@ protected: // Integration kernel for line-of-sight integral class halo_kernel_los : public GFunction { - public : + public : halo_kernel_los(const double& scale_radius , const double& halo_distance, - const double& theta ) : + const double& theta) : m_scale_radius(scale_radius), m_halo_distance(halo_distance), m_theta(theta) {} - double eval( const double& los ) ; - protected : - double m_scale_radius ; - double m_halo_distance ; - double m_theta ; - } ; + double eval(const double& los); + protected : + double m_scale_radius; + double m_halo_distance; + double m_theta; + }; // Protected members - GModelPar m_theta_min ; //!< maximum theta angle - GModelPar m_theta_max ; //!< maximum theta angle - GModelPar m_scale_radius ; //!< scale radius of halo profile - GModelPar m_halo_distance ; //!< distance from earth to halo center + GModelPar m_theta_min; //!< Minimum theta angle + GModelPar m_theta_max; //!< Maximum theta angle + GModelPar m_scale_radius; //!< Scale radius of halo profile + GModelPar m_halo_distance; //!< Distance from Earth to halo center // Cached members used for precomputation - mutable double m_last_scale_radius ; - mutable double m_mass_radius ; + mutable double m_last_scale_radius; + mutable double m_mass_radius; }; diff --git a/include/GModelSpatialRadialProfileDMEinasto.hpp b/include/GModelSpatialRadialProfileDMEinasto.hpp index 30d04af..0aa2a58 100644 --- a/include/GModelSpatialRadialProfileDMEinasto.hpp +++ b/include/GModelSpatialRadialProfileDMEinasto.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GModelSpatialRadialProfileDMEinasto.hpp - DM Einasto profile class * * ----------------------------------------------------------------------- * - * copyright (C) 2016 by Juergen Knoedlseder * + * copyright (C) 2016 by Nathan Kelley-Hoskins * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -59,21 +59,21 @@ public: virtual GModelSpatialRadialProfileDMEinasto& operator=(const GModelSpatialRadialProfileDMEinasto& model); // Implemented pure virtual base class methods - virtual void clear(void); + virtual void clear(void); virtual GModelSpatialRadialProfileDMEinasto* clone(void) const; - virtual std::string classname(void) const; - virtual std::string type(void) const; - virtual double theta_min(void) const; - virtual double theta_max(void) const; - virtual void read(const GXmlElement& xml); - virtual void write(GXmlElement& xml) const; - virtual std::string print(const GChatter& chatter = NORMAL) const; + virtual std::string classname(void) const; + virtual std::string type(void) const; + virtual double theta_min(void) const; + virtual double theta_max(void) const; + virtual void read(const GXmlElement& xml); + virtual void write(GXmlElement& xml) const; + virtual std::string print(const GChatter& chatter = NORMAL) const; // Other methods - double scale_radius(void) const ; - void scale_radius(const double& scale_radius ) ; - //double prof_val(const double& theta) ; - double j_factor( const double& minangle, const double& maxangle, const int& npoints ) const ; + double scale_radius(void) const; + void scale_radius(const double& scale_radius); + //double prof_val(const double& theta); + double j_factor( const double& minangle, const double& maxangle, const int& npoints) const; protected: // Protected methods @@ -85,33 +85,33 @@ protected: // Integration kernel for line-of-sight integral class halo_kernel_los : public GFunction { - public : - halo_kernel_los(const double& scale_radius , + public : + halo_kernel_los(const double& scale_radius, const double& halo_distance, - const double& alpha , - const double& theta ) : + const double& alpha, + const double& theta) : m_scale_radius(scale_radius), m_halo_distance(halo_distance), m_alpha(alpha), m_theta(theta) {} - double eval( const double& los ) ; - protected : - double m_scale_radius ; - double m_halo_distance ; - double m_alpha ; - double m_theta ; - } ; + double eval(const double& los); + protected : + double m_scale_radius; + double m_halo_distance; + double m_alpha; + double m_theta; + }; // Protected members - GModelPar m_scale_radius ; //!< scale radius of halo profile - GModelPar m_halo_distance ; //!< distance from earth to halo center - GModelPar m_alpha ; //!< einasto spatial power index - GModelPar m_theta_min ; //!< minimum theta angle - GModelPar m_theta_max ; //!< maximum theta angle + GModelPar m_theta_min; //!< Minimum theta angle + GModelPar m_theta_max; //!< Maximum theta angle + GModelPar m_scale_radius; //!< Scale radius of halo profile + GModelPar m_halo_distance; //!< Distance from earth to halo center + GModelPar m_alpha; //!< Einasto spatial power index // Cached members used for pre-computation - mutable double m_last_scale_radius ; - mutable double m_mass_radius ; + mutable double m_last_scale_radius; + mutable double m_mass_radius; }; diff --git a/include/GModelSpatialRadialProfileDMZhao.hpp b/include/GModelSpatialRadialProfileDMZhao.hpp index b156f69..2951822 100644 --- a/include/GModelSpatialRadialProfileDMZhao.hpp +++ b/include/GModelSpatialRadialProfileDMZhao.hpp @@ -1,7 +1,7 @@ /*************************************************************************** - * GModelSpatialRadialProfileDMZhao.hpp - DM Einasto profile class * + * GModelSpatialRadialProfileDMZhao.hpp - DM Zhao profile class * * ----------------------------------------------------------------------- * - * copyright (C) 2016 by Juergen Knoedlseder * + * copyright (C) 2016 by Nathan Kelley-Hoskins * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -20,7 +20,7 @@ ***************************************************************************/ /** * @file GModelSpatialRadialProfileDMZhao.hpp - * @brief Dark Matter Einasto profile model class interface definition + * @brief Dark Matter Zhao profile model class interface definition * @author Nathan Kelley-Hoskins */ @@ -41,10 +41,10 @@ class GXmlElement; /************************************************************************** * @class GModelSpatialRadialProfileDMZhao * - * @brief Radial Dark Matter Einasto profile source model class + * @brief Radial Dark Matter Zhao profile source model class * * This class implements the spatial component of the factorised source - * model for a radial Dark Matter profile, using an Einasto density halo. + * model for a radial Dark Matter profile, using an Zhao density halo. ***************************************************************************/ class GModelSpatialRadialProfileDMZhao : public GModelSpatialRadialProfile { @@ -59,20 +59,20 @@ public: virtual GModelSpatialRadialProfileDMZhao& operator=(const GModelSpatialRadialProfileDMZhao& model); // Implemented pure virtual base class methods - virtual void clear(void); + virtual void clear(void); virtual GModelSpatialRadialProfileDMZhao* clone(void) const; - virtual std::string classname(void) const; - virtual std::string type(void) const; - virtual double theta_min(void) const; - virtual double theta_max(void) const; - virtual void read(const GXmlElement& xml); - virtual void write(GXmlElement& xml) const; - virtual std::string print(const GChatter& chatter = NORMAL) const; + virtual std::string classname(void) const; + virtual std::string type(void) const; + virtual double theta_min(void) const; + virtual double theta_max(void) const; + virtual void read(const GXmlElement& xml); + virtual void write(GXmlElement& xml) const; + virtual std::string print(const GChatter& chatter = NORMAL) const; // Other methods - double scale_radius(void) const ; - void scale_radius(const double& scale_radius ) ; - double prof_val(const double& theta ) ; + double scale_radius(void) const; + void scale_radius(const double& scale_radius); + double prof_val(const double& theta); protected: // Protected methods @@ -84,38 +84,39 @@ protected: // Integration kernel for line-of-sight integral class halo_kernel_los : public GFunction { - public : - halo_kernel_los(const double& scale_radius , + public : + halo_kernel_los(const double& scale_radius, const double& halo_distance, - const double& alpha , - const double& beta , - const double& gamma , - const double& theta ) : + const double& alpha, + const double& beta, + const double& gamma, + const double& theta) : m_scale_radius(scale_radius), m_halo_distance(halo_distance), m_alpha(alpha), m_beta(beta), m_gamma(gamma), m_theta(theta) {} - double eval( const double& los ) ; - protected : - double m_scale_radius ; - double m_halo_distance ; - double m_alpha ; - double m_beta ; - double m_gamma ; - double m_theta ; + double eval(const double& los); + protected : + double m_scale_radius; + double m_halo_distance; + double m_alpha; + double m_beta; + double m_gamma; + double m_theta; } ; // Protected members - GModelPar m_scale_radius ; //!< scale radius of halo profile - GModelPar m_halo_distance ; //!< distance from earth to halo center - GModelPar m_alpha ; //!< power index, inverse transition region width - GModelPar m_beta ; //!< power index, slope at >> m_scale_radius - GModelPar m_gamma ; //!< power index, slope at << m_scale_radius - GModelPar m_theta_min ; //!< minimum theta angle - GModelPar m_theta_max ; //!< maximum theta angle - + GModelPar m_theta_min; //!< Minimum theta angle + GModelPar m_theta_max; //!< Maximum theta angle + GModelPar m_scale_radius; //!< Scale radius of halo profile + GModelPar m_halo_distance; //!< Distance from earth to halo center + GModelPar m_alpha; //!< Power index, inverse transition region width + GModelPar m_beta; //!< Power index, slope at >> m_scale_radius + GModelPar m_gamma; //!< Power index, slope at << m_scale_radius + + // Protected cached members mutable double m_last_scale_radius ; mutable double m_mass_radius ; }; diff --git a/pyext/GModelSpatialRadialProfile.i b/pyext/GModelSpatialRadialProfile.i index 8503dec..413d3c3 100644 --- a/pyext/GModelSpatialRadialProfile.i +++ b/pyext/GModelSpatialRadialProfile.i @@ -55,17 +55,16 @@ public: virtual double theta_max(void) const = 0; // Implemented pure virtual base class methods - virtual double eval(const double& theta, - const GEnergy& energy, - const GTime& time) const; - virtual double eval_gradients(const double& theta, - const GEnergy& energy, - const GTime& time) const; - virtual GSkyDir mc(const GEnergy& energy, - const GTime& time, - GRan& ran) const; - virtual bool contains(const GSkyDir& dir, - const double& margin = 0.0) const; + virtual double eval(const double& theta, + const GEnergy& energy, + const GTime& time, + const bool& gradients = false) const; + virtual GSkyDir mc(const GEnergy& energy, + const GTime& time, + GRan& ran) const; + virtual bool contains(const GSkyDir& dir, + const double& margin = 0.0) const; + virtual GSkyRegion* region(void) const; // Implement other methods int num_nodes(void) const; diff --git a/pyext/GModelSpatialRadialProfileDMBurkert.i b/pyext/GModelSpatialRadialProfileDMBurkert.i index 1f4b7d2..7ced7dd 100644 --- a/pyext/GModelSpatialRadialProfileDMBurkert.i +++ b/pyext/GModelSpatialRadialProfileDMBurkert.i @@ -1,7 +1,7 @@ /*************************************************************************** - * GModelSpatialRadialProfileDMBurkert.i - Einasto radial profile class * + * GModelSpatialRadialProfileDMBurkert.i - Burkert radial profile class * * ----------------------------------------------------------------------- * - * copyright (C) 2016 by Juergen Knoedlseder * + * copyright (C) 2016 by Nathan Kelley-Hoskins * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -20,7 +20,7 @@ ***************************************************************************/ /** * @file GModelSpatialRadialProfileDMBurkert.hpp - * @brief Radial Dark Matter halo for Einasto Density Profile + * @brief Radial Dark Matter halo for Burkert Density Profile * @author Nathan Kelley-Hoskins */ %{ @@ -32,10 +32,10 @@ /************************************************************************** * @class GModelSpatialRadialProfileDMBurkert * - * @brief Radial DM Einasto profile source model class + * @brief Radial DM Burkert profile source model class * * This class implements the spatial component of the factorised source - * model for a Dark Matter Einasto halo radial profile. + * model for a Dark Matter Burkert halo radial profile. ***************************************************************************/ class GModelSpatialRadialProfileDMBurkert : public GModelSpatialRadialProfile { @@ -47,25 +47,23 @@ public: virtual ~GModelSpatialRadialProfileDMBurkert(void); // Implemented pure virtual base class methods - virtual void clear(void); + virtual void clear(void); virtual GModelSpatialRadialProfileDMBurkert* clone(void) const; - virtual std::string classname(void) const; - virtual std::string type(void) const; - virtual double theta_min(void) const; - virtual double theta_max(void) const; - virtual void read(const GXmlElement& xml); - virtual void write(GXmlElement& xml) const; + virtual std::string classname(void) const; + virtual std::string type(void) const; + virtual double theta_min(void) const; + virtual double theta_max(void) const; + virtual void read(const GXmlElement& xml); + virtual void write(GXmlElement& xml) const; + + // Other methods + double scale_radius(void) const; + void scale_radius(const double& scale_radius); }; /***********************************************************************//** * @brief GModelSpatialRadialGauss class extension - * - * The eval(GSkyDir&) and eval_gradients(GSkyDir&) need to be defined in the - * extension to force swig to build also the interface for these methods that - * are implemented in the base class only. It's not clear to me why these - * methods are not inherited automatically. Maybe this could also be handled - * by a %typemap(typecheck) construct. ***************************************************************************/ %extend GModelSpatialRadialProfileDMBurkert { GModelSpatialRadialProfileDMBurkert copy() { diff --git a/pyext/GModelSpatialRadialProfileDMEinasto.i b/pyext/GModelSpatialRadialProfileDMEinasto.i index fdc8758..d7452be 100644 --- a/pyext/GModelSpatialRadialProfileDMEinasto.i +++ b/pyext/GModelSpatialRadialProfileDMEinasto.i @@ -1,7 +1,7 @@ /*************************************************************************** * GModelSpatialRadialProfileDMEinasto.i - Einasto radial profile class * * ----------------------------------------------------------------------- * - * copyright (C) 2016 by Juergen Knoedlseder * + * copyright (C) 2016 by Nathan Kelley-Hoskins * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -47,31 +47,25 @@ public: virtual ~GModelSpatialRadialProfileDMEinasto(void); // Implemented pure virtual base class methods - virtual void clear(void); + virtual void clear(void); virtual GModelSpatialRadialProfileDMEinasto* clone(void) const; - virtual std::string classname(void) const; - virtual std::string type(void) const; - virtual double theta_min(void) const; - virtual double theta_max(void) const; - virtual void read(const GXmlElement& xml); - virtual void write(GXmlElement& xml) const; + virtual std::string classname(void) const; + virtual std::string type(void) const; + virtual double theta_min(void) const; + virtual double theta_max(void) const; + virtual void read(const GXmlElement& xml); + virtual void write(GXmlElement& xml) const; - // other methods - double scale_radius(void) const ; - void scale_radius(const double& scale_radius ) ; - //double prof_val(const double& theta) ; - + // Other methods + double scale_radius(void) const; + void scale_radius(const double& scale_radius); + //double prof_val(const double& theta); + double j_factor( const double& minangle, const double& maxangle, const int& npoints) const; }; /***********************************************************************//** * @brief GModelSpatialRadialGauss class extension - * - * The eval(GSkyDir&) and eval_gradients(GSkyDir&) need to be defined in the - * extension to force swig to build also the interface for these methods that - * are implemented in the base class only. It's not clear to me why these - * methods are not inherited automatically. Maybe this could also be handled - * by a %typemap(typecheck) construct. ***************************************************************************/ %extend GModelSpatialRadialProfileDMEinasto { GModelSpatialRadialProfileDMEinasto copy() { diff --git a/pyext/GModelSpatialRadialProfileDMZhao.i b/pyext/GModelSpatialRadialProfileDMZhao.i index c26e071..e643cd9 100644 --- a/pyext/GModelSpatialRadialProfileDMZhao.i +++ b/pyext/GModelSpatialRadialProfileDMZhao.i @@ -1,7 +1,7 @@ /*************************************************************************** - * GModelSpatialRadialProfileDMZhao.i - Einasto radial profile class * + * GModelSpatialRadialProfileDMZhao.i - Zhao radial profile class * * ----------------------------------------------------------------------- * - * copyright (C) 2016 by Juergen Knoedlseder * + * copyright (C) 2016 by Nathan Kelley-Hoskins * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -20,7 +20,7 @@ ***************************************************************************/ /** * @file GModelSpatialRadialProfileDMZhao.hpp - * @brief Radial Dark Matter halo for Einasto Density Profile + * @brief Radial Dark Matter halo for Zhao Density Profile * @author Nathan Kelley-Hoskins */ %{ @@ -32,7 +32,7 @@ /************************************************************************** * @class GModelSpatialRadialProfileDMZhao * - * @brief Radial DM Einasto profile source model class + * @brief Radial DM Zhao profile source model class * * This class implements the spatial component of the factorised source * model for a Dark Matter Zhao halo radial profile. @@ -47,26 +47,24 @@ public: virtual ~GModelSpatialRadialProfileDMZhao(void); // Implemented pure virtual base class methods - virtual void clear(void); + virtual void clear(void); virtual GModelSpatialRadialProfileDMZhao* clone(void) const; - virtual std::string classname(void) const; - virtual std::string type(void) const; - virtual double theta_min(void) const; - virtual double theta_max(void) const; - virtual void read(const GXmlElement& xml); - virtual void write(GXmlElement& xml) const; - double prof_val(const double& theta ) ; + virtual std::string classname(void) const; + virtual std::string type(void) const; + virtual double theta_min(void) const; + virtual double theta_max(void) const; + virtual void read(const GXmlElement& xml); + virtual void write(GXmlElement& xml) const; + + // Other methods + double scale_radius(void) const; + void scale_radius(const double& scale_radius); + double prof_val(const double& theta); }; /***********************************************************************//** * @brief GModelSpatialRadialGauss class extension - * - * The eval(GSkyDir&) and eval_gradients(GSkyDir&) need to be defined in the - * extension to force swig to build also the interface for these methods that - * are implemented in the base class only. It's not clear to me why these - * methods are not inherited automatically. Maybe this could also be handled - * by a %typemap(typecheck) construct. ***************************************************************************/ %extend GModelSpatialRadialProfileDMZhao { GModelSpatialRadialProfileDMZhao copy() { diff --git a/src/model/GModelSpatialRadialProfile.cpp b/src/model/GModelSpatialRadialProfile.cpp index c71ca9a..815e37d 100644 --- a/src/model/GModelSpatialRadialProfile.cpp +++ b/src/model/GModelSpatialRadialProfile.cpp @@ -180,6 +180,7 @@ GModelSpatialRadialProfile& GModelSpatialRadialProfile::operator=(const GModelSp * @param[in] theta Angular distance from model centre (radians). * @param[in] energy Photon energy. * @param[in] time Photon arrival time. + * @param[in] gradients Compute gradients? * @return Model value. * * Evaluate the radial profile model for a given angular distance @p theta @@ -188,7 +189,8 @@ GModelSpatialRadialProfile& GModelSpatialRadialProfile::operator=(const GModelSp ***************************************************************************/ double GModelSpatialRadialProfile::eval(const double& theta, const GEnergy& energy, - const GTime& time) const + const GTime& time, + const bool& gradients) const { // Get pre-computation cache index int icache = cache_index(); @@ -217,29 +219,6 @@ double GModelSpatialRadialProfile::eval(const double& theta, /***********************************************************************//** - * @brief Evaluate function and gradients (in units of sr^-1) - * - * @param[in] theta Angular distance from model centre (radians). - * @param[in] energy Photon energy. - * @param[in] time Photon arrival time. - * @return Model value. - * - * Evaluates the function value. No gradient computation is implemented as - * radial models will be convolved with the instrument response and thus - * require the numerical computation of the derivatives. - * - * See the eval() method for more information. - ***************************************************************************/ -double GModelSpatialRadialProfile::eval_gradients(const double& theta, - const GEnergy& energy, - const GTime& time) const -{ - // Return value - return (eval(theta, energy, time)); -} - - -/***********************************************************************//** * @brief Return MC sky direction * * @param[in] energy Photon energy. @@ -315,6 +294,7 @@ void GModelSpatialRadialProfile::init_members(void) { // Initialise members m_num_nodes = 100; + m_region.clear(); // Initialise pre-computation cache m_profile.clear(); @@ -334,6 +314,7 @@ void GModelSpatialRadialProfile::copy_members(const GModelSpatialRadialProfile& // Copy members m_num_nodes = model.m_num_nodes; m_profile = model.m_profile; + m_region = model.m_region; // Return return; @@ -466,3 +447,19 @@ int GModelSpatialRadialProfile::cache_index(void) const // Return index return (index); } + + +/***********************************************************************//** + * @brief Set boundary sky region + ***************************************************************************/ +void GModelSpatialRadialProfile::set_region(void) const +{ + // Set sky region centre to disk centre + m_region.centre(m_ra.value(), m_dec.value()); + + // Set sky region radius to maximum theta angle + m_region.radius(theta_max()*gammalib::rad2deg); + + // Return + return; +} diff --git a/src/model/GModelSpatialRadialProfileDMBurkert.cpp b/src/model/GModelSpatialRadialProfileDMBurkert.cpp index 598a114..4926f58 100644 --- a/src/model/GModelSpatialRadialProfileDMBurkert.cpp +++ b/src/model/GModelSpatialRadialProfileDMBurkert.cpp @@ -1,7 +1,7 @@ /*************************************************************************** - * GModelSpatialRadialProfileDMBurkert.cpp - DMBurkert radial profile class * + * GModelSpatialRadialProfileDMBurkert.cpp - Burkert radial profile class * * ----------------------------------------------------------------------- * - * copyright (C) 2016 by Juergen Knoedlseder * + * copyright (C) 2016 by Nathan Kelley-Hoskins * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -20,7 +20,7 @@ ***************************************************************************/ /** * @file GModelSpatialRadialProfileDMBurkert.cpp - * @brief Radial DM Einasto profile model class implementation + * @brief Radial DM Burkert profile model class implementation * @author Nathan Kelley-Hoskins */ @@ -40,7 +40,7 @@ /* __ Globals ____________________________________________________________ */ const GModelSpatialRadialProfileDMBurkert g_radial_disk_seed; -const GModelSpatialRegistry g_radial_disk_registry(&g_radial_disk_seed); +const GModelSpatialRegistry g_radial_disk_registry(&g_radial_disk_seed); /* __ Method name definitions ____________________________________________ */ #define G_READ "GModelSpatialRadialProfileDMBurkert::read(GXmlElement&)" @@ -65,7 +65,7 @@ const GModelSpatialRegistry g_radial_disk_registry(&g_radial_disk_seed * Constructs empty radial DMBurkert profile ***************************************************************************/ GModelSpatialRadialProfileDMBurkert::GModelSpatialRadialProfileDMBurkert(void) : - GModelSpatialRadialProfile() + GModelSpatialRadialProfile() { // Initialise members init_members(); @@ -85,7 +85,7 @@ GModelSpatialRadialProfileDMBurkert::GModelSpatialRadialProfileDMBurkert(void) : * expected structure of the XML element. ***************************************************************************/ GModelSpatialRadialProfileDMBurkert::GModelSpatialRadialProfileDMBurkert(const GXmlElement& xml) : - GModelSpatialRadialProfile() + GModelSpatialRadialProfile() { // Initialise members init_members(); @@ -106,7 +106,7 @@ GModelSpatialRadialProfileDMBurkert::GModelSpatialRadialProfileDMBurkert(const G * Copies radial DMBurkert profile model from another radial profile model. ***************************************************************************/ GModelSpatialRadialProfileDMBurkert::GModelSpatialRadialProfileDMBurkert(const GModelSpatialRadialProfileDMBurkert& model) : - GModelSpatialRadialProfile(model) + GModelSpatialRadialProfile(model) { // Initialise members init_members(); @@ -227,7 +227,7 @@ double GModelSpatialRadialProfileDMBurkert::theta_min(void) const update(); // Return value - return m_theta_min.value() ; + return m_theta_min.value(); } /***********************************************************************//** @@ -237,31 +237,32 @@ double GModelSpatialRadialProfileDMBurkert::theta_min(void) const ***************************************************************************/ double GModelSpatialRadialProfileDMBurkert::theta_max(void) const { - + // Update update(); // initialize value - double theta = 0.0 ; + double theta = 0.0; - // if earth is within the significant radius, then we must integrate - // the entire profile + // If Earth is within the significant radius, then we must integrate + // the entire profile ... if ( m_halo_distance.value() < m_mass_radius ) { - theta = gammalib::pi ; + theta = gammalib::pi; + } - // if the halo is far enough away (further than the significant radius) - // then we just need to deal with the angles within the sphere of the - // significant radius. - } else { - theta = std::atan( m_mass_radius / m_halo_distance.value() ) ; + // ... otherwise, if the halo is far enough away (further than the + // significant radius) then we just need to deal with the angles within + // the sphere of the significant radius + else { + theta = std::atan( m_mass_radius / m_halo_distance.value() ) ; } - // always chose the lesser of ( mass_radius theta, theta_max ) - if ( m_theta_max.value() * gammalib::deg2rad < theta ) { - theta = m_theta_max.value() * gammalib::deg2rad ; + // Always chose the lesser of ( mass_radius theta, theta_max ) + if (m_theta_max.value() * gammalib::deg2rad < theta ) { + theta = m_theta_max.value() * gammalib::deg2rad; } // Return value - return theta ; + return theta; } @@ -274,23 +275,23 @@ double GModelSpatialRadialProfileDMBurkert::theta_max(void) const * The XML element shall have either the format * * - * - * - * + * + * + * * * - * + * * * * or * * - * - * - * + * + * + * * * - * + * * ***************************************************************************/ void GModelSpatialRadialProfileDMBurkert::read(const GXmlElement& xml) @@ -323,9 +324,12 @@ void GModelSpatialRadialProfileDMBurkert::read(const GXmlElement& xml) * The XML element will have the format * * - * - * - * + * + * + * + * + * + * * ***************************************************************************/ void GModelSpatialRadialProfileDMBurkert::write(GXmlElement& xml) const @@ -419,9 +423,9 @@ void GModelSpatialRadialProfileDMBurkert::init_members(void) m_theta_min.clear(); m_theta_min.name("Theta Min"); m_theta_min.unit("degrees"); - m_theta_min.value( 1.0e-6 ); // - m_theta_min.min(1.0e-10); // arbitrarily chosen - m_theta_min.fix(); // should always be fixed! + m_theta_min.value(1.0e-6); + m_theta_min.min(1.0e-10); // arbitrarily chosen + m_theta_min.fix(); // should always be fixed! m_theta_min.scale(1.0); m_theta_min.gradient(0.0); m_theta_min.has_grad(false); // Radial components never have gradients @@ -430,18 +434,18 @@ void GModelSpatialRadialProfileDMBurkert::init_members(void) m_theta_max.clear(); m_theta_max.name("Theta Max"); m_theta_max.unit("degrees"); - m_theta_max.value( 180.0 ); // can only go from halo center to opposite halo center - m_theta_max.min(1.0e-6); // arbitrarily chosen - m_theta_max.fix(); // should always be fixed! + m_theta_max.value(180.0); // can only go from halo center to opposite halo center + m_theta_max.min(1.0e-6); // arbitrarily chosen + m_theta_max.fix(); // should always be fixed! m_theta_max.scale(1.0); m_theta_max.gradient(0.0); m_theta_max.has_grad(false); // Radial components never have gradients // Set parameter pointer(s) - m_pars.push_back(&m_scale_radius ); + m_pars.push_back(&m_scale_radius); m_pars.push_back(&m_halo_distance); - m_pars.push_back(&m_theta_max ); - m_pars.push_back(&m_theta_min ); + m_pars.push_back(&m_theta_max); + m_pars.push_back(&m_theta_min); // Initialize precomputation cache. Note that zero values flag // uninitialised, as a zero radius is not meaningful @@ -465,14 +469,14 @@ void GModelSpatialRadialProfileDMBurkert::copy_members(const GModelSpatialRadial // Copy members. We do not have to push back the members on the parameter // stack as this should have been done by init_members() that was called // before. Otherwise we would have sigma twice on the stack. - m_scale_radius = model.m_scale_radius ; - m_halo_distance = model.m_halo_distance ; - m_theta_min = model.m_theta_min ; - m_theta_max = model.m_theta_max ; + m_scale_radius = model.m_scale_radius; + m_halo_distance = model.m_halo_distance; + m_theta_min = model.m_theta_min; + m_theta_max = model.m_theta_max; // copy cache values - m_last_scale_radius = model.m_last_scale_radius ; - m_mass_radius = model.m_mass_radius ; + m_last_scale_radius = model.m_last_scale_radius; + m_mass_radius = model.m_mass_radius; // Return return; @@ -503,25 +507,25 @@ double GModelSpatialRadialProfileDMBurkert::profile_value(const double& theta) c update(); // initialize integral value - double value = 0.0 ; + double value = 0.0; // Set up integration limits - double los_min = m_halo_distance.value() - m_mass_radius ; - double los_max = m_halo_distance.value() + m_mass_radius ; + double los_min = m_halo_distance.value() - m_mass_radius; + double los_max = m_halo_distance.value() + m_mass_radius; - // handle case where observer is within halo mass radius + // Handle case where observer is within halo mass radius if ( los_min < 0.0 ) { - los_min = 0.0 ; + los_min = 0.0 ; } // Set up integral - halo_kernel_los integrand( m_scale_radius.value(), - m_halo_distance.value(), - theta ) ; - GIntegral integral(&integrand) ; + halo_kernel_los integrand(m_scale_radius.value(), + m_halo_distance.value(), + theta); + GIntegral integral(&integrand); // Compute value - value = integral.romberg( los_min, los_max ) ; + value = integral.romberg(los_min, los_max); //double tm = theta_max() ; //std::cout << "burkert theta_max=" << tm << std::setprecision(10) << " theta=" << theta << " d=" << m_halo_distance.value() << " mr=" << m_mass_radius << " rs=" << m_scale_radius.value() << " value=" << value << std::endl; @@ -530,6 +534,7 @@ double GModelSpatialRadialProfileDMBurkert::profile_value(const double& theta) c return value; } + /***********************************************************************//** * @brief Kernel for halo density profile squared * @@ -558,52 +563,55 @@ double GModelSpatialRadialProfileDMBurkert::profile_value(const double& theta) c * @return unit * ***************************************************************************/ -double GModelSpatialRadialProfileDMBurkert::halo_kernel_los::eval( const double &los ) +double GModelSpatialRadialProfileDMBurkert::halo_kernel_los::eval(const double &los) { + // PLEASE ADD COMMENTS + double r = 0.0; + r = los * los; + r += m_halo_distance * m_halo_distance; + r -= 2.0 * los * m_halo_distance * std::cos(m_theta); + r = std::sqrt(r); - double r = 0.0 ; - r = los * los ; - r += m_halo_distance * m_halo_distance ; - r -= 2.0 * los * m_halo_distance * std::cos(m_theta) ; - r = sqrt(r) ; - - double bot = ( m_scale_radius + r ) ; - bot *= (m_scale_radius*m_scale_radius) + (r*r) ; + double bot = (m_scale_radius + r); + bot *= (m_scale_radius*m_scale_radius) + (r*r); - double f = m_scale_radius * m_scale_radius * m_scale_radius ; - f /= bot ; + double f = m_scale_radius * m_scale_radius * m_scale_radius; + f /= bot; - // squared, for annihilating dm - // would just be f if it was decaying dm - f = f * f ; + // squared, for annihilating dm + // would just be f if it was decaying dm + f = f * f; - //std::cout << "kernel::eval los=" << los << " d=" << m_halo_distance << " theta=" << m_theta << " r=" << r << " rs=" << m_scale_radius << " f=" << f << " bot=" << bot << std::endl; - return f; - + //std::cout << "kernel::eval los=" << los << " d=" << m_halo_distance << " theta=" << m_theta << " r=" << r << " rs=" << m_scale_radius << " f=" << f << " bot=" << bot << std::endl; + + // Return function value + return f; } + /***********************************************************************//** * @brief Update precomputation cache * * Computes the mass_radius calculation, determining the radius around - * the halo that contains 99.99% of the mass. For an einasto halo profile, - * this is just 80.0 * scale_radius . - * + * the halo that contains 99.99% of the mass. For an Burket halo profile, + * this is just 80.0 * scale_radius. ***************************************************************************/ void GModelSpatialRadialProfileDMBurkert::update() const { - // Update if scale radius has changed - if ( m_last_scale_radius != scale_radius() ) { + // Update if scale radius has changed + if (m_last_scale_radius != scale_radius()) { - // Store last values - m_last_scale_radius = scale_radius() ; - - // perform precomputations - // set the mass radius to 80*scale_radius, meaning - // 99.99% of the mass is contained within the mass radius, - // and integration only needs to worry about whats inside this radius. - m_mass_radius = 80.0 * scale_radius() ; + // Store last values + m_last_scale_radius = scale_radius(); + + // perform precomputations + // set the mass radius to 80*scale_radius, meaning + // 99.99% of the mass is contained within the mass radius, + // and integration only needs to worry about whats inside this radius. + m_mass_radius = 80.0 * scale_radius(); - } + } + + // Return + return; } - diff --git a/src/model/GModelSpatialRadialProfileDMEinasto.cpp b/src/model/GModelSpatialRadialProfileDMEinasto.cpp index 762fe3f..8010a77 100644 --- a/src/model/GModelSpatialRadialProfileDMEinasto.cpp +++ b/src/model/GModelSpatialRadialProfileDMEinasto.cpp @@ -1,7 +1,7 @@ /*************************************************************************** - * GModelSpatialRadialProfileDMEinasto.cpp - DMEinasto radial profile class * + * GModelSpatialRadialProfileDMEinasto.cpp - Einasto radial profile class * * ----------------------------------------------------------------------- * - * copyright (C) 2016 by Juergen Knoedlseder * + * copyright (C) 2016 by Nathan Kelley-Hoskins * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -40,7 +40,7 @@ /* __ Globals ____________________________________________________________ */ const GModelSpatialRadialProfileDMEinasto g_radial_disk_seed; -const GModelSpatialRegistry g_radial_disk_registry(&g_radial_disk_seed); +const GModelSpatialRegistry g_radial_disk_registry(&g_radial_disk_seed); /* __ Method name definitions ____________________________________________ */ #define G_READ "GModelSpatialRadialProfileDMEinasto::read(GXmlElement&)" @@ -65,7 +65,7 @@ const GModelSpatialRegistry g_radial_disk_registry(&g_radial_disk_seed * Constructs empty radial DMEinasto profile ***************************************************************************/ GModelSpatialRadialProfileDMEinasto::GModelSpatialRadialProfileDMEinasto(void) : - GModelSpatialRadialProfile() + GModelSpatialRadialProfile() { // Initialise members init_members(); @@ -85,7 +85,7 @@ GModelSpatialRadialProfileDMEinasto::GModelSpatialRadialProfileDMEinasto(void) : * expected structure of the XML element. ***************************************************************************/ GModelSpatialRadialProfileDMEinasto::GModelSpatialRadialProfileDMEinasto(const GXmlElement& xml) : - GModelSpatialRadialProfile() + GModelSpatialRadialProfile() { // Initialise members init_members(); @@ -106,7 +106,7 @@ GModelSpatialRadialProfileDMEinasto::GModelSpatialRadialProfileDMEinasto(const G * Copies radial DMEinasto profile model from another radial profile model. ***************************************************************************/ GModelSpatialRadialProfileDMEinasto::GModelSpatialRadialProfileDMEinasto(const GModelSpatialRadialProfileDMEinasto& model) : - GModelSpatialRadialProfile(model) + GModelSpatialRadialProfile(model) { // Initialise members init_members(); @@ -227,7 +227,7 @@ double GModelSpatialRadialProfileDMEinasto::theta_min(void) const update(); // Return value - return m_theta_min.value() ; + return m_theta_min.value(); } /***********************************************************************//** @@ -241,27 +241,28 @@ double GModelSpatialRadialProfileDMEinasto::theta_max(void) const // update precomputation cache update(); - double theta = 0.0 ; + double theta = 0.0; - // if earth is within the significant radius, then theta_max must - // contain the entire profile (180deg) - if ( m_halo_distance.value() < m_mass_radius ) { - theta = gammalib::pi ; + // If Earth is within the significant radius, then theta_max must + // contain the entire profile (180deg) ... + if (m_halo_distance.value() < m_mass_radius) { + theta = gammalib::pi; + } - // if the halo is far enough away (further than the mass radius) - // then we just need to deal with the angles within the sphere of the - // significant radius. - } else { - theta = std::atan( m_mass_radius / m_halo_distance.value() ) ; + // ... otherwise if the halo is far enough away (further than the mass + // radius) then we just need to deal with the angles within the sphere + // of the significant radius. + else { + theta = std::atan(m_mass_radius / m_halo_distance.value()); } - // always chose the lesser of ( mass_radius theta, theta_max ) - if ( m_theta_max.value() * gammalib::deg2rad < theta ) { - theta = m_theta_max.value() * gammalib::deg2rad ; + // Always chose the lesser of ( mass_radius theta, theta_max ) + if (m_theta_max.value() * gammalib::deg2rad < theta) { + theta = m_theta_max.value() * gammalib::deg2rad; } // Return value - return theta ; + return theta; } @@ -436,7 +437,7 @@ void GModelSpatialRadialProfileDMEinasto::init_members(void) m_theta_min.clear(); m_theta_min.name("Theta Min"); m_theta_min.unit("degrees"); - m_theta_min.value( 180.0 ); // can only go from halo center to opposite halo center + m_theta_min.value(180.0); // can only go from halo center to opposite halo center m_theta_min.min(1.0e-6); // arbitrarily chosen m_theta_min.fix(); // should always be fixed! m_theta_min.scale(1.0); @@ -447,7 +448,7 @@ void GModelSpatialRadialProfileDMEinasto::init_members(void) m_theta_max.clear(); m_theta_max.name("Theta Max"); m_theta_max.unit("degrees"); - m_theta_max.value( 1.0e-6 ); // can only go from halo center to opposite halo center + m_theta_max.value(1.0e-6); // can only go from halo center to opposite halo center m_theta_max.min(1.0e-10); // arbitrarily chosen m_theta_max.fix(); // should always be fixed! m_theta_max.scale(1.0); @@ -455,16 +456,16 @@ void GModelSpatialRadialProfileDMEinasto::init_members(void) m_theta_max.has_grad(false); // Radial components never have gradients // Set parameter pointer(s) - m_pars.push_back(&m_scale_radius ); + m_pars.push_back(&m_scale_radius); m_pars.push_back(&m_halo_distance); - m_pars.push_back(&m_alpha ); - m_pars.push_back(&m_theta_min ); - m_pars.push_back(&m_theta_max ); + m_pars.push_back(&m_alpha); + m_pars.push_back(&m_theta_min); + m_pars.push_back(&m_theta_max); // Initialize precomputation cache. Note that zero values flag // uninitialised, as a zero radius is not meaningful - m_last_scale_radius = 0.0 ; - m_mass_radius = 0.0 ; + m_last_scale_radius = 0.0; + m_mass_radius = 0.0; // Return return; @@ -483,15 +484,15 @@ void GModelSpatialRadialProfileDMEinasto::copy_members(const GModelSpatialRadial // Copy members. We do not have to push back the members on the parameter // stack as this should have been done by init_members() that was called // before. Otherwise we would have sigma twice on the stack. - m_scale_radius = model.m_scale_radius ; - m_halo_distance = model.m_halo_distance ; - m_alpha = model.m_alpha ; - m_theta_min = model.m_theta_min ; - m_theta_max = model.m_theta_max ; + m_scale_radius = model.m_scale_radius; + m_halo_distance = model.m_halo_distance; + m_alpha = model.m_alpha; + m_theta_min = model.m_theta_min; + m_theta_max = model.m_theta_max; // copy cache values - m_last_scale_radius = model.m_last_scale_radius ; - m_mass_radius = model.m_mass_radius ; + m_last_scale_radius = model.m_last_scale_radius; + m_mass_radius = model.m_mass_radius; // Return return; @@ -521,37 +522,37 @@ double GModelSpatialRadialProfileDMEinasto::profile_value(const double& theta) c update(); // initialize integral value - double value = 0.0 ; + double value = 0.0; // Set up integration limits - double los_min = m_halo_distance.value() - m_mass_radius ; - double los_max = m_halo_distance.value() + m_mass_radius ; + double los_min = m_halo_distance.value() - m_mass_radius; + double los_max = m_halo_distance.value() + m_mass_radius; // handle case where observer is within halo mass radius - if ( los_min < 0.0 ) { - los_min = 0.0 ; + if (los_min < 0.0) { + los_min = 0.0; } // Set up integral - halo_kernel_los integrand( m_scale_radius.value(), - m_halo_distance.value(), - m_alpha.value(), - theta ) ; - GIntegral integral(&integrand) ; - integral.max_iter( 30 ) ; + halo_kernel_los integrand(m_scale_radius.value(), + m_halo_distance.value(), + m_alpha.value(), + theta); + GIntegral integral(&integrand); + integral.max_iter(30); // also play with eps() or max_iter(), bounds() // Set up integration boundaries // As there is usually an infinity at the halo center, this splits // the integral at the m_halo_distance. - std::vector bounds ; - bounds.push_back( los_min ) ; - bounds.push_back( los_max ) ; - bounds.push_back( m_halo_distance.value() ); + std::vector bounds; + bounds.push_back(los_min); + bounds.push_back(los_max); + bounds.push_back(m_halo_distance.value()); // Compute value - value = integral.romberg( bounds ) ; + value = integral.romberg(bounds); std::cout << "profile_value=" << value << std::endl; @@ -559,6 +560,7 @@ double GModelSpatialRadialProfileDMEinasto::profile_value(const double& theta) c return value; } + /***********************************************************************//** * @brief Kernel for halo density profile squared * @@ -593,53 +595,57 @@ double GModelSpatialRadialProfileDMEinasto::profile_value(const double& theta) c ***************************************************************************/ double GModelSpatialRadialProfileDMEinasto::halo_kernel_los::eval( const double &los ) { + // PLEASE ADD COMMENTS + double g = 0.0; + g = los * los; + g += m_halo_distance * m_halo_distance; + g -= 2.0 * los * m_halo_distance * std::cos(m_theta); + g = std::sqrt(g); + g /= m_scale_radius; - double g = 0.0 ; - g = los * los ; - g += m_halo_distance * m_halo_distance ; - g -= 2.0 * los * m_halo_distance * std::cos(m_theta) ; - g = sqrt(g) ; - g /= m_scale_radius ; + double f = 0.0; + f = std::pow(g, m_alpha); + f -= 1.0; + f *= -2.0 / m_alpha; + f = std::exp(f); - double f = 0.0 ; - f = pow( g , m_alpha ) ; - f -= 1.0 ; - f *= -2.0 / m_alpha ; - f = std::exp( f ) ; - - // squared, for annihilating dm - // would just be f if it was decaying dm - f = f * f ; + // squared, for annihilating dm + // would just be f if it was decaying dm + f = f * f ; - //std::cout << "kernel::eval los=" << los << " theta=" << m_theta << " d=" << m_halo_distance << " rs=" << m_scale_radius << " alpha=" << m_alpha << " g=" << g << " f=" << std::setprecision(12)<< f << std::endl; - - return f; + //std::cout << "kernel::eval los=" << los << " theta=" << m_theta << " d=" << m_halo_distance << " rs=" << m_scale_radius << " alpha=" << m_alpha << " g=" << g << " f=" << std::setprecision(12)<< f << std::endl; + // Return function value + return f; } + /***********************************************************************//** * @brief Update precomputation cache * * Computes the m_mass_radius calculation, determining the radius around - * the halo that contains 99.99% of the mass. For an einasto halo profile, + * the halo that contains 99.99% of the mass. For an Einasto halo profile, * this is just 10.0 * scale_radius . - * ***************************************************************************/ void GModelSpatialRadialProfileDMEinasto::update() const { - // Update if scale radius has changed - if ( m_last_scale_radius != scale_radius() ) { + // Update if scale radius has changed + if (m_last_scale_radius != scale_radius()) { - // Store last values - m_last_scale_radius = scale_radius() ; + // Store last values + m_last_scale_radius = scale_radius(); - // perform precomputations - m_mass_radius = 10.0 * scale_radius() ; + // perform precomputations + m_mass_radius = 10.0 * scale_radius(); + + } - } + // Return + return; } + /***********************************************************************//** * @brief Calculate Halo J Factor * @@ -649,30 +655,31 @@ void GModelSpatialRadialProfileDMEinasto::update() const * Calculates the halo's j-factor of the halo by integrating across the * profile between two angles. * + * @todo This should be done using a proper numerical integration ***************************************************************************/ -double GModelSpatialRadialProfileDMEinasto::j_factor( const double& minangle, const double& maxangle, const int& npoints) const +double GModelSpatialRadialProfileDMEinasto::j_factor(const double& minangle, + const double& maxangle, + const int& npoints) const { + // init variables + double minradian = minangle * gammalib::deg2rad; + double maxradian = maxangle * gammalib::deg2rad; + double dr = (maxradian - minradian) / npoints; + double r = 0.0; + double total = 0.0; - // init variables - double minradian = minangle * gammalib::deg2rad ; - double maxradian = maxangle * gammalib::deg2rad ; - double dr = ( maxradian - minradian ) / npoints ; - double r = 0.0 ; - double total = 0.0 ; - - // loop over different radii in the profile - for (int i = 0; i #endif +//#include #include "GException.hpp" #include "GTools.hpp" #include "GMath.hpp" #include "GXmlElement.hpp" #include "GModelSpatialRadialProfileDMZhao.hpp" #include "GModelSpatialRegistry.hpp" -#include /* __ Constants __________________________________________________________ */ /* __ Globals ____________________________________________________________ */ const GModelSpatialRadialProfileDMZhao g_radial_disk_seed; -const GModelSpatialRegistry g_radial_disk_registry(&g_radial_disk_seed); +const GModelSpatialRegistry g_radial_disk_registry(&g_radial_disk_seed); /* __ Method name definitions ____________________________________________ */ #define G_READ "GModelSpatialRadialProfileDMZhao::read(GXmlElement&)" @@ -65,7 +65,7 @@ const GModelSpatialRegistry g_radial_disk_registry(&g_radial_disk_seed * Constructs empty radial DMZhao profile ***************************************************************************/ GModelSpatialRadialProfileDMZhao::GModelSpatialRadialProfileDMZhao(void) : - GModelSpatialRadialProfile() + GModelSpatialRadialProfile() { // Initialise members init_members(); @@ -85,7 +85,7 @@ GModelSpatialRadialProfileDMZhao::GModelSpatialRadialProfileDMZhao(void) : * expected structure of the XML element. ***************************************************************************/ GModelSpatialRadialProfileDMZhao::GModelSpatialRadialProfileDMZhao(const GXmlElement& xml) : - GModelSpatialRadialProfile() + GModelSpatialRadialProfile() { // Initialise members init_members(); @@ -106,7 +106,7 @@ GModelSpatialRadialProfileDMZhao::GModelSpatialRadialProfileDMZhao(const GXmlEle * Copies radial DMZhao profile model from another radial profile model. ***************************************************************************/ GModelSpatialRadialProfileDMZhao::GModelSpatialRadialProfileDMZhao(const GModelSpatialRadialProfileDMZhao& model) : - GModelSpatialRadialProfile(model) + GModelSpatialRadialProfile(model) { // Initialise members init_members(); @@ -227,9 +227,10 @@ double GModelSpatialRadialProfileDMZhao::theta_min(void) const update(); // Return value - return m_theta_min.value() ; + return m_theta_min.value(); } + /***********************************************************************//** * @brief Return maximum model radius (in radians) * @@ -241,27 +242,28 @@ double GModelSpatialRadialProfileDMZhao::theta_max(void) const // update precomputation cache update(); - double theta = 0.0 ; + double theta = 0.0; - // if earth is within the significant radius, then theta_max must - // contain the entire profile (180deg) - if ( m_halo_distance.value() < m_mass_radius ) { - theta = gammalib::pi ; + // If Earth is within the significant radius, then theta_max must + // contain the entire profile (180deg) ... + if (m_halo_distance.value() < m_mass_radius) { + theta = gammalib::pi; + } - // if the halo is far enough away (further than the significant radius) - // then we just need to deal with the angles within the sphere of the - // significant radius. - } else { - theta = std::atan( m_mass_radius / m_halo_distance.value() ) ; + // ... otherwise, if the halo is far enough away (further than the + // significant radius) then we just need to deal with the angles within + // the sphere of the significant radius. + else { + theta = std::atan(m_mass_radius / m_halo_distance.value()); } - // always chose the lesser of ( mass_radius theta, theta_max ) - if ( m_theta_max.value() * gammalib::deg2rad < theta ) { - theta = m_theta_max.value() * gammalib::deg2rad ; + // Always chose the lesser of ( mass_radius theta, theta_max ) + if (m_theta_max.value() * gammalib::deg2rad < theta) { + theta = m_theta_max.value() * gammalib::deg2rad; } // Return value - return theta ; + return theta; } @@ -507,13 +509,13 @@ void GModelSpatialRadialProfileDMZhao::init_members(void) m_theta_max.has_grad(false); // Radial components never have gradients // Set parameter pointer(s) - m_pars.push_back(&m_scale_radius ); + m_pars.push_back(&m_scale_radius); m_pars.push_back(&m_halo_distance); - m_pars.push_back(&m_alpha ); - m_pars.push_back(&m_beta ); - m_pars.push_back(&m_gamma ); - m_pars.push_back(&m_theta_min ); - m_pars.push_back(&m_theta_max ); + m_pars.push_back(&m_alpha); + m_pars.push_back(&m_beta); + m_pars.push_back(&m_gamma); + m_pars.push_back(&m_theta_min); + m_pars.push_back(&m_theta_max); // Return return; @@ -532,13 +534,13 @@ void GModelSpatialRadialProfileDMZhao::copy_members(const GModelSpatialRadialPro // Copy members. We do not have to push back the members on the parameter // stack as this should have been done by init_members() that was called // before. Otherwise we would have sigma twice on the stack. - m_scale_radius = model.m_scale_radius ; - m_halo_distance = model.m_halo_distance ; - m_alpha = model.m_alpha ; - m_beta = model.m_beta ; - m_gamma = model.m_gamma ; - m_theta_max = model.m_theta_min ; - m_theta_max = model.m_theta_max ; + m_scale_radius = model.m_scale_radius; + m_halo_distance = model.m_halo_distance; + m_alpha = model.m_alpha; + m_beta = model.m_beta; + m_gamma = model.m_gamma; + m_theta_max = model.m_theta_min; + m_theta_max = model.m_theta_max; // Return return; @@ -568,32 +570,32 @@ double GModelSpatialRadialProfileDMZhao::profile_value(const double& theta) cons update(); // initialize integral value - double value = 0.0 ; + double value = 0.0; // Set up integration limits - double los_min = m_halo_distance.value() - m_mass_radius ; - double los_max = m_halo_distance.value() + m_mass_radius ; + double los_min = m_halo_distance.value() - m_mass_radius; + double los_max = m_halo_distance.value() + m_mass_radius; // Set up integral - halo_kernel_los integrand( m_scale_radius.value(), - m_halo_distance.value(), - m_alpha.value(), - m_beta.value(), - m_gamma.value(), - theta ) ; - GIntegral integral(&integrand) ; - integral.max_iter( 25 ) ; + halo_kernel_los integrand(m_scale_radius.value(), + m_halo_distance.value(), + m_alpha.value(), + m_beta.value(), + m_gamma.value(), + theta); + GIntegral integral(&integrand); + integral.max_iter(25); // Set up integration boundaries // As there is usually an infinity at the halo center, this splits // the integral at the m_halo_distance. - std::vector bounds ; - bounds.push_back( los_min ) ; - bounds.push_back( los_max ) ; - bounds.push_back( m_halo_distance.value() ); + std::vector bounds; + bounds.push_back(los_min); + bounds.push_back(los_max); + bounds.push_back(m_halo_distance.value()); // Compute value - value = integral.romberg( bounds ) ; + value = integral.romberg(bounds); //std::cout << " theta=" << theta << " profile_value=" << value << std::endl; @@ -601,6 +603,7 @@ double GModelSpatialRadialProfileDMZhao::profile_value(const double& theta) cons return value; } + /***********************************************************************//** * @brief Kernel for zhao halo density profile squared * @@ -636,55 +639,64 @@ double GModelSpatialRadialProfileDMZhao::profile_value(const double& theta) cons ***************************************************************************/ double GModelSpatialRadialProfileDMZhao::halo_kernel_los::eval( const double &los ) { + // ADD COMMENTS + double g = 0.0; + g = los * los; + g += m_halo_distance * m_halo_distance; + g -= 2.0 * los * m_halo_distance * std::cos(m_theta); + g = std::sqrt(g) ; + g /= m_scale_radius ; - double g = 0.0 ; - g = los * los ; - g += m_halo_distance * m_halo_distance ; - g -= 2.0 * los * m_halo_distance * std::cos(m_theta) ; - g = sqrt(g) ; - g /= m_scale_radius ; - - double f = 0.0 ; - f = pow( g , m_alpha ) ; - f += 1 ; - f = pow( f, (m_beta-m_gamma)/m_alpha ) ; - f *= pow( g, m_gamma ) ; - f = 1 / f ; - - // squared, for annihilating dm - // would just be f if it was decaying dm - f = f * f ; + double f = 0.0 ; + f = std::pow(g, m_alpha); + f += 1 ; + f = std::pow(f, (m_beta-m_gamma)/m_alpha); + f *= std::pow(g, m_gamma); + f = 1 / f ; + + // squared, for annihilating dm + // would just be f if it was decaying dm + f = f * f; - //std::cout << std::setprecision(10) << "kern_los::eval los=" << los << " hd=" << m_halo_distance << " theta=" << m_theta << " alpha=" << m_alpha << " beta=" << m_beta << " gamma=" << m_gamma << " g=" << g << " f=" << f << std::endl ; - - return f; + //std::cout << std::setprecision(10) << "kern_los::eval los=" << los << " hd=" << m_halo_distance << " theta=" << m_theta << " alpha=" << m_alpha << " beta=" << m_beta << " gamma=" << m_gamma << " g=" << g << " f=" << f << std::endl ; + // Return function value + return f; } + /***********************************************************************//** * @brief Update precomputation cache * * Computes the m_mass_radius calculation, determining the radius around - * the halo that contains 99.99% of the mass. For a zhao halo profile, + * the halo that contains 99.99% of the mass. For a Zhao halo profile, * this is just 10.0 * scale_radius . - * ***************************************************************************/ void GModelSpatialRadialProfileDMZhao::update() const { - // Update if scale radius has changed - if ( m_last_scale_radius != scale_radius() ) { + // Update if scale radius has changed + if (m_last_scale_radius != scale_radius()) { - // Store last values - m_last_scale_radius = scale_radius() ; + // Store last values + m_last_scale_radius = scale_radius(); - // perform precomputations - m_mass_radius = 10.0 * scale_radius() ; + // perform precomputations + m_mass_radius = 10.0 * scale_radius(); - } + } + + // Return + return; } -double GModelSpatialRadialProfileDMZhao::prof_val( const double& theta ) + +/***********************************************************************//** + * @brief Compute profile value + * + * @return Profile value + ***************************************************************************/ +double GModelSpatialRadialProfileDMZhao::prof_val(const double& theta) { - return this->profile_value(theta) ; + return this->profile_value(theta); }