Maria Haupt / gammalib

Commit 045c40f4f63b336fedcad716201189da491db487

Authored by Jürgen Knödlseder
1 parent f5b90eec

Add GModelSpatialRadialProfile base class and GModelSpatialRadialProfileGauss derived class (#1520) 

The GModelSpatialRadialProfile class is the base class for an arbitrary radial model. The GModelSpatialRadialProfileGauss class basically illustrates its usage, and provides a comparison with the analytical Gaussian model, implemented by GModelSpatialRadialGauss.
Inline Side-by-side
Showing 12 changed files with 1267 additions and 3 deletions
ChangeLog
  View file @045c40f
1   -2016-02-28 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
  1 +2016-03-01 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
2 2  
3 3 * Version 1.1.0 released
4 4 =========================
5 5  
  6 + Add GModelSpatialRadialProfile and GModelSpatialRadialProfileGauss classes (#1520)
6 7 Make GFilename::is_fits thread safe (#1717).
7 8 Revise GTime operators (#1714)
8 9 Add GSkyMap::publish() method to publish map on VO Hub (#1710)
... ...
... ... @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 1.1.0
2 2  
3 3 Author(s) : Juergen Knoedlseder
4 4  
5   -28 February 2016
  5 +1 March 2016
6 6  
7 7  
8 8 1. Introduction
... ... @@ -142,6 +142,14 @@ GModelSpatialDiffuseMap class now recompute only the Monte Carlo acceptance
142 142 cone in case that the cone centre or radius have changed since the last
143 143 computation. This speeds up the computations (#1673).
144 144  
  145 +Add GModelSpatialRadialProfile class to handle arbitrary radial profiles.
  146 +The method deals with pre-computation and normalization of the profiles,
  147 +and implements the eval(), eval_grandients() and mc() method. The class
  148 +GModelSpatialRadialProfileGauss has been added to show how to implement
  149 +radial profile for the case of a Gaussian, and also to cross-check the
  150 +results with the GModelSpatialRadialGauss class that should return basically
  151 +the same results (#1520).
  152 +
145 153  
146 154 9. Numerics module
147 155 ------------------
... ...
include/GModelSpatialRadialProfile.hpp 0 → 100644
  View file @045c40f
  1 +/***************************************************************************
  2 + * GModelSpatialRadialProfile.hpp - Radial profile source model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2016 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GModelSpatialRadialProfile.hpp
  23 + * @brief Radial profile model class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GMODELSPATIALRADIALPROFILE_HPP
  28 +#define GMODELSPATIALRADIALPROFILE_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include "GModelSpatialRadial.hpp"
  33 +#include "GNodeArray.hpp"
  34 +
  35 +/* __ Forward declaration ________________________________________________ */
  36 +class GXmlElement;
  37 +class GSkyDir;
  38 +
  39 +
  40 +/**************************************************************************
  41 + * @class GModelSpatialRadialProfile
  42 + *
  43 + * @brief Radial profile source model class
  44 + *
  45 + * This class implements the spatial component of the factorised source
  46 + * model for an arbitrary radial profile.
  47 + ***************************************************************************/
  48 +class GModelSpatialRadialProfile : public GModelSpatialRadial {
  49 +
  50 +public:
  51 + // Constructors and destructors
  52 + GModelSpatialRadialProfile(void);
  53 + explicit GModelSpatialRadialProfile(const GXmlElement& xml);
  54 + GModelSpatialRadialProfile(const GModelSpatialRadialProfile& model);
  55 + virtual ~GModelSpatialRadialProfile(void);
  56 +
  57 + // Operators
  58 + virtual GModelSpatialRadialProfile& operator=(const GModelSpatialRadialProfile& model);
  59 +
  60 + // Pure virtual methods
  61 + virtual void clear(void) = 0;
  62 + virtual GModelSpatialRadialProfile* clone(void) const = 0;
  63 + virtual std::string classname(void) const = 0;
  64 + virtual std::string type(void) const = 0;
  65 + virtual double theta_max(void) const = 0;
  66 + virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
  67 +
  68 + // Implemented pure virtual base class methods
  69 + virtual double eval(const double& theta,
  70 + const GEnergy& energy,
  71 + const GTime& time) const;
  72 + virtual double eval_gradients(const double& theta,
  73 + const GEnergy& energy,
  74 + const GTime& time) const;
  75 + virtual GSkyDir mc(const GEnergy& energy,
  76 + const GTime& time,
  77 + GRan& ran) const;
  78 + virtual bool contains(const GSkyDir& dir,
  79 + const double& margin = 0.0) const;
  80 +
  81 + // Implement other methods
  82 + int num_nodes(void) const;
  83 + void num_nodes(const int& number);
  84 +
  85 +protected:
  86 + // Protected methods
  87 + void init_members(void);
  88 + void copy_members(const GModelSpatialRadialProfile& model);
  89 + void free_members(void);
  90 + int cache_index(void) const;
  91 + virtual double profile_value(const double& theta) const = 0;
  92 +
  93 + // Protected members
  94 + int m_num_nodes; //!< Number of profile nodes
  95 +
  96 + // Pre-computed radial profile
  97 + struct profile {
  98 + std::vector<double> pars; //!< Profile parameters
  99 + GNodeArray nodes; //!< Profile nodes
  100 + std::vector<double> values; //!< Profile values
  101 + std::vector<double> mc; //!< Profile for MC
  102 + double mc_max; //!< Maximum of profile for MC
  103 + };
  104 +
  105 + // Pre-computation cache
  106 + mutable std::vector<profile> m_profile; // Pre-computation cache
  107 +};
  108 +
  109 +
  110 +/***********************************************************************//**
  111 + * @brief Return number of nodes
  112 + *
  113 + * @return Number of nodes.
  114 + *
  115 + * Returns the number of nodes in the radial profile.
  116 + ***************************************************************************/
  117 +inline
  118 +int GModelSpatialRadialProfile::num_nodes(void) const
  119 +{
  120 + return (m_num_nodes);
  121 +}
  122 +
  123 +
  124 +/***********************************************************************//**
  125 + * @brief Set number of nodes
  126 + *
  127 + * @param[in] number Number of nodes.
  128 + *
  129 + * Sets the number of nodes in the radial profile.
  130 + ***************************************************************************/
  131 +inline
  132 +void GModelSpatialRadialProfile::num_nodes(const int& number)
  133 +{
  134 + m_num_nodes = number;
  135 + return;
  136 +}
  137 +
  138 +#endif /* GMODELSPATIALRADIALPROFILE_HPP */
... ...
include/GModelSpatialRadialProfileGauss.hpp 0 → 100644
  View file @045c40f
  1 +/***************************************************************************
  2 + * GModelSpatialRadialProfileGauss.hpp - Gaussian radial profile class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2016 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GModelSpatialRadialProfileGauss.hpp
  23 + * @brief Radial Gaussian profile model class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GMODELSPATIALRADIALPROFILEGAUSS_HPP
  28 +#define GMODELSPATIALRADIALPROFILEGAUSS_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include "GModelSpatialRadialProfile.hpp"
  33 +#include "GModelPar.hpp"
  34 +
  35 +/* __ Forward declaration ________________________________________________ */
  36 +class GXmlElement;
  37 +
  38 +
  39 +/**************************************************************************
  40 + * @class GModelSpatialRadialProfileGauss
  41 + *
  42 + * @brief Radial Gaussian profile source model class
  43 + *
  44 + * This class implements the spatial component of the factorised source
  45 + * model for a Gaussian radial profile.
  46 + ***************************************************************************/
  47 +class GModelSpatialRadialProfileGauss : public GModelSpatialRadialProfile {
  48 +
  49 +public:
  50 + // Constructors and destructors
  51 + GModelSpatialRadialProfileGauss(void);
  52 + explicit GModelSpatialRadialProfileGauss(const GXmlElement& xml);
  53 + GModelSpatialRadialProfileGauss(const GModelSpatialRadialProfileGauss& model);
  54 + virtual ~GModelSpatialRadialProfileGauss(void);
  55 +
  56 + // Operators
  57 + virtual GModelSpatialRadialProfileGauss& operator=(const GModelSpatialRadialProfileGauss& model);
  58 +
  59 + // Implemented pure virtual base class methods
  60 + virtual void clear(void);
  61 + virtual GModelSpatialRadialProfileGauss* clone(void) const;
  62 + virtual std::string classname(void) const;
  63 + virtual std::string type(void) const;
  64 + virtual double theta_max(void) const;
  65 + virtual void read(const GXmlElement& xml);
  66 + virtual void write(GXmlElement& xml) const;
  67 + virtual std::string print(const GChatter& chatter = NORMAL) const;
  68 +
  69 +protected:
  70 + // Protected methods
  71 + void init_members(void);
  72 + void copy_members(const GModelSpatialRadialProfileGauss& model);
  73 + void free_members(void);
  74 + virtual double profile_value(const double& theta) const;
  75 +
  76 + // Protected members
  77 + GModelPar m_sigma; //!< Gaussian width (deg)
  78 +};
  79 +
  80 +
  81 +/***********************************************************************//**
  82 + * @brief Return class name
  83 + *
  84 + * @return String containing the class name ("GModelSpatialRadialProfileGauss").
  85 + ***************************************************************************/
  86 +inline
  87 +std::string GModelSpatialRadialProfileGauss::classname(void) const
  88 +{
  89 + return ("GModelSpatialRadialProfileGauss");
  90 +}
  91 +
  92 +
  93 +/***********************************************************************//**
  94 + * @brief Return model type
  95 + *
  96 + * @return "GaussianProfile".
  97 + *
  98 + * Returns the type of the radial profile model.
  99 + ***************************************************************************/
  100 +inline
  101 +std::string GModelSpatialRadialProfileGauss::type(void) const
  102 +{
  103 + return "GaussianProfile";
  104 +}
  105 +
  106 +#endif /* GMODELSPATIALRADIALPROFILEGAUSS_HPP */
... ...
include/GammaLib.hpp
  View file @045c40f
1 1 /***************************************************************************
2 2 * GammaLib.hpp - Gamma-Ray Astronomy Library Header file *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2006-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2006-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -196,6 +196,8 @@
196 196 #include "GModelSpatialRadialDisk.hpp"
197 197 #include "GModelSpatialRadialGauss.hpp"
198 198 #include "GModelSpatialRadialShell.hpp"
  199 +#include "GModelSpatialRadialProfile.hpp"
  200 +#include "GModelSpatialRadialProfileGauss.hpp"
199 201 #include "GModelSpatialElliptical.hpp"
200 202 #include "GModelSpatialEllipticalDisk.hpp"
201 203 #include "GModelSpatialEllipticalGauss.hpp"
... ...
include/Makefile.am
  View file @045c40f
... ... @@ -133,6 +133,8 @@ pkginclude_HEADERS = GammaLib.hpp \
133 133 GModelSpatialRadialDisk.hpp \
134 134 GModelSpatialRadialGauss.hpp \
135 135 GModelSpatialRadialShell.hpp \
  136 + GModelSpatialRadialProfile.hpp \
  137 + GModelSpatialRadialProfileGauss.hpp \
136 138 GModelSpatialElliptical.hpp \
137 139 GModelSpatialEllipticalDisk.hpp \
138 140 GModelSpatialEllipticalGauss.hpp \
... ...
pyext/GModelSpatialRadialProfile.i 0 → 100644
  View file @045c40f
  1 +/***************************************************************************
  2 + * GModelSpatialRadialProfile.i - Radial profile source model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2016 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GModelSpatialRadialProfile.i
  23 + * @brief Radial profile model class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GModelSpatialRadialProfile.hpp"
  29 +%}
  30 +
  31 +
  32 +/**************************************************************************
  33 + * @class GModelSpatialRadialProfile
  34 + *
  35 + * @brief Radial profile source model class
  36 + *
  37 + * This class implements the spatial component of the factorised source
  38 + * model for an arbitrary radial profile.
  39 + ***************************************************************************/
  40 +class GModelSpatialRadialProfile : public GModelSpatialRadial {
  41 +
  42 +public:
  43 + // Constructors and destructors
  44 + GModelSpatialRadialProfile(void);
  45 + explicit GModelSpatialRadialProfile(const GXmlElement& xml);
  46 + GModelSpatialRadialProfile(const GModelSpatialRadialProfile& model);
  47 + virtual ~GModelSpatialRadialProfile(void);
  48 +
  49 + // Pure virtual methods
  50 + virtual void clear(void) = 0;
  51 + virtual GModelSpatialRadialProfile* clone(void) const = 0;
  52 + virtual std::string classname(void) const = 0;
  53 + virtual std::string type(void) const = 0;
  54 + virtual double theta_max(void) const = 0;
  55 +
  56 + // Implemented pure virtual base class methods
  57 + virtual double eval(const double& theta,
  58 + const GEnergy& energy,
  59 + const GTime& time) const;
  60 + virtual double eval_gradients(const double& theta,
  61 + const GEnergy& energy,
  62 + const GTime& time) const;
  63 + virtual GSkyDir mc(const GEnergy& energy,
  64 + const GTime& time,
  65 + GRan& ran) const;
  66 + virtual bool contains(const GSkyDir& dir,
  67 + const double& margin = 0.0) const;
  68 +
  69 + // Implement other methods
  70 + int num_nodes(void) const;
  71 + void num_nodes(const int& number);
  72 +};
... ...
pyext/GModelSpatialRadialProfileGauss.i 0 → 100644
  View file @045c40f
  1 +/***************************************************************************
  2 + * GModelSpatialRadialProfileGauss.i - Gaussian radial profile class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2016 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GModelSpatialRadialProfileGauss.hpp
  23 + * @brief Radial Gaussian profile model class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GModelSpatialRadialProfileGauss.hpp"
  29 +%}
  30 +
  31 +
  32 +/**************************************************************************
  33 + * @class GModelSpatialRadialProfileGauss
  34 + *
  35 + * @brief Radial Gaussian profile source model class
  36 + *
  37 + * This class implements the spatial component of the factorised source
  38 + * model for a Gaussian radial profile.
  39 + ***************************************************************************/
  40 +class GModelSpatialRadialProfileGauss : public GModelSpatialRadialProfile {
  41 +
  42 +public:
  43 + // Constructors and destructors
  44 + GModelSpatialRadialProfileGauss(void);
  45 + explicit GModelSpatialRadialProfileGauss(const GXmlElement& xml);
  46 + GModelSpatialRadialProfileGauss(const GModelSpatialRadialProfileGauss& model);
  47 + virtual ~GModelSpatialRadialProfileGauss(void);
  48 +
  49 + // Implemented pure virtual base class methods
  50 + virtual void clear(void);
  51 + virtual GModelSpatialRadialProfileGauss* clone(void) const;
  52 + virtual std::string classname(void) const;
  53 + virtual std::string type(void) const;
  54 + virtual double theta_max(void) const;
  55 + virtual void read(const GXmlElement& xml);
  56 + virtual void write(GXmlElement& xml) const;
  57 +};
  58 +
  59 +
  60 +/***********************************************************************//**
  61 + * @brief GModelSpatialRadialGauss class extension
  62 + *
  63 + * The eval(GSkyDir&) and eval_gradients(GSkyDir&) need to be defined in the
  64 + * extension to force swig to build also the interface for these methods that
  65 + * are implemented in the base class only. It's not clear to me why these
  66 + * methods are not inherited automatically. Maybe this could also be handled
  67 + * by a %typemap(typecheck) construct.
  68 + ***************************************************************************/
  69 +%extend GModelSpatialRadialProfileGauss {
  70 + double eval(const GPhoton& photon) const {
  71 + return self->GModelSpatialRadial::eval(photon);
  72 + }
  73 + double eval_gradients(const GPhoton& photon) const {
  74 + return self->GModelSpatialRadial::eval_gradients(photon);
  75 + }
  76 +};
... ...
pyext/gammalib/model.i
  View file @045c40f
... ... @@ -134,6 +134,8 @@
134 134 %include "GModelSpatialRadialDisk.i"
135 135 %include "GModelSpatialRadialGauss.i"
136 136 %include "GModelSpatialRadialShell.i"
  137 +%include "GModelSpatialRadialProfile.i"
  138 +%include "GModelSpatialRadialProfileGauss.i"
137 139 %include "GModelSpatialElliptical.i"
138 140 %include "GModelSpatialEllipticalDisk.i"
139 141 %include "GModelSpatialEllipticalGauss.i"
... ...
src/model/GModelSpatialRadialProfile.cpp 0 → 100644
  View file @045c40f
  1 +/***************************************************************************
  2 + * GModelSpatialRadialProfile.cpp - Radial profile source model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2016 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GModelSpatialRadialProfile.cpp
  23 + * @brief Radial profile model class implementation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GException.hpp"
  32 +#include "GTools.hpp"
  33 +#include "GMath.hpp"
  34 +#include "GSkyDir.hpp"
  35 +#include "GXmlElement.hpp"
  36 +#include "GModelSpatialRadialProfile.hpp"
  37 +
  38 +/* __ Constants __________________________________________________________ */
  39 +
  40 +/* __ Globals ____________________________________________________________ */
  41 +
  42 +/* __ Method name definitions ____________________________________________ */
  43 +
  44 +/* __ Macros _____________________________________________________________ */
  45 +
  46 +/* __ Coding definitions _________________________________________________ */
  47 +
  48 +/* __ Debug definitions __________________________________________________ */
  49 +//#define G_DEBUG_PRECOMPUTATION
  50 +
  51 +
  52 +/*==========================================================================
  53 + = =
  54 + = Constructors/destructors =
  55 + = =
  56 + ==========================================================================*/
  57 +
  58 +/***********************************************************************//**
  59 + * @brief Void constructor
  60 + *
  61 + * Constructs empty radial profile
  62 + ***************************************************************************/
  63 +GModelSpatialRadialProfile::GModelSpatialRadialProfile(void) :
  64 + GModelSpatialRadial()
  65 +{
  66 + // Initialise members
  67 + init_members();
  68 +
  69 + // Return
  70 + return;
  71 +}
  72 +
  73 +
  74 +/***********************************************************************//**
  75 + * @brief XML constructor
  76 + *
  77 + * @param[in] xml XML element.
  78 + *
  79 + * Constructs radial profile model by extracting information from an XML
  80 + * elements.
  81 + ***************************************************************************/
  82 +GModelSpatialRadialProfile::GModelSpatialRadialProfile(const GXmlElement& xml) :
  83 + GModelSpatialRadial()
  84 +{
  85 + // Initialise members
  86 + init_members();
  87 +
  88 + // Read information from XML element
  89 + read(xml);
  90 +
  91 + // Return
  92 + return;
  93 +}
  94 +
  95 +
  96 +/***********************************************************************//**
  97 + * @brief Copy constructor
  98 + *
  99 + * @param[in] model Radial profile model.
  100 + *
  101 + * Copies radial profile model from another radial profile model.
  102 + ***************************************************************************/
  103 +GModelSpatialRadialProfile::GModelSpatialRadialProfile(const GModelSpatialRadialProfile& model) :
  104 + GModelSpatialRadial(model)
  105 +{
  106 + // Initialise members
  107 + init_members();
  108 +
  109 + // Copy members
  110 + copy_members(model);
  111 +
  112 + // Return
  113 + return;
  114 +}
  115 +
  116 +
  117 +/***********************************************************************//**
  118 + * @brief Destructor
  119 + *
  120 + * Destructs radial profile model.
  121 + ***************************************************************************/
  122 +GModelSpatialRadialProfile::~GModelSpatialRadialProfile(void)
  123 +{
  124 + // Free members
  125 + free_members();
  126 +
  127 + // Return
  128 + return;
  129 +}
  130 +
  131 +
  132 +/*==========================================================================
  133 + = =
  134 + = Operators =
  135 + = =
  136 + ==========================================================================*/
  137 +
  138 +/***********************************************************************//**
  139 + * @brief Assignment operator
  140 + *
  141 + * @param[in] model Radial profile model.
  142 + * @return Radial profile model.
  143 + *
  144 + * Assigns radial profile model.
  145 + ***************************************************************************/
  146 +GModelSpatialRadialProfile& GModelSpatialRadialProfile::operator=(const GModelSpatialRadialProfile& model)
  147 +{
  148 + // Execute only if object is not identical
  149 + if (this != &model) {
  150 +
  151 + // Copy base class members
  152 + this->GModelSpatialRadial::operator=(model);
  153 +
  154 + // Free members
  155 + free_members();
  156 +
  157 + // Initialise members
  158 + init_members();
  159 +
  160 + // Copy members
  161 + copy_members(model);
  162 +
  163 + } // endif: object was not identical
  164 +
  165 + // Return
  166 + return *this;
  167 +}
  168 +
  169 +
  170 +/*==========================================================================
  171 + = =
  172 + = Public methods =
  173 + = =
  174 + ==========================================================================*/
  175 +
  176 +/***********************************************************************//**
  177 + * @brief Evaluate function (in units of sr^-1)
  178 + *
  179 + * @param[in] theta Angular distance from model centre (radians).
  180 + * @param[in] energy Photon energy.
  181 + * @param[in] time Photon arrival time.
  182 + * @return Model value.
  183 + *
  184 + * Evaluate the radial profile model for a given angular distance @p theta
  185 + * the model centre, a given photon @p energy and a given @p time. The
  186 + * method evaluates the model by linear interpolation.
  187 + ***************************************************************************/
  188 +double GModelSpatialRadialProfile::eval(const double& theta,
  189 + const GEnergy& energy,
  190 + const GTime& time) const
  191 +{
  192 + // Get pre-computation cache index
  193 + int icache = cache_index();
  194 +
  195 + // Get interpolation value
  196 + double value = m_profile[icache].nodes.interpolate(theta, m_profile[icache].values);
  197 +
  198 + // Make sure that value is not-negative
  199 + if (value < 0.0) {
  200 + value = 0.0;
  201 + }
  202 +
  203 + // Compile option: Check for NaN/Inf
  204 + #if defined(G_NAN_CHECK)
  205 + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
  206 + std::cout << "*** ERROR: GModelSpatialRadialProfile::eval";
  207 + std::cout << "(theta=" << theta << "): NaN/Inf encountered";
  208 + std::cout << " (value=" << value;
  209 + std::cout << ")" << std::endl;
  210 + }
  211 + #endif
  212 +
  213 + // Return value
  214 + return value;
  215 +}
  216 +
  217 +
  218 +/***********************************************************************//**
  219 + * @brief Evaluate function and gradients (in units of sr^-1)
  220 + *
  221 + * @param[in] theta Angular distance from model centre (radians).
  222 + * @param[in] energy Photon energy.
  223 + * @param[in] time Photon arrival time.
  224 + * @return Model value.
  225 + *
  226 + * Evaluates the function value. No gradient computation is implemented as
  227 + * radial models will be convolved with the instrument response and thus
  228 + * require the numerical computation of the derivatives.
  229 + *
  230 + * See the eval() method for more information.
  231 + ***************************************************************************/
  232 +double GModelSpatialRadialProfile::eval_gradients(const double& theta,
  233 + const GEnergy& energy,
  234 + const GTime& time) const
  235 +{
  236 + // Return value
  237 + return (eval(theta, energy, time));
  238 +}
  239 +
  240 +
  241 +/***********************************************************************//**
  242 + * @brief Return MC sky direction
  243 + *
  244 + * @param[in] energy Photon energy.
  245 + * @param[in] time Photon arrival time.
  246 + * @param[in,out] ran Random number generator.
  247 + * @return Sky direction.
  248 + *
  249 + * Draws an arbitrary sky direction from the radial profile using a rejection
  250 + * method.
  251 + ***************************************************************************/
  252 +GSkyDir GModelSpatialRadialProfile::mc(const GEnergy& energy,
  253 + const GTime& time,
  254 + GRan& ran) const
  255 +{
  256 + // Get pre-computation cache index
  257 + int icache = cache_index();
  258 +
  259 + // Simulate theta angle (degrees) using a rejection method
  260 + double theta = 0.0;
  261 + while (true) {
  262 + theta = ran.uniform() * theta_max();
  263 + double mc = m_profile[icache].nodes.interpolate(theta, m_profile[icache].mc);
  264 + if ((ran.uniform() * m_profile[icache].mc_max) <= mc) {
  265 + break;
  266 + }
  267 + }
  268 + theta *= gammalib::rad2deg;
  269 +
  270 + // Simulate uniform azmiuth angle (degrees)
  271 + double phi = 360.0 * ran.uniform();
  272 +
  273 + // Rotate sky direction by offset
  274 + GSkyDir sky_dir = dir();
  275 + sky_dir.rotate_deg(phi, theta);
  276 +
  277 + // Return sky direction
  278 + return sky_dir;
  279 +}
  280 +
  281 +
  282 +/***********************************************************************//**
  283 + * @brief Checks where model contains specified sky direction
  284 + *
  285 + * @param[in] dir Sky direction.
  286 + * @param[in] margin Margin to be added to sky direction (degrees)
  287 + * @return True if the model contains the sky direction.
  288 + *
  289 + * Signals whether a sky direction is contained in the radial disk model.
  290 + ***************************************************************************/
  291 +bool GModelSpatialRadialProfile::contains(const GSkyDir& dir,
  292 + const double& margin) const
  293 +{
  294 + // Compute distance to centre (radians)
  295 + double distance = dir.dist(this->dir());
  296 +
  297 + // Return flag
  298 + return (distance <= theta_max() + margin*gammalib::deg2rad);
  299 +}
  300 +
  301 +
  302 +/*==========================================================================
  303 + = =
  304 + = Private methods =
  305 + = =
  306 + ==========================================================================*/
  307 +
  308 +/***********************************************************************//**
  309 + * @brief Initialise class members
  310 + ***************************************************************************/
  311 +void GModelSpatialRadialProfile::init_members(void)
  312 +{
  313 + // Initialise members
  314 + m_num_nodes = 100;
  315 +
  316 + // Initialise pre-computation cache
  317 + m_profile.clear();
  318 +
  319 + // Return
  320 + return;
  321 +}
  322 +
  323 +
  324 +/***********************************************************************//**
  325 + * @brief Copy class members
  326 + *
  327 + * @param[in] model Radial disk model.
  328 + ***************************************************************************/
  329 +void GModelSpatialRadialProfile::copy_members(const GModelSpatialRadialProfile& model)
  330 +{
  331 + // Copy members
  332 + m_num_nodes = model.m_num_nodes;
  333 + m_profile = model.m_profile;
  334 +
  335 + // Return
  336 + return;
  337 +}
  338 +
  339 +
  340 +/***********************************************************************//**
  341 + * @brief Delete class members
  342 + ***************************************************************************/
  343 +void GModelSpatialRadialProfile::free_members(void)
  344 +{
  345 + // Return
  346 + return;
  347 +}
  348 +
  349 +
  350 +/***********************************************************************//**
  351 + * @brief Return index to pre-computation cache
  352 + *
  353 + * @return Index to pre-computation cache
  354 + *
  355 + * Returns the index to the pre-computation cache. If no pre-computation
  356 + * cache was found the method will create one and return the index to that
  357 + * cache.
  358 + ***************************************************************************/
  359 +int GModelSpatialRadialProfile::cache_index(void) const
  360 +{
  361 + // Initialise index
  362 + int index = -1;
  363 +
  364 + // Search for index
  365 + for (int i = 0; i < m_profile.size(); ++i) {
  366 +
  367 + // Initialise found flag
  368 + bool found = true;
  369 +
  370 + // Loop over all spatial parameters and check if they are equal
  371 + // to the parameters in the pre-computation cache. Break and set
  372 + // the found flag to false on non-equality.
  373 + for (int k = 0; k < m_pars.size(); ++k) {
  374 + if (m_pars[k]->value() != m_profile[i].pars[k]) {
  375 + found = false;
  376 + break;
  377 + }
  378 + }
  379 +
  380 + // If profile was found then set index
  381 + if (found) {
  382 + index = i;
  383 + break;
  384 + }
  385 +
  386 + } // endfor: looped over index
  387 +
  388 + // If index was not found then allocate a new pre-computation cache
  389 + if (index == -1) {
  390 +
  391 + // Initialise profile
  392 + profile prf;
  393 + prf.nodes.clear();
  394 + prf.values.clear();
  395 + prf.mc.clear();
  396 + prf.mc_max = 0.0;
  397 +
  398 + // Set profile parameters
  399 + for (int k = 0; k < m_pars.size(); ++k) {
  400 + prf.pars.push_back(m_pars[k]->value());
  401 + }
  402 +
  403 + // Pre-compute radial values, MC values, and mc_max. Compute
  404 + // also the normalization.
  405 + double rmax = theta_max();
  406 + double dr = rmax / m_num_nodes;
  407 + double r = 0.0;
  408 + double norm = 0.0;
  409 + for (int j = 0; j < m_num_nodes; ++j) {
  410 + double value = profile_value(r);
  411 + double mc = value * std::sin(r) * dr;
  412 + norm += mc;
  413 + if (mc > prf.mc_max) {
  414 + prf.mc_max = mc;
  415 + }
  416 + prf.nodes.append(r);
  417 + prf.values.push_back(value);
  418 + prf.mc.push_back(mc);
  419 + r += dr;
  420 + }
  421 + norm *= gammalib::twopi;
  422 +
  423 + // Normalize radial profile
  424 + if (norm > 0.0) {
  425 + double inv_norm = 1.0 / norm;
  426 + for (int j = 0; j < m_num_nodes; ++j) {
  427 + prf.values[j] *= inv_norm;
  428 + }
  429 + }
  430 +
  431 + // Push profile in pre-computation cache
  432 + m_profile.push_back(prf);
  433 +
  434 + // Get index to last profile
  435 + index = m_profile.size() - 1;
  436 +
  437 + // Log pre-computation
  438 + #if defined(G_DEBUG_PRECOMPUTATION)
  439 + std::cout << "GModelSpatialRadialProfile::cache_index";
  440 + std::cout << std::endl;
  441 + std::cout << " dr=" << dr << std::endl;
  442 + for (int k = 0; k < m_pars.size(); ++k) {
  443 + std::cout << " par[" << k << "]=" << m_pars[k]->value();
  444 + std::cout << " (" << m_pars[k]->name();
  445 + std::cout << ")" << std::endl;
  446 + }
  447 + std::cout << " norm=" << norm << std::endl;
  448 + #endif
  449 +
  450 + } // endif: allocated new pre-computation cache
  451 +
  452 + // Return index
  453 + return (index);
  454 +}
... ...
src/model/GModelSpatialRadialProfileGauss.cpp 0 → 100644
  View file @045c40f
  1 +/***************************************************************************
  2 + * GModelSpatialRadialProfileGauss.cpp - Gaussian radial profile class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2016 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GModelSpatialRadialProfileGauss.cpp
  23 + * @brief Radial Gaussian profile model class implementation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GException.hpp"
  32 +#include "GTools.hpp"
  33 +#include "GMath.hpp"
  34 +#include "GXmlElement.hpp"
  35 +#include "GModelSpatialRadialProfileGauss.hpp"
  36 +#include "GModelSpatialRegistry.hpp"
  37 +
  38 +/* __ Constants __________________________________________________________ */
  39 +
  40 +/* __ Globals ____________________________________________________________ */
  41 +const GModelSpatialRadialProfileGauss g_radial_disk_seed;
  42 +const GModelSpatialRegistry g_radial_disk_registry(&g_radial_disk_seed);
  43 +
  44 +/* __ Method name definitions ____________________________________________ */
  45 +#define G_READ "GModelSpatialRadialProfileGauss::read(GXmlElement&)"
  46 +#define G_WRITE "GModelSpatialRadialProfileGauss::write(GXmlElement&)"
  47 +
  48 +/* __ Macros _____________________________________________________________ */
  49 +
  50 +/* __ Coding definitions _________________________________________________ */
  51 +
  52 +/* __ Debug definitions __________________________________________________ */
  53 +
  54 +
  55 +/*==========================================================================
  56 + = =
  57 + = Constructors/destructors =
  58 + = =
  59 + ==========================================================================*/
  60 +
  61 +/***********************************************************************//**
  62 + * @brief Void constructor
  63 + *
  64 + * Constructs empty radial Gaussian profile
  65 + ***************************************************************************/
  66 +GModelSpatialRadialProfileGauss::GModelSpatialRadialProfileGauss(void) :
  67 + GModelSpatialRadialProfile()
  68 +{
  69 + // Initialise members
  70 + init_members();
  71 +
  72 + // Return
  73 + return;
  74 +}
  75 +
  76 +
  77 +/***********************************************************************//**
  78 + * @brief XML constructor
  79 + *
  80 + * @param[in] xml XML element.
  81 + *
  82 + * Constructs radial Gaussian profile model by extracting information from
  83 + * an XML element. See the read() method for more information about the
  84 + * expected structure of the XML element.
  85 + ***************************************************************************/
  86 +GModelSpatialRadialProfileGauss::GModelSpatialRadialProfileGauss(const GXmlElement& xml) :
  87 + GModelSpatialRadialProfile()
  88 +{
  89 + // Initialise members
  90 + init_members();
  91 +
  92 + // Read information from XML element
  93 + read(xml);
  94 +
  95 + // Return
  96 + return;
  97 +}
  98 +
  99 +
  100 +/***********************************************************************//**
  101 + * @brief Copy constructor
  102 + *
  103 + * @param[in] model Radial Gaussian profile model.
  104 + *
  105 + * Copies radial Gaussian profile model from another radial profile model.
  106 + ***************************************************************************/
  107 +GModelSpatialRadialProfileGauss::GModelSpatialRadialProfileGauss(const GModelSpatialRadialProfileGauss& model) :
  108 + GModelSpatialRadialProfile(model)
  109 +{
  110 + // Initialise members
  111 + init_members();
  112 +
  113 + // Copy members
  114 + copy_members(model);
  115 +
  116 + // Return
  117 + return;
  118 +}
  119 +
  120 +
  121 +/***********************************************************************//**
  122 + * @brief Destructor
  123 + *
  124 + * Destructs radial Gaussian profile model.
  125 + ***************************************************************************/
  126 +GModelSpatialRadialProfileGauss::~GModelSpatialRadialProfileGauss(void)
  127 +{
  128 + // Free members
  129 + free_members();
  130 +
  131 + // Return
  132 + return;
  133 +}
  134 +
  135 +
  136 +/*==========================================================================
  137 + = =
  138 + = Operators =
  139 + = =
  140 + ==========================================================================*/
  141 +
  142 +/***********************************************************************//**
  143 + * @brief Assignment operator
  144 + *
  145 + * @param[in] model Radial Gaussian profile model.
  146 + * @return Radial Gaussian profile model.
  147 + *
  148 + * Assigns radial Gaussian profile model.
  149 + ***************************************************************************/
  150 +GModelSpatialRadialProfileGauss& GModelSpatialRadialProfileGauss::operator=(const GModelSpatialRadialProfileGauss& model)
  151 +{
  152 + // Execute only if object is not identical
  153 + if (this != &model) {
  154 +
  155 + // Copy base class members
  156 + this->GModelSpatialRadialProfile::operator=(model);
  157 +
  158 + // Free members
  159 + free_members();
  160 +
  161 + // Initialise members
  162 + init_members();
  163 +
  164 + // Copy members
  165 + copy_members(model);
  166 +
  167 + } // endif: object was not identical
  168 +
  169 + // Return
  170 + return *this;
  171 +}
  172 +
  173 +
  174 +/*==========================================================================
  175 + = =
  176 + = Public methods =
  177 + = =
  178 + ==========================================================================*/
  179 +
  180 +/***********************************************************************//**
  181 + * @brief Clear radial Gaussian profile model
  182 + *
  183 + * Clears radial Gaussian profile model.
  184 + ***************************************************************************/
  185 +void GModelSpatialRadialProfileGauss::clear(void)
  186 +{
  187 + // Free class members
  188 + free_members();
  189 + this->GModelSpatialRadialProfile::free_members();
  190 + this->GModelSpatialRadial::free_members();
  191 + this->GModelSpatial::free_members();
  192 +
  193 + // Initialise members
  194 + this->GModelSpatial::init_members();
  195 + this->GModelSpatialRadial::init_members();
  196 + this->GModelSpatialRadialProfile::init_members();
  197 + init_members();
  198 +
  199 + // Return
  200 + return;
  201 +}
  202 +
  203 +
  204 +/***********************************************************************//**
  205 + * @brief Clone radial Gaussian profile model
  206 + *
  207 + * @return Pointer to deep copy of radial Gaussian profile model.
  208 + *
  209 + * Returns a deep copy of the radial Gaussian profile model.
  210 + ***************************************************************************/
  211 +GModelSpatialRadialProfileGauss* GModelSpatialRadialProfileGauss::clone(void) const
  212 +{
  213 + // Clone radial disk model
  214 + return new GModelSpatialRadialProfileGauss(*this);
  215 +}
  216 +
  217 +
  218 +/***********************************************************************//**
  219 + * @brief Return maximum model radius (in radians)
  220 + *
  221 + * @return Maximum model radius (in radians).
  222 + ***************************************************************************/
  223 +double GModelSpatialRadialProfileGauss::theta_max(void) const
  224 +{
  225 + // Return value
  226 + return (m_sigma.value() * gammalib::deg2rad * 5.0);
  227 +}
  228 +
  229 +
  230 +/***********************************************************************//**
  231 + * @brief Read model from XML element
  232 + *
  233 + * @param[in] xml XML element.
  234 + *
  235 + * Reads the Gaussian radial profile model information from an XML element.
  236 + * The XML element shall have either the format
  237 + *
  238 + * <spatialModel type="GaussianProfile">
  239 + * <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
  240 + * <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
  241 + * <parameter name="Sigma" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
  242 + * </spatialModel>
  243 + *
  244 + * or
  245 + *
  246 + * <spatialModel type="GaussianProfile">
  247 + * <parameter name="GLON" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
  248 + * <parameter name="GLAT" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
  249 + * <parameter name="Sigma" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
  250 + * </spatialModel>
  251 + ***************************************************************************/
  252 +void GModelSpatialRadialProfileGauss::read(const GXmlElement& xml)
  253 +{
  254 + // Read Gaussian location
  255 + GModelSpatialRadial::read(xml);
  256 +
  257 + // Read Sigma parameter
  258 + const GXmlElement* par = gammalib::xml_get_par(G_READ, xml, "Sigma");
  259 + m_sigma.read(*par);
  260 +
  261 + // Return
  262 + return;
  263 +}
  264 +
  265 +
  266 +/***********************************************************************//**
  267 + * @brief Write model into XML element
  268 + *
  269 + * @param[in] xml XML element into which model information is written.
  270 + *
  271 + * Writes the Gaussian radial profile model information into an XML element.
  272 + * The XML element will have the format
  273 + *
  274 + * <spatialModel type="GaussianProfile">
  275 + * <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/>
  276 + * <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/>
  277 + * <parameter name="Sigma" scale="1.0" value="0.45" min="0.01" max="10" free="1"/>
  278 + * </spatialModel>
  279 + ***************************************************************************/
  280 +void GModelSpatialRadialProfileGauss::write(GXmlElement& xml) const
  281 +{
  282 + // Write Gaussian location
  283 + GModelSpatialRadial::write(xml);
  284 +
  285 + // Write Sigma parameter
  286 + GXmlElement* par = gammalib::xml_need_par(G_WRITE, xml, "Sigma");
  287 + m_sigma.write(*par);
  288 +
  289 + // Return
  290 + return;
  291 +}
  292 +
  293 +
  294 +/***********************************************************************//**
  295 + * @brief Print information
  296 + *
  297 + * @param[in] chatter Chattiness (defaults to NORMAL).
  298 + * @return String containing model information.
  299 + ***************************************************************************/
  300 +std::string GModelSpatialRadialProfileGauss::print(const GChatter& chatter) const
  301 +{
  302 + // Initialise result string
  303 + std::string result;
  304 +
  305 + // Continue only if chatter is not silent
  306 + if (chatter != SILENT) {
  307 +
  308 + // Append header
  309 + result.append("=== GModelSpatialRadialProfileGauss ===");
  310 +
  311 + // Append parameters
  312 + result.append("\n"+gammalib::parformat("Number of parameters"));
  313 + result.append(gammalib::str(size()));
  314 + for (int i = 0; i < size(); ++i) {
  315 + result.append("\n"+m_pars[i]->print(chatter));
  316 + }
  317 +
  318 + } // endif: chatter was not silent
  319 +
  320 + // Return result
  321 + return result;
  322 +}
  323 +
  324 +
  325 +/*==========================================================================
  326 + = =
  327 + = Private methods =
  328 + = =
  329 + ==========================================================================*/
  330 +
  331 +/***********************************************************************//**
  332 + * @brief Initialise class members
  333 + ***************************************************************************/
  334 +void GModelSpatialRadialProfileGauss::init_members(void)
  335 +{
  336 + // Initialise Gaussian sigma
  337 + m_sigma.clear();
  338 + m_sigma.name("Sigma");
  339 + m_sigma.unit("deg");
  340 + m_sigma.value(2.778e-4); // 1 arcsec
  341 + m_sigma.min(2.778e-4); // 1 arcsec
  342 + m_sigma.free();
  343 + m_sigma.scale(1.0);
  344 + m_sigma.gradient(0.0);
  345 + m_sigma.has_grad(false); // Radial components never have gradients
  346 +
  347 + // Set parameter pointer(s)
  348 + m_pars.push_back(&m_sigma);
  349 +
  350 + // Return
  351 + return;
  352 +}
  353 +
  354 +
  355 +/***********************************************************************//**
  356 + * @brief Copy class members
  357 + *
  358 + * @param[in] model Radial Gaussian model.
  359 + *
  360 + * Copies class members from another radial profile model.
  361 + ***************************************************************************/
  362 +void GModelSpatialRadialProfileGauss::copy_members(const GModelSpatialRadialProfileGauss& model)
  363 +{
  364 + // Copy members. We do not have to push back the members on the parameter
  365 + // stack as this should have been done by init_members() that was called
  366 + // before. Otherwise we would have sigma twice on the stack.
  367 + m_sigma = model.m_sigma;
  368 +
  369 + // Return
  370 + return;
  371 +}
  372 +
  373 +
  374 +/***********************************************************************//**
  375 + * @brief Delete class members
  376 + ***************************************************************************/
  377 +void GModelSpatialRadialProfileGauss::free_members(void)
  378 +{
  379 + // Return
  380 + return;
  381 +}
  382 +
  383 +
  384 +/***********************************************************************//**
  385 + * @brief Radial profile
  386 + *
  387 + * @param[in] theta Angular distance from Gaussian centre (radians).
  388 + * @return Profile value.
  389 + ***************************************************************************/
  390 +double GModelSpatialRadialProfileGauss::profile_value(const double& theta) const
  391 +{
  392 + // Compute value
  393 + double sigma_rad = m_sigma.value() * gammalib::deg2rad;
  394 + double sigma2 = sigma_rad * sigma_rad;
  395 + double theta2 = theta * theta;
  396 + double value = std::exp(-0.5 * theta2 / sigma2) /
  397 + (gammalib::twopi * sigma2);
  398 +
  399 + // Return value
  400 + return value;
  401 +}
... ...
src/model/Makefile.am
  View file @045c40f
... ... @@ -19,6 +19,8 @@ sources = GModelPar.cpp \
19 19 GModelSpatialRadialDisk.cpp \
20 20 GModelSpatialRadialGauss.cpp \
21 21 GModelSpatialRadialShell.cpp \
  22 + GModelSpatialRadialProfile.cpp \
  23 + GModelSpatialRadialProfileGauss.cpp \
22 24 GModelSpatialElliptical.cpp \
23 25 GModelSpatialEllipticalDisk.cpp \
24 26 GModelSpatialEllipticalGauss.cpp \
... ...