diff --git a/ChangeLog b/ChangeLog index 6a29de2..61168d3 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,12 +1,14 @@ -2021-12-13 +2022-01-10 * Version 2.0.0 released ======================== Add GCOMObservation::tim() and GCOMObservation::oads() setter methods + Add GModelSpatialRadialGeneralGauss and + GModelSpatialEllipticalGeneralGauss classes (#3924) Remove GCOMSelection::zeta_min() and GCOMSelection::zeta_max() methods - Add GCOMSelection::fpmtflag(), GCOMSelection::use_d1() and + Add GCOMSelection::fpmtflag(), GCOMSelection::use_d1() and GCOMSelection::use_d2() methods (#3917) Add GCOMModelDRBPhibarBins class (#3913) Compute spectral model gradients for negative factor values diff --git a/NEWS b/NEWS index 713d527..a1b8b91 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in GammaLib 2.0.0 -13 December 2021 +10 January 2022 1. Introduction @@ -20,6 +20,8 @@ The following classes have been added: - GModelAssociation - GModelAssociations - GModelSpectralBins +- GModelSpatialRadialGeneralGauss +- GModelSpatialEllipticalGeneralGauss - GCTAModelSkyCube - GCTABackground2D - GCOMModelDRM @@ -207,6 +209,9 @@ Implemented GMatrixSparse::row() and GMatrixSparse::add_to_row() methods. 9. Model module --------------- +Added GModelSpatialRadialGeneralGauss and GModelSpatialEllipticalGeneralGauss +classes to model general Gaussian functions as spatial components (#3924). + Gradients for spectral model parameters are now also computed for all models in case that the factor values are negative. This provides now error bars for photon and energy flux power laws that are allowed to take negative flux values. diff --git a/README.md b/README.md index e3092f0..134432d 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ GammaLib information ==================== -* Version: 2.0.0.dev (13 December 2021) +* Version: 2.0.0.dev (10 January 2022) [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=gammalib-integrate-os)](https://cta-jenkins.irap.omp.eu/job/gammalib-integrate-os/) diff --git a/doc/source/admin/release_history/2.0.rst b/doc/source/admin/release_history/2.0.rst index c8d78b7..dba18db 100644 --- a/doc/source/admin/release_history/2.0.rst +++ b/doc/source/admin/release_history/2.0.rst @@ -1,6 +1,6 @@ .. _2.0: -GammaLib 2.0 release (3 December 2021) +GammaLib 2.0 release (10 January 2022) ====================================== Overview @@ -96,6 +96,8 @@ Improvements New features ------------ +* [`3924 `_] - + Add ``GModelSpatialRadialGeneralGauss`` and ``GModelSpatialEllipticalGeneralGauss`` classes * [`3910 `_] - Add ``GCOMModelDRM`` class * [`3606 `_] - diff --git a/include/GModelSpatialEllipticalGeneralGauss.hpp b/include/GModelSpatialEllipticalGeneralGauss.hpp new file mode 100644 index 0000000..61d016b --- /dev/null +++ b/include/GModelSpatialEllipticalGeneralGauss.hpp @@ -0,0 +1,160 @@ +/*************************************************************************** + * GModelSpatialEllipticalGeneralGauss.hpp - Generalised elliptical * + * Gaussian source model class * + * ----------------------------------------------------------------------- * + * copyright (C) 2022 by Luigi Tibaldo * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GModelSpatialEllipticalGeneralGauss.hpp + * @brief Generalised elliptical gaussian model class interface definition + * @author Luigi Tibaldo + */ + +#ifndef GMODELSPATIALELLIPTICALGENERALGAUSS_HPP +#define GMODELSPATIALELLIPTICALGENERALGAUSS_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include "GModelSpatialElliptical.hpp" +#include "GModelPar.hpp" +#include "GSkyDir.hpp" +#include "GSkyRegionCircle.hpp" +#include "GXmlElement.hpp" + + +/************************************************************************** + * @class GModelSpatialEllipticalGeneralGauss + * + * @brief Generalised elliptical gaussian source model class + * + * This class implements the spatial component of the factorised source + * model for an generalised elliptical Gaussian source, i.e. surface + * brightness according to an asymmetric generalised Gaussian. + ***************************************************************************/ +class GModelSpatialEllipticalGeneralGauss : public GModelSpatialElliptical { + +public: + // Constructors and destructors + GModelSpatialEllipticalGeneralGauss(void); + GModelSpatialEllipticalGeneralGauss(const GSkyDir& dir, + const double& major, + const double& minor, + const double& posangle, + const double& ridx); + explicit GModelSpatialEllipticalGeneralGauss(const GXmlElement& xml); + GModelSpatialEllipticalGeneralGauss(const GModelSpatialEllipticalGeneralGauss& model); + virtual ~GModelSpatialEllipticalGeneralGauss(void); + + // Operators + virtual GModelSpatialEllipticalGeneralGauss& operator=(const GModelSpatialEllipticalGeneralGauss& model); + + // Implemented pure virtual base class methods + virtual void clear(void); + virtual GModelSpatialEllipticalGeneralGauss* clone(void) const; + virtual std::string classname(void) const; + virtual double eval(const double& theta, + const double& posangle, + const GEnergy& energy, + const GTime& time, + const bool& gradients = false) const; + virtual GSkyDir mc(const GEnergy& energy, + const GTime& time, + GRan& ran) const; + virtual bool contains(const GSkyDir& dir, + const double& margin = 0.0) const; + virtual double theta_max(void) const; + virtual void read(const GXmlElement& xml); + virtual void write(GXmlElement& xml) const; + virtual std::string print(const GChatter& chatter = NORMAL) const; + + // Other methods + double ridx(void) const; + void ridx(const double& ridx); + + +protected: + // Protected methods + void init_members(void); + void copy_members(const GModelSpatialEllipticalGeneralGauss& model); + void free_members(void); + void update(void) const; + virtual void set_region(void) const; + + // Protected members + GModelPar m_ridx; //!< Reciprocal of exponent of the radial profile + + // Cached members used for pre-computations + mutable double m_last_minor; //!< Last semi-minor axis + mutable double m_last_major; //!< Last semi-major axis + mutable double m_last_posangle; //!< Last position angle + mutable double m_cospos2; //!< squared cosine of position angle + mutable double m_sinpos2; //!< squared sine of position angle + mutable double m_sin2pos; //!< sine of twice the position angle + mutable double m_minor2; //!< square of minor axis + mutable double m_major2; //!< square of major axis + mutable double m_minor_rad; //!< Minor axis in radians + mutable double m_major_rad; //!< Major axis in radians + mutable double m_norm; //!< Normalization + mutable double m_term1; //!< Help term 1 + mutable double m_term2; //!< Help term 2 + mutable double m_term3; //!< Help term 3 + mutable double m_last_ridx; //!< Last reciprocal radial index + mutable double m_inv_ridx; //!< Spatial profile index + +}; + + +/***********************************************************************//** + * @brief Return class name + * + * @return String containing the class name ("GModelSpatialEllipticalGeneralGauss"). + ***************************************************************************/ +inline +std::string GModelSpatialEllipticalGeneralGauss::classname(void) const +{ + return ("GModelSpatialEllipticalGeneralGauss"); +} + +/***********************************************************************//** + * @brief Return reciprocal of the elliptical profile index + * + * @return Reciprocal of the elliptical profile index. + * + * Returns the reciprocal of the elliptical profile index. + ***************************************************************************/ +inline +double GModelSpatialEllipticalGeneralGauss::ridx(void) const +{ + return (m_ridx.value()); +} + +/***********************************************************************//** + * @brief Set reciprocal of the elliptical profile index + * + * @param[in] ridx Reciprocal of the elliptical profile index. + * + * Sets the reciprocal of the elliptical profile index. + ***************************************************************************/ +inline +void GModelSpatialEllipticalGeneralGauss::ridx(const double& ridx) +{ + m_ridx.value(ridx); + return; +} + +#endif /* GMODELSPATIALELLIPTICALGENERALGAUSS_HPP */ diff --git a/include/GModelSpatialRadialGeneralGauss.hpp b/include/GModelSpatialRadialGeneralGauss.hpp new file mode 100644 index 0000000..a3b0c0b --- /dev/null +++ b/include/GModelSpatialRadialGeneralGauss.hpp @@ -0,0 +1,176 @@ +/*************************************************************************** + * GModelSpatialRadialGeneralGauss.hpp - Generalised radial Gaussian * + * source model class * + * ----------------------------------------------------------------------- * + * copyright (C) 2021-2022 by Luigi Tibaldo * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GModelSpatialRadialGeneralGauss.hpp + * @brief Generalized radial Gaussian model class interface definition + * @author Luigi Tibaldo + */ + +#ifndef GMODELSPATIALRADIALGENERALGAUSS_HPP +#define GMODELSPATIALRADIALGENERALGAUSS_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include "GModelSpatialRadial.hpp" +#include "GModelPar.hpp" +#include "GSkyDir.hpp" +#include "GSkyRegionCircle.hpp" +#include "GEnergy.hpp" +#include "GTime.hpp" +#include "GXmlElement.hpp" + + +/***********************************************************************//** + * @class GModelSpatialRadialGeneralGauss + * + * @brief Generalized radial Gaussian model class + * + * This class implements the spatial component of the factorised source + * model for a generalised radial Gaussian source. + ***************************************************************************/ +class GModelSpatialRadialGeneralGauss : public GModelSpatialRadial { + +public: + // Constructors and destructors + GModelSpatialRadialGeneralGauss(void); + GModelSpatialRadialGeneralGauss(const GSkyDir& dir, + const double& radius, + const double& ridx); + explicit GModelSpatialRadialGeneralGauss(const GXmlElement& xml); + GModelSpatialRadialGeneralGauss(const GModelSpatialRadialGeneralGauss& model); + virtual ~GModelSpatialRadialGeneralGauss(void); + + // Operators + virtual GModelSpatialRadialGeneralGauss& operator=(const GModelSpatialRadialGeneralGauss& model); + + // Implemented pure virtual methods + virtual void clear(void); + virtual GModelSpatialRadialGeneralGauss* clone(void) const; + virtual std::string classname(void) const; + virtual double eval(const double& theta, + const GEnergy& energy, + const GTime& time, + const bool& gradients = false) const; + virtual GSkyDir mc(const GEnergy& energy, + const GTime& time, + GRan& ran) const; + virtual bool contains(const GSkyDir& dir, + const double& margin = 0.0) const; + virtual double theta_max(void) const; + virtual void read(const GXmlElement& xml); + virtual void write(GXmlElement& xml) const; + virtual std::string print(const GChatter& chatter = NORMAL) const; + + // Other methods + double radius(void) const; + void radius(const double& radius); + double ridx(void) const; + void ridx(const double& ridx); + +protected: + // Protected methods + void init_members(void); + void copy_members(const GModelSpatialRadialGeneralGauss& model); + void free_members(void); + void update(void) const; + virtual void set_region(void) const; + + // Protected members + GModelPar m_radius; //!< Gaussian width (deg) + GModelPar m_ridx; //!< Reciprocal of exponent of the radial profile + + // Cached members used for pre-computations + mutable double m_last_radius; //!< Last radius + mutable double m_inv_radius_rad; //!< radius(rad)^-1 + mutable double m_last_ridx; //!< Last reciprocal radial index + mutable double m_inv_ridx; //!< Spatial profile index + mutable double m_value_norm; //!< 1/(2pi radius(rad)^2 ridx Gamma(ridx)) +}; + + +/***********************************************************************//** + * @brief Return class name + * + * @return String containing the class name ("GModelSpatialRadialGeneralGauss"). + ***************************************************************************/ +inline +std::string GModelSpatialRadialGeneralGauss::classname(void) const +{ + return ("GModelSpatialRadialGeneralGauss"); +} + + +/***********************************************************************//** + * @brief Return radius + * + * @return Radius (deg). + * + * Returns the radius in degrees. + ***************************************************************************/ +inline +double GModelSpatialRadialGeneralGauss::radius(void) const +{ + return (m_radius.value()); +} + +/***********************************************************************//** + * @brief Set radius + * + * @param[in] radius (deg). + * + * Sets the radius in degrees. + ***************************************************************************/ +inline +void GModelSpatialRadialGeneralGauss::radius(const double& radius) +{ + m_radius.value(radius); + return; +} + +/***********************************************************************//** + * @brief Return ridx + * + * @return Reciprocal of the radial profile index. + * + * Returns the reciprocal of the radial profile index. + ***************************************************************************/ +inline +double GModelSpatialRadialGeneralGauss::ridx(void) const +{ + return (m_ridx.value()); +} + +/***********************************************************************//** + * @brief Set reciprocal index + * + * @param[in] ridx Reciprocal of the radial profile index. + * + * Sets the reciprocal index of the radial profile. + ***************************************************************************/ +inline +void GModelSpatialRadialGeneralGauss::ridx(const double& ridx) +{ + m_ridx.value(ridx); + return; +} + +#endif /* GMODELSPATIALRADIALGENERALGAUSS_HPP */ diff --git a/include/GammaLib.hpp b/include/GammaLib.hpp index 2e75ffa..19915ed 100644 --- a/include/GammaLib.hpp +++ b/include/GammaLib.hpp @@ -213,6 +213,7 @@ #include "GModelSpatialRadial.hpp" #include "GModelSpatialRadialDisk.hpp" #include "GModelSpatialRadialGauss.hpp" +#include "GModelSpatialRadialGeneralGauss.hpp" #include "GModelSpatialRadialShell.hpp" #include "GModelSpatialRadialRing.hpp" #include "GModelSpatialRadialProfile.hpp" @@ -223,6 +224,7 @@ #include "GModelSpatialElliptical.hpp" #include "GModelSpatialEllipticalDisk.hpp" #include "GModelSpatialEllipticalGauss.hpp" +#include "GModelSpatialEllipticalGeneralGauss.hpp" #include "GModelSpatialDiffuse.hpp" #include "GModelSpatialDiffuseConst.hpp" #include "GModelSpatialDiffuseCube.hpp" diff --git a/include/Makefile.am b/include/Makefile.am index f7cbeed..ff83756 100644 --- a/include/Makefile.am +++ b/include/Makefile.am @@ -171,6 +171,7 @@ pkginclude_HEADERS = GammaLib.hpp \ GModelSpatialRadial.hpp \ GModelSpatialRadialDisk.hpp \ GModelSpatialRadialGauss.hpp \ + GModelSpatialRadialGeneralGauss.hpp \ GModelSpatialRadialShell.hpp \ GModelSpatialRadialRing.hpp \ GModelSpatialRadialProfile.hpp \ @@ -181,6 +182,7 @@ pkginclude_HEADERS = GammaLib.hpp \ GModelSpatialElliptical.hpp \ GModelSpatialEllipticalDisk.hpp \ GModelSpatialEllipticalGauss.hpp \ + GModelSpatialEllipticalGeneralGauss.hpp \ GModelSpatialDiffuse.hpp \ GModelSpatialDiffuseConst.hpp \ GModelSpatialDiffuseCube.hpp \ diff --git a/pyext/GModelSpatialEllipticalGeneralGauss.i b/pyext/GModelSpatialEllipticalGeneralGauss.i new file mode 100644 index 0000000..9bd722b --- /dev/null +++ b/pyext/GModelSpatialEllipticalGeneralGauss.i @@ -0,0 +1,98 @@ +/*************************************************************************** + * GModelSpatialEllipticalGeneralGauss.i - Generalised elliptical * + * Gaussian source model class * + * ----------------------------------------------------------------------- * + * copyright (C) 2022 by Luigi Tibaldo * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GModelSpatialEllipticalGeneralGauss.i + * @brief Elliptical generalised gauss model class interface definition + * @author Luigi Tibaldo + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GModelSpatialEllipticalGeneralGauss.hpp" +%} + + +/***********************************************************************//** + * @class GModelSpatialEllipticalGeneralGauss + * + * @brief Elliptical generalised gauss model class + ***************************************************************************/ +class GModelSpatialEllipticalGeneralGauss : public GModelSpatialElliptical { + +public: + // Constructors and destructors + GModelSpatialEllipticalGeneralGauss(void); + GModelSpatialEllipticalGeneralGauss(const GSkyDir& dir, + const double& major, + const double& minor, + const double& posangle, + const double& ridx); + explicit GModelSpatialEllipticalGeneralGauss(const GXmlElement& xml); + GModelSpatialEllipticalGeneralGauss(const GModelSpatialEllipticalGeneralGauss& model); + virtual ~GModelSpatialEllipticalGeneralGauss(void); + + // Implemented pure virtual base class methods + virtual void clear(void); + virtual GModelSpatialEllipticalGeneralGauss* clone(void) const; + virtual std::string classname(void) const; + virtual double eval(const double& theta, + const double& posangle, + const GEnergy& energy, + const GTime& time, + const bool& gradients = false) const; + virtual GSkyDir mc(const GEnergy& energy, + const GTime& time, + GRan& ran) const; + virtual bool contains(const GSkyDir& dir, + const double& margin = 0.0) const; + virtual double theta_max(void) const; + virtual void read(const GXmlElement& xml); + virtual void write(GXmlElement& xml) const; + + // Other methods + double ridx(void) const; + void ridx(const double& ridx); +}; + + +/***********************************************************************//** + * @brief GModelSpatialEllipticalGeneralGauss class extension + ***************************************************************************/ +%extend GModelSpatialEllipticalGeneralGauss { + GModelSpatialEllipticalGeneralGauss copy() { + return (*self); + } + double eval(const GPhoton& photon) const { + return self->GModelSpatialEllipticalGeneralGauss::eval(photon); + } + double eval(const GPhoton& photon, const bool& gradients) const { + return self->GModelSpatialEllipticalGeneralGauss::eval(photon, gradients); + } +%pythoncode { + def __getstate__(self): + xml = gammalib.GXmlElement() + self.write(xml) + state = (xml,) + return state + def __setstate__(self, state): + self.__init__(state[0]) +} +}; diff --git a/pyext/GModelSpatialRadialGeneralGauss.i b/pyext/GModelSpatialRadialGeneralGauss.i new file mode 100644 index 0000000..313dbcc --- /dev/null +++ b/pyext/GModelSpatialRadialGeneralGauss.i @@ -0,0 +1,103 @@ +/*************************************************************************** + * GModelSpatialRadialGeneralGauss.i - Generalised radial Gaussian * + * source model class * + * ----------------------------------------------------------------------- * + * copyright (C) 2021-2022 by Luigi Tibaldo * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GModelSpatialRadialGeneralGauss.i + * @brief Generalised radial Gaussian model class Python interface definition + * @author Luigi Tibaldo + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GModelSpatialRadialGeneralGauss.hpp" +%} + + +/***********************************************************************//** + * @class GModelSpatialRadialGeneralGauss + * + * @brief Generalised radial Gaussian model class + ***************************************************************************/ +class GModelSpatialRadialGeneralGauss : public GModelSpatialRadial { + +public: + // Constructors and destructors + GModelSpatialRadialGeneralGauss(void); + GModelSpatialRadialGeneralGauss(const GSkyDir& dir, + const double& radius, + const double& ridx); + explicit GModelSpatialRadialGeneralGauss(const GXmlElement& xml); + GModelSpatialRadialGeneralGauss(const GModelSpatialRadialGeneralGauss& model); + virtual ~GModelSpatialRadialGeneralGauss(void); + + // Implemented pure virtual methods + virtual void clear(void); + virtual GModelSpatialRadialGeneralGauss* clone(void) const; + virtual std::string classname(void) const; + virtual double eval(const double& theta, + const GEnergy& energy, + const GTime& time, + const bool& gradients = false) const; + virtual GSkyDir mc(const GEnergy& energy, + const GTime& time, + GRan& ran) const; + virtual bool contains(const GSkyDir& dir, + const double& margin = 0.0) const; + virtual double theta_max(void) const; + virtual void read(const GXmlElement& xml); + virtual void write(GXmlElement& xml) const; + + // Other methods + double radius(void) const; + double ridx(void) const; + void radius(const double& radius); + void ridx(const double& ridx); +}; + + +/***********************************************************************//** + * @brief GModelSpatialRadialGeneralGauss class extension + * + * The eval(GPhoton&, bool&) method need to be defined in the extension to + * force swig to build also the interface for these methods that are + * implemented in the base class only. It's not clear to me why these + * methods are not inherited automatically. Maybe this could also be handled + * by a %typemap(typecheck) construct. + ***************************************************************************/ +%extend GModelSpatialRadialGeneralGauss { + GModelSpatialRadialGeneralGauss copy() { + return (*self); + } + double eval(const GPhoton& photon) const { + return self->GModelSpatialRadialGeneralGauss::eval(photon); + } + double eval(const GPhoton& photon, const bool& gradients) const { + return self->GModelSpatialRadialGeneralGauss::eval(photon, gradients); + } +%pythoncode { + def __getstate__(self): + xml = gammalib.GXmlElement() + self.write(xml) + state = (xml,) + return state + def __setstate__(self, state): + self.__init__(state[0]) +} +}; diff --git a/pyext/gammalib/model.i b/pyext/gammalib/model.i index d941f43..2063e22 100644 --- a/pyext/gammalib/model.i +++ b/pyext/gammalib/model.i @@ -180,6 +180,7 @@ %include "GModelSpatialRadialDisk.i" %include "GModelSpatialRadialRing.i" %include "GModelSpatialRadialGauss.i" +%include "GModelSpatialRadialGeneralGauss.i" %include "GModelSpatialRadialShell.i" %include "GModelSpatialRadialProfile.i" %include "GModelSpatialRadialProfileGauss.i" @@ -189,6 +190,7 @@ %include "GModelSpatialElliptical.i" %include "GModelSpatialEllipticalDisk.i" %include "GModelSpatialEllipticalGauss.i" +%include "GModelSpatialEllipticalGeneralGauss.i" %include "GModelSpatialDiffuse.i" %include "GModelSpatialDiffuseConst.i" %include "GModelSpatialDiffuseCube.i" diff --git a/src/model/GModelSpatialEllipticalGeneralGauss.cpp b/src/model/GModelSpatialEllipticalGeneralGauss.cpp new file mode 100644 index 0000000..cbcf3f1 --- /dev/null +++ b/src/model/GModelSpatialEllipticalGeneralGauss.cpp @@ -0,0 +1,799 @@ +/*************************************************************************** + * GModelSpatialRadialGeneralGauss.cpp - Generalised radial Gaussian * + * source model class * + * ----------------------------------------------------------------------- * + * copyright (C) 2022 by Luigi Tibaldo * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GModelSpatialEllipticalGeneralGauss.cpp + * @brief Generalised elliptical Gaussian model class implementation + * @author Luigi Tibaldo + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include "GException.hpp" +#include "GTools.hpp" +#include "GMath.hpp" +#include "GModelSpatialEllipticalGeneralGauss.hpp" +#include "GModelSpatialRadialDisk.hpp" +#include "GModelSpatialRegistry.hpp" + +/* __ Constants __________________________________________________________ */ +namespace { + const double c_theta_max = 5.0; //!< max distance from center in units of major semiaxis +} + +/* __ Globals ____________________________________________________________ */ +const GModelSpatialEllipticalGeneralGauss g_elliptical_general_gauss_seed; +const GModelSpatialRegistry g_elliptical_general_gauss_registry(&g_elliptical_general_gauss_seed); + +/* __ Method name definitions ____________________________________________ */ +#define G_READ "GModelSpatialEllipticalGeneralGauss::read(GXmlElement&)" +#define G_WRITE "GModelSpatialEllipticalGeneralGauss::write(GXmlElement&)" +#define G_MC "GModelSpatialEllipticalGeneralGauss::mc(GEnergy&, GTime&, "\ + "GRan&)" + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ +#define G_SMALL_ANGLE_APPROXIMATION //!< Use small angle approximation + +/* __ Debug definitions __________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +GModelSpatialEllipticalGeneralGauss::GModelSpatialEllipticalGeneralGauss(void) : + GModelSpatialElliptical() +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Elliptical Gaussian constructor + * + * @param[in] dir Centre of elliptical Gaussian. + * @param[in] semimajor Semi-major axis (deg). + * @param[in] semiminor Semi-minor axis (deg). + * @param[in] posangle Position angle of semi-major axis (deg). + * @param[in] ridx Reciprocal of radial profile index. + * + * Construct elliptical Gaussian model from the ellipse centre direction + * (@p dir), the @p semimajor and @p semiminor axes, the position + * angle (@p posangle), and the reciprocal of the radial profile index + * (@p ridx). + ***************************************************************************/ +GModelSpatialEllipticalGeneralGauss::GModelSpatialEllipticalGeneralGauss(const GSkyDir& dir, + const double& semimajor, + const double& semiminor, + const double& posangle, + const double& ridx) : + GModelSpatialElliptical() +{ + // Initialise members + init_members(); + + // Assign parameters + this->dir(dir); + this->semiminor(semiminor); + this->semimajor(semimajor); + this->posangle(posangle); + this->ridx(ridx); + + // Return + return; +} + + +/***********************************************************************//** + * @brief XML constructor + * + * @param[in] xml XML element. + * + * Constructs generalised elliptical Gaussian model by extracting information + * from an XML element. See the read() method for more information about the + * expected structure of the XML element. + ***************************************************************************/ +GModelSpatialEllipticalGeneralGauss::GModelSpatialEllipticalGeneralGauss(const GXmlElement& xml) : + GModelSpatialElliptical() +{ + // Initialise members + init_members(); + + // Read information from XML element + read(xml); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] model Generalised elliptical Gaussian model. + ***************************************************************************/ +GModelSpatialEllipticalGeneralGauss::GModelSpatialEllipticalGeneralGauss(const GModelSpatialEllipticalGeneralGauss& model) : + GModelSpatialElliptical(model) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(model); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +GModelSpatialEllipticalGeneralGauss::~GModelSpatialEllipticalGeneralGauss(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] model Generalised elliptical Gaussian model. + * @return Generalised elliptical Gaussian model. + ***************************************************************************/ +GModelSpatialEllipticalGeneralGauss& GModelSpatialEllipticalGeneralGauss::operator=(const GModelSpatialEllipticalGeneralGauss& model) +{ + // Execute only if object is not identical + if (this != &model) { + + // Copy base class members + this->GModelSpatialElliptical::operator=(model); + + // Free members + free_members(); + + // Initialise members + init_members(); + + // Copy members + copy_members(model); + + } // endif: object was not identical + + // Return + return *this; +} + + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Clear generalised elliptical Gaussian model + ***************************************************************************/ +void GModelSpatialEllipticalGeneralGauss::clear(void) +{ + // Free class members (base and derived classes, derived class first) + free_members(); + this->GModelSpatialElliptical::free_members(); + this->GModelSpatial::free_members(); + + // Initialise members + this->GModelSpatial::init_members(); + this->GModelSpatialElliptical::init_members(); + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Clone generalised elliptical Gaussian model + * + * @return Pointer to deep copy of generalised elliptical Gaussian model. + ***************************************************************************/ +GModelSpatialEllipticalGeneralGauss* GModelSpatialEllipticalGeneralGauss::clone(void) const +{ + // Clone generalised elliptical Gaussian model + return new GModelSpatialEllipticalGeneralGauss(*this); +} + + +/***********************************************************************//** + * @brief Evaluate function (in units of sr^-1) + * + * @param[in] theta Angular distance from disk centre (radians). + * @param[in] posangle Position angle (counterclockwise from North) (radians). + * @param[in] energy Photon energy. + * @param[in] time Photon arrival time. + * @param[in] gradients Compute gradients? + * @return Model value. + * + * Evaluates the spatial component for a generalised elliptical Gaussian + * source model. + * + * The generalised elliptical Gaussian function is defined by + * \f$S_{\rm p}(\theta, \phi | E, t)\f$, where + * \f$\theta\f$ is the angular separation between elliptical Gaussian centre + * and the actual location and \f$\phi\f$ the position angle with respect to + * the model centre, counted counterclockwise from North. + * + * The function \f$S_{\rm p}(\theta, \phi | E, t)\f$ is given by + * + * \f[ + * S_{\rm p}(\theta, \phi | E, t) = {\tt m\_norm} \times + * \exp \left( -\left(\frac{\theta}{r_{\rm eff}}\right)^1/\eta \right) + * \f] + * + * where the effective ellipse radius \f$r_{\rm eff}\f$ towards a given + * position angle is given by + * + * \f[ + * r_{\rm eff} = \frac{ab} + * {\sqrt{\left( a \sin (\phi - \phi_0) \right)^2} + + * \sqrt{\left( b \cos (\phi - \phi_0) \right)^2}} + * \f] + * and + * \f$a\f$ is the semi-major axis of the ellipse, + * \f$b\f$ is the semi-minor axis, + * \f$\phi_0\f$ is the position angle of the ellipse, counted + * counterclockwise from North, and + * \f$\et\f$ is the reciprocal of the radial profile index + * \f${\tt m\_norm}\f$ is a normalization constant given by + * + * \f[ + * {\tt m\_norm} = \frac{1}{2 \pi a b \eta \Gamma(\eta)} + * \f] + * + * The method will not compute analytical parameter gradients, even if the + * @p gradients argument is set to true. Parameter gradients + * need to be computed numerically. + * + * @warning + * The above formula for an elliptical generalised Gaussian are accurate for small + * angles, with semimajor and semiminor axes below a few degrees. This + * should be fine for most use cases, but you should be aware that the + * ellipse will get distorted for larger angles. + * Note that the model normalisation is only correct in the small angle + * approximation and for \f$\eta\f$ order unity or smaller. + ***************************************************************************/ +double GModelSpatialEllipticalGeneralGauss::eval(const double& theta, + const double& posangle, + const GEnergy& energy, + const GTime& time, + const bool& gradients) const +{ + // Initialise value + double value = 0.0; + + // Continue only if we're inside circle enclosing the ellipse + if (theta <= theta_max()) { + + // Update precomputation cache + update(); + + // Perform computations to compute exponent + #if defined(G_SMALL_ANGLE_APPROXIMATION) + double rel_posangle = posangle - this->posangle() * gammalib::deg2rad; + double cosinus = std::cos(rel_posangle); + double sinus = std::sin(rel_posangle); + double arg1 = m_minor_rad * cosinus; + double arg2 = m_major_rad * sinus; + double r_ellipse = m_minor_rad * m_major_rad / + std::sqrt(arg1*arg1 + arg2*arg2); //!< small angle + double r_relative = theta/r_ellipse; + double exponent = std::pow(r_relative,m_inv_ridx); + #else + // Perform computations + double sinphi = std::sin(posangle); + double cosphi = std::cos(posangle); + + // Compute help term + double term1 = m_term1 * cosphi * cosphi; + double term2 = m_term2 * sinphi * sinphi; + double term3 = m_term3 * sinphi * cosphi; + + // Compute exponent + double exponent = std::pow(theta * std::sqrt(term1 + term2 + term3), m_inv_ridx) ; + #endif + + // Set value + value = m_norm * std::exp(-exponent); + + // Compile option: Check for NaN/Inf + #if defined(G_NAN_CHECK) + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { + std::cout << "*** ERROR: GModelSpatialEllipticalGeneralGauss::eval"; + std::cout << "(theta=" << theta << "): NaN/Inf encountered"; + std::cout << "(posangle=" << posangle << "): NaN/Inf encountered"; + std::cout << " (value=" << value; + std::cout << ", MinorAxis^2=" << m_minor2; + std::cout << ", MaxjorAxis^2=" << m_major2; + std::cout << ", m_norm=" << m_norm; + std::cout << ")" << std::endl; + } + #endif + + } // endif: position was inside enclosing circle + + // Return value + return value; +} + + +/***********************************************************************//** + * @brief Returns MC sky direction + * + * @param[in] energy Photon energy. + * @param[in] time Photon arrival time. + * @param[in,out] ran Random number generator. + * @return Sky direction. + * + * Draws an arbitrary sky direction from the elliptical generalised + * Gaussian model. + ***************************************************************************/ +GSkyDir GModelSpatialEllipticalGeneralGauss::mc(const GEnergy& energy, + const GTime& time, + GRan& ran) const +{ + // Update precomputation cache + update(); + + // Initialise simulation + #if defined(G_DEBUG_MC) + int n_samples = 0; + #endif + double u_max = m_norm * std::sin(theta_max()); + double value = 0.0; + double u = 1.0; + double theta = 0.0; + double phi = 0.0; + + // Throw an exception if the maximum value is zero + if (u_max <= 0.0) { + std::string msg = "Non positive maximum ellpitical general Gauss model value."; + throw GException::invalid_return_value(G_MC, msg); + } + + // Simulate offset from photon arrival direction using a rejection method + do { + theta = ran.uniform() * theta_max(); + phi = ran.uniform() * gammalib::twopi; + value = eval(theta, phi, energy, time) * std::sin(theta); + u = ran.uniform() * u_max; + #if defined(G_DEBUG_MC) + n_samples++; + #endif + } while (u > value); + #if defined(G_DEBUG_MC) + std::cout << "#=" << n_samples << " "; + #endif + + // Rotate pointing direction by offset and azimuth angle + GSkyDir sky_dir = dir(); + sky_dir.rotate_deg(phi * gammalib::rad2deg, theta * gammalib::rad2deg); + + // Return sky direction + return sky_dir; +} + + +/***********************************************************************//** + * @brief Checks whether model contains specified sky direction + * + * @param[in] dir Sky direction. + * @param[in] margin Margin to be added to sky direction (deg). + * @return True if the model contains the sky direction. + * + * Signals whether a sky direction is contained in the elliptical gauss + * model. + * + * @todo Implement correct evaluation of effective ellipse radius. + ***************************************************************************/ +bool GModelSpatialEllipticalGeneralGauss::contains(const GSkyDir& dir, + const double& margin) const +{ + // Compute distance to centre (radian) + double distance = dir.dist(this->dir()); + + // Return flag + return (distance <= theta_max() + margin*gammalib::deg2rad); +} + + +/***********************************************************************//** + * @brief Return maximum model radius (in radians) + * + * @return Returns maximum model radius (radians). + * + * Returns the maximum of 3.0 semimajor() and 3.0 semiminor() as + * approximate edge of the Gaussian. This limit is of course arbitrary, but + * allows to limit the integration region for response computation. The value + * of 3.0 has been determined by experiment (#1561). + ***************************************************************************/ +double GModelSpatialEllipticalGeneralGauss::theta_max(void) const +{ + // Set maximum model radius + double theta_max = (semimajor() > semiminor()) + ? semimajor() * gammalib::deg2rad * c_theta_max + : semiminor() * gammalib::deg2rad * c_theta_max; + + // Return value + return theta_max; +} + + +/***********************************************************************//** + * @brief Read model from XML element + * + * @param[in] xml XML element. + * + * Reads the elliptical gauss model information from an XML element. The XML + * element shall have either the format + * + * + * + * + * + * + * + * + * + * + * or + * + * + * + * + * + * + * + * + * + * + * @todo Implement a test of the ellipse boundary. The axes + * and axes minimum should be >0. + ***************************************************************************/ +void GModelSpatialEllipticalGeneralGauss::read(const GXmlElement& xml) +{ + // Verify number of model parameters + gammalib::xml_check_parnum(G_READ, xml, 6); + + // Read gauss location + GModelSpatialElliptical::read(xml); + + // Get parameters + const GXmlElement* minor = gammalib::xml_get_par(G_READ, xml, m_semiminor.name()); + const GXmlElement* major = gammalib::xml_get_par(G_READ, xml, m_semimajor.name()); + const GXmlElement* ridx = gammalib::xml_get_par(G_READ, xml, m_ridx.name()); + + // Read parameters + m_semiminor.read(*minor); + m_semimajor.read(*major); + m_ridx.read(*ridx); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Write model into XML element + * + * @param[in] xml XML element into which model information is written. + * + * Write the elliptical gauss model information into an XML element. The XML + * element will have the format + * + * + * + * + * + * + * + * + * + * + ***************************************************************************/ +void GModelSpatialEllipticalGeneralGauss::write(GXmlElement& xml) const +{ + // Verify model type + gammalib::xml_check_type(G_WRITE, xml, type()); + + // Write disk location + GModelSpatialElliptical::write(xml); + + // Get or create parameters + GXmlElement* minor = gammalib::xml_need_par(G_WRITE, xml, m_semiminor.name()); + GXmlElement* major = gammalib::xml_need_par(G_WRITE, xml, m_semimajor.name()); + GXmlElement* ridx = gammalib::xml_need_par(G_WRITE, xml, m_ridx.name()); + + // Write parameters + m_semiminor.write(*minor); + m_semimajor.write(*major); + m_ridx.write(*ridx); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Print information + * + * @param[in] chatter Chattiness. + * @return String containing model information. + ***************************************************************************/ +std::string GModelSpatialEllipticalGeneralGauss::print(const GChatter& chatter) const +{ + // Initialise result string + std::string result; + + // Continue only if chatter is not silent + if (chatter != SILENT) { + + // Append header + result.append("=== GModelSpatialEllipticalGeneralGauss ==="); + + // Append parameters + result.append("\n"+gammalib::parformat("Number of parameters")); + result.append(gammalib::str(size())); + for (int i = 0; i < size(); ++i) { + result.append("\n"+m_pars[i]->print(chatter)); + } + + } // endif: chatter was not silent + + // Return result + return result; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void GModelSpatialEllipticalGeneralGauss::init_members(void) +{ + // Initialise model type + m_type = "EllipticalGeneralGaussian"; + + // Initialise Reciprocal index + m_ridx.clear(); + m_ridx.name("R_Index"); + m_ridx.value(0.5); // Gaussian + m_ridx.min(1.e-15); // need to be > 0 + m_ridx.free(); + m_ridx.scale(1.0); + m_ridx.gradient(0.0); + m_ridx.has_grad(false); + + // Set parameter pointer(s) + m_pars.push_back(&m_ridx); + + // Initialise precomputation cache. Note that zero values flag + // uninitialised as a zero radius is not meaningful + m_last_minor = 0.0; + m_last_major = 0.0; + m_minor_rad = 0.0; + m_major_rad = 0.0; + m_norm = 0.0; + m_last_posangle = 9999.0; // Signals that has not been initialised + m_sin2pos = 0.0; + m_cospos2 = 0.0; + m_sinpos2 = 0.0; + m_minor2 = 0.0; + m_major2 = 0.0; + m_term1 = 0.0; + m_term2 = 0.0; + m_term3 = 0.0; + m_last_ridx = 0.0; + m_inv_ridx = 0.0; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] model Elliptical Gaussian model. + ***************************************************************************/ +void GModelSpatialEllipticalGeneralGauss::copy_members(const GModelSpatialEllipticalGeneralGauss& model) +{ + // Copy members + m_type = model.m_type; // Needed to conserve model type + m_ridx = model.m_ridx; + + // Copy precomputation cache + m_last_minor = model.m_last_minor; + m_last_major = model.m_last_major; + m_minor_rad = model.m_minor_rad; + m_major_rad = model.m_major_rad; + m_norm = model.m_norm; + m_last_posangle = model.m_last_posangle; + m_sin2pos = model.m_sin2pos; + m_cospos2 = model.m_cospos2; + m_sinpos2 = model.m_sinpos2; + m_minor2 = model.m_minor2; + m_major2 = model.m_major2; + m_term1 = model.m_term1; + m_term2 = model.m_term2; + m_term3 = model.m_term3; + m_last_ridx = model.m_last_ridx; + m_inv_ridx = model.m_inv_ridx; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void GModelSpatialEllipticalGeneralGauss::free_members(void) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief Update precomputation cache + * + * Computes the normalization + * \f[ + * {\tt m\_norm} = \frac{1}{2 \pi \times a \times b} + * \f] + * where + * \f$a\f$ is the semi-major axis and + * \f$b\f$ is the semi-minor axis of the ellipse. + * + * @warning + * The normalization of the elliptical Gaussian is only valid in the small + * angle approximation. + ***************************************************************************/ +void GModelSpatialEllipticalGeneralGauss::update() const +{ + // Initialise flag if something has changed + bool changed = false; + + // Update if one axis has changed or radial profile index has changed + if (m_last_minor != semiminor() || m_last_major != semimajor() || m_last_ridx != ridx()) { + + // Signal parameter changes + changed = true; + + // Store last values + m_last_minor = semiminor(); + m_last_major = semimajor(); + m_last_ridx = ridx(); + + // Compute axes in radians + m_minor_rad = m_last_minor * gammalib::deg2rad; + m_major_rad = m_last_major * gammalib::deg2rad; + + // Perform precomputations. Note also + // that the model normalisation is only correct in the small angle + // approximation and for ridx order unity or smaller. + if (m_minor_rad > 0.0 and m_major_rad > 0.0 and ridx() > 0.0) { + m_inv_ridx = 1.0 / ridx() ; + m_norm = 1.0 / (gammalib::twopi * m_minor_rad * m_major_rad * ridx() * + std::exp(gammalib::gammln(2.0 * ridx())) ); + } + else { + m_inv_ridx = 1.0; + m_norm = 0.0; + } + + } // endif: update required + + // Update chache if position angle changed + #if !defined(G_SMALL_ANGLE_APPROXIMATION) + if (m_last_posangle != posangle()) { + + // Signal parameter changes + changed = true; + + // Store last value + m_last_posangle = posangle(); + + // Compute angle in radians + double posangle_rad = m_last_posangle * gammalib::deg2rad; + + // Compute sine and cosine + double cospos = std::cos(posangle_rad); + double sinpos = std::sin(posangle_rad); + + // Cache important values for further computations + m_cospos2 = cospos * cospos; + m_sinpos2 = sinpos * sinpos; + m_sin2pos = std::sin(2.0 * posangle_rad); + + } // endif: position angle update required + + // Perform precomputations in case anything has changed + if (changed) { + + // Compute help terms + m_term1 = 0.5 * (m_cospos2 / m_minor2 + m_sinpos2 / m_major2); + m_term2 = 0.5 * (m_cospos2 / m_major2 + m_sinpos2 / m_minor2); + m_term3 = 0.5 * (m_sin2pos / m_major2 - m_sin2pos / m_minor2); + + } // endif: something has changed + #endif + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set boundary sky region + ***************************************************************************/ +void GModelSpatialEllipticalGeneralGauss::set_region(void) const +{ + // Set maximum model radius + double max_radius = (semimajor() > semiminor()) ? semimajor() : semiminor(); + + // Set sky region circle (maximum Gaussian sigma times a scaling + // factor (actually 5)) + GSkyRegionCircle region(m_ra.value(), m_dec.value(), max_radius * c_theta_max); + + // Set region (circumvent const correctness) + const_cast(this)->m_region = region; + + // Return + return; +} diff --git a/src/model/GModelSpatialRadialGeneralGauss.cpp b/src/model/GModelSpatialRadialGeneralGauss.cpp new file mode 100644 index 0000000..fe53b1c --- /dev/null +++ b/src/model/GModelSpatialRadialGeneralGauss.cpp @@ -0,0 +1,644 @@ +/*************************************************************************** + * GModelSpatialRadialGeneralGauss.cpp - Generalised radial Gaussian * + * source model class * + * ----------------------------------------------------------------------- * + * copyright (C) 2021-2022 by Luigi Tibaldo * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GModelSpatialRadialGeneralGauss.cpp + * @brief Radial Generallized Gaussian model class implementation + * @author Luigi Tibaldo + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include "GException.hpp" +#include "GTools.hpp" +#include "GMath.hpp" +#include "GModelSpatialRadialGeneralGauss.hpp" +#include "GModelSpatialRegistry.hpp" + +/* __ Constants __________________________________________________________ */ + +/* __ Globals ____________________________________________________________ */ +const GModelSpatialRadialGeneralGauss g_radial_general_gauss_seed; +const GModelSpatialRegistry g_radial_general_gauss_registry(&g_radial_general_gauss_seed); + +/* __ Method name definitions ____________________________________________ */ +#define G_EVAL "GModelSpatialRadialGeneralGauss::eval(double&, GEnergy&, "\ + "GTime&, bool&)" +#define G_READ "GModelSpatialRadialGeneralGauss::read(GXmlElement&)" +#define G_WRITE "GModelSpatialRadialGeneralGauss::write(GXmlElement&)" +#define G_MC "GModelSpatialRadialGeneralGauss::mc(GEnergy&, GTime&, GRan&)" + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + +/* __ Debug definitions __________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + * + * Constructs empty radial Gaussian model. + ***************************************************************************/ +GModelSpatialRadialGeneralGauss::GModelSpatialRadialGeneralGauss(void) : + GModelSpatialRadial() +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Constructor + * + * @param[in] dir Sky position of Gaussian. + * @param[in] radius Width of generalised Gaussian (deg). + * @param[in] ridx Reciprocal of exponential index of radial profile. + * + * Constructs a generalised Gaussian spatial model using a sky direction + * (@p dir), a Gaussian width parameter @p radius in degrees and a reciprocal + * of the exponential index @p ridx. + ***************************************************************************/ +GModelSpatialRadialGeneralGauss::GModelSpatialRadialGeneralGauss(const GSkyDir& dir, + const double& radius, + const double& ridx) : + GModelSpatialRadial() +{ + // Initialise members + init_members(); + + // Assign direction and radius + this->dir(dir); + this->radius(radius); + this->ridx(ridx); + + // Return + return; +} + + +/***********************************************************************//** + * @brief XML constructor + * + * @param[in] xml XML element. + * + * Constructs a generalised Gaussian spatial model by extracting information + * from an XML element. See the method read() for more information about the + * expected structure of the XML element. + ***************************************************************************/ +GModelSpatialRadialGeneralGauss::GModelSpatialRadialGeneralGauss(const GXmlElement& xml) : + GModelSpatialRadial() +{ + // Initialise members + init_members(); + + // Read information from XML element + read(xml); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] model Generalised radial Gaussian model. + ***************************************************************************/ +GModelSpatialRadialGeneralGauss::GModelSpatialRadialGeneralGauss(const GModelSpatialRadialGeneralGauss& model) : + GModelSpatialRadial(model) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(model); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +GModelSpatialRadialGeneralGauss::~GModelSpatialRadialGeneralGauss(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] model Generalised radial Gaussian model. + * @return Radial Generalised radial Gaussian model. + ***************************************************************************/ +GModelSpatialRadialGeneralGauss& GModelSpatialRadialGeneralGauss::operator=(const GModelSpatialRadialGeneralGauss& model) +{ + // Execute only if object is not identical + if (this != &model) { + + // Copy base class members + this->GModelSpatialRadial::operator=(model); + + // Free members + free_members(); + + // Initialise members + init_members(); + + // Copy members + copy_members(model); + + } // endif: object was not identical + + // Return + return *this; +} + + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Clear radial Gauss model + ***************************************************************************/ +void GModelSpatialRadialGeneralGauss::clear(void) +{ + // Free class members (base and derived classes, derived class first) + free_members(); + this->GModelSpatialRadial::free_members(); + this->GModelSpatial::free_members(); + + // Initialise members + this->GModelSpatial::init_members(); + this->GModelSpatialRadial::init_members(); + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Clone radial Gauss model + ***************************************************************************/ +GModelSpatialRadialGeneralGauss* GModelSpatialRadialGeneralGauss::clone(void) const +{ + // Clone radial Gauss model + return new GModelSpatialRadialGeneralGauss(*this); +} + + +/***********************************************************************//** + * @brief Evaluate Generalised Gaussian source model + * + * @param[in] theta Angular distance from source centre (radians). + * @param[in] energy Photon energy (not used). + * @param[in] time Photon arrival time (not used). + * @param[in] gradients Compute gradients? + * @return Model value (\f${\rm sr}^{-1}\f$). + * + * Evaluates the spatial component for a Generalised Gaussian source model using + * + * \f[ + * M_{\rm S}({\bf p} | E, t) = M_{\rm S}(\theta) = + * \frac{1}{2 \pi r^2 \eta \Gamma(\eta)} \exp + * \left(-(\frac{\theta}{r})^1/\eta \right) + * \f] + * + * where + * - \f$\theta\f$ is the angular separation from the centre of the model, and + * - \f$r\f$ is the Generalised Gaussian radius. + * - \f$\eta\f$ is the reciprocal of the radial profile exponent. + * + * If @p gradients is `true` the method will also compute parameter + * gradients + * + * + * Note that the model normalisation is only correct in the small angle + * approximation and for \f$\eta\f$ order unity or smaller. + ***************************************************************************/ +double GModelSpatialRadialGeneralGauss::eval(const double& theta, + const GEnergy& energy, + const GTime& time, + const bool& gradients) const +{ + // Update + update(); + + // Compute value + double arg = -1.0 * std::pow(m_inv_radius_rad * theta,m_inv_ridx); + double value = m_value_norm * std::exp(arg); + + // Compile option: Check for NaN/Inf + #if defined(G_NAN_CHECK) + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { + std::string msg = "Model value not a number:"; + for (int i = 0; i < m_pars.size(); ++i) { + msg += " " + m_pars[i]->name() + "="; + msg += gammalib::str(m_pars[i]->value()); + } + msg += " energy=" + energy.print(); + msg += " time=" + time.print(); + gammalib::warning(G_EVAL, msg); + } + #endif + + // Return value + return value; +} + + +/***********************************************************************//** + * @brief Returns MC sky direction + * + * @param[in] energy Photon energy. + * @param[in] time Photon arrival time. + * @param[in,out] ran Random number generator. + * @return Sky direction. + * + * Draws an arbitrary sky direction from the 2D Generalised Gaussian + * distribution as function of the photon @p energy and arrival @p time. + * + * @todo This method is only valid in the small angle approximation. + ***************************************************************************/ +GSkyDir GModelSpatialRadialGeneralGauss::mc(const GEnergy& energy, + const GTime& time, + GRan& ran) const +{ + + // Update precomputation cache + update(); + + // Initialise simulation + #if defined(G_DEBUG_MC) + int n_samples = 0; + #endif + double u_max = m_value_norm * std::sin(theta_max()); + double value = 0.0; + double u = 1.0; + double theta = 0.0; + + // Throw an exception if the maximum value is zero + if (u_max <= 0.0) { + std::string msg = "Non positive maximum radial genral Gauss model value."; + throw GException::invalid_return_value(G_MC, msg); + } + + // Simulate offset from photon arrival direction using a rejection method + do { + theta = ran.uniform() * theta_max(); + value = eval(theta, energy, time) * std::sin(theta); + u = ran.uniform() * u_max; + #if defined(G_DEBUG_MC) + n_samples++; + #endif + } while (u > value); + #if defined(G_DEBUG_MC) + std::cout << "#=" << n_samples << " "; + #endif + + // Simulate azimuth angle + double phi = 360.0 * ran.uniform(); + + // Rotate pointing direction by offset and azimuth angle + GSkyDir sky_dir = dir(); + sky_dir.rotate_deg(phi, theta * gammalib::rad2deg); + + // Return sky direction + return sky_dir; + +} + + +/***********************************************************************//** + * @brief Checks where model contains specified sky direction + * + * @param[in] dir Sky direction. + * @param[in] margin Margin to be added to sky direction (deg) + * @return True if the model contains the sky direction. + * + * Signals whether a sky direction is contained in the radial gauss model. + ***************************************************************************/ +bool GModelSpatialRadialGeneralGauss::contains(const GSkyDir& dir, + const double& margin) const +{ + // Compute distance to centre (radians) + double distance = dir.dist(this->dir()); + + // Return flag + return (distance <= theta_max() + margin*gammalib::deg2rad); +} + + +/***********************************************************************//** + * @brief Return maximum model radius (in radians) + * + * Returns \f$5 \sigma\f$ as approximate edge of the Gaussian. This limit + * is of course arbitrary, but allows to limit the integration region for + * response computation. + ***************************************************************************/ +double GModelSpatialRadialGeneralGauss::theta_max(void) const +{ + // Return value + return (radius() * gammalib::deg2rad * 5.0); +} + + +/***********************************************************************//** + * @brief Read model from XML element + * + * @param[in] xml XML element. + * + * Reads the generalised radial Gaussian model information from an XML + * element. The XML element shall have either the format + * + * + * + * + * + * + * + * + * or + * + * + * + * + * + * + * + * + * @todo Implement a test of the radius and radius boundary. The sigma + * and sigma minimum should be >0. + ***************************************************************************/ +void GModelSpatialRadialGeneralGauss::read(const GXmlElement& xml) +{ + // Verify number of model parameters + gammalib::xml_check_parnum(G_READ, xml, 4); + + // Read location + GModelSpatialRadial::read(xml); + + // Get parameters + const GXmlElement* radius = gammalib::xml_get_par(G_READ, xml, m_radius.name()); + const GXmlElement* ridx = gammalib::xml_get_par(G_READ, xml, m_ridx.name()); + + // Read parameters + m_radius.read(*radius); + m_ridx.read(*ridx); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Write model into XML element + * + * @param[in] xml XML element into which model information is written. + * + * Writes the radial disk model information into an XML element. The XML + * element will have the format + * + * + * + * + * + * + * + * + ***************************************************************************/ +void GModelSpatialRadialGeneralGauss::write(GXmlElement& xml) const +{ + // Write Gaussian location + GModelSpatialRadial::write(xml); + + // Get or create parameters + GXmlElement* radius = gammalib::xml_need_par(G_WRITE, xml, m_radius.name()); + GXmlElement* ridx = gammalib::xml_need_par(G_WRITE, xml, m_ridx.name()); + + // Write parameters + m_radius.write(*radius); + m_ridx.write(*ridx); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Print generalised radial Gaussian source information + * + * @param[in] chatter Chattiness. + * @return String containing model information. + ***************************************************************************/ +std::string GModelSpatialRadialGeneralGauss::print(const GChatter& chatter) const +{ + // Initialise result string + std::string result; + + // Continue only if chatter is not silent + if (chatter != SILENT) { + + // Append header + result.append("=== GModelSpatialRadialGeneralGauss ==="); + + // Append parameters + result.append("\n"+gammalib::parformat("Number of parameters")); + result.append(gammalib::str(size())); + for (int i = 0; i < size(); ++i) { + result.append("\n"+m_pars[i]->print(chatter)); + } + + } // endif: chatter was not silent + + // Return result + return result; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + * + * Note that the minimum Gaussian width is set to 1 arcsec + * and w set the minimum reciprocal index to 1.e-15 + ***************************************************************************/ +void GModelSpatialRadialGeneralGauss::init_members(void) +{ + // Initialise model type + m_type = "RadialGeneralGaussian"; + + // Initialise Radius + m_radius.clear(); + m_radius.name("Radius"); + m_radius.unit("deg"); + m_radius.value(2.778e-4); // 1 arcsec + m_radius.min(2.778e-4); // 1 arcsec + m_radius.free(); + m_radius.scale(1.0); + m_radius.gradient(0.0); + m_radius.has_grad(false); + + // Initialise Reciprocal index + m_ridx.clear(); + m_ridx.name("R_Index"); + m_ridx.value(0.5); // Gaussian + m_ridx.min(1.e-15); // need to be > 0 + m_ridx.free(); + m_ridx.scale(1.0); + m_ridx.gradient(0.0); + m_ridx.has_grad(false); // Radial components never have gradients + + // Set parameter pointer(s) + m_pars.push_back(&m_radius); + m_pars.push_back(&m_ridx); + + // Initialise precomputation cache. Note that zero values flag + // uninitialised as a zero radius and width shell is not meaningful + m_last_radius = 0.0; + m_inv_radius_rad = 0.0; + m_last_ridx = 0.0; + m_inv_ridx = 0.0; + m_value_norm = 0.0; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] model Generalised radial Gaussian spatial model. + * + * We do not have to push back the members on the parameter stack as this + * should have been done by init_members() that was called before. Otherwise + * we would have sigma twice on the stack. + ***************************************************************************/ +void GModelSpatialRadialGeneralGauss::copy_members(const GModelSpatialRadialGeneralGauss& model) +{ + // Copy members + m_type = model.m_type; // Needed to conserve model type + m_radius = model.m_radius; + m_ridx = model.m_ridx; + + // Copy pre-computation cache + m_last_radius = model.m_last_radius; + m_inv_radius_rad = model.m_inv_radius_rad; + m_last_ridx = model.m_last_ridx; + m_inv_ridx = model.m_inv_ridx; + m_value_norm = model.m_value_norm; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void GModelSpatialRadialGeneralGauss::free_members(void) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief Update precomputation cache + ***************************************************************************/ +void GModelSpatialRadialGeneralGauss::update() const +{ + // Update if radius has changed + if (m_last_radius != radius() or m_last_ridx != ridx()) { + + // Store last values + m_last_radius = radius(); + m_last_ridx = ridx(); + + // Compute radius in radians + double radius_rad = radius() * gammalib::deg2rad; + if (radius_rad > 0.0 and ridx() > 0.0) { + m_inv_radius_rad = 1.0 / radius_rad; + m_inv_ridx = 1.0 / ridx() ; + m_value_norm = 1.0 / (gammalib::twopi * radius_rad * radius_rad * ridx() * + std::exp(gammalib::gammln(2.0 * ridx()))); + } + else { + m_inv_radius_rad = 0.0; + m_inv_ridx = 1.0; + m_value_norm = 0.0; + } + + } // endif: update required + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set boundary sky region + ***************************************************************************/ +void GModelSpatialRadialGeneralGauss::set_region(void) const +{ + // Set sky region circle (5 times radius) + GSkyRegionCircle region(m_ra.value(), m_dec.value(), m_radius.value() * 5.0); + + // Set region (circumvent const correctness) + const_cast(this)->m_region = region; + + // Return + return; +} diff --git a/src/model/Makefile.am b/src/model/Makefile.am index 912a8cd..b95a8b4 100644 --- a/src/model/Makefile.am +++ b/src/model/Makefile.am @@ -40,6 +40,7 @@ sources = GModelPar.cpp \ GModelSpatialRadialDisk.cpp \ GModelSpatialRadialRing.cpp \ GModelSpatialRadialGauss.cpp \ + GModelSpatialRadialGeneralGauss.cpp \ GModelSpatialRadialShell.cpp \ GModelSpatialRadialProfile.cpp \ GModelSpatialRadialProfileGauss.cpp \ @@ -49,6 +50,7 @@ sources = GModelPar.cpp \ GModelSpatialElliptical.cpp \ GModelSpatialEllipticalDisk.cpp \ GModelSpatialEllipticalGauss.cpp \ + GModelSpatialEllipticalGeneralGauss.cpp \ GModelSpatialDiffuse.cpp \ GModelSpatialDiffuseConst.cpp \ GModelSpatialDiffuseCube.cpp \ diff --git a/test/data/model_elliptical_general_gauss.xml b/test/data/model_elliptical_general_gauss.xml new file mode 100644 index 0000000..8cacd9f --- /dev/null +++ b/test/data/model_elliptical_general_gauss.xml @@ -0,0 +1,18 @@ + + + + + + + + + + + + + + + + + + diff --git a/test/data/model_radial_general_gauss.xml b/test/data/model_radial_general_gauss.xml new file mode 100644 index 0000000..eee156f --- /dev/null +++ b/test/data/model_radial_general_gauss.xml @@ -0,0 +1,16 @@ + + + + + + + + + + + + + + + + diff --git a/test/test_GModel.cpp b/test/test_GModel.cpp index 7a02a77..c31ef15 100644 --- a/test/test_GModel.cpp +++ b/test/test_GModel.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * test_GModel.cpp - test GModel class * * ----------------------------------------------------------------------- * - * copyright (C) 2009-2021 by Juergen Knoedlseder * + * copyright (C) 2009-2022 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -84,37 +84,39 @@ void TestGModel::set(void) m_assoc_file = datadir + "/crab_associations.xml"; // Set model definiton XML files - m_xml_model_point_const = datadir + "/model_point_const.xml"; - m_xml_model_point_gauss = datadir + "/model_point_gauss.xml"; - m_xml_model_point_plaw = datadir + "/model_point_plaw.xml"; - m_xml_model_point_plaw_phflux = datadir + "/model_point_plaw_phflux.xml"; - m_xml_model_point_plaw_eflux = datadir + "/model_point_plaw_eflux.xml"; - m_xml_model_point_eplaw = datadir + "/model_point_eplaw.xml"; - m_xml_model_point_einvplaw = datadir + "/model_point_einvplaw.xml"; - m_xml_model_point_bplaw = datadir + "/model_point_bplaw.xml"; - m_xml_model_point_smoothbplaw = datadir + "/model_point_smoothbplaw.xml"; - m_xml_model_point_supeplaw = datadir + "/model_point_supeplaw.xml"; - m_xml_model_point_logparabola = datadir + "/model_point_logparabola.xml"; - m_xml_model_point_bins = datadir + "/model_point_bins.xml"; - m_xml_model_point_nodes = datadir + "/model_point_nodes.xml"; - m_xml_model_point_filefct = datadir + "/model_point_filefct.xml"; - m_xml_model_point_table = datadir + "/model_point_table.xml"; - m_xml_point_multiplicative = datadir + "/model_point_multiplicative.xml"; - m_xml_point_exponential = datadir + "/model_point_exponential.xml"; - m_xml_model_spectral_composite = datadir + "/model_spectral_composite.xml"; - m_xml_model_diffuse_const = datadir + "/model_diffuse_const.xml"; - m_xml_model_diffuse_cube = datadir + "/model_diffuse_cube.xml"; - m_xml_model_diffuse_map = datadir + "/model_diffuse_map.xml"; - m_xml_model_radial_disk = datadir + "/model_radial_disk.xml"; - m_xml_model_radial_ring = datadir + "/model_radial_ring.xml"; - m_xml_model_radial_gauss = datadir + "/model_radial_gauss.xml"; - m_xml_model_radial_shell = datadir + "/model_radial_shell.xml"; - m_xml_model_profile_dmburkert = datadir + "/model_profile_dmburkert.xml"; - m_xml_model_profile_dmeinasto = datadir + "/model_profile_dmeinasto.xml"; - m_xml_model_profile_dmzhao = datadir + "/model_profile_dmzhao.xml"; - m_xml_model_elliptical_disk = datadir + "/model_elliptical_disk.xml"; - m_xml_model_elliptical_gauss = datadir + "/model_elliptical_gauss.xml"; - m_xml_model_spatial_composite = datadir + "/model_spatial_composite.xml"; + m_xml_model_point_const = datadir + "/model_point_const.xml"; + m_xml_model_point_gauss = datadir + "/model_point_gauss.xml"; + m_xml_model_point_plaw = datadir + "/model_point_plaw.xml"; + m_xml_model_point_plaw_phflux = datadir + "/model_point_plaw_phflux.xml"; + m_xml_model_point_plaw_eflux = datadir + "/model_point_plaw_eflux.xml"; + m_xml_model_point_eplaw = datadir + "/model_point_eplaw.xml"; + m_xml_model_point_einvplaw = datadir + "/model_point_einvplaw.xml"; + m_xml_model_point_bplaw = datadir + "/model_point_bplaw.xml"; + m_xml_model_point_smoothbplaw = datadir + "/model_point_smoothbplaw.xml"; + m_xml_model_point_supeplaw = datadir + "/model_point_supeplaw.xml"; + m_xml_model_point_logparabola = datadir + "/model_point_logparabola.xml"; + m_xml_model_point_bins = datadir + "/model_point_bins.xml"; + m_xml_model_point_nodes = datadir + "/model_point_nodes.xml"; + m_xml_model_point_filefct = datadir + "/model_point_filefct.xml"; + m_xml_model_point_table = datadir + "/model_point_table.xml"; + m_xml_point_multiplicative = datadir + "/model_point_multiplicative.xml"; + m_xml_point_exponential = datadir + "/model_point_exponential.xml"; + m_xml_model_spectral_composite = datadir + "/model_spectral_composite.xml"; + m_xml_model_diffuse_const = datadir + "/model_diffuse_const.xml"; + m_xml_model_diffuse_cube = datadir + "/model_diffuse_cube.xml"; + m_xml_model_diffuse_map = datadir + "/model_diffuse_map.xml"; + m_xml_model_radial_disk = datadir + "/model_radial_disk.xml"; + m_xml_model_radial_ring = datadir + "/model_radial_ring.xml"; + m_xml_model_radial_gauss = datadir + "/model_radial_gauss.xml"; + m_xml_model_radial_general_gauss = datadir + "/model_radial_general_gauss.xml"; + m_xml_model_radial_shell = datadir + "/model_radial_shell.xml"; + m_xml_model_profile_dmburkert = datadir + "/model_profile_dmburkert.xml"; + m_xml_model_profile_dmeinasto = datadir + "/model_profile_dmeinasto.xml"; + m_xml_model_profile_dmzhao = datadir + "/model_profile_dmzhao.xml"; + m_xml_model_elliptical_disk = datadir + "/model_elliptical_disk.xml"; + m_xml_model_elliptical_gauss = datadir + "/model_elliptical_gauss.xml"; + m_xml_model_elliptical_general_gauss = datadir + "/model_elliptical_general_gauss.xml"; + m_xml_model_spatial_composite = datadir + "/model_spatial_composite.xml"; // Set temporal model definition XML files m_xml_model_point_temp_lightcurve = datadir + "/model_temporal_lightcurve.xml"; @@ -155,6 +157,8 @@ void TestGModel::set(void) "Test GModelSpatialRadialRing"); append(static_cast(&TestGModel::test_radial_gauss), "Test GModelSpatialRadialGauss"); + append(static_cast(&TestGModel::test_radial_general_gauss), + "Test GModelSpatialRadialGeneralGauss"); append(static_cast(&TestGModel::test_radial_shell), "Test GModelSpatialRadialShell"); append(static_cast(&TestGModel::test_profile_dmburkert), @@ -167,6 +171,8 @@ void TestGModel::set(void) "Test GModelSpatialEllipticalDisk"); append(static_cast(&TestGModel::test_elliptical_gauss), "Test GModelSpatialEllipticalGauss"); + append(static_cast(&TestGModel::test_elliptical_general_gauss), + "Test GModelSpatialEllipticalGeneralGauss"); append(static_cast(&TestGModel::test_diffuse_const), "Test GModelSpatialDiffuseConst"); append(static_cast(&TestGModel::test_diffuse_cube), @@ -1398,6 +1404,80 @@ void TestGModel::test_radial_gauss(void) return; } +/***********************************************************************//** + * @brief Test GModelSpatialRadialGauss class + ***************************************************************************/ +void TestGModel::test_radial_general_gauss(void) +{ + // Test void constructor + GModelSpatialRadialGeneralGauss model1; + test_value(model1.type(), "RadialGeneralGaussian"); + + // Test value constructor + GSkyDir dir1; + dir1.radec_deg(83.6331, +22.0145); + GModelSpatialRadialGeneralGauss model2(dir1, 3.0, 0.5); + test_value(model2.ra(), 83.6331); + test_value(model2.dec(), 22.0145); + test_value(model2.radius(), 3.0); + test_value(model2.ridx(), 0.5); + + // Test XML constructor + GXml xml(m_xml_model_radial_general_gauss); + GXmlElement* element = xml.element(0)->element(0)->element("spatialModel", 0); + GModelSpatialRadialGeneralGauss model3(*element); + test_value(model3.size(), 4); + test_value(model3.type(), "RadialGeneralGaussian"); + test_value(model3.ra(), 83.6331); + test_value(model3.dec(), 22.0145); + test_value(model3.radius(), 0.20); + test_value(model3.ridx(), 0.3); + + // Test ra method + model3.ra(100.0); + test_value(model3.ra(), 100.0); + + // Test dec method + model3.dec(10.0); + test_value(model3.dec(), 10.0); + + // Test dir method + GSkyDir dir2; + dir2.radec_deg(83.6331, +22.0145); + model3.dir(dir2); + test_assert(model3.dir() == dir2, "Test sky direction"); + + // Test radius method + model3.radius(3.9); + test_value(model3.radius(), 3.9); + + // Test ridx method + model3.ridx(0.76); + test_value(model3.ridx(), 0.76); + + // Test region method + GSkyDir dir3; + dir3.radec_deg(83.6331, +22.0145); + test_assert(model3.region()->contains(dir3), "Test region() method (inside)"); + dir3.radec_deg(83.6331, -70.0); + test_assert(!model3.region()->contains(dir3), "Test region() method (outside)"); + + // Test operator access + const char* strarray[] = {"RA", "DEC", "Radius", "R_Index"}; + for (int i = 0; i < 4; ++i) { + std::string keyname(strarray[i]); + model3[keyname].value(2.1); + model3[keyname].error(1.9); + model3[keyname].gradient(0.8); + test_value(model3[keyname].value(), 2.1); + test_value(model3[keyname].error(), 1.9); + test_value(model3[keyname].gradient(), 0.8); + } + + // Exit test + return; +} + /***********************************************************************//** * @brief Test GModelSpatialRadialShell class @@ -1716,6 +1796,92 @@ void TestGModel::test_elliptical_gauss(void) return; } +/***********************************************************************//** + * @brief Test GModelSpatialEllipticalGeneralGauss class + ***************************************************************************/ +void TestGModel::test_elliptical_general_gauss(void) +{ + // Test void constructor + GModelSpatialEllipticalGeneralGauss model1; + test_value(model1.type(), "EllipticalGeneralGaussian"); + + // Test value constructor + GSkyDir dir1; + dir1.radec_deg(83.6331, +22.0145); + GModelSpatialEllipticalGeneralGauss model2(dir1, 3.0, 2.0, 45.0,0.5); + test_value(model2.ra(), 83.6331); + test_value(model2.dec(), 22.0145); + test_value(model2.posangle(), 45.0); + test_value(model2.semimajor(), 3.0); + test_value(model2.semiminor(), 2.0); + test_value(model2.semiminor(), 2.0); + test_value(model2.ridx(), 0.5); + + // Test XML constructor + GXml xml(m_xml_model_elliptical_general_gauss); + GXmlElement* element = xml.element(0)->element(0)->element("spatialModel", 0); + GModelSpatialEllipticalGeneralGauss model3(*element); + test_value(model3.size(), 6); + test_value(model3.type(), "EllipticalGeneralGaussian"); + test_value(model3.ra(), 83.6331); + test_value(model3.dec(), 22.0145); + test_value(model3.posangle(), 45.0); + test_value(model3.semimajor(), 0.3); + test_value(model3.semiminor(), 0.1); + test_value(model3.ridx(), 0.3); + + // Test ra method + model3.ra(100.0); + test_value(model3.ra(), 100.0); + + // Test dec method + model3.dec(10.0); + test_value(model3.dec(), 10.0); + + // Test dir method + GSkyDir dir2; + dir2.radec_deg(83.6331, +22.0145); + model3.dir(dir2); + test_assert(model3.dir() == dir2, "Test sky direction"); + + // Test posangle method + model3.posangle(3.9); + test_value(model3.posangle(), 3.9); + + // Test semimajor method + model3.semimajor(3.9); + test_value(model3.semimajor(), 3.9); + + // Test semiminor method + model3.semiminor(3.9); + test_value(model3.semiminor(), 3.9); + + // Test ridx method + model3.ridx(0.7); + test_value(model3.ridx(), 0.7); + + // Test region method + GSkyDir dir3; + dir3.radec_deg(83.6331, +22.0145); + test_assert(model3.region()->contains(dir3), "Test region() method (inside)"); + dir3.radec_deg(83.6331, +61.5); + test_assert(!model3.region()->contains(dir3), "Test region() method (outside)"); + + // Test operator access + const char* strarray[] = {"RA", "DEC", "PA", "MinorRadius", "MajorRadius", "R_Index"}; + for (int i = 0; i < 6; ++i) { + std::string keyname(strarray[i]); + model3[keyname].value(2.1); + model3[keyname].error(1.9); + model3[keyname].gradient(0.8); + test_value(model3[keyname].value(), 2.1); + test_value(model3[keyname].error(), 1.9); + test_value(model3[keyname].gradient(), 0.8); + } + + // Exit test + return; +} /***********************************************************************//** * @brief Test GModelSpatialComposite class diff --git a/test/test_GModel.hpp b/test/test_GModel.hpp index a3e6a57..942bdae 100644 --- a/test/test_GModel.hpp +++ b/test/test_GModel.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * test_GModel.hpp - Test model module * * ----------------------------------------------------------------------- * - * copyright (C) 2012-2021 by Jean-Baptiste Cayrou * + * copyright (C) 2012-2022 by Jean-Baptiste Cayrou * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -58,12 +58,14 @@ public: void test_radial_disk(void); void test_radial_ring(void); void test_radial_gauss(void); + void test_radial_general_gauss(); void test_radial_shell(void); void test_profile_dmburkert(void); void test_profile_dmeinasto(void); void test_profile_dmzhao(void); void test_elliptical_disk(void); void test_elliptical_gauss(void); + void test_elliptical_general_gauss(void); void test_spatial_composite(void); void test_spatial_model(void); void test_const(void); @@ -147,12 +149,14 @@ private: std::string m_xml_model_radial_disk; std::string m_xml_model_radial_ring; std::string m_xml_model_radial_gauss; + std::string m_xml_model_radial_general_gauss; std::string m_xml_model_radial_shell; std::string m_xml_model_profile_dmburkert; std::string m_xml_model_profile_dmeinasto; std::string m_xml_model_profile_dmzhao; std::string m_xml_model_elliptical_disk; std::string m_xml_model_elliptical_gauss; + std::string m_xml_model_elliptical_general_gauss; std::string m_xml_model_spatial_composite; std::string m_xml_model_point_temp_lightcurve; std::string m_xml_model_point_temp_phasecurve;