GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpatialElliptical.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpatialElliptical.cpp - Abstract elliptical spatial model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2013-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 GModelSpatialElliptical.cpp
23  * @brief Abstract elliptical 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"
33 
34 /* __ Method name definitions ____________________________________________ */
35 #define G_READ "GModelSpatialElliptical::read(GXmlElement&)"
36 #define G_WRITE "GModelSpatialElliptical::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 an elliptical spatial model component by extracting information
70  * from an XML element. See the read() method for more information about the
71  * expected 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  * @brief Copy constructor
88  *
89  * @param[in] model Elliptical spatial model.
90  ***************************************************************************/
93 {
94  // Initialise members
95  init_members();
96 
97  // Copy members
98  copy_members(model);
99 
100  // Return
101  return;
102 }
103 
104 
105 /***********************************************************************//**
106  * @brief Destructor
107  ***************************************************************************/
109 {
110  // Free members
111  free_members();
112 
113  // Return
114  return;
115 }
116 
117 
118 /*==========================================================================
119  = =
120  = Operators =
121  = =
122  ==========================================================================*/
123 
124 /***********************************************************************//**
125  * @brief Assignment operator
126  *
127  * @param[in] model Elliptical spatial model.
128  * @return Elliptical spatial model.
129  ***************************************************************************/
131 {
132  // Execute only if object is not identical
133  if (this != &model) {
134 
135  // Copy base class members
136  this->GModelSpatial::operator=(model);
137 
138  // Free members
139  free_members();
140 
141  // Initialise members
142  init_members();
143 
144  // Copy members
145  copy_members(model);
146 
147  } // endif: object was not identical
148 
149  // Return
150  return *this;
151 }
152 
153 
154 /*==========================================================================
155  = =
156  = Public methods =
157  = =
158  ==========================================================================*/
159 
160 /***********************************************************************//**
161  * @brief Return model value
162  *
163  * @param[in] photon Incident Photon.
164  * @param[in] gradients Compute gradients?
165  * @return Value of spatial elliptical model.
166  *
167  * Evaluates the elliptical spatial model value for a specific incident
168  * @p photon.
169  *
170  * If the @p gradients flag is true the method will also compute the
171  * parameter gradients for all model parameters.
172  ***************************************************************************/
174  const bool& gradients) const
175 {
176  // Compute distance from source and position angle (in radians)
177  const GSkyDir& srcDir = photon.dir();
178  double theta = dir().dist(srcDir);
179  double posang = dir().posang(srcDir);
180 
181  // Evaluate model
182  double value = eval(theta, posang, photon.energy(), photon.time(), gradients);
183 
184  // Return result
185  return value;
186 }
187 
188 
189 /***********************************************************************//**
190  * @brief Read model from XML element
191  *
192  * @param[in] xml XML element.
193  *
194  * @exception GException::model_invalid_parnum
195  * Invalid number of model parameters found in XML element.
196  * @exception GException::model_invalid_parnames
197  * Invalid model parameter names found in XML element.
198  *
199  * Reads the elliptical source location and position angle information from
200  * an XML element in the following format
201  *
202  * <spatialModel type="...">
203  * <parameter name="RA" scale="1" value="83.63" min="-360" max="360" free="1"/>
204  * <parameter name="DEC" scale="1" value="22.01" min="-90" max="90" free="1"/>
205  * <parameter name="PA" scale="1" value="45.0" min="-360" max="360" free="1"/>
206  * ...
207  * </spatialModel>
208  *
209  * or
210  *
211  * <spatialModel type="...">
212  * <parameter name="GLON" scale="1" value="83.63" min="-360" max="360" free="1"/>
213  * <parameter name="GLAT" scale="1" value="22.01" min="-90" max="90" free="1"/>
214  * <parameter name="PA" scale="1" value="45.0" min="-360" max="360" free="1"/>
215  * ...
216  * </spatialModel>
217  *
218  ***************************************************************************/
220 {
221  // Read RA/DEC parameters
222  if (gammalib::xml_has_par(xml, "RA") && gammalib::xml_has_par(xml, "DEC")) {
223 
224  // Get parameters
225  const GXmlElement* ra = gammalib::xml_get_par(G_READ, xml, "RA");
226  const GXmlElement* dec = gammalib::xml_get_par(G_READ, xml, "DEC");
227 
228  // Read parameters
229  m_ra.read(*ra);
230  m_dec.read(*dec);
231 
232  }
233 
234  // ... otherwise read GLON/GLAT parameters
235  else {
236 
237  // Get parameters
238  const GXmlElement* glon = gammalib::xml_get_par(G_READ, xml, "GLON");
239  const GXmlElement* glat = gammalib::xml_get_par(G_READ, xml, "GLAT");
240 
241  // Read parameters
242  m_ra.read(*glon);
243  m_dec.read(*glat);
244 
245  // Convert into RA/DEC
246  GSkyDir dir;
247  dir.lb_deg(ra(), dec()),
248  m_ra.value(dir.ra_deg());
249  m_dec.value(dir.dec_deg());
250 
251  // Set names to RA/DEC
252  m_ra.name("RA");
253  m_dec.name("DEC");
254 
255  }
256 
257  // Get other parameters
259 
260  // Read other parameters
261  m_posangle.read(*pa);
262 
263  // Return
264  return;
265 }
266 
267 
268 /***********************************************************************//**
269  * @brief Write model into XML element
270  *
271  * @param[in] xml XML element into which model information is written.
272  *
273  * @exception GException::model_invalid_spatial
274  * Existing XML element is not of requested type.
275  * @exception GException::model_invalid_parnum
276  * Invalid number of model parameters found in XML element.
277  * @exception GException::model_invalid_parnames
278  * Invalid model parameter names found in XML element.
279  *
280  * Writes the elliptical source location and position angle information into
281  * an XML element in the following format
282  *
283  * <spatialModel type="...">
284  * <parameter name="RA" scale="1" value="83.63" min="-360" max="360" free="1"/>
285  * <parameter name="DEC" scale="1" value="22.01" min="-90" max="90" free="1"/>
286  * <parameter name="PA" scale="1" value="45.0" min="-360" max="360" free="1"/>
287  * ...
288  * </spatialModel>
289  *
290  * @todo The case that an existing spatial XML element with "GLON" and "GLAT"
291  * as coordinates is not supported.
292  ***************************************************************************/
294 {
295  // Set model type
296  if (xml.attribute("type") == "") {
297  xml.attribute("type", type());
298  }
299 
300  // Verify model type
301  if (xml.attribute("type") != type()) {
303  "Elliptical model is not of type \""+type()+"\".");
304  }
305 
306  // Get or create parameters
310 
311  // Write parameters
312  m_ra.write(*ra);
313  m_dec.write(*dec);
314  m_posangle.write(*pa);
315 
316  // Return
317  return;
318 }
319 
320 
321 /***********************************************************************//**
322  * @brief Return position of elliptical spatial model
323  ***************************************************************************/
325 {
326  // Allocate sky direction
327  GSkyDir srcDir;
328 
329  // Set sky direction
330  srcDir.radec_deg(ra(), dec());
331 
332  // Return direction
333  return srcDir;
334 }
335 
336 
337 /***********************************************************************//**
338  * @brief Set position of elliptical spatial model
339  ***************************************************************************/
341 {
342  // Assign Right Ascension and Declination
343  m_ra.value(dir.ra_deg());
344  m_dec.value(dir.dec_deg());
345 
346  // Return
347  return;
348 }
349 
350 
351 /*==========================================================================
352  = =
353  = Private methods =
354  = =
355  ==========================================================================*/
356 
357 /***********************************************************************//**
358  * @brief Initialise class members
359  ***************************************************************************/
361 {
362  // Initialise Right Ascension
363  m_ra.clear();
364  m_ra.name("RA");
365  m_ra.unit("deg");
366  m_ra.fix();
367  m_ra.scale(1.0);
368  m_ra.gradient(0.0);
369  m_ra.has_grad(false);
370 
371  // Initialise Declination
372  m_dec.clear();
373  m_dec.name("DEC");
374  m_dec.unit("deg");
375  m_dec.fix();
376  m_dec.scale(1.0);
377  m_dec.gradient(0.0);
378  m_dec.has_grad(false);
379 
380  // Initialise Position Angle
381  m_posangle.clear();
382  m_posangle.name("PA");
383  m_posangle.unit("deg");
384  m_posangle.fix();
385  m_posangle.scale(1.0);
386  m_posangle.gradient(0.0);
387  m_posangle.has_grad(false);
388 
389  // Initialise semi-major axis
390  m_semimajor.clear();
391  m_semimajor.name("MajorRadius");
392  m_semimajor.unit("deg");
393  m_semimajor.value(2.778e-4); // 1 arcsec
394  m_semimajor.min(2.778e-4); // 1 arcsec
395  m_semimajor.free();
396  m_semimajor.scale(1.0);
397  m_semimajor.gradient(0.0);
398  m_semimajor.has_grad(false); // Elliptical components never have gradients
399 
400  // Initialise semi-minor axis
401  m_semiminor.clear();
402  m_semiminor.name("MinorRadius");
403  m_semiminor.unit("deg");
404  m_semiminor.value(2.778e-4); // 1 arcsec
405  m_semiminor.min(2.778e-4); // 1 arcsec
406  m_semiminor.free();
407  m_semiminor.scale(1.0);
408  m_semiminor.gradient(0.0);
409  m_semiminor.has_grad(false); // Elliptical components never have gradients
410 
411  // Set parameter pointer(s)
412  m_pars.clear();
413  m_pars.push_back(&m_ra);
414  m_pars.push_back(&m_dec);
415  m_pars.push_back(&m_posangle);
416  m_pars.push_back(&m_semimajor);
417  m_pars.push_back(&m_semiminor);
418 
419  // Return
420  return;
421 }
422 
423 
424 /***********************************************************************//**
425  * @brief Copy class members
426  *
427  * @param[in] model Radial spatial model.
428  ***************************************************************************/
430 {
431  // Copy members
432  m_ra = model.m_ra;
433  m_dec = model.m_dec;
434  m_posangle = model.m_posangle;
435  m_semiminor = model.m_semiminor;
436  m_semimajor = model.m_semimajor;
437 
438  // Set parameter pointer(s)
439  m_pars.clear();
440  m_pars.push_back(&m_ra);
441  m_pars.push_back(&m_dec);
442  m_pars.push_back(&m_posangle);
443  m_pars.push_back(&m_semimajor);
444  m_pars.push_back(&m_semiminor);
445 
446  // Return
447  return;
448 }
449 
450 
451 /***********************************************************************//**
452  * @brief Delete class members
453  ***************************************************************************/
455 {
456  // Return
457  return;
458 }
GModelPar m_semimajor
Semi-major axis of ellipse (deg)
GSkyDir dir(void) const
Return position of elliptical spatial model.
virtual void write(GXmlElement &xml) const
Write model into XML element.
void free_members(void)
Delete class members.
double dec_deg(void) const
Returns Declination in degrees.
Definition: GSkyDir.hpp:250
Abstract elliptical spatial model base class.
const std::string & name(void) const
Return parameter name.
GModelPar m_ra
Right Ascension (deg)
#define G_WRITE
double gradient(void) const
Return parameter gradient.
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:353
virtual GModelSpatialElliptical & operator=(const GModelSpatialElliptical &model)
Assignment operator.
virtual std::string type(void) const =0
XML element node class.
Definition: GXmlElement.hpp:47
virtual GModelSpatial & operator=(const GModelSpatial &model)
Assignment operator.
GModelPar m_dec
Declination (deg)
GModelPar m_posangle
Position angle from North, counterclockwise (deg)
double min(void) const
Return parameter minimum boundary.
void init_members(void)
Initialise class members.
Class that handles photons.
Definition: GPhoton.hpp:47
const double & scale(void) const
Return parameter scale.
GModelSpatialElliptical(void)
Void constructor.
virtual ~GModelSpatialElliptical(void)
Destructor.
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
const GXmlAttribute * attribute(const int &index) const
Return attribute.
double ra(void) const
Return Right Ascencion of model centre.
void free(void)
Free a parameter.
std::vector< GModelPar * > m_pars
Parameter pointers.
double dec(void) const
Return Declination of model centre.
void fix(void)
Fix a parameter.
void copy_members(const GModelSpatialElliptical &model)
Copy class members.
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
void clear(void)
Clear parameter.
double ra_deg(void) const
Returns Right Ascension in degrees.
Definition: GSkyDir.hpp:223
const GTime & time(void) const
Return photon time.
Definition: GPhoton.hpp:134
virtual void read(const GXmlElement &xml)
Read model from XML element.
GModelPar m_semiminor
Semi-minor axis of ellipse (deg)
void radec_deg(const double &ra, const double &dec)
Set equatorial sky direction (degrees)
Definition: GSkyDir.cpp:202
bool xml_has_par(const GXmlElement &xml, const std::string &name)
Checks if parameter with given name in XML element exists.
Definition: GTools.cpp:1475
Abstract elliptical spatial model base class interface definition.
void read(const GXmlElement &xml)
Extract parameter attributes from XML element.
Definition: GModelPar.cpp:229
double value(void) const
Return parameter value.
#define G_READ
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:256
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:265
Abstract spatial model base class.
virtual double eval(const double &theta, const double &posangle, const GEnergy &energy, const GTime &time, const bool &gradients=false) const =0
double posang(const GSkyDir &dir) const
Compute position angle between sky directions in radians.
Definition: GSkyDir.cpp:1022
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:1562