Simon Bonnefoy / gammalib

Commit 4e0202a5ec3d218297c9eb1981c4f97d3e780423

Authored by Jürgen Knödlseder
1 parent eb61a120

Integrate code

Inline Side-by-side
Showing 13 changed files with 479 additions and 521 deletions
include/GModelSpatialRadialProfileDMBurkert.hpp
  View file @4e0202a
... ... @@ -74,10 +74,10 @@ public:
74 74 void scale_radius(const double& scale_radius);
75 75 double scale_density(void) const;
76 76 void scale_density(const double& scale_density);
77   - double halo_distance(void) const ;
  77 + double halo_distance(void) const;
78 78 void halo_distance(const double& halo_distance);
79   - double mass_density( const double& radius ) const;
80   - double jfactor( const double& angle ) const;
  79 + double mass_density(const double& radius ) const;
  80 + double jfactor(const double& angle) const;
81 81  
82 82 protected:
83 83 // Protected methods
... ... @@ -85,12 +85,12 @@ protected:
85 85 void copy_members(const GModelSpatialRadialProfileDMBurkert& model);
86 86 void free_members(void);
87 87 virtual double profile_value(const double& theta) const;
88   - void update(void) const ;
  88 + void update(void) const;
89 89  
90 90 // Integration kernel for line-of-sight integral
91 91 class halo_kernel_los : public GFunction {
92 92 public :
93   - halo_kernel_los(const double& scale_radius ,
  93 + halo_kernel_los(const double& scale_radius,
94 94 const double& halo_distance,
95 95 const double& theta,
96 96 const double& core_radius) :
... ...
include/GModelSpatialRadialProfileDMEinasto.hpp
  View file @4e0202a
... ... @@ -78,8 +78,8 @@ public:
78 78 void halo_distance(const double& halo_distance);
79 79 double alpha(void) const;
80 80 void alpha(const double& alpha);
81   - double mass_density( const double& radius ) const ;
82   - double jfactor( const double& angle ) const ;
  81 + double mass_density(const double& radius) const;
  82 + double jfactor(const double& angle) const;
83 83  
84 84 protected:
85 85 // Protected methods
... ... @@ -263,4 +263,3 @@ void GModelSpatialRadialProfileDMEinasto::alpha(const double& alpha)
263 263 }
264 264  
265 265 #endif /* GMODELSPATIALRADIALPROFILEDMEINASTO_HPP */
266   -
... ...
include/GModelSpatialRadialProfileDMZhao.hpp
  View file @4e0202a
... ... @@ -72,19 +72,18 @@ public:
72 72 // Other methods
73 73 double scale_radius(void) const;
74 74 void scale_radius(const double& scale_radius);
75   - //double prof_val(const double& theta);
76 75 double scale_density(void) const;
77 76 void scale_density(const double& scale_density);
78   - double halo_distance(void) const ;
  77 + double halo_distance(void) const;
79 78 void halo_distance(const double& halo_distance);
80   - double alpha(void) const ;
81   - void alpha(const double& alpha) ;
82   - double beta(void) const ;
83   - void beta(const double& beta) ;
84   - double gamma(void) const ;
85   - void gamma(const double& gamma) ;
  79 + double alpha(void) const;
  80 + void alpha(const double& alpha);
  81 + double beta(void) const;
  82 + void beta(const double& beta);
  83 + double gamma(void) const;
  84 + void gamma(const double& gamma);
86 85 double mass_density(const double& radius) const;
87   - double jfactor(const double& angle) const ;
  86 + double jfactor(const double& angle) const;
88 87  
89 88 protected:
90 89 // Protected methods
... ... @@ -134,10 +133,10 @@ protected:
134 133 GModelPar m_core_radius; //!< Core radius
135 134  
136 135 // Protected cached members
137   - mutable double m_last_scale_radius ;
138   - mutable double m_last_scale_density ;
139   - mutable double m_mass_radius ;
140   - mutable double m_scale_density_squared ;
  136 + mutable double m_last_scale_radius;
  137 + mutable double m_last_scale_density;
  138 + mutable double m_mass_radius;
  139 + mutable double m_scale_density_squared;
141 140 };
142 141  
143 142  
... ...
include/GammaLib.hpp
  View file @4e0202a
... ... @@ -205,6 +205,9 @@
205 205 #include "GModelSpatialRadialShell.hpp"
206 206 #include "GModelSpatialRadialProfile.hpp"
207 207 #include "GModelSpatialRadialProfileGauss.hpp"
  208 +#include "GModelSpatialRadialProfileDMBurkert.hpp"
  209 +#include "GModelSpatialRadialProfileDMEinasto.hpp"
  210 +#include "GModelSpatialRadialProfileDMZhao.hpp"
208 211 #include "GModelSpatialElliptical.hpp"
209 212 #include "GModelSpatialEllipticalDisk.hpp"
210 213 #include "GModelSpatialEllipticalGauss.hpp"
... ...
include/Makefile.am
  View file @4e0202a
... ... @@ -163,6 +163,9 @@ pkginclude_HEADERS = GammaLib.hpp \
163 163 GModelSpatialRadialShell.hpp \
164 164 GModelSpatialRadialProfile.hpp \
165 165 GModelSpatialRadialProfileGauss.hpp \
  166 + GModelSpatialRadialProfileDMBurkert.hpp \
  167 + GModelSpatialRadialProfileDMEinasto.hpp \
  168 + GModelSpatialRadialProfileDMZhao.hpp \
166 169 GModelSpatialElliptical.hpp \
167 170 GModelSpatialEllipticalDisk.hpp \
168 171 GModelSpatialEllipticalGauss.hpp \
... ...
pyext/GModelSpatialRadialProfileDMBurkert.i
  View file @4e0202a
... ... @@ -61,10 +61,10 @@ public:
61 61 void scale_radius(const double& scale_radius);
62 62 double scale_density(void) const;
63 63 void scale_density(const double& scale_density);
64   - double halo_distance(void) const ;
  64 + double halo_distance(void) const;
65 65 void halo_distance(const double& halo_distance);
66   - double mass_density( const double& radius ) const;
67   - double jfactor( const double& angle ) const;
  66 + double mass_density(const double& radius) const;
  67 + double jfactor(const double& angle) const;
68 68 };
69 69  
70 70  
... ...
pyext/GModelSpatialRadialProfileDMEinasto.i
  View file @4e0202a
... ... @@ -65,8 +65,8 @@ public:
65 65 void halo_distance(const double& halo_distance);
66 66 double alpha(void) const;
67 67 void alpha(const double& alpha);
68   - double mass_density( const double& radius ) const ;
69   - double jfactor( const double& angle ) const ;
  68 + double mass_density(const double& radius) const;
  69 + double jfactor(const double& angle) const;
70 70 };
71 71  
72 72  
... ...
pyext/GModelSpatialRadialProfileDMZhao.i
  View file @4e0202a
... ... @@ -61,17 +61,16 @@ public:
61 61 void scale_radius(const double& scale_radius);
62 62 double scale_density(void) const;
63 63 void scale_density(const double& scale_density);
64   - double halo_distance(void) const ;
  64 + double halo_distance(void) const;
65 65 void halo_distance(const double& halo_distance);
66   - double alpha(void) const ;
67   - void alpha(const double& alpha) ;
68   - double beta(void) const ;
69   - void beta(const double& beta) ;
70   - double gamma(void) const ;
71   - void gamma(const double& gamma) ;
72   - //double prof_val(const double& theta);
  66 + double alpha(void) const;
  67 + void alpha(const double& alpha);
  68 + double beta(void) const;
  69 + void beta(const double& beta);
  70 + double gamma(void) const;
  71 + void gamma(const double& gamma);
73 72 double mass_density(const double& radius) const;
74   - double jfactor(const double& angle) const ;
  73 + double jfactor(const double& angle) const;
75 74 };
76 75  
77 76  
... ...
pyext/gammalib/model.i
  View file @4e0202a
... ... @@ -154,10 +154,10 @@
154 154 %include "GModelSpatialRadialGauss.i"
155 155 %include "GModelSpatialRadialShell.i"
156 156 %include "GModelSpatialRadialProfile.i"
  157 +%include "GModelSpatialRadialProfileGauss.i"
157 158 %include "GModelSpatialRadialProfileDMBurkert.i"
158 159 %include "GModelSpatialRadialProfileDMEinasto.i"
159 160 %include "GModelSpatialRadialProfileDMZhao.i"
160   -%include "GModelSpatialRadialProfileGauss.i"
161 161 %include "GModelSpatialElliptical.i"
162 162 %include "GModelSpatialEllipticalDisk.i"
163 163 %include "GModelSpatialEllipticalGauss.i"
... ...
src/model/GModelSpatialRadialProfileDMBurkert.cpp
  View file @4e0202a
... ... @@ -34,7 +34,7 @@
34 34 #include "GXmlElement.hpp"
35 35 #include "GModelSpatialRadialProfileDMBurkert.hpp"
36 36 #include "GModelSpatialRegistry.hpp"
37   -#include <iomanip>
  37 +//#include <iomanip>
38 38  
39 39 /* __ Constants __________________________________________________________ */
40 40  
... ... @@ -43,8 +43,8 @@ const GModelSpatialRadialProfileDMBurkert g_radial_disk_seed;
43 43 const GModelSpatialRegistry g_radial_disk_registry(&g_radial_disk_seed);
44 44  
45 45 /* __ Method name definitions ____________________________________________ */
46   -#define G_READ "GModelSpatialRadialProfileDMBurkert::read(GXmlElement&)"
47   -#define G_WRITE "GModelSpatialRadialProfileDMBurkert::write(GXmlElement&)"
  46 +#define G_READ "GModelSpatialRadialProfileDMBurkert::read(GXmlElement&)"
  47 +#define G_WRITE "GModelSpatialRadialProfileDMBurkert::write(GXmlElement&)"
48 48  
49 49 /* __ Macros _____________________________________________________________ */
50 50  
... ... @@ -245,7 +245,7 @@ double GModelSpatialRadialProfileDMBurkert::theta_max(void) const
245 245  
246 246 // If Earth is within the significant radius, then we must integrate
247 247 // the entire profile ...
248   - if ( m_halo_distance.value() < m_mass_radius ) {
  248 + if (m_halo_distance.value() < m_mass_radius) {
249 249 theta = gammalib::pi;
250 250 }
251 251  
... ... @@ -253,12 +253,13 @@ double GModelSpatialRadialProfileDMBurkert::theta_max(void) const
253 253 // significant radius) then we just need to deal with the angles within
254 254 // the sphere of the significant radius
255 255 else {
256   - theta = std::atan( m_mass_radius / m_halo_distance.value() ) ;
  256 + theta = std::atan(m_mass_radius / m_halo_distance.value());
257 257 }
258 258  
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;
  259 + // Always chose the lesser of mass_radius theta and theta_max
  260 + double theta_max = m_theta_max.value() * gammalib::deg2rad;
  261 + if (theta > theta_max) {
  262 + theta = theta_max;
262 263 }
263 264  
264 265 // Return value
... ... @@ -272,49 +273,41 @@ double GModelSpatialRadialProfileDMBurkert::theta_max(void) const
272 273 * @param[in] xml XML element.
273 274 *
274 275 * Reads the DMBurkert radial profile model information from an XML element.
275   - * The XML element shall have either the format
  276 + * The XML element shall have the format
276 277 *
277 278 * <spatialModel type="DMBurkertProfile">
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"/>
281   - * <parameter name="Scale Radius" scale="1.0" value="21.5" min="1.0e-6" free="1"/>
282   - * <parameter name="Halo Distance" scale="1.0" value="7.94" min="1.0e-6" free="1"/>
283   - * <parameter name="Theta Max" scale="1.0" value="0.17" min="0.01" max="10" free="1"/>
  279 + * <parameter name="RA" scale=.. value=.. min=.. max=.. free=../>
  280 + * <parameter name="DEC" scale=.. value=.. min=.. max=.. free=../>
  281 + * <parameter name="ScaleRadius" scale=.. value=.. min=.. max=.. free=../>
  282 + * <parameter name="ScaleDensity" scale=.. value=.. min=.. max=.. free=../>
  283 + * <parameter name="HaloDistance" scale=.. value=.. min=.. max=.. free=../>
  284 + * <parameter name="ThetaMin" scale=.. value=.. min=.. max=.. free=../>
  285 + * <parameter name="ThetaMax" scale=.. value=.. min=.. max=.. free=../>
  286 + * <parameter name="CoreRadius" scale=.. value=.. min=.. max=.. free=../>
284 287 * </spatialModel>
285 288 *
286   - * or
287   - *
288   - * <spatialModel type="DMBurkertProfile">
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"/>
292   - * <parameter name="Scale Radius" scale="1.0" value="21.5" min="1.0e-6" free="1"/>
293   - * <parameter name="Halo Distance" scale="1.0" value="7.94" min="1.0e-6" free="1"/>
294   - * <parameter name="Theta Max" scale="1.0" value="0.17" min="0.01" max="10" free="1"/>
295   - * </spatialModel>
296 289 ***************************************************************************/
297 290 void GModelSpatialRadialProfileDMBurkert::read(const GXmlElement& xml)
298 291 {
299 292 // Read DMBurkert location
300 293 GModelSpatialRadial::read(xml);
301 294  
302   - const GXmlElement* par1 = gammalib::xml_get_par(G_READ, xml, "Scale Radius");
  295 + const GXmlElement* par1 = gammalib::xml_get_par(G_READ, xml, "ScaleRadius");
303 296 m_scale_radius.read(*par1);
304 297  
305   - const GXmlElement* par2 = gammalib::xml_get_par(G_READ, xml, "Scale Density");
  298 + const GXmlElement* par2 = gammalib::xml_get_par(G_READ, xml, "ScaleDensity");
306 299 m_scale_density.read(*par2);
307 300  
308   - const GXmlElement* par3 = gammalib::xml_get_par(G_READ, xml, "Halo Distance");
  301 + const GXmlElement* par3 = gammalib::xml_get_par(G_READ, xml, "HaloDistance");
309 302 m_halo_distance.read(*par3);
310 303  
311   - const GXmlElement* par4 = gammalib::xml_get_par(G_READ, xml, "Theta Min");
  304 + const GXmlElement* par4 = gammalib::xml_get_par(G_READ, xml, "ThetaMin");
312 305 m_theta_min.read(*par4);
313 306  
314   - const GXmlElement* par5 = gammalib::xml_get_par(G_READ, xml, "Theta Max");
  307 + const GXmlElement* par5 = gammalib::xml_get_par(G_READ, xml, "ThetaMax");
315 308 m_theta_max.read(*par5);
316 309  
317   - const GXmlElement* par6 = gammalib::xml_get_par(G_READ, xml, "Core Radius");
  310 + const GXmlElement* par6 = gammalib::xml_get_par(G_READ, xml, "CoreRadius");
318 311 m_core_radius.read(*par6);
319 312  
320 313 // Return
... ... @@ -331,12 +324,14 @@ void GModelSpatialRadialProfileDMBurkert::read(const GXmlElement&amp; xml)
331 324 * The XML element will have the format
332 325 *
333 326 * <spatialModel type="DMBurkertProfile">
334   - * <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
335   - * <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
336   - * <parameter name="Sigma" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
337   - * <parameter name="Scale Radius" scale="1.0" value="21.5" min="1.0e-6" free="1"/>
338   - * <parameter name="Halo Distance" scale="1.0" value="7.94" min="1.0e-6" free="1"/>
339   - * <parameter name="Theta Max" scale="1.0" value="0.17" min="0.01" max="10" free="1"/>
  327 + * <parameter name="RA" scale=.. value=.. min=.. max=.. free=../>
  328 + * <parameter name="DEC" scale=.. value=.. min=.. max=.. free=../>
  329 + * <parameter name="ScaleRadius" scale=.. value=.. min=.. max=.. free=../>
  330 + * <parameter name="ScaleDensity" scale=.. value=.. min=.. max=.. free=../>
  331 + * <parameter name="HaloDistance" scale=.. value=.. min=.. max=.. free=../>
  332 + * <parameter name="ThetaMin" scale=.. value=.. min=.. max=.. free=../>
  333 + * <parameter name="ThetaMax" scale=.. value=.. min=.. max=.. free=../>
  334 + * <parameter name="CoreRadius" scale=.. value=.. min=.. max=.. free=../>
340 335 * </spatialModel>
341 336 ***************************************************************************/
342 337 void GModelSpatialRadialProfileDMBurkert::write(GXmlElement& xml) const
... ... @@ -345,27 +340,27 @@ void GModelSpatialRadialProfileDMBurkert::write(GXmlElement&amp; xml) const
345 340 GModelSpatialRadial::write(xml);
346 341  
347 342 // Write Scale Radius parameter
348   - GXmlElement* par1 = gammalib::xml_need_par(G_WRITE, xml, "Scale Radius");
  343 + GXmlElement* par1 = gammalib::xml_need_par(G_WRITE, xml, "ScaleRadius");
349 344 m_scale_radius.write(*par1);
350 345  
351 346 // Write Scale Radius parameter
352   - GXmlElement* par2 = gammalib::xml_need_par(G_WRITE, xml, "Scale Density");
  347 + GXmlElement* par2 = gammalib::xml_need_par(G_WRITE, xml, "ScaleDensity");
353 348 m_scale_density.write(*par2);
354 349  
355 350 // Write Halo Distance parameter
356   - GXmlElement* par3 = gammalib::xml_need_par(G_WRITE, xml, "Halo Distance");
  351 + GXmlElement* par3 = gammalib::xml_need_par(G_WRITE, xml, "HaloDistance");
357 352 m_halo_distance.write(*par3);
358 353  
359 354 // Write Halo Distance parameter
360   - GXmlElement* par4 = gammalib::xml_need_par(G_WRITE, xml, "Theta Min");
  355 + GXmlElement* par4 = gammalib::xml_need_par(G_WRITE, xml, "ThetaMin");
361 356 m_theta_min.write(*par4);
362 357  
363 358 // Write Halo Distance parameter
364   - GXmlElement* par5 = gammalib::xml_need_par(G_WRITE, xml, "Theta Max");
  359 + GXmlElement* par5 = gammalib::xml_need_par(G_WRITE, xml, "ThetaMax");
365 360 m_theta_max.write(*par5);
366 361  
367 362 // Write Core Radius parameter
368   - GXmlElement* par6 = gammalib::xml_need_par(G_WRITE, xml, "Core Radius");
  363 + GXmlElement* par6 = gammalib::xml_need_par(G_WRITE, xml, "CoreRadius");
369 364 m_core_radius.write(*par6);
370 365  
371 366 // Return
... ... @@ -376,7 +371,7 @@ void GModelSpatialRadialProfileDMBurkert::write(GXmlElement&amp; xml) const
376 371 /***********************************************************************//**
377 372 * @brief Print information
378 373 *
379   - * @param[in] chatter Chattiness (defaults to NORMAL).
  374 + * @param[in] chatter Chattiness.
380 375 * @return String containing model information.
381 376 ***************************************************************************/
382 377 std::string GModelSpatialRadialProfileDMBurkert::print(const GChatter& chatter) const
... ... @@ -418,10 +413,10 @@ void GModelSpatialRadialProfileDMBurkert::init_members(void)
418 413  
419 414 // Initialise scale radius
420 415 m_scale_radius.clear();
421   - m_scale_radius.name("Scale Radius");
  416 + m_scale_radius.name("ScaleRadius");
422 417 m_scale_radius.unit("kpc");
423   - m_scale_radius.value(21.5); // default to GC scale radius
424   - m_scale_radius.min(1.0e-6); // arbitrarily chosen :/
  418 + m_scale_radius.value(21.5); // default to GC scale radius
  419 + m_scale_radius.min(1.0e-6); // arbitrarily chosen :/
425 420 m_scale_radius.free();
426 421 m_scale_radius.scale(1.0);
427 422 m_scale_radius.gradient(0.0);
... ... @@ -429,47 +424,47 @@ void GModelSpatialRadialProfileDMBurkert::init_members(void)
429 424  
430 425 // Initialise halo distance
431 426 m_halo_distance.clear();
432   - m_halo_distance.name("Halo Distance");
  427 + m_halo_distance.name("HaloDistance");
433 428 m_halo_distance.unit("kpc");
434   - m_halo_distance.value(7.94); // default to GC halo distance
435   - m_halo_distance.min(1.0e-6); // arbitrarily chosen
  429 + m_halo_distance.value(7.94); // default to GC halo distance
  430 + m_halo_distance.min(1.0e-6); // arbitrarily chosen
436 431 m_halo_distance.free();
437 432 m_halo_distance.scale(1.0);
438 433 m_halo_distance.gradient(0.0);
439   - m_halo_distance.has_grad(false); // Radial components never have gradients
  434 + m_halo_distance.has_grad(false); // Radial components never have gradients
440 435  
441 436 // Initialise theta max
442 437 m_theta_min.clear();
443   - m_theta_min.name("Theta Min");
  438 + m_theta_min.name("ThetaMin");
444 439 m_theta_min.unit("degrees");
445 440 m_theta_min.value(1.0e-6);
446   - m_theta_min.min(1.0e-10); // arbitrarily chosen
447   - m_theta_min.fix(); // should always be fixed!
  441 + m_theta_min.min(1.0e-10); // arbitrarily chosen
  442 + m_theta_min.fix(); // should always be fixed!
448 443 m_theta_min.scale(1.0);
449 444 m_theta_min.gradient(0.0);
450   - m_theta_min.has_grad(false); // Radial components never have gradients
  445 + m_theta_min.has_grad(false); // Radial components never have gradients
451 446  
452 447 // Initialise theta max
453 448 m_theta_max.clear();
454   - m_theta_max.name("Theta Max");
  449 + m_theta_max.name("ThetaMax");
455 450 m_theta_max.unit("degrees");
456   - m_theta_max.value(180.0); // can only go from halo center to opposite halo center
457   - m_theta_max.min(1.0e-6); // arbitrarily chosen
458   - m_theta_max.fix(); // should always be fixed!
  451 + m_theta_max.value(180.0); // can only go from halo center to opposite halo center
  452 + m_theta_max.min(1.0e-6); // arbitrarily chosen
  453 + m_theta_max.fix(); // should always be fixed!
459 454 m_theta_max.scale(1.0);
460 455 m_theta_max.gradient(0.0);
461   - m_theta_max.has_grad(false); // Radial components never have gradients
  456 + m_theta_max.has_grad(false); // Radial components never have gradients
462 457  
463 458 // Initialise core radius
464 459 m_core_radius.clear();
465   - m_core_radius.name("Core Radius");
  460 + m_core_radius.name("CoreRadius");
466 461 m_core_radius.unit("kpc");
467   - m_core_radius.value(0.5); // example: galactic center core
  462 + m_core_radius.value(0.5); // example: galactic center core
468 463 m_core_radius.min(0.0);
469 464 m_core_radius.fix();
470 465 m_core_radius.scale(1.0);
471 466 m_core_radius.gradient(0.0);
472   - m_core_radius.has_grad(false);
  467 + m_core_radius.has_grad(false); // Radial components never have gradients
473 468  
474 469 // Set parameter pointer(s)
475 470 m_pars.push_back(&m_scale_radius);
... ... @@ -512,7 +507,7 @@ void GModelSpatialRadialProfileDMBurkert::copy_members(const GModelSpatialRadial
512 507 m_last_scale_radius = model.m_last_scale_radius;
513 508 m_last_scale_density = model.m_last_scale_density;
514 509 m_mass_radius = model.m_mass_radius;
515   - m_scale_density_squared = model.m_scale_density_squared ;
  510 + m_scale_density_squared = model.m_scale_density_squared;
516 511  
517 512 // Return
518 513 return;
... ... @@ -537,12 +532,10 @@ void GModelSpatialRadialProfileDMBurkert::free_members(void)
537 532 ***************************************************************************/
538 533 double GModelSpatialRadialProfileDMBurkert::profile_value(const double& theta) const
539 534 {
540   -
541   -
542   - // update precompuation cache
  535 + // Update precompuation cache
543 536 update();
544 537  
545   - // initialize integral value
  538 + // Initialize integral value
546 539 double value = 0.0;
547 540  
548 541 // Set up integration limits
... ... @@ -550,8 +543,8 @@ double GModelSpatialRadialProfileDMBurkert::profile_value(const double&amp; theta) c
550 543 double los_max = m_halo_distance.value() + m_mass_radius;
551 544  
552 545 // Handle case where observer is within halo mass radius
553   - if ( los_min < 0.0 ) {
554   - los_min = 0.0 ;
  546 + if (los_min < 0.0) {
  547 + los_min = 0.0;
555 548 }
556 549  
557 550 // Set up integral
... ... @@ -562,9 +555,7 @@ double GModelSpatialRadialProfileDMBurkert::profile_value(const double&amp; theta) c
562 555 GIntegral integral(&integrand);
563 556  
564 557 // Compute value
565   - value = integral.romberg(los_min, los_max);
566   -
567   - value *= m_scale_density_squared ;
  558 + value = integral.romberg(los_min, los_max) * m_scale_density_squared;
568 559  
569 560 // Return value
570 561 return value;
... ... @@ -575,12 +566,14 @@ double GModelSpatialRadialProfileDMBurkert::profile_value(const double&amp; theta) c
575 566 * @brief Kernel for halo density profile squared
576 567 *
577 568 * @param[in] distance from observer to point in space (meters)
  569 + * @return Density profile value
578 570 *
579 571 * Computes the value of an einasto halo density profile squared,
580 572 * at distance l from observer, at angle \f[\theta\f] from the halo center:
581 573 *
582 574 * \f[
583   - * f(\theta, l) = \frac{r_{scale}^3}{ \left( r + r_{scale} \right) \left( r^2 + r_{scale}^2 \right)}
  575 + * f(\theta, l) = \frac{r_{scale}^3}{ \left( r + r_{scale} \right)
  576 + * \left( r^2 + r_{scale}^2 \right)}
584 577 * \f]
585 578 *
586 579 * where
... ... @@ -596,40 +589,34 @@ double GModelSpatialRadialProfileDMBurkert::profile_value(const double&amp; theta) c
596 589 * http://iopscience.iop.org/article/10.1086/309560/pdf
597 590 * Equation 2
598 591 *
599   - * @return unit
600   - *
601 592 ***************************************************************************/
602 593 double GModelSpatialRadialProfileDMBurkert::halo_kernel_los::eval(const double &los)
603 594 {
604   - // calculate the scale distance r, the ( distance from integration point
  595 + // Calculate the scale distance r, the ( distance from integration point
605 596 // to the halo center ) divided by ( the halo scale radius )
606 597  
607 598 // first calculate the distance of the integration point from the halo
608 599 // center via the law of cosines
609   - double r = 0.0;
610   - r = los * los;
611   - r += m_halo_distance * m_halo_distance;
612   - r -= 2.0 * los * m_halo_distance * std::cos(m_theta);
613   - r = std::sqrt(r);
  600 + double r = los * los;
  601 + r += m_halo_distance * m_halo_distance;
  602 + r -= 2.0 * los * m_halo_distance * std::cos(m_theta);
  603 + r = std::sqrt(r);
614 604  
615   - // if we have a core radius specified, all halo values inside this core
  605 + // If we have a core radius specified, all halo values inside this core
616 606 // radius should be the same as at the core radius itself.
617   - if ( r < m_core_radius )
618   - {
619   - r = m_core_radius ;
  607 + if (r < m_core_radius) {
  608 + r = m_core_radius;
620 609 }
621 610  
622   - // now that we've found out how far we are from the halo center,
  611 + // Now that we've found out how far we are from the halo center,
623 612 // feed this radius to the burkert formula
624 613 double bot = (m_scale_radius + r);
625   - bot *= (m_scale_radius*m_scale_radius) + (r*r);
626   -
627   - double f = m_scale_radius * m_scale_radius * m_scale_radius;
628   - f /= bot;
  614 + bot *= (m_scale_radius*m_scale_radius) + (r*r);
  615 + double f = m_scale_radius * m_scale_radius * m_scale_radius;
  616 + f /= bot;
629 617  
630   - // squared, for annihilating dm
631   - // would just be f if it was decaying dm
632   - f = f * f;
  618 + // Squared, for annihilating DM, would just be f if it was decaying dm
  619 + f = f * f;
633 620  
634 621 // Return function value
635 622 return f;
... ... @@ -646,7 +633,8 @@ double GModelSpatialRadialProfileDMBurkert::halo_kernel_los::eval(const double &amp;
646 633 void GModelSpatialRadialProfileDMBurkert::update() const
647 634 {
648 635 // Update if scale radius has changed
649   - if (m_last_scale_radius != scale_radius() || m_last_scale_density != scale_density() ) {
  636 + if (m_last_scale_radius != scale_radius() ||
  637 + m_last_scale_density != scale_density() ) {
650 638  
651 639 // Store last values
652 640 m_last_scale_radius = scale_radius();
... ... @@ -656,75 +644,76 @@ void GModelSpatialRadialProfileDMBurkert::update() const
656 644 // set the mass radius to 80*scale_radius, meaning
657 645 // 99.99% of the mass is contained within the mass radius,
658 646 // and integration only needs to worry about whats inside this radius.
659   - m_mass_radius = 80.0 * scale_radius();
  647 + m_mass_radius = 80.0 * scale_radius();
660 648 m_scale_density_squared = scale_density() * scale_density() ;
661   -
  649 +
662 650 }
663 651  
664 652 // Return
665 653 return;
666 654 }
667 655  
  656 +
668 657 /***********************************************************************//**
669 658 * @brief Calculate Halo Mass Density
670 659 *
671 660 * @param[in] distance from halo center (kpc)
  661 + * @return Halo mass density
672 662 *
673 663 * Calculates the halo's mass density at a given radial distance from the halo
674 664 * center.
675   - *
676 665 ***************************************************************************/
677 666 double GModelSpatialRadialProfileDMBurkert::mass_density(const double& radius) const
678 667 {
679   - double density = 0.0 ;
680   -
681   - // initialize halo kernel with stored values
682   - // we use theta=0 to effectivly sample the halo radially
683   - halo_kernel_los halo_shape( m_scale_radius.value(),
684   - m_halo_distance.value(),
685   - 0.0,
686   - m_core_radius.value());
  668 + // Initialize halo kernel with stored values. We use theta=0 to effectivly
  669 + // sample the halo radially
  670 + halo_kernel_los halo_shape(m_scale_radius.value(),
  671 + m_halo_distance.value(),
  672 + 0.0,
  673 + m_core_radius.value());
687 674  
688 675 // eval produces a unitless density^2, so we must take its square root
689   - density = std::sqrt( halo_shape.eval( m_halo_distance.value() + radius ) ) ;
  676 + double density = std::sqrt(halo_shape.eval(m_halo_distance.value() + radius));
690 677  
691   - // multiply in the missing scale density
692   - density *= m_scale_density.value() ;
  678 + // Multiply in the missing scale density
  679 + density *= m_scale_density.value();
693 680  
694   - return density ;
  681 + // Return halo mass density
  682 + return density;
695 683 }
696 684  
  685 +
697 686 /***********************************************************************//**
698   - * @brief Calculate J Factor
699   - *
700   - * @param[in] angle from halo center (radians)
  687 + * @brief Calculate J-factor
701 688 *
702   - * Calculates the halo's J-Factor at an angle from the halo center.
  689 + * @param[in] angle Angle from halo center (radians).
  690 + * @return J-factor.
703 691 *
  692 + * Calculates the halo's J-factor at an angle from the halo center.
704 693 ***************************************************************************/
705   -double GModelSpatialRadialProfileDMBurkert::jfactor( const double& angle ) const
  694 +double GModelSpatialRadialProfileDMBurkert::jfactor(const double& angle) const
706 695 {
707   - // Integration settings
708   - double minradian = 0.0 ;
709   - int npoints = 200 ;
  696 + // Integration settings
  697 + double minradian = 0.0;
  698 + int npoints = 200;
710 699  
711   - // initialize other variables
712   - double jfactor = 0.0 ;
713   - double dr = (angle - minradian) / npoints ;
714   - double r = 0.0;
  700 + // Initialize other variables
  701 + double jfactor = 0.0 ;
  702 + double dr = (angle - minradian) / npoints ;
  703 + double r = 0.0;
715 704  
716   - // loop over different radii in the profile
717   - for (int i = 0; i < npoints; ++i) {
  705 + // Loop over different radii in the profile
  706 + for (int i = 0; i < npoints; ++i) {
718 707  
719   - // integration: Int[ profile(r) * r * dr ]
720   - r = minradian + (i * dr);
721   - jfactor += profile_value(r) * r * dr;
  708 + // integration: Int[ profile(r) * r * dr ]
  709 + r = minradian + (i * dr);
  710 + jfactor += profile_value(r) * r * dr;
722 711  
723   - }
724   -
725   - // J-Factor = 2 * pi * Int[ profile(r) * r * dr , {r,minradian,angle} ]
726   - jfactor *= gammalib::twopi;
  712 + }
727 713  
728   - return jfactor ;
729   -}
  714 + // J-Factor = 2 * pi * Int[ profile(r) * r * dr , {r,minradian,angle} ]
  715 + jfactor *= gammalib::twopi;
730 716  
  717 + // Return J-factor
  718 + return jfactor;
  719 +}
... ...
src/model/GModelSpatialRadialProfileDMEinasto.cpp
  View file @4e0202a
... ... @@ -34,7 +34,7 @@
34 34 #include "GXmlElement.hpp"
35 35 #include "GModelSpatialRadialProfileDMEinasto.hpp"
36 36 #include "GModelSpatialRegistry.hpp"
37   -#include <iomanip>
  37 +//#include <iomanip>
38 38  
39 39 /* __ Constants __________________________________________________________ */
40 40  
... ... @@ -43,8 +43,8 @@ const GModelSpatialRadialProfileDMEinasto g_radial_disk_seed;
43 43 const GModelSpatialRegistry g_radial_disk_registry(&g_radial_disk_seed);
44 44  
45 45 /* __ Method name definitions ____________________________________________ */
46   -#define G_READ "GModelSpatialRadialProfileDMEinasto::read(GXmlElement&)"
47   -#define G_WRITE "GModelSpatialRadialProfileDMEinasto::write(GXmlElement&)"
  46 +#define G_READ "GModelSpatialRadialProfileDMEinasto::read(GXmlElement&)"
  47 +#define G_WRITE "GModelSpatialRadialProfileDMEinasto::write(GXmlElement&)"
48 48  
49 49 /* __ Macros _____________________________________________________________ */
50 50  
... ... @@ -215,6 +215,7 @@ GModelSpatialRadialProfileDMEinasto* GModelSpatialRadialProfileDMEinasto::clone(
215 215 return new GModelSpatialRadialProfileDMEinasto(*this);
216 216 }
217 217  
  218 +
218 219 /***********************************************************************//**
219 220 * @brief Return minimum model radius (in radians)
220 221 *
... ... @@ -230,6 +231,7 @@ double GModelSpatialRadialProfileDMEinasto::theta_min(void) const
230 231 return m_theta_min.value();
231 232 }
232 233  
  234 +
233 235 /***********************************************************************//**
234 236 * @brief Return maximum model radius (in radians)
235 237 *
... ... @@ -238,9 +240,10 @@ double GModelSpatialRadialProfileDMEinasto::theta_min(void) const
238 240 double GModelSpatialRadialProfileDMEinasto::theta_max(void) const
239 241 {
240 242  
241   - // update precomputation cache
  243 + // Update precomputation cache
242 244 update();
243   -
  245 +
  246 + // Initialise maximum theta angle
244 247 double theta = 0.0;
245 248  
246 249 // If Earth is within the significant radius, then theta_max must
... ... @@ -256,12 +259,13 @@ double GModelSpatialRadialProfileDMEinasto::theta_max(void) const
256 259 theta = std::atan(m_mass_radius / m_halo_distance.value());
257 260 }
258 261  
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;
  262 + // Always chose the lesser of mass_radius theta and theta_max
  263 + double theta_max = m_theta_max.value() * gammalib::deg2rad;
  264 + if (theta > theta_max) {
  265 + theta = theta_max;
262 266 }
263 267  
264   - // Return value
  268 + // Return maximum theta value
265 269 return theta;
266 270 }
267 271  
... ... @@ -272,30 +276,18 @@ double GModelSpatialRadialProfileDMEinasto::theta_max(void) const
272 276 * @param[in] xml XML element.
273 277 *
274 278 * Reads the DMEinasto radial profile model information from an XML element.
275   - * The XML element shall have either the format
  279 + * The XML element shall have the format
276 280 *
277 281 * <spatialModel type="DMEinastoProfile">
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"/>
281   - * <parameter name="Scale Radius" scale="1.0" value="21.5" min="1.0e-6" free="1"/>
282   - * <parameter name="Scale Density" scale="1.0" value="0.2" free="1"/>
283   - * <parameter name="Halo Distance" scale="1.0" value="7.94" min="1.0e-6" free="1"/>
284   - * <parameter name="Alpha" scale="1.0" value="0.17" min="0.01" max="10" free="1"/>
285   - * <parameter name='Core Radius" scale="1.0" value="10.0" min="1.0e-6" free="1"/>
286   - * </spatialModel>
287   - *
288   - * or
289   - *
290   - * <spatialModel type="DMEinastoProfile">
291   - * <parameter name="GLON" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
292   - * <parameter name="GLAT" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
293   - * <parameter name="Sigma" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
294   - * <parameter name="Scale Radius" scale="1.0" value="21.5" min="1.0e-6" free="1"/>
295   - * <parameter name="Scale Density" scale="1.0" value="0.2" free="1"/>
296   - * <parameter name="Halo Distance" scale="1.0" value="7.94" min="1.0e-6" free="1"/>
297   - * <parameter name="Alpha" scale="1.0" value="0.17" min="0.01" max="10" free="1"/>
298   - * <parameter name='Core Radius" scale="1.0" value="10.0" min="1.0e-6" free="1"/>
  282 + * <parameter name="RA" scale=.. value=.. min=.. max=.. free=../>
  283 + * <parameter name="DEC" scale=.. value=.. min=.. max=.. free=../>
  284 + * <parameter name="ScaleRadius" scale=.. value=.. min=.. max=.. free=../>
  285 + * <parameter name="ScaleDensity" scale=.. value=.. min=.. max=.. free=../>
  286 + * <parameter name="HaloDistance" scale=.. value=.. min=.. max=.. free=../>
  287 + * <parameter name="Alpha" scale=.. value=.. min=.. max=.. free=../>
  288 + * <parameter name="ThetaMin" scale=.. value=.. min=.. max=.. free=../>
  289 + * <parameter name="ThetaMax" scale=.. value=.. min=.. max=.. free=../>
  290 + * <parameter name="CoreRadius" scale=.. value=.. min=.. max=.. free=../>
299 291 * </spatialModel>
300 292 ***************************************************************************/
301 293 void GModelSpatialRadialProfileDMEinasto::read(const GXmlElement& xml)
... ... @@ -303,25 +295,32 @@ void GModelSpatialRadialProfileDMEinasto::read(const GXmlElement&amp; xml)
303 295 // Read DMEinasto location
304 296 GModelSpatialRadial::read(xml);
305 297  
306   - const GXmlElement* par1 = gammalib::xml_get_par(G_READ, xml, "Scale Radius");
  298 + // Read ScaleRadius parameter
  299 + const GXmlElement* par1 = gammalib::xml_get_par(G_READ, xml, "ScaleRadius");
307 300 m_scale_radius.read(*par1);
308 301  
309   - const GXmlElement* par2 = gammalib::xml_get_par(G_READ, xml, "Scale Density");
  302 + // Read ScaleDensity parameter
  303 + const GXmlElement* par2 = gammalib::xml_get_par(G_READ, xml, "ScaleDensity");
310 304 m_scale_density.read(*par2);
311 305  
312   - const GXmlElement* par3 = gammalib::xml_get_par(G_READ, xml, "Halo Distance");
  306 + // Read HaloDistance parameter
  307 + const GXmlElement* par3 = gammalib::xml_get_par(G_READ, xml, "HaloDistance");
313 308 m_halo_distance.read(*par3);
314 309  
  310 + // Read Alpha parameter
315 311 const GXmlElement* par4 = gammalib::xml_get_par(G_READ, xml, "Alpha");
316 312 m_alpha.read(*par4);
317 313  
318   - const GXmlElement* par5 = gammalib::xml_get_par(G_READ, xml, "Theta Min");
  314 + // Read ThetaMin parameter
  315 + const GXmlElement* par5 = gammalib::xml_get_par(G_READ, xml, "ThetaMin");
319 316 m_theta_min.read(*par5);
320 317  
321   - const GXmlElement* par6 = gammalib::xml_get_par(G_READ, xml, "Theta Max");
  318 + // Read ThetaMax parameter
  319 + const GXmlElement* par6 = gammalib::xml_get_par(G_READ, xml, "ThetaMax");
322 320 m_theta_max.read(*par6);
323 321  
324   - const GXmlElement* par7 = gammalib::xml_get_par(G_READ, xml, "Core Radius");
  322 + // Read CoreRadius parameter
  323 + const GXmlElement* par7 = gammalib::xml_get_par(G_READ, xml, "CoreRadius");
325 324 m_core_radius.read(*par7);
326 325  
327 326 // Return
... ... @@ -338,9 +337,15 @@ void GModelSpatialRadialProfileDMEinasto::read(const GXmlElement&amp; xml)
338 337 * The XML element will have the format
339 338 *
340 339 * <spatialModel type="DMEinastoProfile">
341   - * <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
342   - * <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
343   - * <parameter name="Sigma" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
  340 + * <parameter name="RA" scale=.. value=.. min=.. max=.. free=../>
  341 + * <parameter name="DEC" scale=.. value=.. min=.. max=.. free=../>
  342 + * <parameter name="ScaleRadius" scale=.. value=.. min=.. max=.. free=../>
  343 + * <parameter name="ScaleDensity" scale=.. value=.. min=.. max=.. free=../>
  344 + * <parameter name="HaloDistance" scale=.. value=.. min=.. max=.. free=../>
  345 + * <parameter name="Alpha" scale=.. value=.. min=.. max=.. free=../>
  346 + * <parameter name="ThetaMin" scale=.. value=.. min=.. max=.. free=../>
  347 + * <parameter name="ThetaMax" scale=.. value=.. min=.. max=.. free=../>
  348 + * <parameter name="CoreRadius" scale=.. value=.. min=.. max=.. free=../>
344 349 * </spatialModel>
345 350 ***************************************************************************/
346 351 void GModelSpatialRadialProfileDMEinasto::write(GXmlElement& xml) const
... ... @@ -348,32 +353,32 @@ void GModelSpatialRadialProfileDMEinasto::write(GXmlElement&amp; xml) const
348 353 // Write DMEinasto location
349 354 GModelSpatialRadial::write(xml);
350 355  
351   - // Write Scale Radius parameter
352   - GXmlElement* par1 = gammalib::xml_need_par(G_WRITE, xml, "Scale Radius");
  356 + // Write ScaleRadius parameter
  357 + GXmlElement* par1 = gammalib::xml_need_par(G_WRITE, xml, "ScaleRadius");
353 358 m_scale_radius.write(*par1);
354 359  
355   - // Write Scale Density parameter
356   - GXmlElement* par2 = gammalib::xml_need_par(G_WRITE, xml, "Scale Density");
  360 + // Write ScaleDensity parameter
  361 + GXmlElement* par2 = gammalib::xml_need_par(G_WRITE, xml, "ScaleDensity");
357 362 m_scale_density.write(*par2);
358 363  
359   - // Write Halo Distance parameter
360   - GXmlElement* par3 = gammalib::xml_need_par(G_WRITE, xml, "Halo Distance");
  364 + // Write HaloDistance parameter
  365 + GXmlElement* par3 = gammalib::xml_need_par(G_WRITE, xml, "HaloDistance");
361 366 m_halo_distance.write(*par3);
362 367  
363 368 // Write Alpha parameter
364 369 GXmlElement* par4 = gammalib::xml_need_par(G_WRITE, xml, "Alpha");
365 370 m_alpha.write(*par4);
366 371  
367   - // Write Theta Min parameter
368   - GXmlElement* par5 = gammalib::xml_need_par(G_WRITE, xml, "Theta Min");
  372 + // Write ThetaMin parameter
  373 + GXmlElement* par5 = gammalib::xml_need_par(G_WRITE, xml, "ThetaMin");
369 374 m_theta_min.write(*par5);
370 375  
371   - // Write Theta Max parameter
372   - GXmlElement* par6 = gammalib::xml_need_par(G_WRITE, xml, "Theta Max");
  376 + // Write ThetaMax parameter
  377 + GXmlElement* par6 = gammalib::xml_need_par(G_WRITE, xml, "ThetaMax");
373 378 m_theta_max.write(*par6);
374 379  
375   - // Write Core Radius parameter
376   - GXmlElement* par7 = gammalib::xml_need_par(G_WRITE, xml, "Core Radius");
  380 + // Write CoreRadius parameter
  381 + GXmlElement* par7 = gammalib::xml_need_par(G_WRITE, xml, "CoreRadius");
377 382 m_core_radius.write(*par7);
378 383  
379 384 // Return
... ... @@ -384,7 +389,7 @@ void GModelSpatialRadialProfileDMEinasto::write(GXmlElement&amp; xml) const
384 389 /***********************************************************************//**
385 390 * @brief Print information
386 391 *
387   - * @param[in] chatter Chattiness (defaults to NORMAL).
  392 + * @param[in] chatter Chattiness.
388 393 * @return String containing model information.
389 394 ***************************************************************************/
390 395 std::string GModelSpatialRadialProfileDMEinasto::print(const GChatter& chatter) const
... ... @@ -426,20 +431,20 @@ void GModelSpatialRadialProfileDMEinasto::init_members(void)
426 431  
427 432 // Initialise scale radius
428 433 m_scale_radius.clear();
429   - m_scale_radius.name("Scale Radius");
  434 + m_scale_radius.name("ScaleRadius");
430 435 m_scale_radius.unit("kpc");
431   - m_scale_radius.value(21.5); // default to GC scale radius
432   - m_scale_radius.min(1.0e-6); // arbitrarily chosen :/
  436 + m_scale_radius.value(21.5); // default to GC scale radius
  437 + m_scale_radius.min(1.0e-6); // arbitrarily chosen :/
433 438 m_scale_radius.free();
434 439 m_scale_radius.scale(1.0);
435 440 m_scale_radius.gradient(0.0);
436   - m_scale_radius.has_grad(false); // Radial components never have gradients
  441 + m_scale_radius.has_grad(false); // Radial components never have gradients
437 442  
438 443 // Initialise scale density
439 444 m_scale_density.clear();
440   - m_scale_density.name("Scale Density");
  445 + m_scale_density.name("ScaleDensity");
441 446 m_scale_density.unit("GeV/cm^3");
442   - m_scale_density.value(0.2); // GeV/cm3, default to GC scale density
  447 + m_scale_density.value(0.2); // GeV/cm3, default to GC scale density
443 448 m_scale_density.min(1.0e-6);
444 449 m_scale_density.fix();
445 450 m_scale_density.scale(1.0);
... ... @@ -448,54 +453,54 @@ void GModelSpatialRadialProfileDMEinasto::init_members(void)
448 453  
449 454 // Initialise halo distance
450 455 m_halo_distance.clear();
451   - m_halo_distance.name("Halo Distance");
  456 + m_halo_distance.name("HaloDistance");
452 457 m_halo_distance.unit("kpc");
453   - m_halo_distance.value(7.94); // default to GC halo distance
454   - m_halo_distance.min(1.0e-6); // arbitrarily chosen
  458 + m_halo_distance.value(7.94); // default to GC halo distance
  459 + m_halo_distance.min(1.0e-6); // arbitrarily chosen
455 460 m_halo_distance.free();
456 461 m_halo_distance.scale(1.0);
457 462 m_halo_distance.gradient(0.0);
458   - m_halo_distance.has_grad(false); // Radial components never have gradients
  463 + m_halo_distance.has_grad(false); // Radial components never have gradients
459 464  
460 465 // Initialise alpha
461 466 m_alpha.clear();
462 467 m_alpha.name("Alpha");
463 468 m_alpha.unit("unitless");
464   - m_alpha.value(0.17); // default to GC alpha
465   - m_alpha.min(0.01); // arbitrarily chosen
466   - m_alpha.max(10.0); // arbitrarily chosen
  469 + m_alpha.value(0.17); // default to GC alpha
  470 + m_alpha.min(0.01); // arbitrarily chosen
  471 + m_alpha.max(10.0); // arbitrarily chosen
467 472 m_alpha.free();
468 473 m_alpha.scale(1.0);
469 474 m_alpha.gradient(0.0);
470   - m_alpha.has_grad(false); // Radial components never have gradients
  475 + m_alpha.has_grad(false); // Radial components never have gradients
471 476  
472 477 // Initialise theta max
473 478 m_theta_min.clear();
474   - m_theta_min.name("Theta Min");
  479 + m_theta_min.name("ThetaMin");
475 480 m_theta_min.unit("degrees");
476   - m_theta_min.value(180.0); // can only go from halo center to opposite halo center
477   - m_theta_min.min(1.0e-6); // arbitrarily chosen
478   - m_theta_min.fix(); // should always be fixed!
  481 + m_theta_min.value(180.0); // can only go from halo center to opposite halo center
  482 + m_theta_min.min(1.0e-6); // arbitrarily chosen
  483 + m_theta_min.fix(); // should always be fixed!
479 484 m_theta_min.scale(1.0);
480 485 m_theta_min.gradient(0.0);
481   - m_theta_min.has_grad(false); // Radial components never have gradients
  486 + m_theta_min.has_grad(false); // Radial components never have gradients
482 487  
483 488 // Initialise theta max
484 489 m_theta_max.clear();
485   - m_theta_max.name("Theta Max");
  490 + m_theta_max.name("ThetaMax");
486 491 m_theta_max.unit("degrees");
487   - m_theta_max.value(1.0e-6); // can only go from halo center to opposite halo center
488   - m_theta_max.min(1.0e-10); // arbitrarily chosen
489   - m_theta_max.fix(); // should always be fixed!
  492 + m_theta_max.value(1.0e-6); // can only go from halo center to opposite halo center
  493 + m_theta_max.min(1.0e-10); // arbitrarily chosen
  494 + m_theta_max.fix(); // should always be fixed!
490 495 m_theta_max.scale(1.0);
491 496 m_theta_max.gradient(0.0);
492   - m_theta_max.has_grad(false); // Radial components never have gradients
  497 + m_theta_max.has_grad(false); // Radial components never have gradients
493 498  
494 499 // Initialise core radius
495 500 m_core_radius.clear();
496   - m_core_radius.name("Core Radius");
  501 + m_core_radius.name("CoreRadius");
497 502 m_core_radius.unit("kpc");
498   - m_core_radius.value(0.5); // example: galactic center core
  503 + m_core_radius.value(0.5); // example: galactic center core
499 504 m_core_radius.min(0.0);
500 505 m_core_radius.fix();
501 506 m_core_radius.scale(1.0);
... ... @@ -513,10 +518,10 @@ void GModelSpatialRadialProfileDMEinasto::init_members(void)
513 518  
514 519 // Initialize precomputation cache. Note that zero values flag
515 520 // uninitialised, as a zero radius is not meaningful
516   - m_last_scale_radius = 0.0 ;
517   - m_last_scale_density = 0.0 ;
518   - m_mass_radius = 0.0 ;
519   - m_scale_density_squared = 0.0 ;
  521 + m_last_scale_radius = 0.0;
  522 + m_last_scale_density = 0.0;
  523 + m_mass_radius = 0.0;
  524 + m_scale_density_squared = 0.0;
520 525  
521 526 // Return
522 527 return;
... ... @@ -534,7 +539,7 @@ void GModelSpatialRadialProfileDMEinasto::copy_members(const GModelSpatialRadial
534 539 {
535 540 // Copy members. We do not have to push back the members on the parameter
536 541 // stack as this should have been done by init_members() that was called
537   - // before. Otherwise we would have sigma twice on the stack.
  542 + // before.
538 543 m_scale_radius = model.m_scale_radius;
539 544 m_scale_density = model.m_scale_density;
540 545 m_halo_distance = model.m_halo_distance;
... ... @@ -572,11 +577,10 @@ void GModelSpatialRadialProfileDMEinasto::free_members(void)
572 577 ***************************************************************************/
573 578 double GModelSpatialRadialProfileDMEinasto::profile_value(const double& theta) const
574 579 {
575   -
576   - // update precomputation cache
  580 + // Update precomputation cache
577 581 update();
578 582  
579   - // initialize integral value
  583 + // Initialize integral value
580 584 double value = 0.0;
581 585  
582 586 // Set up integration limits
... ... @@ -597,11 +601,8 @@ double GModelSpatialRadialProfileDMEinasto::profile_value(const double&amp; theta) c
597 601 GIntegral integral(&integrand);
598 602 integral.max_iter(30);
599 603  
600   - // also play with eps() or max_iter(), bounds()
601   -
602   - // Set up integration boundaries
603   - // As there is usually an infinity at the halo center, this splits
604   - // the integral at the m_halo_distance.
  604 + // Set up integration boundaries. As there is usually an infinity at the
  605 + // halo center, this splits the integral at the m_halo_distance.
605 606 std::vector<double> bounds;
606 607 bounds.push_back(los_min);
607 608 bounds.push_back(los_max);
... ... @@ -609,11 +610,9 @@ double GModelSpatialRadialProfileDMEinasto::profile_value(const double&amp; theta) c
609 610  
610 611 // Compute value
611 612 value = integral.romberg(bounds);
612   -
613   - //std::cout << "profile_value=" << value << std::endl;
614 613  
615   - // multiply in the density^2
616   - value *= m_scale_density_squared ;
  614 + // Multiply in the density^2
  615 + value *= m_scale_density_squared;
617 616  
618 617 // Return value
619 618 return value;
... ... @@ -654,40 +653,33 @@ double GModelSpatialRadialProfileDMEinasto::profile_value(const double&amp; theta) c
654 653 ***************************************************************************/
655 654 double GModelSpatialRadialProfileDMEinasto::halo_kernel_los::eval( const double &los )
656 655 {
657   - // calculate the scale distance g, the ( distance from integration point
  656 + // Calculate the scale distance g, the ( distance from integration point
658 657 // to the halo center ) divided by ( the halo scale radius )
659 658  
660   - // first calculate the distance of the integration point from the halo
  659 + // First calculate the distance of the integration point from the halo
661 660 // center via the law of cosines
662   - double g = 0.0;
663   - g = los * los;
664   - g += m_halo_distance * m_halo_distance;
665   - g -= 2.0 * los * m_halo_distance * std::cos(m_theta);
666   - g = std::sqrt(g);
  661 + double g = los * los;
  662 + g += m_halo_distance * m_halo_distance;
  663 + g -= 2.0 * los * m_halo_distance * std::cos(m_theta);
  664 + g = std::sqrt(g);
667 665  
668   - // if we have a core radius specified, all halo values inside this core
  666 + // If we have a core radius specified, all halo values inside this core
669 667 // radius should be the same as at the core radius itself.
670   - if ( g < m_core_radius )
671   - {
672   - //std::cout << "los=" << los << " is inside core radius " << m_core_radius << std::endl;
673   - g = m_core_radius ;
  668 + if (g < m_core_radius) {
  669 + g = m_core_radius;
674 670 }
675 671  
676   - // finish scaling the integration point by the halo's scale radius
  672 + // Finish scaling the integration point by the halo's scale radius
677 673 g /= m_scale_radius;
678 674  
679   - // calculate the halo density f at the scale distance g
680   - double f = 0.0;
681   - f = std::pow(g, m_alpha);
682   - f -= 1.0;
683   - f *= -2.0 / m_alpha;
684   - f = std::exp(f);
  675 + // Calculate the halo density f at the scale distance g
  676 + double f = std::pow(g, m_alpha);
  677 + f -= 1.0;
  678 + f *= -2.0 / m_alpha;
  679 + f = std::exp(f);
685 680  
686   - // squared, for annihilating dm
687   - // would just be f if it was decaying dm
688   - f = f * f ;
689   -
690   - //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;
  681 + // Squared, for annihilating DM, would just be f if it was decaying DM
  682 + f = f * f;
691 683  
692 684 // Return function value
693 685 return f;
... ... @@ -703,16 +695,16 @@ double GModelSpatialRadialProfileDMEinasto::halo_kernel_los::eval( const double
703 695 ***************************************************************************/
704 696 void GModelSpatialRadialProfileDMEinasto::update() const
705 697 {
706   -
707 698 // Update if scale radius has changed
708   - if (m_last_scale_radius != scale_radius() || m_last_scale_density != scale_density()) {
  699 + if (m_last_scale_radius != scale_radius() ||
  700 + m_last_scale_density != scale_density()) {
709 701  
710 702 // Store last values
711 703 m_last_scale_radius = scale_radius();
712 704 m_last_scale_density = scale_density();
713 705  
714   - // perform precomputations
715   - m_mass_radius = 10.0 * scale_radius();
  706 + // Perform precomputations
  707 + m_mass_radius = 10.0 * scale_radius();
716 708 m_scale_density_squared = scale_density() * scale_density();
717 709  
718 710 }
... ... @@ -723,68 +715,63 @@ void GModelSpatialRadialProfileDMEinasto::update() const
723 715  
724 716  
725 717 /***********************************************************************//**
726   - * @brief Calculate Halo Mass Density
  718 + * @brief Calculate halo mass density
727 719 *
728 720 * @param[in] radius from halo center (kpc)
  721 + * @return Halo mass density.
729 722 *
730   - * Calculates the halo's mass density at a given radial distance from the halo
731   - * center.
732   - *
  723 + * Calculates the halo's mass density at a given radial distance from the
  724 + * halo center.
733 725 ***************************************************************************/
734 726 double GModelSpatialRadialProfileDMEinasto::mass_density( const double& radius ) const
735 727 {
736   - double density = 0.0 ;
737   -
738   - //m_halo_distance ;
739   - halo_kernel_los halo_shape( m_scale_radius.value(),
740   - m_halo_distance.value(),
741   - m_alpha.value(),
742   - 0.0,
743   - m_core_radius.value() ) ;
  728 + // Set-up kernel
  729 + halo_kernel_los halo_shape(m_scale_radius.value(),
  730 + m_halo_distance.value(),
  731 + m_alpha.value(),
  732 + 0.0,
  733 + m_core_radius.value());
744 734  
745 735  
746   - // eval produces a unitless density^2, so we must take its square root
747   - density = std::sqrt( halo_shape.eval( m_halo_distance.value() + radius ) ) ;
  736 + // eval produces a unitless density^2, so we must take its square root
  737 + double density = std::sqrt(halo_shape.eval(m_halo_distance.value() + radius));
748 738  
749   - // multiply in the missing scale density
750   - density *= m_scale_density.value() ;
  739 + // Multiply in the missing scale density
  740 + density *= m_scale_density.value();
751 741  
752   - return density ;
  742 + // Return halo mass density
  743 + return density;
753 744 }
754 745  
  746 +
755 747 /***********************************************************************//**
756   - * @brief Calculate J Factor
  748 + * @brief Calculate J-factor
757 749 *
758 750 * @param[in] angle from halo center (radians)
  751 + * @return J-factor.
759 752 *
760   - * Calculates the halo's J-Factor at an angle from the halo center.
761   - *
  753 + * Calculates the halo's J-factor at an angle from the halo center.
762 754 ***************************************************************************/
763   -double GModelSpatialRadialProfileDMEinasto::jfactor( const double& angle ) const
  755 +double GModelSpatialRadialProfileDMEinasto::jfactor(const double& angle) const
764 756 {
765   -
766   - // Integration settings
767   - double minradian = 0.0 ;
768   - int npoints = 200 ;
769   -
770   - // initialize other variables
771   - double jfactor = 0.0 ;
772   - double dr = (angle - minradian) / npoints ;
773   - double r = 0.0;
774   -
775   - // loop over different radii in the profile
776   - for (int i = 0; i < npoints; ++i) {
777   -
778   - // integration: Int[ profile(r) * r * dr ]
779   - r = minradian + (i * dr);
780   - jfactor += profile_value(r) * r * dr;
781   -
782   - }
  757 + // Integration settings
  758 + double minradian = 0.0;
  759 + int npoints = 200;
  760 +
  761 + // Initialize other variables
  762 + double jfactor = 0.0 ;
  763 + double dr = (angle - minradian) / npoints;
  764 + double r = 0.0;
  765 +
  766 + // Loop over different radii in the profile
  767 + for (int i = 0; i < npoints; ++i) {
  768 + r = minradian + (i * dr);
  769 + jfactor += profile_value(r) * r * dr;
  770 + }
783 771  
784   - // J-Factor = 2 * pi * Int[ profile(r) * r * dr , {r,minradian,angle} ]
785   - jfactor *= gammalib::twopi;
  772 + // J-factor = 2 * pi * Int[ profile(r) * r * dr , {r,minradian,angle} ]
  773 + jfactor *= gammalib::twopi;
786 774  
787   - return jfactor ;
  775 + // Return J-factor
  776 + return jfactor;
788 777 }
789   -
790   -
... ...
src/model/GModelSpatialRadialProfileDMZhao.cpp
  View file @4e0202a
... ... @@ -28,7 +28,6 @@
28 28 #ifdef HAVE_CONFIG_H
29 29 #include <config.h>
30 30 #endif
31   -//#include <iomanip>
32 31 #include "GException.hpp"
33 32 #include "GTools.hpp"
34 33 #include "GMath.hpp"
... ... @@ -43,8 +42,8 @@ const GModelSpatialRadialProfileDMZhao g_radial_disk_seed;
43 42 const GModelSpatialRegistry g_radial_disk_registry(&g_radial_disk_seed);
44 43  
45 44 /* __ Method name definitions ____________________________________________ */
46   -#define G_READ "GModelSpatialRadialProfileDMZhao::read(GXmlElement&)"
47   -#define G_WRITE "GModelSpatialRadialProfileDMZhao::write(GXmlElement&)"
  45 +#define G_READ "GModelSpatialRadialProfileDMZhao::read(GXmlElement&)"
  46 +#define G_WRITE "GModelSpatialRadialProfileDMZhao::write(GXmlElement&)"
48 47  
49 48 /* __ Macros _____________________________________________________________ */
50 49  
... ... @@ -238,18 +237,18 @@ double GModelSpatialRadialProfileDMZhao::theta_min(void) const
238 237 ***************************************************************************/
239 238 double GModelSpatialRadialProfileDMZhao::theta_max(void) const
240 239 {
241   -
242   - // update precomputation cache
  240 + // Update precomputation cache
243 241 update();
244   -
  242 +
  243 + // Initialize maximum theta angle
245 244 double theta = 0.0;
246   -
  245 +
247 246 // If Earth is within the significant radius, then theta_max must
248 247 // contain the entire profile (180deg) ...
249 248 if (m_halo_distance.value() < m_mass_radius) {
250 249 theta = gammalib::pi;
251 250 }
252   -
  251 +
253 252 // ... otherwise, if the halo is far enough away (further than the
254 253 // significant radius) then we just need to deal with the angles within
255 254 // the sphere of the significant radius.
... ... @@ -258,10 +257,11 @@ double GModelSpatialRadialProfileDMZhao::theta_max(void) const
258 257 }
259 258  
260 259 // 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;
  260 + double theta_max = m_theta_max.value() * gammalib::deg2rad;
  261 + if (theta_max < theta) {
  262 + theta = theta_max;
263 263 }
264   -
  264 +
265 265 // Return value
266 266 return theta;
267 267 }
... ... @@ -273,26 +273,20 @@ double GModelSpatialRadialProfileDMZhao::theta_max(void) const
273 273 * @param[in] xml XML element.
274 274 *
275 275 * Reads the DMZhao radial profile model information from an XML element.
276   - * The XML element shall have either the format
  276 + * The XML element shall have the format
277 277 *
278 278 * <spatialModel type="DMZhaoProfile">
279   - * <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
280   - * <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
281   - * <parameter name="Sigma" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
282   - * <parameter name="Scale Radius" scale="1.0" value="21.5" min="1.0e-6" free="1"/>
283   - * <parameter name="Halo Distance" scale="1.0" value="7.94" min="1.0e-6" free="1"/>
284   - * <parameter name="Alpha" scale="1.0" value="0.17" min="0.01" max="10" free="1"/>
285   - * </spatialModel>
286   - *
287   - * or
288   - *
289   - * <spatialModel type="DMZhaoProfile">
290   - * <parameter name="GLON" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
291   - * <parameter name="GLAT" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
292   - * <parameter name="Sigma" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
293   - * <parameter name="Scale Radius" scale="1.0" value="21.5" min="1.0e-6" free="1"/>
294   - * <parameter name="Halo Distance" scale="1.0" value="7.94" min="1.0e-6" free="1"/>
295   - * <parameter name="Alpha" scale="1.0" value="0.17" min="0.01" max="10" free="1"/>
  279 + * <parameter name="RA" scale=.. value=.. min=.. max=.. free=../>
  280 + * <parameter name="DEC" scale=.. value=.. min=.. max=.. free=../>
  281 + * <parameter name="ScaleRadius" scale=.. value=.. min=.. max=.. free=../>
  282 + * <parameter name="ScaleDensity" scale=.. value=.. min=.. max=.. free=../>
  283 + * <parameter name="HaloDistance" scale=.. value=.. min=.. max=.. free=../>
  284 + * <parameter name="Alpha" scale=.. value=.. min=.. max=.. free=../>
  285 + * <parameter name="Beta" scale=.. value=.. min=.. max=.. free=../>
  286 + * <parameter name="Gamma" scale=.. value=.. min=.. max=.. free=../>
  287 + * <parameter name="ThetaMin" scale=.. value=.. min=.. max=.. free=../>
  288 + * <parameter name="ThetaMax" scale=.. value=.. min=.. max=.. free=../>
  289 + * <parameter name="CoreRadius" scale=.. value=.. min=.. max=.. free=../>
296 290 * </spatialModel>
297 291 ***************************************************************************/
298 292 void GModelSpatialRadialProfileDMZhao::read(const GXmlElement& xml)
... ... @@ -300,16 +294,16 @@ void GModelSpatialRadialProfileDMZhao::read(const GXmlElement&amp; xml)
300 294 // Read halo sky coordinate
301 295 GModelSpatialRadial::read(xml);
302 296  
303   - // Read Scale Radius parameter
304   - const GXmlElement* par1 = gammalib::xml_get_par(G_READ, xml, "Scale Radius");
  297 + // Read ScaleRadius parameter
  298 + const GXmlElement* par1 = gammalib::xml_get_par(G_READ, xml, "ScaleRadius");
305 299 m_scale_radius.read(*par1);
306 300  
307   - // Read Scale Density parameter
308   - const GXmlElement* par2 = gammalib::xml_get_par(G_READ, xml, "Scale Density");
  301 + // Read ScaleDensity parameter
  302 + const GXmlElement* par2 = gammalib::xml_get_par(G_READ, xml, "ScaleDensity");
309 303 m_scale_density.read(*par2);
310 304  
311   - // Read Halo Distance parameter
312   - const GXmlElement* par3 = gammalib::xml_get_par(G_READ, xml, "Halo Distance");
  305 + // Read HaloDistance parameter
  306 + const GXmlElement* par3 = gammalib::xml_get_par(G_READ, xml, "HaloDistance");
313 307 m_halo_distance.read(*par3);
314 308  
315 309 // Read Alpha parameter
... ... @@ -324,16 +318,16 @@ void GModelSpatialRadialProfileDMZhao::read(const GXmlElement&amp; xml)
324 318 const GXmlElement* par6 = gammalib::xml_get_par(G_READ, xml, "Gamma");
325 319 m_gamma.read(*par6);
326 320  
327   - // Read Theta Min parameter
328   - const GXmlElement* par7 = gammalib::xml_get_par(G_READ, xml, "Theta Min");
  321 + // Read ThetaMin parameter
  322 + const GXmlElement* par7 = gammalib::xml_get_par(G_READ, xml, "ThetaMin");
329 323 m_theta_min.read(*par7);
330 324  
331   - // Read Theta Max parameter
332   - const GXmlElement* par8 = gammalib::xml_get_par(G_READ, xml, "Theta Max");
  325 + // Read ThetaMax parameter
  326 + const GXmlElement* par8 = gammalib::xml_get_par(G_READ, xml, "ThetaMax");
333 327 m_theta_max.read(*par8);
334 328  
335   - // Read Core Radius parameter
336   - const GXmlElement* par9 = gammalib::xml_get_par(G_READ, xml, "Core Radius");
  329 + // Read CoreRadius parameter
  330 + const GXmlElement* par9 = gammalib::xml_get_par(G_READ, xml, "CoreRadius");
337 331 m_core_radius.read(*par9);
338 332  
339 333 // Return
... ... @@ -350,9 +344,17 @@ void GModelSpatialRadialProfileDMZhao::read(const GXmlElement&amp; xml)
350 344 * The XML element will have the format
351 345 *
352 346 * <spatialModel type="DMZhaoProfile">
353   - * <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
354   - * <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
355   - * <parameter name="Sigma" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
  347 + * <parameter name="RA" scale=.. value=.. min=.. max=.. free=../>
  348 + * <parameter name="DEC" scale=.. value=.. min=.. max=.. free=../>
  349 + * <parameter name="ScaleRadius" scale=.. value=.. min=.. max=.. free=../>
  350 + * <parameter name="ScaleDensity" scale=.. value=.. min=.. max=.. free=../>
  351 + * <parameter name="HaloDistance" scale=.. value=.. min=.. max=.. free=../>
  352 + * <parameter name="Alpha" scale=.. value=.. min=.. max=.. free=../>
  353 + * <parameter name="Beta" scale=.. value=.. min=.. max=.. free=../>
  354 + * <parameter name="Gamma" scale=.. value=.. min=.. max=.. free=../>
  355 + * <parameter name="ThetaMin" scale=.. value=.. min=.. max=.. free=../>
  356 + * <parameter name="ThetaMax" scale=.. value=.. min=.. max=.. free=../>
  357 + * <parameter name="CoreRadius" scale=.. value=.. min=.. max=.. free=../>
356 358 * </spatialModel>
357 359 ***************************************************************************/
358 360 void GModelSpatialRadialProfileDMZhao::write(GXmlElement& xml) const
... ... @@ -360,16 +362,16 @@ void GModelSpatialRadialProfileDMZhao::write(GXmlElement&amp; xml) const
360 362 // Write DMZhao location
361 363 GModelSpatialRadial::write(xml);
362 364  
363   - // Write Scale Radius parameter
364   - GXmlElement* par1 = gammalib::xml_need_par(G_WRITE, xml, "Scale Radius");
  365 + // Write ScaleRadius parameter
  366 + GXmlElement* par1 = gammalib::xml_need_par(G_WRITE, xml, "ScaleRadius");
365 367 m_scale_radius.write(*par1);
366 368  
367   - // Write Scale Density parameter
368   - GXmlElement* par2 = gammalib::xml_need_par(G_WRITE, xml, "Scale Density");
  369 + // Write ScaleDensity parameter
  370 + GXmlElement* par2 = gammalib::xml_need_par(G_WRITE, xml, "ScaleDensity");
369 371 m_scale_density.write(*par2);
370 372  
371   - // Write Halo Distance parameter
372   - GXmlElement* par3 = gammalib::xml_need_par(G_WRITE, xml, "Halo Distance");
  373 + // Write HaloDistance parameter
  374 + GXmlElement* par3 = gammalib::xml_need_par(G_WRITE, xml, "HaloDistance");
373 375 m_halo_distance.write(*par3);
374 376  
375 377 // Write Alpha parameter
... ... @@ -384,16 +386,16 @@ void GModelSpatialRadialProfileDMZhao::write(GXmlElement&amp; xml) const
384 386 GXmlElement* par6 = gammalib::xml_need_par(G_WRITE, xml, "Gamma");
385 387 m_gamma.write(*par6);
386 388  
387   - // Write Theta Max parameter
388   - GXmlElement* par7 = gammalib::xml_need_par(G_WRITE, xml, "Theta Min");
  389 + // Write ThetaMax parameter
  390 + GXmlElement* par7 = gammalib::xml_need_par(G_WRITE, xml, "ThetaMin");
389 391 m_theta_min.write(*par7);
390 392  
391   - // Write Theta Max parameter
392   - GXmlElement* par8 = gammalib::xml_need_par(G_WRITE, xml, "Theta Max");
  393 + // Write ThetaMax parameter
  394 + GXmlElement* par8 = gammalib::xml_need_par(G_WRITE, xml, "ThetaMax");
393 395 m_theta_max.write(*par8);
394 396  
395   - // Write Core Radius parameter
396   - GXmlElement* par9 = gammalib::xml_need_par(G_WRITE, xml, "Core Radius");
  397 + // Write CoreRadius parameter
  398 + GXmlElement* par9 = gammalib::xml_need_par(G_WRITE, xml, "CoreRadius");
397 399 m_core_radius.write(*par9);
398 400  
399 401 // Return
... ... @@ -404,7 +406,7 @@ void GModelSpatialRadialProfileDMZhao::write(GXmlElement&amp; xml) const
404 406 /***********************************************************************//**
405 407 * @brief Print information
406 408 *
407   - * @param[in] chatter Chattiness (defaults to NORMAL).
  409 + * @param[in] chatter Chattiness.
408 410 * @return String containing model information.
409 411 ***************************************************************************/
410 412 std::string GModelSpatialRadialProfileDMZhao::print(const GChatter& chatter) const
... ... @@ -446,20 +448,20 @@ void GModelSpatialRadialProfileDMZhao::init_members(void)
446 448  
447 449 // Initialise scale radius
448 450 m_scale_radius.clear();
449   - m_scale_radius.name("Scale Radius");
  451 + m_scale_radius.name("ScaleRadius");
450 452 m_scale_radius.unit("kpc");
451   - m_scale_radius.value(21.5); // default to GC scale radius
452   - m_scale_radius.min(1.0e-6); // arbitrarily chosen :/
  453 + m_scale_radius.value(21.5); // default to GC scale radius
  454 + m_scale_radius.min(1.0e-6); // arbitrarily chosen :/
453 455 m_scale_radius.free();
454 456 m_scale_radius.scale(1.0);
455 457 m_scale_radius.gradient(0.0);
456   - m_scale_radius.has_grad(false); // Radial components never have gradients
  458 + m_scale_radius.has_grad(false); // Radial components never have gradients
457 459  
458 460 // Initialise scale density
459 461 m_scale_density.clear();
460   - m_scale_density.name("Scale Density");
  462 + m_scale_density.name("ScaleDensity");
461 463 m_scale_density.unit("GeV/cm^3");
462   - m_scale_density.value(0.2); // GeV/cm^3 default to GC scale density
  464 + m_scale_density.value(0.2); // GeV/cm^3 default to GC scale density
463 465 m_scale_density.free();
464 466 m_scale_density.scale(1.0);
465 467 m_scale_density.gradient(0.0);
... ... @@ -467,76 +469,76 @@ void GModelSpatialRadialProfileDMZhao::init_members(void)
467 469  
468 470 // Initialise halo distance
469 471 m_halo_distance.clear();
470   - m_halo_distance.name("Halo Distance");
  472 + m_halo_distance.name("HaloDistance");
471 473 m_halo_distance.unit("kpc");
472   - m_halo_distance.value(7.94); // default to GC halo distance
473   - m_halo_distance.min(1.0e-6); // arbitrarily chosen
  474 + m_halo_distance.value(7.94); // default to GC halo distance
  475 + m_halo_distance.min(1.0e-6); // arbitrarily chosen
474 476 m_halo_distance.free();
475 477 m_halo_distance.scale(1.0);
476 478 m_halo_distance.gradient(0.0);
477   - m_halo_distance.has_grad(false); // Radial components never have gradients
  479 + m_halo_distance.has_grad(false); // Radial components never have gradients
478 480  
479 481 // Initialise alpha
480 482 m_alpha.clear();
481 483 m_alpha.name("Alpha");
482 484 m_alpha.unit("unitless");
483   - m_alpha.value(1.0); // default to NFW alpha
484   - m_alpha.min(0.01); // arbitrarily chosen
485   - m_alpha.max(10.0); // arbitrarily chosen
  485 + m_alpha.value(1.0); // default to NFW alpha
  486 + m_alpha.min(0.01); // arbitrarily chosen
  487 + m_alpha.max(10.0); // arbitrarily chosen
486 488 m_alpha.free();
487 489 m_alpha.scale(1.0);
488 490 m_alpha.gradient(0.0);
489   - m_alpha.has_grad(false); // Radial components never have gradients
  491 + m_alpha.has_grad(false); // Radial components never have gradients
490 492  
491 493 // Initialise
492 494 m_beta.clear();
493 495 m_beta.name("Beta");
494 496 m_beta.unit("unitless");
495   - m_beta.value(3.0); // default to NFW beta
496   - m_beta.min(0.01); // arbitrarily chosen
497   - m_beta.max(10.0); // arbitrarily chosen
  497 + m_beta.value(3.0); // default to NFW beta
  498 + m_beta.min(0.01); // arbitrarily chosen
  499 + m_beta.max(10.0); // arbitrarily chosen
498 500 m_beta.free();
499 501 m_beta.scale(1.0);
500 502 m_beta.gradient(0.0);
501   - m_beta.has_grad(false); // Radial components never have gradients
  503 + m_beta.has_grad(false); // Radial components never have gradients
502 504  
503 505 // Initialise
504 506 m_gamma.clear();
505 507 m_gamma.name("Gamma");
506 508 m_gamma.unit("unitless");
507   - m_gamma.value(1.0); // default to NFW gamma
508   - m_gamma.min(0.01); // arbitrarily chosen
509   - m_gamma.max(10.0); // arbitrarily chosen
  509 + m_gamma.value(1.0); // default to NFW gamma
  510 + m_gamma.min(0.01); // arbitrarily chosen
  511 + m_gamma.max(10.0); // arbitrarily chosen
510 512 m_gamma.free();
511 513 m_gamma.scale(1.0);
512 514 m_gamma.gradient(0.0);
513   - m_gamma.has_grad(false); // Radial components never have gradients
  515 + m_gamma.has_grad(false); // Radial components never have gradients
514 516  
515 517 // Initialise theta max
516 518 m_theta_min.clear();
517   - m_theta_min.name("Theta Min");
  519 + m_theta_min.name("ThetaMin");
518 520 m_theta_min.unit("degrees");
519   - m_theta_min.value( 1.0e-6 ); // can only go from halo center to opposite halo center
520   - m_theta_min.min(1.0e-10); // arbitrarily chosen, some halos don't converge at theta=0.0
521   - m_theta_min.fix(); // should always be fixed!
  521 + m_theta_min.value(1.0e-6); // can only go from halo center to opposite halo center
  522 + m_theta_min.min(1.0e-10); // arbitrarily chosen, some halos don't converge at theta=0.0
  523 + m_theta_min.fix(); // should always be fixed!
522 524 m_theta_min.scale(1.0);
523 525 m_theta_min.gradient(0.0);
524   - m_theta_min.has_grad(false); // Radial components never have gradients
  526 + m_theta_min.has_grad(false); // Radial components never have gradients
525 527  
526 528 // Initialise theta max
527 529 m_theta_max.clear();
528   - m_theta_max.name("Theta Max");
  530 + m_theta_max.name("ThetaMax");
529 531 m_theta_max.unit("degrees");
530   - m_theta_max.value( 180.0 ); // can only go from halo center to opposite halo center
531   - m_theta_max.min(1.0e-6); // arbitrarily chosen
532   - m_theta_max.fix(); // should always be fixed!
  532 + m_theta_max.value(180.0); // can only go from halo center to opposite halo center
  533 + m_theta_max.min(1.0e-6); // arbitrarily chosen
  534 + m_theta_max.fix(); // should always be fixed!
533 535 m_theta_max.scale(1.0);
534 536 m_theta_max.gradient(0.0);
535   - m_theta_max.has_grad(false); // Radial components never have gradients
  537 + m_theta_max.has_grad(false); // Radial components never have gradients
536 538  
537 539 // Initialise core radius
538 540 m_core_radius.clear();
539   - m_core_radius.name("Core Radius");
  541 + m_core_radius.name("CoreRadius");
540 542 m_core_radius.unit("kpc");
541 543 m_core_radius.value(10.0);
542 544 m_core_radius.min(0.0);
... ... @@ -558,10 +560,10 @@ void GModelSpatialRadialProfileDMZhao::init_members(void)
558 560  
559 561 // Initialize precomputation cache. Note that zero values flag
560 562 // uninitialised, as a zero radius is not meaningful
561   - m_last_scale_radius = 0.0 ;
562   - m_last_scale_density = 0.0 ;
563   - m_mass_radius = 0.0 ;
564   - m_scale_density_squared = 0.0 ;
  563 + m_last_scale_radius = 0.0;
  564 + m_last_scale_density = 0.0;
  565 + m_mass_radius = 0.0;
  566 + m_scale_density_squared = 0.0;
565 567  
566 568 // Return
567 569 return;
... ... @@ -579,7 +581,7 @@ void GModelSpatialRadialProfileDMZhao::copy_members(const GModelSpatialRadialPro
579 581 {
580 582 // Copy members. We do not have to push back the members on the parameter
581 583 // stack as this should have been done by init_members() that was called
582   - // before. Otherwise we would have sigma twice on the stack.
  584 + // before.
583 585 m_scale_radius = model.m_scale_radius;
584 586 m_scale_density = model.m_scale_density;
585 587 m_halo_distance = model.m_halo_distance;
... ... @@ -619,17 +621,16 @@ void GModelSpatialRadialProfileDMZhao::free_members(void)
619 621 ***************************************************************************/
620 622 double GModelSpatialRadialProfileDMZhao::profile_value(const double& theta) const
621 623 {
622   -
623   - // update precomputation cache
  624 + // Update precomputation cache
624 625 update();
625   -
626   - // initialize integral value
  626 +
  627 + // Initialize integral value
627 628 double value = 0.0;
628   -
  629 +
629 630 // Set up integration limits
630 631 double los_min = m_halo_distance.value() - m_mass_radius;
631 632 double los_max = m_halo_distance.value() + m_mass_radius;
632   -
  633 +
633 634 // Set up integral
634 635 halo_kernel_los integrand(m_scale_radius.value(),
635 636 m_halo_distance.value(),
... ... @@ -640,7 +641,7 @@ double GModelSpatialRadialProfileDMZhao::profile_value(const double&amp; theta) cons
640 641 m_core_radius.value());
641 642 GIntegral integral(&integrand);
642 643 integral.max_iter(25);
643   -
  644 +
644 645 // Set up integration boundaries
645 646 // As there is usually an infinity at the halo center, this splits
646 647 // the integral at the m_halo_distance.
... ... @@ -651,10 +652,10 @@ double GModelSpatialRadialProfileDMZhao::profile_value(const double&amp; theta) cons
651 652  
652 653 // Compute value
653 654 value = integral.romberg(bounds);
654   -
  655 +
655 656 // apply scale density squared
656 657 // TODO: must be multiplied by the particle physics factor
657   - value *= m_scale_density_squared ;
  658 + value *= m_scale_density_squared;
658 659  
659 660 // Return value
660 661 return value;
... ... @@ -665,6 +666,7 @@ double GModelSpatialRadialProfileDMZhao::profile_value(const double&amp; theta) cons
665 666 * @brief Kernel for zhao halo density profile squared
666 667 *
667 668 * @param[in] distance from observer to point in space (meters)
  669 + * @return Jalo density.
668 670 *
669 671 * Computes the value of a zhao halo density profile squared,
670 672 * at distance l from observer, at angle \f[\theta\f] from the halo center,
... ... @@ -691,46 +693,38 @@ double GModelSpatialRadialProfileDMZhao::profile_value(const double&amp; theta) cons
691 693 * Monthly Notices of the Royal Astronomical Society, 278, 488-49
692 694 * http://mnras.oxfordjournals.org/content/278/2/488.short
693 695 *
694   - * @return unit
695 696 *
696 697 ***************************************************************************/
697   -double GModelSpatialRadialProfileDMZhao::halo_kernel_los::eval( const double &los )
  698 +double GModelSpatialRadialProfileDMZhao::halo_kernel_los::eval(const double &los)
698 699 {
699   - // calculate the scale distance g, the ( distance from integration point
  700 + // Calculate the scale distance g, the ( distance from integration point
700 701 // to the halo center ) divided by ( the halo scale radius )
701 702  
702   - // first calculate the distance of the integration point from the halo
  703 + // First calculate the distance of the integration point from the halo
703 704 // center via the law of cosines
704   - double g = 0.0;
705   - g = los * los;
706   - g += m_halo_distance * m_halo_distance;
707   - g -= 2.0 * los * m_halo_distance * std::cos(m_theta);
708   - g = std::sqrt(g) ;
  705 + double g = los * los;
  706 + g += m_halo_distance * m_halo_distance;
  707 + g -= 2.0 * los * m_halo_distance * std::cos(m_theta);
  708 + g = std::sqrt(g);
709 709  
710   - // if we have a core radius specified, all halo values inside this core
  710 + // If we have a core radius specified, all halo values inside this core
711 711 // radius should be the same as at the core radius itself.
712   - if ( g < m_core_radius )
713   - {
714   - //std::cout << "los=" << los << " is inside core radius " << m_core_radius << std::endl;
715   - g = m_core_radius ;
  712 + if (g < m_core_radius) {
  713 + g = m_core_radius;
716 714 }
717 715  
718   - // finish scaling the integration point by the halo's scale radius
  716 + // Finish scaling the integration point by the halo's scale radius
719 717 g /= m_scale_radius ;
720 718  
721 719 // apply the zhao profile from the scaled distance to the halo center
722   - double f = 0.0 ;
723   - f = std::pow(g, m_alpha);
724   - f += 1 ;
725   - f = std::pow(f, (m_beta-m_gamma)/m_alpha);
726   - f *= std::pow(g, m_gamma);
727   - f = 1 / f ;
728   -
729   - // squared, for annihilating dm
730   - // would just be f if it was decaying dm
  720 + double f = std::pow(g, m_alpha);
  721 + f += 1.0;
  722 + f = std::pow(f, (m_beta-m_gamma)/m_alpha);
  723 + f *= std::pow(g, m_gamma);
  724 + f = 1.0 / f;
  725 +
  726 + // Squared, for annihilating DM, would just be f if it was decaying DM
731 727 f = f * f;
732   -
733   - //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 ;
734 728  
735 729 // Return function value
736 730 return f;
... ... @@ -746,14 +740,14 @@ double GModelSpatialRadialProfileDMZhao::halo_kernel_los::eval( const double &amp;lo
746 740 ***************************************************************************/
747 741 void GModelSpatialRadialProfileDMZhao::update() const
748 742 {
749   -
750 743 // Update if scale radius has changed
751   - if (m_last_scale_radius != scale_radius() || m_last_scale_density != scale_density() ) {
752   -
  744 + if (m_last_scale_radius != scale_radius() ||
  745 + m_last_scale_density != scale_density() ) {
  746 +
753 747 // Store last values
754 748 m_last_scale_radius = scale_radius();
755 749 m_last_scale_density = scale_density();
756   -
  750 +
757 751 // perform precomputations
758 752 m_mass_radius = 10.0 * scale_radius();
759 753 m_scale_density_squared = scale_density() * scale_density() ;
... ... @@ -766,18 +760,6 @@ void GModelSpatialRadialProfileDMZhao::update() const
766 760  
767 761  
768 762 /***********************************************************************//**
769   - * @brief Compute profile value
770   - *
771   - * @return Profile value
772   - ***************************************************************************/
773   -/*
774   -double GModelSpatialRadialProfileDMZhao::prof_val(const double& theta)
775   -{
776   - return this->profile_value(theta);
777   -}
778   -*/
779   -
780   -/***********************************************************************//**
781 763 * @brief Calculate Halo Mass Density
782 764 *
783 765 * @param[in] distance from halo center (kpc)
... ... @@ -788,57 +770,54 @@ double GModelSpatialRadialProfileDMZhao::prof_val(const double&amp; theta)
788 770 ***************************************************************************/
789 771 double GModelSpatialRadialProfileDMZhao::mass_density(const double& radius) const
790 772 {
791   - double density = 0.0 ;
792   -
793   - // initialize halo kernel with stored values
794   - // we use theta=0 to effectivly sample the halo radially
795   - halo_kernel_los halo_shape( m_scale_radius.value(),
796   - m_halo_distance.value(),
797   - m_alpha.value(),
798   - m_beta.value(),
799   - m_gamma.value(),
800   - 0.0,
801   - m_core_radius.value());
  773 + // Initialize halo kernel with stored values. We use theta=0 to effectivly
  774 + // sample the halo radially
  775 + halo_kernel_los halo_shape(m_scale_radius.value(),
  776 + m_halo_distance.value(),
  777 + m_alpha.value(),
  778 + m_beta.value(),
  779 + m_gamma.value(),
  780 + 0.0,
  781 + m_core_radius.value());
802 782  
803 783 // eval produces a unitless density^2, so we must take its square root
804   - density = std::sqrt( halo_shape.eval( m_halo_distance.value() + radius ) ) ;
  784 + double density = std::sqrt(halo_shape.eval(m_halo_distance.value() + radius));
805 785  
806   - // multiply in the missing scale density
807   - density *= m_scale_density.value() ;
  786 + // Multiply in the missing scale density
  787 + density *= m_scale_density.value();
808 788  
809 789 return density ;
810 790 }
811 791  
  792 +
812 793 /***********************************************************************//**
813   - * @brief Calculate J Factor
  794 + * @brief Calculate J-factor
814 795 *
815 796 * @param[in] angle from halo center (radians)
  797 + * @return J-factor
816 798 *
817   - * Calculates the halo's J-Factor at an angle from the halo center.
818   - *
  799 + * Calculates the halo's J-factor at an angle from the halo center.
819 800 ***************************************************************************/
820   -double GModelSpatialRadialProfileDMZhao::jfactor( const double& angle ) const
  801 +double GModelSpatialRadialProfileDMZhao::jfactor(const double& angle) const
821 802 {
822   - // Integration settings
823   - double minradian = 0.0 ;
824   - int npoints = 200 ;
  803 + // Integration settings
  804 + double minradian = 0.0;
  805 + int npoints = 200;
825 806  
826   - // initialize other variables
827   - double jfactor = 0.0 ;
828   - double dr = (angle - minradian) / npoints ;
829   - double r = 0.0;
  807 + // Initialize other variables
  808 + double jfactor = 0.0;
  809 + double dr = (angle - minradian) / npoints;
  810 + double r = 0.0;
830 811  
831   - // loop over different radii in the profile
832   - for (int i = 0; i < npoints; ++i) {
833   -
834   - // integration: Int[ profile(r) * r * dr ]
835   - r = minradian + (i * dr);
836   - jfactor += profile_value(r) * r * dr;
837   -
838   - }
839   -
840   - // J-Factor = 2 * pi * Int[ profile(r) * r * dr , {r,minradian,angle} ]
841   - jfactor *= gammalib::twopi;
  812 + // Loop over different radii in the profile
  813 + for (int i = 0; i < npoints; ++i) {
  814 + r = minradian + (i * dr);
  815 + jfactor += profile_value(r) * r * dr;
  816 + }
842 817  
843   - return jfactor ;
  818 + // J-factor = 2 * pi * Int[ profile(r) * r * dr , {r,minradian,angle} ]
  819 + jfactor *= gammalib::twopi;
  820 +
  821 + // Return J-factor
  822 + return jfactor;
844 823 }
... ...
src/model/Makefile.am
  View file @4e0202a
... ... @@ -39,10 +39,10 @@ sources = GModelPar.cpp \
39 39 GModelSpatialRadialGauss.cpp \
40 40 GModelSpatialRadialShell.cpp \
41 41 GModelSpatialRadialProfile.cpp \
  42 + GModelSpatialRadialProfileGauss.cpp \
42 43 GModelSpatialRadialProfileDMZhao.cpp \
43 44 GModelSpatialRadialProfileDMBurkert.cpp \
44 45 GModelSpatialRadialProfileDMEinasto.cpp \
45   - GModelSpatialRadialProfileGauss.cpp \
46 46 GModelSpatialElliptical.cpp \
47 47 GModelSpatialEllipticalDisk.cpp \
48 48 GModelSpatialEllipticalGauss.cpp \
... ...