Jürgen Knödlseder / gammalib

Commit e9681ff1b66b2dc273d4591200d75c12f8a002aa

Authored by Jürgen Knödlseder
2 parents b8a93dd1 4e43a6db

Merge branch '3924-GeneralGauss' into integration 

# Conflicts:
#	ChangeLog
#	NEWS
#	README.md
Inline Side-by-side
Showing 18 changed files with 2239 additions and 38 deletions
ChangeLog
  View file @e9681ff
1   -2021-12-13
  1 +2022-01-10
2 2  
3 3 * Version 2.0.0 released
4 4 ========================
5 5  
6 6 Add GCOMObservation::tim() and GCOMObservation::oads() setter methods
  7 + Add GModelSpatialRadialGeneralGauss and
  8 + GModelSpatialEllipticalGeneralGauss classes (#3924)
7 9 Remove GCOMSelection::zeta_min() and GCOMSelection::zeta_max()
8 10 methods
9   - Add GCOMSelection::fpmtflag(), GCOMSelection::use_d1() and
  11 + Add GCOMSelection::fpmtflag(), GCOMSelection::use_d1() and
10 12 GCOMSelection::use_d2() methods (#3917)
11 13 Add GCOMModelDRBPhibarBins class (#3913)
12 14 Compute spectral model gradients for negative factor values
... ...
1 1 New Features and Important Changes in GammaLib 2.0.0
2 2  
3   -13 December 2021
  3 +10 January 2022
4 4  
5 5  
6 6 1. Introduction
... ... @@ -20,6 +20,8 @@ The following classes have been added:
20 20 - GModelAssociation
21 21 - GModelAssociations
22 22 - GModelSpectralBins
  23 +- GModelSpatialRadialGeneralGauss
  24 +- GModelSpatialEllipticalGeneralGauss
23 25 - GCTAModelSkyCube
24 26 - GCTABackground2D
25 27 - GCOMModelDRM
... ... @@ -207,6 +209,9 @@ Implemented GMatrixSparse::row() and GMatrixSparse::add_to_row() methods.
207 209  
208 210 9. Model module
209 211 ---------------
  212 +Added GModelSpatialRadialGeneralGauss and GModelSpatialEllipticalGeneralGauss
  213 +classes to model general Gaussian functions as spatial components (#3924).
  214 +
210 215 Gradients for spectral model parameters are now also computed for all models in
211 216 case that the factor values are negative. This provides now error bars for photon
212 217 and energy flux power laws that are allowed to take negative flux values.
... ...
README.md
  View file @e9681ff
1 1 GammaLib information
2 2 ====================
3   -* Version: 2.0.0.dev (13 December 2021)
  3 +* Version: 2.0.0.dev (10 January 2022)
4 4  
5 5 [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=gammalib-integrate-os)](https://cta-jenkins.irap.omp.eu/job/gammalib-integrate-os/)
6 6  
... ...
doc/source/admin/release_history/2.0.rst
  View file @e9681ff
1 1 .. _2.0:
2 2  
3   -GammaLib 2.0 release (3 December 2021)
  3 +GammaLib 2.0 release (10 January 2022)
4 4 ======================================
5 5  
6 6 Overview
... ... @@ -96,6 +96,8 @@ Improvements
96 96 New features
97 97 ------------
98 98  
  99 +* [`3924 <https://cta-redmine.irap.omp.eu/issues/3924>`_] -
  100 + Add ``GModelSpatialRadialGeneralGauss`` and ``GModelSpatialEllipticalGeneralGauss`` classes
99 101 * [`3910 <https://cta-redmine.irap.omp.eu/issues/3910>`_] -
100 102 Add ``GCOMModelDRM`` class
101 103 * [`3606 <https://cta-redmine.irap.omp.eu/issues/3606>`_] -
... ...
include/GModelSpatialEllipticalGeneralGauss.hpp 0 → 100644
  View file @e9681ff
  1 +/***************************************************************************
  2 + * GModelSpatialEllipticalGeneralGauss.hpp - Generalised elliptical *
  3 + * Gaussian 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.hpp
  24 + * @brief Generalised elliptical gaussian model class interface definition
  25 + * @author Luigi Tibaldo
  26 + */
  27 +
  28 +#ifndef GMODELSPATIALELLIPTICALGENERALGAUSS_HPP
  29 +#define GMODELSPATIALELLIPTICALGENERALGAUSS_HPP
  30 +
  31 +/* __ Includes ___________________________________________________________ */
  32 +#include <string>
  33 +#include "GModelSpatialElliptical.hpp"
  34 +#include "GModelPar.hpp"
  35 +#include "GSkyDir.hpp"
  36 +#include "GSkyRegionCircle.hpp"
  37 +#include "GXmlElement.hpp"
  38 +
  39 +
  40 +/**************************************************************************
  41 + * @class GModelSpatialEllipticalGeneralGauss
  42 + *
  43 + * @brief Generalised elliptical gaussian source model class
  44 + *
  45 + * This class implements the spatial component of the factorised source
  46 + * model for an generalised elliptical Gaussian source, i.e. surface
  47 + * brightness according to an asymmetric generalised Gaussian.
  48 + ***************************************************************************/
  49 +class GModelSpatialEllipticalGeneralGauss : public GModelSpatialElliptical {
  50 +
  51 +public:
  52 + // Constructors and destructors
  53 + GModelSpatialEllipticalGeneralGauss(void);
  54 + GModelSpatialEllipticalGeneralGauss(const GSkyDir& dir,
  55 + const double& major,
  56 + const double& minor,
  57 + const double& posangle,
  58 + const double& ridx);
  59 + explicit GModelSpatialEllipticalGeneralGauss(const GXmlElement& xml);
  60 + GModelSpatialEllipticalGeneralGauss(const GModelSpatialEllipticalGeneralGauss& model);
  61 + virtual ~GModelSpatialEllipticalGeneralGauss(void);
  62 +
  63 + // Operators
  64 + virtual GModelSpatialEllipticalGeneralGauss& operator=(const GModelSpatialEllipticalGeneralGauss& model);
  65 +
  66 + // Implemented pure virtual base class methods
  67 + virtual void clear(void);
  68 + virtual GModelSpatialEllipticalGeneralGauss* clone(void) const;
  69 + virtual std::string classname(void) const;
  70 + virtual double eval(const double& theta,
  71 + const double& posangle,
  72 + const GEnergy& energy,
  73 + const GTime& time,
  74 + const bool& gradients = false) const;
  75 + virtual GSkyDir mc(const GEnergy& energy,
  76 + const GTime& time,
  77 + GRan& ran) const;
  78 + virtual bool contains(const GSkyDir& dir,
  79 + const double& margin = 0.0) const;
  80 + virtual double theta_max(void) const;
  81 + virtual void read(const GXmlElement& xml);
  82 + virtual void write(GXmlElement& xml) const;
  83 + virtual std::string print(const GChatter& chatter = NORMAL) const;
  84 +
  85 + // Other methods
  86 + double ridx(void) const;
  87 + void ridx(const double& ridx);
  88 +
  89 +
  90 +protected:
  91 + // Protected methods
  92 + void init_members(void);
  93 + void copy_members(const GModelSpatialEllipticalGeneralGauss& model);
  94 + void free_members(void);
  95 + void update(void) const;
  96 + virtual void set_region(void) const;
  97 +
  98 + // Protected members
  99 + GModelPar m_ridx; //!< Reciprocal of exponent of the radial profile
  100 +
  101 + // Cached members used for pre-computations
  102 + mutable double m_last_minor; //!< Last semi-minor axis
  103 + mutable double m_last_major; //!< Last semi-major axis
  104 + mutable double m_last_posangle; //!< Last position angle
  105 + mutable double m_cospos2; //!< squared cosine of position angle
  106 + mutable double m_sinpos2; //!< squared sine of position angle
  107 + mutable double m_sin2pos; //!< sine of twice the position angle
  108 + mutable double m_minor2; //!< square of minor axis
  109 + mutable double m_major2; //!< square of major axis
  110 + mutable double m_minor_rad; //!< Minor axis in radians
  111 + mutable double m_major_rad; //!< Major axis in radians
  112 + mutable double m_norm; //!< Normalization
  113 + mutable double m_term1; //!< Help term 1
  114 + mutable double m_term2; //!< Help term 2
  115 + mutable double m_term3; //!< Help term 3
  116 + mutable double m_last_ridx; //!< Last reciprocal radial index
  117 + mutable double m_inv_ridx; //!< Spatial profile index
  118 +
  119 +};
  120 +
  121 +
  122 +/***********************************************************************//**
  123 + * @brief Return class name
  124 + *
  125 + * @return String containing the class name ("GModelSpatialEllipticalGeneralGauss").
  126 + ***************************************************************************/
  127 +inline
  128 +std::string GModelSpatialEllipticalGeneralGauss::classname(void) const
  129 +{
  130 + return ("GModelSpatialEllipticalGeneralGauss");
  131 +}
  132 +
  133 +/***********************************************************************//**
  134 + * @brief Return reciprocal of the elliptical profile index
  135 + *
  136 + * @return Reciprocal of the elliptical profile index.
  137 + *
  138 + * Returns the reciprocal of the elliptical profile index.
  139 + ***************************************************************************/
  140 +inline
  141 +double GModelSpatialEllipticalGeneralGauss::ridx(void) const
  142 +{
  143 + return (m_ridx.value());
  144 +}
  145 +
  146 +/***********************************************************************//**
  147 + * @brief Set reciprocal of the elliptical profile index
  148 + *
  149 + * @param[in] ridx Reciprocal of the elliptical profile index.
  150 + *
  151 + * Sets the reciprocal of the elliptical profile index.
  152 + ***************************************************************************/
  153 +inline
  154 +void GModelSpatialEllipticalGeneralGauss::ridx(const double& ridx)
  155 +{
  156 + m_ridx.value(ridx);
  157 + return;
  158 +}
  159 +
  160 +#endif /* GMODELSPATIALELLIPTICALGENERALGAUSS_HPP */
... ...
include/GModelSpatialRadialGeneralGauss.hpp 0 → 100644
  View file @e9681ff
  1 +/***************************************************************************
  2 + * GModelSpatialRadialGeneralGauss.hpp - Generalised radial Gaussian *
  3 + * source model class *
  4 + * ----------------------------------------------------------------------- *
  5 + * copyright (C) 2021-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 GModelSpatialRadialGeneralGauss.hpp
  24 + * @brief Generalized radial Gaussian model class interface definition
  25 + * @author Luigi Tibaldo
  26 + */
  27 +
  28 +#ifndef GMODELSPATIALRADIALGENERALGAUSS_HPP
  29 +#define GMODELSPATIALRADIALGENERALGAUSS_HPP
  30 +
  31 +/* __ Includes ___________________________________________________________ */
  32 +#include <string>
  33 +#include "GModelSpatialRadial.hpp"
  34 +#include "GModelPar.hpp"
  35 +#include "GSkyDir.hpp"
  36 +#include "GSkyRegionCircle.hpp"
  37 +#include "GEnergy.hpp"
  38 +#include "GTime.hpp"
  39 +#include "GXmlElement.hpp"
  40 +
  41 +
  42 +/***********************************************************************//**
  43 + * @class GModelSpatialRadialGeneralGauss
  44 + *
  45 + * @brief Generalized radial Gaussian model class
  46 + *
  47 + * This class implements the spatial component of the factorised source
  48 + * model for a generalised radial Gaussian source.
  49 + ***************************************************************************/
  50 +class GModelSpatialRadialGeneralGauss : public GModelSpatialRadial {
  51 +
  52 +public:
  53 + // Constructors and destructors
  54 + GModelSpatialRadialGeneralGauss(void);
  55 + GModelSpatialRadialGeneralGauss(const GSkyDir& dir,
  56 + const double& radius,
  57 + const double& ridx);
  58 + explicit GModelSpatialRadialGeneralGauss(const GXmlElement& xml);
  59 + GModelSpatialRadialGeneralGauss(const GModelSpatialRadialGeneralGauss& model);
  60 + virtual ~GModelSpatialRadialGeneralGauss(void);
  61 +
  62 + // Operators
  63 + virtual GModelSpatialRadialGeneralGauss& operator=(const GModelSpatialRadialGeneralGauss& model);
  64 +
  65 + // Implemented pure virtual methods
  66 + virtual void clear(void);
  67 + virtual GModelSpatialRadialGeneralGauss* clone(void) const;
  68 + virtual std::string classname(void) const;
  69 + virtual double eval(const double& theta,
  70 + const GEnergy& energy,
  71 + const GTime& time,
  72 + const bool& gradients = false) const;
  73 + virtual GSkyDir mc(const GEnergy& energy,
  74 + const GTime& time,
  75 + GRan& ran) const;
  76 + virtual bool contains(const GSkyDir& dir,
  77 + const double& margin = 0.0) const;
  78 + virtual double theta_max(void) const;
  79 + virtual void read(const GXmlElement& xml);
  80 + virtual void write(GXmlElement& xml) const;
  81 + virtual std::string print(const GChatter& chatter = NORMAL) const;
  82 +
  83 + // Other methods
  84 + double radius(void) const;
  85 + void radius(const double& radius);
  86 + double ridx(void) const;
  87 + void ridx(const double& ridx);
  88 +
  89 +protected:
  90 + // Protected methods
  91 + void init_members(void);
  92 + void copy_members(const GModelSpatialRadialGeneralGauss& model);
  93 + void free_members(void);
  94 + void update(void) const;
  95 + virtual void set_region(void) const;
  96 +
  97 + // Protected members
  98 + GModelPar m_radius; //!< Gaussian width (deg)
  99 + GModelPar m_ridx; //!< Reciprocal of exponent of the radial profile
  100 +
  101 + // Cached members used for pre-computations
  102 + mutable double m_last_radius; //!< Last radius
  103 + mutable double m_inv_radius_rad; //!< radius(rad)^-1
  104 + mutable double m_last_ridx; //!< Last reciprocal radial index
  105 + mutable double m_inv_ridx; //!< Spatial profile index
  106 + mutable double m_value_norm; //!< 1/(2pi radius(rad)^2 ridx Gamma(ridx))
  107 +};
  108 +
  109 +
  110 +/***********************************************************************//**
  111 + * @brief Return class name
  112 + *
  113 + * @return String containing the class name ("GModelSpatialRadialGeneralGauss").
  114 + ***************************************************************************/
  115 +inline
  116 +std::string GModelSpatialRadialGeneralGauss::classname(void) const
  117 +{
  118 + return ("GModelSpatialRadialGeneralGauss");
  119 +}
  120 +
  121 +
  122 +/***********************************************************************//**
  123 + * @brief Return radius
  124 + *
  125 + * @return Radius (deg).
  126 + *
  127 + * Returns the radius in degrees.
  128 + ***************************************************************************/
  129 +inline
  130 +double GModelSpatialRadialGeneralGauss::radius(void) const
  131 +{
  132 + return (m_radius.value());
  133 +}
  134 +
  135 +/***********************************************************************//**
  136 + * @brief Set radius
  137 + *
  138 + * @param[in] radius (deg).
  139 + *
  140 + * Sets the radius in degrees.
  141 + ***************************************************************************/
  142 +inline
  143 +void GModelSpatialRadialGeneralGauss::radius(const double& radius)
  144 +{
  145 + m_radius.value(radius);
  146 + return;
  147 +}
  148 +
  149 +/***********************************************************************//**
  150 + * @brief Return ridx
  151 + *
  152 + * @return Reciprocal of the radial profile index.
  153 + *
  154 + * Returns the reciprocal of the radial profile index.
  155 + ***************************************************************************/
  156 +inline
  157 +double GModelSpatialRadialGeneralGauss::ridx(void) const
  158 +{
  159 + return (m_ridx.value());
  160 +}
  161 +
  162 +/***********************************************************************//**
  163 + * @brief Set reciprocal index
  164 + *
  165 + * @param[in] ridx Reciprocal of the radial profile index.
  166 + *
  167 + * Sets the reciprocal index of the radial profile.
  168 + ***************************************************************************/
  169 +inline
  170 +void GModelSpatialRadialGeneralGauss::ridx(const double& ridx)
  171 +{
  172 + m_ridx.value(ridx);
  173 + return;
  174 +}
  175 +
  176 +#endif /* GMODELSPATIALRADIALGENERALGAUSS_HPP */
... ...
include/GammaLib.hpp
  View file @e9681ff
... ... @@ -213,6 +213,7 @@
213 213 #include "GModelSpatialRadial.hpp"
214 214 #include "GModelSpatialRadialDisk.hpp"
215 215 #include "GModelSpatialRadialGauss.hpp"
  216 +#include "GModelSpatialRadialGeneralGauss.hpp"
216 217 #include "GModelSpatialRadialShell.hpp"
217 218 #include "GModelSpatialRadialRing.hpp"
218 219 #include "GModelSpatialRadialProfile.hpp"
... ... @@ -223,6 +224,7 @@
223 224 #include "GModelSpatialElliptical.hpp"
224 225 #include "GModelSpatialEllipticalDisk.hpp"
225 226 #include "GModelSpatialEllipticalGauss.hpp"
  227 +#include "GModelSpatialEllipticalGeneralGauss.hpp"
226 228 #include "GModelSpatialDiffuse.hpp"
227 229 #include "GModelSpatialDiffuseConst.hpp"
228 230 #include "GModelSpatialDiffuseCube.hpp"
... ...
include/Makefile.am
  View file @e9681ff
... ... @@ -171,6 +171,7 @@ pkginclude_HEADERS = GammaLib.hpp \
171 171 GModelSpatialRadial.hpp \
172 172 GModelSpatialRadialDisk.hpp \
173 173 GModelSpatialRadialGauss.hpp \
  174 + GModelSpatialRadialGeneralGauss.hpp \
174 175 GModelSpatialRadialShell.hpp \
175 176 GModelSpatialRadialRing.hpp \
176 177 GModelSpatialRadialProfile.hpp \
... ... @@ -181,6 +182,7 @@ pkginclude_HEADERS = GammaLib.hpp \
181 182 GModelSpatialElliptical.hpp \
182 183 GModelSpatialEllipticalDisk.hpp \
183 184 GModelSpatialEllipticalGauss.hpp \
  185 + GModelSpatialEllipticalGeneralGauss.hpp \
184 186 GModelSpatialDiffuse.hpp \
185 187 GModelSpatialDiffuseConst.hpp \
186 188 GModelSpatialDiffuseCube.hpp \
... ...
pyext/GModelSpatialEllipticalGeneralGauss.i 0 → 100644
  View file @e9681ff
  1 +/***************************************************************************
  2 + * GModelSpatialEllipticalGeneralGauss.i - Generalised elliptical *
  3 + * Gaussian 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.i
  24 + * @brief Elliptical generalised gauss model class interface definition
  25 + * @author Luigi Tibaldo
  26 + */
  27 +%{
  28 +/* Put headers and other declarations here that are needed for compilation */
  29 +#include "GModelSpatialEllipticalGeneralGauss.hpp"
  30 +%}
  31 +
  32 +
  33 +/***********************************************************************//**
  34 + * @class GModelSpatialEllipticalGeneralGauss
  35 + *
  36 + * @brief Elliptical generalised gauss model class
  37 + ***************************************************************************/
  38 +class GModelSpatialEllipticalGeneralGauss : public GModelSpatialElliptical {
  39 +
  40 +public:
  41 + // Constructors and destructors
  42 + GModelSpatialEllipticalGeneralGauss(void);
  43 + GModelSpatialEllipticalGeneralGauss(const GSkyDir& dir,
  44 + const double& major,
  45 + const double& minor,
  46 + const double& posangle,
  47 + const double& ridx);
  48 + explicit GModelSpatialEllipticalGeneralGauss(const GXmlElement& xml);
  49 + GModelSpatialEllipticalGeneralGauss(const GModelSpatialEllipticalGeneralGauss& model);
  50 + virtual ~GModelSpatialEllipticalGeneralGauss(void);
  51 +
  52 + // Implemented pure virtual base class methods
  53 + virtual void clear(void);
  54 + virtual GModelSpatialEllipticalGeneralGauss* clone(void) const;
  55 + virtual std::string classname(void) const;
  56 + virtual double eval(const double& theta,
  57 + const double& posangle,
  58 + const GEnergy& energy,
  59 + const GTime& time,
  60 + const bool& gradients = false) const;
  61 + virtual GSkyDir mc(const GEnergy& energy,
  62 + const GTime& time,
  63 + GRan& ran) const;
  64 + virtual bool contains(const GSkyDir& dir,
  65 + const double& margin = 0.0) const;
  66 + virtual double theta_max(void) const;
  67 + virtual void read(const GXmlElement& xml);
  68 + virtual void write(GXmlElement& xml) const;
  69 +
  70 + // Other methods
  71 + double ridx(void) const;
  72 + void ridx(const double& ridx);
  73 +};
  74 +
  75 +
  76 +/***********************************************************************//**
  77 + * @brief GModelSpatialEllipticalGeneralGauss class extension
  78 + ***************************************************************************/
  79 +%extend GModelSpatialEllipticalGeneralGauss {
  80 + GModelSpatialEllipticalGeneralGauss copy() {
  81 + return (*self);
  82 + }
  83 + double eval(const GPhoton& photon) const {
  84 + return self->GModelSpatialEllipticalGeneralGauss::eval(photon);
  85 + }
  86 + double eval(const GPhoton& photon, const bool& gradients) const {
  87 + return self->GModelSpatialEllipticalGeneralGauss::eval(photon, gradients);
  88 + }
  89 +%pythoncode {
  90 + def __getstate__(self):
  91 + xml = gammalib.GXmlElement()
  92 + self.write(xml)
  93 + state = (xml,)
  94 + return state
  95 + def __setstate__(self, state):
  96 + self.__init__(state[0])
  97 +}
  98 +};
... ...
pyext/GModelSpatialRadialGeneralGauss.i 0 → 100644
  View file @e9681ff
  1 +/***************************************************************************
  2 + * GModelSpatialRadialGeneralGauss.i - Generalised radial Gaussian *
  3 + * source model class *
  4 + * ----------------------------------------------------------------------- *
  5 + * copyright (C) 2021-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 GModelSpatialRadialGeneralGauss.i
  24 + * @brief Generalised radial Gaussian model class Python interface definition
  25 + * @author Luigi Tibaldo
  26 + */
  27 +%{
  28 +/* Put headers and other declarations here that are needed for compilation */
  29 +#include "GModelSpatialRadialGeneralGauss.hpp"
  30 +%}
  31 +
  32 +
  33 +/***********************************************************************//**
  34 + * @class GModelSpatialRadialGeneralGauss
  35 + *
  36 + * @brief Generalised radial Gaussian model class
  37 + ***************************************************************************/
  38 +class GModelSpatialRadialGeneralGauss : public GModelSpatialRadial {
  39 +
  40 +public:
  41 + // Constructors and destructors
  42 + GModelSpatialRadialGeneralGauss(void);
  43 + GModelSpatialRadialGeneralGauss(const GSkyDir& dir,
  44 + const double& radius,
  45 + const double& ridx);
  46 + explicit GModelSpatialRadialGeneralGauss(const GXmlElement& xml);
  47 + GModelSpatialRadialGeneralGauss(const GModelSpatialRadialGeneralGauss& model);
  48 + virtual ~GModelSpatialRadialGeneralGauss(void);
  49 +
  50 + // Implemented pure virtual methods
  51 + virtual void clear(void);
  52 + virtual GModelSpatialRadialGeneralGauss* clone(void) const;
  53 + virtual std::string classname(void) const;
  54 + virtual double eval(const double& theta,
  55 + const GEnergy& energy,
  56 + const GTime& time,
  57 + const bool& gradients = false) const;
  58 + virtual GSkyDir mc(const GEnergy& energy,
  59 + const GTime& time,
  60 + GRan& ran) const;
  61 + virtual bool contains(const GSkyDir& dir,
  62 + const double& margin = 0.0) const;
  63 + virtual double theta_max(void) const;
  64 + virtual void read(const GXmlElement& xml);
  65 + virtual void write(GXmlElement& xml) const;
  66 +
  67 + // Other methods
  68 + double radius(void) const;
  69 + double ridx(void) const;
  70 + void radius(const double& radius);
  71 + void ridx(const double& ridx);
  72 +};
  73 +
  74 +
  75 +/***********************************************************************//**
  76 + * @brief GModelSpatialRadialGeneralGauss class extension
  77 + *
  78 + * The eval(GPhoton&, bool&) method need to be defined in the extension to
  79 + * force swig to build also the interface for these methods that are
  80 + * implemented in the base class only. It's not clear to me why these
  81 + * methods are not inherited automatically. Maybe this could also be handled
  82 + * by a %typemap(typecheck) construct.
  83 + ***************************************************************************/
  84 +%extend GModelSpatialRadialGeneralGauss {
  85 + GModelSpatialRadialGeneralGauss copy() {
  86 + return (*self);
  87 + }
  88 + double eval(const GPhoton& photon) const {
  89 + return self->GModelSpatialRadialGeneralGauss::eval(photon);
  90 + }
  91 + double eval(const GPhoton& photon, const bool& gradients) const {
  92 + return self->GModelSpatialRadialGeneralGauss::eval(photon, gradients);
  93 + }
  94 +%pythoncode {
  95 + def __getstate__(self):
  96 + xml = gammalib.GXmlElement()
  97 + self.write(xml)
  98 + state = (xml,)
  99 + return state
  100 + def __setstate__(self, state):
  101 + self.__init__(state[0])
  102 +}
  103 +};
... ...
pyext/gammalib/model.i
  View file @e9681ff
... ... @@ -180,6 +180,7 @@
180 180 %include "GModelSpatialRadialDisk.i"
181 181 %include "GModelSpatialRadialRing.i"
182 182 %include "GModelSpatialRadialGauss.i"
  183 +%include "GModelSpatialRadialGeneralGauss.i"
183 184 %include "GModelSpatialRadialShell.i"
184 185 %include "GModelSpatialRadialProfile.i"
185 186 %include "GModelSpatialRadialProfileGauss.i"
... ... @@ -189,6 +190,7 @@
189 190 %include "GModelSpatialElliptical.i"
190 191 %include "GModelSpatialEllipticalDisk.i"
191 192 %include "GModelSpatialEllipticalGauss.i"
  193 +%include "GModelSpatialEllipticalGeneralGauss.i"
192 194 %include "GModelSpatialDiffuse.i"
193 195 %include "GModelSpatialDiffuseConst.i"
194 196 %include "GModelSpatialDiffuseCube.i"
... ...
src/model/GModelSpatialEllipticalGeneralGauss.cpp 0 → 100644
  View file @e9681ff
  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"
  35 +#include "GModelSpatialEllipticalGeneralGauss.hpp"
  36 +#include "GModelSpatialRadialDisk.hpp"
  37 +#include "GModelSpatialRegistry.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 ____________________________________________________________ */
  45 +const GModelSpatialEllipticalGeneralGauss g_elliptical_general_gauss_seed;
  46 +const GModelSpatialRegistry g_elliptical_general_gauss_registry(&g_elliptical_general_gauss_seed);
  47 +
  48 +/* __ Method name definitions ____________________________________________ */
  49 +#define G_READ "GModelSpatialEllipticalGeneralGauss::read(GXmlElement&)"
  50 +#define G_WRITE "GModelSpatialEllipticalGeneralGauss::write(GXmlElement&)"
  51 +#define G_MC "GModelSpatialEllipticalGeneralGauss::mc(GEnergy&, GTime&, "\
  52 + "GRan&)"
  53 +
  54 +/* __ Macros _____________________________________________________________ */
  55 +
  56 +/* __ Coding definitions _________________________________________________ */
  57 +#define G_SMALL_ANGLE_APPROXIMATION //!< Use small angle approximation
  58 +
  59 +/* __ Debug definitions __________________________________________________ */
  60 +
  61 +
  62 +/*==========================================================================
  63 + = =
  64 + = Constructors/destructors =
  65 + = =
  66 + ==========================================================================*/
  67 +
  68 +/***********************************************************************//**
  69 + * @brief Void constructor
  70 + ***************************************************************************/
  71 +GModelSpatialEllipticalGeneralGauss::GModelSpatialEllipticalGeneralGauss(void) :
  72 + GModelSpatialElliptical()
  73 +{
  74 + // Initialise members
  75 + init_members();
  76 +
  77 + // Return
  78 + return;
  79 +}
  80 +
  81 +
  82 +/***********************************************************************//**
  83 + * @brief Elliptical Gaussian constructor
  84 + *
  85 + * @param[in] dir Centre of elliptical Gaussian.
  86 + * @param[in] semimajor Semi-major axis (deg).
  87 + * @param[in] semiminor Semi-minor axis (deg).
  88 + * @param[in] posangle Position angle of semi-major axis (deg).
  89 + * @param[in] ridx Reciprocal of radial profile index.
  90 + *
  91 + * Construct elliptical Gaussian model from the ellipse centre direction
  92 + * (@p dir), the @p semimajor and @p semiminor axes, the position
  93 + * angle (@p posangle), and the reciprocal of the radial profile index
  94 + * (@p ridx).
  95 + ***************************************************************************/
  96 +GModelSpatialEllipticalGeneralGauss::GModelSpatialEllipticalGeneralGauss(const GSkyDir& dir,
  97 + const double& semimajor,
  98 + const double& semiminor,
  99 + const double& posangle,
  100 + const double& ridx) :
  101 + GModelSpatialElliptical()
  102 +{
  103 + // Initialise members
  104 + init_members();
  105 +
  106 + // Assign parameters
  107 + this->dir(dir);
  108 + this->semiminor(semiminor);
  109 + this->semimajor(semimajor);
  110 + this->posangle(posangle);
  111 + this->ridx(ridx);
  112 +
  113 + // Return
  114 + return;
  115 +}
  116 +
  117 +
  118 +/***********************************************************************//**
  119 + * @brief XML constructor
  120 + *
  121 + * @param[in] xml XML element.
  122 + *
  123 + * Constructs generalised elliptical Gaussian model by extracting information
  124 + * from an XML element. See the read() method for more information about the
  125 + * expected structure of the XML element.
  126 + ***************************************************************************/
  127 +GModelSpatialEllipticalGeneralGauss::GModelSpatialEllipticalGeneralGauss(const GXmlElement& xml) :
  128 + GModelSpatialElliptical()
  129 +{
  130 + // Initialise members
  131 + init_members();
  132 +
  133 + // Read information from XML element
  134 + read(xml);
  135 +
  136 + // Return
  137 + return;
  138 +}
  139 +
  140 +
  141 +/***********************************************************************//**
  142 + * @brief Copy constructor
  143 + *
  144 + * @param[in] model Generalised elliptical Gaussian model.
  145 + ***************************************************************************/
  146 +GModelSpatialEllipticalGeneralGauss::GModelSpatialEllipticalGeneralGauss(const GModelSpatialEllipticalGeneralGauss& model) :
  147 + GModelSpatialElliptical(model)
  148 +{
  149 + // Initialise members
  150 + init_members();
  151 +
  152 + // Copy members
  153 + copy_members(model);
  154 +
  155 + // Return
  156 + return;
  157 +}
  158 +
  159 +
  160 +/***********************************************************************//**
  161 + * @brief Destructor
  162 + ***************************************************************************/
  163 +GModelSpatialEllipticalGeneralGauss::~GModelSpatialEllipticalGeneralGauss(void)
  164 +{
  165 + // Free members
  166 + free_members();
  167 +
  168 + // Return
  169 + return;
  170 +}
  171 +
  172 +
  173 +/*==========================================================================
  174 + = =
  175 + = Operators =
  176 + = =
  177 + ==========================================================================*/
  178 +
  179 +/***********************************************************************//**
  180 + * @brief Assignment operator
  181 + *
  182 + * @param[in] model Generalised elliptical Gaussian model.
  183 + * @return Generalised elliptical Gaussian model.
  184 + ***************************************************************************/
  185 +GModelSpatialEllipticalGeneralGauss& GModelSpatialEllipticalGeneralGauss::operator=(const GModelSpatialEllipticalGeneralGauss& model)
  186 +{
  187 + // Execute only if object is not identical
  188 + if (this != &model) {
  189 +
  190 + // Copy base class members
  191 + this->GModelSpatialElliptical::operator=(model);
  192 +
  193 + // Free members
  194 + free_members();
  195 +
  196 + // Initialise members
  197 + init_members();
  198 +
  199 + // Copy members
  200 + copy_members(model);
  201 +
  202 + } // endif: object was not identical
  203 +
  204 + // Return
  205 + return *this;
  206 +}
  207 +
  208 +
  209 +/*==========================================================================
  210 + = =
  211 + = Public methods =
  212 + = =
  213 + ==========================================================================*/
  214 +
  215 +/***********************************************************************//**
  216 + * @brief Clear generalised elliptical Gaussian model
  217 + ***************************************************************************/
  218 +void GModelSpatialEllipticalGeneralGauss::clear(void)
  219 +{
  220 + // Free class members (base and derived classes, derived class first)
  221 + free_members();
  222 + this->GModelSpatialElliptical::free_members();
  223 + this->GModelSpatial::free_members();
  224 +
  225 + // Initialise members
  226 + this->GModelSpatial::init_members();
  227 + this->GModelSpatialElliptical::init_members();
  228 + init_members();
  229 +
  230 + // Return
  231 + return;
  232 +}
  233 +
  234 +
  235 +/***********************************************************************//**
  236 + * @brief Clone generalised elliptical Gaussian model
  237 + *
  238 + * @return Pointer to deep copy of generalised elliptical Gaussian model.
  239 + ***************************************************************************/
  240 +GModelSpatialEllipticalGeneralGauss* GModelSpatialEllipticalGeneralGauss::clone(void) const
  241 +{
  242 + // Clone generalised elliptical Gaussian model
  243 + return new GModelSpatialEllipticalGeneralGauss(*this);
  244 +}
  245 +
  246 +
  247 +/***********************************************************************//**
  248 + * @brief Evaluate function (in units of sr^-1)
  249 + *
  250 + * @param[in] theta Angular distance from disk centre (radians).
  251 + * @param[in] posangle Position angle (counterclockwise from North) (radians).
  252 + * @param[in] energy Photon energy.
  253 + * @param[in] time Photon arrival time.
  254 + * @param[in] gradients Compute gradients?
  255 + * @return Model value.
  256 + *
  257 + * Evaluates the spatial component for a generalised elliptical Gaussian
  258 + * source model.
  259 + *
  260 + * The generalised elliptical Gaussian function is defined by
  261 + * \f$S_{\rm p}(\theta, \phi | E, t)\f$, where
  262 + * \f$\theta\f$ is the angular separation between elliptical Gaussian centre
  263 + * and the actual location and \f$\phi\f$ the position angle with respect to
  264 + * the model centre, counted counterclockwise from North.
  265 + *
  266 + * The function \f$S_{\rm p}(\theta, \phi | E, t)\f$ is given by
  267 + *
  268 + * \f[
  269 + * S_{\rm p}(\theta, \phi | E, t) = {\tt m\_norm} \times
  270 + * \exp \left( -\left(\frac{\theta}{r_{\rm eff}}\right)^1/\eta \right)
  271 + * \f]
  272 + *
  273 + * where the effective ellipse radius \f$r_{\rm eff}\f$ towards a given
  274 + * position angle is given by
  275 + *
  276 + * \f[
  277 + * r_{\rm eff} = \frac{ab}
  278 + * {\sqrt{\left( a \sin (\phi - \phi_0) \right)^2} +
  279 + * \sqrt{\left( b \cos (\phi - \phi_0) \right)^2}}
  280 + * \f]
  281 + * and
  282 + * \f$a\f$ is the semi-major axis of the ellipse,
  283 + * \f$b\f$ is the semi-minor axis,
  284 + * \f$\phi_0\f$ is the position angle of the ellipse, counted
  285 + * counterclockwise from North, and
  286 + * \f$\et\f$ is the reciprocal of the radial profile index
  287 + * \f${\tt m\_norm}\f$ is a normalization constant given by
  288 + *
  289 + * \f[
  290 + * {\tt m\_norm} = \frac{1}{2 \pi a b \eta \Gamma(\eta)}
  291 + * \f]
  292 + *
  293 + * The method will not compute analytical parameter gradients, even if the
  294 + * @p gradients argument is set to true. Parameter gradients
  295 + * need to be computed numerically.
  296 + *
  297 + * @warning
  298 + * The above formula for an elliptical generalised Gaussian are accurate for small
  299 + * angles, with semimajor and semiminor axes below a few degrees. This
  300 + * should be fine for most use cases, but you should be aware that the
  301 + * ellipse will get distorted for larger angles.
  302 + * Note that the model normalisation is only correct in the small angle
  303 + * approximation and for \f$\eta\f$ order unity or smaller.
  304 + ***************************************************************************/
  305 +double GModelSpatialEllipticalGeneralGauss::eval(const double& theta,
  306 + const double& posangle,
  307 + const GEnergy& energy,
  308 + const GTime& time,
  309 + const bool& gradients) const
  310 +{
  311 + // Initialise value
  312 + double value = 0.0;
  313 +
  314 + // Continue only if we're inside circle enclosing the ellipse
  315 + if (theta <= theta_max()) {
  316 +
  317 + // Update precomputation cache
  318 + update();
  319 +
  320 + // Perform computations to compute exponent
  321 + #if defined(G_SMALL_ANGLE_APPROXIMATION)
  322 + double rel_posangle = posangle - this->posangle() * gammalib::deg2rad;
  323 + double cosinus = std::cos(rel_posangle);
  324 + double sinus = std::sin(rel_posangle);
  325 + double arg1 = m_minor_rad * cosinus;
  326 + double arg2 = m_major_rad * sinus;
  327 + double r_ellipse = m_minor_rad * m_major_rad /
  328 + std::sqrt(arg1*arg1 + arg2*arg2); //!< small angle
  329 + double r_relative = theta/r_ellipse;
  330 + double exponent = std::pow(r_relative,m_inv_ridx);
  331 + #else
  332 + // Perform computations
  333 + double sinphi = std::sin(posangle);
  334 + double cosphi = std::cos(posangle);
  335 +
  336 + // Compute help term
  337 + double term1 = m_term1 * cosphi * cosphi;
  338 + double term2 = m_term2 * sinphi * sinphi;
  339 + double term3 = m_term3 * sinphi * cosphi;
  340 +
  341 + // Compute exponent
  342 + double exponent = std::pow(theta * std::sqrt(term1 + term2 + term3), m_inv_ridx) ;
  343 + #endif
  344 +
  345 + // Set value
  346 + value = m_norm * std::exp(-exponent);
  347 +
  348 + // Compile option: Check for NaN/Inf
  349 + #if defined(G_NAN_CHECK)
  350 + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
  351 + std::cout << "*** ERROR: GModelSpatialEllipticalGeneralGauss::eval";
  352 + std::cout << "(theta=" << theta << "): NaN/Inf encountered";
  353 + std::cout << "(posangle=" << posangle << "): NaN/Inf encountered";
  354 + std::cout << " (value=" << value;
  355 + std::cout << ", MinorAxis^2=" << m_minor2;
  356 + std::cout << ", MaxjorAxis^2=" << m_major2;
  357 + std::cout << ", m_norm=" << m_norm;
  358 + std::cout << ")" << std::endl;
  359 + }
  360 + #endif
  361 +
  362 + } // endif: position was inside enclosing circle
  363 +
  364 + // Return value
  365 + return value;
  366 +}
  367 +
  368 +
  369 +/***********************************************************************//**
  370 + * @brief Returns MC sky direction
  371 + *
  372 + * @param[in] energy Photon energy.
  373 + * @param[in] time Photon arrival time.
  374 + * @param[in,out] ran Random number generator.
  375 + * @return Sky direction.
  376 + *
  377 + * Draws an arbitrary sky direction from the elliptical generalised
  378 + * Gaussian model.
  379 + ***************************************************************************/
  380 +GSkyDir GModelSpatialEllipticalGeneralGauss::mc(const GEnergy& energy,
  381 + const GTime& time,
  382 + GRan& ran) const
  383 +{
  384 + // Update precomputation cache
  385 + update();
  386 +
  387 + // Initialise simulation
  388 + #if defined(G_DEBUG_MC)
  389 + int n_samples = 0;
  390 + #endif
  391 + double u_max = m_norm * std::sin(theta_max());
  392 + double value = 0.0;
  393 + double u = 1.0;
  394 + double theta = 0.0;
  395 + double phi = 0.0;
  396 +
  397 + // Throw an exception if the maximum value is zero
  398 + if (u_max <= 0.0) {
  399 + std::string msg = "Non positive maximum ellpitical general Gauss model value.";
  400 + throw GException::invalid_return_value(G_MC, msg);
  401 + }
  402 +
  403 + // Simulate offset from photon arrival direction using a rejection method
  404 + do {
  405 + theta = ran.uniform() * theta_max();
  406 + phi = ran.uniform() * gammalib::twopi;
  407 + value = eval(theta, phi, energy, time) * std::sin(theta);
  408 + u = ran.uniform() * u_max;
  409 + #if defined(G_DEBUG_MC)
  410 + n_samples++;
  411 + #endif
  412 + } while (u > value);
  413 + #if defined(G_DEBUG_MC)
  414 + std::cout << "#=" << n_samples << " ";
  415 + #endif
  416 +
  417 + // Rotate pointing direction by offset and azimuth angle
  418 + GSkyDir sky_dir = dir();
  419 + sky_dir.rotate_deg(phi * gammalib::rad2deg, theta * gammalib::rad2deg);
  420 +
  421 + // Return sky direction
  422 + return sky_dir;
  423 +}
  424 +
  425 +
  426 +/***********************************************************************//**
  427 + * @brief Checks whether model contains specified sky direction
  428 + *
  429 + * @param[in] dir Sky direction.
  430 + * @param[in] margin Margin to be added to sky direction (deg).
  431 + * @return True if the model contains the sky direction.
  432 + *
  433 + * Signals whether a sky direction is contained in the elliptical gauss
  434 + * model.
  435 + *
  436 + * @todo Implement correct evaluation of effective ellipse radius.
  437 + ***************************************************************************/
  438 +bool GModelSpatialEllipticalGeneralGauss::contains(const GSkyDir& dir,
  439 + const double& margin) const
  440 +{
  441 + // Compute distance to centre (radian)
  442 + double distance = dir.dist(this->dir());
  443 +
  444 + // Return flag
  445 + return (distance <= theta_max() + margin*gammalib::deg2rad);
  446 +}
  447 +
  448 +
  449 +/***********************************************************************//**
  450 + * @brief Return maximum model radius (in radians)
  451 + *
  452 + * @return Returns maximum model radius (radians).
  453 + *
  454 + * Returns the maximum of 3.0 semimajor() and 3.0 semiminor() as
  455 + * approximate edge of the Gaussian. This limit is of course arbitrary, but
  456 + * allows to limit the integration region for response computation. The value
  457 + * of 3.0 has been determined by experiment (#1561).
  458 + ***************************************************************************/
  459 +double GModelSpatialEllipticalGeneralGauss::theta_max(void) const
  460 +{
  461 + // Set maximum model radius
  462 + double theta_max = (semimajor() > semiminor())
  463 + ? semimajor() * gammalib::deg2rad * c_theta_max
  464 + : semiminor() * gammalib::deg2rad * c_theta_max;
  465 +
  466 + // Return value
  467 + return theta_max;
  468 +}
  469 +
  470 +
  471 +/***********************************************************************//**
  472 + * @brief Read model from XML element
  473 + *
  474 + * @param[in] xml XML element.
  475 + *
  476 + * Reads the elliptical gauss model information from an XML element. The XML
  477 + * element shall have either the format
  478 + *
  479 + * <spatialModel type="EllipticalGeneralGaussian">
  480 + * <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
  481 + * <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
  482 + * <parameter name="PA" scale="1.0" value="45.0" min="-360" max="360" free="1"/>
  483 + * <parameter name="MinorRadius" scale="1.0" value="0.5" min="0.001" max="10" free="1"/>
  484 + * <parameter name="MajorRadius" scale="1.0" value="2.0" min="0.001" max="10" free="1"/>
  485 + * <parameter name="R_Index" scale="1.0" value="0.5" min="0.01" max="10" free="1"/>
  486 + * </spatialModel>
  487 + *
  488 + * or
  489 + *
  490 + * <spatialModel type="EllipticalGeneralGaussian">
  491 + * <parameter name="GLON" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
  492 + * <parameter name="GLAT" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
  493 + * <parameter name="PA" scale="1.0" value="45.0" min="-360" max="360" free="1"/>
  494 + * <parameter name="MinorRadius" scale="1.0" value="0.5" min="0.001" max="10" free="1"/>
  495 + * <parameter name="MajorRadius" scale="1.0" value="2.0" min="0.001" max="10" free="1"/>
  496 + * <parameter name="R_Index" scale="1.0" value="0.5" min="0.01" max="10" free="1"/>
  497 + * </spatialModel>
  498 + *
  499 + * @todo Implement a test of the ellipse boundary. The axes
  500 + * and axes minimum should be >0.
  501 + ***************************************************************************/
  502 +void GModelSpatialEllipticalGeneralGauss::read(const GXmlElement& xml)
  503 +{
  504 + // Verify number of model parameters
  505 + gammalib::xml_check_parnum(G_READ, xml, 6);
  506 +
  507 + // Read gauss location
  508 + GModelSpatialElliptical::read(xml);
  509 +
  510 + // Get parameters
  511 + const GXmlElement* minor = gammalib::xml_get_par(G_READ, xml, m_semiminor.name());
  512 + const GXmlElement* major = gammalib::xml_get_par(G_READ, xml, m_semimajor.name());
  513 + const GXmlElement* ridx = gammalib::xml_get_par(G_READ, xml, m_ridx.name());
  514 +
  515 + // Read parameters
  516 + m_semiminor.read(*minor);
  517 + m_semimajor.read(*major);
  518 + m_ridx.read(*ridx);
  519 +
  520 + // Return
  521 + return;
  522 +}
  523 +
  524 +
  525 +/***********************************************************************//**
  526 + * @brief Write model into XML element
  527 + *
  528 + * @param[in] xml XML element into which model information is written.
  529 + *
  530 + * Write the elliptical gauss model information into an XML element. The XML
  531 + * element will have the format
  532 + *
  533 + * <spatialModel type="EllipticalGeneralGaussian">
  534 + * <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
  535 + * <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
  536 + * <parameter name="PA" scale="1.0" value="45.0" min="-360" max="360" free="1"/>
  537 + * <parameter name="MinorRadius" scale="1.0" value="0.5" min="0.001" max="10" free="1"/>
  538 + * <parameter name="MajorRadius" scale="1.0" value="2.0" min="0.001" max="10" free="1"/>
  539 + * <parameter name="R_Index" scale="1.0" value="0.5" min="0.01" max="10" free="1"/>
  540 + * </spatialModel>
  541 + *
  542 + ***************************************************************************/
  543 +void GModelSpatialEllipticalGeneralGauss::write(GXmlElement& xml) const
  544 +{
  545 + // Verify model type
  546 + gammalib::xml_check_type(G_WRITE, xml, type());
  547 +
  548 + // Write disk location
  549 + GModelSpatialElliptical::write(xml);
  550 +
  551 + // Get or create parameters
  552 + GXmlElement* minor = gammalib::xml_need_par(G_WRITE, xml, m_semiminor.name());
  553 + GXmlElement* major = gammalib::xml_need_par(G_WRITE, xml, m_semimajor.name());
  554 + GXmlElement* ridx = gammalib::xml_need_par(G_WRITE, xml, m_ridx.name());
  555 +
  556 + // Write parameters
  557 + m_semiminor.write(*minor);
  558 + m_semimajor.write(*major);
  559 + m_ridx.write(*ridx);
  560 +
  561 + // Return
  562 + return;
  563 +}
  564 +
  565 +
  566 +/***********************************************************************//**
  567 + * @brief Print information
  568 + *
  569 + * @param[in] chatter Chattiness.
  570 + * @return String containing model information.
  571 + ***************************************************************************/
  572 +std::string GModelSpatialEllipticalGeneralGauss::print(const GChatter& chatter) const
  573 +{
  574 + // Initialise result string
  575 + std::string result;
  576 +
  577 + // Continue only if chatter is not silent
  578 + if (chatter != SILENT) {
  579 +
  580 + // Append header
  581 + result.append("=== GModelSpatialEllipticalGeneralGauss ===");
  582 +
  583 + // Append parameters
  584 + result.append("\n"+gammalib::parformat("Number of parameters"));
  585 + result.append(gammalib::str(size()));
  586 + for (int i = 0; i < size(); ++i) {
  587 + result.append("\n"+m_pars[i]->print(chatter));
  588 + }
  589 +
  590 + } // endif: chatter was not silent
  591 +
  592 + // Return result
  593 + return result;
  594 +}
  595 +
  596 +
  597 +/*==========================================================================
  598 + = =
  599 + = Private methods =
  600 + = =
  601 + ==========================================================================*/
  602 +
  603 +/***********************************************************************//**
  604 + * @brief Initialise class members
  605 + ***************************************************************************/
  606 +void GModelSpatialEllipticalGeneralGauss::init_members(void)
  607 +{
  608 + // Initialise model type
  609 + m_type = "EllipticalGeneralGaussian";
  610 +
  611 + // Initialise Reciprocal index
  612 + m_ridx.clear();
  613 + m_ridx.name("R_Index");
  614 + m_ridx.value(0.5); // Gaussian
  615 + m_ridx.min(1.e-15); // need to be > 0
  616 + m_ridx.free();
  617 + m_ridx.scale(1.0);
  618 + m_ridx.gradient(0.0);
  619 + m_ridx.has_grad(false);
  620 +
  621 + // Set parameter pointer(s)
  622 + m_pars.push_back(&m_ridx);
  623 +
  624 + // Initialise precomputation cache. Note that zero values flag
  625 + // uninitialised as a zero radius is not meaningful
  626 + m_last_minor = 0.0;
  627 + m_last_major = 0.0;
  628 + m_minor_rad = 0.0;
  629 + m_major_rad = 0.0;
  630 + m_norm = 0.0;
  631 + m_last_posangle = 9999.0; // Signals that has not been initialised
  632 + m_sin2pos = 0.0;
  633 + m_cospos2 = 0.0;
  634 + m_sinpos2 = 0.0;
  635 + m_minor2 = 0.0;
  636 + m_major2 = 0.0;
  637 + m_term1 = 0.0;
  638 + m_term2 = 0.0;
  639 + m_term3 = 0.0;
  640 + m_last_ridx = 0.0;
  641 + m_inv_ridx = 0.0;
  642 +
  643 + // Return
  644 + return;
  645 +}
  646 +
  647 +
  648 +/***********************************************************************//**
  649 + * @brief Copy class members
  650 + *
  651 + * @param[in] model Elliptical Gaussian model.
  652 + ***************************************************************************/
  653 +void GModelSpatialEllipticalGeneralGauss::copy_members(const GModelSpatialEllipticalGeneralGauss& model)
  654 +{
  655 + // Copy members
  656 + m_type = model.m_type; // Needed to conserve model type
  657 + m_ridx = model.m_ridx;
  658 +
  659 + // Copy precomputation cache
  660 + m_last_minor = model.m_last_minor;
  661 + m_last_major = model.m_last_major;
  662 + m_minor_rad = model.m_minor_rad;
  663 + m_major_rad = model.m_major_rad;
  664 + m_norm = model.m_norm;
  665 + m_last_posangle = model.m_last_posangle;
  666 + m_sin2pos = model.m_sin2pos;
  667 + m_cospos2 = model.m_cospos2;
  668 + m_sinpos2 = model.m_sinpos2;
  669 + m_minor2 = model.m_minor2;
  670 + m_major2 = model.m_major2;
  671 + m_term1 = model.m_term1;
  672 + m_term2 = model.m_term2;
  673 + m_term3 = model.m_term3;
  674 + m_last_ridx = model.m_last_ridx;
  675 + m_inv_ridx = model.m_inv_ridx;
  676 +
  677 + // Return
  678 + return;
  679 +}
  680 +
  681 +
  682 +/***********************************************************************//**
  683 + * @brief Delete class members
  684 + ***************************************************************************/
  685 +void GModelSpatialEllipticalGeneralGauss::free_members(void)
  686 +{
  687 + // Return
  688 + return;
  689 +}
  690 +
  691 +
  692 +/***********************************************************************//**
  693 + * @brief Update precomputation cache
  694 + *
  695 + * Computes the normalization
  696 + * \f[
  697 + * {\tt m\_norm} = \frac{1}{2 \pi \times a \times b}
  698 + * \f]
  699 + * where
  700 + * \f$a\f$ is the semi-major axis and
  701 + * \f$b\f$ is the semi-minor axis of the ellipse.
  702 + *
  703 + * @warning
  704 + * The normalization of the elliptical Gaussian is only valid in the small
  705 + * angle approximation.
  706 + ***************************************************************************/
  707 +void GModelSpatialEllipticalGeneralGauss::update() const
  708 +{
  709 + // Initialise flag if something has changed
  710 + bool changed = false;
  711 +
  712 + // Update if one axis has changed or radial profile index has changed
  713 + if (m_last_minor != semiminor() || m_last_major != semimajor() || m_last_ridx != ridx()) {
  714 +
  715 + // Signal parameter changes
  716 + changed = true;
  717 +
  718 + // Store last values
  719 + m_last_minor = semiminor();
  720 + m_last_major = semimajor();
  721 + m_last_ridx = ridx();
  722 +
  723 + // Compute axes in radians
  724 + m_minor_rad = m_last_minor * gammalib::deg2rad;
  725 + m_major_rad = m_last_major * gammalib::deg2rad;
  726 +
  727 + // Perform precomputations. Note also
  728 + // that the model normalisation is only correct in the small angle
  729 + // approximation and for ridx order unity or smaller.
  730 + if (m_minor_rad > 0.0 and m_major_rad > 0.0 and ridx() > 0.0) {
  731 + m_inv_ridx = 1.0 / ridx() ;
  732 + m_norm = 1.0 / (gammalib::twopi * m_minor_rad * m_major_rad * ridx() *
  733 + std::exp(gammalib::gammln(2.0 * ridx())) );
  734 + }
  735 + else {
  736 + m_inv_ridx = 1.0;
  737 + m_norm = 0.0;
  738 + }
  739 +
  740 + } // endif: update required
  741 +
  742 + // Update chache if position angle changed
  743 + #if !defined(G_SMALL_ANGLE_APPROXIMATION)
  744 + if (m_last_posangle != posangle()) {
  745 +
  746 + // Signal parameter changes
  747 + changed = true;
  748 +
  749 + // Store last value
  750 + m_last_posangle = posangle();
  751 +
  752 + // Compute angle in radians
  753 + double posangle_rad = m_last_posangle * gammalib::deg2rad;
  754 +
  755 + // Compute sine and cosine
  756 + double cospos = std::cos(posangle_rad);
  757 + double sinpos = std::sin(posangle_rad);
  758 +
  759 + // Cache important values for further computations
  760 + m_cospos2 = cospos * cospos;
  761 + m_sinpos2 = sinpos * sinpos;
  762 + m_sin2pos = std::sin(2.0 * posangle_rad);
  763 +
  764 + } // endif: position angle update required
  765 +
  766 + // Perform precomputations in case anything has changed
  767 + if (changed) {
  768 +
  769 + // Compute help terms
  770 + m_term1 = 0.5 * (m_cospos2 / m_minor2 + m_sinpos2 / m_major2);
  771 + m_term2 = 0.5 * (m_cospos2 / m_major2 + m_sinpos2 / m_minor2);
  772 + m_term3 = 0.5 * (m_sin2pos / m_major2 - m_sin2pos / m_minor2);
  773 +
  774 + } // endif: something has changed
  775 + #endif
  776 +
  777 + // Return
  778 + return;
  779 +}
  780 +
  781 +
  782 +/***********************************************************************//**
  783 + * @brief Set boundary sky region
  784 + ***************************************************************************/
  785 +void GModelSpatialEllipticalGeneralGauss::set_region(void) const
  786 +{
  787 + // Set maximum model radius
  788 + double max_radius = (semimajor() > semiminor()) ? semimajor() : semiminor();
  789 +
  790 + // Set sky region circle (maximum Gaussian sigma times a scaling
  791 + // factor (actually 5))
  792 + GSkyRegionCircle region(m_ra.value(), m_dec.value(), max_radius * c_theta_max);
  793 +
  794 + // Set region (circumvent const correctness)
  795 + const_cast<GModelSpatialEllipticalGeneralGauss*>(this)->m_region = region;
  796 +
  797 + // Return
  798 + return;
  799 +}
... ...
src/model/GModelSpatialRadialGeneralGauss.cpp 0 → 100644
  View file @e9681ff
  1 +/***************************************************************************
  2 + * GModelSpatialRadialGeneralGauss.cpp - Generalised radial Gaussian *
  3 + * source model class *
  4 + * ----------------------------------------------------------------------- *
  5 + * copyright (C) 2021-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 GModelSpatialRadialGeneralGauss.cpp
  24 + * @brief Radial Generallized 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"
  35 +#include "GModelSpatialRadialGeneralGauss.hpp"
  36 +#include "GModelSpatialRegistry.hpp"
  37 +
  38 +/* __ Constants __________________________________________________________ */
  39 +
  40 +/* __ Globals ____________________________________________________________ */
  41 +const GModelSpatialRadialGeneralGauss g_radial_general_gauss_seed;
  42 +const GModelSpatialRegistry g_radial_general_gauss_registry(&g_radial_general_gauss_seed);
  43 +
  44 +/* __ Method name definitions ____________________________________________ */
  45 +#define G_EVAL "GModelSpatialRadialGeneralGauss::eval(double&, GEnergy&, "\
  46 + "GTime&, bool&)"
  47 +#define G_READ "GModelSpatialRadialGeneralGauss::read(GXmlElement&)"
  48 +#define G_WRITE "GModelSpatialRadialGeneralGauss::write(GXmlElement&)"
  49 +#define G_MC "GModelSpatialRadialGeneralGauss::mc(GEnergy&, GTime&, GRan&)"
  50 +
  51 +/* __ Macros _____________________________________________________________ */
  52 +
  53 +/* __ Coding definitions _________________________________________________ */
  54 +
  55 +/* __ Debug definitions __________________________________________________ */
  56 +
  57 +
  58 +/*==========================================================================
  59 + = =
  60 + = Constructors/destructors =
  61 + = =
  62 + ==========================================================================*/
  63 +
  64 +/***********************************************************************//**
  65 + * @brief Void constructor
  66 + *
  67 + * Constructs empty radial Gaussian model.
  68 + ***************************************************************************/
  69 +GModelSpatialRadialGeneralGauss::GModelSpatialRadialGeneralGauss(void) :
  70 + GModelSpatialRadial()
  71 +{
  72 + // Initialise members
  73 + init_members();
  74 +
  75 + // Return
  76 + return;
  77 +}
  78 +
  79 +
  80 +/***********************************************************************//**
  81 + * @brief Constructor
  82 + *
  83 + * @param[in] dir Sky position of Gaussian.
  84 + * @param[in] radius Width of generalised Gaussian (deg).
  85 + * @param[in] ridx Reciprocal of exponential index of radial profile.
  86 + *
  87 + * Constructs a generalised Gaussian spatial model using a sky direction
  88 + * (@p dir), a Gaussian width parameter @p radius in degrees and a reciprocal
  89 + * of the exponential index @p ridx.
  90 + ***************************************************************************/
  91 +GModelSpatialRadialGeneralGauss::GModelSpatialRadialGeneralGauss(const GSkyDir& dir,
  92 + const double& radius,
  93 + const double& ridx) :
  94 + GModelSpatialRadial()
  95 +{
  96 + // Initialise members
  97 + init_members();
  98 +
  99 + // Assign direction and radius
  100 + this->dir(dir);
  101 + this->radius(radius);
  102 + this->ridx(ridx);
  103 +
  104 + // Return
  105 + return;
  106 +}
  107 +
  108 +
  109 +/***********************************************************************//**
  110 + * @brief XML constructor
  111 + *
  112 + * @param[in] xml XML element.
  113 + *
  114 + * Constructs a generalised Gaussian spatial model by extracting information
  115 + * from an XML element. See the method read() for more information about the
  116 + * expected structure of the XML element.
  117 + ***************************************************************************/
  118 +GModelSpatialRadialGeneralGauss::GModelSpatialRadialGeneralGauss(const GXmlElement& xml) :
  119 + GModelSpatialRadial()
  120 +{
  121 + // Initialise members
  122 + init_members();
  123 +
  124 + // Read information from XML element
  125 + read(xml);
  126 +
  127 + // Return
  128 + return;
  129 +}
  130 +
  131 +
  132 +/***********************************************************************//**
  133 + * @brief Copy constructor
  134 + *
  135 + * @param[in] model Generalised radial Gaussian model.
  136 + ***************************************************************************/
  137 +GModelSpatialRadialGeneralGauss::GModelSpatialRadialGeneralGauss(const GModelSpatialRadialGeneralGauss& model) :
  138 + GModelSpatialRadial(model)
  139 +{
  140 + // Initialise members
  141 + init_members();
  142 +
  143 + // Copy members
  144 + copy_members(model);
  145 +
  146 + // Return
  147 + return;
  148 +}
  149 +
  150 +
  151 +/***********************************************************************//**
  152 + * @brief Destructor
  153 + ***************************************************************************/
  154 +GModelSpatialRadialGeneralGauss::~GModelSpatialRadialGeneralGauss(void)
  155 +{
  156 + // Free members
  157 + free_members();
  158 +
  159 + // Return
  160 + return;
  161 +}
  162 +
  163 +
  164 +/*==========================================================================
  165 + = =
  166 + = Operators =
  167 + = =
  168 + ==========================================================================*/
  169 +
  170 +/***********************************************************************//**
  171 + * @brief Assignment operator
  172 + *
  173 + * @param[in] model Generalised radial Gaussian model.
  174 + * @return Radial Generalised radial Gaussian model.
  175 + ***************************************************************************/
  176 +GModelSpatialRadialGeneralGauss& GModelSpatialRadialGeneralGauss::operator=(const GModelSpatialRadialGeneralGauss& model)
  177 +{
  178 + // Execute only if object is not identical
  179 + if (this != &model) {
  180 +
  181 + // Copy base class members
  182 + this->GModelSpatialRadial::operator=(model);
  183 +
  184 + // Free members
  185 + free_members();
  186 +
  187 + // Initialise members
  188 + init_members();
  189 +
  190 + // Copy members
  191 + copy_members(model);
  192 +
  193 + } // endif: object was not identical
  194 +
  195 + // Return
  196 + return *this;
  197 +}
  198 +
  199 +
  200 +/*==========================================================================
  201 + = =
  202 + = Public methods =
  203 + = =
  204 + ==========================================================================*/
  205 +
  206 +/***********************************************************************//**
  207 + * @brief Clear radial Gauss model
  208 + ***************************************************************************/
  209 +void GModelSpatialRadialGeneralGauss::clear(void)
  210 +{
  211 + // Free class members (base and derived classes, derived class first)
  212 + free_members();
  213 + this->GModelSpatialRadial::free_members();
  214 + this->GModelSpatial::free_members();
  215 +
  216 + // Initialise members
  217 + this->GModelSpatial::init_members();
  218 + this->GModelSpatialRadial::init_members();
  219 + init_members();
  220 +
  221 + // Return
  222 + return;
  223 +}
  224 +
  225 +
  226 +/***********************************************************************//**
  227 + * @brief Clone radial Gauss model
  228 + ***************************************************************************/
  229 +GModelSpatialRadialGeneralGauss* GModelSpatialRadialGeneralGauss::clone(void) const
  230 +{
  231 + // Clone radial Gauss model
  232 + return new GModelSpatialRadialGeneralGauss(*this);
  233 +}
  234 +
  235 +
  236 +/***********************************************************************//**
  237 + * @brief Evaluate Generalised Gaussian source model
  238 + *
  239 + * @param[in] theta Angular distance from source centre (radians).
  240 + * @param[in] energy Photon energy (not used).
  241 + * @param[in] time Photon arrival time (not used).
  242 + * @param[in] gradients Compute gradients?
  243 + * @return Model value (\f${\rm sr}^{-1}\f$).
  244 + *
  245 + * Evaluates the spatial component for a Generalised Gaussian source model using
  246 + *
  247 + * \f[
  248 + * M_{\rm S}({\bf p} | E, t) = M_{\rm S}(\theta) =
  249 + * \frac{1}{2 \pi r^2 \eta \Gamma(\eta)} \exp
  250 + * \left(-(\frac{\theta}{r})^1/\eta \right)
  251 + * \f]
  252 + *
  253 + * where
  254 + * - \f$\theta\f$ is the angular separation from the centre of the model, and
  255 + * - \f$r\f$ is the Generalised Gaussian radius.
  256 + * - \f$\eta\f$ is the reciprocal of the radial profile exponent.
  257 + *
  258 + * If @p gradients is `true` the method will also compute parameter
  259 + * gradients
  260 + *
  261 + *
  262 + * Note that the model normalisation is only correct in the small angle
  263 + * approximation and for \f$\eta\f$ order unity or smaller.
  264 + ***************************************************************************/
  265 +double GModelSpatialRadialGeneralGauss::eval(const double& theta,
  266 + const GEnergy& energy,
  267 + const GTime& time,
  268 + const bool& gradients) const
  269 +{
  270 + // Update
  271 + update();
  272 +
  273 + // Compute value
  274 + double arg = -1.0 * std::pow(m_inv_radius_rad * theta,m_inv_ridx);
  275 + double value = m_value_norm * std::exp(arg);
  276 +
  277 + // Compile option: Check for NaN/Inf
  278 + #if defined(G_NAN_CHECK)
  279 + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
  280 + std::string msg = "Model value not a number:";
  281 + for (int i = 0; i < m_pars.size(); ++i) {
  282 + msg += " " + m_pars[i]->name() + "=";
  283 + msg += gammalib::str(m_pars[i]->value());
  284 + }
  285 + msg += " energy=" + energy.print();
  286 + msg += " time=" + time.print();
  287 + gammalib::warning(G_EVAL, msg);
  288 + }
  289 + #endif
  290 +
  291 + // Return value
  292 + return value;
  293 +}
  294 +
  295 +
  296 +/***********************************************************************//**
  297 + * @brief Returns MC sky direction
  298 + *
  299 + * @param[in] energy Photon energy.
  300 + * @param[in] time Photon arrival time.
  301 + * @param[in,out] ran Random number generator.
  302 + * @return Sky direction.
  303 + *
  304 + * Draws an arbitrary sky direction from the 2D Generalised Gaussian
  305 + * distribution as function of the photon @p energy and arrival @p time.
  306 + *
  307 + * @todo This method is only valid in the small angle approximation.
  308 + ***************************************************************************/
  309 +GSkyDir GModelSpatialRadialGeneralGauss::mc(const GEnergy& energy,
  310 + const GTime& time,
  311 + GRan& ran) const
  312 +{
  313 +
  314 + // Update precomputation cache
  315 + update();
  316 +
  317 + // Initialise simulation
  318 + #if defined(G_DEBUG_MC)
  319 + int n_samples = 0;
  320 + #endif
  321 + double u_max = m_value_norm * std::sin(theta_max());
  322 + double value = 0.0;
  323 + double u = 1.0;
  324 + double theta = 0.0;
  325 +
  326 + // Throw an exception if the maximum value is zero
  327 + if (u_max <= 0.0) {
  328 + std::string msg = "Non positive maximum radial genral Gauss model value.";
  329 + throw GException::invalid_return_value(G_MC, msg);
  330 + }
  331 +
  332 + // Simulate offset from photon arrival direction using a rejection method
  333 + do {
  334 + theta = ran.uniform() * theta_max();
  335 + value = eval(theta, energy, time) * std::sin(theta);
  336 + u = ran.uniform() * u_max;
  337 + #if defined(G_DEBUG_MC)
  338 + n_samples++;
  339 + #endif
  340 + } while (u > value);
  341 + #if defined(G_DEBUG_MC)
  342 + std::cout << "#=" << n_samples << " ";
  343 + #endif
  344 +
  345 + // Simulate azimuth angle
  346 + double phi = 360.0 * ran.uniform();
  347 +
  348 + // Rotate pointing direction by offset and azimuth angle
  349 + GSkyDir sky_dir = dir();
  350 + sky_dir.rotate_deg(phi, theta * gammalib::rad2deg);
  351 +
  352 + // Return sky direction
  353 + return sky_dir;
  354 +
  355 +}
  356 +
  357 +
  358 +/***********************************************************************//**
  359 + * @brief Checks where model contains specified sky direction
  360 + *
  361 + * @param[in] dir Sky direction.
  362 + * @param[in] margin Margin to be added to sky direction (deg)
  363 + * @return True if the model contains the sky direction.
  364 + *
  365 + * Signals whether a sky direction is contained in the radial gauss model.
  366 + ***************************************************************************/
  367 +bool GModelSpatialRadialGeneralGauss::contains(const GSkyDir& dir,
  368 + const double& margin) const
  369 +{
  370 + // Compute distance to centre (radians)
  371 + double distance = dir.dist(this->dir());
  372 +
  373 + // Return flag
  374 + return (distance <= theta_max() + margin*gammalib::deg2rad);
  375 +}
  376 +
  377 +
  378 +/***********************************************************************//**
  379 + * @brief Return maximum model radius (in radians)
  380 + *
  381 + * Returns \f$5 \sigma\f$ as approximate edge of the Gaussian. This limit
  382 + * is of course arbitrary, but allows to limit the integration region for
  383 + * response computation.
  384 + ***************************************************************************/
  385 +double GModelSpatialRadialGeneralGauss::theta_max(void) const
  386 +{
  387 + // Return value
  388 + return (radius() * gammalib::deg2rad * 5.0);
  389 +}
  390 +
  391 +
  392 +/***********************************************************************//**
  393 + * @brief Read model from XML element
  394 + *
  395 + * @param[in] xml XML element.
  396 + *
  397 + * Reads the generalised radial Gaussian model information from an XML
  398 + * element. The XML element shall have either the format
  399 + *
  400 + * <spatialModel type="RadialGeneralGaussian">
  401 + * <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
  402 + * <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
  403 + * <parameter name="Radius" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
  404 + * <parameter name="R_Index" scale="1.0" value="0.5" min="0.01" max="10" free="1"/>
  405 + * </spatialModel>
  406 + *
  407 + * or
  408 + *
  409 + * <spatialModel type="RadialGeneralGaussian">
  410 + * <parameter name="GLON" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
  411 + * <parameter name="GLAT" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
  412 + * <parameter name="Radius" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
  413 + * <parameter name="R_Index" scale="1.0" value="0.5" min="0.01" max="10" free="1"/>
  414 + * </spatialModel>
  415 + *
  416 + * @todo Implement a test of the radius and radius boundary. The sigma
  417 + * and sigma minimum should be >0.
  418 + ***************************************************************************/
  419 +void GModelSpatialRadialGeneralGauss::read(const GXmlElement& xml)
  420 +{
  421 + // Verify number of model parameters
  422 + gammalib::xml_check_parnum(G_READ, xml, 4);
  423 +
  424 + // Read location
  425 + GModelSpatialRadial::read(xml);
  426 +
  427 + // Get parameters
  428 + const GXmlElement* radius = gammalib::xml_get_par(G_READ, xml, m_radius.name());
  429 + const GXmlElement* ridx = gammalib::xml_get_par(G_READ, xml, m_ridx.name());
  430 +
  431 + // Read parameters
  432 + m_radius.read(*radius);
  433 + m_ridx.read(*ridx);
  434 +
  435 + // Return
  436 + return;
  437 +}
  438 +
  439 +
  440 +/***********************************************************************//**
  441 + * @brief Write model into XML element
  442 + *
  443 + * @param[in] xml XML element into which model information is written.
  444 + *
  445 + * Writes the radial disk model information into an XML element. The XML
  446 + * element will have the format
  447 + *
  448 + * <spatialModel type="RadialGeneralGaussian">
  449 + * <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
  450 + * <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
  451 + * <parameter name="Sigma" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
  452 + * <parameter name="R_Index" scale="1.0" value="0.5" min="0.01" max="10" free="1"/>
  453 + * </spatialModel>
  454 + *
  455 + ***************************************************************************/
  456 +void GModelSpatialRadialGeneralGauss::write(GXmlElement& xml) const
  457 +{
  458 + // Write Gaussian location
  459 + GModelSpatialRadial::write(xml);
  460 +
  461 + // Get or create parameters
  462 + GXmlElement* radius = gammalib::xml_need_par(G_WRITE, xml, m_radius.name());
  463 + GXmlElement* ridx = gammalib::xml_need_par(G_WRITE, xml, m_ridx.name());
  464 +
  465 + // Write parameters
  466 + m_radius.write(*radius);
  467 + m_ridx.write(*ridx);
  468 +
  469 + // Return
  470 + return;
  471 +}
  472 +
  473 +
  474 +/***********************************************************************//**
  475 + * @brief Print generalised radial Gaussian source information
  476 + *
  477 + * @param[in] chatter Chattiness.
  478 + * @return String containing model information.
  479 + ***************************************************************************/
  480 +std::string GModelSpatialRadialGeneralGauss::print(const GChatter& chatter) const
  481 +{
  482 + // Initialise result string
  483 + std::string result;
  484 +
  485 + // Continue only if chatter is not silent
  486 + if (chatter != SILENT) {
  487 +
  488 + // Append header
  489 + result.append("=== GModelSpatialRadialGeneralGauss ===");
  490 +
  491 + // Append parameters
  492 + result.append("\n"+gammalib::parformat("Number of parameters"));
  493 + result.append(gammalib::str(size()));
  494 + for (int i = 0; i < size(); ++i) {
  495 + result.append("\n"+m_pars[i]->print(chatter));
  496 + }
  497 +
  498 + } // endif: chatter was not silent
  499 +
  500 + // Return result
  501 + return result;
  502 +}
  503 +
  504 +
  505 +/*==========================================================================
  506 + = =
  507 + = Private methods =
  508 + = =
  509 + ==========================================================================*/
  510 +
  511 +/***********************************************************************//**
  512 + * @brief Initialise class members
  513 + *
  514 + * Note that the minimum Gaussian width is set to 1 arcsec
  515 + * and w set the minimum reciprocal index to 1.e-15
  516 + ***************************************************************************/
  517 +void GModelSpatialRadialGeneralGauss::init_members(void)
  518 +{
  519 + // Initialise model type
  520 + m_type = "RadialGeneralGaussian";
  521 +
  522 + // Initialise Radius
  523 + m_radius.clear();
  524 + m_radius.name("Radius");
  525 + m_radius.unit("deg");
  526 + m_radius.value(2.778e-4); // 1 arcsec
  527 + m_radius.min(2.778e-4); // 1 arcsec
  528 + m_radius.free();
  529 + m_radius.scale(1.0);
  530 + m_radius.gradient(0.0);
  531 + m_radius.has_grad(false);
  532 +
  533 + // Initialise Reciprocal index
  534 + m_ridx.clear();
  535 + m_ridx.name("R_Index");
  536 + m_ridx.value(0.5); // Gaussian
  537 + m_ridx.min(1.e-15); // need to be > 0
  538 + m_ridx.free();
  539 + m_ridx.scale(1.0);
  540 + m_ridx.gradient(0.0);
  541 + m_ridx.has_grad(false); // Radial components never have gradients
  542 +
  543 + // Set parameter pointer(s)
  544 + m_pars.push_back(&m_radius);
  545 + m_pars.push_back(&m_ridx);
  546 +
  547 + // Initialise precomputation cache. Note that zero values flag
  548 + // uninitialised as a zero radius and width shell is not meaningful
  549 + m_last_radius = 0.0;
  550 + m_inv_radius_rad = 0.0;
  551 + m_last_ridx = 0.0;
  552 + m_inv_ridx = 0.0;
  553 + m_value_norm = 0.0;
  554 +
  555 + // Return
  556 + return;
  557 +}
  558 +
  559 +
  560 +/***********************************************************************//**
  561 + * @brief Copy class members
  562 + *
  563 + * @param[in] model Generalised radial Gaussian spatial model.
  564 + *
  565 + * We do not have to push back the members on the parameter stack as this
  566 + * should have been done by init_members() that was called before. Otherwise
  567 + * we would have sigma twice on the stack.
  568 + ***************************************************************************/
  569 +void GModelSpatialRadialGeneralGauss::copy_members(const GModelSpatialRadialGeneralGauss& model)
  570 +{
  571 + // Copy members
  572 + m_type = model.m_type; // Needed to conserve model type
  573 + m_radius = model.m_radius;
  574 + m_ridx = model.m_ridx;
  575 +
  576 + // Copy pre-computation cache
  577 + m_last_radius = model.m_last_radius;
  578 + m_inv_radius_rad = model.m_inv_radius_rad;
  579 + m_last_ridx = model.m_last_ridx;
  580 + m_inv_ridx = model.m_inv_ridx;
  581 + m_value_norm = model.m_value_norm;
  582 +
  583 + // Return
  584 + return;
  585 +}
  586 +
  587 +
  588 +/***********************************************************************//**
  589 + * @brief Delete class members
  590 + ***************************************************************************/
  591 +void GModelSpatialRadialGeneralGauss::free_members(void)
  592 +{
  593 + // Return
  594 + return;
  595 +}
  596 +
  597 +
  598 +/***********************************************************************//**
  599 + * @brief Update precomputation cache
  600 + ***************************************************************************/
  601 +void GModelSpatialRadialGeneralGauss::update() const
  602 +{
  603 + // Update if radius has changed
  604 + if (m_last_radius != radius() or m_last_ridx != ridx()) {
  605 +
  606 + // Store last values
  607 + m_last_radius = radius();
  608 + m_last_ridx = ridx();
  609 +
  610 + // Compute radius in radians
  611 + double radius_rad = radius() * gammalib::deg2rad;
  612 + if (radius_rad > 0.0 and ridx() > 0.0) {
  613 + m_inv_radius_rad = 1.0 / radius_rad;
  614 + m_inv_ridx = 1.0 / ridx() ;
  615 + m_value_norm = 1.0 / (gammalib::twopi * radius_rad * radius_rad * ridx() *
  616 + std::exp(gammalib::gammln(2.0 * ridx())));
  617 + }
  618 + else {
  619 + m_inv_radius_rad = 0.0;
  620 + m_inv_ridx = 1.0;
  621 + m_value_norm = 0.0;
  622 + }
  623 +
  624 + } // endif: update required
  625 +
  626 + // Return
  627 + return;
  628 +}
  629 +
  630 +
  631 +/***********************************************************************//**
  632 + * @brief Set boundary sky region
  633 + ***************************************************************************/
  634 +void GModelSpatialRadialGeneralGauss::set_region(void) const
  635 +{
  636 + // Set sky region circle (5 times radius)
  637 + GSkyRegionCircle region(m_ra.value(), m_dec.value(), m_radius.value() * 5.0);
  638 +
  639 + // Set region (circumvent const correctness)
  640 + const_cast<GModelSpatialRadialGeneralGauss*>(this)->m_region = region;
  641 +
  642 + // Return
  643 + return;
  644 +}
... ...
src/model/Makefile.am
  View file @e9681ff
... ... @@ -40,6 +40,7 @@ sources = GModelPar.cpp \
40 40 GModelSpatialRadialDisk.cpp \
41 41 GModelSpatialRadialRing.cpp \
42 42 GModelSpatialRadialGauss.cpp \
  43 + GModelSpatialRadialGeneralGauss.cpp \
43 44 GModelSpatialRadialShell.cpp \
44 45 GModelSpatialRadialProfile.cpp \
45 46 GModelSpatialRadialProfileGauss.cpp \
... ... @@ -49,6 +50,7 @@ sources = GModelPar.cpp \
49 50 GModelSpatialElliptical.cpp \
50 51 GModelSpatialEllipticalDisk.cpp \
51 52 GModelSpatialEllipticalGauss.cpp \
  53 + GModelSpatialEllipticalGeneralGauss.cpp \
52 54 GModelSpatialDiffuse.cpp \
53 55 GModelSpatialDiffuseConst.cpp \
54 56 GModelSpatialDiffuseCube.cpp \
... ...
test/data/model_elliptical_general_gauss.xml 0 → 100644
  View file @e9681ff
  1 +<?xml version="1.0" standalone="no"?>
  2 +<source_library title="source library">
  3 + <source name="Cen A lobes" type="ExtendedSource">
  4 + <spectrum type="PowerLaw">
  5 + <parameter scale="1e-07" name="Prefactor" min="1e-07" max="1000.0" value="1.73" free="1"/>
  6 + <parameter scale="1.0" name="Index" min="-5.0" max="+5.0" value="-2.1" free="1"/>
  7 + <parameter scale="1.0" name="PivotEnergy" min="10.0" max="1000000.0" value="100.0" free="0"/>
  8 + </spectrum>
  9 + <spatialModel type="EllipticalGeneralGaussian">
  10 + <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
  11 + <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
  12 + <parameter name="PA" scale="1.0" value="45.0" min="-360" max="360" free="1"/>
  13 + <parameter name="MinorRadius" scale="1.0" value="0.1" min="0.001" max="10" free="1"/>
  14 + <parameter name="MajorRadius" scale="1.0" value="0.3" min="0.001" max="10" free="1"/>
  15 + <parameter name="R_Index" scale="1.0" value="0.3" min="0.01" max="10" free="1"/>
  16 + </spatialModel>
  17 + </source>
  18 +</source_library>
... ...
test/data/model_radial_general_gauss.xml 0 → 100644
  View file @e9681ff
  1 +<?xml version="1.0" standalone="no"?>
  2 +<source_library title="source library">
  3 + <source name="Cen A lobes" type="ExtendedSource">
  4 + <spectrum type="PowerLaw">
  5 + <parameter scale="1e-07" name="Prefactor" min="1e-07" max="1000.0" value="1.73" free="1"/>
  6 + <parameter scale="1.0" name="Index" min="-5.0" max="+5.0" value="-2.1" free="1"/>
  7 + <parameter scale="1.0" name="PivotEnergy" min="10.0" max="1000000.0" value="100.0" free="0"/>
  8 + </spectrum>
  9 + <spatialModel type="RadialGeneralGaussian">
  10 + <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
  11 + <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
  12 + <parameter name="Radius" scale="1.0" value="0.20" min="0.01" max="10" free="1"/>
  13 + <parameter name="R_Index" scale="1.0" value="0.3" min="0.01" max="10" free="1"/>
  14 + </spatialModel>
  15 + </source>
  16 +</source_library>
... ...
test/test_GModel.cpp
  View file @e9681ff
1 1 /***************************************************************************
2 2 * test_GModel.cpp - test GModel class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2021 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2022 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -84,37 +84,39 @@ void TestGModel::set(void)
84 84 m_assoc_file = datadir + "/crab_associations.xml";
85 85  
86 86 // Set model definiton XML files
87   - m_xml_model_point_const = datadir + "/model_point_const.xml";
88   - m_xml_model_point_gauss = datadir + "/model_point_gauss.xml";
89   - m_xml_model_point_plaw = datadir + "/model_point_plaw.xml";
90   - m_xml_model_point_plaw_phflux = datadir + "/model_point_plaw_phflux.xml";
91   - m_xml_model_point_plaw_eflux = datadir + "/model_point_plaw_eflux.xml";
92   - m_xml_model_point_eplaw = datadir + "/model_point_eplaw.xml";
93   - m_xml_model_point_einvplaw = datadir + "/model_point_einvplaw.xml";
94   - m_xml_model_point_bplaw = datadir + "/model_point_bplaw.xml";
95   - m_xml_model_point_smoothbplaw = datadir + "/model_point_smoothbplaw.xml";
96   - m_xml_model_point_supeplaw = datadir + "/model_point_supeplaw.xml";
97   - m_xml_model_point_logparabola = datadir + "/model_point_logparabola.xml";
98   - m_xml_model_point_bins = datadir + "/model_point_bins.xml";
99   - m_xml_model_point_nodes = datadir + "/model_point_nodes.xml";
100   - m_xml_model_point_filefct = datadir + "/model_point_filefct.xml";
101   - m_xml_model_point_table = datadir + "/model_point_table.xml";
102   - m_xml_point_multiplicative = datadir + "/model_point_multiplicative.xml";
103   - m_xml_point_exponential = datadir + "/model_point_exponential.xml";
104   - m_xml_model_spectral_composite = datadir + "/model_spectral_composite.xml";
105   - m_xml_model_diffuse_const = datadir + "/model_diffuse_const.xml";
106   - m_xml_model_diffuse_cube = datadir + "/model_diffuse_cube.xml";
107   - m_xml_model_diffuse_map = datadir + "/model_diffuse_map.xml";
108   - m_xml_model_radial_disk = datadir + "/model_radial_disk.xml";
109   - m_xml_model_radial_ring = datadir + "/model_radial_ring.xml";
110   - m_xml_model_radial_gauss = datadir + "/model_radial_gauss.xml";
111   - m_xml_model_radial_shell = datadir + "/model_radial_shell.xml";
112   - m_xml_model_profile_dmburkert = datadir + "/model_profile_dmburkert.xml";
113   - m_xml_model_profile_dmeinasto = datadir + "/model_profile_dmeinasto.xml";
114   - m_xml_model_profile_dmzhao = datadir + "/model_profile_dmzhao.xml";
115   - m_xml_model_elliptical_disk = datadir + "/model_elliptical_disk.xml";
116   - m_xml_model_elliptical_gauss = datadir + "/model_elliptical_gauss.xml";
117   - m_xml_model_spatial_composite = datadir + "/model_spatial_composite.xml";
  87 + m_xml_model_point_const = datadir + "/model_point_const.xml";
  88 + m_xml_model_point_gauss = datadir + "/model_point_gauss.xml";
  89 + m_xml_model_point_plaw = datadir + "/model_point_plaw.xml";
  90 + m_xml_model_point_plaw_phflux = datadir + "/model_point_plaw_phflux.xml";
  91 + m_xml_model_point_plaw_eflux = datadir + "/model_point_plaw_eflux.xml";
  92 + m_xml_model_point_eplaw = datadir + "/model_point_eplaw.xml";
  93 + m_xml_model_point_einvplaw = datadir + "/model_point_einvplaw.xml";
  94 + m_xml_model_point_bplaw = datadir + "/model_point_bplaw.xml";
  95 + m_xml_model_point_smoothbplaw = datadir + "/model_point_smoothbplaw.xml";
  96 + m_xml_model_point_supeplaw = datadir + "/model_point_supeplaw.xml";
  97 + m_xml_model_point_logparabola = datadir + "/model_point_logparabola.xml";
  98 + m_xml_model_point_bins = datadir + "/model_point_bins.xml";
  99 + m_xml_model_point_nodes = datadir + "/model_point_nodes.xml";
  100 + m_xml_model_point_filefct = datadir + "/model_point_filefct.xml";
  101 + m_xml_model_point_table = datadir + "/model_point_table.xml";
  102 + m_xml_point_multiplicative = datadir + "/model_point_multiplicative.xml";
  103 + m_xml_point_exponential = datadir + "/model_point_exponential.xml";
  104 + m_xml_model_spectral_composite = datadir + "/model_spectral_composite.xml";
  105 + m_xml_model_diffuse_const = datadir + "/model_diffuse_const.xml";
  106 + m_xml_model_diffuse_cube = datadir + "/model_diffuse_cube.xml";
  107 + m_xml_model_diffuse_map = datadir + "/model_diffuse_map.xml";
  108 + m_xml_model_radial_disk = datadir + "/model_radial_disk.xml";
  109 + m_xml_model_radial_ring = datadir + "/model_radial_ring.xml";
  110 + m_xml_model_radial_gauss = datadir + "/model_radial_gauss.xml";
  111 + m_xml_model_radial_general_gauss = datadir + "/model_radial_general_gauss.xml";
  112 + m_xml_model_radial_shell = datadir + "/model_radial_shell.xml";
  113 + m_xml_model_profile_dmburkert = datadir + "/model_profile_dmburkert.xml";
  114 + m_xml_model_profile_dmeinasto = datadir + "/model_profile_dmeinasto.xml";
  115 + m_xml_model_profile_dmzhao = datadir + "/model_profile_dmzhao.xml";
  116 + m_xml_model_elliptical_disk = datadir + "/model_elliptical_disk.xml";
  117 + m_xml_model_elliptical_gauss = datadir + "/model_elliptical_gauss.xml";
  118 + m_xml_model_elliptical_general_gauss = datadir + "/model_elliptical_general_gauss.xml";
  119 + m_xml_model_spatial_composite = datadir + "/model_spatial_composite.xml";
118 120  
119 121 // Set temporal model definition XML files
120 122 m_xml_model_point_temp_lightcurve = datadir + "/model_temporal_lightcurve.xml";
... ... @@ -155,6 +157,8 @@ void TestGModel::set(void)
155 157 "Test GModelSpatialRadialRing");
156 158 append(static_cast<pfunction>(&TestGModel::test_radial_gauss),
157 159 "Test GModelSpatialRadialGauss");
  160 + append(static_cast<pfunction>(&TestGModel::test_radial_general_gauss),
  161 + "Test GModelSpatialRadialGeneralGauss");
158 162 append(static_cast<pfunction>(&TestGModel::test_radial_shell),
159 163 "Test GModelSpatialRadialShell");
160 164 append(static_cast<pfunction>(&TestGModel::test_profile_dmburkert),
... ... @@ -167,6 +171,8 @@ void TestGModel::set(void)
167 171 "Test GModelSpatialEllipticalDisk");
168 172 append(static_cast<pfunction>(&TestGModel::test_elliptical_gauss),
169 173 "Test GModelSpatialEllipticalGauss");
  174 + append(static_cast<pfunction>(&TestGModel::test_elliptical_general_gauss),
  175 + "Test GModelSpatialEllipticalGeneralGauss");
170 176 append(static_cast<pfunction>(&TestGModel::test_diffuse_const),
171 177 "Test GModelSpatialDiffuseConst");
172 178 append(static_cast<pfunction>(&TestGModel::test_diffuse_cube),
... ... @@ -1398,6 +1404,80 @@ void TestGModel::test_radial_gauss(void)
1398 1404 return;
1399 1405 }
1400 1406  
  1407 +/***********************************************************************//**
  1408 + * @brief Test GModelSpatialRadialGauss class
  1409 + ***************************************************************************/
  1410 +void TestGModel::test_radial_general_gauss(void)
  1411 +{
  1412 + // Test void constructor
  1413 + GModelSpatialRadialGeneralGauss model1;
  1414 + test_value(model1.type(), "RadialGeneralGaussian");
  1415 +
  1416 + // Test value constructor
  1417 + GSkyDir dir1;
  1418 + dir1.radec_deg(83.6331, +22.0145);
  1419 + GModelSpatialRadialGeneralGauss model2(dir1, 3.0, 0.5);
  1420 + test_value(model2.ra(), 83.6331);
  1421 + test_value(model2.dec(), 22.0145);
  1422 + test_value(model2.radius(), 3.0);
  1423 + test_value(model2.ridx(), 0.5);
  1424 +
  1425 + // Test XML constructor
  1426 + GXml xml(m_xml_model_radial_general_gauss);
  1427 + GXmlElement* element = xml.element(0)->element(0)->element("spatialModel", 0);
  1428 + GModelSpatialRadialGeneralGauss model3(*element);
  1429 + test_value(model3.size(), 4);
  1430 + test_value(model3.type(), "RadialGeneralGaussian");
  1431 + test_value(model3.ra(), 83.6331);
  1432 + test_value(model3.dec(), 22.0145);
  1433 + test_value(model3.radius(), 0.20);
  1434 + test_value(model3.ridx(), 0.3);
  1435 +
  1436 + // Test ra method
  1437 + model3.ra(100.0);
  1438 + test_value(model3.ra(), 100.0);
  1439 +
  1440 + // Test dec method
  1441 + model3.dec(10.0);
  1442 + test_value(model3.dec(), 10.0);
  1443 +
  1444 + // Test dir method
  1445 + GSkyDir dir2;
  1446 + dir2.radec_deg(83.6331, +22.0145);
  1447 + model3.dir(dir2);
  1448 + test_assert(model3.dir() == dir2, "Test sky direction");
  1449 +
  1450 + // Test radius method
  1451 + model3.radius(3.9);
  1452 + test_value(model3.radius(), 3.9);
  1453 +
  1454 + // Test ridx method
  1455 + model3.ridx(0.76);
  1456 + test_value(model3.ridx(), 0.76);
  1457 +
  1458 + // Test region method
  1459 + GSkyDir dir3;
  1460 + dir3.radec_deg(83.6331, +22.0145);
  1461 + test_assert(model3.region()->contains(dir3), "Test region() method (inside)");
  1462 + dir3.radec_deg(83.6331, -70.0);
  1463 + test_assert(!model3.region()->contains(dir3), "Test region() method (outside)");
  1464 +
  1465 + // Test operator access
  1466 + const char* strarray[] = {"RA", "DEC", "Radius", "R_Index"};
  1467 + for (int i = 0; i < 4; ++i) {
  1468 + std::string keyname(strarray[i]);
  1469 + model3[keyname].value(2.1);
  1470 + model3[keyname].error(1.9);
  1471 + model3[keyname].gradient(0.8);
  1472 + test_value(model3[keyname].value(), 2.1);
  1473 + test_value(model3[keyname].error(), 1.9);
  1474 + test_value(model3[keyname].gradient(), 0.8);
  1475 + }
  1476 +
  1477 + // Exit test
  1478 + return;
  1479 +}
  1480 +
1401 1481  
1402 1482 /***********************************************************************//**
1403 1483 * @brief Test GModelSpatialRadialShell class
... ... @@ -1716,6 +1796,92 @@ void TestGModel::test_elliptical_gauss(void)
1716 1796 return;
1717 1797 }
1718 1798  
  1799 +/***********************************************************************//**
  1800 + * @brief Test GModelSpatialEllipticalGeneralGauss class
  1801 + ***************************************************************************/
  1802 +void TestGModel::test_elliptical_general_gauss(void)
  1803 +{
  1804 + // Test void constructor
  1805 + GModelSpatialEllipticalGeneralGauss model1;
  1806 + test_value(model1.type(), "EllipticalGeneralGaussian");
  1807 +
  1808 + // Test value constructor
  1809 + GSkyDir dir1;
  1810 + dir1.radec_deg(83.6331, +22.0145);
  1811 + GModelSpatialEllipticalGeneralGauss model2(dir1, 3.0, 2.0, 45.0,0.5);
  1812 + test_value(model2.ra(), 83.6331);
  1813 + test_value(model2.dec(), 22.0145);
  1814 + test_value(model2.posangle(), 45.0);
  1815 + test_value(model2.semimajor(), 3.0);
  1816 + test_value(model2.semiminor(), 2.0);
  1817 + test_value(model2.semiminor(), 2.0);
  1818 + test_value(model2.ridx(), 0.5);
  1819 +
  1820 + // Test XML constructor
  1821 + GXml xml(m_xml_model_elliptical_general_gauss);
  1822 + GXmlElement* element = xml.element(0)->element(0)->element("spatialModel", 0);
  1823 + GModelSpatialEllipticalGeneralGauss model3(*element);
  1824 + test_value(model3.size(), 6);
  1825 + test_value(model3.type(), "EllipticalGeneralGaussian");
  1826 + test_value(model3.ra(), 83.6331);
  1827 + test_value(model3.dec(), 22.0145);
  1828 + test_value(model3.posangle(), 45.0);
  1829 + test_value(model3.semimajor(), 0.3);
  1830 + test_value(model3.semiminor(), 0.1);
  1831 + test_value(model3.ridx(), 0.3);
  1832 +
  1833 + // Test ra method
  1834 + model3.ra(100.0);
  1835 + test_value(model3.ra(), 100.0);
  1836 +
  1837 + // Test dec method
  1838 + model3.dec(10.0);
  1839 + test_value(model3.dec(), 10.0);
  1840 +
  1841 + // Test dir method
  1842 + GSkyDir dir2;
  1843 + dir2.radec_deg(83.6331, +22.0145);
  1844 + model3.dir(dir2);
  1845 + test_assert(model3.dir() == dir2, "Test sky direction");
  1846 +
  1847 + // Test posangle method
  1848 + model3.posangle(3.9);
  1849 + test_value(model3.posangle(), 3.9);
  1850 +
  1851 + // Test semimajor method
  1852 + model3.semimajor(3.9);
  1853 + test_value(model3.semimajor(), 3.9);
  1854 +
  1855 + // Test semiminor method
  1856 + model3.semiminor(3.9);
  1857 + test_value(model3.semiminor(), 3.9);
  1858 +
  1859 + // Test ridx method
  1860 + model3.ridx(0.7);
  1861 + test_value(model3.ridx(), 0.7);
  1862 +
  1863 + // Test region method
  1864 + GSkyDir dir3;
  1865 + dir3.radec_deg(83.6331, +22.0145);
  1866 + test_assert(model3.region()->contains(dir3), "Test region() method (inside)");
  1867 + dir3.radec_deg(83.6331, +61.5);
  1868 + test_assert(!model3.region()->contains(dir3), "Test region() method (outside)");
  1869 +
  1870 + // Test operator access
  1871 + const char* strarray[] = {"RA", "DEC", "PA", "MinorRadius", "MajorRadius", "R_Index"};
  1872 + for (int i = 0; i < 6; ++i) {
  1873 + std::string keyname(strarray[i]);
  1874 + model3[keyname].value(2.1);
  1875 + model3[keyname].error(1.9);
  1876 + model3[keyname].gradient(0.8);
  1877 + test_value(model3[keyname].value(), 2.1);
  1878 + test_value(model3[keyname].error(), 1.9);
  1879 + test_value(model3[keyname].gradient(), 0.8);
  1880 + }
  1881 +
  1882 + // Exit test
  1883 + return;
  1884 +}
1719 1885  
1720 1886 /***********************************************************************//**
1721 1887 * @brief Test GModelSpatialComposite class
... ...
test/test_GModel.hpp
  View file @e9681ff
1 1 /***************************************************************************
2 2 * test_GModel.hpp - Test model module *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2021 by Jean-Baptiste Cayrou *
  4 + * copyright (C) 2012-2022 by Jean-Baptiste Cayrou *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -58,12 +58,14 @@ public:
58 58 void test_radial_disk(void);
59 59 void test_radial_ring(void);
60 60 void test_radial_gauss(void);
  61 + void test_radial_general_gauss();
61 62 void test_radial_shell(void);
62 63 void test_profile_dmburkert(void);
63 64 void test_profile_dmeinasto(void);
64 65 void test_profile_dmzhao(void);
65 66 void test_elliptical_disk(void);
66 67 void test_elliptical_gauss(void);
  68 + void test_elliptical_general_gauss(void);
67 69 void test_spatial_composite(void);
68 70 void test_spatial_model(void);
69 71 void test_const(void);
... ... @@ -147,12 +149,14 @@ private:
147 149 std::string m_xml_model_radial_disk;
148 150 std::string m_xml_model_radial_ring;
149 151 std::string m_xml_model_radial_gauss;
  152 + std::string m_xml_model_radial_general_gauss;
150 153 std::string m_xml_model_radial_shell;
151 154 std::string m_xml_model_profile_dmburkert;
152 155 std::string m_xml_model_profile_dmeinasto;
153 156 std::string m_xml_model_profile_dmzhao;
154 157 std::string m_xml_model_elliptical_disk;
155 158 std::string m_xml_model_elliptical_gauss;
  159 + std::string m_xml_model_elliptical_general_gauss;
156 160 std::string m_xml_model_spatial_composite;
157 161 std::string m_xml_model_point_temp_lightcurve;
158 162 std::string m_xml_model_point_temp_phasecurve;
... ...