GammaLib  2.1.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-2023 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::invalid_argument
140  * Parameter with specified name not found.
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  std::string msg = "Model parameter \""+name+"\" not found in model. "
157  "Please specify a valid model parameter name.";
159  }
160 
161  // Return reference
162  return *(m_pars[index]);
163 }
164 
165 
166 /***********************************************************************//**
167  * @brief Returns model parameter (const version)
168  *
169  * @param[in] name Parameter name.
170  * @return Model parameter reference.
171  *
172  * @exception GException::invalid_argument
173  * Parameter with specified name not found.
174  *
175  * Returns reference to the model parameter of specified @p name.
176  ***************************************************************************/
177 const GModelPar& GModelSpatial::operator[](const std::string& name) const
178 {
179  // Get parameter index
180  int index = 0;
181  for (; index < size(); ++index) {
182  if (m_pars[index]->name() == name) {
183  break;
184  }
185  }
186 
187  // Throw exception if parameter name was not found
188  if (index >= size()) {
189  std::string msg = "Model parameter \""+name+"\" not found in model. "
190  "Please specify a valid model parameter name.";
192  }
193 
194  // Return reference
195  return *(m_pars[index]);
196 }
197 
198 
199 /*==========================================================================
200  = =
201  = Public methods =
202  = =
203  ==========================================================================*/
204 
205 /***********************************************************************//**
206  * @brief Returns model parameter
207  *
208  * @param[in] index Parameter index [0,...,size()[.
209  * @return Model parameter.
210  *
211  * @exception GException::out_of_range
212  * Parameter index is out of range.
213  *
214  * Returns model parameter with @p index range checking.
215  ***************************************************************************/
216 GModelPar& GModelSpatial::at(const int& index)
217 {
218  // Compile option: raise exception if index is out of range
219  if (index < 0 || index >= size()) {
220  throw GException::out_of_range(G_AT, "Parameter index", index, size());
221  }
222 
223  // Return reference
224  return *(m_pars[index]);
225 }
226 
227 
228 /***********************************************************************//**
229  * @brief Returns model parameter (const version)
230  *
231  * @param[in] index Parameter index [0,...,size()[.
232  * @return Model parameter.
233  *
234  * @exception GException::out_of_range
235  * Parameter index is out of range.
236  *
237  * Returns model parameter with @p index range checking.
238  ***************************************************************************/
239 const GModelPar& GModelSpatial::at(const int& index) const
240 {
241  // Compile option: raise exception if index is out of range
242  if (index < 0 || index >= size()) {
243  throw GException::out_of_range(G_AT, "Parameter index", index, size());
244  }
245 
246  // Return reference
247  return *(m_pars[index]);
248 }
249 
250 
251 /***********************************************************************//**
252  * @brief Checks if parameter name exists
253  *
254  * @param[in] name Parameter name.
255  * @return True if parameter with specified @p name exists.
256  *
257  * Searches all parameter names for a match with the specified @p name. If
258  * the specified name has been found, true is returned.
259  ***************************************************************************/
260 bool GModelSpatial::has_par(const std::string& name) const
261 {
262  // Default found flag to false
263  bool found = false;
264 
265  // Search for parameter name
266  for (int i = 0; i < size(); ++i) {
267  if (m_pars[i]->name() == name) {
268  found = true;
269  break;
270  }
271  }
272 
273  // Return
274  return found;
275 }
276 
277 
278 /***********************************************************************//**
279  * @brief Checks if the spatial model has free parameters
280  *
281  * @return True if spatial model has free parameters.
282  ***************************************************************************/
284 {
285  // Initialise result flag
286  bool has_free_pars = false;
287 
288  // Search for free parameters
289  for (int i = 0; i < size(); ++i) {
290  if (m_pars[i]->is_free()) {
291  has_free_pars = true;
292  break;
293  }
294  }
295 
296  // Return result
297  return has_free_pars;
298 }
299 
300 
301 /***********************************************************************//**
302  * @brief Autoscale parameters
303  *
304  * Sets the scale factors for all parameters so that the values are unity.
305  ***************************************************************************/
307 {
308  // Loop over all parameters
309  for (int i = 0; i < m_pars.size(); ++i) {
310  if (m_pars[i] != NULL) {
311  m_pars[i]->autoscale();
312  }
313  }
314 
315  // Return
316  return;
317 }
318 
319 
320 /***********************************************************************//**
321  * @brief Returns model flux within sky region
322  *
323  * @param[in] region Sky region.
324  * @param[in] srcEng Energy.
325  * @param[in] srcTime Time.
326  * @return Flux (adimensional or ph/cm2/s).
327  *
328  * Computes
329  *
330  * \f[
331  * \int_{\rho_{\rm min}}^{\rho_{\rm max}} \sin \rho \times
332  * \int_{\omega} M(\rho, \omega | E, t) d\omega d\rho
333  * \f]
334  *
335  * where
336  * \f$M(\rho, \omega | E, t)\f$ is the spatial model,
337  * \f$\rho\f$ is the distance from the model centre, and
338  * \f$\omega\f$ is the azimuth angle around the model centre.
339  *
340  * The integration is performed over the circular region that bounds the
341  * spatial model.
342  ***************************************************************************/
343 double GModelSpatial::flux(const GSkyRegion& region,
344  const GEnergy& srcEng,
345  const GTime& srcTime) const
346 {
347  // Initialise flux
348  double flux = 0.0;
349 
350  // Make sure that m_region is set
351  set_region();
352 
353  // Continue only if region overlaps with model
354  if (m_region.overlaps(region)) {
355 
356  // Retrieve model centre and radial integration boundaries in radians
358  double rho_min = 0.0;
359  double rho_max = m_region.radius() * gammalib::deg2rad;
360 
361  // Setup integration kernel
362  GModelSpatial::circle_int_kern_rho integrand(this,
363  region,
364  centre,
365  srcEng,
366  srcTime);
367  GIntegral integral(&integrand);
368 
369  // Suppress integration warnings
370  integral.silent(true);
371 
372  // Perform integration
373  flux = integral.romberg(rho_min, rho_max);
374 
375  } // endif: model overlaps with region
376 
377  // Return
378  return flux;
379 }
380 
381 
382 /*==========================================================================
383  = =
384  = Private methods =
385  = =
386  ==========================================================================*/
387 
388 /***********************************************************************//**
389  * @brief Initialise class members
390  ***************************************************************************/
392 {
393  // Initialise members
394  m_type.clear();
395  m_region.clear();
396  m_pars.clear();
397 
398  // Return
399  return;
400 }
401 
402 
403 /***********************************************************************//**
404  * @brief Copy class members
405  *
406  * @param[in] model Spatial model.
407  ***************************************************************************/
409 {
410  // Copy members
411  m_type = model.m_type;
412  m_region = model.m_region;
413  m_pars = model.m_pars;
414 
415  // Return
416  return;
417 }
418 
419 
420 /***********************************************************************//**
421  * @brief Delete class members
422  ***************************************************************************/
424 {
425  // Return
426  return;
427 }
428 
429 
430 /***********************************************************************//**
431  * @brief Kernel for spatial model radial integration
432  *
433  * @param[in] rho Radial distance from model centre (radians).
434  * @return Integration kernel defined as
435  *
436  * \f[
437  * \int_{\rho_{\rm min}}^{\rho_{\rm max}} K(\rho | E, t) d\rho
438  * \f]
439  *
440  * of a spatial model over a circular region. The eval() method computes
441  *
442  * \f[
443  * K(\rho | E, t) = \sin \rho \times
444  * \int_{\omega} M(\rho, \omega | E, t) d\omega
445  * \f]
446  *
447  * where
448  * \f$M(\rho, \omega | E, t)\f$ is the spatial model,
449  * \f$\rho\f$ is the distance from the model centre, and
450  * \f$\omega\f$ is the azimuth angle around the model centre.
451  ***************************************************************************/
453 {
454  // Initialise model value
455  double flux = 0.0;
456 
457  // Continue only if rho is positive
458  if (rho > 0.0) {
459 
460  // Set omega integration range
461  const double omega_min = 0.0;
462  const double omega_max = gammalib::twopi;
463 
464  // Setup integration kernel for azimuth integration
466  m_region,
467  m_centre,
468  rho,
469  m_srcEng,
470  m_srcTime);
471 
472  // Setup integrator
473  GIntegral integral(&integrand);
474 
475  // Compute sine term for radial integration
476  double sin_rho = std::sin(rho);
477 
478  // Integrate over omega
479  flux = integral.romberg(omega_min, omega_max) * sin_rho;
480 
481  } // endif: rho was positive
482 
483  // Return
484  return flux;
485 }
486 
487 
488 /***********************************************************************//**
489  * @brief Kernel for spatial model azimuth angle integration
490  *
491  * @param[in] omega Azimuth angle (radians).
492  * @return Integration kernel defined as
493  *
494  * \f[
495  * K(\omega | \rho, E, t) = M(\omega | \rho)
496  * \f]
497  *
498  * where
499  * \f$M(\omega | \rho)\f$ is the spatial model,
500  * \f$\rho\f$ is the distance from the model centre, and
501  * \f$\omega\f$ is the azimuth angle around the model centre.
502  ***************************************************************************/
504 {
505  // Initialise flux
506  double flux = 0.0;
507 
508  // Compute sky direction corresponding to given rho, omega by rotating
509  // model centre
510  GSkyDir dir = m_centre;
511  dir.rotate(omega, m_rho);
512 
513  // Evaluate model if sky direction is contained in region
514  if (m_region.contains(dir)) {
515 
516  // Set photon for this sky direction
517  GPhoton photon = GPhoton(dir, m_srcEng, m_srcTime);
518 
519  // Evaluate model for this sky direction
520  flux = m_model->eval(photon);
521 
522  } // endif: sky direction was contained in region
523 
524  // Return
525  return flux;
526 }
const GSkyDir & m_centre
Model centre.
bool overlaps(const GSkyRegion &reg) const
Checks if region is overlapping with this region.
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:382
Point source spatial model class interface definition.
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.
double eval(const double &omega)
Kernel for spatial model azimuth angle integration.
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:55
GIntegral class interface definition.
Definition: GIntegral.hpp:46
std::string centre(const std::string &s, const int &n, const char &c= ' ')
Centre string to achieve a length of n characters.
Definition: GTools.cpp:1118
GSkyRegionCircle m_region
Bounding circle.
void copy_members(const GModelSpatial &model)
Copy class members.
const double & radius(void) const
Return circular region radius (in degrees)
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 spatial model radial integration.
const GModelSpatial * m_model
Spatial model.
std::vector< GModelPar * > m_pars
Parameter pointers.
#define G_AT
void rotate(const double &phi, const double &theta)
Rotate sky direction by zenith and azimuth angle.
Definition: GSkyDir.cpp:378
const double g_kludge_radius
Tiny angle (radians)
#define G_ACCESS
void silent(const bool &silent)
Set silence flag.
Definition: GIntegral.hpp:251
const GSkyDir & centre(void) const
Return circular region centre.
void free_members(void)
Delete class members.
std::string m_type
Spatial model type.
GModelPar & at(const int &index)
Returns model parameter.
bool has_par(const std::string &name) const
Checks if parameter name exists.
const GSkyRegion & m_region
Sky region.
virtual double flux(const GSkyRegion &region, const GEnergy &srcEng=GEnergy(), const GTime &srcTime=GTime()) const
Returns model flux within sky region.
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1316
bool contains(const GSkyDir &dir) const
Checks if sky direction lies within region.
Exception handler interface definition.
Abstract spatial model base class.
virtual void set_region(void) const =0
void init_members(void)
Initialise class members.
const double twopi
Definition: GMath.hpp:36
Sky direction class.
Definition: GSkyDir.hpp:62
void autoscale(void)
Autoscale parameters.
void clear(void)
Clear instance.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48