GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpatial.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpatial.cpp - Abstract spatial model base class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2009-2020 by Juergen Knoedlseder *
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 GModelSpatial.cpp
23  * @brief Abstract spatial model base class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GException.hpp"
32 #include "GModelSpatial.hpp"
34 
35 /* __ Method name definitions ____________________________________________ */
36 #define G_ACCESS "GModelSpatial::operator[](std::string&)"
37 #define G_AT "GModelPar& GModelSpatial::at(int&)"
38 #define G_FLUX "GModelSpatial::flux(GSkyRegionCircle&, GEnergy&, GTime&)"
39 
40 /* __ Constants __________________________________________________________ */
41 const double g_kludge_radius = 1.0e-12; //!< Tiny angle (radians)
42 
43 /* __ Macros _____________________________________________________________ */
44 
45 /* __ Coding definitions _________________________________________________ */
46 
47 /* __ Debug definitions __________________________________________________ */
48 
49 
50 /*==========================================================================
51  = =
52  = Constructors/destructors =
53  = =
54  ==========================================================================*/
55 
56 /***********************************************************************//**
57  * @brief Void constructor
58  ***************************************************************************/
60 {
61  // Initialise members
62  init_members();
63 
64  // Return
65  return;
66 }
67 
68 
69 /***********************************************************************//**
70  * @brief Copy constructor
71  *
72  * @param[in] model Spatial model.
73  ***************************************************************************/
75 {
76  // Initialise members
77  init_members();
78 
79  // Copy members
80  copy_members(model);
81 
82  // Return
83  return;
84 }
85 
86 
87 /***********************************************************************//**
88  * @brief Destructor
89  ***************************************************************************/
91 {
92  // Free members
93  free_members();
94 
95  // Return
96  return;
97 }
98 
99 
100 /*==========================================================================
101  = =
102  = Operators =
103  = =
104  ==========================================================================*/
105 
106 /***********************************************************************//**
107  * @brief Assignment operator
108  *
109  * @param[in] model Spatial model.
110  * @return Spatial model.
111  ***************************************************************************/
113 {
114  // Execute only if object is not identical
115  if (this != &model) {
116 
117  // Free members
118  free_members();
119 
120  // Initialise members
121  init_members();
122 
123  // Copy members
124  copy_members(model);
125 
126  } // endif: object was not identical
127 
128  // Return
129  return *this;
130 }
131 
132 
133 /***********************************************************************//**
134  * @brief Returns model parameter
135  *
136  * @param[in] name Parameter name.
137  * @return Model parameter reference.
138  *
139  * @exception GException::par_not_found
140  * Parameter with specified name not found in container.
141  *
142  * Returns reference to the model parameter of specified @p name.
143  ***************************************************************************/
144 GModelPar& GModelSpatial::operator[](const std::string& name)
145 {
146  // Get parameter index
147  int index = 0;
148  for (; index < size(); ++index) {
149  if (m_pars[index]->name() == name) {
150  break;
151  }
152  }
153 
154  // Throw exception if parameter name was not found
155  if (index >= size()) {
156  throw GException::par_not_found(G_ACCESS, name);
157  }
158 
159  // Return reference
160  return *(m_pars[index]);
161 }
162 
163 
164 /***********************************************************************//**
165  * @brief Returns model parameter (const version)
166  *
167  * @param[in] name Parameter name.
168  * @return Model parameter reference.
169  *
170  * @exception GException::par_not_found
171  * Parameter with specified name not found in container.
172  *
173  * Returns reference to the model parameter of specified @p name.
174  ***************************************************************************/
175 const GModelPar& GModelSpatial::operator[](const std::string& name) const
176 {
177  // Get parameter index
178  int index = 0;
179  for (; index < size(); ++index) {
180  if (m_pars[index]->name() == name) {
181  break;
182  }
183  }
184 
185  // Throw exception if parameter name was not found
186  if (index >= size()) {
187  throw GException::par_not_found(G_ACCESS, name);
188  }
189 
190  // Return reference
191  return *(m_pars[index]);
192 }
193 
194 
195 /*==========================================================================
196  = =
197  = Public methods =
198  = =
199  ==========================================================================*/
200 
201 /***********************************************************************//**
202  * @brief Returns model parameter
203  *
204  * @param[in] index Parameter index [0,...,size()-1].
205  * @return Model parameter.
206  *
207  * @exception GException::out_of_range
208  * Parameter index is out of range.
209  *
210  * Returns model parameter with @p index range checking.
211  ***************************************************************************/
212 GModelPar& GModelSpatial::at(const int& index)
213 {
214  // Compile option: raise exception if index is out of range
215  if (index < 0 || index >= size()) {
216  throw GException::out_of_range(G_AT, index, 0, size()-1);
217  }
218 
219  // Return reference
220  return *(m_pars[index]);
221 }
222 
223 
224 /***********************************************************************//**
225  * @brief Returns model parameter (const version)
226  *
227  * @param[in] index Parameter index [0,...,size()-1].
228  * @return Model parameter.
229  *
230  * @exception GException::out_of_range
231  * Parameter index is out of range.
232  *
233  * Returns model parameter with @p index range checking.
234  ***************************************************************************/
235 const GModelPar& GModelSpatial::at(const int& index) const
236 {
237  // Compile option: raise exception if index is out of range
238  if (index < 0 || index >= size()) {
239  throw GException::out_of_range(G_AT, index, 0, size()-1);
240  }
241 
242  // Return reference
243  return *(m_pars[index]);
244 }
245 
246 
247 /***********************************************************************//**
248  * @brief Checks if parameter name exists
249  *
250  * @param[in] name Parameter name.
251  * @return True if parameter with specified @p name exists.
252  *
253  * Searches all parameter names for a match with the specified @p name. If
254  * the specified name has been found, true is returned.
255  ***************************************************************************/
256 bool GModelSpatial::has_par(const std::string& name) const
257 {
258  // Default found flag to false
259  bool found = false;
260 
261  // Search for parameter name
262  for (int i = 0; i < size(); ++i) {
263  if (m_pars[i]->name() == name) {
264  found = true;
265  break;
266  }
267  }
268 
269  // Return
270  return found;
271 }
272 
273 
274 /***********************************************************************//**
275  * @brief Checks if the spatial model has free parameters
276  *
277  * @return True if spatial model has free parameters.
278  ***************************************************************************/
280 {
281  // Initialise result flag
282  bool has_free_pars = false;
283 
284  // Search for free parameters
285  for (int i = 0; i < size(); ++i) {
286  if (m_pars[i]->is_free()) {
287  has_free_pars = true;
288  break;
289  }
290  }
291 
292  // Return result
293  return has_free_pars;
294 }
295 
296 
297 /***********************************************************************//**
298  * @brief Autoscale parameters
299  *
300  * Sets the scale factors for all parameters so that the values are unity.
301  ***************************************************************************/
303 {
304  // Loop over all parameters
305  for (int i = 0; i < m_pars.size(); ++i) {
306  if (m_pars[i] != NULL) {
307  m_pars[i]->autoscale();
308  }
309  }
310 
311  // Return
312  return;
313 }
314 
315 
316 /***********************************************************************//**
317  * @brief Returns model flux integrated in circular sky region
318  *
319  * @param[in] reg Sky region.
320  * @param[in] srcEng Energy.
321  * @param[in] srcTime Time.
322  * @return Flux (adimensional or ph/cm2/s).
323  *
324  * @exception GException::feature_not_implemented
325  * Regions are not circular.
326  *
327  * Computes
328  *
329  * \f[
330  * \int_{\rho_{\rm min}}^{\rho_{\rm max}} \sin \rho \times
331  * \int_{\omega} M(\rho, \omega | E, t) d\omega d\rho
332  * \f]
333  *
334  * where
335  * \f$M(\rho, \omega | E, t)\f$ is the spatial model,
336  * \f$\rho\f$ is the distance from the region centre, and
337  * \f$\omega\f$ is the position angle with respect to the connecting line
338  * between the region centre and the direction on the sky.
339  ***************************************************************************/
340 double GModelSpatial::flux(const GSkyRegion* reg,
341  const GEnergy& srcEng,
342  const GTime& srcTime) const
343 {
344  // Initialise flux
345  double flux = 0.0;
346 
347  // Continue only if region overlaps with model
348  const GSkyRegion* model_reg = this->region();
349  if (model_reg->overlaps(*reg)) {
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 regions overlap
356  if (ps_model == NULL) {
357 
358  // Thrown an exception if region is not a sky circle
359  const GSkyRegionCircle* reg_circle =
360  dynamic_cast<const GSkyRegionCircle*>(reg);
361  if (reg_circle == NULL) {
362  std::string msg = "Flux can only be computed for a circular "
363  "region.";
365  }
366 
367  // Thrown an exception if model region is not a sky circle
368  const GSkyRegionCircle* model_circle =
369  dynamic_cast<const GSkyRegionCircle*>(model_reg);
370  if (model_circle == NULL) {
371  std::string msg = "Flux can only be computed for spatial model "
372  "with region defined as circle.";
374  }
375 
376  // Model centre and radius in radians
377  GSkyDir model_centre = model_circle->centre();
378  double model_radius = model_circle->radius() * gammalib::deg2rad;
379 
380  // Distance between region and model centres in radians
381  double distance = model_centre.dist(reg_circle->centre());
382 
383  // Take the distance between centers minus radius of model region as
384  // the minimum radial integration boundary (in radians)
385  double rho_min = distance - model_radius;
386 
387  // Make sure that rho_min does not become negative
388  if (rho_min < 0.0) {
389  rho_min = 0.0;
390  }
391 
392  // Take the region radius as the maximum radial integration boundary
393  // (in radians)
394  double rho_max = reg_circle->radius() * gammalib::deg2rad;
395 
396  // Pre-compute some quantities for arclength
397  double cosdist = std::cos(distance);
398  double sindist = std::sin(distance);
399  double cosmodrad = std::cos(model_radius);
400 
401  // Setup integration kernel
402  GModelSpatial::circle_int_kern_rho integrand(this,
403  reg_circle,
404  srcEng,
405  srcTime,
406  distance,
407  cosdist,
408  sindist,
409  model_radius,
410  cosmodrad);
411  GIntegral integral(&integrand);
412 
413  // Perform integration
414  flux = integral.romberg(rho_min, rho_max);
415 
416  } // endif: source model was not a point source
417 
418  // ... otherwise, if model is point source since it overlaps
419  // with region set flux to 1
420  else {
421  flux = 1.0;
422  }
423 
424  } // endif: model overlaps with region
425 
426  // Return
427  return flux;
428 }
429 
430 
431 /*==========================================================================
432  = =
433  = Private methods =
434  = =
435  ==========================================================================*/
436 
437 /***********************************************************************//**
438  * @brief Initialise class members
439  ***************************************************************************/
441 {
442  // Initialise members
443  m_pars.clear();
444 
445  // Return
446  return;
447 }
448 
449 
450 /***********************************************************************//**
451  * @brief Copy class members
452  *
453  * @param[in] model Spatial model.
454  ***************************************************************************/
456 {
457  // Copy members
458  m_pars = model.m_pars;
459 
460  // Return
461  return;
462 }
463 
464 
465 /***********************************************************************//**
466  * @brief Delete class members
467  ***************************************************************************/
469 {
470  // Return
471  return;
472 }
473 
474 
475 /***********************************************************************//**
476  * @brief Kernel for circular sky region radial integration
477  *
478  * @param[in] rho Radial distance from region centre (radians).
479  * @return Integration kernel defined as
480  *
481  * \f[
482  * \int_{\rho_{\rm min}}^{\rho_{\rm max}} K(\rho | E, t) d\rho
483  * \f]
484  *
485  * of a spatial model over a circular region. The eval() method computes
486  *
487  * \f[
488  * K(\rho | E, t) = \sin \rho \times
489  * \int_{\omega} M(\rho, \omega | E, t) d\omega
490  * \f]
491  *
492  * where
493  * \f$M(\rho, \omega | E, t)\f$ is the spatial model,
494  * \f$\rho\f$ is the distance from the region centre, and
495  * \f$\omega\f$ is the position angle with respect to the connecting line
496  * between the region centre and the direction on the sky.
497  ***************************************************************************/
499 {
500  // Initialise model value
501  double flux = 0.0;
502 
503  // Continue only if rho is positive
504  if (rho > 0.0) {
505 
506  // Compute half length of the arc (in radians) from a circle with
507  // radius rho that intersects with model circle
508  double domega = 0.5 * gammalib::roi_arclength(rho,
509  m_dist,
510  m_cosdist,
511  m_sindist,
512  m_modrad,
513  m_cosmodrad);
514 
515  // Continue only if arc length is positive
516  if (domega > 0.0) {
517 
518  // Reduce rho by an infinite amount to avoid rounding errors
519  // at the boundary of a sharp edged model
520  double rho_kludge = rho - g_kludge_radius;
521  if (rho_kludge < 0.0) {
522  rho_kludge = 0.0;
523  }
524 
525  // Compute omega integration range
526  double omega_min = -domega;
527  double omega_max = +domega;
528 
529  // Setup integration kernel for azimuth integration
531  m_reg,
532  rho_kludge,
533  m_srcEng,
534  m_srcTime);
535 
536  // Setup integrator
537  GIntegral integral(&integrand);
538 
539  // Compute sine term for radial integration
540  double sin_rho = std::sin(rho);
541 
542  // Integrate over omega
543  flux = integral.romberg(omega_min, omega_max) * sin_rho;
544 
545  } //endif: domega was positive
546 
547  } // endif: rho was positive
548 
549  // Return
550  return flux;
551 }
552 
553 
554 /***********************************************************************//**
555  * @brief Kernel for circular sky region azimuth angle integration
556  *
557  * @param[in] omega Azimuth angle (radians).
558  * @return Integration kernel defined as
559  *
560  * \f[
561  * K(\omega | \rho, E, t) = M(\omega | \rho)
562  * \f]
563  *
564  * where
565  * \fM(\omega | \rho)\f$ is the spatial model,
566  * \f$\rho\f$ is the distance from the region centre, and
567  * \f$\omega\f$ is the position angle with respect to the connecting line
568  * between the region centre and the direction on the sky.
569  ***************************************************************************/
570 double GModelSpatial::circle_int_kern_omega::eval(const double& omega)
571 {
572  // Initialise model value
573  double flux = 0.0;
574 
575  // Precompute rho and omega in degrees
576  double omega_deg = omega * gammalib::rad2deg;
577  double rho_deg = m_rho * gammalib::rad2deg;
578 
579  // Compute sky direction corresponding to given rho, omega
580  GSkyDir dir = GSkyDir(m_reg->centre());
581 
582  // Rotate to obtain requested sky direction
583  dir.rotate_deg(omega_deg, rho_deg);
584 
585  // Set photon for this sky direction
586  GPhoton photon = GPhoton(dir, m_srcEng, m_srcTime);
587 
588  // Evaluate model for this sky direction
589  flux = m_model->eval(photon);
590 
591  // Return
592  return flux;
593 }
double m_modrad
Model radius (rad)
double roi_arclength(const double &rad, const double &dist, const double &cosdist, const double &sindist, const double &roi, const double &cosroi)
Returns length of circular arc within circular ROI.
Definition: GTools.cpp:1860
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:380
Point source spatial model class interface definition.
double m_dist
Distance model-region (rad)
int size(void) const
Return number of parameters.
GModelSpatial(void)
Void constructor.
bool has_free_pars(void) const
Checks if the spatial model has free parameters.
virtual ~GModelSpatial(void)
Destructor.
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
double eval(const double &omega)
virtual GModelSpatial & operator=(const GModelSpatial &model)
Assignment operator.
Abstract spatial model base class interface definition.
virtual GModelPar & operator[](const int &index)
Returns model parameter.
Time class.
Definition: GTime.hpp:54
GIntegral class interface definition.
Definition: GIntegral.hpp:45
Interface for the circular sky region class.
void copy_members(const GModelSpatial &model)
Copy class members.
const double & radius(void) const
Return circular region radius (in degrees)
double m_sindist
Sin of distance model-region.
Class that handles photons.
Definition: GPhoton.hpp:47
Model parameter class.
Definition: GModelPar.hpp:87
const double deg2rad
Definition: GMath.hpp:43
Abstract interface for the sky region class.
Definition: GSkyRegion.hpp:57
double eval(const double &rho)
Kernel for circular sky region radial integration.
const GModelSpatial * m_model
Spatial model.
std::vector< GModelPar * > m_pars
Parameter pointers.
void rotate_deg(const double &phi, const double &theta)
Rotate sky direction by zenith and azimuth angle.
Definition: GSkyDir.cpp:321
#define G_AT
#define G_FLUX
double flux(const GSkyRegion *reg, const GEnergy &srcEng=GEnergy(), const GTime &srcTime=GTime()) const
Returns model flux integrated in circular sky region.
double m_cosmodrad
Cos of model radius.
const GSkyRegionCircle * m_reg
Interation region.
virtual bool overlaps(const GSkyRegion &reg) const =0
const double g_kludge_radius
Tiny angle (radians)
#define G_ACCESS
Point source spatial model.
const GSkyDir & centre(void) const
Return circular region centre.
void free_members(void)
Delete class members.
GModelPar & at(const int &index)
Returns model parameter.
bool has_par(const std::string &name) const
Checks if parameter name exists.
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1226
Exception handler interface definition.
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
Definition: GSkyDir.hpp:265
Abstract spatial model base class.
void init_members(void)
Initialise class members.
const double rad2deg
Definition: GMath.hpp:44
Sky direction class.
Definition: GSkyDir.hpp:62
virtual GSkyRegion * region(void) const =0
void autoscale(void)
Autoscale parameters.
double m_cosdist
Cos of distance model-region.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48