gammalib / gammalib

Commit 58af4ad623c1cb118ddb34f4bac9c95415e1998d

Authored by Jürgen Knödlseder
1 parent 9cc600e0

Adapt DM radial profile code to new spatial model interface

Inline Side-by-side
Showing 12 changed files with 501 additions and 470 deletions
include/GModelSpatialRadialProfile.hpp
  View file @58af4ad
... ... @@ -31,6 +31,7 @@
31 31 #include <string>
32 32 #include "GModelSpatialRadial.hpp"
33 33 #include "GNodeArray.hpp"
  34 +#include "GSkyRegionCircle.hpp"
34 35  
35 36 /* __ Forward declaration ________________________________________________ */
36 37 class GXmlElement;
... ... @@ -67,17 +68,16 @@ public:
67 68 virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
68 69  
69 70 // Implemented pure virtual base class methods
70   - virtual double eval(const double& theta,
71   - const GEnergy& energy,
72   - const GTime& time) const;
73   - virtual double eval_gradients(const double& theta,
74   - const GEnergy& energy,
75   - const GTime& time) const;
76   - virtual GSkyDir mc(const GEnergy& energy,
77   - const GTime& time,
78   - GRan& ran) const;
79   - virtual bool contains(const GSkyDir& dir,
80   - const double& margin = 0.0) const;
  71 + virtual double eval(const double& theta,
  72 + const GEnergy& energy,
  73 + const GTime& time,
  74 + const bool& gradients = false) const;
  75 + virtual GSkyDir mc(const GEnergy& energy,
  76 + const GTime& time,
  77 + GRan& ran) const;
  78 + virtual bool contains(const GSkyDir& dir,
  79 + const double& margin = 0.0) const;
  80 + virtual GSkyRegion* region(void) const;
81 81  
82 82 // Implement other methods
83 83 int num_nodes(void) const;
... ... @@ -90,9 +90,11 @@ protected:
90 90 void free_members(void);
91 91 int cache_index(void) const;
92 92 virtual double profile_value(const double& theta) const = 0;
  93 + void set_region(void) const;
93 94  
94 95 // Protected members
95   - int m_num_nodes; //!< Number of profile nodes
  96 + int m_num_nodes; //!< Number of profile nodes
  97 + mutable GSkyRegionCircle m_region; //!< Bounding circle
96 98  
97 99 // Pre-computed radial profile
98 100 struct profile {
... ... @@ -136,4 +138,19 @@ void GModelSpatialRadialProfile::num_nodes(const int&amp; number)
136 138 return;
137 139 }
138 140  
  141 +
  142 +/***********************************************************************//**
  143 + * @brief Return boundary sky region
  144 + *
  145 + * @return Boundary sky region.
  146 + *
  147 + * Returns a sky region that fully encloses the spatial model component.
  148 + ***************************************************************************/
  149 +inline
  150 +GSkyRegion* GModelSpatialRadialProfile::region(void) const
  151 +{
  152 + set_region();
  153 + return (&m_region);
  154 +}
  155 +
139 156 #endif /* GMODELSPATIALRADIALPROFILE_HPP */
... ...
include/GModelSpatialRadialProfileDMBurkert.hpp
  View file @58af4ad
1 1 /***************************************************************************
2   - * GModelSpatialRadialProfileDMBurkert.hpp - DM Einasto profile class *
  2 + * GModelSpatialRadialProfileDMBurkert.hpp - DM Burkert profile class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2016 by Nathan Kelley-Hoskins *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -20,7 +20,7 @@
20 20 ***************************************************************************/
21 21 /**
22 22 * @file GModelSpatialRadialProfileDMBurkert.hpp
23   - * @brief Dark Matter Einasto profile model class interface definition
  23 + * @brief Dark Matter Burkert profile model class interface definition
24 24 * @author Nathan Kelley-Hoskins
25 25 */
26 26  
... ... @@ -41,10 +41,10 @@ class GXmlElement;
41 41 /**************************************************************************
42 42 * @class GModelSpatialRadialProfileDMBurkert
43 43 *
44   - * @brief Radial Dark Matter Einasto profile source model class
  44 + * @brief Radial Dark Matter Burkert profile source model class
45 45 *
46 46 * This class implements the spatial component of the factorised source
47   - * model for a radial Dark Matter profile, using an Einasto density halo.
  47 + * model for a radial Dark Matter profile, using an Burkert density halo.
48 48 ***************************************************************************/
49 49 class GModelSpatialRadialProfileDMBurkert : public GModelSpatialRadialProfile {
50 50  
... ... @@ -59,19 +59,19 @@ public:
59 59 virtual GModelSpatialRadialProfileDMBurkert& operator=(const GModelSpatialRadialProfileDMBurkert& model);
60 60  
61 61 // Implemented pure virtual base class methods
62   - virtual void clear(void);
  62 + virtual void clear(void);
63 63 virtual GModelSpatialRadialProfileDMBurkert* clone(void) const;
64   - virtual std::string classname(void) const;
65   - virtual std::string type(void) const;
66   - virtual double theta_min(void) const;
67   - virtual double theta_max(void) const;
68   - virtual void read(const GXmlElement& xml);
69   - virtual void write(GXmlElement& xml) const;
70   - virtual std::string print(const GChatter& chatter = NORMAL) const;
  64 + virtual std::string classname(void) const;
  65 + virtual std::string type(void) const;
  66 + virtual double theta_min(void) const;
  67 + virtual double theta_max(void) const;
  68 + virtual void read(const GXmlElement& xml);
  69 + virtual void write(GXmlElement& xml) const;
  70 + virtual std::string print(const GChatter& chatter = NORMAL) const;
71 71  
72 72 // Other methods
73   - double scale_radius(void) const ;
74   - void scale_radius(const double& scale_radius ) ;
  73 + double scale_radius(void) const;
  74 + void scale_radius(const double& scale_radius);
75 75  
76 76 protected:
77 77 // Protected methods
... ... @@ -83,29 +83,29 @@ protected:
83 83  
84 84 // Integration kernel for line-of-sight integral
85 85 class halo_kernel_los : public GFunction {
86   - public :
  86 + public :
87 87 halo_kernel_los(const double& scale_radius ,
88 88 const double& halo_distance,
89   - const double& theta ) :
  89 + const double& theta) :
90 90 m_scale_radius(scale_radius),
91 91 m_halo_distance(halo_distance),
92 92 m_theta(theta) {}
93   - double eval( const double& los ) ;
94   - protected :
95   - double m_scale_radius ;
96   - double m_halo_distance ;
97   - double m_theta ;
98   - } ;
  93 + double eval(const double& los);
  94 + protected :
  95 + double m_scale_radius;
  96 + double m_halo_distance;
  97 + double m_theta;
  98 + };
99 99  
100 100 // Protected members
101   - GModelPar m_theta_min ; //!< maximum theta angle
102   - GModelPar m_theta_max ; //!< maximum theta angle
103   - GModelPar m_scale_radius ; //!< scale radius of halo profile
104   - GModelPar m_halo_distance ; //!< distance from earth to halo center
  101 + GModelPar m_theta_min; //!< Minimum theta angle
  102 + GModelPar m_theta_max; //!< Maximum theta angle
  103 + GModelPar m_scale_radius; //!< Scale radius of halo profile
  104 + GModelPar m_halo_distance; //!< Distance from Earth to halo center
105 105  
106 106 // Cached members used for precomputation
107   - mutable double m_last_scale_radius ;
108   - mutable double m_mass_radius ;
  107 + mutable double m_last_scale_radius;
  108 + mutable double m_mass_radius;
109 109 };
110 110  
111 111  
... ...
include/GModelSpatialRadialProfileDMEinasto.hpp
  View file @58af4ad
1 1 /***************************************************************************
2 2 * GModelSpatialRadialProfileDMEinasto.hpp - DM Einasto profile class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2016 by Nathan Kelley-Hoskins *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -59,21 +59,21 @@ public:
59 59 virtual GModelSpatialRadialProfileDMEinasto& operator=(const GModelSpatialRadialProfileDMEinasto& model);
60 60  
61 61 // Implemented pure virtual base class methods
62   - virtual void clear(void);
  62 + virtual void clear(void);
63 63 virtual GModelSpatialRadialProfileDMEinasto* clone(void) const;
64   - virtual std::string classname(void) const;
65   - virtual std::string type(void) const;
66   - virtual double theta_min(void) const;
67   - virtual double theta_max(void) const;
68   - virtual void read(const GXmlElement& xml);
69   - virtual void write(GXmlElement& xml) const;
70   - virtual std::string print(const GChatter& chatter = NORMAL) const;
  64 + virtual std::string classname(void) const;
  65 + virtual std::string type(void) const;
  66 + virtual double theta_min(void) const;
  67 + virtual double theta_max(void) const;
  68 + virtual void read(const GXmlElement& xml);
  69 + virtual void write(GXmlElement& xml) const;
  70 + virtual std::string print(const GChatter& chatter = NORMAL) const;
71 71  
72 72 // Other methods
73   - double scale_radius(void) const ;
74   - void scale_radius(const double& scale_radius ) ;
75   - //double prof_val(const double& theta) ;
76   - double j_factor( const double& minangle, const double& maxangle, const int& npoints ) const ;
  73 + double scale_radius(void) const;
  74 + void scale_radius(const double& scale_radius);
  75 + //double prof_val(const double& theta);
  76 + double j_factor( const double& minangle, const double& maxangle, const int& npoints) const;
77 77  
78 78 protected:
79 79 // Protected methods
... ... @@ -85,33 +85,33 @@ protected:
85 85  
86 86 // Integration kernel for line-of-sight integral
87 87 class halo_kernel_los : public GFunction {
88   - public :
89   - halo_kernel_los(const double& scale_radius ,
  88 + public :
  89 + halo_kernel_los(const double& scale_radius,
90 90 const double& halo_distance,
91   - const double& alpha ,
92   - const double& theta ) :
  91 + const double& alpha,
  92 + const double& theta) :
93 93 m_scale_radius(scale_radius),
94 94 m_halo_distance(halo_distance),
95 95 m_alpha(alpha),
96 96 m_theta(theta) {}
97   - double eval( const double& los ) ;
98   - protected :
99   - double m_scale_radius ;
100   - double m_halo_distance ;
101   - double m_alpha ;
102   - double m_theta ;
103   - } ;
  97 + double eval(const double& los);
  98 + protected :
  99 + double m_scale_radius;
  100 + double m_halo_distance;
  101 + double m_alpha;
  102 + double m_theta;
  103 + };
104 104  
105 105 // Protected members
106   - GModelPar m_scale_radius ; //!< scale radius of halo profile
107   - GModelPar m_halo_distance ; //!< distance from earth to halo center
108   - GModelPar m_alpha ; //!< einasto spatial power index
109   - GModelPar m_theta_min ; //!< minimum theta angle
110   - GModelPar m_theta_max ; //!< maximum theta angle
  106 + GModelPar m_theta_min; //!< Minimum theta angle
  107 + GModelPar m_theta_max; //!< Maximum theta angle
  108 + GModelPar m_scale_radius; //!< Scale radius of halo profile
  109 + GModelPar m_halo_distance; //!< Distance from earth to halo center
  110 + GModelPar m_alpha; //!< Einasto spatial power index
111 111  
112 112 // Cached members used for pre-computation
113   - mutable double m_last_scale_radius ;
114   - mutable double m_mass_radius ;
  113 + mutable double m_last_scale_radius;
  114 + mutable double m_mass_radius;
115 115 };
116 116  
117 117  
... ...
include/GModelSpatialRadialProfileDMZhao.hpp
  View file @58af4ad
1 1 /***************************************************************************
2   - * GModelSpatialRadialProfileDMZhao.hpp - DM Einasto profile class *
  2 + * GModelSpatialRadialProfileDMZhao.hpp - DM Zhao profile class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2016 by Nathan Kelley-Hoskins *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -20,7 +20,7 @@
20 20 ***************************************************************************/
21 21 /**
22 22 * @file GModelSpatialRadialProfileDMZhao.hpp
23   - * @brief Dark Matter Einasto profile model class interface definition
  23 + * @brief Dark Matter Zhao profile model class interface definition
24 24 * @author Nathan Kelley-Hoskins
25 25 */
26 26  
... ... @@ -41,10 +41,10 @@ class GXmlElement;
41 41 /**************************************************************************
42 42 * @class GModelSpatialRadialProfileDMZhao
43 43 *
44   - * @brief Radial Dark Matter Einasto profile source model class
  44 + * @brief Radial Dark Matter Zhao profile source model class
45 45 *
46 46 * This class implements the spatial component of the factorised source
47   - * model for a radial Dark Matter profile, using an Einasto density halo.
  47 + * model for a radial Dark Matter profile, using an Zhao density halo.
48 48 ***************************************************************************/
49 49 class GModelSpatialRadialProfileDMZhao : public GModelSpatialRadialProfile {
50 50  
... ... @@ -59,20 +59,20 @@ public:
59 59 virtual GModelSpatialRadialProfileDMZhao& operator=(const GModelSpatialRadialProfileDMZhao& model);
60 60  
61 61 // Implemented pure virtual base class methods
62   - virtual void clear(void);
  62 + virtual void clear(void);
63 63 virtual GModelSpatialRadialProfileDMZhao* clone(void) const;
64   - virtual std::string classname(void) const;
65   - virtual std::string type(void) const;
66   - virtual double theta_min(void) const;
67   - virtual double theta_max(void) const;
68   - virtual void read(const GXmlElement& xml);
69   - virtual void write(GXmlElement& xml) const;
70   - virtual std::string print(const GChatter& chatter = NORMAL) const;
  64 + virtual std::string classname(void) const;
  65 + virtual std::string type(void) const;
  66 + virtual double theta_min(void) const;
  67 + virtual double theta_max(void) const;
  68 + virtual void read(const GXmlElement& xml);
  69 + virtual void write(GXmlElement& xml) const;
  70 + virtual std::string print(const GChatter& chatter = NORMAL) const;
71 71  
72 72 // Other methods
73   - double scale_radius(void) const ;
74   - void scale_radius(const double& scale_radius ) ;
75   - double prof_val(const double& theta ) ;
  73 + double scale_radius(void) const;
  74 + void scale_radius(const double& scale_radius);
  75 + double prof_val(const double& theta);
76 76  
77 77 protected:
78 78 // Protected methods
... ... @@ -84,38 +84,39 @@ protected:
84 84  
85 85 // Integration kernel for line-of-sight integral
86 86 class halo_kernel_los : public GFunction {
87   - public :
88   - halo_kernel_los(const double& scale_radius ,
  87 + public :
  88 + halo_kernel_los(const double& scale_radius,
89 89 const double& halo_distance,
90   - const double& alpha ,
91   - const double& beta ,
92   - const double& gamma ,
93   - const double& theta ) :
  90 + const double& alpha,
  91 + const double& beta,
  92 + const double& gamma,
  93 + const double& theta) :
94 94 m_scale_radius(scale_radius),
95 95 m_halo_distance(halo_distance),
96 96 m_alpha(alpha),
97 97 m_beta(beta),
98 98 m_gamma(gamma),
99 99 m_theta(theta) {}
100   - double eval( const double& los ) ;
101   - protected :
102   - double m_scale_radius ;
103   - double m_halo_distance ;
104   - double m_alpha ;
105   - double m_beta ;
106   - double m_gamma ;
107   - double m_theta ;
  100 + double eval(const double& los);
  101 + protected :
  102 + double m_scale_radius;
  103 + double m_halo_distance;
  104 + double m_alpha;
  105 + double m_beta;
  106 + double m_gamma;
  107 + double m_theta;
108 108 } ;
109 109  
110 110 // Protected members
111   - GModelPar m_scale_radius ; //!< scale radius of halo profile
112   - GModelPar m_halo_distance ; //!< distance from earth to halo center
113   - GModelPar m_alpha ; //!< power index, inverse transition region width
114   - GModelPar m_beta ; //!< power index, slope at >> m_scale_radius
115   - GModelPar m_gamma ; //!< power index, slope at << m_scale_radius
116   - GModelPar m_theta_min ; //!< minimum theta angle
117   - GModelPar m_theta_max ; //!< maximum theta angle
118   -
  111 + GModelPar m_theta_min; //!< Minimum theta angle
  112 + GModelPar m_theta_max; //!< Maximum theta angle
  113 + GModelPar m_scale_radius; //!< Scale radius of halo profile
  114 + GModelPar m_halo_distance; //!< Distance from earth to halo center
  115 + GModelPar m_alpha; //!< Power index, inverse transition region width
  116 + GModelPar m_beta; //!< Power index, slope at >> m_scale_radius
  117 + GModelPar m_gamma; //!< Power index, slope at << m_scale_radius
  118 +
  119 + // Protected cached members
119 120 mutable double m_last_scale_radius ;
120 121 mutable double m_mass_radius ;
121 122 };
... ...
pyext/GModelSpatialRadialProfile.i
  View file @58af4ad
... ... @@ -55,17 +55,16 @@ public:
55 55 virtual double theta_max(void) const = 0;
56 56  
57 57 // Implemented pure virtual base class methods
58   - virtual double eval(const double& theta,
59   - const GEnergy& energy,
60   - const GTime& time) const;
61   - virtual double eval_gradients(const double& theta,
62   - const GEnergy& energy,
63   - const GTime& time) const;
64   - virtual GSkyDir mc(const GEnergy& energy,
65   - const GTime& time,
66   - GRan& ran) const;
67   - virtual bool contains(const GSkyDir& dir,
68   - const double& margin = 0.0) const;
  58 + virtual double eval(const double& theta,
  59 + const GEnergy& energy,
  60 + const GTime& time,
  61 + const bool& gradients = false) const;
  62 + virtual GSkyDir mc(const GEnergy& energy,
  63 + const GTime& time,
  64 + GRan& ran) const;
  65 + virtual bool contains(const GSkyDir& dir,
  66 + const double& margin = 0.0) const;
  67 + virtual GSkyRegion* region(void) const;
69 68  
70 69 // Implement other methods
71 70 int num_nodes(void) const;
... ...
pyext/GModelSpatialRadialProfileDMBurkert.i
  View file @58af4ad
1 1 /***************************************************************************
2   - * GModelSpatialRadialProfileDMBurkert.i - Einasto radial profile class *
  2 + * GModelSpatialRadialProfileDMBurkert.i - Burkert radial profile class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2016 by Nathan Kelley-Hoskins *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -20,7 +20,7 @@
20 20 ***************************************************************************/
21 21 /**
22 22 * @file GModelSpatialRadialProfileDMBurkert.hpp
23   - * @brief Radial Dark Matter halo for Einasto Density Profile
  23 + * @brief Radial Dark Matter halo for Burkert Density Profile
24 24 * @author Nathan Kelley-Hoskins
25 25 */
26 26 %{
... ... @@ -32,10 +32,10 @@
32 32 /**************************************************************************
33 33 * @class GModelSpatialRadialProfileDMBurkert
34 34 *
35   - * @brief Radial DM Einasto profile source model class
  35 + * @brief Radial DM Burkert profile source model class
36 36 *
37 37 * This class implements the spatial component of the factorised source
38   - * model for a Dark Matter Einasto halo radial profile.
  38 + * model for a Dark Matter Burkert halo radial profile.
39 39 ***************************************************************************/
40 40 class GModelSpatialRadialProfileDMBurkert : public GModelSpatialRadialProfile {
41 41  
... ... @@ -47,25 +47,23 @@ public:
47 47 virtual ~GModelSpatialRadialProfileDMBurkert(void);
48 48  
49 49 // Implemented pure virtual base class methods
50   - virtual void clear(void);
  50 + virtual void clear(void);
51 51 virtual GModelSpatialRadialProfileDMBurkert* clone(void) const;
52   - virtual std::string classname(void) const;
53   - virtual std::string type(void) const;
54   - virtual double theta_min(void) const;
55   - virtual double theta_max(void) const;
56   - virtual void read(const GXmlElement& xml);
57   - virtual void write(GXmlElement& xml) const;
  52 + virtual std::string classname(void) const;
  53 + virtual std::string type(void) const;
  54 + virtual double theta_min(void) const;
  55 + virtual double theta_max(void) const;
  56 + virtual void read(const GXmlElement& xml);
  57 + virtual void write(GXmlElement& xml) const;
  58 +
  59 + // Other methods
  60 + double scale_radius(void) const;
  61 + void scale_radius(const double& scale_radius);
58 62 };
59 63  
60 64  
61 65 /***********************************************************************//**
62 66 * @brief GModelSpatialRadialGauss class extension
63   - *
64   - * The eval(GSkyDir&) and eval_gradients(GSkyDir&) need to be defined in the
65   - * extension to force swig to build also the interface for these methods that
66   - * are implemented in the base class only. It's not clear to me why these
67   - * methods are not inherited automatically. Maybe this could also be handled
68   - * by a %typemap(typecheck) construct.
69 67 ***************************************************************************/
70 68 %extend GModelSpatialRadialProfileDMBurkert {
71 69 GModelSpatialRadialProfileDMBurkert copy() {
... ...
pyext/GModelSpatialRadialProfileDMEinasto.i
  View file @58af4ad
1 1 /***************************************************************************
2 2 * GModelSpatialRadialProfileDMEinasto.i - Einasto radial profile class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2016 by Nathan Kelley-Hoskins *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -47,31 +47,25 @@ public:
47 47 virtual ~GModelSpatialRadialProfileDMEinasto(void);
48 48  
49 49 // Implemented pure virtual base class methods
50   - virtual void clear(void);
  50 + virtual void clear(void);
51 51 virtual GModelSpatialRadialProfileDMEinasto* clone(void) const;
52   - virtual std::string classname(void) const;
53   - virtual std::string type(void) const;
54   - virtual double theta_min(void) const;
55   - virtual double theta_max(void) const;
56   - virtual void read(const GXmlElement& xml);
57   - virtual void write(GXmlElement& xml) const;
  52 + virtual std::string classname(void) const;
  53 + virtual std::string type(void) const;
  54 + virtual double theta_min(void) const;
  55 + virtual double theta_max(void) const;
  56 + virtual void read(const GXmlElement& xml);
  57 + virtual void write(GXmlElement& xml) const;
58 58  
59   - // other methods
60   - double scale_radius(void) const ;
61   - void scale_radius(const double& scale_radius ) ;
62   - //double prof_val(const double& theta) ;
63   -
  59 + // Other methods
  60 + double scale_radius(void) const;
  61 + void scale_radius(const double& scale_radius);
  62 + //double prof_val(const double& theta);
  63 + double j_factor( const double& minangle, const double& maxangle, const int& npoints) const;
64 64 };
65 65  
66 66  
67 67 /***********************************************************************//**
68 68 * @brief GModelSpatialRadialGauss class extension
69   - *
70   - * The eval(GSkyDir&) and eval_gradients(GSkyDir&) need to be defined in the
71   - * extension to force swig to build also the interface for these methods that
72   - * are implemented in the base class only. It's not clear to me why these
73   - * methods are not inherited automatically. Maybe this could also be handled
74   - * by a %typemap(typecheck) construct.
75 69 ***************************************************************************/
76 70 %extend GModelSpatialRadialProfileDMEinasto {
77 71 GModelSpatialRadialProfileDMEinasto copy() {
... ...
pyext/GModelSpatialRadialProfileDMZhao.i
  View file @58af4ad
1 1 /***************************************************************************
2   - * GModelSpatialRadialProfileDMZhao.i - Einasto radial profile class *
  2 + * GModelSpatialRadialProfileDMZhao.i - Zhao radial profile class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2016 by Nathan Kelley-Hoskins *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -20,7 +20,7 @@
20 20 ***************************************************************************/
21 21 /**
22 22 * @file GModelSpatialRadialProfileDMZhao.hpp
23   - * @brief Radial Dark Matter halo for Einasto Density Profile
  23 + * @brief Radial Dark Matter halo for Zhao Density Profile
24 24 * @author Nathan Kelley-Hoskins
25 25 */
26 26 %{
... ... @@ -32,7 +32,7 @@
32 32 /**************************************************************************
33 33 * @class GModelSpatialRadialProfileDMZhao
34 34 *
35   - * @brief Radial DM Einasto profile source model class
  35 + * @brief Radial DM Zhao profile source model class
36 36 *
37 37 * This class implements the spatial component of the factorised source
38 38 * model for a Dark Matter Zhao halo radial profile.
... ... @@ -47,26 +47,24 @@ public:
47 47 virtual ~GModelSpatialRadialProfileDMZhao(void);
48 48  
49 49 // Implemented pure virtual base class methods
50   - virtual void clear(void);
  50 + virtual void clear(void);
51 51 virtual GModelSpatialRadialProfileDMZhao* clone(void) const;
52   - virtual std::string classname(void) const;
53   - virtual std::string type(void) const;
54   - virtual double theta_min(void) const;
55   - virtual double theta_max(void) const;
56   - virtual void read(const GXmlElement& xml);
57   - virtual void write(GXmlElement& xml) const;
58   - double prof_val(const double& theta ) ;
  52 + virtual std::string classname(void) const;
  53 + virtual std::string type(void) const;
  54 + virtual double theta_min(void) const;
  55 + virtual double theta_max(void) const;
  56 + virtual void read(const GXmlElement& xml);
  57 + virtual void write(GXmlElement& xml) const;
  58 +
  59 + // Other methods
  60 + double scale_radius(void) const;
  61 + void scale_radius(const double& scale_radius);
  62 + double prof_val(const double& theta);
59 63 };
60 64  
61 65  
62 66 /***********************************************************************//**
63 67 * @brief GModelSpatialRadialGauss class extension
64   - *
65   - * The eval(GSkyDir&) and eval_gradients(GSkyDir&) need to be defined in the
66   - * extension to force swig to build also the interface for these methods that
67   - * are implemented in the base class only. It's not clear to me why these
68   - * methods are not inherited automatically. Maybe this could also be handled
69   - * by a %typemap(typecheck) construct.
70 68 ***************************************************************************/
71 69 %extend GModelSpatialRadialProfileDMZhao {
72 70 GModelSpatialRadialProfileDMZhao copy() {
... ...
src/model/GModelSpatialRadialProfile.cpp
  View file @58af4ad
... ... @@ -180,6 +180,7 @@ GModelSpatialRadialProfile&amp; GModelSpatialRadialProfile::operator=(const GModelSp
180 180 * @param[in] theta Angular distance from model centre (radians).
181 181 * @param[in] energy Photon energy.
182 182 * @param[in] time Photon arrival time.
  183 + * @param[in] gradients Compute gradients?
183 184 * @return Model value.
184 185 *
185 186 * Evaluate the radial profile model for a given angular distance @p theta
... ... @@ -188,7 +189,8 @@ GModelSpatialRadialProfile&amp; GModelSpatialRadialProfile::operator=(const GModelSp
188 189 ***************************************************************************/
189 190 double GModelSpatialRadialProfile::eval(const double& theta,
190 191 const GEnergy& energy,
191   - const GTime& time) const
  192 + const GTime& time,
  193 + const bool& gradients) const
192 194 {
193 195 // Get pre-computation cache index
194 196 int icache = cache_index();
... ... @@ -217,29 +219,6 @@ double GModelSpatialRadialProfile::eval(const double&amp; theta,
217 219  
218 220  
219 221 /***********************************************************************//**
220   - * @brief Evaluate function and gradients (in units of sr^-1)
221   - *
222   - * @param[in] theta Angular distance from model centre (radians).
223   - * @param[in] energy Photon energy.
224   - * @param[in] time Photon arrival time.
225   - * @return Model value.
226   - *
227   - * Evaluates the function value. No gradient computation is implemented as
228   - * radial models will be convolved with the instrument response and thus
229   - * require the numerical computation of the derivatives.
230   - *
231   - * See the eval() method for more information.
232   - ***************************************************************************/
233   -double GModelSpatialRadialProfile::eval_gradients(const double& theta,
234   - const GEnergy& energy,
235   - const GTime& time) const
236   -{
237   - // Return value
238   - return (eval(theta, energy, time));
239   -}
240   -
241   -
242   -/***********************************************************************//**
243 222 * @brief Return MC sky direction
244 223 *
245 224 * @param[in] energy Photon energy.
... ... @@ -315,6 +294,7 @@ void GModelSpatialRadialProfile::init_members(void)
315 294 {
316 295 // Initialise members
317 296 m_num_nodes = 100;
  297 + m_region.clear();
318 298  
319 299 // Initialise pre-computation cache
320 300 m_profile.clear();
... ... @@ -334,6 +314,7 @@ void GModelSpatialRadialProfile::copy_members(const GModelSpatialRadialProfile&amp;
334 314 // Copy members
335 315 m_num_nodes = model.m_num_nodes;
336 316 m_profile = model.m_profile;
  317 + m_region = model.m_region;
337 318  
338 319 // Return
339 320 return;
... ... @@ -466,3 +447,19 @@ int GModelSpatialRadialProfile::cache_index(void) const
466 447 // Return index
467 448 return (index);
468 449 }
  450 +
  451 +
  452 +/***********************************************************************//**
  453 + * @brief Set boundary sky region
  454 + ***************************************************************************/
  455 +void GModelSpatialRadialProfile::set_region(void) const
  456 +{
  457 + // Set sky region centre to disk centre
  458 + m_region.centre(m_ra.value(), m_dec.value());
  459 +
  460 + // Set sky region radius to maximum theta angle
  461 + m_region.radius(theta_max()*gammalib::rad2deg);
  462 +
  463 + // Return
  464 + return;
  465 +}
... ...
src/model/GModelSpatialRadialProfileDMBurkert.cpp
  View file @58af4ad
1 1 /***************************************************************************
2   - * GModelSpatialRadialProfileDMBurkert.cpp - DMBurkert radial profile class *
  2 + * GModelSpatialRadialProfileDMBurkert.cpp - Burkert radial profile class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2016 by Nathan Kelley-Hoskins *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -20,7 +20,7 @@
20 20 ***************************************************************************/
21 21 /**
22 22 * @file GModelSpatialRadialProfileDMBurkert.cpp
23   - * @brief Radial DM Einasto profile model class implementation
  23 + * @brief Radial DM Burkert profile model class implementation
24 24 * @author Nathan Kelley-Hoskins
25 25 */
26 26  
... ... @@ -40,7 +40,7 @@
40 40  
41 41 /* __ Globals ____________________________________________________________ */
42 42 const GModelSpatialRadialProfileDMBurkert g_radial_disk_seed;
43   -const GModelSpatialRegistry g_radial_disk_registry(&g_radial_disk_seed);
  43 +const GModelSpatialRegistry g_radial_disk_registry(&g_radial_disk_seed);
44 44  
45 45 /* __ Method name definitions ____________________________________________ */
46 46 #define G_READ "GModelSpatialRadialProfileDMBurkert::read(GXmlElement&)"
... ... @@ -65,7 +65,7 @@ const GModelSpatialRegistry g_radial_disk_registry(&amp;g_radial_disk_seed
65 65 * Constructs empty radial DMBurkert profile
66 66 ***************************************************************************/
67 67 GModelSpatialRadialProfileDMBurkert::GModelSpatialRadialProfileDMBurkert(void) :
68   - GModelSpatialRadialProfile()
  68 + GModelSpatialRadialProfile()
69 69 {
70 70 // Initialise members
71 71 init_members();
... ... @@ -85,7 +85,7 @@ GModelSpatialRadialProfileDMBurkert::GModelSpatialRadialProfileDMBurkert(void) :
85 85 * expected structure of the XML element.
86 86 ***************************************************************************/
87 87 GModelSpatialRadialProfileDMBurkert::GModelSpatialRadialProfileDMBurkert(const GXmlElement& xml) :
88   - GModelSpatialRadialProfile()
  88 + GModelSpatialRadialProfile()
89 89 {
90 90 // Initialise members
91 91 init_members();
... ... @@ -106,7 +106,7 @@ GModelSpatialRadialProfileDMBurkert::GModelSpatialRadialProfileDMBurkert(const G
106 106 * Copies radial DMBurkert profile model from another radial profile model.
107 107 ***************************************************************************/
108 108 GModelSpatialRadialProfileDMBurkert::GModelSpatialRadialProfileDMBurkert(const GModelSpatialRadialProfileDMBurkert& model) :
109   - GModelSpatialRadialProfile(model)
  109 + GModelSpatialRadialProfile(model)
110 110 {
111 111 // Initialise members
112 112 init_members();
... ... @@ -227,7 +227,7 @@ double GModelSpatialRadialProfileDMBurkert::theta_min(void) const
227 227 update();
228 228  
229 229 // Return value
230   - return m_theta_min.value() ;
  230 + return m_theta_min.value();
231 231 }
232 232  
233 233 /***********************************************************************//**
... ... @@ -237,31 +237,32 @@ double GModelSpatialRadialProfileDMBurkert::theta_min(void) const
237 237 ***************************************************************************/
238 238 double GModelSpatialRadialProfileDMBurkert::theta_max(void) const
239 239 {
240   -
  240 + // Update
241 241 update();
242 242  
243 243 // initialize value
244   - double theta = 0.0 ;
  244 + double theta = 0.0;
245 245  
246   - // if earth is within the significant radius, then we must integrate
247   - // the entire profile
  246 + // If Earth is within the significant radius, then we must integrate
  247 + // the entire profile ...
248 248 if ( m_halo_distance.value() < m_mass_radius ) {
249   - theta = gammalib::pi ;
  249 + theta = gammalib::pi;
  250 + }
250 251  
251   - // if the halo is far enough away (further than the significant radius)
252   - // then we just need to deal with the angles within the sphere of the
253   - // significant radius.
254   - } else {
255   - theta = std::atan( m_mass_radius / m_halo_distance.value() ) ;
  252 + // ... otherwise, if the halo is far enough away (further than the
  253 + // significant radius) then we just need to deal with the angles within
  254 + // the sphere of the significant radius
  255 + else {
  256 + theta = std::atan( m_mass_radius / m_halo_distance.value() ) ;
256 257 }
257 258  
258   - // always chose the lesser of ( mass_radius theta, theta_max )
259   - if ( m_theta_max.value() * gammalib::deg2rad < theta ) {
260   - theta = m_theta_max.value() * gammalib::deg2rad ;
  259 + // Always chose the lesser of ( mass_radius theta, theta_max )
  260 + if (m_theta_max.value() * gammalib::deg2rad < theta ) {
  261 + theta = m_theta_max.value() * gammalib::deg2rad;
261 262 }
262 263  
263 264 // Return value
264   - return theta ;
  265 + return theta;
265 266 }
266 267  
267 268  
... ... @@ -274,23 +275,23 @@ double GModelSpatialRadialProfileDMBurkert::theta_max(void) const
274 275 * The XML element shall have either the format
275 276 *
276 277 * <spatialModel type="DMBurkertProfile">
277   - * <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
278   - * <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
279   - * <parameter name="Sigma" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
  278 + * <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
  279 + * <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
  280 + * <parameter name="Sigma" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
280 281 * <parameter name="Scale Radius" scale="1.0" value="21.5" min="1.0e-6" free="1"/>
281 282 * <parameter name="Halo Distance" scale="1.0" value="7.94" min="1.0e-6" free="1"/>
282   - * <parameter name="Alpha" scale="1.0" value="0.17" min="0.01" max="10" free="1"/>
  283 + * <parameter name="Theta Max" scale="1.0" value="0.17" min="0.01" max="10" free="1"/>
283 284 * </spatialModel>
284 285 *
285 286 * or
286 287 *
287 288 * <spatialModel type="DMBurkertProfile">
288   - * <parameter name="GLON" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
289   - * <parameter name="GLAT" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
290   - * <parameter name="Sigma" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
  289 + * <parameter name="GLON" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
  290 + * <parameter name="GLAT" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
  291 + * <parameter name="Sigma" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
291 292 * <parameter name="Scale Radius" scale="1.0" value="21.5" min="1.0e-6" free="1"/>
292 293 * <parameter name="Halo Distance" scale="1.0" value="7.94" min="1.0e-6" free="1"/>
293   - * <parameter name="Alpha" scale="1.0" value="0.17" min="0.01" max="10" free="1"/>
  294 + * <parameter name="Theta Max" scale="1.0" value="0.17" min="0.01" max="10" free="1"/>
294 295 * </spatialModel>
295 296 ***************************************************************************/
296 297 void GModelSpatialRadialProfileDMBurkert::read(const GXmlElement& xml)
... ... @@ -323,9 +324,12 @@ void GModelSpatialRadialProfileDMBurkert::read(const GXmlElement&amp; xml)
323 324 * The XML element will have the format
324 325 *
325 326 * <spatialModel type="DMBurkertProfile">
326   - * <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
327   - * <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
328   - * <parameter name="Sigma" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
  327 + * <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
  328 + * <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
  329 + * <parameter name="Sigma" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
  330 + * <parameter name="Scale Radius" scale="1.0" value="21.5" min="1.0e-6" free="1"/>
  331 + * <parameter name="Halo Distance" scale="1.0" value="7.94" min="1.0e-6" free="1"/>
  332 + * <parameter name="Theta Max" scale="1.0" value="0.17" min="0.01" max="10" free="1"/>
329 333 * </spatialModel>
330 334 ***************************************************************************/
331 335 void GModelSpatialRadialProfileDMBurkert::write(GXmlElement& xml) const
... ... @@ -419,9 +423,9 @@ void GModelSpatialRadialProfileDMBurkert::init_members(void)
419 423 m_theta_min.clear();
420 424 m_theta_min.name("Theta Min");
421 425 m_theta_min.unit("degrees");
422   - m_theta_min.value( 1.0e-6 ); //
423   - m_theta_min.min(1.0e-10); // arbitrarily chosen
424   - m_theta_min.fix(); // should always be fixed!
  426 + m_theta_min.value(1.0e-6);
  427 + m_theta_min.min(1.0e-10); // arbitrarily chosen
  428 + m_theta_min.fix(); // should always be fixed!
425 429 m_theta_min.scale(1.0);
426 430 m_theta_min.gradient(0.0);
427 431 m_theta_min.has_grad(false); // Radial components never have gradients
... ... @@ -430,18 +434,18 @@ void GModelSpatialRadialProfileDMBurkert::init_members(void)
430 434 m_theta_max.clear();
431 435 m_theta_max.name("Theta Max");
432 436 m_theta_max.unit("degrees");
433   - m_theta_max.value( 180.0 ); // can only go from halo center to opposite halo center
434   - m_theta_max.min(1.0e-6); // arbitrarily chosen
435   - m_theta_max.fix(); // should always be fixed!
  437 + m_theta_max.value(180.0); // can only go from halo center to opposite halo center
  438 + m_theta_max.min(1.0e-6); // arbitrarily chosen
  439 + m_theta_max.fix(); // should always be fixed!
436 440 m_theta_max.scale(1.0);
437 441 m_theta_max.gradient(0.0);
438 442 m_theta_max.has_grad(false); // Radial components never have gradients
439 443  
440 444 // Set parameter pointer(s)
441   - m_pars.push_back(&m_scale_radius );
  445 + m_pars.push_back(&m_scale_radius);
442 446 m_pars.push_back(&m_halo_distance);
443   - m_pars.push_back(&m_theta_max );
444   - m_pars.push_back(&m_theta_min );
  447 + m_pars.push_back(&m_theta_max);
  448 + m_pars.push_back(&m_theta_min);
445 449  
446 450 // Initialize precomputation cache. Note that zero values flag
447 451 // uninitialised, as a zero radius is not meaningful
... ... @@ -465,14 +469,14 @@ void GModelSpatialRadialProfileDMBurkert::copy_members(const GModelSpatialRadial
465 469 // Copy members. We do not have to push back the members on the parameter
466 470 // stack as this should have been done by init_members() that was called
467 471 // before. Otherwise we would have sigma twice on the stack.
468   - m_scale_radius = model.m_scale_radius ;
469   - m_halo_distance = model.m_halo_distance ;
470   - m_theta_min = model.m_theta_min ;
471   - m_theta_max = model.m_theta_max ;
  472 + m_scale_radius = model.m_scale_radius;
  473 + m_halo_distance = model.m_halo_distance;
  474 + m_theta_min = model.m_theta_min;
  475 + m_theta_max = model.m_theta_max;
472 476  
473 477 // copy cache values
474   - m_last_scale_radius = model.m_last_scale_radius ;
475   - m_mass_radius = model.m_mass_radius ;
  478 + m_last_scale_radius = model.m_last_scale_radius;
  479 + m_mass_radius = model.m_mass_radius;
476 480  
477 481 // Return
478 482 return;
... ... @@ -503,25 +507,25 @@ double GModelSpatialRadialProfileDMBurkert::profile_value(const double&amp; theta) c
503 507 update();
504 508  
505 509 // initialize integral value
506   - double value = 0.0 ;
  510 + double value = 0.0;
507 511  
508 512 // Set up integration limits
509   - double los_min = m_halo_distance.value() - m_mass_radius ;
510   - double los_max = m_halo_distance.value() + m_mass_radius ;
  513 + double los_min = m_halo_distance.value() - m_mass_radius;
  514 + double los_max = m_halo_distance.value() + m_mass_radius;
511 515  
512   - // handle case where observer is within halo mass radius
  516 + // Handle case where observer is within halo mass radius
513 517 if ( los_min < 0.0 ) {
514   - los_min = 0.0 ;
  518 + los_min = 0.0 ;
515 519 }
516 520  
517 521 // Set up integral
518   - halo_kernel_los integrand( m_scale_radius.value(),
519   - m_halo_distance.value(),
520   - theta ) ;
521   - GIntegral integral(&integrand) ;
  522 + halo_kernel_los integrand(m_scale_radius.value(),
  523 + m_halo_distance.value(),
  524 + theta);
  525 + GIntegral integral(&integrand);
522 526  
523 527 // Compute value
524   - value = integral.romberg( los_min, los_max ) ;
  528 + value = integral.romberg(los_min, los_max);
525 529  
526 530 //double tm = theta_max() ;
527 531 //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&amp; theta) c
530 534 return value;
531 535 }
532 536  
  537 +
533 538 /***********************************************************************//**
534 539 * @brief Kernel for halo density profile squared
535 540 *
... ... @@ -558,52 +563,55 @@ double GModelSpatialRadialProfileDMBurkert::profile_value(const double&amp; theta) c
558 563 * @return unit
559 564 *
560 565 ***************************************************************************/
561   -double GModelSpatialRadialProfileDMBurkert::halo_kernel_los::eval( const double &los )
  566 +double GModelSpatialRadialProfileDMBurkert::halo_kernel_los::eval(const double &los)
562 567 {
  568 + // PLEASE ADD COMMENTS
  569 + double r = 0.0;
  570 + r = los * los;
  571 + r += m_halo_distance * m_halo_distance;
  572 + r -= 2.0 * los * m_halo_distance * std::cos(m_theta);
  573 + r = std::sqrt(r);
563 574  
564   - double r = 0.0 ;
565   - r = los * los ;
566   - r += m_halo_distance * m_halo_distance ;
567   - r -= 2.0 * los * m_halo_distance * std::cos(m_theta) ;
568   - r = sqrt(r) ;
569   -
570   - double bot = ( m_scale_radius + r ) ;
571   - bot *= (m_scale_radius*m_scale_radius) + (r*r) ;
  575 + double bot = (m_scale_radius + r);
  576 + bot *= (m_scale_radius*m_scale_radius) + (r*r);
572 577  
573   - double f = m_scale_radius * m_scale_radius * m_scale_radius ;
574   - f /= bot ;
  578 + double f = m_scale_radius * m_scale_radius * m_scale_radius;
  579 + f /= bot;
575 580  
576   - // squared, for annihilating dm
577   - // would just be f if it was decaying dm
578   - f = f * f ;
  581 + // squared, for annihilating dm
  582 + // would just be f if it was decaying dm
  583 + f = f * f;
579 584  
580   - //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;
581   - return f;
582   -
  585 + //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;
  586 +
  587 + // Return function value
  588 + return f;
583 589 }
584 590  
  591 +
585 592 /***********************************************************************//**
586 593 * @brief Update precomputation cache
587 594 *
588 595 * Computes the mass_radius calculation, determining the radius around
589   - * the halo that contains 99.99% of the mass. For an einasto halo profile,
590   - * this is just 80.0 * scale_radius .
591   - *
  596 + * the halo that contains 99.99% of the mass. For an Burket halo profile,
  597 + * this is just 80.0 * scale_radius.
592 598 ***************************************************************************/
593 599 void GModelSpatialRadialProfileDMBurkert::update() const
594 600 {
595   - // Update if scale radius has changed
596   - if ( m_last_scale_radius != scale_radius() ) {
  601 + // Update if scale radius has changed
  602 + if (m_last_scale_radius != scale_radius()) {
597 603  
598   - // Store last values
599   - m_last_scale_radius = scale_radius() ;
600   -
601   - // perform precomputations
602   - // set the mass radius to 80*scale_radius, meaning
603   - // 99.99% of the mass is contained within the mass radius,
604   - // and integration only needs to worry about whats inside this radius.
605   - m_mass_radius = 80.0 * scale_radius() ;
  604 + // Store last values
  605 + m_last_scale_radius = scale_radius();
  606 +
  607 + // perform precomputations
  608 + // set the mass radius to 80*scale_radius, meaning
  609 + // 99.99% of the mass is contained within the mass radius,
  610 + // and integration only needs to worry about whats inside this radius.
  611 + m_mass_radius = 80.0 * scale_radius();
606 612  
607   - }
  613 + }
  614 +
  615 + // Return
  616 + return;
608 617 }
609   -
... ...
src/model/GModelSpatialRadialProfileDMEinasto.cpp
  View file @58af4ad
1 1 /***************************************************************************
2   - * GModelSpatialRadialProfileDMEinasto.cpp - DMEinasto radial profile class *
  2 + * GModelSpatialRadialProfileDMEinasto.cpp - Einasto radial profile class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2016 by Nathan Kelley-Hoskins *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -40,7 +40,7 @@
40 40  
41 41 /* __ Globals ____________________________________________________________ */
42 42 const GModelSpatialRadialProfileDMEinasto g_radial_disk_seed;
43   -const GModelSpatialRegistry g_radial_disk_registry(&g_radial_disk_seed);
  43 +const GModelSpatialRegistry g_radial_disk_registry(&g_radial_disk_seed);
44 44  
45 45 /* __ Method name definitions ____________________________________________ */
46 46 #define G_READ "GModelSpatialRadialProfileDMEinasto::read(GXmlElement&)"
... ... @@ -65,7 +65,7 @@ const GModelSpatialRegistry g_radial_disk_registry(&amp;g_radial_disk_seed
65 65 * Constructs empty radial DMEinasto profile
66 66 ***************************************************************************/
67 67 GModelSpatialRadialProfileDMEinasto::GModelSpatialRadialProfileDMEinasto(void) :
68   - GModelSpatialRadialProfile()
  68 + GModelSpatialRadialProfile()
69 69 {
70 70 // Initialise members
71 71 init_members();
... ... @@ -85,7 +85,7 @@ GModelSpatialRadialProfileDMEinasto::GModelSpatialRadialProfileDMEinasto(void) :
85 85 * expected structure of the XML element.
86 86 ***************************************************************************/
87 87 GModelSpatialRadialProfileDMEinasto::GModelSpatialRadialProfileDMEinasto(const GXmlElement& xml) :
88   - GModelSpatialRadialProfile()
  88 + GModelSpatialRadialProfile()
89 89 {
90 90 // Initialise members
91 91 init_members();
... ... @@ -106,7 +106,7 @@ GModelSpatialRadialProfileDMEinasto::GModelSpatialRadialProfileDMEinasto(const G
106 106 * Copies radial DMEinasto profile model from another radial profile model.
107 107 ***************************************************************************/
108 108 GModelSpatialRadialProfileDMEinasto::GModelSpatialRadialProfileDMEinasto(const GModelSpatialRadialProfileDMEinasto& model) :
109   - GModelSpatialRadialProfile(model)
  109 + GModelSpatialRadialProfile(model)
110 110 {
111 111 // Initialise members
112 112 init_members();
... ... @@ -227,7 +227,7 @@ double GModelSpatialRadialProfileDMEinasto::theta_min(void) const
227 227 update();
228 228  
229 229 // Return value
230   - return m_theta_min.value() ;
  230 + return m_theta_min.value();
231 231 }
232 232  
233 233 /***********************************************************************//**
... ... @@ -241,27 +241,28 @@ double GModelSpatialRadialProfileDMEinasto::theta_max(void) const
241 241 // update precomputation cache
242 242 update();
243 243  
244   - double theta = 0.0 ;
  244 + double theta = 0.0;
245 245  
246   - // if earth is within the significant radius, then theta_max must
247   - // contain the entire profile (180deg)
248   - if ( m_halo_distance.value() < m_mass_radius ) {
249   - theta = gammalib::pi ;
  246 + // If Earth is within the significant radius, then theta_max must
  247 + // contain the entire profile (180deg) ...
  248 + if (m_halo_distance.value() < m_mass_radius) {
  249 + theta = gammalib::pi;
  250 + }
250 251  
251   - // if the halo is far enough away (further than the mass radius)
252   - // then we just need to deal with the angles within the sphere of the
253   - // significant radius.
254   - } else {
255   - theta = std::atan( m_mass_radius / m_halo_distance.value() ) ;
  252 + // ... otherwise if the halo is far enough away (further than the mass
  253 + // radius) then we just need to deal with the angles within the sphere
  254 + // of the significant radius.
  255 + else {
  256 + theta = std::atan(m_mass_radius / m_halo_distance.value());
256 257 }
257 258  
258   - // always chose the lesser of ( mass_radius theta, theta_max )
259   - if ( m_theta_max.value() * gammalib::deg2rad < theta ) {
260   - theta = m_theta_max.value() * gammalib::deg2rad ;
  259 + // Always chose the lesser of ( mass_radius theta, theta_max )
  260 + if (m_theta_max.value() * gammalib::deg2rad < theta) {
  261 + theta = m_theta_max.value() * gammalib::deg2rad;
261 262 }
262 263  
263 264 // Return value
264   - return theta ;
  265 + return theta;
265 266 }
266 267  
267 268  
... ... @@ -436,7 +437,7 @@ void GModelSpatialRadialProfileDMEinasto::init_members(void)
436 437 m_theta_min.clear();
437 438 m_theta_min.name("Theta Min");
438 439 m_theta_min.unit("degrees");
439   - m_theta_min.value( 180.0 ); // can only go from halo center to opposite halo center
  440 + m_theta_min.value(180.0); // can only go from halo center to opposite halo center
440 441 m_theta_min.min(1.0e-6); // arbitrarily chosen
441 442 m_theta_min.fix(); // should always be fixed!
442 443 m_theta_min.scale(1.0);
... ... @@ -447,7 +448,7 @@ void GModelSpatialRadialProfileDMEinasto::init_members(void)
447 448 m_theta_max.clear();
448 449 m_theta_max.name("Theta Max");
449 450 m_theta_max.unit("degrees");
450   - m_theta_max.value( 1.0e-6 ); // can only go from halo center to opposite halo center
  451 + m_theta_max.value(1.0e-6); // can only go from halo center to opposite halo center
451 452 m_theta_max.min(1.0e-10); // arbitrarily chosen
452 453 m_theta_max.fix(); // should always be fixed!
453 454 m_theta_max.scale(1.0);
... ... @@ -455,16 +456,16 @@ void GModelSpatialRadialProfileDMEinasto::init_members(void)
455 456 m_theta_max.has_grad(false); // Radial components never have gradients
456 457  
457 458 // Set parameter pointer(s)
458   - m_pars.push_back(&m_scale_radius );
  459 + m_pars.push_back(&m_scale_radius);
459 460 m_pars.push_back(&m_halo_distance);
460   - m_pars.push_back(&m_alpha );
461   - m_pars.push_back(&m_theta_min );
462   - m_pars.push_back(&m_theta_max );
  461 + m_pars.push_back(&m_alpha);
  462 + m_pars.push_back(&m_theta_min);
  463 + m_pars.push_back(&m_theta_max);
463 464  
464 465 // Initialize precomputation cache. Note that zero values flag
465 466 // uninitialised, as a zero radius is not meaningful
466   - m_last_scale_radius = 0.0 ;
467   - m_mass_radius = 0.0 ;
  467 + m_last_scale_radius = 0.0;
  468 + m_mass_radius = 0.0;
468 469  
469 470 // Return
470 471 return;
... ... @@ -483,15 +484,15 @@ void GModelSpatialRadialProfileDMEinasto::copy_members(const GModelSpatialRadial
483 484 // Copy members. We do not have to push back the members on the parameter
484 485 // stack as this should have been done by init_members() that was called
485 486 // before. Otherwise we would have sigma twice on the stack.
486   - m_scale_radius = model.m_scale_radius ;
487   - m_halo_distance = model.m_halo_distance ;
488   - m_alpha = model.m_alpha ;
489   - m_theta_min = model.m_theta_min ;
490   - m_theta_max = model.m_theta_max ;
  487 + m_scale_radius = model.m_scale_radius;
  488 + m_halo_distance = model.m_halo_distance;
  489 + m_alpha = model.m_alpha;
  490 + m_theta_min = model.m_theta_min;
  491 + m_theta_max = model.m_theta_max;
491 492  
492 493 // copy cache values
493   - m_last_scale_radius = model.m_last_scale_radius ;
494   - m_mass_radius = model.m_mass_radius ;
  494 + m_last_scale_radius = model.m_last_scale_radius;
  495 + m_mass_radius = model.m_mass_radius;
495 496  
496 497 // Return
497 498 return;
... ... @@ -521,37 +522,37 @@ double GModelSpatialRadialProfileDMEinasto::profile_value(const double&amp; theta) c
521 522 update();
522 523  
523 524 // initialize integral value
524   - double value = 0.0 ;
  525 + double value = 0.0;
525 526  
526 527 // Set up integration limits
527   - double los_min = m_halo_distance.value() - m_mass_radius ;
528   - double los_max = m_halo_distance.value() + m_mass_radius ;
  528 + double los_min = m_halo_distance.value() - m_mass_radius;
  529 + double los_max = m_halo_distance.value() + m_mass_radius;
529 530  
530 531 // handle case where observer is within halo mass radius
531   - if ( los_min < 0.0 ) {
532   - los_min = 0.0 ;
  532 + if (los_min < 0.0) {
  533 + los_min = 0.0;
533 534 }
534 535  
535 536 // Set up integral
536   - halo_kernel_los integrand( m_scale_radius.value(),
537   - m_halo_distance.value(),
538   - m_alpha.value(),
539   - theta ) ;
540   - GIntegral integral(&integrand) ;
541   - integral.max_iter( 30 ) ;
  537 + halo_kernel_los integrand(m_scale_radius.value(),
  538 + m_halo_distance.value(),
  539 + m_alpha.value(),
  540 + theta);
  541 + GIntegral integral(&integrand);
  542 + integral.max_iter(30);
542 543  
543 544 // also play with eps() or max_iter(), bounds()
544 545  
545 546 // Set up integration boundaries
546 547 // As there is usually an infinity at the halo center, this splits
547 548 // the integral at the m_halo_distance.
548   - std::vector<double> bounds ;
549   - bounds.push_back( los_min ) ;
550   - bounds.push_back( los_max ) ;
551   - bounds.push_back( m_halo_distance.value() );
  549 + std::vector<double> bounds;
  550 + bounds.push_back(los_min);
  551 + bounds.push_back(los_max);
  552 + bounds.push_back(m_halo_distance.value());
552 553  
553 554 // Compute value
554   - value = integral.romberg( bounds ) ;
  555 + value = integral.romberg(bounds);
555 556  
556 557 std::cout << "profile_value=" << value << std::endl;
557 558  
... ... @@ -559,6 +560,7 @@ double GModelSpatialRadialProfileDMEinasto::profile_value(const double&amp; theta) c
559 560 return value;
560 561 }
561 562  
  563 +
562 564 /***********************************************************************//**
563 565 * @brief Kernel for halo density profile squared
564 566 *
... ... @@ -593,53 +595,57 @@ double GModelSpatialRadialProfileDMEinasto::profile_value(const double&amp; theta) c
593 595 ***************************************************************************/
594 596 double GModelSpatialRadialProfileDMEinasto::halo_kernel_los::eval( const double &los )
595 597 {
  598 + // PLEASE ADD COMMENTS
  599 + double g = 0.0;
  600 + g = los * los;
  601 + g += m_halo_distance * m_halo_distance;
  602 + g -= 2.0 * los * m_halo_distance * std::cos(m_theta);
  603 + g = std::sqrt(g);
  604 + g /= m_scale_radius;
596 605  
597   - double g = 0.0 ;
598   - g = los * los ;
599   - g += m_halo_distance * m_halo_distance ;
600   - g -= 2.0 * los * m_halo_distance * std::cos(m_theta) ;
601   - g = sqrt(g) ;
602   - g /= m_scale_radius ;
  606 + double f = 0.0;
  607 + f = std::pow(g, m_alpha);
  608 + f -= 1.0;
  609 + f *= -2.0 / m_alpha;
  610 + f = std::exp(f);
603 611  
604   - double f = 0.0 ;
605   - f = pow( g , m_alpha ) ;
606   - f -= 1.0 ;
607   - f *= -2.0 / m_alpha ;
608   - f = std::exp( f ) ;
609   -
610   - // squared, for annihilating dm
611   - // would just be f if it was decaying dm
612   - f = f * f ;
  612 + // squared, for annihilating dm
  613 + // would just be f if it was decaying dm
  614 + f = f * f ;
613 615  
614   - //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;
615   -
616   - return f;
  616 + //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;
617 617  
  618 + // Return function value
  619 + return f;
618 620 }
619 621  
  622 +
620 623 /***********************************************************************//**
621 624 * @brief Update precomputation cache
622 625 *
623 626 * Computes the m_mass_radius calculation, determining the radius around
624   - * the halo that contains 99.99% of the mass. For an einasto halo profile,
  627 + * the halo that contains 99.99% of the mass. For an Einasto halo profile,
625 628 * this is just 10.0 * scale_radius .
626   - *
627 629 ***************************************************************************/
628 630 void GModelSpatialRadialProfileDMEinasto::update() const
629 631 {
630 632  
631   - // Update if scale radius has changed
632   - if ( m_last_scale_radius != scale_radius() ) {
  633 + // Update if scale radius has changed
  634 + if (m_last_scale_radius != scale_radius()) {
633 635  
634   - // Store last values
635   - m_last_scale_radius = scale_radius() ;
  636 + // Store last values
  637 + m_last_scale_radius = scale_radius();
636 638  
637   - // perform precomputations
638   - m_mass_radius = 10.0 * scale_radius() ;
  639 + // perform precomputations
  640 + m_mass_radius = 10.0 * scale_radius();
  641 +
  642 + }
639 643  
640   - }
  644 + // Return
  645 + return;
641 646 }
642 647  
  648 +
643 649 /***********************************************************************//**
644 650 * @brief Calculate Halo J Factor
645 651 *
... ... @@ -649,30 +655,31 @@ void GModelSpatialRadialProfileDMEinasto::update() const
649 655 * Calculates the halo's j-factor of the halo by integrating across the
650 656 * profile between two angles.
651 657 *
  658 + * @todo This should be done using a proper numerical integration
652 659 ***************************************************************************/
653   -double GModelSpatialRadialProfileDMEinasto::j_factor( const double& minangle, const double& maxangle, const int& npoints) const
  660 +double GModelSpatialRadialProfileDMEinasto::j_factor(const double& minangle,
  661 + const double& maxangle,
  662 + const int& npoints) const
654 663 {
  664 + // init variables
  665 + double minradian = minangle * gammalib::deg2rad;
  666 + double maxradian = maxangle * gammalib::deg2rad;
  667 + double dr = (maxradian - minradian) / npoints;
  668 + double r = 0.0;
  669 + double total = 0.0;
655 670  
656   - // init variables
657   - double minradian = minangle * gammalib::deg2rad ;
658   - double maxradian = maxangle * gammalib::deg2rad ;
659   - double dr = ( maxradian - minradian ) / npoints ;
660   - double r = 0.0 ;
661   - double total = 0.0 ;
662   -
663   - // loop over different radii in the profile
664   - for (int i = 0; i<npoints; ++i)
665   - {
  671 + // loop over different radii in the profile
  672 + for (int i = 0; i < npoints; ++i) {
666 673  
667   - // integration: Int[ profile(r) * r * dr ]
668   - r = minradian + ( i * dr ) ;
669   - total += profile_value( r ) * r * dr ;
  674 + // integration: Int[ profile(r) * r * dr ]
  675 + r = minradian + (i * dr);
  676 + total += profile_value(r) * r * dr;
670 677  
671   - }
  678 + }
672 679  
673   - // 2 * pi * Int[ profile(r) * r * dr
674   - total *= gammalib::twopi ;
  680 + // 2 * pi * Int[ profile(r) * r * dr
  681 + total *= gammalib::twopi;
675 682  
676   - return total ;
  683 + // Return J factor
  684 + return total;
677 685 }
678   -
... ...
src/model/GModelSpatialRadialProfileDMZhao.cpp
  View file @58af4ad
1 1 /***************************************************************************
2   - * GModelSpatialRadialProfileDMZhao.cpp - DMZhao radial profile class *
  2 + * GModelSpatialRadialProfileDMZhao.cpp - Zhao radial profile class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2016 by Nathan Kelley-Hoskins *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -28,19 +28,19 @@
28 28 #ifdef HAVE_CONFIG_H
29 29 #include <config.h>
30 30 #endif
  31 +//#include <iomanip>
31 32 #include "GException.hpp"
32 33 #include "GTools.hpp"
33 34 #include "GMath.hpp"
34 35 #include "GXmlElement.hpp"
35 36 #include "GModelSpatialRadialProfileDMZhao.hpp"
36 37 #include "GModelSpatialRegistry.hpp"
37   -#include <iomanip>
38 38  
39 39 /* __ Constants __________________________________________________________ */
40 40  
41 41 /* __ Globals ____________________________________________________________ */
42 42 const GModelSpatialRadialProfileDMZhao g_radial_disk_seed;
43   -const GModelSpatialRegistry g_radial_disk_registry(&g_radial_disk_seed);
  43 +const GModelSpatialRegistry g_radial_disk_registry(&g_radial_disk_seed);
44 44  
45 45 /* __ Method name definitions ____________________________________________ */
46 46 #define G_READ "GModelSpatialRadialProfileDMZhao::read(GXmlElement&)"
... ... @@ -65,7 +65,7 @@ const GModelSpatialRegistry g_radial_disk_registry(&amp;g_radial_disk_seed
65 65 * Constructs empty radial DMZhao profile
66 66 ***************************************************************************/
67 67 GModelSpatialRadialProfileDMZhao::GModelSpatialRadialProfileDMZhao(void) :
68   - GModelSpatialRadialProfile()
  68 + GModelSpatialRadialProfile()
69 69 {
70 70 // Initialise members
71 71 init_members();
... ... @@ -85,7 +85,7 @@ GModelSpatialRadialProfileDMZhao::GModelSpatialRadialProfileDMZhao(void) :
85 85 * expected structure of the XML element.
86 86 ***************************************************************************/
87 87 GModelSpatialRadialProfileDMZhao::GModelSpatialRadialProfileDMZhao(const GXmlElement& xml) :
88   - GModelSpatialRadialProfile()
  88 + GModelSpatialRadialProfile()
89 89 {
90 90 // Initialise members
91 91 init_members();
... ... @@ -106,7 +106,7 @@ GModelSpatialRadialProfileDMZhao::GModelSpatialRadialProfileDMZhao(const GXmlEle
106 106 * Copies radial DMZhao profile model from another radial profile model.
107 107 ***************************************************************************/
108 108 GModelSpatialRadialProfileDMZhao::GModelSpatialRadialProfileDMZhao(const GModelSpatialRadialProfileDMZhao& model) :
109   - GModelSpatialRadialProfile(model)
  109 + GModelSpatialRadialProfile(model)
110 110 {
111 111 // Initialise members
112 112 init_members();
... ... @@ -227,9 +227,10 @@ double GModelSpatialRadialProfileDMZhao::theta_min(void) const
227 227 update();
228 228  
229 229 // Return value
230   - return m_theta_min.value() ;
  230 + return m_theta_min.value();
231 231 }
232 232  
  233 +
233 234 /***********************************************************************//**
234 235 * @brief Return maximum model radius (in radians)
235 236 *
... ... @@ -241,27 +242,28 @@ double GModelSpatialRadialProfileDMZhao::theta_max(void) const
241 242 // update precomputation cache
242 243 update();
243 244  
244   - double theta = 0.0 ;
  245 + double theta = 0.0;
245 246  
246   - // if earth is within the significant radius, then theta_max must
247   - // contain the entire profile (180deg)
248   - if ( m_halo_distance.value() < m_mass_radius ) {
249   - theta = gammalib::pi ;
  247 + // If Earth is within the significant radius, then theta_max must
  248 + // contain the entire profile (180deg) ...
  249 + if (m_halo_distance.value() < m_mass_radius) {
  250 + theta = gammalib::pi;
  251 + }
250 252  
251   - // if the halo is far enough away (further than the significant radius)
252   - // then we just need to deal with the angles within the sphere of the
253   - // significant radius.
254   - } else {
255   - theta = std::atan( m_mass_radius / m_halo_distance.value() ) ;
  253 + // ... otherwise, if the halo is far enough away (further than the
  254 + // significant radius) then we just need to deal with the angles within
  255 + // the sphere of the significant radius.
  256 + else {
  257 + theta = std::atan(m_mass_radius / m_halo_distance.value());
256 258 }
257 259  
258   - // always chose the lesser of ( mass_radius theta, theta_max )
259   - if ( m_theta_max.value() * gammalib::deg2rad < theta ) {
260   - theta = m_theta_max.value() * gammalib::deg2rad ;
  260 + // Always chose the lesser of ( mass_radius theta, theta_max )
  261 + if (m_theta_max.value() * gammalib::deg2rad < theta) {
  262 + theta = m_theta_max.value() * gammalib::deg2rad;
261 263 }
262 264  
263 265 // Return value
264   - return theta ;
  266 + return theta;
265 267 }
266 268  
267 269  
... ... @@ -507,13 +509,13 @@ void GModelSpatialRadialProfileDMZhao::init_members(void)
507 509 m_theta_max.has_grad(false); // Radial components never have gradients
508 510  
509 511 // Set parameter pointer(s)
510   - m_pars.push_back(&m_scale_radius );
  512 + m_pars.push_back(&m_scale_radius);
511 513 m_pars.push_back(&m_halo_distance);
512   - m_pars.push_back(&m_alpha );
513   - m_pars.push_back(&m_beta );
514   - m_pars.push_back(&m_gamma );
515   - m_pars.push_back(&m_theta_min );
516   - m_pars.push_back(&m_theta_max );
  514 + m_pars.push_back(&m_alpha);
  515 + m_pars.push_back(&m_beta);
  516 + m_pars.push_back(&m_gamma);
  517 + m_pars.push_back(&m_theta_min);
  518 + m_pars.push_back(&m_theta_max);
517 519  
518 520 // Return
519 521 return;
... ... @@ -532,13 +534,13 @@ void GModelSpatialRadialProfileDMZhao::copy_members(const GModelSpatialRadialPro
532 534 // Copy members. We do not have to push back the members on the parameter
533 535 // stack as this should have been done by init_members() that was called
534 536 // before. Otherwise we would have sigma twice on the stack.
535   - m_scale_radius = model.m_scale_radius ;
536   - m_halo_distance = model.m_halo_distance ;
537   - m_alpha = model.m_alpha ;
538   - m_beta = model.m_beta ;
539   - m_gamma = model.m_gamma ;
540   - m_theta_max = model.m_theta_min ;
541   - m_theta_max = model.m_theta_max ;
  537 + m_scale_radius = model.m_scale_radius;
  538 + m_halo_distance = model.m_halo_distance;
  539 + m_alpha = model.m_alpha;
  540 + m_beta = model.m_beta;
  541 + m_gamma = model.m_gamma;
  542 + m_theta_max = model.m_theta_min;
  543 + m_theta_max = model.m_theta_max;
542 544  
543 545 // Return
544 546 return;
... ... @@ -568,32 +570,32 @@ double GModelSpatialRadialProfileDMZhao::profile_value(const double&amp; theta) cons
568 570 update();
569 571  
570 572 // initialize integral value
571   - double value = 0.0 ;
  573 + double value = 0.0;
572 574  
573 575 // Set up integration limits
574   - double los_min = m_halo_distance.value() - m_mass_radius ;
575   - double los_max = m_halo_distance.value() + m_mass_radius ;
  576 + double los_min = m_halo_distance.value() - m_mass_radius;
  577 + double los_max = m_halo_distance.value() + m_mass_radius;
576 578  
577 579 // Set up integral
578   - halo_kernel_los integrand( m_scale_radius.value(),
579   - m_halo_distance.value(),
580   - m_alpha.value(),
581   - m_beta.value(),
582   - m_gamma.value(),
583   - theta ) ;
584   - GIntegral integral(&integrand) ;
585   - integral.max_iter( 25 ) ;
  580 + halo_kernel_los integrand(m_scale_radius.value(),
  581 + m_halo_distance.value(),
  582 + m_alpha.value(),
  583 + m_beta.value(),
  584 + m_gamma.value(),
  585 + theta);
  586 + GIntegral integral(&integrand);
  587 + integral.max_iter(25);
586 588  
587 589 // Set up integration boundaries
588 590 // As there is usually an infinity at the halo center, this splits
589 591 // the integral at the m_halo_distance.
590   - std::vector<double> bounds ;
591   - bounds.push_back( los_min ) ;
592   - bounds.push_back( los_max ) ;
593   - bounds.push_back( m_halo_distance.value() );
  592 + std::vector<double> bounds;
  593 + bounds.push_back(los_min);
  594 + bounds.push_back(los_max);
  595 + bounds.push_back(m_halo_distance.value());
594 596  
595 597 // Compute value
596   - value = integral.romberg( bounds ) ;
  598 + value = integral.romberg(bounds);
597 599  
598 600 //std::cout << " theta=" << theta << " profile_value=" << value << std::endl;
599 601  
... ... @@ -601,6 +603,7 @@ double GModelSpatialRadialProfileDMZhao::profile_value(const double&amp; theta) cons
601 603 return value;
602 604 }
603 605  
  606 +
604 607 /***********************************************************************//**
605 608 * @brief Kernel for zhao halo density profile squared
606 609 *
... ... @@ -636,55 +639,64 @@ double GModelSpatialRadialProfileDMZhao::profile_value(const double&amp; theta) cons
636 639 ***************************************************************************/
637 640 double GModelSpatialRadialProfileDMZhao::halo_kernel_los::eval( const double &los )
638 641 {
  642 + // ADD COMMENTS
  643 + double g = 0.0;
  644 + g = los * los;
  645 + g += m_halo_distance * m_halo_distance;
  646 + g -= 2.0 * los * m_halo_distance * std::cos(m_theta);
  647 + g = std::sqrt(g) ;
  648 + g /= m_scale_radius ;
639 649  
640   - double g = 0.0 ;
641   - g = los * los ;
642   - g += m_halo_distance * m_halo_distance ;
643   - g -= 2.0 * los * m_halo_distance * std::cos(m_theta) ;
644   - g = sqrt(g) ;
645   - g /= m_scale_radius ;
646   -
647   - double f = 0.0 ;
648   - f = pow( g , m_alpha ) ;
649   - f += 1 ;
650   - f = pow( f, (m_beta-m_gamma)/m_alpha ) ;
651   - f *= pow( g, m_gamma ) ;
652   - f = 1 / f ;
653   -
654   - // squared, for annihilating dm
655   - // would just be f if it was decaying dm
656   - f = f * f ;
  650 + double f = 0.0 ;
  651 + f = std::pow(g, m_alpha);
  652 + f += 1 ;
  653 + f = std::pow(f, (m_beta-m_gamma)/m_alpha);
  654 + f *= std::pow(g, m_gamma);
  655 + f = 1 / f ;
  656 +
  657 + // squared, for annihilating dm
  658 + // would just be f if it was decaying dm
  659 + f = f * f;
657 660  
658   - //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 ;
659   -
660   - return f;
  661 + //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 ;
661 662  
  663 + // Return function value
  664 + return f;
662 665 }
663 666  
  667 +
664 668 /***********************************************************************//**
665 669 * @brief Update precomputation cache
666 670 *
667 671 * Computes the m_mass_radius calculation, determining the radius around
668   - * the halo that contains 99.99% of the mass. For a zhao halo profile,
  672 + * the halo that contains 99.99% of the mass. For a Zhao halo profile,
669 673 * this is just 10.0 * scale_radius .
670   - *
671 674 ***************************************************************************/
672 675 void GModelSpatialRadialProfileDMZhao::update() const
673 676 {
674 677  
675   - // Update if scale radius has changed
676   - if ( m_last_scale_radius != scale_radius() ) {
  678 + // Update if scale radius has changed
  679 + if (m_last_scale_radius != scale_radius()) {
677 680  
678   - // Store last values
679   - m_last_scale_radius = scale_radius() ;
  681 + // Store last values
  682 + m_last_scale_radius = scale_radius();
680 683  
681   - // perform precomputations
682   - m_mass_radius = 10.0 * scale_radius() ;
  684 + // perform precomputations
  685 + m_mass_radius = 10.0 * scale_radius();
683 686  
684   - }
  687 + }
  688 +
  689 + // Return
  690 + return;
685 691 }
686 692  
687   -double GModelSpatialRadialProfileDMZhao::prof_val( const double& theta )
  693 +
  694 +/***********************************************************************//**
  695 + * @brief Compute profile value
  696 + *
  697 + * @return Profile value
  698 + ***************************************************************************/
  699 +double GModelSpatialRadialProfileDMZhao::prof_val(const double& theta)
688 700 {
689   - return this->profile_value(theta) ;
  701 + return this->profile_value(theta);
690 702 }
... ...