GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpatialEllipticalGeneralGauss.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpatialRadialGeneralGauss.cpp - Generalised radial Gaussian *
3  * source model class *
4  * ----------------------------------------------------------------------- *
5  * copyright (C) 2022 by Luigi Tibaldo *
6  * ----------------------------------------------------------------------- *
7  * *
8  * This program is free software: you can redistribute it and/or modify *
9  * it under the terms of the GNU General Public License as published by *
10  * the Free Software Foundation, either version 3 of the License, or *
11  * (at your option) any later version. *
12  * *
13  * This program is distributed in the hope that it will be useful, *
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
16  * GNU General Public License for more details. *
17  * *
18  * You should have received a copy of the GNU General Public License *
19  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
20  * *
21  ***************************************************************************/
22 /**
23  * @file GModelSpatialEllipticalGeneralGauss.cpp
24  * @brief Generalised elliptical Gaussian model class implementation
25  * @author Luigi Tibaldo
26  */
27 
28 /* __ Includes ___________________________________________________________ */
29 #ifdef HAVE_CONFIG_H
30 #include <config.h>
31 #endif
32 #include "GException.hpp"
33 #include "GTools.hpp"
34 #include "GMath.hpp"
38 
39 /* __ Constants __________________________________________________________ */
40 namespace {
41  const double c_theta_max = 5.0; //!< max distance from center in units of major semiaxis
42 }
43 
44 /* __ Globals ____________________________________________________________ */
46 const GModelSpatialRegistry g_elliptical_general_gauss_registry(&g_elliptical_general_gauss_seed);
47 
48 /* __ Method name definitions ____________________________________________ */
49 #define G_CONSTRUCTOR "GModelSpatialEllipticalGeneralGauss::"\
50  "GModelSpatialEllipticalGeneralGauss(GSkyDir&, double&, double&, "\
51  "double&, double&, std::string&)"
52 #define G_READ "GModelSpatialEllipticalGeneralGauss::read(GXmlElement&)"
53 #define G_WRITE "GModelSpatialEllipticalGeneralGauss::write(GXmlElement&)"
54 #define G_MC "GModelSpatialEllipticalGeneralGauss::mc(GEnergy&, GTime&, "\
55  "GRan&)"
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 Elliptical Gaussian constructor
87  *
88  * @param[in] dir Centre of elliptical Gaussian.
89  * @param[in] semimajor Semi-major axis (deg).
90  * @param[in] semiminor Semi-minor axis (deg).
91  * @param[in] posangle Position angle of semi-major axis (deg).
92  * @param[in] ridx Reciprocal of radial profile index.
93  * @param[in] coordsys Coordinate system (either "CEL" or "GAL")
94  *
95  * @exception GException::invalid_argument
96  * Invalid @p coordsys argument specified.
97  *
98  * Construct elliptical Gaussian model from the ellipse centre direction
99  * (@p dir), the @p semimajor and @p semiminor axes, the position
100  * angle (@p posangle), and the reciprocal of the radial profile index
101  * (@p ridx). The @p coordsys parameter specifies whether the sky direction
102  * should be interpreted in the celestial or Galactic coordinate system.
103  ***************************************************************************/
105  const double& semimajor,
106  const double& semiminor,
107  const double& posangle,
108  const double& ridx,
109  const std::string& coordsys) :
111 {
112  // Throw an exception if the coordinate system is invalid
113  if ((coordsys != "CEL") && (coordsys != "GAL")) {
114  std::string msg = "Invalid coordinate system \""+coordsys+"\" "
115  "specified. Please specify either \"CEL\" or "
116  "\"GAL\".";
118  }
119 
120  // Initialise members
121  init_members();
122 
123  // Set parameter names
124  if (coordsys == "CEL") {
125  m_lon.name("RA");
126  m_lat.name("DEC");
127  }
128  else {
129  m_lon.name("GLON");
130  m_lat.name("GLAT");
131  }
132 
133  // Assign parameters
134  this->dir(dir);
135  this->semiminor(semiminor);
136  this->semimajor(semimajor);
137  this->posangle(posangle);
138  this->ridx(ridx);
139 
140  // Return
141  return;
142 }
143 
144 
145 /***********************************************************************//**
146  * @brief XML constructor
147  *
148  * @param[in] xml XML element.
149  *
150  * Constructs generalised elliptical Gaussian model by extracting information
151  * from an XML element. See the read() method for more information about the
152  * expected structure of the XML element.
153  ***************************************************************************/
156 {
157  // Initialise members
158  init_members();
159 
160  // Read information from XML element
161  read(xml);
162 
163  // Return
164  return;
165 }
166 
167 
168 /***********************************************************************//**
169  * @brief Copy constructor
170  *
171  * @param[in] model Generalised elliptical Gaussian model.
172  ***************************************************************************/
175 {
176  // Initialise members
177  init_members();
178 
179  // Copy members
180  copy_members(model);
181 
182  // Return
183  return;
184 }
185 
186 
187 /***********************************************************************//**
188  * @brief Destructor
189  ***************************************************************************/
191 {
192  // Free members
193  free_members();
194 
195  // Return
196  return;
197 }
198 
199 
200 /*==========================================================================
201  = =
202  = Operators =
203  = =
204  ==========================================================================*/
205 
206 /***********************************************************************//**
207  * @brief Assignment operator
208  *
209  * @param[in] model Generalised elliptical Gaussian model.
210  * @return Generalised elliptical Gaussian model.
211  ***************************************************************************/
213 {
214  // Execute only if object is not identical
215  if (this != &model) {
216 
217  // Copy base class members
219 
220  // Free members
221  free_members();
222 
223  // Initialise members
224  init_members();
225 
226  // Copy members
227  copy_members(model);
228 
229  } // endif: object was not identical
230 
231  // Return
232  return *this;
233 }
234 
235 
236 /*==========================================================================
237  = =
238  = Public methods =
239  = =
240  ==========================================================================*/
241 
242 /***********************************************************************//**
243  * @brief Clear generalised elliptical Gaussian model
244  ***************************************************************************/
246 {
247  // Free class members (base and derived classes, derived class first)
248  free_members();
251 
252  // Initialise members
255  init_members();
256 
257  // Return
258  return;
259 }
260 
261 
262 /***********************************************************************//**
263  * @brief Clone generalised elliptical Gaussian model
264  *
265  * @return Pointer to deep copy of generalised elliptical Gaussian model.
266  ***************************************************************************/
268 {
269  // Clone generalised elliptical Gaussian model
270  return new GModelSpatialEllipticalGeneralGauss(*this);
271 }
272 
273 
274 /***********************************************************************//**
275  * @brief Evaluate function (in units of sr^-1)
276  *
277  * @param[in] theta Angular distance from disk centre (radians).
278  * @param[in] posangle Position angle (counterclockwise from North) (radians).
279  * @param[in] energy Photon energy.
280  * @param[in] time Photon arrival time.
281  * @param[in] gradients Compute gradients?
282  * @return Model value.
283  *
284  * Evaluates the spatial component for a generalised elliptical Gaussian
285  * source model.
286  *
287  * The generalised elliptical Gaussian function is defined by
288  * \f$S_{\rm p}(\theta, \phi | E, t)\f$, where
289  * \f$\theta\f$ is the angular separation between elliptical Gaussian centre
290  * and the actual location and \f$\phi\f$ the position angle with respect to
291  * the model centre, counted counterclockwise from North.
292  *
293  * The function \f$S_{\rm p}(\theta, \phi | E, t)\f$ is given by
294  *
295  * \f[
296  * S_{\rm p}(\theta, \phi | E, t) = {\tt m\_norm} \times
297  * \exp \left( -\left(\frac{\theta}{r_{\rm eff}}\right)^1/\eta \right)
298  * \f]
299  *
300  * where the effective ellipse radius \f$r_{\rm eff}\f$ towards a given
301  * position angle is given by
302  *
303  * \f[
304  * r_{\rm eff} = \frac{ab}
305  * {\sqrt{\left( a \sin (\phi - \phi_0) \right)^2} +
306  * \sqrt{\left( b \cos (\phi - \phi_0) \right)^2}}
307  * \f]
308  * and
309  * \f$a\f$ is the semi-major axis of the ellipse,
310  * \f$b\f$ is the semi-minor axis,
311  * \f$\phi_0\f$ is the position angle of the ellipse, counted
312  * counterclockwise from North, and
313  * \f$\et\f$ is the reciprocal of the radial profile index
314  * \f${\tt m\_norm}\f$ is a normalization constant given by
315  *
316  * \f[
317  * {\tt m\_norm} = \frac{1}{2 \pi a b \eta \Gamma(\eta)}
318  * \f]
319  *
320  * The method will not compute analytical parameter gradients, even if the
321  * @p gradients argument is set to true. Parameter gradients
322  * need to be computed numerically.
323  *
324  * @warning
325  * The above formula for an elliptical generalised Gaussian are accurate for small
326  * angles, with semimajor and semiminor axes below a few degrees. This
327  * should be fine for most use cases, but you should be aware that the
328  * ellipse will get distorted for larger angles.
329  * Note that the model normalisation is only correct in the small angle
330  * approximation and for \f$\eta\f$ order unity or smaller.
331  ***************************************************************************/
332 double GModelSpatialEllipticalGeneralGauss::eval(const double& theta,
333  const double& posangle,
334  const GEnergy& energy,
335  const GTime& time,
336  const bool& gradients) const
337 {
338  // Initialise value
339  double value = 0.0;
340 
341  // Continue only if we're inside circle enclosing the ellipse
342  if (theta <= theta_max()) {
343 
344  // Update precomputation cache
345  update();
346 
347  // Perform computations to compute exponent
348  #if defined(G_SMALL_ANGLE_APPROXIMATION)
349  double rel_posangle = posangle - this->posangle() * gammalib::deg2rad;
350  double cosinus = std::cos(rel_posangle);
351  double sinus = std::sin(rel_posangle);
352  double arg1 = m_minor_rad * cosinus;
353  double arg2 = m_major_rad * sinus;
354  double r_ellipse = m_minor_rad * m_major_rad /
355  std::sqrt(arg1*arg1 + arg2*arg2); //!< small angle
356  double r_relative = theta/r_ellipse;
357  double exponent = std::pow(r_relative,m_inv_ridx);
358  #else
359  // Perform computations
360  double sinphi = std::sin(posangle);
361  double cosphi = std::cos(posangle);
362 
363  // Compute help term
364  double term1 = m_term1 * cosphi * cosphi;
365  double term2 = m_term2 * sinphi * sinphi;
366  double term3 = m_term3 * sinphi * cosphi;
367 
368  // Compute exponent
369  double exponent = std::pow(theta * std::sqrt(term1 + term2 + term3), m_inv_ridx) ;
370  #endif
371 
372  // Set value
373  value = m_norm * std::exp(-exponent);
374 
375  // Compile option: Check for NaN/Inf
376  #if defined(G_NAN_CHECK)
377  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
378  std::cout << "*** ERROR: GModelSpatialEllipticalGeneralGauss::eval";
379  std::cout << "(theta=" << theta << "): NaN/Inf encountered";
380  std::cout << "(posangle=" << posangle << "): NaN/Inf encountered";
381  std::cout << " (value=" << value;
382  std::cout << ", MinorAxis^2=" << m_minor2;
383  std::cout << ", MaxjorAxis^2=" << m_major2;
384  std::cout << ", m_norm=" << m_norm;
385  std::cout << ")" << std::endl;
386  }
387  #endif
388 
389  } // endif: position was inside enclosing circle
390 
391  // Return value
392  return value;
393 }
394 
395 
396 /***********************************************************************//**
397  * @brief Returns MC sky direction
398  *
399  * @param[in] energy Photon energy.
400  * @param[in] time Photon arrival time.
401  * @param[in,out] ran Random number generator.
402  * @return Sky direction.
403  *
404  * Draws an arbitrary sky direction from the elliptical generalised
405  * Gaussian model.
406  ***************************************************************************/
408  const GTime& time,
409  GRan& ran) const
410 {
411  // Update precomputation cache
412  update();
413 
414  // Initialise simulation
415  #if defined(G_DEBUG_MC)
416  int n_samples = 0;
417  #endif
418  double u_max = m_norm * std::sin(theta_max());
419  double value = 0.0;
420  double u = 1.0;
421  double theta = 0.0;
422  double phi = 0.0;
423 
424  // Throw an exception if the maximum value is zero
425  if (u_max <= 0.0) {
426  std::string msg = "Non positive maximum ellpitical general Gauss model value.";
428  }
429 
430  // Simulate offset from photon arrival direction using a rejection method
431  do {
432  theta = ran.uniform() * theta_max();
433  phi = ran.uniform() * gammalib::twopi;
434  value = eval(theta, phi, energy, time) * std::sin(theta);
435  u = ran.uniform() * u_max;
436  #if defined(G_DEBUG_MC)
437  n_samples++;
438  #endif
439  } while (u > value);
440  #if defined(G_DEBUG_MC)
441  std::cout << "#=" << n_samples << " ";
442  #endif
443 
444  // Rotate pointing direction by offset and azimuth angle
445  GSkyDir sky_dir = dir();
446  sky_dir.rotate_deg(phi * gammalib::rad2deg, theta * gammalib::rad2deg);
447 
448  // Return sky direction
449  return sky_dir;
450 }
451 
452 
453 /***********************************************************************//**
454  * @brief Checks whether model contains specified sky direction
455  *
456  * @param[in] dir Sky direction.
457  * @param[in] margin Margin to be added to sky direction (deg).
458  * @return True if the model contains the sky direction.
459  *
460  * Signals whether a sky direction is contained in the elliptical gauss
461  * model.
462  *
463  * @todo Implement correct evaluation of effective ellipse radius.
464  ***************************************************************************/
466  const double& margin) const
467 {
468  // Compute distance to centre (radian)
469  double distance = dir.dist(this->dir());
470 
471  // Return flag
472  return (distance <= theta_max() + margin*gammalib::deg2rad);
473 }
474 
475 
476 /***********************************************************************//**
477  * @brief Return maximum model radius (in radians)
478  *
479  * @return Returns maximum model radius (radians).
480  *
481  * Returns the maximum of 3.0 semimajor() and 3.0 semiminor() as
482  * approximate edge of the Gaussian. This limit is of course arbitrary, but
483  * allows to limit the integration region for response computation. The value
484  * of 3.0 has been determined by experiment (#1561).
485  ***************************************************************************/
487 {
488  // Set maximum model radius
489  double theta_max = (semimajor() > semiminor())
490  ? semimajor() * gammalib::deg2rad * c_theta_max
491  : semiminor() * gammalib::deg2rad * c_theta_max;
492 
493  // Return value
494  return theta_max;
495 }
496 
497 
498 /***********************************************************************//**
499  * @brief Read model from XML element
500  *
501  * @param[in] xml XML element.
502  *
503  * Reads the elliptical gauss model information from an XML element. The XML
504  * element shall have either the format
505  *
506  * <spatialModel type="EllipticalGeneralGaussian">
507  * <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
508  * <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
509  * <parameter name="PA" scale="1.0" value="45.0" min="-360" max="360" free="1"/>
510  * <parameter name="MinorRadius" scale="1.0" value="0.5" min="0.001" max="10" free="1"/>
511  * <parameter name="MajorRadius" scale="1.0" value="2.0" min="0.001" max="10" free="1"/>
512  * <parameter name="R_Index" scale="1.0" value="0.5" min="0.01" max="10" free="1"/>
513  * </spatialModel>
514  *
515  * or
516  *
517  * <spatialModel type="EllipticalGeneralGaussian">
518  * <parameter name="GLON" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
519  * <parameter name="GLAT" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
520  * <parameter name="PA" scale="1.0" value="45.0" min="-360" max="360" free="1"/>
521  * <parameter name="MinorRadius" scale="1.0" value="0.5" min="0.001" max="10" free="1"/>
522  * <parameter name="MajorRadius" scale="1.0" value="2.0" min="0.001" max="10" free="1"/>
523  * <parameter name="R_Index" scale="1.0" value="0.5" min="0.01" max="10" free="1"/>
524  * </spatialModel>
525  *
526  * @todo Implement a test of the ellipse boundary. The axes
527  * and axes minimum should be >0.
528  ***************************************************************************/
530 {
531  // Verify number of model parameters
533 
534  // Read gauss location
536 
537  // Get parameters
538  const GXmlElement* minor = gammalib::xml_get_par(G_READ, xml, m_semiminor.name());
539  const GXmlElement* major = gammalib::xml_get_par(G_READ, xml, m_semimajor.name());
541 
542  // Read parameters
543  m_semiminor.read(*minor);
544  m_semimajor.read(*major);
545  m_ridx.read(*ridx);
546 
547  // Return
548  return;
549 }
550 
551 
552 /***********************************************************************//**
553  * @brief Write model into XML element
554  *
555  * @param[in] xml XML element into which model information is written.
556  *
557  * Write the elliptical gauss model information into an XML element. The XML
558  * element will have the format
559  *
560  * <spatialModel type="EllipticalGeneralGaussian">
561  * <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
562  * <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
563  * <parameter name="PA" scale="1.0" value="45.0" min="-360" max="360" free="1"/>
564  * <parameter name="MinorRadius" scale="1.0" value="0.5" min="0.001" max="10" free="1"/>
565  * <parameter name="MajorRadius" scale="1.0" value="2.0" min="0.001" max="10" free="1"/>
566  * <parameter name="R_Index" scale="1.0" value="0.5" min="0.01" max="10" free="1"/>
567  * </spatialModel>
568  *
569  ***************************************************************************/
571 {
572  // Verify model type
574 
575  // Write disk location
577 
578  // Get or create parameters
582 
583  // Write parameters
584  m_semiminor.write(*minor);
585  m_semimajor.write(*major);
586  m_ridx.write(*ridx);
587 
588  // Return
589  return;
590 }
591 
592 
593 /***********************************************************************//**
594  * @brief Print information
595  *
596  * @param[in] chatter Chattiness.
597  * @return String containing model information.
598  ***************************************************************************/
599 std::string GModelSpatialEllipticalGeneralGauss::print(const GChatter& chatter) const
600 {
601  // Initialise result string
602  std::string result;
603 
604  // Continue only if chatter is not silent
605  if (chatter != SILENT) {
606 
607  // Append header
608  result.append("=== GModelSpatialEllipticalGeneralGauss ===");
609 
610  // Append parameters
611  result.append("\n"+gammalib::parformat("Number of parameters"));
612  result.append(gammalib::str(size()));
613  for (int i = 0; i < size(); ++i) {
614  result.append("\n"+m_pars[i]->print(chatter));
615  }
616 
617  } // endif: chatter was not silent
618 
619  // Return result
620  return result;
621 }
622 
623 
624 /*==========================================================================
625  = =
626  = Private methods =
627  = =
628  ==========================================================================*/
629 
630 /***********************************************************************//**
631  * @brief Initialise class members
632  ***************************************************************************/
634 {
635  // Initialise model type
636  m_type = "EllipticalGeneralGaussian";
637 
638  // Initialise Reciprocal index
639  m_ridx.clear();
640  m_ridx.name("R_Index");
641  m_ridx.value(0.5); // Gaussian
642  m_ridx.min(1.e-15); // need to be > 0
643  m_ridx.free();
644  m_ridx.scale(1.0);
645  m_ridx.gradient(0.0);
646  m_ridx.has_grad(false);
647 
648  // Set parameter pointer(s)
649  m_pars.push_back(&m_ridx);
650 
651  // Initialise precomputation cache. Note that zero values flag
652  // uninitialised as a zero radius is not meaningful
653  m_last_minor = 0.0;
654  m_last_major = 0.0;
655  m_minor_rad = 0.0;
656  m_major_rad = 0.0;
657  m_norm = 0.0;
658  m_last_posangle = 9999.0; // Signals that has not been initialised
659  m_sin2pos = 0.0;
660  m_cospos2 = 0.0;
661  m_sinpos2 = 0.0;
662  m_minor2 = 0.0;
663  m_major2 = 0.0;
664  m_term1 = 0.0;
665  m_term2 = 0.0;
666  m_term3 = 0.0;
667  m_last_ridx = 0.0;
668  m_inv_ridx = 0.0;
669 
670  // Return
671  return;
672 }
673 
674 
675 /***********************************************************************//**
676  * @brief Copy class members
677  *
678  * @param[in] model Elliptical Gaussian model.
679  ***************************************************************************/
681 {
682  // Copy members
683  m_type = model.m_type; // Needed to conserve model type
684  m_ridx = model.m_ridx;
685 
686  // Copy precomputation cache
687  m_last_minor = model.m_last_minor;
688  m_last_major = model.m_last_major;
689  m_minor_rad = model.m_minor_rad;
690  m_major_rad = model.m_major_rad;
691  m_norm = model.m_norm;
693  m_sin2pos = model.m_sin2pos;
694  m_cospos2 = model.m_cospos2;
695  m_sinpos2 = model.m_sinpos2;
696  m_minor2 = model.m_minor2;
697  m_major2 = model.m_major2;
698  m_term1 = model.m_term1;
699  m_term2 = model.m_term2;
700  m_term3 = model.m_term3;
701  m_last_ridx = model.m_last_ridx;
702  m_inv_ridx = model.m_inv_ridx;
703 
704  // Return
705  return;
706 }
707 
708 
709 /***********************************************************************//**
710  * @brief Delete class members
711  ***************************************************************************/
713 {
714  // Return
715  return;
716 }
717 
718 
719 /***********************************************************************//**
720  * @brief Update precomputation cache
721  *
722  * Computes the normalization
723  * \f[
724  * {\tt m\_norm} = \frac{1}{2 \pi \times a \times b}
725  * \f]
726  * where
727  * \f$a\f$ is the semi-major axis and
728  * \f$b\f$ is the semi-minor axis of the ellipse.
729  *
730  * @warning
731  * The normalization of the elliptical Gaussian is only valid in the small
732  * angle approximation.
733  ***************************************************************************/
735 {
736  // Initialise flag if something has changed
737  bool changed = false;
738 
739  // Update if one axis has changed or radial profile index has changed
740  if (m_last_minor != semiminor() || m_last_major != semimajor() || m_last_ridx != ridx()) {
741 
742  // Signal parameter changes
743  changed = true;
744 
745  // Store last values
748  m_last_ridx = ridx();
749 
750  // Compute axes in radians
753 
754  // Perform precomputations. Note also
755  // that the model normalisation is only correct in the small angle
756  // approximation and for ridx order unity or smaller.
757  if ((m_minor_rad > 0.0) && (m_major_rad > 0.0) && (ridx() > 0.0)) {
758  m_inv_ridx = 1.0 / ridx() ;
760  std::exp(gammalib::gammln(2.0 * ridx())) );
761  }
762  else {
763  m_inv_ridx = 1.0;
764  m_norm = 0.0;
765  }
766 
767  } // endif: update required
768 
769  // Update chache if position angle changed
770  #if !defined(G_SMALL_ANGLE_APPROXIMATION)
771  if (m_last_posangle != posangle()) {
772 
773  // Signal parameter changes
774  changed = true;
775 
776  // Store last value
778 
779  // Compute angle in radians
780  double posangle_rad = m_last_posangle * gammalib::deg2rad;
781 
782  // Compute sine and cosine
783  double cospos = std::cos(posangle_rad);
784  double sinpos = std::sin(posangle_rad);
785 
786  // Cache important values for further computations
787  m_cospos2 = cospos * cospos;
788  m_sinpos2 = sinpos * sinpos;
789  m_sin2pos = std::sin(2.0 * posangle_rad);
790 
791  } // endif: position angle update required
792 
793  // Perform precomputations in case anything has changed
794  if (changed) {
795 
796  // Compute help terms
797  m_term1 = 0.5 * (m_cospos2 / m_minor2 + m_sinpos2 / m_major2);
798  m_term2 = 0.5 * (m_cospos2 / m_major2 + m_sinpos2 / m_minor2);
799  m_term3 = 0.5 * (m_sin2pos / m_major2 - m_sin2pos / m_minor2);
800 
801  } // endif: something has changed
802  #endif
803 
804  // Return
805  return;
806 }
807 
808 
809 /***********************************************************************//**
810  * @brief Set boundary sky region
811  ***************************************************************************/
813 {
814  // Set maximum model radius
815  double max_radius = (semimajor() > semiminor()) ? semimajor() : semiminor();
816 
817  // Set sky region circle (maximum Gaussian sigma times a scaling
818  // factor (actually 5))
819  GSkyRegionCircle region(dir(), max_radius * c_theta_max);
820 
821  // Set region (circumvent const correctness)
822  const_cast<GModelSpatialEllipticalGeneralGauss*>(this)->m_region = region;
823 
824  // Return
825  return;
826 }
GModelPar m_semimajor
Semi-major axis of ellipse (deg)
virtual void write(GXmlElement &xml) const
Write model into XML element.
void free_members(void)
Delete class members.
Abstract elliptical spatial model base class.
const std::string & name(void) const
Return parameter name.
virtual void clear(void)
Clear generalised elliptical Gaussian model.
GModelPar m_lon
Right Ascension or Galactic longitude (deg)
double gradient(void) const
Return parameter gradient.
int size(void) const
Return number of parameters.
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
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.
virtual GModelSpatialElliptical & operator=(const GModelSpatialElliptical &model)
Assignment operator.
virtual GSkyDir mc(const GEnergy &energy, const GTime &time, GRan &ran) const
Returns MC sky direction.
XML element node class.
Definition: GXmlElement.hpp:48
const GSkyDir & dir(void) const
Return position of elliptical spatial model.
Random number generator class.
Definition: GRan.hpp:44
Time class.
Definition: GTime.hpp:55
Gammalib tools definition.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print information.
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.
virtual void set_region(void) const
Set boundary sky region.
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
double semiminor(void) const
Return semi-minor axis of ellipse.
const GModelSpatialEllipticalGeneralGauss g_elliptical_general_gauss_seed
const double deg2rad
Definition: GMath.hpp:43
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1358
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
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
GModelPar m_ridx
Reciprocal of exponent of the radial profile.
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
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.
virtual double theta_max(void) const
Return maximum model radius (in radians)
GChatter
Definition: GTypemaps.hpp:33
virtual GModelSpatialEllipticalGeneralGauss & operator=(const GModelSpatialEllipticalGeneralGauss &model)
Assignment operator.
virtual void read(const GXmlElement &xml)
Read model from XML element.
GModelPar m_semiminor
Semi-minor axis of ellipse (deg)
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
double gammln(const double &arg)
Computes logarithm of gamma function.
Definition: GMath.cpp:383
void update(void) const
Update precomputation cache.
const GSkyRegion * region(void) const
Return boundary sky region.
void free_members(void)
Delete class members.
void copy_members(const GModelSpatialEllipticalGeneralGauss &model)
Copy class members.
double ridx(void) const
Return reciprocal of the elliptical profile index.
std::string m_type
Spatial model type.
double posangle(void) const
Return Position Angle of model.
Generalised elliptical gaussian model 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.
GModelPar m_lat
Declination or Galactic latitude (deg)
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1422
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1316
double m_sin2pos
sine of twice the position angle
Exception handler interface definition.
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
Definition: GSkyDir.hpp:271
double m_cospos2
squared cosine of position angle
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1232
void init_members(void)
Initialise class members.
const double twopi
Definition: GMath.hpp:36
const double rad2deg
Definition: GMath.hpp:44
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
virtual bool contains(const GSkyDir &dir, const double &margin=0.0) const
Checks whether model contains specified sky direction.
Sky direction class.
Definition: GSkyDir.hpp:62
virtual GModelSpatialEllipticalGeneralGauss * clone(void) const
Clone generalised elliptical Gaussian model.
double semimajor(void) const
Return semi-major axis of ellipse.
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
virtual void write(GXmlElement &xml) const
Write model into XML element.
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