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