Justine Devin / gammalib

Commit b1be222c8b2e3fe5f0df46031cb3a0e5cea7218a

Authored by Jürgen Knödlseder
2 parents 01fb6dd4 7d33bb30

Merge branch '1561-fix-ellipse-gauss' into integration

Inline Side-by-side
Showing 2 changed files with 95 additions and 56 deletions
inst/cta/src/GCTAResponse_helpers.cpp
  View file @b1be222
... ... @@ -864,11 +864,15 @@ double cta_irf_elliptical_kern_omega::eval(const double& omega)
864 864 double arg2 = m_semimajor_rad * sinus;
865 865 double r_ellipse = m_semiminor_rad * m_semimajor_rad /
866 866 std::sqrt(arg1*arg1 + arg2*arg2);
867   - std::cout << "*** WARNING: cta_irf_elliptical_kern_omega::eval";
  867 + std::cout << "*** WARNING: cta_irf_elliptical_kern_omega::eval:";
868 868 std::cout << " zero model for (rho,omega)=(";
869 869 std::cout << m_rho*gammalib::rad2deg << ",";
870 870 std::cout << omega*gammalib::rad2deg << ")";
871   - std::cout << " rho-r_ellipse=" << (m_rho-r_ellipse) << " radians";
  871 + std::cout << " semiminor=" << m_semiminor_rad;
  872 + std::cout << " semimajor=" << m_semimajor_rad;
  873 + std::cout << " posang=" << m_model->posangle() * gammalib::deg2rad;
  874 + std::cout << " rel_posang=" << diff_angle * gammalib::deg2rad;
  875 + std::cout << " rho/r_ellipse=" << (m_rho/r_ellipse);
872 876 std::cout << std::endl;
873 877 }
874 878 #endif
... ...
src/model/GModelSpatialEllipticalGauss.cpp
  View file @b1be222
... ... @@ -36,6 +36,11 @@
36 36 #include "GModelSpatialRegistry.hpp"
37 37  
38 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 +}
39 44  
40 45 /* __ Globals ____________________________________________________________ */
41 46 const GModelSpatialEllipticalGauss g_elliptical_gauss_seed;
... ... @@ -48,6 +53,7 @@ const GModelSpatialRegistry g_elliptical_gauss_registry(&amp;g_elliptical_gauss_seed
48 53 /* __ Macros _____________________________________________________________ */
49 54  
50 55 /* __ Coding definitions _________________________________________________ */
  56 +#define G_SMALL_ANGLE_APPROXIMATION //!< Use small angle approximation
51 57  
52 58 /* __ Debug definitions __________________________________________________ */
53 59  
... ... @@ -73,14 +79,14 @@ GModelSpatialEllipticalGauss::GModelSpatialEllipticalGauss(void) :
73 79  
74 80  
75 81 /***********************************************************************//**
76   - * @brief Elliptical Gauss constructor
  82 + * @brief Elliptical Gaussian constructor
77 83 *
78   - * @param[in] dir Sky position of gauss centre.
  84 + * @param[in] dir Centre of elliptical Gaussian.
79 85 * @param[in] semimajor Semi-major axis (degrees).
80 86 * @param[in] semiminor Semi-minor axis (degrees).
81 87 * @param[in] posangle Position angle of semi-major axis (degrees).
82 88 *
83   - * Construct elliptical gauss model from sky position of the ellipse centre
  89 + * Construct elliptical Gaussian model from the ellipse centre direction
84 90 * (@p dir), the @p semimajor and @p semiminor axes, and the position
85 91 * angle (@p posangle).
86 92 ***************************************************************************/
... ... @@ -109,7 +115,7 @@ GModelSpatialEllipticalGauss::GModelSpatialEllipticalGauss(const GSkyDir&amp; dir,
109 115 *
110 116 * @param[in] xml XML element.
111 117 *
112   - * Constructs elliptical gauss model by extracting information from an XML
  118 + * Constructs elliptical Gaussian model by extracting information from an XML
113 119 * element. See the read() method for more information about the expected
114 120 * structure of the XML element.
115 121 ***************************************************************************/
... ... @@ -130,7 +136,7 @@ GModelSpatialEllipticalGauss::GModelSpatialEllipticalGauss(const GXmlElement&amp; xm
130 136 /***********************************************************************//**
131 137 * @brief Copy constructor
132 138 *
133   - * @param[in] model Elliptical gauss model.
  139 + * @param[in] model Elliptical Gaussian model.
134 140 ***************************************************************************/
135 141 GModelSpatialEllipticalGauss::GModelSpatialEllipticalGauss(const GModelSpatialEllipticalGauss& model) :
136 142 GModelSpatialElliptical(model)
... ... @@ -169,7 +175,7 @@ GModelSpatialEllipticalGauss::~GModelSpatialEllipticalGauss(void)
169 175 * @brief Assignment operator
170 176 *
171 177 * @param[in] model Elliptical gauss model.
172   - * @return Elliptical gauss model.
  178 + * @return Elliptical Gaussian model.
173 179 ***************************************************************************/
174 180 GModelSpatialEllipticalGauss& GModelSpatialEllipticalGauss::operator=(const GModelSpatialEllipticalGauss& model)
175 181 {
... ... @@ -202,7 +208,7 @@ GModelSpatialEllipticalGauss&amp; GModelSpatialEllipticalGauss::operator=(const GMod
202 208 ==========================================================================*/
203 209  
204 210 /***********************************************************************//**
205   - * @brief Clear elliptical gauss model
  211 + * @brief Clear elliptical Gaussian model
206 212 ***************************************************************************/
207 213 void GModelSpatialEllipticalGauss::clear(void)
208 214 {
... ... @@ -222,13 +228,13 @@ void GModelSpatialEllipticalGauss::clear(void)
222 228  
223 229  
224 230 /***********************************************************************//**
225   - * @brief Clone elliptical gauss model
  231 + * @brief Clone elliptical Gaussian model
226 232 *
227   - * @return Pointer to deep copy of elliptical gauss model.
  233 + * @return Pointer to deep copy of elliptical Gaussian model.
228 234 ***************************************************************************/
229 235 GModelSpatialEllipticalGauss* GModelSpatialEllipticalGauss::clone(void) const
230 236 {
231   - // Clone elliptical gauss model
  237 + // Clone elliptical Gaussian model
232 238 return new GModelSpatialEllipticalGauss(*this);
233 239 }
234 240  
... ... @@ -242,47 +248,50 @@ GModelSpatialEllipticalGauss* GModelSpatialEllipticalGauss::clone(void) const
242 248 * @param[in] time Photon arrival time.
243 249 * @return Model value.
244 250 *
245   - * Evaluates the spatial component for an elliptical gauss source model. The
246   - * source model is an elliptical gaussian function
  251 + * Evaluates the spatial component for an elliptical Gaussian source model.
  252 + *
  253 + * The elliptical Gaussian function is defined by
247 254 * \f$S_{\rm p}(\theta, \phi | E, t)\f$, where
248   - * \f$\theta\f$ is the angular separation between elliptical gauss centre and
249   - * the actual location and \f$\phi\f$ the position angle with respect to the
250   - * model centre, counted counterclockwise from North.
  255 + * \f$\theta\f$ is the angular separation between elliptical Gaussian centre
  256 + * and the actual location and \f$\phi\f$ the position angle with respect to
  257 + * the model centre, counted counterclockwise from North.
251 258 *
252   - * The function \f$f(\theta, \phi)\f$ is given by
  259 + * The function \f$S_{\rm p}(\theta, \phi | E, t)\f$ is given by
253 260 *
254 261 * \f[
255   - * S_{\rm p}(\theta, \phi | E, t) = \left \{
256   - * \begin{array}{l l}
257   - * {\tt m\_norm}
258   - * & \mbox{if} \, \, \theta \le \theta_0 \\
259   - * \\
260   - * 0 & \mbox{else}
261   - * \end{array}
262   - * \right .
  262 + * S_{\rm p}(\theta, \phi | E, t) = {\tt m\_norm} \times
  263 + * \exp \left( -\frac{\theta^2}{2 r_{rm eff}} \right)
263 264 * \f]
264 265 *
265   - * where \f$\theta_0\f$ is the effective radius of the ellipse on the sphere
266   - * given by
  266 + * where the effective ellipse radius \f$r_{rm eff}\f$ towards a given
  267 + * position angle is given by
267 268 *
268   - * \f[\theta_0\ =
269   - * \frac{ab}{\sqrt{b^2 \cos^2(\phi-\phi_0) + a^2 \sin^2(\phi-\phi_0)}}\f]
270   - *
271   - * and
  269 + * \f[
  270 + * r_{rm eff} = \frac{ab}
  271 + * {\sqrt{\left( a \sin (\phi - \phi_0) \right)^2 +
  272 + * \sqrt{\left( b \cos (\phi - \phi_0) \right)^2}
  273 + * \f]
  274 + * and
272 275 * \f$a\f$ is the semi-major axis of the ellipse,
273 276 * \f$b\f$ is the semi-minor axis, and
274 277 * \f$\phi_0\f$ is the position angle of the ellipse, counted
275 278 * counterclockwise from North.
  279 + * \f${\tt m\_norm}\f$ is a normalization constant given by
276 280 *
277   - * The normalisation constant \f${\tt m\_norm}\f$ which is the inverse of the
278   - * solid angle subtended by an ellipse is given by
279 281 * \f[
280 282 * {\tt m\_norm} = \frac{1}{2 \pi \times a \times b}
281 283 * \f]
282 284 *
283 285 * @warning
284   - * The normalization of the elliptical Gaussian is only valid in the small
285   - * angle approximation.
  286 + * The above formula for an elliptical Gaussian are accurate for small
  287 + * angles, with semimajor and semiminor axes below a few degrees. This
  288 + * should be fine for most use cases, but you should be aware that the
  289 + * ellipse will get distorted for larger angles.
  290 + *
  291 + * @warning
  292 + * For numerical reasons the elliptical Gaussian will be truncated for
  293 + * \f$\theta\f$ angles that correspond to 3.0 times the effective ellipse
  294 + * radius.
286 295 ***************************************************************************/
287 296 double GModelSpatialEllipticalGauss::eval(const double& theta,
288 297 const double& posangle,
... ... @@ -298,6 +307,18 @@ double GModelSpatialEllipticalGauss::eval(const double&amp; theta,
298 307 // Update precomputation cache
299 308 update();
300 309  
  310 + // Perform computations to compute exponent
  311 + #if defined(G_SMALL_ANGLE_APPROXIMATION)
  312 + double rel_posangle = posangle - this->posangle() * gammalib::deg2rad;
  313 + double cosinus = std::cos(rel_posangle);
  314 + double sinus = std::sin(rel_posangle);
  315 + double arg1 = m_minor_rad * cosinus;
  316 + double arg2 = m_major_rad * sinus;
  317 + double r_ellipse = m_minor_rad * m_major_rad /
  318 + std::sqrt(arg1*arg1 + arg2*arg2); //!< small angle
  319 + double r_relative = theta/r_ellipse;
  320 + double exponent = 0.5*r_relative*r_relative;
  321 + #else
301 322 // Perform computations
302 323 double sinphi = std::sin(posangle);
303 324 double cosphi = std::cos(posangle);
... ... @@ -308,10 +329,13 @@ double GModelSpatialEllipticalGauss::eval(const double&amp; theta,
308 329 double term3 = m_term3 * sinphi * cosphi;
309 330  
310 331 // Compute exponent
311   - double exponent = theta * theta * (term1 + term2 + term3);
  332 + double exponent = theta * theta * (term1 + term2 + term3);
  333 + #endif
312 334  
313   - // Set value
314   - value = m_norm * std::exp(-exponent);
  335 + // Set value if the exponent is within the truncation region
  336 + if (exponent <= c_max_exponent) {
  337 + value = m_norm * std::exp(-exponent);
  338 + }
315 339  
316 340 // Compile option: Check for NaN/Inf
317 341 #if defined(G_NAN_CHECK)
... ... @@ -344,7 +368,7 @@ double GModelSpatialEllipticalGauss::eval(const double&amp; theta,
344 368 * @return Model value.
345 369 *
346 370 * Evaluates the function value. No gradient computation is implemented as
347   - * Elliptical models will be convolved with the instrument response and thus
  371 + * elliptical models will be convolved with the instrument response and thus
348 372 * require the numerical computation of the derivatives.
349 373 *
350 374 * See the eval() method for more information.
... ... @@ -367,7 +391,12 @@ double GModelSpatialEllipticalGauss::eval_gradients(const double&amp; theta,
367 391 * @param[in,out] ran Random number generator.
368 392 * @return Sky direction.
369 393 *
370   - * Draws an arbitrary sky position from the 2D gauss distribution.
  394 + * Draws an arbitrary sky direction from the elliptical Gaussian model.
  395 + *
  396 + * @warning
  397 + * For numerical reasons the elliptical Gaussian will be truncated for
  398 + * \f$\theta\f$ angles that correspond to 3.0 times the effective ellipse
  399 + * radius.
371 400 ***************************************************************************/
372 401 GSkyDir GModelSpatialEllipticalGauss::mc(const GEnergy& energy,
373 402 const GTime& time,
... ... @@ -382,10 +411,18 @@ GSkyDir GModelSpatialEllipticalGauss::mc(const GEnergy&amp; energy,
382 411 photon.time(time);
383 412  
384 413 // Draw gaussian offset from each axis
385   - double theta1 = semimajor() * ran.normal();
386   - double theta2 = semiminor() * ran.normal();
387   -
388   - // Compute total offset from model centre
  414 + double ran_major;
  415 + double ran_minor;
  416 + do {
  417 + ran_major = ran.normal();
  418 + } while (ran_major > c_theta_max);
  419 + do {
  420 + ran_minor = ran.normal();
  421 + } while (ran_minor > c_theta_max);
  422 + double theta1 = semimajor() * ran_major;
  423 + double theta2 = semiminor() * ran_minor;
  424 +
  425 + // Compute total offset from model centre in small angle approximation
389 426 double theta = std::sqrt(theta1 * theta1 + theta2 * theta2);
390 427  
391 428 // Compute rotation angle, taking into account given position angle
... ... @@ -431,19 +468,17 @@ bool GModelSpatialEllipticalGauss::contains(const GSkyDir&amp; dir,
431 468 *
432 469 * @return Returns maximum model radius.
433 470 *
434   - * Returns the maximum of \f$2.5\f$ semimajor() and \f$2.5\f$ semiminor() as
  471 + * Returns the maximum of 3.0 semimajor() and 3.0 semiminor() as
435 472 * approximate edge of the Gaussian. This limit is of course arbitrary, but
436 473 * allows to limit the integration region for response computation. The value
437   - * of 2.5 has been determined by experiment. Ideally, a large value should
438   - * be used to properly take into account the tails of the distribution, but
439   - * this apparently leads to fit convergence problems (#1561).
  474 + * of 3.0 has been determined by experiment (#1561).
440 475 ***************************************************************************/
441 476 double GModelSpatialEllipticalGauss::theta_max(void) const
442 477 {
443 478 // Set maximum model radius
444 479 double theta_max = (semimajor() > semiminor())
445   - ? semimajor() * gammalib::deg2rad * 2.5
446   - : semiminor() * gammalib::deg2rad * 2.5;
  480 + ? semimajor() * gammalib::deg2rad * c_theta_max
  481 + : semiminor() * gammalib::deg2rad * c_theta_max;
447 482  
448 483 // Return value
449 484 return theta_max;
... ... @@ -690,7 +725,7 @@ void GModelSpatialEllipticalGauss::init_members(void)
690 725 /***********************************************************************//**
691 726 * @brief Copy class members
692 727 *
693   - * @param[in] model Elliptical gauss model.
  728 + * @param[in] model Elliptical Gaussian model.
694 729 ***************************************************************************/
695 730 void GModelSpatialEllipticalGauss::copy_members(const GModelSpatialEllipticalGauss& model)
696 731 {
... ... @@ -739,9 +774,6 @@ void GModelSpatialEllipticalGauss::free_members(void)
739 774 * @warning
740 775 * The normalization of the elliptical Gaussian is only valid in the small
741 776 * angle approximation.
742   - *
743   - * @todo Generalize normalization to the sphere (for the moment only accurate
744   - * in the small angle approximation).
745 777 ***************************************************************************/
746 778 void GModelSpatialEllipticalGauss::update() const
747 779 {
... ... @@ -768,12 +800,14 @@ void GModelSpatialEllipticalGauss::update() const
768 800 // approximation.
769 801 m_minor2 = m_minor_rad * m_minor_rad;
770 802 m_major2 = m_major_rad * m_major_rad;
771   - double denom = gammalib::twopi * m_minor_rad * m_major_rad;
  803 + double denom = gammalib::twopi * m_minor_rad * m_major_rad *
  804 + c_fraction;
772 805 m_norm = (denom > 0.0) ? 1.0 / denom : 0.0;
773 806  
774 807 } // endif: update required
775 808  
776 809 // Update chache if position angle changed
  810 + #if !defined(G_SMALL_ANGLE_APPROXIMATION)
777 811 if (m_last_posangle != posangle()) {
778 812  
779 813 // Signal parameter changes
... ... @@ -805,6 +839,7 @@ void GModelSpatialEllipticalGauss::update() const
805 839 m_term3 = 0.5 * (m_sin2pos / m_major2 - m_sin2pos / m_minor2);
806 840  
807 841 } // endif: something has changed
  842 + #endif
808 843  
809 844 // Return
810 845 return;
... ...