GammaLib  2.1.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpatialRadial.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpatialRadial.cpp - Abstract radial spatial model base class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2011-2022 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 GModelSpatialRadial.cpp
23  * @brief Abstract radial 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 "GModelSpatialRadial.hpp"
33 
34 /* __ Method name definitions ____________________________________________ */
35 #define G_READ "GModelSpatialRadial::read(GXmlElement&)"
36 #define G_WRITE "GModelSpatialRadial::write(GXmlElement&)"
37 
38 /* __ Macros _____________________________________________________________ */
39 
40 /* __ Coding definitions _________________________________________________ */
41 
42 /* __ Debug definitions __________________________________________________ */
43 
44 
45 /*==========================================================================
46  = =
47  = Constructors/destructors =
48  = =
49  ==========================================================================*/
50 
51 /***********************************************************************//**
52  * @brief Void constructor
53  ***************************************************************************/
55 {
56  // Initialise members
57  init_members();
58 
59  // Return
60  return;
61 }
62 
63 
64 /***********************************************************************//**
65  * @brief XML constructor
66  *
67  * @param[in] xml XML element.
68  *
69  * Constructs a radial spatial model by extracting information from an XML
70  * element. See the read() method for more information about the expected
71  * structure of the XML element.
72  ***************************************************************************/
75 {
76  // Initialise members
77  init_members();
78 
79  // Read information from XML element
80  read(xml);
81 
82  // Return
83  return;
84 }
85 
86 
87 /***********************************************************************//**
88  * @brief Copy constructor
89  *
90  * @param[in] model Radial spatial model.
91  ***************************************************************************/
93  GModelSpatial(model)
94 {
95  // Initialise members
96  init_members();
97 
98  // Copy members
99  copy_members(model);
100 
101  // Return
102  return;
103 }
104 
105 
106 /***********************************************************************//**
107  * @brief Destructor
108  ***************************************************************************/
110 {
111  // Free members
112  free_members();
113 
114  // Return
115  return;
116 }
117 
118 
119 /*==========================================================================
120  = =
121  = Operators =
122  = =
123  ==========================================================================*/
124 
125 /***********************************************************************//**
126  * @brief Assignment operator
127  *
128  * @param[in] model Radial spatial model.
129  * @return Radial spatial model.
130  ***************************************************************************/
132 {
133  // Execute only if object is not identical
134  if (this != &model) {
135 
136  // Copy base class members
137  this->GModelSpatial::operator=(model);
138 
139  // Free members
140  free_members();
141 
142  // Initialise members
143  init_members();
144 
145  // Copy members
146  copy_members(model);
147 
148  } // endif: object was not identical
149 
150  // Return
151  return *this;
152 }
153 
154 
155 /*==========================================================================
156  = =
157  = Public methods =
158  = =
159  ==========================================================================*/
160 
161 /***********************************************************************//**
162  * @brief Return model value
163  *
164  * @param[in] photon Incident Photon.
165  * @param[in] gradients Compute gradients?
166  * @return Value of spatial radial model.
167  *
168  * Evaluates the radial spatial model value for a specific incident
169  * @p photon.
170  *
171  * If the @p gradients flag is true the method will also compute the
172  * parameter gradients for all model parameters.
173  ***************************************************************************/
174 double GModelSpatialRadial::eval(const GPhoton& photon,
175  const bool& gradients) const
176 {
177  // Compute distance from source (in radians)
178  double theta = photon.dir().dist(dir());
179 
180  // Evaluate model
181  double value = eval(theta, photon.energy(), photon.time(), gradients);
182 
183  // Return result
184  return value;
185 }
186 
187 
188 /***********************************************************************//**
189  * @brief Read model from XML element
190  *
191  * @param[in] xml XML element.
192  *
193  * @exception GException::model_invalid_parnum
194  * Invalid number of model parameters found in XML element.
195  * @exception GException::model_invalid_parnames
196  * Invalid model parameter names found in XML element.
197  *
198  * Reads the radial source location and position angle information from an
199  * XML element in the following format
200  *
201  * <spatialModel type="...">
202  * <parameter name="RA" scale="1" value="83.6331" min="-360" max="360" free="0" />
203  * <parameter name="DEC" scale="1" value="22.0145" min="-90" max="90" free="0" />
204  * ...
205  * </spatialModel>
206  *
207  * or
208  *
209  * <spatialModel type="...">
210  * <parameter name="GLON" scale="1" value="184.5575" min="-360" max="360" free="0" />
211  * <parameter name="GLAT" scale="1" value="-5.7843" min="-90" max="90" free="0" />
212  * ...
213  * </spatialModel>
214  *
215  ***************************************************************************/
217 {
218  // Read RA/DEC parameters
219  if (gammalib::xml_has_par(xml, "RA") && gammalib::xml_has_par(xml, "DEC")) {
220 
221  // Get parameters
222  const GXmlElement* ra = gammalib::xml_get_par(G_READ, xml, "RA");
223  const GXmlElement* dec = gammalib::xml_get_par(G_READ, xml, "DEC");
224 
225  // Read parameters
226  m_lon.read(*ra);
227  m_lat.read(*dec);
228 
229  }
230 
231  // ... otherwise read GLON/GLAT parameters
232  else {
233 
234  // Get parameters
235  const GXmlElement* glon = gammalib::xml_get_par(G_READ, xml, "GLON");
236  const GXmlElement* glat = gammalib::xml_get_par(G_READ, xml, "GLAT");
237 
238  // Read parameters
239  m_lon.read(*glon);
240  m_lat.read(*glat);
241 
242  }
243 
244  // Return
245  return;
246 }
247 
248 
249 /***********************************************************************//**
250  * @brief Write model into XML element
251  *
252  * @param[in] xml XML element into which model information is written.
253  *
254  * Depending on the coordinate system, write the radial source information
255  * into an XML element with the following format
256  *
257  * <spatialModel type="...">
258  * <parameter name="RA" scale="1" value="83.6331" min="-360" max="360" free="0" />
259  * <parameter name="DEC" scale="1" value="22.0145" min="-90" max="90" free="0" />
260  * ...
261  * </spatialModel>
262  *
263  * or
264  *
265  * <spatialModel type="PointSource">
266  * <parameter name="GLON" scale="1" value="184.5575" min="-360" max="360" free="0" />
267  * <parameter name="GLAT" scale="1" value="-5.7843" min="-90" max="90" free="0" />
268  * ...
269  * </spatialModel>
270  *
271  ***************************************************************************/
273 {
274  // Check model type
276 
277  // Get or create parameters
280 
281  // Write parameters
282  m_lon.write(*lon);
283  m_lat.write(*lat);
284 
285  // Return
286  return;
287 }
288 
289 
290 /***********************************************************************//**
291  * @brief Return position of radial spatial model
292  *
293  * @return Radial spatial model sky direction.
294  *
295  * Returns the sky direction of the radial spatial model.
296  ***************************************************************************/
298 {
299  // Get longitude and latitude values
300  double lon = m_lon.value();
301  double lat = m_lat.value();
302 
303  // If longitude or latitude values have changed then update sky
304  // direction cache
305  if ((lon != m_last_lon) || (lat != m_last_lat)) {
306 
307  // Update last values
308  m_last_lon = lon;
309  m_last_lat = lat;
310 
311  // Update sky direction dependent on model coordinate system
312  if (is_celestial()) {
314  }
315  else {
317  }
318 
319  } // endif: update of sky direction cache required
320 
321  // Return sky direction
322  return (m_dir);
323 }
324 
325 
326 /***********************************************************************//**
327  * @brief Set position of radial spatial model
328  *
329  * @param[in] dir Sky direction of radial spatial model.
330  *
331  * Sets the sky direction of the radial spatial model.
332  ***************************************************************************/
334 {
335  // Assign sky direction depending on the model coordinate system
336  if (is_celestial()) {
337  m_lon.value(dir.ra_deg());
338  m_lat.value(dir.dec_deg());
339  }
340  else {
341  m_lon.value(dir.l_deg());
342  m_lat.value(dir.b_deg());
343  }
344 
345  // Return
346  return;
347 }
348 
349 
350 /*==========================================================================
351  = =
352  = Private methods =
353  = =
354  ==========================================================================*/
355 
356 /***********************************************************************//**
357  * @brief Initialise class members
358  ***************************************************************************/
360 {
361  // Initialise Right Ascension
362  m_lon.clear();
363  m_lon.name("RA");
364  m_lon.unit("deg");
365  m_lon.fix();
366  m_lon.scale(1.0);
367  m_lon.gradient(0.0);
368  m_lon.has_grad(false);
369 
370  // Initialise Declination
371  m_lat.clear();
372  m_lat.name("DEC");
373  m_lat.unit("deg");
374  m_lat.fix();
375  m_lat.scale(1.0);
376  m_lat.gradient(0.0);
377  m_lat.has_grad(false);
378 
379  // Set parameter pointer(s)
380  m_pars.clear();
381  m_pars.push_back(&m_lon);
382  m_pars.push_back(&m_lat);
383 
384  // Initialise cache
385  m_dir.clear();
386  m_last_lon = -9999.0;
387  m_last_lat = -9999.0;
388 
389  // Return
390  return;
391 }
392 
393 
394 /***********************************************************************//**
395  * @brief Copy class members
396  *
397  * @param[in] model Radial spatial model.
398  ***************************************************************************/
400 {
401  // Copy members
402  m_lon = model.m_lon;
403  m_lat = model.m_lat;
404 
405  // Set parameter pointer(s)
406  m_pars.clear();
407  m_pars.push_back(&m_lon);
408  m_pars.push_back(&m_lat);
409 
410  // Copy cache
411  m_dir = model.m_dir;
412  m_last_lon = model.m_last_lon;
413  m_last_lat = model.m_last_lat;
414 
415  // Return
416  return;
417 }
418 
419 
420 /***********************************************************************//**
421  * @brief Delete class members
422  ***************************************************************************/
424 {
425  // Return
426  return;
427 }
double dec_deg(void) const
Returns Declination in degrees.
Definition: GSkyDir.hpp:256
double m_last_lon
Last longitude.
const std::string & name(void) const
Return parameter name.
double gradient(void) const
Return parameter gradient.
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:351
Abstract radial spatial model base class interface definition.
void copy_members(const GModelSpatialRadial &model)
Copy class members.
XML element node class.
Definition: GXmlElement.hpp:48
virtual GModelSpatial & operator=(const GModelSpatial &model)
Assignment operator.
GModelSpatialRadial(void)
Void constructor.
void init_members(void)
Initialise class members.
Class that handles photons.
Definition: GPhoton.hpp:47
const double & scale(void) const
Return parameter scale.
virtual ~GModelSpatialRadial(void)
Destructor.
virtual double eval(const double &theta, const GEnergy &energy, const GTime &time, const bool &gradients=false) const =0
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
GSkyDir m_dir
Sky direction representing parameters.
const GSkyDir & dir(void) const
Return position of radial spatial model.
void free_members(void)
Delete class members.
std::vector< GModelPar * > m_pars
Parameter pointers.
void fix(void)
Fix a parameter.
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.
void clear(void)
Clear parameter.
double m_last_lat
Last latitude.
double ra_deg(void) const
Returns Right Ascension in degrees.
Definition: GSkyDir.hpp:229
const GTime & time(void) const
Return photon time.
Definition: GPhoton.hpp:134
bool is_celestial(void) const
Check if model holds celestial coordinates.
std::string type(void) const
Return model type.
void radec_deg(const double &ra, const double &dec)
Set equatorial sky direction (degrees)
Definition: GSkyDir.cpp:224
bool xml_has_par(const GXmlElement &xml, const std::string &name)
Checks if parameter with given name in XML element exists.
Definition: GTools.cpp:1596
virtual void write(GXmlElement &xml) const
Write model into XML element.
virtual GModelSpatialRadial & operator=(const GModelSpatialRadial &model)
Assignment operator.
void clear(void)
Clear sky direction.
Definition: GSkyDir.cpp:164
GModelPar m_lat
Declination or Galactic latitude (deg)
#define G_WRITE
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.
void lb_deg(const double &l, const double &b)
Set galactic sky direction (degrees)
Definition: GSkyDir.cpp:278
const GEnergy & energy(void) const
Return photon energy.
Definition: GPhoton.hpp:122
Exception handler interface definition.
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
Definition: GSkyDir.hpp:271
Abstract spatial model base class.
Abstract radial spatial model base class.
double l_deg(void) const
Return galactic longitude in degrees.
Definition: GSkyDir.hpp:175
double b_deg(void) const
Returns galactic latitude in degrees.
Definition: GSkyDir.hpp:202
#define G_READ
Sky direction class.
Definition: GSkyDir.hpp:62
const GSkyDir & dir(void) const
Return photon sky direction.
Definition: GPhoton.hpp:110
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
void xml_check_type(const std::string &origin, GXmlElement &xml, const std::string &type)
Checks the model type.
Definition: GTools.cpp:1819