GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpatialEllipticalGauss.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpatialEllipticalGauss.cpp - Elliptical gauss source model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2015-2018 by Michael Mayer *
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 GModelSpatialEllipticalGauss.cpp
23  * @brief Elliptical gauss model class implementation
24  * @author Michael Mayer
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"
37 
38 /* __ Constants __________________________________________________________ */
39 namespace {
40  const double c_theta_max = 3.0; //!< semiaxis multiplied for theta_max()
41  const double c_max_exponent = 0.5 * c_theta_max * c_theta_max;
42  const double c_fraction = 1.0 - std::exp(-c_max_exponent);
43 }
44 
45 /* __ Globals ____________________________________________________________ */
47 const GModelSpatialRegistry g_elliptical_gauss_registry(&g_elliptical_gauss_seed);
48 #if defined(G_LEGACY_XML_FORMAT)
49 const GModelSpatialEllipticalGauss g_elliptical_gauss_legacy_seed(true, "EllipticalGauss");
50 const GModelSpatialRegistry g_elliptical_gauss_legacy_registry(&g_elliptical_gauss_legacy_seed);
51 #endif
52 
53 /* __ Method name definitions ____________________________________________ */
54 #define G_READ "GModelSpatialEllipticalGauss::read(GXmlElement&)"
55 #define G_WRITE "GModelSpatialEllipticalGauss::write(GXmlElement&)"
56 
57 /* __ Macros _____________________________________________________________ */
58 
59 /* __ Coding definitions _________________________________________________ */
60 #define G_SMALL_ANGLE_APPROXIMATION //!< Use small angle approximation
61 
62 /* __ Debug definitions __________________________________________________ */
63 
64 
65 /*==========================================================================
66  = =
67  = Constructors/destructors =
68  = =
69  ==========================================================================*/
70 
71 /***********************************************************************//**
72  * @brief Void constructor
73  ***************************************************************************/
76 {
77  // Initialise members
78  init_members();
79 
80  // Return
81  return;
82 }
83 
84 
85 /***********************************************************************//**
86  * @brief Model type constructor
87  *
88  * @param[in] dummy Dummy flag.
89  * @param[in] type Model type.
90  *
91  * Constructs empty elliptical Gaussian model by specifying a model @p type.
92  ***************************************************************************/
94  const std::string& type) :
96 {
97  // Initialise members
98  init_members();
99 
100  // Set model type
101  m_type = type;
102 
103  // Return
104  return;
105 }
106 
107 
108 /***********************************************************************//**
109  * @brief Elliptical Gaussian constructor
110  *
111  * @param[in] dir Centre of elliptical Gaussian.
112  * @param[in] semimajor Semi-major axis (degrees).
113  * @param[in] semiminor Semi-minor axis (degrees).
114  * @param[in] posangle Position angle of semi-major axis (degrees).
115  *
116  * Construct elliptical Gaussian model from the ellipse centre direction
117  * (@p dir), the @p semimajor and @p semiminor axes, and the position
118  * angle (@p posangle).
119  ***************************************************************************/
121  const double& semimajor,
122  const double& semiminor,
123  const double& posangle) :
125 {
126  // Initialise members
127  init_members();
128 
129  // Assign parameters
130  this->dir(dir);
131  this->semiminor(semiminor);
132  this->semimajor(semimajor);
133  this->posangle(posangle);
134 
135  // Return
136  return;
137 }
138 
139 
140 /***********************************************************************//**
141  * @brief XML constructor
142  *
143  * @param[in] xml XML element.
144  *
145  * Constructs elliptical Gaussian model by extracting information from an XML
146  * element. See the read() method for more information about the expected
147  * structure of the XML element.
148  ***************************************************************************/
151 {
152  // Initialise members
153  init_members();
154 
155  // Read information from XML element
156  read(xml);
157 
158  // Return
159  return;
160 }
161 
162 
163 /***********************************************************************//**
164  * @brief Copy constructor
165  *
166  * @param[in] model Elliptical Gaussian model.
167  ***************************************************************************/
170 {
171  // Initialise members
172  init_members();
173 
174  // Copy members
175  copy_members(model);
176 
177  // Return
178  return;
179 }
180 
181 
182 /***********************************************************************//**
183  * @brief Destructor
184  ***************************************************************************/
186 {
187  // Free members
188  free_members();
189 
190  // Return
191  return;
192 }
193 
194 
195 /*==========================================================================
196  = =
197  = Operators =
198  = =
199  ==========================================================================*/
200 
201 /***********************************************************************//**
202  * @brief Assignment operator
203  *
204  * @param[in] model Elliptical gauss model.
205  * @return Elliptical Gaussian model.
206  ***************************************************************************/
208 {
209  // Execute only if object is not identical
210  if (this != &model) {
211 
212  // Copy base class members
214 
215  // Free members
216  free_members();
217 
218  // Initialise members
219  init_members();
220 
221  // Copy members
222  copy_members(model);
223 
224  } // endif: object was not identical
225 
226  // Return
227  return *this;
228 }
229 
230 
231 /*==========================================================================
232  = =
233  = Public methods =
234  = =
235  ==========================================================================*/
236 
237 /***********************************************************************//**
238  * @brief Clear elliptical Gaussian model
239  ***************************************************************************/
241 {
242  // Free class members (base and derived classes, derived class first)
243  free_members();
246 
247  // Initialise members
250  init_members();
251 
252  // Return
253  return;
254 }
255 
256 
257 /***********************************************************************//**
258  * @brief Clone elliptical Gaussian model
259  *
260  * @return Pointer to deep copy of elliptical Gaussian model.
261  ***************************************************************************/
263 {
264  // Clone elliptical Gaussian model
265  return new GModelSpatialEllipticalGauss(*this);
266 }
267 
268 
269 /***********************************************************************//**
270  * @brief Evaluate function (in units of sr^-1)
271  *
272  * @param[in] theta Angular distance from disk centre (radians).
273  * @param[in] posangle Position angle (counterclockwise from North) (radians).
274  * @param[in] energy Photon energy.
275  * @param[in] time Photon arrival time.
276  * @param[in] gradients Compute gradients?
277  * @return Model value.
278  *
279  * Evaluates the spatial component for an elliptical Gaussian source model.
280  *
281  * The elliptical Gaussian function is defined by
282  * \f$S_{\rm p}(\theta, \phi | E, t)\f$, where
283  * \f$\theta\f$ is the angular separation between elliptical Gaussian centre
284  * and the actual location and \f$\phi\f$ the position angle with respect to
285  * the model centre, counted counterclockwise from North.
286  *
287  * The function \f$S_{\rm p}(\theta, \phi | E, t)\f$ is given by
288  *
289  * \f[
290  * S_{\rm p}(\theta, \phi | E, t) = {\tt m\_norm} \times
291  * \exp \left( -\frac{\theta^2}{2 r_{\rm eff}^2} \right)
292  * \f]
293  *
294  * where the effective ellipse radius \f$r_{\rm eff}\f$ towards a given
295  * position angle is given by
296  *
297  * \f[
298  * r_{\rm eff} = \frac{ab}
299  * {\sqrt{\left( a \sin (\phi - \phi_0) \right)^2} +
300  * \sqrt{\left( b \cos (\phi - \phi_0) \right)^2}}
301  * \f]
302  * and
303  * \f$a\f$ is the semi-major axis of the ellipse,
304  * \f$b\f$ is the semi-minor axis, and
305  * \f$\phi_0\f$ is the position angle of the ellipse, counted
306  * counterclockwise from North.
307  * \f${\tt m\_norm}\f$ is a normalization constant given by
308  *
309  * \f[
310  * {\tt m\_norm} = \frac{1}{2 \pi \times a \times b}
311  * \f]
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  * @warning
318  * The above formula for an elliptical Gaussian are accurate for small
319  * angles, with semimajor and semiminor axes below a few degrees. This
320  * should be fine for most use cases, but you should be aware that the
321  * ellipse will get distorted for larger angles.
322  *
323  * @warning
324  * For numerical reasons the elliptical Gaussian will be truncated for
325  * \f$\theta\f$ angles that correspond to 3.0 times the effective ellipse
326  * radius.
327  ***************************************************************************/
328 double GModelSpatialEllipticalGauss::eval(const double& theta,
329  const double& posangle,
330  const GEnergy& energy,
331  const GTime& time,
332  const bool& gradients) const
333 {
334  // Initialise value
335  double value = 0.0;
336 
337  // Continue only if we're inside circle enclosing the ellipse
338  if (theta <= theta_max()) {
339 
340  // Update precomputation cache
341  update();
342 
343  // Perform computations to compute exponent
344  #if defined(G_SMALL_ANGLE_APPROXIMATION)
345  double rel_posangle = posangle - this->posangle() * gammalib::deg2rad;
346  double cosinus = std::cos(rel_posangle);
347  double sinus = std::sin(rel_posangle);
348  double arg1 = m_minor_rad * cosinus;
349  double arg2 = m_major_rad * sinus;
350  double r_ellipse = m_minor_rad * m_major_rad /
351  std::sqrt(arg1*arg1 + arg2*arg2); //!< small angle
352  double r_relative = theta/r_ellipse;
353  double exponent = 0.5*r_relative*r_relative;
354  #else
355  // Perform computations
356  double sinphi = std::sin(posangle);
357  double cosphi = std::cos(posangle);
358 
359  // Compute help term
360  double term1 = m_term1 * cosphi * cosphi;
361  double term2 = m_term2 * sinphi * sinphi;
362  double term3 = m_term3 * sinphi * cosphi;
363 
364  // Compute exponent
365  double exponent = theta * theta * (term1 + term2 + term3);
366  #endif
367 
368  // Set value if the exponent is within the truncation region
369  if (exponent <= c_max_exponent) {
370  value = m_norm * std::exp(-exponent);
371  }
372 
373  // Compile option: Check for NaN/Inf
374  #if defined(G_NAN_CHECK)
375  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
376  std::cout << "*** ERROR: GModelSpatialEllipticalGauss::eval";
377  std::cout << "(theta=" << theta << "): NaN/Inf encountered";
378  std::cout << "(posangle=" << posangle << "): NaN/Inf encountered";
379  std::cout << " (value=" << value;
380  std::cout << ", MinorAxis^2=" << m_minor2;
381  std::cout << ", MaxjorAxis^2=" << m_major2;
382  std::cout << ", m_norm=" << m_norm;
383  std::cout << ")" << std::endl;
384  }
385  #endif
386 
387  } // endif: position was inside enclosing circle
388 
389  // Return value
390  return value;
391 }
392 
393 
394 /***********************************************************************//**
395  * @brief Returns MC sky direction
396  *
397  * @param[in] energy Photon energy.
398  * @param[in] time Photon arrival time.
399  * @param[in,out] ran Random number generator.
400  * @return Sky direction.
401  *
402  * Draws an arbitrary sky direction from the elliptical Gaussian model.
403  *
404  * @warning
405  * For numerical reasons the elliptical Gaussian will be truncated for
406  * \f$\theta\f$ angles that correspond to 3.0 times the effective ellipse
407  * radius.
408  ***************************************************************************/
410  const GTime& time,
411  GRan& ran) const
412 {
413  // Update precomputation cache
414  update();
415 
416  // Initialise photon
417  GPhoton photon;
418  photon.energy(energy);
419  photon.time(time);
420 
421  // Draw gaussian offset from each axis
422  double ran_major;
423  double ran_minor;
424  do {
425  ran_major = ran.normal();
426  } while (ran_major > c_theta_max);
427  do {
428  ran_minor = ran.normal();
429  } while (ran_minor > c_theta_max);
430  double theta1 = semimajor() * ran_major;
431  double theta2 = semiminor() * ran_minor;
432 
433  // Compute total offset from model centre in small angle approximation
434  double theta = std::sqrt(theta1 * theta1 + theta2 * theta2);
435 
436  // Compute rotation angle, taking into account given position angle
437  double phi = gammalib::atan2d(theta2, theta1) + posangle();
438 
439  // Rotate sky direction by offset
440  GSkyDir sky_dir = dir();
441  sky_dir.rotate_deg(phi , theta);
442 
443  // Set photon sky direction
444  photon.dir(sky_dir);
445 
446  // Return photon direction
447  return (photon.dir());
448 }
449 
450 
451 /***********************************************************************//**
452  * @brief Checks whether model contains specified sky direction
453  *
454  * @param[in] dir Sky direction.
455  * @param[in] margin Margin to be added to sky direction (degrees)
456  * @return True if the model contains the sky direction.
457  *
458  * Signals whether a sky direction is contained in the elliptical gauss
459  * model.
460  *
461  * @todo Implement correct evaluation of effective ellipse radius.
462  ***************************************************************************/
464  const double& margin) const
465 {
466  // Compute distance to centre (radian)
467  double distance = dir.dist(this->dir());
468 
469  // Return flag
470  return (distance <= theta_max() + margin*gammalib::deg2rad);
471 }
472 
473 
474 /***********************************************************************//**
475  * @brief Return maximum model radius (in radians)
476  *
477  * @return Returns maximum model radius.
478  *
479  * Returns the maximum of 3.0 semimajor() and 3.0 semiminor() as
480  * approximate edge of the Gaussian. This limit is of course arbitrary, but
481  * allows to limit the integration region for response computation. The value
482  * of 3.0 has been determined by experiment (#1561).
483  ***************************************************************************/
485 {
486  // Set maximum model radius
487  double theta_max = (semimajor() > semiminor())
488  ? semimajor() * gammalib::deg2rad * c_theta_max
489  : semiminor() * gammalib::deg2rad * c_theta_max;
490 
491  // Return value
492  return theta_max;
493 }
494 
495 
496 /***********************************************************************//**
497  * @brief Read model from XML element
498  *
499  * @param[in] xml XML element.
500  *
501  * @exception GException::model_invalid_parnum
502  * Invalid number of model parameters found in XML element.
503  * @exception GException::model_invalid_parnames
504  * Invalid model parameter names found in XML element.
505  *
506  * Reads the elliptical gauss model information from an XML element. The XML
507  * element shall have either the format
508  *
509  * <spatialModel type="EllipticalGaussian">
510  * <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
511  * <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
512  * <parameter name="PA" scale="1.0" value="45.0" min="-360" max="360" free="1"/>
513  * <parameter name="MinorRadius" scale="1.0" value="0.5" min="0.001" max="10" free="1"/>
514  * <parameter name="MajorRadius" scale="1.0" value="2.0" min="0.001" max="10" free="1"/>
515  * </spatialModel>
516  *
517  * or
518  *
519  * <spatialModel type="EllipticalGaussian">
520  * <parameter name="GLON" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
521  * <parameter name="GLAT" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
522  * <parameter name="PA" scale="1.0" value="45.0" min="-360" max="360" free="1"/>
523  * <parameter name="MinorRadius" scale="1.0" value="0.5" min="0.001" max="10" free="1"/>
524  * <parameter name="MajorRadius" scale="1.0" value="2.0" min="0.001" max="10" free="1"/>
525  * </spatialModel>
526  *
527  * @todo Implement a test of the ellipse boundary. The axes
528  * and axes minimum should be >0.
529  ***************************************************************************/
531 {
532  // Determine number of parameter nodes in XML element
533  int npars = xml.elements("parameter");
534 
535  // Verify that XML element has exactly 5 parameters
536  if (xml.elements() != 5 || npars != 5) {
538  "Elliptical gauss model requires exactly 5 parameters.");
539  }
540 
541  // Read gauss location
543 
544  // Get parameters
545  const GXmlElement* minor = gammalib::xml_get_par(G_READ, xml, m_semiminor.name());
546  const GXmlElement* major = gammalib::xml_get_par(G_READ, xml, m_semimajor.name());
547 
548  // Read parameters
549  m_semiminor.read(*minor);
550  m_semimajor.read(*major);
551 
552  // Return
553  return;
554 }
555 
556 
557 /***********************************************************************//**
558  * @brief Write model into XML element
559  *
560  * @param[in] xml XML element into which model information is written.
561  *
562  * @exception GException::model_invalid_parnum
563  * Invalid number of model parameters found in XML element.
564  * @exception GException::model_invalid_parnames
565  * Invalid model parameter names found in XML element.
566  *
567  * Write the elliptical gauss model information into an XML element. The XML
568  * element will have the format
569  *
570  * <spatialModel type="EllipticalGaussian">
571  * <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
572  * <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
573  * <parameter name="PA" scale="1.0" value="45.0" min="-360" max="360" free="1"/>
574  * <parameter name="MinorRadius" scale="1.0" value="0.5" min="0.001" max="10" free="1"/>
575  * <parameter name="MajorRadius" scale="1.0" value="2.0" min="0.001" max="10" free="1"/>
576  * </spatialModel>
577  *
578  ***************************************************************************/
580 {
581  // Write disk location
583 
584  // Get or create parameters
587 
588  // Write parameters
589  m_semiminor.write(*minor);
590  m_semimajor.write(*major);
591 
592  // Return
593  return;
594 }
595 
596 
597 /***********************************************************************//**
598  * @brief Print information
599  *
600  * @param[in] chatter Chattiness.
601  * @return String containing model information.
602  ***************************************************************************/
603 std::string GModelSpatialEllipticalGauss::print(const GChatter& chatter) const
604 {
605  // Initialise result string
606  std::string result;
607 
608  // Continue only if chatter is not silent
609  if (chatter != SILENT) {
610 
611  // Append header
612  result.append("=== GModelSpatialEllipticalGauss ===");
613 
614  // Append parameters
615  result.append("\n"+gammalib::parformat("Number of parameters"));
616  result.append(gammalib::str(size()));
617  for (int i = 0; i < size(); ++i) {
618  result.append("\n"+m_pars[i]->print(chatter));
619  }
620 
621  } // endif: chatter was not silent
622 
623  // Return result
624  return result;
625 }
626 
627 
628 /*==========================================================================
629  = =
630  = Private methods =
631  = =
632  ==========================================================================*/
633 
634 /***********************************************************************//**
635  * @brief Initialise class members
636  ***************************************************************************/
638 {
639  // Initialise model type
640  m_type = "EllipticalGaussian";
641 
642  // Initialise precomputation cache. Note that zero values flag
643  // uninitialised as a zero radius is not meaningful
644  m_last_minor = 0.0;
645  m_last_major = 0.0;
646  m_minor_rad = 0.0;
647  m_major_rad = 0.0;
648  m_norm = 0.0;
649  m_last_posangle = 9999.0; // Signals that has not been initialised
650  m_sin2pos = 0.0;
651  m_cospos2 = 0.0;
652  m_sinpos2 = 0.0;
653  m_minor2 = 0.0;
654  m_major2 = 0.0;
655  m_term1 = 0.0;
656  m_term2 = 0.0;
657  m_term3 = 0.0;
658 
659  // Initialise other members
660  m_region.clear();
661 
662  // Return
663  return;
664 }
665 
666 
667 /***********************************************************************//**
668  * @brief Copy class members
669  *
670  * @param[in] model Elliptical Gaussian model.
671  ***************************************************************************/
673 {
674  // Copy members
675  m_type = model.m_type;
676  m_region = model.m_region;
677 
678  // Copy precomputation cache
679  m_last_minor = model.m_last_minor;
680  m_last_major = model.m_last_major;
681  m_minor_rad = model.m_minor_rad;
682  m_major_rad = model.m_major_rad;
683  m_norm = model.m_norm;
685  m_sin2pos = model.m_sin2pos;
686  m_cospos2 = model.m_cospos2;
687  m_sinpos2 = model.m_sinpos2;
688  m_minor2 = model.m_minor2;
689  m_major2 = model.m_major2;
690  m_term1 = model.m_term1;
691  m_term2 = model.m_term2;
692  m_term3 = model.m_term3;
693 
694  // Return
695  return;
696 }
697 
698 
699 /***********************************************************************//**
700  * @brief Delete class members
701  ***************************************************************************/
703 {
704  // Return
705  return;
706 }
707 
708 
709 /***********************************************************************//**
710  * @brief Update precomputation cache
711  *
712  * Computes the normalization
713  * \f[
714  * {\tt m\_norm} = \frac{1}{2 \pi \times a \times b}
715  * \f]
716  * where
717  * \f$a\f$ is the semi-major axis and
718  * \f$b\f$ is the semi-minor axis of the ellipse.
719  *
720  * @warning
721  * The normalization of the elliptical Gaussian is only valid in the small
722  * angle approximation.
723  ***************************************************************************/
725 {
726  // Initialise flag if something has changed
727  bool changed = false;
728 
729  // Update if one axis has changed
730  if (m_last_minor != semiminor() || m_last_major != semimajor()) {
731 
732  // Signal parameter changes
733  changed = true;
734 
735  // Store last values
738 
739  // Compute axes in radians
742 
743  // Perform precomputations. Note that I verified the normalization
744  // by numerical integration of the resulting Gaussian. Note also
745  // also that the normalization is only correct in the small angle
746  // approximation.
749  double denom = gammalib::twopi * m_minor_rad * m_major_rad *
750  c_fraction;
751  m_norm = (denom > 0.0) ? 1.0 / denom : 0.0;
752 
753  } // endif: update required
754 
755  // Update chache if position angle changed
756  #if !defined(G_SMALL_ANGLE_APPROXIMATION)
757  if (m_last_posangle != posangle()) {
758 
759  // Signal parameter changes
760  changed = true;
761 
762  // Store last value
764 
765  // Compute angle in radians
766  double posangle_rad = m_last_posangle * gammalib::deg2rad;
767 
768  // Compute sine and cosine
769  double cospos = std::cos(posangle_rad);
770  double sinpos = std::sin(posangle_rad);
771 
772  // Cache important values for further computations
773  m_cospos2 = cospos * cospos;
774  m_sinpos2 = sinpos * sinpos;
775  m_sin2pos = std::sin(2.0 * posangle_rad);
776 
777  } // endif: position angle update required
778 
779  // Perform precomputations in case anything has changed
780  if (changed) {
781 
782  // Compute help terms
783  m_term1 = 0.5 * (m_cospos2 / m_minor2 + m_sinpos2 / m_major2);
784  m_term2 = 0.5 * (m_cospos2 / m_major2 + m_sinpos2 / m_minor2);
785  m_term3 = 0.5 * (m_sin2pos / m_major2 - m_sin2pos / m_minor2);
786 
787  } // endif: something has changed
788  #endif
789 
790  // Return
791  return;
792 }
793 
794 
795 /***********************************************************************//**
796  * @brief Set boundary sky region
797  ***************************************************************************/
799 {
800  // Set sky region centre to Gaussian centre
802 
803  // Set maximum model radius
804  double max_radius = (semimajor() > semiminor()) ? semimajor() : semiminor();
805 
806  // Set sky region radius to maximum Gaussian sigma times a scaling
807  // factor (actually 3)
808  m_region.radius(max_radius * c_theta_max);
809 
810  // Return
811  return;
812 }
GModelPar m_semimajor
Semi-major axis of ellipse (deg)
double m_last_major
Last semi-major axis.
GSkyDir dir(void) const
Return position of elliptical spatial model.
virtual void write(GXmlElement &xml) const
Write model into XML element.
double m_sinpos2
squared sine of position angle
void free_members(void)
Delete class members.
GSkyRegionCircle m_region
Bounding circle.
virtual GModelSpatialEllipticalGauss * clone(void) const
Clone elliptical Gaussian model.
virtual bool contains(const GSkyDir &dir, const double &margin=0.0) const
Checks whether model contains specified sky direction.
virtual void write(GXmlElement &xml) const
Write model into XML element.
Abstract elliptical spatial model base class.
const std::string & name(void) const
Return parameter name.
GModelPar m_ra
Right Ascension (deg)
GModelSpatialEllipticalGauss(void)
Void constructor.
virtual GModelSpatialEllipticalGauss & operator=(const GModelSpatialEllipticalGauss &model)
Assignment operator.
virtual void clear(void)
Clear elliptical Gaussian model.
int size(void) const
Return number of parameters.
virtual double eval(const double &theta, const double &posangle, 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 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
double m_minor_rad
Minor axis in radians.
virtual GModelSpatialElliptical & operator=(const GModelSpatialElliptical &model)
Assignment operator.
XML element node class.
Definition: GXmlElement.hpp:47
Random number generator class.
Definition: GRan.hpp:44
virtual int elements(void) const
Return number of GXMLElement children of node.
Definition: GXmlNode.cpp:580
GModelPar m_dec
Declination (deg)
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
Interface definition for the spatial model registry class.
double m_cospos2
squared cosine of position angle
void init_members(void)
Initialise class members.
void set_region(void) const
Set boundary sky region.
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)
Class that handles photons.
Definition: GPhoton.hpp:47
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
void update(void) const
Update precomputation cache.
double semiminor(void) const
Return semi-minor axis of ellipse.
const double deg2rad
Definition: GMath.hpp:43
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
double m_last_minor
Last semi-minor axis.
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
double m_sin2pos
sine of twice the position angle
Spatial model registry class definition.
const GTime & time(void) const
Return photon time.
Definition: GPhoton.hpp:134
GChatter
Definition: GTypemaps.hpp:33
void init_members(void)
Initialise class members.
virtual void read(const GXmlElement &xml)
Read model from XML element.
GModelPar m_semiminor
Semi-minor axis of ellipse (deg)
virtual GSkyDir mc(const GEnergy &energy, const GTime &time, GRan &ran) const
Returns MC sky direction.
const GModelSpatialEllipticalGauss g_elliptical_gauss_seed
double normal(void)
Returns normal deviates.
Definition: GRan.cpp:257
virtual std::string type(void) const
Return model type.
const GSkyDir & centre(void) const
Return circular region centre.
void free_members(void)
Delete class members.
double atan2d(const double &y, const double &x)
Compute arc tangens in degrees.
Definition: GMath.cpp:308
double posangle(void) const
Return Position Angle of model.
void read(const GXmlElement &xml)
Extract parameter attributes from XML element.
Definition: GModelPar.cpp:229
double value(void) const
Return parameter value.
double m_last_posangle
Last position angle.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print information.
void copy_members(const GModelSpatialEllipticalGauss &model)
Copy class members.
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1226
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
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
Elliptical gauss model class interface definition.
void init_members(void)
Initialise class members.
const double twopi
Definition: GMath.hpp:36
virtual ~GModelSpatialEllipticalGauss(void)
Destructor.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
void free_members(void)
Delete class members.
Sky direction class.
Definition: GSkyDir.hpp:62
double m_major_rad
Major axis in radians.
double semimajor(void) const
Return semi-major axis of ellipse.
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
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
virtual double theta_max(void) const
Return maximum model radius (in radians)