Karim Sabri / gammalib

Commit 61df6481cf48e730ab9e9146038f57cb92fa8ac1

Authored by Jürgen Knödlseder
1 parent 7a14eaae

Reimplement GModelSpatial::flux() method for diffuse and combined models (#1785)

Inline Side-by-side
Showing 21 changed files with 347 additions and 96 deletions
ChangeLog
  View file @61df648
... ... @@ -42,6 +42,7 @@
42 42 Add GSkyModel flux_error() and eflux_error() methods (#3074)
43 43 Remove COMPTEL internal response caching (#2973)
44 44 Add GSkyMap::flux(GRegion(s)) and GSkyMap::solidangle(GRegion(s)) methods (#1785)
  45 + Reimplement GModelSpatial::flux() method for diffuse and combined models (#1785)
45 46 Implement pre-computation cache in GModelSpatialGauss (#1300)
46 47  
47 48  
... ...
... ... @@ -47,6 +47,11 @@ The following methods have been added:
47 47 - GModelSky::flux_error(GSkyRegion&, GEnergy&, GEnergy&)
48 48 - GModelSky::eflux_error(GEnergy&, GEnergy&)
49 49 - GModelSky::eflux_error(GSkyRegion&, GEnergy&, GEnergy&)
  50 +- GModelSpatialComposite::flux(GSkyRegion&, GEnergy&, GTime&)
  51 +- GModelSpatialDiffuseConst::flux(GSkyRegion&, GEnergy&, GTime&)
  52 +- GModelSpatialDiffuseCube::flux(GSkyRegion&, GEnergy&, GTime&)
  53 +- GModelSpatialDiffuseMap::flux(GSkyRegion&, GEnergy&, GTime&)
  54 +- GModelSpatialPointSource::flux(GSkyRegion&, GEnergy&, GTime&)
50 55 - GModelSpatialRadial::is_energy_dependent()
51 56 - GModelSpatialRadial::is_time_dependent()
52 57 - GModelSpectralFunc::nodes()
... ... @@ -164,6 +169,17 @@ Added analytical gradients to the GModelSpatialRadialGauss class (#3203).
164 169 Added is_energy_dependent() and is_time_dependent() methods to GModelSpatialRadial
165 170 to allow for checking of energy or time dependence of spatial model components.
166 171  
  172 +Reimplemented the GModelSpatial::flux(GSkyRegion&, GEnergy&, GTime&) method at
  173 +the level of the following derived classes so that flux computation works for all
  174 +kind of regions:
  175 +- GModelSpatialComposite
  176 +- GModelSpatialDiffuseConst
  177 +- GModelSpatialDiffuseCube
  178 +- GModelSpatialDiffuseMap
  179 +- GModelSpatialPointSource
  180 +(#1785)
  181 +
  182 +
167 183 A pre-computation cache was added to the GModelSpatialGauss class (#1300).
168 184  
169 185  
... ...
doc/source/admin/release_history/2.0.rst
  View file @61df648
... ... @@ -58,6 +58,8 @@ Improvements
58 58 the radial Gaussian model in a stacked analysis
59 59 * [`1785 <https://cta-redmine.irap.omp.eu/issues/1785>`_] -
60 60 Add ``GSkyMap::flux(GRegion)``, ``GSkyMap::flux(GRegions)``, ``GSkyMap::solidangle(GRegion)`` and ``GSkyMap::solidangle(GRegions)`` methods
  61 +* [`1785 <https://cta-redmine.irap.omp.eu/issues/1785>`_] -
  62 + Reimplement ``GModelSpatial::flux()`` method for diffuse and combined models
61 63 * [`1300 <https://cta-redmine.irap.omp.eu/issues/1300>`_] -
62 64 Implement pre-computation cache in GModelSpatialGauss
63 65  
... ...
include/GModelSpatial.hpp
  View file @61df648
1 1 /***************************************************************************
2 2 * GModelSpatial.hpp - Spatial model abstract base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2021 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -104,6 +104,11 @@ public:
104 104 virtual void write(GXmlElement& xml) const = 0;
105 105 virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
106 106  
  107 + // Virtual methods
  108 + virtual double flux(const GSkyRegion& region,
  109 + const GEnergy& srcEng = GEnergy(),
  110 + const GTime& srcTime = GTime()) const;
  111 +
107 112 // Methods
108 113 std::string type(void) const;
109 114 void type(const std::string& type);
... ... @@ -113,9 +118,6 @@ public:
113 118 bool has_free_pars(void) const;
114 119 int size(void) const;
115 120 void autoscale(void);
116   - double flux(const GSkyRegion& region,
117   - const GEnergy& srcEng = GEnergy(),
118   - const GTime& srcTime = GTime()) const;
119 121 const GSkyRegion* region(void) const;
120 122  
121 123 protected:
... ...
include/GModelSpatialComposite.hpp
  View file @61df648
1 1 /***************************************************************************
2 2 * GModelSpatialComposite.hpp - Spatial composite model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016-2020 by Domenico Tiziani *
  4 + * copyright (C) 2016-2021 by Domenico Tiziani *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -75,6 +75,11 @@ public:
75 75 virtual void write(GXmlElement& xml) const;
76 76 virtual std::string print(const GChatter& chatter = NORMAL) const;
77 77  
  78 + // Overloaded base class methods
  79 + virtual double flux(const GSkyRegion& region,
  80 + const GEnergy& srcEng = GEnergy(),
  81 + const GTime& srcTime = GTime()) const;
  82 +
78 83 // Other methods
79 84 int components(void) const;
80 85 void append(const GModelSpatial& component,
... ...
include/GModelSpatialDiffuseConst.hpp
  View file @61df648
1 1 /***************************************************************************
2 2 * GModelSpatialDiffuseConst.hpp - Spatial isotropic model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2021 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -76,6 +76,11 @@ public:
76 76 virtual void write(GXmlElement& xml) const;
77 77 virtual std::string print(const GChatter& chatter = NORMAL) const;
78 78  
  79 + // Overloaded base class methods
  80 + virtual double flux(const GSkyRegion& region,
  81 + const GEnergy& srcEng = GEnergy(),
  82 + const GTime& srcTime = GTime()) const;
  83 +
79 84 // Other methods
80 85 double value(void) const;
81 86 void value(const double& value);
... ...
include/GModelSpatialDiffuseCube.hpp
  View file @61df648
1 1 /***************************************************************************
2 2 * GModelSpatialDiffuseCube.hpp - Spatial map cube model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2021 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -91,6 +91,11 @@ public:
91 91 virtual void write(GXmlElement& xml) const;
92 92 virtual std::string print(const GChatter& chatter = NORMAL) const;
93 93  
  94 + // Overloaded base class methods
  95 + virtual double flux(const GSkyRegion& region,
  96 + const GEnergy& srcEng = GEnergy(),
  97 + const GTime& srcTime = GTime()) const;
  98 +
94 99 // Other methods
95 100 double value(void) const;
96 101 void value(const double& value);
... ...
include/GModelSpatialDiffuseMap.hpp
  View file @61df648
1 1 /***************************************************************************
2 2 * GModelSpatialDiffuseMap.hpp - Spatial map model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2021 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -83,6 +83,11 @@ public:
83 83 virtual void write(GXmlElement& xml) const;
84 84 virtual std::string print(const GChatter& chatter = NORMAL) const;
85 85  
  86 + // Overloaded base class methods
  87 + virtual double flux(const GSkyRegion& region,
  88 + const GEnergy& srcEng = GEnergy(),
  89 + const GTime& srcTime = GTime()) const;
  90 +
86 91 // Other methods
87 92 double value(void) const;
88 93 void value(const double& value);
... ...
include/GModelSpatialPointSource.hpp
  View file @61df648
1 1 /***************************************************************************
2 2 * GModelSpatialPointSource.hpp - Spatial point source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2021 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -78,6 +78,11 @@ public:
78 78 virtual void write(GXmlElement& xml) const;
79 79 virtual std::string print(const GChatter& chatter = NORMAL) const;
80 80  
  81 + // Overloaded base class methods
  82 + virtual double flux(const GSkyRegion& region,
  83 + const GEnergy& srcEng = GEnergy(),
  84 + const GTime& srcTime = GTime()) const;
  85 +
81 86 // Other methods
82 87 double ra(void) const;
83 88 double dec(void) const;
... ...
pyext/GModelSpatial.i
  View file @61df648
1 1 /***************************************************************************
2 2 * GModelSpatial.i - Spatial model abstract base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2021 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -57,6 +57,11 @@ public:
57 57 virtual void read(const GXmlElement& xml) = 0;
58 58 virtual void write(GXmlElement& xml) const = 0;
59 59  
  60 + // Virtual methods
  61 + virtual double flux(const GSkyRegion& region,
  62 + const GEnergy& srcEng = GEnergy(),
  63 + const GTime& srcTime = GTime()) const;
  64 +
60 65 // Methods
61 66 std::string type(void) const;
62 67 void type(const std::string& type);
... ... @@ -65,9 +70,6 @@ public:
65 70 bool has_free_pars(void) const;
66 71 int size(void) const;
67 72 void autoscale(void);
68   - double flux(const GSkyRegion& region,
69   - const GEnergy& srcEng = GEnergy(),
70   - const GTime& srcTime = GTime()) const;
71 73 const GSkyRegion* region(void) const;
72 74 };
73 75  
... ...
pyext/GModelSpatialComposite.i
  View file @61df648
1 1 /***************************************************************************
2 2 * GModelSpatialComposite.i - Spatial point source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016-2020 by Domenico Tiziani *
  4 + * copyright (C) 2016-2021 by Domenico Tiziani *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -59,6 +59,11 @@ public:
59 59 virtual void read(const GXmlElement& xml);
60 60 virtual void write(GXmlElement& xml) const;
61 61  
  62 + // Overloaded base class methods
  63 + virtual double flux(const GSkyRegion& region,
  64 + const GEnergy& srcEng = GEnergy(),
  65 + const GTime& srcTime = GTime()) const;
  66 +
62 67 // Other methods
63 68 int components(void) const;
64 69 void append(const GModelSpatial& component,
... ...
pyext/GModelSpatialDiffuseConst.i
  View file @61df648
1 1 /***************************************************************************
2 2 * GModelSpatialDiffuseConst.i - Spatial isotropic model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2021 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -59,6 +59,11 @@ public:
59 59 virtual void read(const GXmlElement& xml);
60 60 virtual void write(GXmlElement& xml) const;
61 61  
  62 + // Overloaded base class methods
  63 + virtual double flux(const GSkyRegion& region,
  64 + const GEnergy& srcEng = GEnergy(),
  65 + const GTime& srcTime = GTime()) const;
  66 +
62 67 // Other methods
63 68 double value(void) const;
64 69 void value(const double& value);
... ...
pyext/GModelSpatialDiffuseCube.i
  View file @61df648
1 1 /***************************************************************************
2 2 * GModelSpatialDiffuseCube.i - Spatial map cube model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2021 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -64,6 +64,11 @@ public:
64 64 virtual void read(const GXmlElement& xml);
65 65 virtual void write(GXmlElement& xml) const;
66 66  
  67 + // Overloaded base class methods
  68 + virtual double flux(const GSkyRegion& region,
  69 + const GEnergy& srcEng = GEnergy(),
  70 + const GTime& srcTime = GTime()) const;
  71 +
67 72 // Other methods
68 73 double value(void) const;
69 74 void value(const double& value);
... ...
pyext/GModelSpatialDiffuseMap.i
  View file @61df648
1 1 /***************************************************************************
2 2 * GModelSpatialDiffuseMap.i - Spatial map model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2021 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -64,6 +64,11 @@ public:
64 64 virtual void read(const GXmlElement& xml);
65 65 virtual void write(GXmlElement& xml) const;
66 66  
  67 + // Overloaded base class methods
  68 + virtual double flux(const GSkyRegion& region,
  69 + const GEnergy& srcEng = GEnergy(),
  70 + const GTime& srcTime = GTime()) const;
  71 +
67 72 // Other methods
68 73 double value(void) const;
69 74 void value(const double& value);
... ...
pyext/GModelSpatialPointSource.i
  View file @61df648
1 1 /***************************************************************************
2 2 * GModelSpatialPointSource.i - Spatial point source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2021 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -61,6 +61,11 @@ public:
61 61 virtual void read(const GXmlElement& xml);
62 62 virtual void write(GXmlElement& xml) const;
63 63  
  64 + // Overloaded base class methods
  65 + virtual double flux(const GSkyRegion& region,
  66 + const GEnergy& srcEng = GEnergy(),
  67 + const GTime& srcTime = GTime()) const;
  68 +
64 69 // Other methods
65 70 double ra(void) const;
66 71 double dec(void) const;
... ...
src/model/GModelSpatial.cpp
  View file @61df648
... ... @@ -348,83 +348,65 @@ double GModelSpatial::flux(const GSkyRegion&amp; region,
348 348 const GSkyRegion* model_reg = this->region();
349 349 if (model_reg->overlaps(region)) {
350 350  
351   - // Check if the model is a point source
352   - const GModelSpatialPointSource* ps_model =
353   - dynamic_cast<const GModelSpatialPointSource*>(this);
354   -
355   - // If the model is not a point source integrate over circular
356   - // regions overlap
357   - if (ps_model == NULL) {
358   -
359   - // Thrown an exception if region is not a sky circle
360   - const GSkyRegionCircle* reg_circle =
361   - dynamic_cast<const GSkyRegionCircle*>(&region);
362   - if (reg_circle == NULL) {
363   - std::string msg = "Flux can only be computed for a circular "
364   - "region.";
365   - throw GException::feature_not_implemented(G_FLUX,msg);
366   - }
367   -
368   - // Thrown an exception if model region is not a sky circle
369   - const GSkyRegionCircle* model_circle =
370   - dynamic_cast<const GSkyRegionCircle*>(model_reg);
371   - if (model_circle == NULL) {
372   - std::string msg = "Flux can only be computed for spatial model "
373   - "with region defined as circle.";
374   - throw GException::feature_not_implemented(G_FLUX,msg);
375   - }
376   -
377   - // Model centre and radius in radians
378   - GSkyDir model_centre = model_circle->centre();
379   - double model_radius = model_circle->radius() * gammalib::deg2rad;
380   -
381   - // Distance between region and model centres in radians
382   - double distance = model_centre.dist(reg_circle->centre());
383   -
384   - // Take the distance between centers minus radius of model region as
385   - // the minimum radial integration boundary (in radians)
386   - double rho_min = distance - model_radius;
387   -
388   - // Make sure that rho_min does not become negative
389   - if (rho_min < 0.0) {
390   - rho_min = 0.0;
391   - }
  351 + // Throw an exception if region is not a sky circle
  352 + const GSkyRegionCircle* reg_circle = dynamic_cast<const GSkyRegionCircle*>(&region);
  353 + if (reg_circle == NULL) {
  354 + std::string msg = "Flux can only be computed for a circular "
  355 + "region.";
  356 + throw GException::feature_not_implemented(G_FLUX,msg);
  357 + }
392 358  
393   - // Take the region radius as the maximum radial integration boundary
394   - // (in radians)
395   - double rho_max = reg_circle->radius() * gammalib::deg2rad;
396   -
397   - // Pre-compute some quantities for arclength
398   - double cosdist = std::cos(distance);
399   - double sindist = std::sin(distance);
400   - double cosmodrad = std::cos(model_radius);
401   -
402   - // Setup integration kernel
403   - GModelSpatial::circle_int_kern_rho integrand(this,
404   - reg_circle,
405   - srcEng,
406   - srcTime,
407   - distance,
408   - cosdist,
409   - sindist,
410   - model_radius,
411   - cosmodrad);
412   - GIntegral integral(&integrand);
  359 + // Throw an exception if model region is not a sky circle
  360 + const GSkyRegionCircle* model_circle = dynamic_cast<const GSkyRegionCircle*>(model_reg);
  361 + if (model_circle == NULL) {
  362 + std::string msg = "Flux can only be computed for spatial model "
  363 + "with region defined as circle.";
  364 + throw GException::feature_not_implemented(G_FLUX,msg);
  365 + }
413 366  
414   - // Suppress integration warnings
415   - integral.silent(true);
  367 + // Model centre and radius in radians
  368 + GSkyDir model_centre = model_circle->centre();
  369 + double model_radius = model_circle->radius() * gammalib::deg2rad;
416 370  
417   - // Perform integration
418   - flux = integral.romberg(rho_min, rho_max);
  371 + // Distance between region and model centres in radians
  372 + double distance = model_centre.dist(reg_circle->centre());
419 373  
420   - } // endif: source model was not a point source
  374 + // Take the distance between centers minus radius of model region as
  375 + // the minimum radial integration boundary (in radians)
  376 + double rho_min = distance - model_radius;
421 377  
422   - // ... otherwise the model is point source that overlaps with the
423   - // region, hence the flux is set flux to 1.0
424   - else {
425   - flux = 1.0;
  378 + // Make sure that rho_min does not become negative
  379 + if (rho_min < 0.0) {
  380 + rho_min = 0.0;
426 381 }
427 382  
  383 + // Take the region radius as the maximum radial integration boundary
  384 + // (in radians)
  385 + double rho_max = reg_circle->radius() * gammalib::deg2rad;
  386 +
  387 + // Pre-compute some quantities for arclength
  388 + double cosdist = std::cos(distance);
  389 + double sindist = std::sin(distance);
  390 + double cosmodrad = std::cos(model_radius);
  391 +
  392 + // Setup integration kernel
  393 + GModelSpatial::circle_int_kern_rho integrand(this,
  394 + reg_circle,
  395 + srcEng,
  396 + srcTime,
  397 + distance,
  398 + cosdist,
  399 + sindist,
  400 + model_radius,
  401 + cosmodrad);
  402 + GIntegral integral(&integrand);
  403 +
  404 + // Suppress integration warnings
  405 + integral.silent(true);
  406 +
  407 + // Perform integration
  408 + flux = integral.romberg(rho_min, rho_max);
  409 +
428 410 } // endif: model overlaps with region
429 411  
430 412 // Return
... ...
src/model/GModelSpatialComposite.cpp
  View file @61df648
1 1 /***************************************************************************
2 2 * GModelSpatialComposite.cpp - Spatial composite model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016-2020 by Domenico Tiziani *
  4 + * copyright (C) 2016-2021 by Domenico Tiziani *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -713,6 +713,39 @@ double GModelSpatialComposite::sum_of_scales(void) const
713 713  
714 714  
715 715 /***********************************************************************//**
  716 + * @brief Returns flux integrated in sky region
  717 + *
  718 + * @param[in] region Sky region.
  719 + * @param[in] srcEng Energy.
  720 + * @param[in] srcTime Time.
  721 + * @return Flux (adimensional or ph/cm2/s).
  722 + *
  723 + * Returns flux within a sky region.
  724 + ***************************************************************************/
  725 +double GModelSpatialComposite::flux(const GSkyRegion& region,
  726 + const GEnergy& srcEng,
  727 + const GTime& srcTime) const
  728 +{
  729 + // Initialise flux
  730 + double flux = 0.0;
  731 +
  732 + // Sum over all components
  733 + for (int i = 0; i < m_components.size(); ++i) {
  734 + flux += m_components[i]->flux(region, srcEng, srcTime) *
  735 + m_scales[i]->value();
  736 + }
  737 +
  738 + // Normalise sum by sum of scales
  739 + if (sum_of_scales() > 0) {
  740 + flux /= sum_of_scales();
  741 + }
  742 +
  743 + // Return flux
  744 + return flux;
  745 +}
  746 +
  747 +
  748 +/***********************************************************************//**
716 749 * @brief Print composite spatial model information
717 750 *
718 751 * @param[in] chatter Chattiness.
... ...
src/model/GModelSpatialDiffuseConst.cpp
  View file @61df648
1 1 /***************************************************************************
2 2 * GModelSpatialDiffuseConst.cpp - Spatial isotropic model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2021 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -351,7 +351,7 @@ GSkyDir GModelSpatialDiffuseConst::mc(const GEnergy&amp; energy,
351 351  
352 352 // Rotate sky direction by offset angles (phi, theta)
353 353 dir.rotate_deg(phi, theta);
354   -
  354 +
355 355 // Return sky direction
356 356 return dir;
357 357 }
... ... @@ -473,9 +473,42 @@ void GModelSpatialDiffuseConst::write(GXmlElement&amp; xml) const
473 473  
474 474  
475 475 /***********************************************************************//**
  476 + * @brief Returns isotropic flux integrated in sky region
  477 + *
  478 + * @param[in] region Sky region.
  479 + * @param[in] srcEng Energy.
  480 + * @param[in] srcTime Time.
  481 + * @return Flux (adimensional or ph/cm2/s).
  482 + *
  483 + * Returns isotropic flux within a sky region. The flux \f$F\f$ is computed
  484 + * using
  485 + *
  486 + * \f[F = \frac{{\tt value}}{4 \pi} \times \Omega\f]
  487 + *
  488 + * where
  489 + * - \f${\tt value}\f$ is the normalisation factor returned by the value()
  490 + * method, and
  491 + * - \f$\Omega\f$ is the solid angle of the sky region.
  492 + ***************************************************************************/
  493 +double GModelSpatialDiffuseConst::flux(const GSkyRegion& region,
  494 + const GEnergy& srcEng,
  495 + const GTime& srcTime) const
  496 +{
  497 + // Set normalization constant
  498 + const double norm = 1.0 / gammalib::fourpi;
  499 +
  500 + // Compute flux in sky region
  501 + double flux = norm * m_value.value() * region.solidangle();
  502 +
  503 + // Return flux
  504 + return flux;
  505 +}
  506 +
  507 +
  508 +/***********************************************************************//**
476 509 * @brief Print isotropic source model information
477 510 *
478   - * @param[in] chatter Chattiness (defaults to NORMAL).
  511 + * @param[in] chatter Chattiness.
479 512 * @return String containing model information.
480 513 ***************************************************************************/
481 514 std::string GModelSpatialDiffuseConst::print(const GChatter& chatter) const
... ...
src/model/GModelSpatialDiffuseCube.cpp
  View file @61df648
1 1 /***************************************************************************
2 2 * GModelSpatialDiffuseCube.cpp - Spatial map cube model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2021 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -1040,6 +1040,83 @@ void GModelSpatialDiffuseCube::write(GFits&amp; fits) const
1040 1040  
1041 1041  
1042 1042 /***********************************************************************//**
  1043 + * @brief Returns diffuse cube flux integrated in sky region
  1044 + *
  1045 + * @param[in] region Sky region.
  1046 + * @param[in] srcEng Energy.
  1047 + * @param[in] srcTime Time.
  1048 + * @return Flux (adimensional or ph/cm2/s).
  1049 + *
  1050 + * Returns diffuse cube flux within a sky region. The flux is computed by
  1051 + * log-log interpolation of the flux of the sky maps that encompass the
  1052 + * specified @p srcEng energy. Specifically
  1053 + *
  1054 + * \f[
  1055 + * F = {\tt value} \times \left \{
  1056 + * \begin{array}{l l}
  1057 + * 0 & \mbox{if $F_l=0$ and $F_r=0$} \\
  1058 + * F_l & \mbox{if $F_l\ne0$ and $F_r=0$} \\
  1059 + * F_r & \mbox{if $F_l=0$ and $F_r\ne0$} \\
  1060 + * \exp \left( w_l \log F_l + w_r \log F_r \right) & \mbox{else}
  1061 + * \end{array}
  1062 + * \right .
  1063 + * \f]
  1064 + *
  1065 + * where
  1066 + * - \f$F\f$ is the diffuse cube flux
  1067 + * - \f${\tt value}\f$ is the cube normalisation that is returned by the
  1068 + * value() method,
  1069 + * - \f$w_l\f$ and \f$w_r\f$ are the weighting factors from the logarithmic
  1070 + * interpolation in @p srcEng, with \f$w_l+w_r=1\f$, and
  1071 + * - \f$F_l\f$ and \f$F_r\f$ are the fluxes in the sky maps that encompass
  1072 + * the @p srcEng energy.
  1073 + *
  1074 + * This method fetches the diffuse cube if it was not yet fetched before.
  1075 + ***************************************************************************/
  1076 +double GModelSpatialDiffuseCube::flux(const GSkyRegion& region,
  1077 + const GEnergy& srcEng,
  1078 + const GTime& srcTime) const
  1079 +{
  1080 + // Fetch cube
  1081 + fetch_cube();
  1082 +
  1083 + // Initialise flux value
  1084 + double flux = 0.0;
  1085 +
  1086 + // Continue only if there is energy information for the map cube
  1087 + if (m_logE.size() > 0) {
  1088 +
  1089 + // Set energy for interpolation in log-energy
  1090 + m_logE.set_value(srcEng.log10MeV());
  1091 +
  1092 + // Compute map cube flux for the left and right map
  1093 + double left_flux = m_cube.flux(region, m_logE.inx_left());
  1094 + double right_flux = m_cube.flux(region, m_logE.inx_right());
  1095 +
  1096 + // Perform log-log interpolation
  1097 + if (left_flux > 0.0 && right_flux > 0.0) {
  1098 + double log_flux = m_logE.wgt_left() * std::log(left_flux) +
  1099 + m_logE.wgt_right() * std::log(right_flux);
  1100 + flux = std::exp(log_flux);
  1101 + }
  1102 + else if (left_flux > 0.0) {
  1103 + flux = left_flux;
  1104 + }
  1105 + else if (right_flux > 0.0) {
  1106 + flux = right_flux;
  1107 + }
  1108 +
  1109 + // Multiply by cube normalisation
  1110 + flux *= m_value.value();
  1111 +
  1112 + } // endif: energy information was available
  1113 +
  1114 + // Return flux
  1115 + return flux;
  1116 +}
  1117 +
  1118 +
  1119 +/***********************************************************************//**
1043 1120 * @brief Print map cube information
1044 1121 *
1045 1122 * @param[in] chatter Chattiness.
... ...
src/model/GModelSpatialDiffuseMap.cpp
  View file @61df648
1 1 /***************************************************************************
2 2 * GModelSpatialDiffuseMap.cpp - Spatial map model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2021 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -701,7 +701,7 @@ void GModelSpatialDiffuseMap::mc_cone(const GSkyRegionCircle&amp; cone) const
701 701 }
702 702 sum_map += flux; // total flux
703 703 }
704   -
  704 +
705 705 // Update maximum intensity
706 706 if (distance <= m_mc_cone.radius()) {
707 707 if (intensity > m_mc_max) {
... ... @@ -746,6 +746,30 @@ void GModelSpatialDiffuseMap::mc_cone(const GSkyRegionCircle&amp; cone) const
746 746  
747 747  
748 748 /***********************************************************************//**
  749 + * @brief Returns diffuse map flux integrated in sky region
  750 + *
  751 + * @param[in] region Sky region.
  752 + * @param[in] srcEng Energy.
  753 + * @param[in] srcTime Time.
  754 + * @return Flux (adimensional or ph/cm2/s).
  755 + *
  756 + * Returns diffuse map flux within a sky region. The flux is defined as the
  757 + * sum of the flux in all diffuse map pixels that are contained within the
  758 + * sky region multiplied by the diffuse map normalisation factor.
  759 + ***************************************************************************/
  760 +double GModelSpatialDiffuseMap::flux(const GSkyRegion& region,
  761 + const GEnergy& srcEng,
  762 + const GTime& srcTime) const
  763 +{
  764 + // Compute flux in sky region
  765 + double flux = m_map.flux(region) * m_value.value();
  766 +
  767 + // Return flux
  768 + return flux;
  769 +}
  770 +
  771 +
  772 +/***********************************************************************//**
749 773 * @brief Print map information
750 774 *
751 775 * @param[in] chatter Chattiness.
... ...
src/model/GModelSpatialPointSource.cpp
  View file @61df648
1 1 /***************************************************************************
2 2 * GModelSpatialPointSource.cpp - Spatial point source model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2021 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -462,6 +462,35 @@ void GModelSpatialPointSource::write(GXmlElement&amp; xml) const
462 462  
463 463  
464 464 /***********************************************************************//**
  465 + * @brief Returns model flux integrated in sky region
  466 + *
  467 + * @param[in] region Sky region.
  468 + * @param[in] srcEng Energy.
  469 + * @param[in] srcTime Time.
  470 + * @return Flux (adimensional or ph/cm2/s).
  471 + *
  472 + * Returns point source flux within a sky region. If the point source is
  473 + * contained within the sky region the flux will be 1, otherwise the flux
  474 + * will be 0.
  475 + ***************************************************************************/
  476 +double GModelSpatialPointSource::flux(const GSkyRegion& region,
  477 + const GEnergy& srcEng,
  478 + const GTime& srcTime) const
  479 +{
  480 + // Initialise flux
  481 + double flux = 0.0;
  482 +
  483 + // If point source region overlaps with sky region then set flux to 1
  484 + if (this->region()->overlaps(region)) {
  485 + flux = 1.0;
  486 + }
  487 +
  488 + // Return flux
  489 + return flux;
  490 +}
  491 +
  492 +
  493 +/***********************************************************************//**
465 494 * @brief Print point source information
466 495 *
467 496 * @param[in] chatter Chattiness.
... ...