GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpatialRadialDisk.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpatialRadialDisk.cpp - Radial disk source model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2011-2022 by Christoph Deil *
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 GModelSpatialRadialDisk.cpp
23  * @brief Radial disk model class implementation
24  * @author Christoph Deil
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"
36 
37 /* __ Constants __________________________________________________________ */
38 
39 /* __ Globals ____________________________________________________________ */
41 const GModelSpatialRegistry g_radial_disk_registry(&g_radial_disk_seed);
42 #if defined(G_LEGACY_XML_FORMAT)
43 const GModelSpatialRadialDisk g_radial_disk_legacy_seed(true, "DiskFunction");
44 const GModelSpatialRegistry g_radial_disk_legacy_registry(&g_radial_disk_legacy_seed);
45 #endif
46 
47 /* __ Method name definitions ____________________________________________ */
48 #define G_CONSTRUCTOR "GModelSpatialRadialDisk::GModelSpatialRadialDisk("\
49  "GSkyDir&, double&, std::string&)"
50 #define G_READ "GModelSpatialRadialDisk::read(GXmlElement&)"
51 #define G_WRITE "GModelSpatialRadialDisk::write(GXmlElement&)"
52 
53 /* __ Macros _____________________________________________________________ */
54 
55 /* __ Coding definitions _________________________________________________ */
56 
57 /* __ Debug definitions __________________________________________________ */
58 
59 
60 /*==========================================================================
61  = =
62  = Constructors/destructors =
63  = =
64  ==========================================================================*/
65 
66 /***********************************************************************//**
67  * @brief Void constructor
68  *
69  * Constructs empty radial disk model.
70  ***************************************************************************/
72 {
73  // Initialise members
74  init_members();
75 
76  // Return
77  return;
78 }
79 
80 
81 /***********************************************************************//**
82  * @brief Model type constructor
83  *
84  * @param[in] dummy Dummy flag.
85  * @param[in] type Model type.
86  *
87  * Constructs empty radial disk model by specifying a model @p type.
88  ***************************************************************************/
90  const std::string& type) :
92 {
93  // Initialise members
94  init_members();
95 
96  // Set model type
97  m_type = type;
98 
99  // Return
100  return;
101 }
102 
103 
104 /***********************************************************************//**
105  * @brief Disk constructor
106  *
107  * @param[in] dir Sky position of disk centre.
108  * @param[in] radius Disk radius (degrees).
109  * @param[in] coordsys Coordinate system (either "CEL" or "GAL")
110  *
111  * @exception GException::invalid_argument
112  * Invalid @p coordsys argument specified.
113  *
114  * Constructs radial disk model from the sky position of the disk centre
115  * (@p dir) and the disk @p radius in degrees. The @p coordsys parameter
116  * specifies whether the sky direction should be interpreted in the
117  * celestial or Galactic coordinate system.
118  ***************************************************************************/
120  const double& radius,
121  const std::string& coordsys) :
123 {
124  // Throw an exception if the coordinate system is invalid
125  if ((coordsys != "CEL") && (coordsys != "GAL")) {
126  std::string msg = "Invalid coordinate system \""+coordsys+"\" "
127  "specified. Please specify either \"CEL\" or "
128  "\"GAL\".";
130  }
131 
132  // Initialise members
133  init_members();
134 
135  // Set parameter names
136  if (coordsys == "CEL") {
137  m_lon.name("RA");
138  m_lat.name("DEC");
139  }
140  else {
141  m_lon.name("GLON");
142  m_lat.name("GLAT");
143  }
144 
145  // Assign parameters
146  this->dir(dir);
147  this->radius(radius);
148 
149  // Return
150  return;
151 }
152 
153 
154 /***********************************************************************//**
155  * @brief XML constructor
156  *
157  * @param[in] xml XML element.
158  *
159  * Creates instance of radial disk model by extracting information from an
160  * XML element. See the read() method for more information about the expected
161  * structure of the XML element.
162  ***************************************************************************/
165 {
166  // Initialise members
167  init_members();
168 
169  // Read information from XML element
170  read(xml);
171 
172  // Return
173  return;
174 }
175 
176 
177 /***********************************************************************//**
178  * @brief Copy constructor
179  *
180  * @param[in] model Radial disk model.
181  ***************************************************************************/
183  GModelSpatialRadial(model)
184 {
185  // Initialise members
186  init_members();
187 
188  // Copy members
189  copy_members(model);
190 
191  // Return
192  return;
193 }
194 
195 
196 /***********************************************************************//**
197  * @brief Destructor
198  ***************************************************************************/
200 {
201  // Free members
202  free_members();
203 
204  // Return
205  return;
206 }
207 
208 
209 /*==========================================================================
210  = =
211  = Operators =
212  = =
213  ==========================================================================*/
214 
215 /***********************************************************************//**
216  * @brief Assignment operator
217  *
218  * @param[in] model Radial disk model.
219  * @return Radial disk model.
220  ***************************************************************************/
222 {
223  // Execute only if object is not identical
224  if (this != &model) {
225 
226  // Copy base class members
227  this->GModelSpatialRadial::operator=(model);
228 
229  // Free members
230  free_members();
231 
232  // Initialise members
233  init_members();
234 
235  // Copy members
236  copy_members(model);
237 
238  } // endif: object was not identical
239 
240  // Return
241  return *this;
242 }
243 
244 
245 /*==========================================================================
246  = =
247  = Public methods =
248  = =
249  ==========================================================================*/
250 
251 /***********************************************************************//**
252  * @brief Clear radial disk model
253  ***************************************************************************/
255 {
256  // Free class members (base and derived classes, derived class first)
257  free_members();
260 
261  // Initialise members
264  init_members();
265 
266  // Return
267  return;
268 }
269 
270 
271 /***********************************************************************//**
272  * @brief Clone radial disk model
273  *
274  * @return Pointer to deep copy of radial disk model.
275  ***************************************************************************/
277 {
278  // Clone radial disk model
279  return new GModelSpatialRadialDisk(*this);
280 }
281 
282 
283 /***********************************************************************//**
284  * @brief Evaluate function (in units of sr^-1)
285  *
286  * @param[in] theta Angular distance from disk centre (radians).
287  * @param[in] energy Photon energy.
288  * @param[in] time Photon arrival time.
289  * @param[in] gradients Compute gradients?
290  * @return Model value.
291  *
292  * Evaluates the spatial component for a disk source model using
293  *
294  * \f[
295  * S_{\rm p}(\vec{p} | E, t) = \left \{
296  * \begin{array}{l l}
297  * \displaystyle
298  * {\tt m\_norm}
299  * & \mbox{if $\theta \le $ radius} \\
300  * \\
301  * \displaystyle
302  * 0 & \mbox{if $\theta > $ radius}
303  * \end{array}
304  * \right .
305  * \f]
306  *
307  * where
308  * - \f$\theta\f$ is the angular separation between disk centre and the
309  * sky direction of interest, and
310  * - \f${\tt m\_norm} = \frac{1}{2 \pi (1 - \cos r)} \f$ is a normalization
311  * constant (see the update() method for details).
312  *
313  * The method will not compute analytical parameter gradients, even if the
314  * @p gradients argument is set to true. Radial disk parameter gradients
315  * need to be computed numerically.
316  ***************************************************************************/
317 double GModelSpatialRadialDisk::eval(const double& theta,
318  const GEnergy& energy,
319  const GTime& time,
320  const bool& gradients) const
321 {
322  // Update precomputation cache
323  update();
324 
325  // Set value
326  double value = (theta <= m_radius_rad) ? m_norm : 0.0;
327 
328  // Compile option: Check for NaN/Inf
329  #if defined(G_NAN_CHECK)
330  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
331  std::cout << "*** ERROR: GModelSpatialRadialDisk::eval";
332  std::cout << "(theta=" << theta << "): NaN/Inf encountered";
333  std::cout << " (value=" << value;
334  std::cout << ", m_radius_rad=" << m_radius_rad;
335  std::cout << ", m_norm=" << m_norm;
336  std::cout << ")" << std::endl;
337  }
338  #endif
339 
340  // Return value
341  return value;
342 }
343 
344 
345 /***********************************************************************//**
346  * @brief Return MC sky direction
347  *
348  * @param[in] energy Photon energy.
349  * @param[in] time Photon arrival time.
350  * @param[in,out] ran Random number generator.
351  * @return Sky direction.
352  *
353  * Draws an arbitrary sky direction from the 2D disk distribution as function
354  * of the photon @p energy and arrival @p time.
355  ***************************************************************************/
357  const GTime& time,
358  GRan& ran) const
359 {
360  // Simulate offset from photon arrival direction
361  double cosrad = std::cos(radius() * gammalib::deg2rad);
362  double theta = std::acos(1.0 - ran.uniform() * (1.0 - cosrad)) *
364  double phi = 360.0 * ran.uniform();
365 
366  // Rotate sky direction by offset
367  GSkyDir sky_dir = dir();
368  sky_dir.rotate_deg(phi, theta);
369 
370  // Return sky direction
371  return sky_dir;
372 }
373 
374 
375 /***********************************************************************//**
376  * @brief Checks where model contains specified sky direction
377  *
378  * @param[in] dir Sky direction.
379  * @param[in] margin Margin to be added to sky direction (degrees)
380  * @return True if the model contains the sky direction.
381  *
382  * Signals whether a sky direction is contained in the radial disk model.
383  ***************************************************************************/
385  const double& margin) const
386 {
387  // Compute distance to centre (radians)
388  double distance = dir.dist(this->dir());
389 
390  // Return flag
391  return (distance <= theta_max() + margin*gammalib::deg2rad);
392 }
393 
394 
395 /***********************************************************************//**
396  * @brief Return maximum model radius (in radians)
397  *
398  * @return Maximum model radius (in radians).
399  ***************************************************************************/
401 {
402  // Return value
403  return (radius() * gammalib::deg2rad);
404 }
405 
406 
407 /***********************************************************************//**
408  * @brief Read model from XML element
409  *
410  * @param[in] xml XML element.
411  *
412  * Reads the radial disk model information from an XML element. The XML
413  * element shall have either the format
414  *
415  * <spatialModel type="RadialDisk">
416  * <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
417  * <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
418  * <parameter name="Radius" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
419  * </spatialModel>
420  *
421  * or
422  *
423  * <spatialModel type="RadialDisk">
424  * <parameter name="GLON" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
425  * <parameter name="GLAT" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
426  * <parameter name="Radius" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
427  * </spatialModel>
428  *
429  * @todo Implement a test of the radius and radius boundary. The radius
430  * and radius minimum should be >0.
431  ***************************************************************************/
433 {
434  // Verify that XML element has exactly 3 parameters
436 
437  // Read disk location
439 
440  // Get parameters
442 
443  // Read parameters
444  m_radius.read(*radius);
445 
446  // Return
447  return;
448 }
449 
450 
451 /***********************************************************************//**
452  * @brief Write model into XML element
453  *
454  * @param[in] xml XML element into which model information is written.
455  *
456  * @exception GException::model_invalid_parnum
457  * Invalid number of model parameters found in XML element.
458  * @exception GException::model_invalid_parnames
459  * Invalid model parameter names found in XML element.
460  *
461  * Writes the radial disk model information into an XML element. The XML
462  * element will have the format
463  *
464  * <spatialModel type="RadialDisk">
465  * <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
466  * <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
467  * <parameter name="Radius" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
468  * </spatialModel>
469  *
470  ***************************************************************************/
472 {
473  // Verify model type
475 
476  // Write disk location
478 
479  // Get or create parameters
481 
482  // Write parameters
483  m_radius.write(*radius);
484 
485  // Return
486  return;
487 }
488 
489 
490 /***********************************************************************//**
491  * @brief Print information
492  *
493  * @param[in] chatter Chattiness.
494  * @return String containing model information.
495  ***************************************************************************/
496 std::string GModelSpatialRadialDisk::print(const GChatter& chatter) const
497 {
498  // Initialise result string
499  std::string result;
500 
501  // Continue only if chatter is not silent
502  if (chatter != SILENT) {
503 
504  // Append header
505  result.append("=== GModelSpatialRadialDisk ===");
506 
507  // Append parameters
508  result.append("\n"+gammalib::parformat("Number of parameters"));
509  result.append(gammalib::str(size()));
510  for (int i = 0; i < size(); ++i) {
511  result.append("\n"+m_pars[i]->print(chatter));
512  }
513 
514  } // endif: chatter was not silent
515 
516  // Return result
517  return result;
518 }
519 
520 
521 /*==========================================================================
522  = =
523  = Private methods =
524  = =
525  ==========================================================================*/
526 
527 /***********************************************************************//**
528  * @brief Initialise class members
529  ***************************************************************************/
531 {
532  // Initialise model type
533  m_type = "RadialDisk";
534 
535  // Initialise Radius
536  m_radius.clear();
537  m_radius.name("Radius");
538  m_radius.unit("deg");
539  m_radius.value(2.778e-4); // 1 arcsec
540  m_radius.min(2.778e-4); // 1 arcsec
541  m_radius.free();
542  m_radius.scale(1.0);
543  m_radius.gradient(0.0);
544  m_radius.has_grad(false); // Radial components never have gradients
545 
546  // Set parameter pointer(s)
547  m_pars.push_back(&m_radius);
548 
549  // Initialise precomputation cache. Note that zero values flag
550  // uninitialised as a zero radius is not meaningful
551  m_last_radius = 0.0;
552  m_radius_rad = 0.0;
553  m_norm = 0.0;
554 
555  // Return
556  return;
557 }
558 
559 
560 /***********************************************************************//**
561  * @brief Copy class members
562  *
563  * @param[in] model Radial disk model.
564  *
565  * We do not have to push back the members on the parameter stack as this
566  * should have been done by init_members() that was called before. Otherwise
567  * we would have the radius twice on the stack.
568  ***************************************************************************/
570 {
571  // Copy members
572  m_type = model.m_type; // Needed to conserve model type
573  m_radius = model.m_radius;
574 
575  // Copy precomputation cache
577  m_radius_rad = model.m_radius_rad;
578  m_norm = model.m_norm;
579 
580  // Return
581  return;
582 }
583 
584 
585 /***********************************************************************//**
586  * @brief Delete class members
587  ***************************************************************************/
589 {
590  // Return
591  return;
592 }
593 
594 
595 /***********************************************************************//**
596  * @brief Update precomputation cache
597  *
598  * Computes the normalization
599  * \f[
600  * {\tt m\_norm} = \frac{1}{2 \pi (1 - \cos r)}
601  * \f]
602  *
603  * Note that this is the correct normalization on the sphere for any
604  * disk radius r. For small r it is very similar to the cartesian
605  * approximation you might have expected:
606  * \f[{\tt m\_norm} = \frac{1}{\pi r ^ 2}\f]
607  ***************************************************************************/
609 {
610  // Update if radius has changed
611  if (m_last_radius != radius()) {
612 
613  // Store last values
614  m_last_radius = radius();
615 
616  // Compute disk radius in radians
618 
619  // Perform precomputations
620  double denom = gammalib::twopi * (1 - std::cos(m_radius_rad));
621  m_norm = (denom > 0.0) ? 1.0 / denom : 0.0;
622 
623  } // endif: update required
624 
625  // Return
626  return;
627 }
628 
629 
630 /***********************************************************************//**
631  * @brief Set boundary sky region
632  ***************************************************************************/
634 {
635  // Set sky region circle
637 
638  // Set region (circumvent const correctness)
639  const_cast<GModelSpatialRadialDisk*>(this)->m_region = region;
640 
641  // Return
642  return;
643 }
virtual GModelSpatialRadialDisk & operator=(const GModelSpatialRadialDisk &model)
Assignment operator.
virtual double theta_max(void) const
Return maximum model radius (in radians)
virtual void write(GXmlElement &xml) const
Write model into XML element.
void update(void) const
Update precomputation cache.
const std::string & name(void) const
Return parameter name.
double gradient(void) const
Return parameter gradient.
int size(void) const
Return number of parameters.
#define G_READ
double radius(void) const
Return disk radius.
virtual void clear(void)
Clear radial disk model.
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:351
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1190
XML element node class.
Definition: GXmlElement.hpp:48
const GModelSpatialRadialDisk g_radial_disk_seed
double acos(const double &arg)
Computes acos by avoiding NaN due to rounding errors.
Definition: GMath.cpp:69
Random number generator class.
Definition: GRan.hpp:44
#define G_WRITE
Time class.
Definition: GTime.hpp:55
Gammalib tools definition.
Interface for the circular sky region class.
Interface definition for the spatial model registry class.
double min(void) const
Return parameter minimum boundary.
void init_members(void)
Initialise class members.
GModelPar m_radius
Disk radius (degrees)
virtual ~GModelSpatialRadialDisk(void)
Destructor.
GSkyRegionCircle m_region
Bounding circle.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:201
const double & scale(void) const
Return parameter scale.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:184
const double deg2rad
Definition: GMath.hpp:43
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
virtual double eval(const double &theta, const GEnergy &energy, const GTime &time, const bool &gradients=false) const
Evaluate function (in units of sr^-1)
const GSkyDir & dir(void) const
Return position of radial spatial model.
virtual void read(const GXmlElement &xml)
Read model from XML element.
void free_members(void)
Delete class members.
void free(void)
Free a parameter.
Radial disk model class interface definition.
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:424
GXmlElement * xml_need_par(const std::string &origin, GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1637
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
void clear(void)
Clear parameter.
Spatial model registry class definition.
GModelSpatialRadialDisk(void)
Void constructor.
GChatter
Definition: GTypemaps.hpp:33
virtual GSkyDir mc(const GEnergy &energy, const GTime &time, GRan &ran) const
Return MC sky direction.
std::string type(void) const
Return model type.
void xml_check_parnum(const std::string &origin, const GXmlElement &xml, const int &number)
Checks number of parameters.
Definition: GTools.cpp:1777
virtual std::string print(const GChatter &chatter=NORMAL) const
Print information.
virtual void set_region(void) const
Set boundary sky region.
virtual void write(GXmlElement &xml) const
Write model into XML element.
const GSkyRegion * region(void) const
Return boundary sky region.
virtual GModelSpatialRadial & operator=(const GModelSpatialRadial &model)
Assignment operator.
void free_members(void)
Delete class members.
std::string m_type
Spatial model type.
GModelPar m_lat
Declination or Galactic latitude (deg)
virtual GModelSpatialRadialDisk * clone(void) const
Clone radial disk model.
#define G_CONSTRUCTOR
void read(const GXmlElement &xml)
Extract parameter attributes from XML element.
Definition: GModelPar.cpp:229
double value(void) const
Return parameter value.
GModelPar m_lon
Right Ascension or Galactic longitude (deg)
const std::string & unit(void) const
Return parameter unit.
Exception handler interface definition.
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
Definition: GSkyDir.hpp:271
void free_members(void)
Delete class members.
Abstract radial spatial model base class.
void init_members(void)
Initialise class members.
const double twopi
Definition: GMath.hpp:36
double m_radius_rad
Radius in radians.
const double rad2deg
Definition: GMath.hpp:44
double m_last_radius
Last disk radius.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
void copy_members(const GModelSpatialRadialDisk &model)
Copy class members.
void init_members(void)
Initialise class members.
Sky direction class.
Definition: GSkyDir.hpp:62
virtual bool contains(const GSkyDir &dir, const double &margin=0.0) const
Checks where model contains specified sky direction.
const GXmlElement * xml_get_par(const std::string &origin, const GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1689
Mathematical function definitions.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489
void xml_check_type(const std::string &origin, GXmlElement &xml, const std::string &type)
Checks the model type.
Definition: GTools.cpp:1819