Josh Cardenzana / gammalib

Commit 68de274ece7318f41c25c5a5a65ad4a12da17098

Authored by Jürgen Knödlseder
1 parent 1315a996

Add GModelSpatial::flux() method (#2961)

Inline Side-by-side
Showing 14 changed files with 259 additions and 306 deletions
ChangeLog
  View file @68de274
1   -2020-01-29
  1 +2020-01-31
2 2  
3 3 * Version 1.7.0 released
4 4 ========================
... ... @@ -11,6 +11,7 @@
11 11 Add GCOMObservation::drm(GModels&) method (#2973)
12 12 Implement methods to remove models from a response cache (#2973)
13 13 Improve COMPTEL classes (#2973)
  14 + Add GModelSpatial::flux() method (#2961)
14 15 Change order of Python module import (#2961)
15 16 Add GGti::overlap() method
16 17 Correctly write model attributes in GCOMModelDRBFitting class
... ...
1 1 New Features and Important Changes in GammaLib 1.7.0
2 2  
3   -29 January 2020
  3 +31 January 2020
4 4  
5 5  
6 6 1. Introduction
... ... @@ -35,6 +35,7 @@ The following methods have been added:
35 35 - GSkyDir::sun() - Set sky position of Sun
36 36 - GSkyDir::moon() - Set sky position of Moon
37 37 - GSkyMap::extract() - Extract section of sky map
  38 +- GModelSpatial::flux() - Extract flux from sky region
38 39 - GCTAObservation::remove_response_cache(string&) - Remove response cache for model
39 40 - GCTAResponse::remove_response_cache(string&) - Remove response cache for model
40 41 - GCTAResponseCache::remove(string&) - Remove response cache for model
... ... @@ -86,7 +87,9 @@ None
86 87  
87 88 8. Model module
88 89 ---------------
89   -None
  90 +The method GModelSpatial::flux() was added that extracts the flux by integrating
  91 +over a specific sky region. For the moment only circular sky regions are
  92 +implemented (#2961).
90 93  
91 94  
92 95 9. Numerics module
... ...
doc/source/admin/release_history/1.7.rst
  View file @68de274
... ... @@ -36,6 +36,8 @@ Improvements
36 36 Implement methods to remove models from a response cache.
37 37 Add ``GCOMObservation::drm(GModels&)`` method.
38 38 Add ``GCOMDri::cone_content()`` method.
  39 +* [`2961 <https://cta-redmine.irap.omp.eu/issues/2961>`_] -
  40 + Add ``GModelSpatial::flux()`` method.
39 41 * Add ``GGti::overlap()`` method.
40 42 * [`2909 <https://cta-redmine.irap.omp.eu/issues/2909>`_] -
41 43 Add ``GTime.datetime()`` method and add optional precision to ``GTime.utc()`` method.
... ... @@ -48,6 +50,9 @@ Improvements
48 50 New features
49 51 ------------
50 52  
  53 +* [`2961 <https://cta-redmine.irap.omp.eu/issues/2961>`_] -
  54 + Change order of Python module import.
  55 +
51 56  
52 57 Clean up
53 58 --------
... ...
include/GModelSpatial.hpp
  View file @68de274
1 1 /***************************************************************************
2 2 * GModelSpatial.hpp - Spatial model abstract base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2018 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -105,7 +105,6 @@ public:
105 105 virtual void write(GXmlElement& xml) const = 0;
106 106 virtual GSkyRegion* region(void) const = 0;
107 107 virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
108   -
109 108  
110 109 // Methods
111 110 GModelPar& at(const int& index);
... ... @@ -114,9 +113,9 @@ public:
114 113 bool has_free_pars(void) const;
115 114 int size(void) const;
116 115 void autoscale(void);
117   - double flux(GSkyRegion* reg,
118   - const GEnergy& srcEng = GEnergy(),
119   - const GTime& srcTime = GTime()) const;
  116 + double flux(const GSkyRegion* reg,
  117 + const GEnergy& srcEng = GEnergy(),
  118 + const GTime& srcTime = GTime()) const;
120 119  
121 120 protected:
122 121 // Protected methods
... ... @@ -124,101 +123,61 @@ protected:
124 123 void copy_members(const GModelSpatial& model);
125 124 void free_members(void);
126 125  
127   - /***********************************************************************//**
128   - * @class circle_int_kern_rho
129   - *
130   - * @brief Kernel for circular sky region radial integration
131   - *
132   - * This class implements the integration kernel \f$K(\rho)\f$ for the
133   - * integration
134   - *
135   - * \f[
136   - * \int_{\rho_{\rm min}}^{\rho_{\rm max}} K(\rho | E, t) d\rho
137   - * \f]
138   - *
139   - * of a spatial model over a circular region. The eval() method computes
140   - *
141   - * \f[
142   - * K(\rho | E, t) = \sin \rho \times
143   - * \int_{\omega} M(\rho, \omega | E, t) d\omega
144   - * \f]
145   - *
146   - * where
147   - * \f$M(\rho, \omega | E, t)\f$ is the spatial model,
148   - * \f$\rho\f$ is the distance from the region centre, and
149   - * \f$\omega\f$ is the position angle with respect to the connecting line
150   - * between the region centre and the direction on the sky.
151   - ***************************************************************************/
152   - class circle_int_kern_rho : public GFunction {
153   - public:
154   - circle_int_kern_rho(const GModelSpatial* model,
155   - const GSkyRegionCircle& reg,
156   - const GEnergy& srcEng,
157   - const GTime& srcTime,
158   - const double distance,
159   - const double cosdist,
160   - const double sindist,
161   - const double modrad,
162   - const double cosmodrad) :
163   - m_model(model),
164   - m_reg(reg),
165   - m_srcEng(srcEng),
166   - m_srcTime(srcTime),
167   - m_dist(distance),
168   - m_cosdist(cosdist),
169   - m_sindist(sindist),
170   - m_modrad(modrad),
171   - m_cosmodrad(cosmodrad) { }
172   - double eval(const double& rho);
173   - public:
174   - const GModelSpatial* m_model; //!< Model
175   - GSkyRegionCircle m_reg; //!< Region
176   - GEnergy m_srcEng; //!< Photon energy
177   - GTime m_srcTime; //!< Photon time
178   - double m_dist; //!< Distance model-region (rad)
179   - double m_cosdist; //!< Cos of distance model-region
180   - double m_sindist; //!< Sin of distance model-region
181   - double m_modrad; //!< Model radius (rad)
182   - double m_cosmodrad; //!< Cos of model radius
183   - };
184   -
185   - /***********************************************************************//**
186   - * @class circle_int_kern_omega
187   - *
188   - * @brief Kernel for circular sky region azimuth angle integration
189   - *
190   - * This class implements the computation of
191   - *
192   - * \f[
193   - * K(\omega | \rho, E, t) = M(\omega | \rho)
194   - * \f]
195   - *
196   - * where
197   - * \fM(\omega | \rho)\f$ is the spatial model,
198   - * \f$\rho\f$ is the distance from the region centre, and
199   - * \f$\omega\f$ is the position angle with respect to the connecting line
200   - * between the region centre and the direction on the sky.
201   - ***************************************************************************/
202   - class circle_int_kern_omega : public GFunction {
203   - public:
204   - circle_int_kern_omega(const GModelSpatial* model,
205   - const GSkyRegionCircle& reg,
206   - const double& rho,
207   - const GEnergy& srcEng,
208   - const GTime& srcTime) :
209   - m_model(model),
210   - m_reg(reg),
211   - m_rho(rho),
212   - m_srcEng(srcEng),
213   - m_srcTime(srcTime) { }
214   - double eval(const double& omega);
215   - public:
216   - const GModelSpatial* m_model; //!< Model
217   - GSkyRegionCircle m_reg; //!< Region
218   - double m_rho; //!< Offset from center of the region
219   - GEnergy m_srcEng; //!< Photon energy
220   - GTime m_srcTime; //!< Photon time
221   - };
  126 + // Kernel for circular sky region radial integration
  127 + class circle_int_kern_rho : public GFunction {
  128 + public:
  129 + circle_int_kern_rho(const GModelSpatial* model,
  130 + const GSkyRegionCircle* reg,
  131 + const GEnergy& srcEng,
  132 + const GTime& srcTime,
  133 + const double distance,
  134 + const double cosdist,
  135 + const double sindist,
  136 + const double modrad,
  137 + const double cosmodrad) :
  138 + m_model(model),
  139 + m_reg(reg),
  140 + m_srcEng(srcEng),
  141 + m_srcTime(srcTime),
  142 + m_dist(distance),
  143 + m_cosdist(cosdist),
  144 + m_sindist(sindist),
  145 + m_modrad(modrad),
  146 + m_cosmodrad(cosmodrad) { }
  147 + double eval(const double& rho);
  148 + public:
  149 + const GModelSpatial* m_model; //!< Spatial model
  150 + const GSkyRegionCircle* m_reg; //!< Interation region
  151 + GEnergy m_srcEng; //!< Photon energy
  152 + GTime m_srcTime; //!< Photon time
  153 + double m_dist; //!< Distance model-region (rad)
  154 + double m_cosdist; //!< Cos of distance model-region
  155 + double m_sindist; //!< Sin of distance model-region
  156 + double m_modrad; //!< Model radius (rad)
  157 + double m_cosmodrad; //!< Cos of model radius
  158 + };
  159 +
  160 + // Kernel for circular sky region azimuth angle integration
  161 + class circle_int_kern_omega : public GFunction {
  162 + public:
  163 + circle_int_kern_omega(const GModelSpatial* model,
  164 + const GSkyRegionCircle* reg,
  165 + const double& rho,
  166 + const GEnergy& srcEng,
  167 + const GTime& srcTime) :
  168 + m_model(model),
  169 + m_reg(reg),
  170 + m_rho(rho),
  171 + m_srcEng(srcEng),
  172 + m_srcTime(srcTime) { }
  173 + double eval(const double& omega);
  174 + public:
  175 + const GModelSpatial* m_model; //!< Spatial model
  176 + const GSkyRegionCircle* m_reg; //!< Interation region
  177 + double m_rho; //!< Offset from center of the region
  178 + GEnergy m_srcEng; //!< Photon energy
  179 + GTime m_srcTime; //!< Photon time
  180 + };
222 181  
223 182 // Proteced members
224 183 std::vector<GModelPar*> m_pars; //!< Parameter pointers
... ...
include/GTools.hpp
  View file @68de274
1 1 /***************************************************************************
2 2 * GTools.hpp - GammaLib tools *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2008-2019 by Juergen Knoedlseder *
  4 + * copyright (C) 2008-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ...
inst/cta/src/GCTAModelRadialAcceptance.cpp
  View file @68de274
1 1 /***************************************************************************
2 2 * GCTAModelRadialAcceptance.cpp - Radial acceptance model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2018 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -1052,17 +1052,17 @@ double GCTAModelRadialAcceptance::roi_kern::eval(const double&amp; offset)
1052 1052  
1053 1053 // Initialise Npred value
1054 1054 double value = 0.0;
1055   -
  1055 +
1056 1056 // Continue only if offset > 0
1057 1057 if (offset > 0.0) {
1058 1058  
1059 1059 // Get arclength for given radius in radians.
1060 1060 double phi = gammalib::roi_arclength(offset,
1061   - m_dist,
1062   - m_cosdist,
1063   - m_sindist,
1064   - m_roi,
1065   - m_cosroi);
  1061 + m_dist,
  1062 + m_cosdist,
  1063 + m_sindist,
  1064 + m_roi,
  1065 + m_cosroi);
1066 1066  
1067 1067 // Get kernel value if phi > 0
1068 1068 if (phi > 0.0) {
... ...
inst/cta/src/GCTAResponse_helpers.cpp
  View file @68de274
1 1 /***************************************************************************
2 2 * GCTAResponse_helpers.cpp - CTA response helper classes *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2018 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -633,11 +633,11 @@ double cta_nroi_radial_kern_rho::eval(const double&amp; rho)
633 633  
634 634 // Compute half length of arc that lies within ROI+PSF radius (radians)
635 635 double domega = 0.5 * gammalib::roi_arclength(rho,
636   - m_dist,
637   - m_cos_dist,
638   - m_sin_dist,
639   - m_radius,
640   - m_cos_radius);
  636 + m_dist,
  637 + m_cos_dist,
  638 + m_sin_dist,
  639 + m_radius,
  640 + m_cos_radius);
641 641  
642 642 // Continue only if arc length is positive
643 643 if (domega > 0.0) {
... ... @@ -810,11 +810,11 @@ double cta_irf_elliptical_kern_rho::eval(const double&amp; rho)
810 810 // radius rho that intersects with the point spread function, defined
811 811 // as a circle with maximum radius m_delta_max
812 812 double domega = 0.5 * gammalib::roi_arclength(rho,
813   - m_rho_obs,
814   - m_cos_rho_obs,
815   - m_sin_rho_obs,
816   - m_delta_max,
817   - m_cos_delta_max);
  813 + m_rho_obs,
  814 + m_cos_rho_obs,
  815 + m_sin_rho_obs,
  816 + m_delta_max,
  817 + m_cos_delta_max);
818 818  
819 819 // Continue only if arc length is positive
820 820 if (domega > 0.0) {
... ... @@ -1116,11 +1116,11 @@ double cta_nroi_elliptical_kern_rho::eval(const double&amp; rho)
1116 1116  
1117 1117 // Compute half length of arc that lies within ROI+PSF radius (radians)
1118 1118 double domega = 0.5 * gammalib::roi_arclength(rho,
1119   - m_rho_roi,
1120   - m_cos_rho_roi,
1121   - m_sin_rho_roi,
1122   - m_radius_roi,
1123   - m_cos_radius_roi);
  1119 + m_rho_roi,
  1120 + m_cos_rho_roi,
  1121 + m_sin_rho_roi,
  1122 + m_radius_roi,
  1123 + m_cos_radius_roi);
1124 1124  
1125 1125 // Continue only if arc length is positive
1126 1126 if (domega > 0.0) {
... ... @@ -1748,7 +1748,7 @@ double cta_psf_radial_kern_rho::eval(const double&amp; rho)
1748 1748 m_rho_obs,
1749 1749 m_cos_rho_obs,
1750 1750 m_sin_rho_obs,
1751   - m_delta_max,
  1751 + m_delta_max,
1752 1752 m_cos_delta_max);
1753 1753  
1754 1754 // Continue only if arc length is positive
... ... @@ -1928,11 +1928,11 @@ double cta_psf_radial_kern_delta::eval(const double&amp; delta)
1928 1928 // radius delta that intersects with the model, defined as a circle
1929 1929 // with maximum radius m_theta_max
1930 1930 double dphi = 0.5 * gammalib::roi_arclength(delta,
1931   - m_delta_mod,
1932   - m_cos_delta_mod,
1933   - m_sin_delta_mod,
1934   - m_theta_max,
1935   - m_cos_theta_max);
  1931 + m_delta_mod,
  1932 + m_cos_delta_mod,
  1933 + m_sin_delta_mod,
  1934 + m_theta_max,
  1935 + m_cos_theta_max);
1936 1936  
1937 1937 // Continue only if arc length is positive
1938 1938 if (dphi > 0.0) {
... ... @@ -2085,11 +2085,11 @@ double cta_psf_elliptical_kern_rho::eval(const double&amp; rho)
2085 2085 // radius rho that intersects with the point spread function, defined
2086 2086 // as a circle with maximum radius m_delta_max
2087 2087 double domega = 0.5 * gammalib::roi_arclength(rho,
2088   - m_rho_obs,
2089   - m_cos_rho_obs,
2090   - m_sin_rho_obs,
2091   - m_delta_max,
2092   - m_cos_delta_max);
  2088 + m_rho_obs,
  2089 + m_cos_rho_obs,
  2090 + m_sin_rho_obs,
  2091 + m_delta_max,
  2092 + m_cos_delta_max);
2093 2093  
2094 2094 // Continue only if arc length is positive
2095 2095 if (domega > 0.0) {
... ...
inst/cta/src/GCTASupport.cpp
  View file @68de274
... ... @@ -67,7 +67,6 @@
67 67 /* __ Coding definitions _________________________________________________ */
68 68  
69 69 /* __ Debug definitions __________________________________________________ */
70   -#define G_CHECK_FOR_NAN 0
71 70  
72 71  
73 72 /***********************************************************************//**
... ...
pyext/GModelSpatial.i
  View file @68de274
1 1 /***************************************************************************
2 2 * GModelSpatial.i - Spatial model abstract base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2018 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -64,9 +64,9 @@ public:
64 64 bool has_free_pars(void) const;
65 65 int size(void) const;
66 66 void autoscale(void);
67   - double flux(GSkyRegion* reg,
68   - const GEnergy& srcEng = GEnergy(),
69   - const GTime& srcTime = GTime()) const;
  67 + double flux(const GSkyRegion* reg,
  68 + const GEnergy& srcEng = GEnergy(),
  69 + const GTime& srcTime = GTime()) const;
70 70 };
71 71  
72 72  
... ...
pyext/GTools.i
  View file @68de274
1 1 /***************************************************************************
2 2 * GTools.i - GammaLib tools *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2008-2018 by Juergen Knoedlseder *
  4 + * copyright (C) 2008-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ...
src/model/GModelSpatial.cpp
  View file @68de274
1 1 /***************************************************************************
2 2 * GModelSpatial.cpp - Abstract spatial model base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2018 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -34,9 +34,8 @@
34 34 /* __ Method name definitions ____________________________________________ */
35 35 #define G_ACCESS "GModelSpatial::operator[](std::string&)"
36 36 #define G_AT "GModelPar& GModelSpatial::at(int&)"
37   -#define G_FLUX "double GModelSpatial::flux(const GSkyRegionCircle& reg, "\
38   - "const GEnergy& srcEng = GEnergy(),"\
39   - " const GTime& srcTime = GTime())"
  37 +#define G_FLUX "GModelSpatial::flux(GSkyRegionCircle&, GEnergy&, GTime&)"
  38 +
40 39 /* __ Constants __________________________________________________________ */
41 40 const double g_kludge_radius = 1.0e-12; //!< Tiny angle (radians)
42 41  
... ... @@ -312,13 +311,17 @@ void GModelSpatial::autoscale(void)
312 311 return;
313 312 }
314 313  
  314 +
315 315 /***********************************************************************//**
316 316 * @brief Returns model flux integrated in circular sky region
317 317 *
318 318 * @param[in] reg Sky region.
319 319 * @param[in] srcEng Energy.
320 320 * @param[in] srcTime Time.
321   - * @return flux (adimensional or or ph/cm2/s).
  321 + * @return Flux (adimensional or ph/cm2/s).
  322 + *
  323 + * @exception GException::feature_not_implemented
  324 + * Regions are not circular.
322 325 *
323 326 * Computes
324 327 *
... ... @@ -333,96 +336,76 @@ void GModelSpatial::autoscale(void)
333 336 * \f$\omega\f$ is the position angle with respect to the connecting line
334 337 * between the region centre and the direction on the sky.
335 338 ***************************************************************************/
336   -double GModelSpatial::flux(GSkyRegion* reg,
337   - const GEnergy& srcEng,
338   - const GTime& srcTime) const
339   -{
  339 +double GModelSpatial::flux(const GSkyRegion* reg,
  340 + const GEnergy& srcEng,
  341 + const GTime& srcTime) const
  342 +{
340 343 // Initialise flux
341 344 double flux = 0.0;
342 345  
343   - // Check if region overlaps with model
344   - GSkyRegion* model_reg = this -> region();
345   - if (model_reg -> overlaps(*reg)) {
  346 + // Continue only if region overlaps with model
  347 + const GSkyRegion* model_reg = this->region();
  348 + if (model_reg->overlaps(*reg)) {
  349 +
  350 + // Thrown an exception if region is not a sky circle
  351 + const GSkyRegionCircle* reg_circle =
  352 + dynamic_cast<const GSkyRegionCircle*>(reg);
  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 + }
346 358  
347   - // Extract the region circle
348   - // This only works for a region that results a GSkyRegionCircle
349   - GSkyRegionCircle* reg_circle = dynamic_cast<GSkyRegionCircle*>(reg);
350   -
351   - // If region circle is valid then
352   - // continue
353   - if (reg_circle != NULL) {
354   -
355   - // Extract the model region circle
356   - // This only works for a model that results a GSkyRegionCircle
357   - GSkyRegionCircle* model_circle = dynamic_cast<GSkyRegionCircle*>(model_reg);
358   -
359   - // If model region circle is valid then
360   - // estimate flux in the overlap region between model circle
361   - // and region of interest
362   - if (model_circle != NULL) {
363   -
364   - // Model centre and radius (in rad)
365   - GSkyDir model_centre = model_circle -> centre();
366   - double model_radius = model_circle -> radius();
367   - model_radius *= gammalib::deg2rad;
368   -
369   - // Distance between region and model centres
370   - double distance = model_centre.dist(reg_circle -> centre());
371   -
372   - // Minimum radial integration boundary (in radians)
373   - // Distance between centers minus radius of model region
374   - double rho_min = distance - model_radius;
375   -
376   - // If rho_min is < 0 set it to 0
377   - if (rho_min < 0.) {
378   - rho_min = 0;
379   - }
380   -
381   - // Maximum radial integration boundary corresponds to region radius (in radians)
382   - double rho_max = reg_circle -> radius() * gammalib::deg2rad;
383   -
384   - // Pre-compute some quantities for arclength
385   - double cosdist = std::cos(distance);
386   - double sindist = std::sin(distance);
387   - double cosmodrad = std::cos(model_radius);
388   -
389   - // Setup integration kernel
390   - GModelSpatial::circle_int_kern_rho integrand(this,
391   - *reg_circle,
392   - srcEng,
393   - srcTime,
394   - distance,
395   - cosdist,
396   - sindist,
397   - model_radius,
398   - cosmodrad);
399   - GIntegral integral(&integrand);
400   -
401   - // Perform integration
402   - flux = integral.romberg(rho_min, rho_max);
403   -
404   - }
405   -
406   - // Otherwise, if model region circle is not valid throw exception
407   - else {
408   -
409   - std::string msg = "Flux can only be computed for spatial model"
410   - " with region defined as circle.";
411   - throw GException::feature_not_implemented(G_FLUX,msg);
412   -
413   - }// endif model has a valid region circle
414   -
415   - }
416   -
417   - else {
418   - std::string msg = "Flux can only be computed for a"
419   - " circular region.";
420   - throw GException::feature_not_implemented(G_FLUX,msg);
421   -
422   - }// endif region is a circle
423   -
424   - }// endif model overlaps with region
  359 + // Thrown an exception if model region is not a sky circle
  360 + const GSkyRegionCircle* model_circle =
  361 + dynamic_cast<const GSkyRegionCircle*>(model_reg);
  362 + if (model_circle == NULL) {
  363 + std::string msg = "Flux can only be computed for spatial model "
  364 + "with region defined as circle.";
  365 + throw GException::feature_not_implemented(G_FLUX,msg);
  366 + }
  367 +
  368 + // Model centre and radius in radians
  369 + GSkyDir model_centre = model_circle->centre();
  370 + double model_radius = model_circle->radius() * gammalib::deg2rad;
  371 +
  372 + // Distance between region and model centres in radians
  373 + double distance = model_centre.dist(reg_circle->centre());
425 374  
  375 + // Take the distance between centers minus radius of model region as
  376 + // the minimum radial integration boundary (in radians)
  377 + double rho_min = distance - model_radius;
  378 +
  379 + // Make sure that rho_min does not become negative
  380 + if (rho_min < 0.0) {
  381 + rho_min = 0.0;
  382 + }
  383 +
  384 + // Take the region radius as the maximum radial integration boundary
  385 + // (in radians)
  386 + double rho_max = reg_circle->radius() * gammalib::deg2rad;
  387 +
  388 + // Pre-compute some quantities for arclength
  389 + double cosdist = std::cos(distance);
  390 + double sindist = std::sin(distance);
  391 + double cosmodrad = std::cos(model_radius);
  392 +
  393 + // Setup integration kernel
  394 + GModelSpatial::circle_int_kern_rho integrand(this,
  395 + reg_circle,
  396 + srcEng,
  397 + srcTime,
  398 + distance,
  399 + cosdist,
  400 + sindist,
  401 + model_radius,
  402 + cosmodrad);
  403 + GIntegral integral(&integrand);
  404 +
  405 + // Perform integration
  406 + flux = integral.romberg(rho_min, rho_max);
  407 +
  408 + } // endif: model overlaps with region
426 409  
427 410 // Return
428 411 return flux;
... ... @@ -473,6 +456,7 @@ void GModelSpatial::free_members(void)
473 456 return;
474 457 }
475 458  
  459 +
476 460 /***********************************************************************//**
477 461 * @brief Kernel for circular sky region radial integration
478 462 *
... ... @@ -501,62 +485,62 @@ double GModelSpatial::circle_int_kern_rho::eval(const double&amp; rho)
501 485 // Initialise model value
502 486 double flux = 0.0;
503 487  
504   - // Continue only if rho is positive
  488 + // Continue only if rho is positive
505 489 if (rho > 0.0) {
506 490  
507   - // Compute half length of the arc (in radians) from a circle with
508   - // radius rho that intersects with model circle
509   - double domega = 0.5 * gammalib::roi_arclength(rho,
510   - m_dist,
511   - m_cosdist,
512   - m_sindist,
513   - m_modrad,
514   - m_cosmodrad);
515   -
516   - // Continue only if arc length is positive
517   - if (domega > 0.0) {
518   -
519   - // Reduce rho by an infinite amount to avoid rounding errors
520   - // at the boundary of a sharp edged model
521   - double rho_kludge = rho - g_kludge_radius;
522   - if (rho_kludge < 0.0) {
523   - rho_kludge = 0.0;
524   - }
  491 + // Compute half length of the arc (in radians) from a circle with
  492 + // radius rho that intersects with model circle
  493 + double domega = 0.5 * gammalib::roi_arclength(rho,
  494 + m_dist,
  495 + m_cosdist,
  496 + m_sindist,
  497 + m_modrad,
  498 + m_cosmodrad);
  499 +
  500 + // Continue only if arc length is positive
  501 + if (domega > 0.0) {
  502 +
  503 + // Reduce rho by an infinite amount to avoid rounding errors
  504 + // at the boundary of a sharp edged model
  505 + double rho_kludge = rho - g_kludge_radius;
  506 + if (rho_kludge < 0.0) {
  507 + rho_kludge = 0.0;
  508 + }
525 509  
526   - // Compute omega integration range
527   - double omega_min = -domega;
528   - double omega_max = +domega;
  510 + // Compute omega integration range
  511 + double omega_min = -domega;
  512 + double omega_max = +domega;
529 513  
530   - // Setup integration kernel for azimuth integration
531   - GModelSpatial::circle_int_kern_omega integrand(m_model,
532   - m_reg,
533   - rho_kludge,
534   - m_srcEng,
535   - m_srcTime);
  514 + // Setup integration kernel for azimuth integration
  515 + GModelSpatial::circle_int_kern_omega integrand(m_model,
  516 + m_reg,
  517 + rho_kludge,
  518 + m_srcEng,
  519 + m_srcTime);
536 520  
537   - // Setup integrator
538   - GIntegral integral(&integrand);
  521 + // Setup integrator
  522 + GIntegral integral(&integrand);
539 523  
540   - // Compute sine term for radial integration
541   - double sin_rho = std::sin(rho);
  524 + // Compute sine term for radial integration
  525 + double sin_rho = std::sin(rho);
542 526  
543   - // Integrate over omega
544   - flux = integral.romberg(omega_min, omega_max) * sin_rho;
  527 + // Integrate over omega
  528 + flux = integral.romberg(omega_min, omega_max) * sin_rho;
545 529  
546   - } //endif: domega was positive
  530 + } //endif: domega was positive
547 531  
548 532 } // endif: rho was positive
549   -
550   -
551   - // Return
  533 +
  534 + // Return
552 535 return flux;
553 536 }
554 537  
  538 +
555 539 /***********************************************************************//**
556 540 * @brief Kernel for circular sky region azimuth angle integration
557 541 *
558 542 * @param[in] omega Azimuth angle (radians).
559   -* Computes
  543 + * @return Integration kernel defined as
560 544 *
561 545 * \f[
562 546 * K(\omega | \rho, E, t) = M(\omega | \rho)
... ... @@ -575,19 +559,20 @@ double GModelSpatial::circle_int_kern_omega::eval(const double&amp; omega)
575 559  
576 560 // Precompute rho and omega in degrees
577 561 double omega_deg = omega * gammalib::rad2deg;
578   - double rho_deg = m_rho * gammalib::rad2deg;
  562 + double rho_deg = m_rho * gammalib::rad2deg;
579 563  
580 564 // Compute sky direction corresponding to given rho, omega
581   - // Deep copy of region centre
582   - GSkyDir dir = GSkyDir(m_reg.centre());
583   - // Rotate to otain requested sky direction
584   - dir.rotate_deg(omega_deg,rho_deg);
  565 + GSkyDir dir = GSkyDir(m_reg->centre());
  566 +
  567 + // Rotate to obtain requested sky direction
  568 + dir.rotate_deg(omega_deg, rho_deg);
  569 +
  570 + // Set photon for this sky direction
  571 + GPhoton photon = GPhoton(dir, m_srcEng, m_srcTime);
585 572  
586 573 // Evaluate model for this sky direction
587   - GPhoton photon = GPhoton(dir,m_srcEng,m_srcTime);
588 574 flux = m_model->eval(photon);
589 575  
590   - // Return
  576 + // Return
591 577 return flux;
592   -
593 578 }
... ...
src/support/GTools.cpp
  View file @68de274
1 1 /***************************************************************************
2 2 * GTools.cpp - GammaLib tools *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2008-2019 by Juergen Knoedlseder *
  4 + * copyright (C) 2008-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -55,6 +55,7 @@
55 55 #include "GXmlText.hpp"
56 56  
57 57 /* __ Compile options ____________________________________________________ */
  58 +//#define G_CHECK_FOR_NAN
58 59  
59 60 /* __ Function name definitions __________________________________________ */
60 61 #define G_XML2STRING "gammalib::xml2string(std::string&)"
... ... @@ -1789,12 +1790,12 @@ int gammalib::recv(int fd, char *buffer, int len, int flags, int timeout)
1789 1790 {
1790 1791 // Initialise error code with time out
1791 1792 int error = -2;
1792   -
  1793 +
1793 1794 // Initialize the set
1794 1795 fd_set readset;
1795 1796 FD_ZERO(&readset);
1796 1797 FD_SET(fd, &readset);
1797   -
  1798 +
1798 1799 // Initialize time out struct
1799 1800 struct timeval tv;
1800 1801 if (timeout >= 1000) {
... ... @@ -1805,7 +1806,7 @@ int gammalib::recv(int fd, char *buffer, int len, int flags, int timeout)
1805 1806 tv.tv_sec = 0;
1806 1807 tv.tv_usec = timeout*1000;
1807 1808 }
1808   -
  1809 +
1809 1810 // select()
1810 1811 int result = select(fd+1, &readset, NULL, NULL, &tv);
1811 1812  
... ... @@ -1822,23 +1823,24 @@ int gammalib::recv(int fd, char *buffer, int len, int flags, int timeout)
1822 1823 if ((iof = fcntl(fd, F_GETFL, 0)) != -1) {
1823 1824 fcntl(fd, F_SETFL, iof | O_NONBLOCK);
1824 1825 }
1825   -
  1826 +
1826 1827 // Receive data
1827 1828 result = ::recv(fd, buffer, len, flags);
1828   -
  1829 +
1829 1830 // Set as before
1830 1831 if (iof != -1) {
1831 1832 fcntl(fd, F_SETFL, iof);
1832 1833 }
1833   -
  1834 +
1834 1835 // Set error
1835 1836 error = result;
1836 1837 }
1837   -
  1838 +
1838 1839 // Return error
1839 1840 return error;
1840 1841 }
1841 1842  
  1843 +
1842 1844 /***********************************************************************//**
1843 1845 * @brief Returns length of circular arc within circular ROI
1844 1846 *
... ... @@ -1894,7 +1896,7 @@ double gammalib::roi_arclength(const double&amp; rad, const double&amp; dist,
1894 1896 double sinrad = std::sin(rad);
1895 1897 double cosang = (cosroi - cosdist*cosrad) / (sindist*sinrad);
1896 1898 arclength = 2.0 * gammalib::acos(cosang);
1897   - #if G_CHECK_FOR_NAN
  1899 + #if defined(G_CHECK_FOR_NAN)
1898 1900 if (gammalib::is_infinite(arclength) ||
1899 1901 gammalib::is_notanumber(arclength)) {
1900 1902 std::cout << "roi_arclength: NaN occured";
... ...
test/test_GModel.cpp
  View file @68de274
1 1 /***************************************************************************
2 2 * test_GModel.cpp - test GModel class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2019 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -3467,29 +3467,28 @@ void TestGModel::test_legacy_model_point_logparabola(void)
3467 3467 ***************************************************************************/
3468 3468 void TestGModel::test_flux(void)
3469 3469 {
3470   -
3471 3470 // Define model for test
3472 3471 GSkyDir centre1;
3473   - centre1.radec_deg(0.,0.5);
3474   - GModelSpatialRadialDisk model = GModelSpatialRadialDisk(centre1,0.25);
  3472 + centre1.radec_deg(0.0, 0.5);
  3473 + GModelSpatialRadialDisk model = GModelSpatialRadialDisk(centre1, 0.25);
3475 3474  
3476 3475 // Define entre of test ROI
3477 3476 GSkyDir centre2;
3478   - centre2.radec_deg(0.,0.);
  3477 + centre2.radec_deg(0.0, 0.0);
3479 3478  
3480 3479 // Test non overlapping region
3481   - GSkyRegionCircle circle1 = GSkyRegionCircle(centre2,0.1);
3482   - GSkyRegion* roi1 = &circle1;
3483   - test_value(model.flux(roi1), 0, 1.e-4, "Check model not overlapping with region");
  3480 + GSkyRegionCircle circle1 = GSkyRegionCircle(centre2, 0.1);
  3481 + GSkyRegion* roi1 = &circle1;
  3482 + test_value(model.flux(roi1), 0.0, 1.e-4, "Check model not overlapping with region");
3484 3483  
3485 3484 // Test fully contained model
3486   - GSkyRegionCircle circle2 = GSkyRegionCircle(centre2,1.);
3487   - GSkyRegion* roi2 = &circle2;
3488   - test_value(model.flux(roi2), 1., 1.e-4, "Check model fully contained in region");
  3485 + GSkyRegionCircle circle2 = GSkyRegionCircle(centre2, 1.0);
  3486 + GSkyRegion* roi2 = &circle2;
  3487 + test_value(model.flux(roi2), 1.0, 1.e-4, "Check model fully contained in region");
3489 3488  
3490 3489 // Test partially contained model
3491   - GSkyRegionCircle circle3 = GSkyRegionCircle(centre2,0.5);
3492   - GSkyRegion* roi3 = &circle3;
  3490 + GSkyRegionCircle circle3 = GSkyRegionCircle(centre2, 0.5);
  3491 + GSkyRegion* roi3 = &circle3;
3493 3492 test_value(model.flux(roi3), 0.44654, 1.e-4, "Check model partially contained in region");
3494 3493  
3495 3494 // Exit test
... ...
test/test_GModel.hpp
  View file @68de274
1 1 /***************************************************************************
2 2 * test_GModel.hpp - Test model module *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2017 by Jean-Baptiste Cayrou *
  4 + * copyright (C) 2012-2020 by Jean-Baptiste Cayrou *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -99,11 +99,11 @@ public:
99 99 void test_legacy_model_point_logparabola(void);
100 100 void test_flux(void);
101 101  
102   -private:
  102 +private:
103 103 // Private methods
104 104 void test_xml_model(const std::string& name, const std::string& filename);
105 105 void test_sky_model_content(const GModelSky& model);
106   -
  106 +
107 107 // Test files
108 108 std::string m_map_file;
109 109 std::string m_cube_file;
... ...