GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpatialRadialProfile.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpatialRadialProfile.cpp - Radial profile source model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2016-2018 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 GModelSpatialRadialProfile.cpp
23  * @brief Radial profile model 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 "GTools.hpp"
33 #include "GMath.hpp"
34 #include "GSkyDir.hpp"
35 #include "GXmlElement.hpp"
37 #include <iomanip>
38 
39 /* __ Constants __________________________________________________________ */
40 
41 /* __ Globals ____________________________________________________________ */
42 
43 /* __ Method name definitions ____________________________________________ */
44 
45 /* __ Macros _____________________________________________________________ */
46 
47 /* __ Coding definitions _________________________________________________ */
48 
49 /* __ Debug definitions __________________________________________________ */
50 
51 
52 /*==========================================================================
53  = =
54  = Constructors/destructors =
55  = =
56  ==========================================================================*/
57 
58 /***********************************************************************//**
59  * @brief Void constructor
60  *
61  * Constructs empty radial profile
62  ***************************************************************************/
65 {
66  // Initialise members
67  init_members();
68 
69  // Return
70  return;
71 }
72 
73 
74 /***********************************************************************//**
75  * @brief XML constructor
76  *
77  * @param[in] xml XML element.
78  *
79  * Constructs radial profile model by extracting information from an XML
80  * elements.
81  ***************************************************************************/
84 {
85  // Initialise members
86  init_members();
87 
88  // Read information from XML element
89  read(xml);
90 
91  // Return
92  return;
93 }
94 
95 
96 /***********************************************************************//**
97  * @brief Copy constructor
98  *
99  * @param[in] model Radial profile model.
100  *
101  * Copies radial profile model from another radial profile model.
102  ***************************************************************************/
104  GModelSpatialRadial(model)
105 {
106  // Initialise members
107  init_members();
108 
109  // Copy members
110  copy_members(model);
111 
112  // Return
113  return;
114 }
115 
116 
117 /***********************************************************************//**
118  * @brief Destructor
119  *
120  * Destructs radial profile model.
121  ***************************************************************************/
123 {
124  // Free members
125  free_members();
126 
127  // Return
128  return;
129 }
130 
131 
132 /*==========================================================================
133  = =
134  = Operators =
135  = =
136  ==========================================================================*/
137 
138 /***********************************************************************//**
139  * @brief Assignment operator
140  *
141  * @param[in] model Radial profile model.
142  * @return Radial profile model.
143  *
144  * Assigns radial profile model.
145  ***************************************************************************/
147 {
148  // Execute only if object is not identical
149  if (this != &model) {
150 
151  // Copy base class members
152  this->GModelSpatialRadial::operator=(model);
153 
154  // Free members
155  free_members();
156 
157  // Initialise members
158  init_members();
159 
160  // Copy members
161  copy_members(model);
162 
163  } // endif: object was not identical
164 
165  // Return
166  return *this;
167 }
168 
169 
170 /*==========================================================================
171  = =
172  = Public methods =
173  = =
174  ==========================================================================*/
175 
176 /***********************************************************************//**
177  * @brief Evaluate function (in units of sr^-1)
178  *
179  * @param[in] theta Angular distance from model centre (radians).
180  * @param[in] energy Photon energy.
181  * @param[in] time Photon arrival time.
182  * @param[in] gradients Compute gradients?
183  * @return Model value.
184  *
185  * Evaluate the radial profile model for a given angular distance @p theta
186  * the model centre, a given photon @p energy and a given @p time. The
187  * method evaluates the model by linear interpolation.
188  ***************************************************************************/
189 double GModelSpatialRadialProfile::eval(const double& theta,
190  const GEnergy& energy,
191  const GTime& time,
192  const bool& gradients) const
193 {
194  // Get pre-computation cache index
195  int icache = cache_index();
196 
197  // Get interpolation value
198  double value = m_profile[icache].nodes.interpolate(theta, m_profile[icache].values);
199 
200  // Make sure that value is not-negative
201  if (value < 0.0) {
202  value = 0.0;
203  }
204 
205  // Compile option: Check for NaN/Inf
206  #if defined(G_NAN_CHECK)
207  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
208  std::cout << "*** ERROR: GModelSpatialRadialProfile::eval";
209  std::cout << "(theta=" << theta << "): NaN/Inf encountered";
210  std::cout << " (value=" << value;
211  std::cout << ")" << std::endl;
212  }
213  #endif
214 
215  // Return value
216  return value;
217 }
218 
219 
220 /***********************************************************************//**
221  * @brief Return MC sky direction
222  *
223  * @param[in] energy Photon energy.
224  * @param[in] time Photon arrival time.
225  * @param[in,out] ran Random number generator.
226  * @return Sky direction.
227  *
228  * Draws an arbitrary sky direction from the radial profile using a rejection
229  * method.
230  ***************************************************************************/
232  const GTime& time,
233  GRan& ran) const
234 {
235  // Get pre-computation cache index
236  int icache = cache_index();
237 
238  // Simulate theta angle (degrees) using a rejection method
239  double theta = 0.0;
240  double theta_width = theta_max() - theta_min() ;
241  while (true) {
242  theta = theta_min() + ( ran.uniform() * theta_width ) ;
243  double mc = m_profile[icache].nodes.interpolate(theta, m_profile[icache].mc);
244  if ((ran.uniform() * m_profile[icache].mc_max) <= mc) {
245  break;
246  }
247  }
248  theta *= gammalib::rad2deg;
249 
250  // Simulate uniform azmiuth angle (degrees)
251  double phi = 360.0 * ran.uniform();
252 
253  // Rotate sky direction by offset
254  GSkyDir sky_dir = dir();
255  sky_dir.rotate_deg(phi, theta);
256 
257  // Return sky direction
258  return sky_dir;
259 }
260 
261 
262 /***********************************************************************//**
263  * @brief Checks where model contains specified sky direction
264  *
265  * @param[in] dir Sky direction.
266  * @param[in] margin Margin to be added to sky direction (degrees)
267  * @return True if the model contains the sky direction.
268  *
269  * Signals whether a sky direction is contained in the radial disk model.
270  ***************************************************************************/
272  const double& margin) const
273 {
274  // Compute distance to centre (radians)
275  double distance = dir.dist(this->dir());
276 
277  // Return flag
278  return ( ( distance <= theta_max() + margin*gammalib::deg2rad ) &&
279  ( distance >= theta_min() + margin*gammalib::deg2rad ) ) ;
280 }
281 
282 
283 /*==========================================================================
284  = =
285  = Private methods =
286  = =
287  ==========================================================================*/
288 
289 /***********************************************************************//**
290  * @brief Initialise class members
291  ***************************************************************************/
293 {
294  // Initialise members
295  m_coord_indep = false;
296  m_num_nodes = 100;
297  m_region.clear();
298 
299  // Initialise pre-computation cache
300  m_profile.clear();
301 
302  // Return
303  return;
304 }
305 
306 
307 /***********************************************************************//**
308  * @brief Copy class members
309  *
310  * @param[in] model Radial disk model.
311  ***************************************************************************/
313 {
314  // Copy members
316  m_num_nodes = model.m_num_nodes;
317  m_profile = model.m_profile;
318  m_region = model.m_region;
319 
320  // Return
321  return;
322 }
323 
324 
325 /***********************************************************************//**
326  * @brief Delete class members
327  ***************************************************************************/
329 {
330  // Return
331  return;
332 }
333 
334 
335 /***********************************************************************//**
336  * @brief Return index to pre-computation cache
337  *
338  * @return Index to pre-computation cache
339  *
340  * Returns the index to the pre-computation cache. If no pre-computation
341  * cache was found the method will create one and return the index to that
342  * cache.
343  ***************************************************************************/
345 {
346  // Initialise index
347  int index = -1;
348 
349  // Search for index
350  for (int i = 0; i < m_profile.size(); ++i) {
351 
352  // Initialise found flag
353  bool found = true;
354 
355  // Loop over all spatial parameters and check if they are equal
356  // to the parameters in the pre-computation cache. Break and set
357  // the found flag to false on non-equality.
358  for (int k = 0; k < m_pars.size(); ++k) {
359 
360  // Skip if model is coordinate independent and par is RA, DEC,
361  // GLON or GLAT
362  if (m_coord_indep && (m_pars[k]->name() == "RA" ||
363  m_pars[k]->name() == "DEC" ||
364  m_pars[k]->name() == "GLON" ||
365  m_pars[k]->name() == "GLAT")) {
366  continue;
367  }
368 
369  // otherwise ...
370  else if (m_pars[k]->value() != m_profile[i].pars[k]) {
371  found = false;
372  break;
373  }
374  }
375 
376  // If profile was found then set index
377  if (found) {
378  index = i;
379  break;
380  }
381 
382  } // endfor: looped over index
383 
384  // If index was not found then allocate a new pre-computation cache
385  if (index == -1) {
386 
387  // Initialise profile
388  profile prf;
389  prf.nodes.clear();
390  prf.values.clear();
391  prf.mc.clear();
392  prf.mc_max = 0.0;
393 
394  // Set profile parameters
395  for (int k = 0; k < m_pars.size(); ++k) {
396  prf.pars.push_back(m_pars[k]->value());
397  }
398 
399  // Pre-compute radial values, MC values, and mc_max. Compute
400  // also the normalization.
401  double rmax = theta_max();
402  double dr = rmax / m_num_nodes;
403  double r = theta_min();
404  double norm = 0.0;
405  for (int j = 0; j < m_num_nodes; ++j) {
406  double value = profile_value(r);
407  double mc = value * std::sin(r) * dr;
408  norm += mc;
409  if (mc > prf.mc_max) {
410  prf.mc_max = mc;
411  }
412  prf.nodes.append(r);
413  prf.values.push_back(value);
414  prf.mc.push_back(mc);
415  r += dr;
416  }
417  norm *= gammalib::twopi;
418 
419  // Normalize radial profile
420  if (norm > 0.0) {
421  double inv_norm = 1.0 / norm;
422  for (int j = 0; j < m_num_nodes; ++j) {
423  prf.values[j] *= inv_norm;
424  }
425  }
426 
427  // Push profile in pre-computation cache
428  m_profile.push_back(prf);
429 
430  // Get index to last profile
431  index = m_profile.size() - 1;
432 
433  } // endif: allocated new pre-computation cache
434 
435  // Return index
436  return (index);
437 }
438 
439 
440 /***********************************************************************//**
441  * @brief Set boundary sky region
442  ***************************************************************************/
444 {
445  // Set sky region centre to disk centre
447 
448  // Set sky region radius to maximum theta angle
450 
451  // Return
452  return;
453 }
GModelSpatialRadialProfile(void)
Void constructor.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:821
double mc_max
Maximum of profile for MC.
XML element node class interface definition.
std::vector< profile > m_profile
Pre-computation cache.
virtual double theta_min(void) const =0
Sky direction class interface definition.
std::vector< double > pars
Profile parameters.
XML element node class.
Definition: GXmlElement.hpp:47
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
Random number generator class.
Definition: GRan.hpp:44
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
GModelPar m_dec
Declination (deg)
void copy_members(const GModelSpatialRadialProfile &model)
Copy class members.
virtual GSkyDir mc(const GEnergy &energy, const GTime &time, GRan &ran) const
Return MC sky direction.
Radial profile model class interface definition.
void init_members(void)
Initialise class members.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:184
const double & radius(void) const
Return circular region radius (in degrees)
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
const double deg2rad
Definition: GMath.hpp:43
GSkyRegionCircle m_region
Bounding circle.
void free_members(void)
Delete class members.
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
virtual double theta_max(void) const =0
virtual void read(const GXmlElement &xml)
Read model from XML element.
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
virtual double eval(const double &theta, const GEnergy &energy, const GTime &time, const bool &gradients=false) const
Evaluate function (in units of sr^-1)
GModelPar m_ra
Right Ascension (deg)
int m_num_nodes
Number of profile nodes.
virtual GModelSpatialRadial & operator=(const GModelSpatialRadial &model)
Assignment operator.
const GSkyDir & centre(void) const
Return circular region centre.
GSkyDir dir(void) const
Return position of radial spatial model.
double value(void) const
Return parameter value.
virtual GModelSpatialRadialProfile & operator=(const GModelSpatialRadialProfile &model)
Assignment operator.
int cache_index(void) const
Return index to pre-computation cache.
void set_region(void) const
Set boundary sky region.
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 radial spatial model base class.
virtual ~GModelSpatialRadialProfile(void)
Destructor.
const double twopi
Definition: GMath.hpp:36
const double rad2deg
Definition: GMath.hpp:44
void append(const double &node)
Append one node to array.
Definition: GNodeArray.cpp:351
virtual bool contains(const GSkyDir &dir, const double &margin=0.0) const
Checks where model contains specified sky direction.
virtual double profile_value(const double &theta) const =0
Sky direction class.
Definition: GSkyDir.hpp:62
std::vector< double > mc
Profile for MC.
std::vector< double > values
Profile values.
Mathematical function definitions.
void clear(void)
Clear instance.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48