diff --git a/ChangeLog b/ChangeLog index 115d362..e989304 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2016-02-28 Juergen Knoedlseder +2016-03-01 Juergen Knoedlseder * Version 1.1.0 released ========================= + Add GModelSpatialRadialProfile and GModelSpatialRadialProfileGauss classes (#1520) Make GFilename::is_fits thread safe (#1717). Revise GTime operators (#1714) Add GSkyMap::publish() method to publish map on VO Hub (#1710) diff --git a/NEWS b/NEWS index 5758b13..1b37e22 100644 --- a/NEWS +++ b/NEWS @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 1.1.0 Author(s) : Juergen Knoedlseder -28 February 2016 +1 March 2016 1. Introduction @@ -142,6 +142,14 @@ GModelSpatialDiffuseMap class now recompute only the Monte Carlo acceptance cone in case that the cone centre or radius have changed since the last computation. This speeds up the computations (#1673). +Add GModelSpatialRadialProfile class to handle arbitrary radial profiles. +The method deals with pre-computation and normalization of the profiles, +and implements the eval(), eval_grandients() and mc() method. The class +GModelSpatialRadialProfileGauss has been added to show how to implement +radial profile for the case of a Gaussian, and also to cross-check the +results with the GModelSpatialRadialGauss class that should return basically +the same results (#1520). + 9. Numerics module ------------------ diff --git a/include/GModelSpatialRadialProfile.hpp b/include/GModelSpatialRadialProfile.hpp new file mode 100644 index 0000000..4d5cfe2 --- /dev/null +++ b/include/GModelSpatialRadialProfile.hpp @@ -0,0 +1,138 @@ +/*************************************************************************** + * GModelSpatialRadialProfile.hpp - Radial profile source model class * + * ----------------------------------------------------------------------- * + * copyright (C) 2016 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * 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 GModelSpatialRadialProfile.hpp + * @brief Radial profile model class interface definition + * @author Juergen Knoedlseder + */ + +#ifndef GMODELSPATIALRADIALPROFILE_HPP +#define GMODELSPATIALRADIALPROFILE_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include "GModelSpatialRadial.hpp" +#include "GNodeArray.hpp" + +/* __ Forward declaration ________________________________________________ */ +class GXmlElement; +class GSkyDir; + + +/************************************************************************** + * @class GModelSpatialRadialProfile + * + * @brief Radial profile source model class + * + * This class implements the spatial component of the factorised source + * model for an arbitrary radial profile. + ***************************************************************************/ +class GModelSpatialRadialProfile : public GModelSpatialRadial { + +public: + // Constructors and destructors + GModelSpatialRadialProfile(void); + explicit GModelSpatialRadialProfile(const GXmlElement& xml); + GModelSpatialRadialProfile(const GModelSpatialRadialProfile& model); + virtual ~GModelSpatialRadialProfile(void); + + // Operators + virtual GModelSpatialRadialProfile& operator=(const GModelSpatialRadialProfile& model); + + // Pure virtual methods + virtual void clear(void) = 0; + virtual GModelSpatialRadialProfile* clone(void) const = 0; + virtual std::string classname(void) const = 0; + virtual std::string type(void) const = 0; + virtual double theta_max(void) const = 0; + virtual std::string print(const GChatter& chatter = NORMAL) const = 0; + + // Implemented pure virtual base class methods + virtual double eval(const double& theta, + const GEnergy& energy, + const GTime& time) const; + virtual double eval_gradients(const double& theta, + const GEnergy& energy, + const GTime& time) 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; + + // Implement other methods + int num_nodes(void) const; + void num_nodes(const int& number); + +protected: + // Protected methods + void init_members(void); + void copy_members(const GModelSpatialRadialProfile& model); + void free_members(void); + int cache_index(void) const; + virtual double profile_value(const double& theta) const = 0; + + // Protected members + int m_num_nodes; //!< Number of profile nodes + + // Pre-computed radial profile + struct profile { + std::vector pars; //!< Profile parameters + GNodeArray nodes; //!< Profile nodes + std::vector values; //!< Profile values + std::vector mc; //!< Profile for MC + double mc_max; //!< Maximum of profile for MC + }; + + // Pre-computation cache + mutable std::vector m_profile; // Pre-computation cache +}; + + +/***********************************************************************//** + * @brief Return number of nodes + * + * @return Number of nodes. + * + * Returns the number of nodes in the radial profile. + ***************************************************************************/ +inline +int GModelSpatialRadialProfile::num_nodes(void) const +{ + return (m_num_nodes); +} + + +/***********************************************************************//** + * @brief Set number of nodes + * + * @param[in] number Number of nodes. + * + * Sets the number of nodes in the radial profile. + ***************************************************************************/ +inline +void GModelSpatialRadialProfile::num_nodes(const int& number) +{ + m_num_nodes = number; + return; +} + +#endif /* GMODELSPATIALRADIALPROFILE_HPP */ diff --git a/include/GModelSpatialRadialProfileGauss.hpp b/include/GModelSpatialRadialProfileGauss.hpp new file mode 100644 index 0000000..db1efb8 --- /dev/null +++ b/include/GModelSpatialRadialProfileGauss.hpp @@ -0,0 +1,106 @@ +/*************************************************************************** + * GModelSpatialRadialProfileGauss.hpp - Gaussian radial profile class * + * ----------------------------------------------------------------------- * + * copyright (C) 2016 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * 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 GModelSpatialRadialProfileGauss.hpp + * @brief Radial Gaussian profile model class interface definition + * @author Juergen Knoedlseder + */ + +#ifndef GMODELSPATIALRADIALPROFILEGAUSS_HPP +#define GMODELSPATIALRADIALPROFILEGAUSS_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include "GModelSpatialRadialProfile.hpp" +#include "GModelPar.hpp" + +/* __ Forward declaration ________________________________________________ */ +class GXmlElement; + + +/************************************************************************** + * @class GModelSpatialRadialProfileGauss + * + * @brief Radial Gaussian profile source model class + * + * This class implements the spatial component of the factorised source + * model for a Gaussian radial profile. + ***************************************************************************/ +class GModelSpatialRadialProfileGauss : public GModelSpatialRadialProfile { + +public: + // Constructors and destructors + GModelSpatialRadialProfileGauss(void); + explicit GModelSpatialRadialProfileGauss(const GXmlElement& xml); + GModelSpatialRadialProfileGauss(const GModelSpatialRadialProfileGauss& model); + virtual ~GModelSpatialRadialProfileGauss(void); + + // Operators + virtual GModelSpatialRadialProfileGauss& operator=(const GModelSpatialRadialProfileGauss& model); + + // Implemented pure virtual base class methods + virtual void clear(void); + virtual GModelSpatialRadialProfileGauss* clone(void) const; + virtual std::string classname(void) const; + virtual std::string type(void) 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; + +protected: + // Protected methods + void init_members(void); + void copy_members(const GModelSpatialRadialProfileGauss& model); + void free_members(void); + virtual double profile_value(const double& theta) const; + + // Protected members + GModelPar m_sigma; //!< Gaussian width (deg) +}; + + +/***********************************************************************//** + * @brief Return class name + * + * @return String containing the class name ("GModelSpatialRadialProfileGauss"). + ***************************************************************************/ +inline +std::string GModelSpatialRadialProfileGauss::classname(void) const +{ + return ("GModelSpatialRadialProfileGauss"); +} + + +/***********************************************************************//** + * @brief Return model type + * + * @return "GaussianProfile". + * + * Returns the type of the radial profile model. + ***************************************************************************/ +inline +std::string GModelSpatialRadialProfileGauss::type(void) const +{ + return "GaussianProfile"; +} + +#endif /* GMODELSPATIALRADIALPROFILEGAUSS_HPP */ diff --git a/include/GammaLib.hpp b/include/GammaLib.hpp index a9756fe..a7a0fdb 100644 --- a/include/GammaLib.hpp +++ b/include/GammaLib.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GammaLib.hpp - Gamma-Ray Astronomy Library Header file * * ----------------------------------------------------------------------- * - * copyright (C) 2006-2015 by Juergen Knoedlseder * + * copyright (C) 2006-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -196,6 +196,8 @@ #include "GModelSpatialRadialDisk.hpp" #include "GModelSpatialRadialGauss.hpp" #include "GModelSpatialRadialShell.hpp" +#include "GModelSpatialRadialProfile.hpp" +#include "GModelSpatialRadialProfileGauss.hpp" #include "GModelSpatialElliptical.hpp" #include "GModelSpatialEllipticalDisk.hpp" #include "GModelSpatialEllipticalGauss.hpp" diff --git a/include/Makefile.am b/include/Makefile.am index bd6fade..cc2374c 100644 --- a/include/Makefile.am +++ b/include/Makefile.am @@ -133,6 +133,8 @@ pkginclude_HEADERS = GammaLib.hpp \ GModelSpatialRadialDisk.hpp \ GModelSpatialRadialGauss.hpp \ GModelSpatialRadialShell.hpp \ + GModelSpatialRadialProfile.hpp \ + GModelSpatialRadialProfileGauss.hpp \ GModelSpatialElliptical.hpp \ GModelSpatialEllipticalDisk.hpp \ GModelSpatialEllipticalGauss.hpp \ diff --git a/pyext/GModelSpatialRadialProfile.i b/pyext/GModelSpatialRadialProfile.i new file mode 100644 index 0000000..5b148fd --- /dev/null +++ b/pyext/GModelSpatialRadialProfile.i @@ -0,0 +1,72 @@ +/*************************************************************************** + * GModelSpatialRadialProfile.i - Radial profile source model class * + * ----------------------------------------------------------------------- * + * copyright (C) 2016 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * 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 GModelSpatialRadialProfile.i + * @brief Radial profile model class interface definition + * @author Juergen Knoedlseder + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GModelSpatialRadialProfile.hpp" +%} + + +/************************************************************************** + * @class GModelSpatialRadialProfile + * + * @brief Radial profile source model class + * + * This class implements the spatial component of the factorised source + * model for an arbitrary radial profile. + ***************************************************************************/ +class GModelSpatialRadialProfile : public GModelSpatialRadial { + +public: + // Constructors and destructors + GModelSpatialRadialProfile(void); + explicit GModelSpatialRadialProfile(const GXmlElement& xml); + GModelSpatialRadialProfile(const GModelSpatialRadialProfile& model); + virtual ~GModelSpatialRadialProfile(void); + + // Pure virtual methods + virtual void clear(void) = 0; + virtual GModelSpatialRadialProfile* clone(void) const = 0; + virtual std::string classname(void) const = 0; + virtual std::string type(void) const = 0; + virtual double theta_max(void) const = 0; + + // Implemented pure virtual base class methods + virtual double eval(const double& theta, + const GEnergy& energy, + const GTime& time) const; + virtual double eval_gradients(const double& theta, + const GEnergy& energy, + const GTime& time) 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; + + // Implement other methods + int num_nodes(void) const; + void num_nodes(const int& number); +}; diff --git a/pyext/GModelSpatialRadialProfileGauss.i b/pyext/GModelSpatialRadialProfileGauss.i new file mode 100644 index 0000000..f74d61c --- /dev/null +++ b/pyext/GModelSpatialRadialProfileGauss.i @@ -0,0 +1,76 @@ +/*************************************************************************** + * GModelSpatialRadialProfileGauss.i - Gaussian radial profile class * + * ----------------------------------------------------------------------- * + * copyright (C) 2016 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * 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 GModelSpatialRadialProfileGauss.hpp + * @brief Radial Gaussian profile model class interface definition + * @author Juergen Knoedlseder + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GModelSpatialRadialProfileGauss.hpp" +%} + + +/************************************************************************** + * @class GModelSpatialRadialProfileGauss + * + * @brief Radial Gaussian profile source model class + * + * This class implements the spatial component of the factorised source + * model for a Gaussian radial profile. + ***************************************************************************/ +class GModelSpatialRadialProfileGauss : public GModelSpatialRadialProfile { + +public: + // Constructors and destructors + GModelSpatialRadialProfileGauss(void); + explicit GModelSpatialRadialProfileGauss(const GXmlElement& xml); + GModelSpatialRadialProfileGauss(const GModelSpatialRadialProfileGauss& model); + virtual ~GModelSpatialRadialProfileGauss(void); + + // Implemented pure virtual base class methods + virtual void clear(void); + virtual GModelSpatialRadialProfileGauss* clone(void) const; + virtual std::string classname(void) const; + virtual std::string type(void) const; + virtual double theta_max(void) const; + virtual void read(const GXmlElement& xml); + virtual void write(GXmlElement& xml) const; +}; + + +/***********************************************************************//** + * @brief GModelSpatialRadialGauss class extension + * + * The eval(GSkyDir&) and eval_gradients(GSkyDir&) 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 GModelSpatialRadialProfileGauss { + double eval(const GPhoton& photon) const { + return self->GModelSpatialRadial::eval(photon); + } + double eval_gradients(const GPhoton& photon) const { + return self->GModelSpatialRadial::eval_gradients(photon); + } +}; diff --git a/pyext/gammalib/model.i b/pyext/gammalib/model.i index d17a4a2..79003a2 100644 --- a/pyext/gammalib/model.i +++ b/pyext/gammalib/model.i @@ -134,6 +134,8 @@ %include "GModelSpatialRadialDisk.i" %include "GModelSpatialRadialGauss.i" %include "GModelSpatialRadialShell.i" +%include "GModelSpatialRadialProfile.i" +%include "GModelSpatialRadialProfileGauss.i" %include "GModelSpatialElliptical.i" %include "GModelSpatialEllipticalDisk.i" %include "GModelSpatialEllipticalGauss.i" diff --git a/src/model/GModelSpatialRadialProfile.cpp b/src/model/GModelSpatialRadialProfile.cpp new file mode 100644 index 0000000..9e898b9 --- /dev/null +++ b/src/model/GModelSpatialRadialProfile.cpp @@ -0,0 +1,454 @@ +/*************************************************************************** + * GModelSpatialRadialProfile.cpp - Radial profile source model class * + * ----------------------------------------------------------------------- * + * copyright (C) 2016 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * 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 GModelSpatialRadialProfile.cpp + * @brief Radial profile model class implementation + * @author Juergen Knoedlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include "GException.hpp" +#include "GTools.hpp" +#include "GMath.hpp" +#include "GSkyDir.hpp" +#include "GXmlElement.hpp" +#include "GModelSpatialRadialProfile.hpp" + +/* __ Constants __________________________________________________________ */ + +/* __ Globals ____________________________________________________________ */ + +/* __ Method name definitions ____________________________________________ */ + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + +/* __ Debug definitions __________________________________________________ */ +//#define G_DEBUG_PRECOMPUTATION + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + * + * Constructs empty radial profile + ***************************************************************************/ +GModelSpatialRadialProfile::GModelSpatialRadialProfile(void) : + GModelSpatialRadial() +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief XML constructor + * + * @param[in] xml XML element. + * + * Constructs radial profile model by extracting information from an XML + * elements. + ***************************************************************************/ +GModelSpatialRadialProfile::GModelSpatialRadialProfile(const GXmlElement& xml) : + GModelSpatialRadial() +{ + // Initialise members + init_members(); + + // Read information from XML element + read(xml); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] model Radial profile model. + * + * Copies radial profile model from another radial profile model. + ***************************************************************************/ +GModelSpatialRadialProfile::GModelSpatialRadialProfile(const GModelSpatialRadialProfile& model) : + GModelSpatialRadial(model) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(model); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + * + * Destructs radial profile model. + ***************************************************************************/ +GModelSpatialRadialProfile::~GModelSpatialRadialProfile(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] model Radial profile model. + * @return Radial profile model. + * + * Assigns radial profile model. + ***************************************************************************/ +GModelSpatialRadialProfile& GModelSpatialRadialProfile::operator=(const GModelSpatialRadialProfile& 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 Evaluate function (in units of sr^-1) + * + * @param[in] theta Angular distance from model centre (radians). + * @param[in] energy Photon energy. + * @param[in] time Photon arrival time. + * @return Model value. + * + * Evaluate the radial profile model for a given angular distance @p theta + * the model centre, a given photon @p energy and a given @p time. The + * method evaluates the model by linear interpolation. + ***************************************************************************/ +double GModelSpatialRadialProfile::eval(const double& theta, + const GEnergy& energy, + const GTime& time) const +{ + // Get pre-computation cache index + int icache = cache_index(); + + // Get interpolation value + double value = m_profile[icache].nodes.interpolate(theta, m_profile[icache].values); + + // Make sure that value is not-negative + if (value < 0.0) { + value = 0.0; + } + + // Compile option: Check for NaN/Inf + #if defined(G_NAN_CHECK) + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { + std::cout << "*** ERROR: GModelSpatialRadialProfile::eval"; + std::cout << "(theta=" << theta << "): NaN/Inf encountered"; + std::cout << " (value=" << value; + std::cout << ")" << std::endl; + } + #endif + + // Return value + return value; +} + + +/***********************************************************************//** + * @brief Evaluate function and gradients (in units of sr^-1) + * + * @param[in] theta Angular distance from model centre (radians). + * @param[in] energy Photon energy. + * @param[in] time Photon arrival time. + * @return Model value. + * + * Evaluates the function value. No gradient computation is implemented as + * radial models will be convolved with the instrument response and thus + * require the numerical computation of the derivatives. + * + * See the eval() method for more information. + ***************************************************************************/ +double GModelSpatialRadialProfile::eval_gradients(const double& theta, + const GEnergy& energy, + const GTime& time) const +{ + // Return value + return (eval(theta, energy, time)); +} + + +/***********************************************************************//** + * @brief Return 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 radial profile using a rejection + * method. + ***************************************************************************/ +GSkyDir GModelSpatialRadialProfile::mc(const GEnergy& energy, + const GTime& time, + GRan& ran) const +{ + // Get pre-computation cache index + int icache = cache_index(); + + // Simulate theta angle (degrees) using a rejection method + double theta = 0.0; + while (true) { + theta = ran.uniform() * theta_max(); + double mc = m_profile[icache].nodes.interpolate(theta, m_profile[icache].mc); + if ((ran.uniform() * m_profile[icache].mc_max) <= mc) { + break; + } + } + theta *= gammalib::rad2deg; + + // Simulate uniform azmiuth angle (degrees) + double phi = 360.0 * ran.uniform(); + + // Rotate sky direction by offset + GSkyDir sky_dir = dir(); + sky_dir.rotate_deg(phi, theta); + + // 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 (degrees) + * @return True if the model contains the sky direction. + * + * Signals whether a sky direction is contained in the radial disk model. + ***************************************************************************/ +bool GModelSpatialRadialProfile::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); +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void GModelSpatialRadialProfile::init_members(void) +{ + // Initialise members + m_num_nodes = 100; + + // Initialise pre-computation cache + m_profile.clear(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] model Radial disk model. + ***************************************************************************/ +void GModelSpatialRadialProfile::copy_members(const GModelSpatialRadialProfile& model) +{ + // Copy members + m_num_nodes = model.m_num_nodes; + m_profile = model.m_profile; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void GModelSpatialRadialProfile::free_members(void) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief Return index to pre-computation cache + * + * @return Index to pre-computation cache + * + * Returns the index to the pre-computation cache. If no pre-computation + * cache was found the method will create one and return the index to that + * cache. + ***************************************************************************/ +int GModelSpatialRadialProfile::cache_index(void) const +{ + // Initialise index + int index = -1; + + // Search for index + for (int i = 0; i < m_profile.size(); ++i) { + + // Initialise found flag + bool found = true; + + // Loop over all spatial parameters and check if they are equal + // to the parameters in the pre-computation cache. Break and set + // the found flag to false on non-equality. + for (int k = 0; k < m_pars.size(); ++k) { + if (m_pars[k]->value() != m_profile[i].pars[k]) { + found = false; + break; + } + } + + // If profile was found then set index + if (found) { + index = i; + break; + } + + } // endfor: looped over index + + // If index was not found then allocate a new pre-computation cache + if (index == -1) { + + // Initialise profile + profile prf; + prf.nodes.clear(); + prf.values.clear(); + prf.mc.clear(); + prf.mc_max = 0.0; + + // Set profile parameters + for (int k = 0; k < m_pars.size(); ++k) { + prf.pars.push_back(m_pars[k]->value()); + } + + // Pre-compute radial values, MC values, and mc_max. Compute + // also the normalization. + double rmax = theta_max(); + double dr = rmax / m_num_nodes; + double r = 0.0; + double norm = 0.0; + for (int j = 0; j < m_num_nodes; ++j) { + double value = profile_value(r); + double mc = value * std::sin(r) * dr; + norm += mc; + if (mc > prf.mc_max) { + prf.mc_max = mc; + } + prf.nodes.append(r); + prf.values.push_back(value); + prf.mc.push_back(mc); + r += dr; + } + norm *= gammalib::twopi; + + // Normalize radial profile + if (norm > 0.0) { + double inv_norm = 1.0 / norm; + for (int j = 0; j < m_num_nodes; ++j) { + prf.values[j] *= inv_norm; + } + } + + // Push profile in pre-computation cache + m_profile.push_back(prf); + + // Get index to last profile + index = m_profile.size() - 1; + + // Log pre-computation + #if defined(G_DEBUG_PRECOMPUTATION) + std::cout << "GModelSpatialRadialProfile::cache_index"; + std::cout << std::endl; + std::cout << " dr=" << dr << std::endl; + for (int k = 0; k < m_pars.size(); ++k) { + std::cout << " par[" << k << "]=" << m_pars[k]->value(); + std::cout << " (" << m_pars[k]->name(); + std::cout << ")" << std::endl; + } + std::cout << " norm=" << norm << std::endl; + #endif + + } // endif: allocated new pre-computation cache + + // Return index + return (index); +} diff --git a/src/model/GModelSpatialRadialProfileGauss.cpp b/src/model/GModelSpatialRadialProfileGauss.cpp new file mode 100644 index 0000000..ecaa362 --- /dev/null +++ b/src/model/GModelSpatialRadialProfileGauss.cpp @@ -0,0 +1,401 @@ +/*************************************************************************** + * GModelSpatialRadialProfileGauss.cpp - Gaussian radial profile class * + * ----------------------------------------------------------------------- * + * copyright (C) 2016 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * 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 GModelSpatialRadialProfileGauss.cpp + * @brief Radial Gaussian profile model class implementation + * @author Juergen Knoedlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include "GException.hpp" +#include "GTools.hpp" +#include "GMath.hpp" +#include "GXmlElement.hpp" +#include "GModelSpatialRadialProfileGauss.hpp" +#include "GModelSpatialRegistry.hpp" + +/* __ Constants __________________________________________________________ */ + +/* __ Globals ____________________________________________________________ */ +const GModelSpatialRadialProfileGauss g_radial_disk_seed; +const GModelSpatialRegistry g_radial_disk_registry(&g_radial_disk_seed); + +/* __ Method name definitions ____________________________________________ */ +#define G_READ "GModelSpatialRadialProfileGauss::read(GXmlElement&)" +#define G_WRITE "GModelSpatialRadialProfileGauss::write(GXmlElement&)" + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + +/* __ Debug definitions __________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + * + * Constructs empty radial Gaussian profile + ***************************************************************************/ +GModelSpatialRadialProfileGauss::GModelSpatialRadialProfileGauss(void) : + GModelSpatialRadialProfile() +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief XML constructor + * + * @param[in] xml XML element. + * + * Constructs radial Gaussian profile model by extracting information from + * an XML element. See the read() method for more information about the + * expected structure of the XML element. + ***************************************************************************/ +GModelSpatialRadialProfileGauss::GModelSpatialRadialProfileGauss(const GXmlElement& xml) : + GModelSpatialRadialProfile() +{ + // Initialise members + init_members(); + + // Read information from XML element + read(xml); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] model Radial Gaussian profile model. + * + * Copies radial Gaussian profile model from another radial profile model. + ***************************************************************************/ +GModelSpatialRadialProfileGauss::GModelSpatialRadialProfileGauss(const GModelSpatialRadialProfileGauss& model) : + GModelSpatialRadialProfile(model) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(model); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + * + * Destructs radial Gaussian profile model. + ***************************************************************************/ +GModelSpatialRadialProfileGauss::~GModelSpatialRadialProfileGauss(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] model Radial Gaussian profile model. + * @return Radial Gaussian profile model. + * + * Assigns radial Gaussian profile model. + ***************************************************************************/ +GModelSpatialRadialProfileGauss& GModelSpatialRadialProfileGauss::operator=(const GModelSpatialRadialProfileGauss& model) +{ + // Execute only if object is not identical + if (this != &model) { + + // Copy base class members + this->GModelSpatialRadialProfile::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 Gaussian profile model + * + * Clears radial Gaussian profile model. + ***************************************************************************/ +void GModelSpatialRadialProfileGauss::clear(void) +{ + // Free class members + free_members(); + this->GModelSpatialRadialProfile::free_members(); + this->GModelSpatialRadial::free_members(); + this->GModelSpatial::free_members(); + + // Initialise members + this->GModelSpatial::init_members(); + this->GModelSpatialRadial::init_members(); + this->GModelSpatialRadialProfile::init_members(); + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Clone radial Gaussian profile model + * + * @return Pointer to deep copy of radial Gaussian profile model. + * + * Returns a deep copy of the radial Gaussian profile model. + ***************************************************************************/ +GModelSpatialRadialProfileGauss* GModelSpatialRadialProfileGauss::clone(void) const +{ + // Clone radial disk model + return new GModelSpatialRadialProfileGauss(*this); +} + + +/***********************************************************************//** + * @brief Return maximum model radius (in radians) + * + * @return Maximum model radius (in radians). + ***************************************************************************/ +double GModelSpatialRadialProfileGauss::theta_max(void) const +{ + // Return value + return (m_sigma.value() * gammalib::deg2rad * 5.0); +} + + +/***********************************************************************//** + * @brief Read model from XML element + * + * @param[in] xml XML element. + * + * Reads the Gaussian radial profile model information from an XML element. + * The XML element shall have either the format + * + * + * + * + * + * + * + * or + * + * + * + * + * + * + ***************************************************************************/ +void GModelSpatialRadialProfileGauss::read(const GXmlElement& xml) +{ + // Read Gaussian location + GModelSpatialRadial::read(xml); + + // Read Sigma parameter + const GXmlElement* par = gammalib::xml_get_par(G_READ, xml, "Sigma"); + m_sigma.read(*par); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Write model into XML element + * + * @param[in] xml XML element into which model information is written. + * + * Writes the Gaussian radial profile model information into an XML element. + * The XML element will have the format + * + * + * + * + * + * + ***************************************************************************/ +void GModelSpatialRadialProfileGauss::write(GXmlElement& xml) const +{ + // Write Gaussian location + GModelSpatialRadial::write(xml); + + // Write Sigma parameter + GXmlElement* par = gammalib::xml_need_par(G_WRITE, xml, "Sigma"); + m_sigma.write(*par); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Print information + * + * @param[in] chatter Chattiness (defaults to NORMAL). + * @return String containing model information. + ***************************************************************************/ +std::string GModelSpatialRadialProfileGauss::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("=== GModelSpatialRadialProfileGauss ==="); + + // 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 GModelSpatialRadialProfileGauss::init_members(void) +{ + // Initialise Gaussian sigma + m_sigma.clear(); + m_sigma.name("Sigma"); + m_sigma.unit("deg"); + m_sigma.value(2.778e-4); // 1 arcsec + m_sigma.min(2.778e-4); // 1 arcsec + m_sigma.free(); + m_sigma.scale(1.0); + m_sigma.gradient(0.0); + m_sigma.has_grad(false); // Radial components never have gradients + + // Set parameter pointer(s) + m_pars.push_back(&m_sigma); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] model Radial Gaussian model. + * + * Copies class members from another radial profile model. + ***************************************************************************/ +void GModelSpatialRadialProfileGauss::copy_members(const GModelSpatialRadialProfileGauss& model) +{ + // Copy members. 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. + m_sigma = model.m_sigma; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void GModelSpatialRadialProfileGauss::free_members(void) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief Radial profile + * + * @param[in] theta Angular distance from Gaussian centre (radians). + * @return Profile value. + ***************************************************************************/ +double GModelSpatialRadialProfileGauss::profile_value(const double& theta) const +{ + // Compute value + double sigma_rad = m_sigma.value() * gammalib::deg2rad; + double sigma2 = sigma_rad * sigma_rad; + double theta2 = theta * theta; + double value = std::exp(-0.5 * theta2 / sigma2) / + (gammalib::twopi * sigma2); + + // Return value + return value; +} diff --git a/src/model/Makefile.am b/src/model/Makefile.am index 8e9a622..1491af0 100644 --- a/src/model/Makefile.am +++ b/src/model/Makefile.am @@ -19,6 +19,8 @@ sources = GModelPar.cpp \ GModelSpatialRadialDisk.cpp \ GModelSpatialRadialGauss.cpp \ GModelSpatialRadialShell.cpp \ + GModelSpatialRadialProfile.cpp \ + GModelSpatialRadialProfileGauss.cpp \ GModelSpatialElliptical.cpp \ GModelSpatialEllipticalDisk.cpp \ GModelSpatialEllipticalGauss.cpp \