diff --git a/include/GModelSpatialRadialRing.hpp b/include/GModelSpatialRadialRing.hpp new file mode 100644 index 0000000..946afa9 --- /dev/null +++ b/include/GModelSpatialRadialRing.hpp @@ -0,0 +1,210 @@ +/*************************************************************************** + * GModelSpatialRadialRing.hpp - Radial ring source model class * + * ----------------------------------------------------------------------- * + * copyright (C) 2011-2016 by Pierrick Martin * + * ----------------------------------------------------------------------- * + * * + * 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 GModelSpatialRadialRing.hpp + * @brief Radial ring model class interface definition + * @author Pierrick Martin + */ + +#ifndef GMODELSPATIALRADIALRING_HPP +#define GMODELSPATIALRADIALRING_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include "GModelSpatialRadial.hpp" +#include "GModelPar.hpp" +#include "GSkyDir.hpp" +#include "GSkyRegionCircle.hpp" +#include "GXmlElement.hpp" + + +/************************************************************************** + * @class GModelSpatialRadialRing + * + * @brief Ring source model class + * + * This class implements the spatial component of the factorised source + * model for a ring source, i.e. constant surface brightness between an + * inner radius and an outer radius, and no emission elsewhere. + ***************************************************************************/ +class GModelSpatialRadialRing : public GModelSpatialRadial { + +public: + // Constructors and destructors + GModelSpatialRadialRing(void); + GModelSpatialRadialRing(const bool& dummy, const std::string& type); + GModelSpatialRadialRing(const GSkyDir& dir, const double& radius, const double& width); + explicit GModelSpatialRadialRing(const GXmlElement& xml); + GModelSpatialRadialRing(const GModelSpatialRadialRing& model); + virtual ~GModelSpatialRadialRing(void); + + // Operators + virtual GModelSpatialRadialRing& operator=(const GModelSpatialRadialRing& model); + + // Implemented pure virtual base class methods + virtual void clear(void); + virtual GModelSpatialRadialRing* clone(void) const; + virtual std::string classname(void) const; + virtual std::string type(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 double theta_min(void) const; + virtual GSkyRegion* region(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 width(void) const; + void width(const double& width); + +protected: + // Protected methods + void init_members(void); + void copy_members(const GModelSpatialRadialRing& model); + void free_members(void); + void update(void) const; + void set_region(void) const; + + // Protected members + std::string m_type; //!< Model type + GModelPar m_radius; //!< Ring inner radius (degrees) + GModelPar m_width; //!< Ring width (degrees) + mutable GSkyRegionCircle m_region; //!< Bounding circle + + // Cached members used for pre-computations + mutable double m_last_radius; //!< Last ring radius + mutable double m_last_width; //!< Last ring width + mutable double m_inner_radius_rad; //!< Inner radius in radians + mutable double m_outer_radius_rad; //!< Outer radius in radians + mutable double m_width_rad; //!< Ring width in radians + mutable double m_norm; //!< Normalization +}; + + +/***********************************************************************//** + * @brief Return class name + * + * @return String containing the class name ("GModelSpatialRadialRing"). + ***************************************************************************/ +inline +std::string GModelSpatialRadialRing::classname(void) const +{ + return ("GModelSpatialRadialRing"); +} + + +/***********************************************************************//** + * @brief Return model type + * + * @return Model type. + * + * Returns the type of the radial ring model. + ***************************************************************************/ +inline +std::string GModelSpatialRadialRing::type(void) const +{ + return (m_type); +} + + +/***********************************************************************//** + * @brief Return ring outer radius + * + * @return Ring outer radius (degrees). + * + * Returns the outer radius of the ring in degrees. + ***************************************************************************/ +inline +double GModelSpatialRadialRing::radius(void) const +{ + return (m_radius.value()); +} + + +/***********************************************************************//** + * @brief Set ring outer radius + * + * @param[in] Ring outer radius (degrees). + * + * Sets the outer radius of the ring in degrees. + ***************************************************************************/ +inline +void GModelSpatialRadialRing::radius(const double& radius) +{ + m_radius.value(radius); + return; +} + + +/***********************************************************************//** + * @brief Return ring width + * + * @return Ring width (degrees). + * + * Returns the width of the ring in degrees. + ***************************************************************************/ +inline +double GModelSpatialRadialRing::width(void) const +{ + return (m_width.value()); +} + + +/***********************************************************************//** + * @brief Set ring width + * + * @param[in] Ring width (degrees). + * + * Sets the width of the ring in degrees. + ***************************************************************************/ +inline +void GModelSpatialRadialRing::width(const double& width) +{ + m_width.value(width); + return; +} + +/***********************************************************************//** + * @brief Return boundary sky region + * + * @return Boundary sky region. + * + * Returns a sky region that fully encloses the spatial model component. + ***************************************************************************/ +inline +GSkyRegion* GModelSpatialRadialRing::region(void) const +{ + set_region(); + return (&m_region); +} + +#endif /* GMODELSPATIALRADIALRING_HPP */ diff --git a/include/GammaLib.hpp b/include/GammaLib.hpp index 11c0343..5a89acd 100644 --- a/include/GammaLib.hpp +++ b/include/GammaLib.hpp @@ -207,6 +207,7 @@ #include "GModelSpatialRadialDisk.hpp" #include "GModelSpatialRadialGauss.hpp" #include "GModelSpatialRadialShell.hpp" +#include "GModelSpatialRadialRing.hpp" #include "GModelSpatialRadialProfile.hpp" #include "GModelSpatialRadialProfileGauss.hpp" #include "GModelSpatialRadialProfileDMBurkert.hpp" diff --git a/include/Makefile.am b/include/Makefile.am index 6a848bc..cd05720 100644 --- a/include/Makefile.am +++ b/include/Makefile.am @@ -165,6 +165,7 @@ pkginclude_HEADERS = GammaLib.hpp \ GModelSpatialRadialDisk.hpp \ GModelSpatialRadialGauss.hpp \ GModelSpatialRadialShell.hpp \ + GModelSpatialRadialRing.hpp \ GModelSpatialRadialProfile.hpp \ GModelSpatialRadialProfileGauss.hpp \ GModelSpatialRadialProfileDMBurkert.hpp \ diff --git a/pyext/GModelSpatialRadialRing.i b/pyext/GModelSpatialRadialRing.i new file mode 100644 index 0000000..ebd4c9a --- /dev/null +++ b/pyext/GModelSpatialRadialRing.i @@ -0,0 +1,106 @@ +/*************************************************************************** + * GModelSpatialRadialRing.i - Radial disk source model class * + * ----------------------------------------------------------------------- * + * copyright (C) 2011-2020 by Pierrick Martin * + * ----------------------------------------------------------------------- * + * * + * 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 GModelSpatialRadialRing.i + * @brief Radial ring model class Python interface definition + * @author Pierrick Martin + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GModelSpatialRadialRing.hpp" +%} + +/************************************************************************** + * @class GModelSpatialRadialRing + * + * @brief Ring source model class + * + * This class implements the spatial component of the factorised source + * model for a ring source, i.e. constant surface brightness between an + * inner radius and an outer radius, and no emission elsewhere. + ***************************************************************************/ +class GModelSpatialRadialRing : public GModelSpatialRadial { + +public: + // Constructors and destructors + GModelSpatialRadialRing(void); + GModelSpatialRadialRing(const GSkyDir& dir, const double& radius, const double& width); + explicit GModelSpatialRadialRing(const GXmlElement& xml); + GModelSpatialRadialRing(const GModelSpatialRadialRing& model); + virtual ~GModelSpatialRadialRing(void); + + // Implemented pure virtual methods + virtual void clear(void); + virtual GModelSpatialRadialRing* clone(void) const; + virtual std::string classname(void) const; + virtual std::string type(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 double theta_min(void) const; + virtual GSkyRegion* region(void) const; + virtual void read(const GXmlElement& xml); + virtual void write(GXmlElement& xml) const; + + // Other methods + double radius(void) const; + void radius(const double& radius); + double width(void) const; + void width(const double& width); +}; + + +/***********************************************************************//** + * @brief GModelSpatialRadialRing 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 GModelSpatialRadialRing { + GModelSpatialRadialRing copy() { + return (*self); + } + double eval(const GPhoton& photon) const { + return self->GModelSpatialRadial::eval(photon); + } + double eval(const GPhoton& photon, const bool& gradients) const { + return self->GModelSpatialRadial::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 10092d5..e096831 100644 --- a/pyext/gammalib/model.i +++ b/pyext/gammalib/model.i @@ -170,6 +170,7 @@ %include "GModelSpatialPointSource.i" %include "GModelSpatialRadial.i" %include "GModelSpatialRadialDisk.i" +%include "GModelSpatialRadialRing.i" %include "GModelSpatialRadialGauss.i" %include "GModelSpatialRadialShell.i" %include "GModelSpatialRadialProfile.i" diff --git a/src/model/GModelSpatialRadialRing.cpp b/src/model/GModelSpatialRadialRing.cpp new file mode 100644 index 0000000..78b56e5 --- /dev/null +++ b/src/model/GModelSpatialRadialRing.cpp @@ -0,0 +1,697 @@ +/*************************************************************************** + * GModelSpatialRadialRing.cpp - Radial ring source model class * + * ----------------------------------------------------------------------- * + * copyright (C) 2011-2018 by Pierrick Martin * + * ----------------------------------------------------------------------- * + * * + * 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 GModelSpatialRadialRing.cpp + * @brief Radial ring model class implementation + * @author Pierrick Martin + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include "GException.hpp" +#include "GTools.hpp" +#include "GMath.hpp" +#include "GModelSpatialRadialRing.hpp" +#include "GModelSpatialRegistry.hpp" + +/* __ Constants __________________________________________________________ */ + +/* __ Globals ____________________________________________________________ */ +const GModelSpatialRadialRing g_radial_ring_seed; +const GModelSpatialRegistry g_radial_ring_registry(&g_radial_ring_seed); +#if defined(G_LEGACY_XML_FORMAT) +const GModelSpatialRadialRing g_radial_ring_legacy_seed(true, "RingFunction"); +const GModelSpatialRegistry g_radial_ring_legacy_registry(&g_radial_ring_legacy_seed); +#endif + +/* __ Method name definitions ____________________________________________ */ +#define G_READ "GModelSpatialRadialRing::read(GXmlElement&)" +#define G_WRITE "GModelSpatialRadialRing::write(GXmlElement&)" + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + +/* __ Debug definitions __________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + * + * Constructs empty radial ring model. + ***************************************************************************/ +GModelSpatialRadialRing::GModelSpatialRadialRing(void) : GModelSpatialRadial() +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Model type constructor + * + * @param[in] dummy Dummy flag. + * @param[in] type Model type. + * + * Constructs empty radial ring model by specifying a model @p type. + ***************************************************************************/ +GModelSpatialRadialRing::GModelSpatialRadialRing(const bool& dummy, + const std::string& type) : + GModelSpatialRadial() +{ + // Initialise members + init_members(); + + // Set model type + m_type = type; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Ring constructor + * + * @param[in] dir Sky position of ring centre. + * @param[in] radius Ring outer radius (degrees). + * @param[in] width Ring width (degrees). + * + * Constructs radial ring model from the sky position of the ring centre + * (@p dir) and the ring outer radius @p radius and width + * @p width in degrees. + ***************************************************************************/ +GModelSpatialRadialRing::GModelSpatialRadialRing(const GSkyDir& dir, + const double& radius, + const double& width) : + GModelSpatialRadial() +{ + // Initialise members + init_members(); + + // Assign center location + this->dir(dir); + + // Assign radius and width after check + if (width <= radius) { + this->radius(radius); + this->width(width); + } else { + this->radius(radius); + this->width(radius); + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief XML constructor + * + * @param[in] xml XML element. + * + * Creates instance of radial ring model by extracting information from an + * XML element. See the read() method for more information about the expected + * structure of the XML element. + ***************************************************************************/ +GModelSpatialRadialRing::GModelSpatialRadialRing(const GXmlElement& xml) : + GModelSpatialRadial() +{ + // Initialise members + init_members(); + + // Read information from XML element + read(xml); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] model Radial ring model. + ***************************************************************************/ +GModelSpatialRadialRing::GModelSpatialRadialRing(const GModelSpatialRadialRing& model) : + GModelSpatialRadial(model) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(model); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +GModelSpatialRadialRing::~GModelSpatialRadialRing(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] model Radial ring model. + * @return Radial ring model. + ***************************************************************************/ +GModelSpatialRadialRing& GModelSpatialRadialRing::operator=(const GModelSpatialRadialRing& 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 ring model + ***************************************************************************/ +void GModelSpatialRadialRing::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 ring model + * + * @return Pointer to deep copy of radial ring model. + ***************************************************************************/ +GModelSpatialRadialRing* GModelSpatialRadialRing::clone(void) const +{ + // Clone radial ring model + return new GModelSpatialRadialRing(*this); +} + + +/***********************************************************************//** + * @brief Evaluate function (in units of sr^-1) + * + * @param[in] theta Angular distance from ring centre (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 ring source model using + * + * \f[ + * S_{\rm p}(\vec{p} | E, t) = \left \{ + * \begin{array}{l l} + * \displaystyle + * {\tt m\_norm} + * & \mbox{if $\theta \le $ outer radius and $\theta \ge $ inner radius} \\ + * \\ + * \displaystyle + * 0 & \mbox{if $\theta > $ outer radius or $\theta < $ inner radius} + * \end{array} + * \right . + * \f] + * + * where + * - \f$\theta\f$ is the angular separation between ring centre and the + * sky direction of interest, and + * - \f${\tt m\_norm} = \frac{1}{2 \pi (\cos r_i - \cos r_o)} \f$ is a normalization + * constant (see the update() method for details). + * + * The method will not compute analytical parameter gradients, even if the + * @p gradients argument is set to true. Radial ring parameter gradients + * need to be computed numerically. + ***************************************************************************/ +double GModelSpatialRadialRing::eval(const double& theta, + const GEnergy& energy, + const GTime& time, + const bool& gradients) const +{ + // Update precomputation cache + update(); + + // Set value + double value = 0.0; + value = (theta <= m_outer_radius_rad) ? 1.0 : 0.0; + value *= (theta >= m_inner_radius_rad) ? m_norm : 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: GModelSpatialRadialRing::eval"; + std::cout << "(theta=" << theta << "): NaN/Inf encountered"; + std::cout << " (value=" << value; + std::cout << ", m_inner_radius_rad=" << m_inner_radius_rad; + std::cout << ", m_outer_radius_rad=" << m_outer_radius_rad; + std::cout << ", m_norm=" << m_norm; + std::cout << ")" << std::endl; + } + #endif + + // Return value + return value; +} + + +/***********************************************************************//** + * @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 2D ring distribution as function + * of the photon @p energy and arrival @p time. + ***************************************************************************/ +GSkyDir GModelSpatialRadialRing::mc(const GEnergy& energy, + const GTime& time, + GRan& ran) const +{ + // Update precomputation cache + update(); + + // Simulate offset from photon arrival direction + double cosirad = std::cos(m_inner_radius_rad); + double cosorad = std::cos(m_outer_radius_rad); + double theta = std::acos(cosirad - ran.uniform() * (cosirad - cosorad)) * + gammalib::rad2deg; + 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 whether 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 ring model. + ***************************************************************************/ +bool GModelSpatialRadialRing::contains(const GSkyDir& dir, + const double& margin) const +{ + // Update precomputation cache + update(); + + // Compute distance to centre (radians) + double distance = dir.dist(this->dir()); + + // Set flag + bool flag; + flag = (distance <= m_outer_radius_rad + margin*gammalib::deg2rad) + *(distance >= m_inner_radius_rad - margin*gammalib::deg2rad); + + // Return flag + return flag; +} + + +/***********************************************************************//** + * @brief Return maximum model radius (in radians) + * + * @return Maximum model radius (in radians). + ***************************************************************************/ +double GModelSpatialRadialRing::theta_max(void) const +{ + // Return value + return m_outer_radius_rad; +} + + +/***********************************************************************//** + * @brief Return minimum model radius (in radians) + * + * @return Minimum model radius (in radians). + ***************************************************************************/ +double GModelSpatialRadialRing::theta_min(void) const +{ + // Return value + return m_inner_radius_rad; +} + +/***********************************************************************//** + * @brief Read model from XML element + * + * @param[in] xml XML element. + * + * @exception GException::model_invalid_parnum + * Invalid number of model parameters found in XML element. + * @exception GException::model_invalid_parnames + * Invalid model parameter names found in XML element. + * + * Reads the radial ring 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 radius + * and radius minimum should be >0. + ***************************************************************************/ +void GModelSpatialRadialRing::read(const GXmlElement& xml) +{ + // Determine number of parameter nodes in XML element + int npars = xml.elements("parameter"); + + // Verify that XML element has exactly 4 parameters + if (xml.elements() != 4 || npars != 4) { + throw GException::model_invalid_parnum(G_READ, xml, + "Ring model requires exactly 4 parameters."); + } + + // Read ring center location + GModelSpatialRadial::read(xml); + + // Get parameters + const GXmlElement* radius = gammalib::xml_get_par(G_READ, xml, m_radius.name()); + const GXmlElement* width = gammalib::xml_get_par(G_READ, xml, m_width.name()); + + // Read parameters + m_radius.read(*radius); + m_width.read(*width); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Write model into XML element + * + * @param[in] xml XML element into which model information is written. + * + * @exception GException::model_invalid_parnum + * Invalid number of model parameters found in XML element. + * @exception GException::model_invalid_parnames + * Invalid model parameter names found in XML element. + * + * Writes the radial ring model information into an XML element. The XML + * element will have the format + * + * + * + * + * + * + * + * + ***************************************************************************/ +void GModelSpatialRadialRing::write(GXmlElement& xml) const +{ + // Write ring center location + GModelSpatialRadial::write(xml); + + // Get or create parameters + GXmlElement* radius = gammalib::xml_need_par(G_WRITE, xml, m_radius.name()); + GXmlElement* width = gammalib::xml_need_par(G_WRITE, xml, m_width.name()); + + // Write parameters + m_radius.write(*radius); + m_width.write(*width); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Print information + * + * @param[in] chatter Chattiness (defaults to NORMAL). + * @return String containing model information. + ***************************************************************************/ +std::string GModelSpatialRadialRing::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("=== GModelSpatialRadialRing ==="); + + // 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 GModelSpatialRadialRing::init_members(void) +{ + // Initialise model type + m_type = "RadialRing"; + + // 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); // Radial components never have gradients + + // Set parameter pointer(s) + m_pars.push_back(&m_radius); + + // Initialise Width + m_width.clear(); + m_width.name("Width"); + m_width.unit("deg"); + m_width.value(2.778e-4); // 1 arcsec + m_width.min(2.778e-4); // 1 arcsec + m_width.free(); + m_width.scale(1.0); + m_width.gradient(0.0); + m_width.has_grad(false); // Radial components never have gradients + + // Set parameter pointer(s) + m_pars.push_back(&m_width); + + // Initialise precomputation cache + m_last_radius = 0.0; + m_last_width = 0.0; + m_inner_radius_rad = 0.0; + m_outer_radius_rad = 0.0; + m_width_rad = 0.0; + m_norm = 0.0; + + // Initialise other members + m_region.clear(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] model Radial ring 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 the radius twice on the stack. + ***************************************************************************/ +void GModelSpatialRadialRing::copy_members(const GModelSpatialRadialRing& model) +{ + // Copy members + m_type = model.m_type; + m_radius = model.m_radius; + m_width = model.m_width; + m_region = model.m_region; + + // Copy precomputation cache + m_last_radius = model.m_last_radius; + m_last_width = model.m_last_width; + m_inner_radius_rad = model.m_inner_radius_rad; + m_outer_radius_rad = model.m_outer_radius_rad; + m_width_rad = model.m_width_rad; + m_norm = model.m_norm; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void GModelSpatialRadialRing::free_members(void) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief Update precomputation cache + * + * Computes the normalization + * \f[ + * {\tt m\_norm} = \frac{1}{2 \pi (\cos r_i - \cos r_o)} + * \f] + * + * This is the correct normalization on the sphere for any radii. + ***************************************************************************/ +void GModelSpatialRadialRing::update() const +{ + // Update if radius or width has changed + if ((m_last_radius != radius()) || (m_last_width != width())) { + + // Store last values + m_last_radius = radius(); + m_last_width = width(); + + // Compute ring parameters in radians + m_outer_radius_rad = radius() * gammalib::deg2rad; + m_inner_radius_rad = (radius()-width()) * gammalib::deg2rad; + m_width_rad = width() * gammalib::deg2rad; + + // Perform precomputations + double denom = gammalib::twopi * + (std::cos(m_inner_radius_rad) - + std::cos(m_outer_radius_rad)); + m_norm = (denom > 0.0) ? 1.0 / denom : 0.0; + + } // endif: update required + + // Also update region + set_region(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set boundary sky region + ***************************************************************************/ +void GModelSpatialRadialRing::set_region(void) const +{ + // Set sky region centre to ring centre + m_region.centre(m_ra.value(), m_dec.value()); + + // Set sky region radius to ring radius + m_region.radius(m_radius.value()); + + // Return + return; +} diff --git a/src/model/Makefile.am b/src/model/Makefile.am index 9450590..297f9f5 100644 --- a/src/model/Makefile.am +++ b/src/model/Makefile.am @@ -36,6 +36,7 @@ sources = GModelPar.cpp \ GModelSpatialPointSource.cpp \ GModelSpatialRadial.cpp \ GModelSpatialRadialDisk.cpp \ + GModelSpatialRadialRing.cpp \ GModelSpatialRadialGauss.cpp \ GModelSpatialRadialShell.cpp \ GModelSpatialRadialProfile.cpp \ diff --git a/test/data/model_radial_ring.xml b/test/data/model_radial_ring.xml new file mode 100644 index 0000000..1e37b41 --- /dev/null +++ b/test/data/model_radial_ring.xml @@ -0,0 +1,16 @@ + + + + + + + + + + + + + + + + diff --git a/test/test_GModel.cpp b/test/test_GModel.cpp index f56ca6f..c355af3 100644 --- a/test/test_GModel.cpp +++ b/test/test_GModel.cpp @@ -76,6 +76,7 @@ void TestGModel::set(void) 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_elliptical_disk = datadir + "/model_elliptical_disk.xml"; @@ -113,6 +114,8 @@ void TestGModel::set(void) "Test GModelSpatialPointSource"); append(static_cast(&TestGModel::test_radial_disk), "Test GModelSpatialRadialDisk"); + append(static_cast(&TestGModel::test_radial_ring), + "Test GModelSpatialRadialRing"); append(static_cast(&TestGModel::test_radial_gauss), "Test GModelSpatialRadialGauss"); append(static_cast(&TestGModel::test_radial_shell), @@ -935,6 +938,88 @@ void TestGModel::test_radial_disk(void) /***********************************************************************//** + * @brief Test GModelSpatialRadialRing class + ***************************************************************************/ +void TestGModel::test_radial_ring(void) +{ + // Test void constructor + GModelSpatialRadialRing model1; + test_value(model1.type(), "RadialRing"); + + // Test value constructor + GSkyDir dir1; + dir1.radec_deg(83.6331, +22.0145); + GModelSpatialRadialRing model2(dir1, 3.0,1.0); + test_value(model2.ra(), 83.6331); + test_value(model2.dec(), 22.0145); + test_value(model2.radius(), 3.0); + test_value(model2.width(), 1.0); + + // Test XML constructor + GXml xml(m_xml_model_radial_ring); + GXmlElement* element = xml.element(0)->element(0)->element("spatialModel", 0); + GModelSpatialRadialRing model3(*element); + test_value(model3.size(), 4); + test_value(model3.type(), "RadialRing"); + test_value(model3.ra(), 83.6331); + test_value(model3.dec(), 22.0145); + test_value(model3.radius(), 0.45); + test_value(model3.width(), 0.15); + + // 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(1.0); + test_value(model3.radius(), 1.0); + + // Test width method + model3.width(0.5); + test_value(model3.width(), 0.5); + + // 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, +26.0); + test_assert(!model3.region()->contains(dir3), "Test region() method (outside)"); + + // Test contains method + dir3.radec_deg(83.6331, +22.75); + test_assert(model3.contains(dir3), "Test contains() method (inside)"); + dir3.radec_deg(83.6331, +22.0145); + test_assert(!model3.contains(dir3), "Test contains() method (outside interior of ring)"); + dir3.radec_deg(83.6331, +26.0); + test_assert(!model3.contains(dir3), "Test contains() method (outside exterior of ring)"); + + // Test operator access + const char* strarray[] = {"RA", "DEC", "Radius", "Width"}; + 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 GModelSpatialRadialGauss class ***************************************************************************/ void TestGModel::test_radial_gauss(void) @@ -2509,6 +2594,7 @@ void TestGModel::test_spatial_model(void) // Test spatial models XML interface test_xml_model("GModelSpatialPointSource", m_xml_model_point_plaw); test_xml_model("GModelSpatialRadialDisk", m_xml_model_radial_disk); + test_xml_model("GModelSpatialRadialRing", m_xml_model_radial_ring); test_xml_model("GModelSpatialRadialGauss", m_xml_model_radial_gauss); test_xml_model("GModelSpatialRadialShell", m_xml_model_radial_shell); test_xml_model("GModelSpatialEllipticalDisk", m_xml_model_elliptical_disk); diff --git a/test/test_GModel.hpp b/test/test_GModel.hpp index 2eb6602..6a831e6 100644 --- a/test/test_GModel.hpp +++ b/test/test_GModel.hpp @@ -54,6 +54,7 @@ public: void test_diffuse_cube(void); void test_diffuse_map(void); void test_radial_disk(void); + void test_radial_ring(void); void test_radial_gauss(void); void test_radial_shell(void); void test_elliptical_disk(void); @@ -133,6 +134,7 @@ private: std::string m_xml_model_diffuse_cube; std::string m_xml_model_diffuse_map; 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_shell; std::string m_xml_model_elliptical_disk; diff --git a/test/test_GModel.py b/test/test_GModel.py index adaa2f7..ddf8ecb 100755 --- a/test/test_GModel.py +++ b/test/test_GModel.py @@ -124,6 +124,7 @@ class Test(gammalib.GPythonTestSuite): test_support.pickeling(self, gammalib.GModelSpatialEllipticalGauss()) test_support.pickeling(self, gammalib.GModelSpatialPointSource()) test_support.pickeling(self, gammalib.GModelSpatialRadialDisk()) + test_support.pickeling(self, gammalib.GModelSpatialRadialRing()) test_support.pickeling(self, gammalib.GModelSpatialRadialGauss()) test_support.pickeling(self, gammalib.GModelSpatialRadialProfileDMBurkert()) test_support.pickeling(self, gammalib.GModelSpatialRadialProfileDMEinasto()) @@ -188,6 +189,7 @@ class Test(gammalib.GPythonTestSuite): test_support.pickeling(self, gammalib.GModelSpatialEllipticalGauss(dir,0.2,0.3,45.0)) test_support.pickeling(self, gammalib.GModelSpatialPointSource(ptsrc)) test_support.pickeling(self, gammalib.GModelSpatialRadialDisk(dir,0.2)) + test_support.pickeling(self, gammalib.GModelSpatialRadialRing(dir,0.3,0.1)) test_support.pickeling(self, gammalib.GModelSpatialRadialGauss(dir,0.2)) #test_support.pickeling(self, gammalib.GModelSpatialRadialProfileDMBurkert()) #test_support.pickeling(self, gammalib.GModelSpatialRadialProfileDMEinasto())