GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpatialRadialProfileDMZhao.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpatialRadialProfileDMZhao.cpp - Zhao radial profile class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2016-2021 by Nathan Kelley-Hoskins *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GModelSpatialRadialProfileDMZhao.cpp
23  * @brief Radial DM Einasto profile model class implementation
24  * @author Nathan Kelley-Hoskins
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GException.hpp"
32 #include "GTools.hpp"
33 #include "GMath.hpp"
34 #include "GXmlElement.hpp"
37 
38 /* __ Constants __________________________________________________________ */
39 
40 /* __ Globals ____________________________________________________________ */
42 const GModelSpatialRegistry g_radial_disk_registry(&g_radial_disk_seed);
43 
44 /* __ Method name definitions ____________________________________________ */
45 #define G_READ "GModelSpatialRadialProfileDMZhao::read(GXmlElement&)"
46 #define G_WRITE "GModelSpatialRadialProfileDMZhao::write(GXmlElement&)"
47 
48 /* __ Macros _____________________________________________________________ */
49 
50 /* __ Coding definitions _________________________________________________ */
51 
52 /* __ Debug definitions __________________________________________________ */
53 
54 
55 /*==========================================================================
56  = =
57  = Constructors/destructors =
58  = =
59  ==========================================================================*/
60 
61 /***********************************************************************//**
62  * @brief Void constructor
63  *
64  * Constructs empty radial DMZhao profile
65  ***************************************************************************/
68 {
69  // Initialise members
70  init_members();
71 
72  // Return
73  return;
74 }
75 
76 
77 /***********************************************************************//**
78  * @brief XML constructor
79  *
80  * @param[in] xml XML element.
81  *
82  * Constructs radial DMZhao profile model by extracting information from
83  * an XML element. See the read() method for more information about the
84  * expected structure of the XML element.
85  ***************************************************************************/
88 {
89  // Initialise members
90  init_members();
91 
92  // Read information from XML element
93  read(xml);
94 
95  // Return
96  return;
97 }
98 
99 
100 /***********************************************************************//**
101  * @brief Copy constructor
102  *
103  * @param[in] model Radial DMZhao profile model.
104  *
105  * Copies radial DMZhao profile model from another radial profile model.
106  ***************************************************************************/
109 {
110  // Initialise members
111  init_members();
112 
113  // Copy members
114  copy_members(model);
115 
116  // Return
117  return;
118 }
119 
120 
121 /***********************************************************************//**
122  * @brief Destructor
123  *
124  * Destructs radial DMZhao profile model.
125  ***************************************************************************/
127 {
128  // Free members
129  free_members();
130 
131  // Return
132  return;
133 }
134 
135 
136 /*==========================================================================
137  = =
138  = Operators =
139  = =
140  ==========================================================================*/
141 
142 /***********************************************************************//**
143  * @brief Assignment operator
144  *
145  * @param[in] model Radial DMZhao profile model.
146  * @return Radial DMZhao profile model.
147  *
148  * Assigns radial DMZhao profile model.
149  ***************************************************************************/
151 {
152  // Execute only if object is not identical
153  if (this != &model) {
154 
155  // Copy base class members
157 
158  // Free members
159  free_members();
160 
161  // Initialise members
162  init_members();
163 
164  // Copy members
165  copy_members(model);
166 
167  } // endif: object was not identical
168 
169  // Return
170  return *this;
171 }
172 
173 
174 /*==========================================================================
175  = =
176  = Public methods =
177  = =
178  ==========================================================================*/
179 
180 /***********************************************************************//**
181  * @brief Clear radial DMZhao profile model
182  *
183  * Clears radial DMZhao profile model.
184  ***************************************************************************/
186 {
187  // Free class members
188  free_members();
192 
193  // Initialise members
197  init_members();
198 
199  // Return
200  return;
201 }
202 
203 
204 /***********************************************************************//**
205  * @brief Clone radial DMZhao profile model
206  *
207  * @return Pointer to deep copy of radial DMZhao profile model.
208  *
209  * Returns a deep copy of the radial DMZhao profile model.
210  ***************************************************************************/
212 {
213  // Clone radial disk model
214  return new GModelSpatialRadialProfileDMZhao(*this);
215 }
216 
217 /***********************************************************************//**
218  * @brief Return minimum model radius (in radians)
219  *
220  * @return Minimum model radius (in radians).
221  ***************************************************************************/
223 {
224 
225  // update precomputation cache
226  update();
227 
228  // Return value
229  return m_theta_min.value();
230 }
231 
232 
233 /***********************************************************************//**
234  * @brief Return maximum model radius (in radians)
235  *
236  * @return Maximum model radius (in radians).
237  ***************************************************************************/
239 {
240  // Update precomputation cache
241  update();
242 
243  // Initialize maximum theta angle
244  double theta = 0.0;
245 
246  // If Earth is within the significant radius, then theta_max must
247  // contain the entire profile (180deg) ...
249  theta = gammalib::pi;
250  }
251 
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 {
257  }
258 
259  // Always chose the lesser of ( mass_radius theta, theta_max )
261  if (theta_max < theta) {
262  theta = theta_max;
263  }
264 
265  // Return value
266  return theta;
267 }
268 
269 
270 /***********************************************************************//**
271  * @brief Read model from XML element
272  *
273  * @param[in] xml XML element.
274  *
275  * Reads the DMZhao radial profile model information from an XML element.
276  * The XML element shall have the format
277  *
278  * <spatialModel type="DMZhaoProfile">
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=../>
290  * </spatialModel>
291  ***************************************************************************/
293 {
294  // Verify number of model parameters
296 
297  // Read halo sky coordinate
299 
300  // Read ScaleRadius parameter
301  const GXmlElement* par1 = gammalib::xml_get_par(G_READ, xml, "ScaleRadius");
302  m_scale_radius.read(*par1);
303 
304  // Read ScaleDensity parameter
305  const GXmlElement* par2 = gammalib::xml_get_par(G_READ, xml, "ScaleDensity");
306  m_scale_density.read(*par2);
307 
308  // Read HaloDistance parameter
309  const GXmlElement* par3 = gammalib::xml_get_par(G_READ, xml, "HaloDistance");
310  m_halo_distance.read(*par3);
311 
312  // Read Alpha parameter
313  const GXmlElement* par4 = gammalib::xml_get_par(G_READ, xml, "Alpha");
314  m_alpha.read(*par4);
315 
316  // Read Beta parameter
317  const GXmlElement* par5 = gammalib::xml_get_par(G_READ, xml, "Beta");
318  m_beta.read(*par5);
319 
320  // Read Gamma parameter
321  const GXmlElement* par6 = gammalib::xml_get_par(G_READ, xml, "Gamma");
322  m_gamma.read(*par6);
323 
324  // Read ThetaMin parameter
325  const GXmlElement* par7 = gammalib::xml_get_par(G_READ, xml, "ThetaMin");
326  m_theta_min.read(*par7);
327 
328  // Read ThetaMax parameter
329  const GXmlElement* par8 = gammalib::xml_get_par(G_READ, xml, "ThetaMax");
330  m_theta_max.read(*par8);
331 
332  // Read CoreRadius parameter
333  const GXmlElement* par9 = gammalib::xml_get_par(G_READ, xml, "CoreRadius");
334  m_core_radius.read(*par9);
335 
336  // Return
337  return;
338 }
339 
340 
341 /***********************************************************************//**
342  * @brief Write model into XML element
343  *
344  * @param[in] xml XML element into which model information is written.
345  *
346  * Writes the DMZhao radial profile model information into an XML element.
347  * The XML element will have the format
348  *
349  * <spatialModel type="DMZhaoProfile">
350  * <parameter name="RA" scale=.. value=.. min=.. max=.. free=../>
351  * <parameter name="DEC" scale=.. value=.. min=.. max=.. free=../>
352  * <parameter name="ScaleRadius" scale=.. value=.. min=.. max=.. free=../>
353  * <parameter name="ScaleDensity" scale=.. value=.. min=.. max=.. free=../>
354  * <parameter name="HaloDistance" scale=.. value=.. min=.. max=.. free=../>
355  * <parameter name="Alpha" scale=.. value=.. min=.. max=.. free=../>
356  * <parameter name="Beta" scale=.. value=.. min=.. max=.. free=../>
357  * <parameter name="Gamma" scale=.. value=.. min=.. max=.. free=../>
358  * <parameter name="ThetaMin" scale=.. value=.. min=.. max=.. free=../>
359  * <parameter name="ThetaMax" scale=.. value=.. min=.. max=.. free=../>
360  * <parameter name="CoreRadius" scale=.. value=.. min=.. max=.. free=../>
361  * </spatialModel>
362  ***************************************************************************/
364 {
365  // Verify model type
367 
368  // Write DMZhao location
370 
371  // Write ScaleRadius parameter
372  GXmlElement* par1 = gammalib::xml_need_par(G_WRITE, xml, "ScaleRadius");
373  m_scale_radius.write(*par1);
374 
375  // Write ScaleDensity parameter
376  GXmlElement* par2 = gammalib::xml_need_par(G_WRITE, xml, "ScaleDensity");
377  m_scale_density.write(*par2);
378 
379  // Write HaloDistance parameter
380  GXmlElement* par3 = gammalib::xml_need_par(G_WRITE, xml, "HaloDistance");
381  m_halo_distance.write(*par3);
382 
383  // Write Alpha parameter
384  GXmlElement* par4 = gammalib::xml_need_par(G_WRITE, xml, "Alpha");
385  m_alpha.write(*par4);
386 
387  // Write Beta parameter
388  GXmlElement* par5 = gammalib::xml_need_par(G_WRITE, xml, "Beta");
389  m_beta.write(*par5);
390 
391  // Write Gamma parameter
392  GXmlElement* par6 = gammalib::xml_need_par(G_WRITE, xml, "Gamma");
393  m_gamma.write(*par6);
394 
395  // Write ThetaMax parameter
396  GXmlElement* par7 = gammalib::xml_need_par(G_WRITE, xml, "ThetaMin");
397  m_theta_min.write(*par7);
398 
399  // Write ThetaMax parameter
400  GXmlElement* par8 = gammalib::xml_need_par(G_WRITE, xml, "ThetaMax");
401  m_theta_max.write(*par8);
402 
403  // Write CoreRadius parameter
404  GXmlElement* par9 = gammalib::xml_need_par(G_WRITE, xml, "CoreRadius");
405  m_core_radius.write(*par9);
406 
407  // Return
408  return;
409 }
410 
411 
412 /***********************************************************************//**
413  * @brief Print information
414  *
415  * @param[in] chatter Chattiness.
416  * @return String containing model information.
417  ***************************************************************************/
418 std::string GModelSpatialRadialProfileDMZhao::print(const GChatter& chatter) const
419 {
420  // Initialise result string
421  std::string result;
422 
423  // Continue only if chatter is not silent
424  if (chatter != SILENT) {
425 
426  // Append header
427  result.append("=== GModelSpatialRadialProfileDMZhao ===");
428 
429  // Append parameters
430  result.append("\n"+gammalib::parformat("Number of parameters"));
431  result.append(gammalib::str(size()));
432  for (int i = 0; i < size(); ++i) {
433  result.append("\n"+m_pars[i]->print(chatter));
434  }
435 
436  } // endif: chatter was not silent
437 
438  // Return result
439  return result;
440 }
441 
442 
443 /*==========================================================================
444  = =
445  = Private methods =
446  = =
447  ==========================================================================*/
448 
449 /***********************************************************************//**
450  * @brief Initialise class members
451  ***************************************************************************/
453 {
454  // Initialise model type
455  m_type = "DMZhaoProfile";
456 
457  // Initialise scale radius
459  m_scale_radius.name("ScaleRadius");
460  m_scale_radius.unit("kpc");
461  m_scale_radius.value(21.5); // default to GC scale radius
462  m_scale_radius.min(1.0e-6); // arbitrarily chosen :/
464  m_scale_radius.scale(1.0);
466  m_scale_radius.has_grad(false); // Radial components never have gradients
467 
468  // Initialise scale density
470  m_scale_density.name("ScaleDensity");
471  m_scale_density.unit("GeV/cm^3");
472  m_scale_density.value(0.2); // GeV/cm^3 default to GC scale density
474  m_scale_density.scale(1.0);
476  m_scale_density.has_grad(false); // Radial components never have gradients
477 
478  // Initialise halo distance
480  m_halo_distance.name("HaloDistance");
481  m_halo_distance.unit("kpc");
482  m_halo_distance.value(7.94); // default to GC halo distance
483  m_halo_distance.min(1.0e-6); // arbitrarily chosen
485  m_halo_distance.scale(1.0);
487  m_halo_distance.has_grad(false); // Radial components never have gradients
488 
489  // Initialise alpha
490  m_alpha.clear();
491  m_alpha.name("Alpha");
492  m_alpha.unit("unitless");
493  m_alpha.value(1.0); // default to NFW alpha
494  m_alpha.min(0.01); // arbitrarily chosen
495  m_alpha.max(10.0); // arbitrarily chosen
496  m_alpha.free();
497  m_alpha.scale(1.0);
498  m_alpha.gradient(0.0);
499  m_alpha.has_grad(false); // Radial components never have gradients
500 
501  // Initialise
502  m_beta.clear();
503  m_beta.name("Beta");
504  m_beta.unit("unitless");
505  m_beta.value(3.0); // default to NFW beta
506  m_beta.min(0.01); // arbitrarily chosen
507  m_beta.max(10.0); // arbitrarily chosen
508  m_beta.free();
509  m_beta.scale(1.0);
510  m_beta.gradient(0.0);
511  m_beta.has_grad(false); // Radial components never have gradients
512 
513  // Initialise
514  m_gamma.clear();
515  m_gamma.name("Gamma");
516  m_gamma.unit("unitless");
517  m_gamma.value(1.0); // default to NFW gamma
518  m_gamma.min(0.01); // arbitrarily chosen
519  m_gamma.max(10.0); // arbitrarily chosen
520  m_gamma.free();
521  m_gamma.scale(1.0);
522  m_gamma.gradient(0.0);
523  m_gamma.has_grad(false); // Radial components never have gradients
524 
525  // Initialise theta max
526  m_theta_min.clear();
527  m_theta_min.name("ThetaMin");
528  m_theta_min.unit("degrees");
529  m_theta_min.value(1.0e-6); // can only go from halo center to opposite halo center
530  m_theta_min.min(1.0e-10); // arbitrarily chosen, some halos don't converge at theta=0.0
531  m_theta_min.fix(); // should always be fixed!
532  m_theta_min.scale(1.0);
533  m_theta_min.gradient(0.0);
534  m_theta_min.has_grad(false); // Radial components never have gradients
535 
536  // Initialise theta max
537  m_theta_max.clear();
538  m_theta_max.name("ThetaMax");
539  m_theta_max.unit("degrees");
540  m_theta_max.value(180.0); // can only go from halo center to opposite halo center
541  m_theta_max.min(1.0e-6); // arbitrarily chosen
542  m_theta_max.fix(); // should always be fixed!
543  m_theta_max.scale(1.0);
544  m_theta_max.gradient(0.0);
545  m_theta_max.has_grad(false); // Radial components never have gradients
546 
547  // Initialise core radius
549  m_core_radius.name("CoreRadius");
550  m_core_radius.unit("kpc");
551  m_core_radius.value(10.0);
552  m_core_radius.min(0.0);
553  m_core_radius.fix();
554  m_core_radius.scale(1.0);
555  m_core_radius.gradient(0.0);
556  m_core_radius.has_grad(false);
557 
558  // Set parameter pointer(s)
559  m_pars.push_back(&m_scale_radius);
560  m_pars.push_back(&m_scale_density);
561  m_pars.push_back(&m_halo_distance);
562  m_pars.push_back(&m_alpha);
563  m_pars.push_back(&m_beta);
564  m_pars.push_back(&m_gamma);
565  m_pars.push_back(&m_theta_min);
566  m_pars.push_back(&m_theta_max);
567  m_pars.push_back(&m_core_radius);
568 
569  // Initialize precomputation cache. Note that zero values flag
570  // uninitialised, as a zero radius is not meaningful
571  m_last_scale_radius = 0.0;
572  m_last_scale_density = 0.0;
573  m_mass_radius = 0.0;
575 
576  // Return
577  return;
578 }
579 
580 
581 /***********************************************************************//**
582  * @brief Copy class members
583  *
584  * @param[in] model Radial DMZhao model.
585  *
586  * Copies class members from another radial profile model.
587  ***************************************************************************/
589 {
590  // Copy members. We do not have to push back the members on the parameter
591  // stack as this should have been done by init_members() that was called
592  // before.
593  m_type = model.m_type; // Needed to conserve model type
597  m_alpha = model.m_alpha;
598  m_beta = model.m_beta;
599  m_gamma = model.m_gamma;
600  m_theta_min = model.m_theta_min;
601  m_theta_max = model.m_theta_max;
603 
604  // copy cache values
609 
610  // Return
611  return;
612 }
613 
614 
615 /***********************************************************************//**
616  * @brief Delete class members
617  ***************************************************************************/
619 {
620  // Return
621  return;
622 }
623 
624 
625 /***********************************************************************//**
626  * @brief Radial profile
627  *
628  * @param[in] theta Angular distance from DMZhao centre (radians).
629  * @return Profile value.
630  ***************************************************************************/
631 double GModelSpatialRadialProfileDMZhao::profile_value(const double& theta) const
632 {
633  // Update precomputation cache
634  update();
635 
636  // Set up integration limits
637  double los_min = m_halo_distance.value() - m_mass_radius;
638  double los_max = m_halo_distance.value() + m_mass_radius;
639 
640  // Set up integral
643  m_alpha.value(),
644  m_beta.value(),
645  m_gamma.value(),
646  theta,
647  m_core_radius.value());
648  GIntegral integral(&integrand);
649  integral.max_iter(25);
650 
651  // Set up integration boundaries
652  // As there is usually an infinity at the halo center, this splits
653  // the integral at the m_halo_distance.
654  std::vector<double> bounds;
655  bounds.push_back(los_min);
656  bounds.push_back(los_max);
657  bounds.push_back(m_halo_distance.value());
658 
659  // Compute value
660  double value = integral.romberg(bounds);
661 
662  // apply scale density squared
663  // TODO: must be multiplied by the particle physics factor
664  value *= m_scale_density_squared;
665 
666  // Return value
667  return value;
668 }
669 
670 
671 /***********************************************************************//**
672  * @brief Kernel for zhao halo density profile squared
673  *
674  * @param[in] los Distance from observer to point in space (meters)
675  * @return Halo density.
676  *
677  * Computes the value of a zhao halo density profile squared,
678  * at distance l from observer, at angle \f[\theta\f] from the halo center,
679  * with a halo posessing a scale radius of \f[r_s\f] :
680  *
681  * \f[
682  * f(\theta, l) = \frac{1}{ g^{\gamma} {\left( g^{\alpha} + 1 \right)}^{ \frac{\beta-\gamma}{\alpha} } }
683  * \f]
684  *
685  * where
686  *
687  * \f[
688  * g = \frac{ \sqrt{l^2+d^2-2ldCos(\theta)} }{r_s}
689  * \f]
690  *
691  * \f[ \beta \f] is the slope of the density at radii much bigger than \f[r_s\f]
692  * \f[ \gamma \f] is the slope of the density profile at radii much smaller than \f[r_s\f]
693  * \f[ \alpha \f] governs how large the transition region is between these two slopes.
694  * larger \f[ \alpha \f] means a smaller transition region.
695  *
696  * This profile is detailed in:
697  * Zhao, 1996
698  * "Analytical models for galactic nuclei"
699  * Monthly Notices of the Royal Astronomical Society, 278, 488-49
700  * http://mnras.oxfordjournals.org/content/278/2/488.short
701  ***************************************************************************/
703 {
704  // Calculate the scale distance g, the ( distance from integration point
705  // to the halo center ) divided by ( the halo scale radius )
706 
707  // First calculate the distance of the integration point from the halo
708  // center via the law of cosines
709  double g = los * los;
711  g -= 2.0 * los * m_halo_distance * std::cos(m_theta);
712  g = std::sqrt(g);
713 
714  // If we have a core radius specified, all halo values inside this core
715  // radius should be the same as at the core radius itself.
716  if (g < m_core_radius) {
717  g = m_core_radius;
718  }
719 
720  // Finish scaling the integration point by the halo's scale radius
721  g /= m_scale_radius ;
722 
723  // apply the zhao profile from the scaled distance to the halo center
724  double f = std::pow(g, m_alpha);
725  f += 1.0;
726  f = std::pow(f, (m_beta-m_gamma)/m_alpha);
727  f *= std::pow(g, m_gamma);
728  f = 1.0 / f;
729 
730  // Squared, for annihilating DM, would just be f if it was decaying DM
731  f = f * f;
732 
733  // Return function value
734  return f;
735 }
736 
737 
738 /***********************************************************************//**
739  * @brief Update precomputation cache
740  *
741  * Computes the m_mass_radius calculation, determining the radius around
742  * the halo that contains 99.99% of the mass. For a Zhao halo profile,
743  * this is just 10.0 * scale_radius .
744  ***************************************************************************/
746 {
747  // Update if scale radius has changed
750 
751  // Store last values
754 
755  // perform precomputations
756  m_mass_radius = 10.0 * scale_radius();
758 
759  }
760 
761  // Return
762  return;
763 }
764 
765 
766 /***********************************************************************//**
767  * @brief Calculate Halo Mass Density
768  *
769  * @param[in] radius Distance from halo center (kpc).
770  * @return Halo mass density.
771  *
772  * Calculates the halo's mass density at a given radial distance from the halo
773  * center.
774  ***************************************************************************/
775 double GModelSpatialRadialProfileDMZhao::mass_density(const double& radius) const
776 {
777  // Initialize halo kernel with stored values. We use theta=0 to effectivly
778  // sample the halo radially
779  halo_kernel_los halo_shape(m_scale_radius.value(),
781  m_alpha.value(),
782  m_beta.value(),
783  m_gamma.value(),
784  0.0,
785  m_core_radius.value());
786 
787  // eval produces a unitless density^2, so we must take its square root
788  double density = std::sqrt(halo_shape.eval(m_halo_distance.value() + radius));
789 
790  // Multiply in the missing scale density
791  density *= m_scale_density.value();
792 
793  return density ;
794 }
795 
796 
797 /***********************************************************************//**
798  * @brief Calculate J-factor
799  *
800  * @param[in] angle from halo center (radians)
801  * @return J-factor
802  *
803  * Calculates the halo's J-factor at an angle from the halo center.
804  ***************************************************************************/
806 {
807  // Integration settings
808  const double minradian = 0.0;
809  const int npoints = 200;
810 
811  // Initialize other variables
812  double jfactor = 0.0;
813  double dr = (angle - minradian) / npoints;
814 
815  // Loop over different radii in the profile
816  for (int i = 0; i < npoints; ++i) {
817  double r = minradian + (i * dr);
818  jfactor += profile_value(r) * r * dr;
819  }
820 
821  // J-factor = 2 * pi * Int[ profile(r) * r * dr , {r,minradian,angle} ]
822  jfactor *= gammalib::twopi;
823 
824  // Return J-factor
825  return jfactor;
826 }
virtual double theta_min(void) const
Return minimum model radius (in radians)
double scale_radius(void) const
Return scale radius.
GModelPar m_alpha
Power index, inverse transition region width.
GModelPar m_beta
Power index, slope at &gt;&gt; m_scale_radius.
const std::string & name(void) const
Return parameter name.
GModelPar m_gamma
Power index, slope at &lt;&lt; m_scale_radius.
virtual void read(const GXmlElement &xml)
Read model from XML element.
XML element node class interface definition.
const double pi
Definition: GMath.hpp:35
double gradient(void) const
Return parameter gradient.
int size(void) const
Return number of parameters.
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:351
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1190
virtual GModelSpatialRadialProfileDMZhao * clone(void) const
Clone radial DMZhao profile model.
GModelPar m_scale_radius
Scale radius of halo profile.
XML element node class.
Definition: GXmlElement.hpp:48
double max(void) const
Return parameter maximum boundary.
const GModelSpatialRadialDisk g_radial_disk_seed
virtual void clear(void)
Clear radial DMZhao profile model.
Gammalib tools definition.
GIntegral class interface definition.
Definition: GIntegral.hpp:46
Interface definition for the spatial model registry class.
double min(void) const
Return parameter minimum boundary.
void init_members(void)
Initialise class members.
virtual void write(GXmlElement &xml) const
Write model into XML element.
void init_members(void)
Initialise class members.
const double & scale(void) const
Return parameter scale.
const double deg2rad
Definition: GMath.hpp:43
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1358
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
void free_members(void)
Delete class members.
void free_members(void)
Delete class members.
void free(void)
Free a parameter.
double mass_density(const double &radius) const
Calculate Halo Mass Density.
std::vector< GModelPar * > m_pars
Parameter pointers.
void fix(void)
Fix a parameter.
GXmlElement * xml_need_par(const std::string &origin, GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1637
virtual void read(const GXmlElement &xml)
Read model from XML element.
void clear(void)
Clear parameter.
Spatial model registry class definition.
virtual GModelSpatialRadialProfileDMZhao & operator=(const GModelSpatialRadialProfileDMZhao &model)
Assignment operator.
Dark Matter Zhao profile model class interface definition.
GChatter
Definition: GTypemaps.hpp:33
double angle(const GVector &a, const GVector &b)
Computes angle between vectors.
Definition: GVector.cpp:974
virtual double profile_value(const double &theta) const
Radial profile.
std::string type(void) const
Return model type.
void xml_check_parnum(const std::string &origin, const GXmlElement &xml, const int &number)
Checks number of parameters.
Definition: GTools.cpp:1777
void copy_members(const GModelSpatialRadialProfileDMZhao &model)
Copy class members.
void init_members(void)
Initialise class members.
void update(void) const
Update precomputation cache.
virtual void write(GXmlElement &xml) const
Write model into XML element.
GModelPar m_scale_density
Scale density of halo profile.
void free_members(void)
Delete class members.
GModelPar m_halo_distance
Distance from earth to halo center.
std::string m_type
Spatial model type.
double jfactor(const double &angle) const
Calculate J-factor.
void read(const GXmlElement &xml)
Extract parameter attributes from XML element.
Definition: GModelPar.cpp:229
double value(void) const
Return parameter value.
virtual GModelSpatialRadialProfile & operator=(const GModelSpatialRadialProfile &model)
Assignment operator.
double scale_density(void) const
Return scale density.
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1422
const std::string & unit(void) const
Return parameter unit.
Exception handler interface definition.
virtual double theta_max(void) const
Return maximum model radius (in radians)
void init_members(void)
Initialise class members.
const double twopi
Definition: GMath.hpp:36
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
GVector atan(const GVector &vector)
Computes arctan of vector elements.
Definition: GVector.cpp:1148
double eval(const double &los)
Kernel for zhao halo density profile squared.
const GXmlElement * xml_get_par(const std::string &origin, const GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1689
virtual std::string print(const GChatter &chatter=NORMAL) const
Print information.
Mathematical function definitions.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489
void xml_check_type(const std::string &origin, GXmlElement &xml, const std::string &type)
Checks the model type.
Definition: GTools.cpp:1819