Jürgen Knödlseder / gammalib

Commit d1cd172c96bc1bc466d973fc4a5a3bbea56b5db8

Authored by Michael Mayer
1 parent 53a83c05

I added Feature 605 - LogParabola spectrum

Inline Side-by-side
Showing 7 changed files with 949 additions and 12 deletions
include/GModelSpectralLogParabola.hpp 0 → 100644
  View file @d1cd172
  1 +/*
  2 + * GModelSpectralLogParabola.hpp
  3 + *
  4 + * Created on: Nov 30, 2012
  5 + * Author: mmayer
  6 + */
  7 +
  8 +#ifndef GMODELSPECTRALLOGPARABOLA_HPP_
  9 +#define GMODELSPECTRALLOGPARABOLA_HPP_
  10 +
  11 +
  12 +/* __ Includes ___________________________________________________________ */
  13 +#include <string>
  14 +#include "GModelPar.hpp"
  15 +#include "GModelSpectral.hpp"
  16 +#include "GEnergy.hpp"
  17 +#include "GXmlElement.hpp"
  18 +#include "GIntegral.hpp"
  19 +
  20 +
  21 +/***********************************************************************//**
  22 + * @class GModelSpectralLogParabola
  23 + *
  24 + * @brief LogParabola spectrall model class
  25 + *
  26 + * This class implements a logparabole as the spectral component of the
  27 + * gamma-ray sky model. The power law is defined as
  28 + * \f[I(E)=norm (E/pivot)^{index-curvature\ln(E/pivot)}\f]
  29 + * where
  30 + * \f$norm\f$ is the normalization or prefactor,
  31 + * \f$pivot\f$ is the pivot energy, and
  32 + * \f$index\f$ is the spectral index.
  33 + * \f$curvature\f$ is the curvature
  34 + ***************************************************************************/
  35 +class GModelSpectralLogParabola : public GModelSpectral {
  36 +
  37 +public:
  38 + // Constructors and destructors
  39 + GModelSpectralLogParabola(void);
  40 + explicit GModelSpectralLogParabola(const double& norm, const double& index, const double &curvature);
  41 + explicit GModelSpectralLogParabola(const GXmlElement& xml);
  42 + GModelSpectralLogParabola(const GModelSpectralLogParabola& model);
  43 + virtual ~GModelSpectralLogParabola(void);
  44 +
  45 + // Operators
  46 + virtual GModelSpectralLogParabola& operator=(const GModelSpectralLogParabola& model);
  47 +
  48 + // Implemented pure virtual methods
  49 + virtual void clear(void);
  50 + virtual GModelSpectralLogParabola* clone(void) const;
  51 + virtual std::string type(void) const { return "LogParabola"; }
  52 + virtual double eval(const GEnergy& srcEng) const;
  53 + virtual double eval_gradients(const GEnergy& srcEng) const;
  54 + virtual double flux(const GEnergy& emin, const GEnergy& emax) const;
  55 + virtual double eflux(const GEnergy& emin, const GEnergy& emax) const;
  56 + virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
  57 + virtual void read(const GXmlElement& xml);
  58 + virtual void write(GXmlElement& xml) const;
  59 + virtual std::string print(void) const;
  60 +
  61 + // Other methods
  62 + void autoscale(void);
  63 + double norm(void) const { return m_norm.real_value(); }
  64 + double index(void) const { return m_index.real_value(); }
  65 + double curvature(void) const { return m_curvature.real_value(); }
  66 + double pivot(void) const { return m_pivot.real_value(); }
  67 +
  68 +protected:
  69 + // Protected methods
  70 + void init_members(void);
  71 + void copy_members(const GModelSpectralLogParabola& model);
  72 + void free_members(void);
  73 +
  74 + //class to determine to the integral photon flux
  75 + class flux_kern: public GIntegrand{
  76 + public:
  77 + //Constructor
  78 + flux_kern(double norm, double index, double curvature, double pivot):
  79 + m_norm(norm),m_index(index),m_curvature(curvature),m_pivot(pivot){};
  80 +
  81 + //Method
  82 + double eval(double x){return m_norm*std::pow(x/m_pivot,m_index+m_curvature*std::log(x/m_pivot));}
  83 +
  84 + //Members
  85 + protected:
  86 + double m_norm; //!< Normalization
  87 + double m_index; //!< Spectral index at pivot
  88 + double m_curvature; //!< Curvature
  89 + double m_pivot; //!< Pivot energy
  90 + };
  91 +
  92 + // class to determine the integral energyflux, derviation of flux_kern
  93 + class eflux_kern: public flux_kern{
  94 + public:
  95 + //Constructor
  96 + eflux_kern(double norm, double index, double curvature, double pivot):flux_kern(norm,index,curvature,pivot){};
  97 +
  98 + //Method
  99 + double eval(double x){return x*flux_kern::eval(x);}
  100 + };
  101 +
  102 +
  103 + // Protected members
  104 + GModelPar m_norm; //!< Normalization factor
  105 + GModelPar m_index; //!< Spectral index
  106 + GModelPar m_curvature; //!< Curvature
  107 + GModelPar m_pivot; //!< Pivot energy
  108 +};
  109 +
  110 +
  111 +
  112 +#endif /* GMODELSPECTRALLOGPARABOLA_HPP_ */
... ...
include/GammaLib.hpp
  View file @d1cd172
... ... @@ -133,9 +133,7 @@
133 133  
134 134 /* __ Observation handling _______________________________________________ */
135 135 #include "GEnergy.hpp"
136   -#include "GTimes.hpp"
137 136 #include "GTime.hpp"
138   -#include "GTimeReference.hpp"
139 137 #include "GCaldb.hpp"
140 138 #include "GObservations.hpp"
141 139 #include "GObservation.hpp"
... ... @@ -183,6 +181,7 @@
183 181 #include "GModelSpectralNodes.hpp"
184 182 #include "GModelSpectralPlaw.hpp"
185 183 #include "GModelSpectralPlaw2.hpp"
  184 +#include "GModelSpectralLogParabola.hpp"
186 185 #include "GModelTemporal.hpp"
187 186 #include "GModelTemporalRegistry.hpp"
188 187 #include "GModelTemporalConst.hpp"
... ...
include/Makefile.am
  View file @d1cd172
1 1 # Define package included headers (to be installed)
2 2 pkginclude_HEADERS = GammaLib.hpp \
3   - GBase.hpp \
4   - GRegistry.hpp \
5 3 GException.hpp \
6 4 GNodeArray.hpp \
7 5 GTools.hpp \
... ... @@ -20,9 +18,7 @@ pkginclude_HEADERS = GammaLib.hpp \
20 18 GWcsSTG.hpp \
21 19 GWcsHPX.hpp \
22 20 GEnergy.hpp \
23   - GTimes.hpp \
24 21 GTime.hpp \
25   - GTimeReference.hpp \
26 22 GCaldb.hpp \
27 23 GObservations.hpp \
28 24 GObservation.hpp \
... ... @@ -67,6 +63,7 @@ pkginclude_HEADERS = GammaLib.hpp \
67 63 GModelSpectralPlaw.hpp \
68 64 GModelSpectralPlaw2.hpp \
69 65 GModelSpectralExpPlaw.hpp \
  66 + GModelSpectralLogParabola.hpp \
70 67 GModelSpectralFunc.hpp \
71 68 GModelSpectralNodes.hpp \
72 69 GModelSpectralConst.hpp \
... ... @@ -125,11 +122,10 @@ pkginclude_HEADERS = GammaLib.hpp \
125 122 GIntegrand.hpp \
126 123 GDerivative.hpp \
127 124 GFunction.hpp \
128   - GNumerics.hpp \
129 125 GOptimizer.hpp \
130 126 GOptimizerLM.hpp \
131 127 GOptimizerPars.hpp \
132 128 GOptimizerFunction.hpp \
133 129 GTestSuite.hpp \
134 130 GTestSuites.hpp \
135   - GTestCase.hpp
136 131 \ No newline at end of file
  132 + GTestCase.hpp
... ...
pyext/GModelSpectralLogParabola.i 0 → 100644
  View file @d1cd172
  1 +/***************************************************************************
  2 + * GModelSpectralLogParabola.i - LogParabola model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2009-2011 by Jurgen Knodlseder *
  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 GModelSpectralLogParabola.i
  23 + * @brief LogParabola spectral model class Python interface definition
  24 + * @author M. Mayer
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GModelSpectralLogParabola.hpp"
  29 +#include "GTools.hpp"
  30 +%}
  31 +
  32 +
  33 +/***********************************************************************//**
  34 + * @class GModelSpectralLogParabola
  35 + *
  36 + * @brief LogParabola spectral model class
  37 + ***************************************************************************/
  38 +class GModelSpectralLogParabola : public GModelSpectral {
  39 +public:
  40 + // Constructors and destructors
  41 + GModelSpectralLogParabola(void);
  42 + explicit GModelSpectralLogParabola(const double& norm, const double& index,const double& curvature);
  43 + explicit GModelSpectralLogParabola(const GXmlElement& xml);
  44 + GModelSpectralLogParabola(const GModelSpectralLogParabola& model);
  45 + virtual ~GModelSpectralLogParabola(void);
  46 +
  47 + // Implemented pure virtual methods
  48 + virtual void clear(void);
  49 + virtual GModelSpectralLogParabola* clone(void) const;
  50 + virtual std::string type(void) const;
  51 + virtual double eval(const GEnergy& srcEng) const;
  52 + virtual double eval_gradients(const GEnergy& srcEng) const;
  53 + virtual double flux(const GEnergy& emin, const GEnergy& emax) const;
  54 + virtual double eflux(const GEnergy& emin, const GEnergy& emax) const;
  55 + virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
  56 + virtual void read(const GXmlElement& xml);
  57 + virtual void write(GXmlElement& xml) const;
  58 +
  59 + // Other methods
  60 + void autoscale(void);
  61 + double norm(void) const;
  62 + double index(void) const;
  63 + double curvature(void) const;
  64 + double pivot(void) const;
  65 +};
  66 +
  67 +
  68 +/***********************************************************************//**
  69 + * @brief GModelSpectralLogParabola class extension
  70 + ***************************************************************************/
  71 +%extend GModelSpectralLogParabola {
  72 + GModelSpectralLogParabola copy() {
  73 + return (*self);
  74 + }
  75 +};
... ...
pyext/gammalib/model.i
  View file @d1cd172
... ... @@ -27,12 +27,10 @@
27 27 /* __ Include standard typemaps for vectors and strings __________________ */
28 28 %include stl.i
29 29  
30   -/* __ Include interface classes __________________________________________ */
31   -%import(module="gammalib.base") "GBase.i";
32   -%import(module="gammalib.base") "GRegistry.i";
33   -
34 30 /* __ Make sure that exceptions are catched ______________________________ */
35 31 %import(module="gammalib.support") "GException.i";
  32 +
  33 +/* __ Inform about base classes __________________________________________ */
36 34 %import(module="gammalib.opt") "GOptimizerPars.i";
37 35  
38 36 /* __ Model handling _____________________________________________________ */
... ... @@ -64,6 +62,7 @@
64 62 %include "GModelSpectralNodes.i"
65 63 %include "GModelSpectralPlaw.i"
66 64 %include "GModelSpectralPlaw2.i"
  65 +%include "GModelSpectralLogParabola.i"
67 66 %include "GModelTemporal.i"
68 67 %include "GModelTemporalRegistry.i"
69 68 %include "GModelTemporalConst.i"
... ...
src/model/GModelSpectralLogParabola.cpp 0 → 100644
  View file @d1cd172
  1 +/*
  2 + * GModelSpectralLogParabola.cpp
  3 + *
  4 + * Created on: Nov 30, 2012
  5 + * Author: mimayer
  6 + */
  7 +
  8 +
  9 +/***************************************************************************
  10 + * GModelSpectralLogParabola.cpp - Spectral LogParabola model class *
  11 + * ----------------------------------------------------------------------- *
  12 + * copyright (C) 2009-2012 by Juergen Knoedlseder *
  13 + * ----------------------------------------------------------------------- *
  14 + * *
  15 + * This program is free software: you can redistribute it and/or modify *
  16 + * it under the terms of the GNU General Public License as published by *
  17 + * the Free Software Foundation, either version 3 of the License, or *
  18 + * (at your option) any later version. *
  19 + * *
  20 + * This program is distributed in the hope that it will be useful, *
  21 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  22 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  23 + * GNU General Public License for more details. *
  24 + * *
  25 + * You should have received a copy of the GNU General Public License *
  26 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  27 + * *
  28 + ***************************************************************************/
  29 +/**
  30 + * @file GModelSpectralPlaw.cpp
  31 + * @brief LogParabola spectral model class implementation
  32 + * @author M. Mayer
  33 + */
  34 +
  35 +/* __ Includes ___________________________________________________________ */
  36 +#ifdef HAVE_CONFIG_H
  37 +#include <config.h>
  38 +#endif
  39 +#include <cmath>
  40 +#include "GException.hpp"
  41 +#include "GTools.hpp"
  42 +#include "GModelSpectralLogParabola.hpp"
  43 +#include "GModelSpectralRegistry.hpp"
  44 +
  45 +/* __ Constants __________________________________________________________ */
  46 +
  47 +/* __ Globals ____________________________________________________________ */
  48 +const GModelSpectralLogParabola g_spectral_logparabola_seed;
  49 +const GModelSpectralRegistry g_spectral_logparabola_registry(&g_spectral_logparabola_seed);
  50 +
  51 +/* __ Method name definitions ____________________________________________ */
  52 +#define G_FLUX "GModelSpectralLogParabola::flux(GEnergy&, GEnergy&)"
  53 +#define G_EFLUX "GModelSpectralLogParabola::eflux(GEnergy&, GEnergy&)"
  54 +#define G_MC "GModelSpectralLogParabola::mc(GEnergy&, GEnergy&, GRan&)"
  55 +#define G_READ "GModelSpectralLogParabola::read(GXmlElement&)"
  56 +#define G_WRITE "GModelSpectralLogParabola::write(GXmlElement&)"
  57 +
  58 +/* __ Macros _____________________________________________________________ */
  59 +
  60 +/* __ Coding definitions _________________________________________________ */
  61 +
  62 +/* __ Debug definitions __________________________________________________ */
  63 +
  64 +
  65 +/*==========================================================================
  66 + = =
  67 + = Constructors/destructors =
  68 + = =
  69 + ==========================================================================*/
  70 +
  71 +/***********************************************************************//**
  72 + * @brief Void constructor
  73 + ***************************************************************************/
  74 +GModelSpectralLogParabola::GModelSpectralLogParabola(void) : GModelSpectral()
  75 +{
  76 + // Initialise private members for clean destruction
  77 + init_members();
  78 +
  79 + // Return
  80 + return;
  81 +}
  82 +
  83 +
  84 +/***********************************************************************//**
  85 + * @brief Constructor
  86 + *
  87 + * @param[in] norm Power law normalization.
  88 + * @param[in] index Power law index.
  89 + * @param[in] curvature
  90 + *
  91 + * Construct a LogParabola model from a normalization value, a spectral
  92 + * index and a curvature
  93 + ***************************************************************************/
  94 +GModelSpectralLogParabola::GModelSpectralLogParabola(const double& norm,
  95 + const double& index, const double& curvature) : GModelSpectral()
  96 +{
  97 + // Initialise members
  98 + init_members();
  99 +
  100 + // Set parameters
  101 + m_norm.real_value(norm);
  102 + m_index.real_value(index);
  103 + m_curvature.real_value(curvature);
  104 +
  105 + // Return
  106 + return;
  107 +}
  108 +
  109 +
  110 +/***********************************************************************//**
  111 + * @brief XML constructor
  112 + *
  113 + * @param[in] xml XML element.
  114 + *
  115 + * Creates instance of a log parabola spectral model by extracting information
  116 + * from an XML element. See GModelSpectralLogParabola::read() for more information
  117 + * about the expected structure of the XML element.
  118 + ***************************************************************************/
  119 +GModelSpectralLogParabola::GModelSpectralLogParabola(const GXmlElement& xml) : GModelSpectral()
  120 +{
  121 + // Initialise members
  122 + init_members();
  123 +
  124 + // Read information from XML element
  125 + read(xml);
  126 +
  127 + // Return
  128 + return;
  129 +}
  130 +
  131 +
  132 +/***********************************************************************//**
  133 + * @brief Copy constructor
  134 + *
  135 + * @param[in] model LogParabola model.
  136 + ***************************************************************************/
  137 +GModelSpectralLogParabola::GModelSpectralLogParabola(const GModelSpectralLogParabola& model)
  138 + : GModelSpectral(model)
  139 +{
  140 + // Initialise members
  141 + init_members();
  142 +
  143 + // Copy members
  144 + copy_members(model);
  145 +
  146 + // Return
  147 + return;
  148 +}
  149 +
  150 +
  151 +/***********************************************************************//**
  152 + * @brief Destructor
  153 + ***************************************************************************/
  154 +GModelSpectralLogParabola::~GModelSpectralLogParabola(void)
  155 +{
  156 + // Free members
  157 + free_members();
  158 +
  159 + // Return
  160 + return;
  161 +}
  162 +
  163 +
  164 +/*==========================================================================
  165 + = =
  166 + = Operators =
  167 + = =
  168 + ==========================================================================*/
  169 +
  170 +/***********************************************************************//**
  171 + * @brief Assignment operator
  172 + *
  173 + * @param[in] model LogParabola model.
  174 + ***************************************************************************/
  175 +GModelSpectralLogParabola& GModelSpectralLogParabola::operator=(const GModelSpectralLogParabola& model)
  176 +{
  177 + // Execute only if object is not identical
  178 + if (this != &model) {
  179 +
  180 + // Copy base class members
  181 + this->GModelSpectral::operator=(model);
  182 +
  183 + // Free members
  184 + free_members();
  185 +
  186 + // Initialise members
  187 + init_members();
  188 +
  189 + // Copy members
  190 + copy_members(model);
  191 +
  192 + } // endif: object was not identical
  193 +
  194 + // Return
  195 + return *this;
  196 +}
  197 +
  198 +
  199 +/*==========================================================================
  200 + = =
  201 + = Public methods =
  202 + = =
  203 + ==========================================================================*/
  204 +
  205 +/***********************************************************************//**
  206 + * @brief Clear instance
  207 +***************************************************************************/
  208 +void GModelSpectralLogParabola::clear(void)
  209 +{
  210 + // Free class members (base and derived classes, derived class first)
  211 + free_members();
  212 + this->GModelSpectral::free_members();
  213 +
  214 + // Initialise members
  215 + this->GModelSpectral::init_members();
  216 + init_members();
  217 +
  218 + // Return
  219 + return;
  220 +}
  221 +
  222 +
  223 +/***********************************************************************//**
  224 + * @brief Clone instance
  225 +***************************************************************************/
  226 +GModelSpectralLogParabola* GModelSpectralLogParabola::clone(void) const
  227 +{
  228 + return new GModelSpectralLogParabola(*this);
  229 +}
  230 +
  231 +
  232 +/***********************************************************************//**
  233 + * @brief Evaluate function
  234 + *
  235 + * @param[in] srcEng True energy of photon.
  236 + *
  237 + * The logparabola function is defined as
  238 + * \f[I(E)=norm (E/pivot)^{index-curvature*log(E/pivot)}\f]
  239 + * where
  240 + * \f$norm\f$ is the normalization or prefactor,
  241 + * \f$pivot\f$ is the pivot energy, and
  242 + * \f$index\f$ is the spectral index.
  243 + * \f$curvature\f$ is the curvature
  244 + *
  245 + * @todo For the moment the pivot energy is fixed to units of MeV. This may
  246 + * not be ideal and should eventually be improved in the future.
  247 + * Furthermore, the method expects that energy!=0. Otherwise Inf or NaN
  248 + * may result.
  249 + ***************************************************************************/
  250 +double GModelSpectralLogParabola::eval(const GEnergy& srcEng) const
  251 +{
  252 + // Compute function value
  253 + double energy = srcEng.MeV() / pivot();
  254 + double exponent = index() + curvature()*log(energy);
  255 + double power = std::pow(energy,exponent);
  256 + double value = norm() * power;
  257 +
  258 + // Compile option: Check for NaN/Inf
  259 + #if defined(G_NAN_CHECK)
  260 + if (isnotanumber(value) || isinfinite(value)) {
  261 + std::cout << "*** ERROR: GModelSpectralPlaw::eval";
  262 + std::cout << "(srcEng=" << srcEng << "):";
  263 + std::cout << " NaN/Inf encountered";
  264 + std::cout << " (value=" << value;
  265 + std::cout << ", energy=" << energy;
  266 + std::cout << ", index=" << index();
  267 + std::cout << ", curvature=" << curvature();
  268 + std::cout << ", pivot=" << pivot();
  269 + std::cout << ", power=" << power;
  270 + std::cout << ")" << std::endl;
  271 + }
  272 + #endif
  273 +
  274 + // Return
  275 + return value;
  276 +}
  277 +
  278 +
  279 +/***********************************************************************//**
  280 + * @brief Evaluate function and gradients
  281 + *
  282 + * @param[in] srcEng True energy of photon.
  283 + *
  284 + * The logparabol function is defined as
  285 + * \f[I(E)=norm (E/pivot)^{index-curvature\ln(E/pivot)}\f]
  286 + * where
  287 + * \f$norm=n_s n_v\f$ is the normalization or prefactor,
  288 + * \f$pivot=p_s p_v\f$ is the pivot energy, and
  289 + * \f$index=i_s i_v\f$ is the spectral index.
  290 + * \f$curvature=i_s i_v\f$ is the curvature
  291 + * Note that each parameter is factorised into a scaling factor and a value
  292 + * and that the method is expected to return the gradient with respect to
  293 + * the parameter value (i.e. n_v, p_v, i_v, and c_v in this case).
  294 + *
  295 + * The partial derivatives of the parameter values are given by
  296 + *@todo implement formulae here
  297 + *
  298 + * @todo For the moment the pivot energy is fixed to units of MeV. This may
  299 + * not be ideal and should eventually be improved in the futur.
  300 + * Furthermore, the method expects that energy!=0. Otherwise Inf or NaN
  301 + * may result.
  302 + ***************************************************************************/
  303 +double GModelSpectralLogParabola::eval_gradients(const GEnergy& srcEng) const
  304 +{
  305 + // Compute function value
  306 + double energy = srcEng.MeV() / pivot();
  307 + double exponent = index() + curvature()*std::log(energy);
  308 + double power = std::pow(energy,exponent);
  309 + double value = norm() * power;
  310 +
  311 + // Compute partial derivatives of the parameter values
  312 + double g_norm = (m_norm.isfree()) ? m_norm.scale() * power : 0.0;
  313 + double g_index = (m_index.isfree()) ? value * m_index.scale() * std::log(energy) : 0.0;
  314 + double g_curvature = (m_curvature.isfree()) ? value * m_curvature.scale() * std::pow(std::log(energy),2) : 0.0;
  315 + double g_pivot = (m_pivot.isfree()) ? -value/m_pivot.value() * (exponent + curvature()*std::log(energy)) : 0.0;
  316 +
  317 + // Set gradients (circumvent const correctness)
  318 + const_cast<GModelSpectralLogParabola*>(this)->m_norm.gradient(g_norm);
  319 + const_cast<GModelSpectralLogParabola*>(this)->m_index.gradient(g_index);
  320 + const_cast<GModelSpectralLogParabola*>(this)->m_curvature.gradient(g_curvature);
  321 + const_cast<GModelSpectralLogParabola*>(this)->m_pivot.gradient(g_pivot);
  322 +
  323 + // Compile option: Check for NaN/Inf
  324 + #if defined(G_NAN_CHECK)
  325 + if (isnotanumber(value) || isinfinite(value)) {
  326 + std::cout << "*** ERROR: GModelSpectralPlaw::eval_gradients";
  327 + std::cout << "(srcEng=" << srcEng << "):";
  328 + std::cout << " NaN/Inf encountered";
  329 + std::cout << " (value=" << value;
  330 + std::cout << ", energy=" << energy;
  331 + std::cout << ", index=" << index();
  332 + std::cout << ", curvature=" << curvature();
  333 + std::cout << ", pivot=" << pivot();
  334 + std::cout << ", power=" << power;
  335 + std::cout << ", g_norm=" << g_norm;
  336 + std::cout << ", g_index=" << g_index;
  337 + std::cout << ", g_curvature=" << g_curvature;
  338 + std::cout << ", g_pivot=" << g_pivot;
  339 + std::cout << ")" << std::endl;
  340 + }
  341 + #endif
  342 +
  343 + // Return
  344 + return value;
  345 +}
  346 +
  347 +
  348 +/***********************************************************************//**
  349 + * @brief Returns model photon flux between [emin, emax] (units: ph/cm2/s)
  350 + *
  351 + * @param[in] emin Minimum photon energy.
  352 + * @param[in] emax Maximum photon energy.
  353 + *
  354 + * Computes
  355 + * \f[\int_{E_{\rm min}}^{E_{\rm max}} I(E) dE\f]
  356 + * where
  357 + * \f$E_{\rm min}\f$ and \f$E_{\rm max}\f$ are the minimum and maximum
  358 + * energy, respectively, and
  359 + * \f$I(E)\f$ is the spectral model (units: ph/cm2/s/MeV).
  360 + * The integration is done numerically.
  361 + ***************************************************************************/
  362 +double GModelSpectralLogParabola::flux(const GEnergy& emin, const GEnergy& emax) const
  363 +{
  364 + //Initialise function to integrate
  365 + flux_kern flux(norm(),index(),curvature(),pivot());
  366 +
  367 + //Initialise integral class with function
  368 + GIntegral integral(&flux);
  369 +
  370 + //Set integration precision
  371 + integral.eps(1.0e-8);
  372 +
  373 + //Calculate integral between emin and emax
  374 + double photonflux = integral.romb(emin.MeV(), emax.MeV());
  375 +
  376 + //Return value
  377 + return photonflux;
  378 +}
  379 +
  380 +
  381 +/***********************************************************************//**
  382 + * @brief Returns model energy flux between [emin, emax] (units: erg/cm2/s)
  383 + *
  384 + * @param[in] emin Minimum photon energy.
  385 + * @param[in] emax Maximum photon energy.
  386 + *
  387 + * Computes
  388 + * \f[\int_{E_{\rm min}}^{E_{\rm max}} I(E) E dE\f]
  389 + * where
  390 + * \f$E_{\rm min}\f$ and \f$E_{\rm max}\f$ are the minimum and maximum
  391 + * energy, respectively, and
  392 + * \f$I(E)\f$ is the spectral model (units: ph/cm2/s/MeV).
  393 + * The integration is done numerically.
  394 + ***************************************************************************/
  395 +double GModelSpectralLogParabola::eflux(const GEnergy& emin, const GEnergy& emax) const
  396 +{
  397 + //Initialise function to integrate
  398 + eflux_kern eflux(norm(),index(),curvature(),pivot());
  399 +
  400 + //Initialise integral class with function
  401 + GIntegral integral(&eflux);
  402 +
  403 + //Set integration precision
  404 + integral.eps(1.0e-8);
  405 +
  406 + //Calculate integral between emin and emax
  407 + double energyflux = integral.romb(emin.MeV(), emax.MeV());
  408 +
  409 + //Return value
  410 + return energyflux;
  411 +}
  412 +
  413 +
  414 +/***********************************************************************//**
  415 + * @brief Returns Monte Carlo energy between [emin, emax]
  416 + *
  417 + * @param[in] emin Minimum photon energy.
  418 + * @param[in] emax Maximum photon energy.
  419 + * @param[in] ran Random number generator.
  420 + *
  421 + * Returns Monte Carlo energy by randomly drawing from a power law.
  422 + ***************************************************************************/
  423 +GEnergy GModelSpectralLogParabola::mc(const GEnergy& emin, const GEnergy& emax,
  424 + GRan& ran) const
  425 +{
  426 + // Allocate energy
  427 + GEnergy energy;
  428 +
  429 + // Dump warning that method is not yet implemented
  430 + throw GException::feature_not_implemented(G_MC);
  431 +
  432 + // Return energy
  433 + return energy;
  434 +}
  435 +
  436 +
  437 +/***********************************************************************//**
  438 + * @brief Autoscale normalization
  439 + *
  440 + * Based on the actual value of the m_norm parameter, set the scale of m_norm
  441 + * so that the value will be 1. If minimum and/or maximum value boundaries
  442 + * exist, the boundaries are also modified accordingly.
  443 + ***************************************************************************/
  444 +void GModelSpectralLogParabola::autoscale(void)
  445 +{
  446 + // Autoscale normalization to a value of 1.0
  447 + if (m_norm.value() != 0.0) {
  448 +
  449 + // Get inverse scaling factor
  450 + double invscale = 1.0 / m_norm.value();
  451 +
  452 + // Set values, error, min and max
  453 + m_norm.value(m_norm.value() * invscale);
  454 + m_norm.error(m_norm.error() * invscale);
  455 + if (m_norm.hasmin()) {
  456 + m_norm.min(m_norm.min() * invscale);
  457 + }
  458 + if (m_norm.hasmax()) {
  459 + m_norm.max(m_norm.max() * invscale);
  460 + }
  461 +
  462 + // Set scale
  463 + m_norm.scale(1.0 / invscale);
  464 +
  465 + }
  466 +
  467 + // Return
  468 + return;
  469 +}
  470 +
  471 +
  472 +/***********************************************************************//**
  473 + * @brief Read model from XML element
  474 + *
  475 + * @param[in] xml XML element containing logparabola model information.
  476 + *
  477 + * @exception GException::model_invalid_parnum
  478 + * Invalid number of model parameters found in XML element.
  479 + * @exception GException::model_invalid_parnames
  480 + * Invalid model parameter names found in XML element.
  481 + *
  482 + * Read the logparabola information from an XML element. The XML
  483 + * element is required to have 4 parameters with names "Prefactor", "Index",
  484 + * "Curvature" and "Scale".
  485 + *
  486 + * @todo Add parameter validity check
  487 + ***************************************************************************/
  488 +void GModelSpectralLogParabola::read(const GXmlElement& xml)
  489 +{
  490 + // Verify that XML element has exactly 3 parameters
  491 + if (xml.elements() != 4 || xml.elements("parameter") != 4) {
  492 + throw GException::model_invalid_parnum(G_READ, xml,
  493 + "LogParabola model requires exactly 4 parameters.");
  494 + }
  495 +
  496 + // Extract model parameters
  497 + int npar[] = {0, 0, 0,0};
  498 + for (int i = 0; i < 4; ++i) {
  499 +
  500 + // Get parameter element
  501 + GXmlElement* par = static_cast<GXmlElement*>(xml.element("parameter", i));
  502 +
  503 + // Handle prefactor
  504 + if (par->attribute("name") == "Prefactor" ||
  505 + par->attribute("name")=="norm") {
  506 + m_norm.read(*par);
  507 + npar[0]++;
  508 + }
  509 +
  510 + // Handle index
  511 + else if (par->attribute("name") == "Index" ||
  512 + par->attribute("name")=="alpha") {
  513 + m_index.read(*par);
  514 + npar[1]++;
  515 + }
  516 +
  517 + // Handle index
  518 + else if (par->attribute("name") == "Curvature" ||
  519 + par->attribute("name")=="beta") {
  520 + m_curvature.read(*par);
  521 + npar[2]++;
  522 + }
  523 +
  524 + // Handle pivot energy
  525 + else if (par->attribute("name") == "Scale" ||
  526 + par->attribute("name") == "Eb") {
  527 + m_pivot.read(*par);
  528 + npar[3]++;
  529 + }
  530 +
  531 + } // endfor: looped over all parameters
  532 +
  533 + // Verify that all parameters were found
  534 + if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1 || npar[3] != 1) {
  535 + throw GException::model_invalid_parnames(G_READ, xml,
  536 + "LogParabola requires \"Prefactor\" or \"norm\", \"Index\" or \"alpha\",\"Curvature\" or \"beta\" and \"Scale\" or \"Eb\" parameters.");
  537 + }
  538 +
  539 + // Return
  540 + return;
  541 +}
  542 +
  543 +
  544 +/***********************************************************************//**
  545 + * @brief Write model into XML element
  546 + *
  547 + * @param[in] xml XML element into which model information is written.
  548 + *
  549 + * @exception GException::model_invalid_spectral
  550 + * Existing XML element is not of type "LogParabola"
  551 + * @exception GException::model_invalid_parnum
  552 + * Invalid number of model parameters or nodes found in XML element.
  553 + * @exception GException::model_invalid_parnames
  554 + * Invalid model parameter names found in XML element.
  555 + *
  556 + * Write the LogParabola model information into an XML element. The XML
  557 + * element has to be of type "LogParabola" and will have 4 parameter leafs
  558 + * named "Prefactor", "Index", "Curvature" and "Scale".
  559 + ***************************************************************************/
  560 +void GModelSpectralLogParabola::write(GXmlElement& xml) const
  561 +{
  562 + // Set model type
  563 + if (xml.attribute("type") == "") {
  564 + xml.attribute("type", "LogParabola");
  565 + }
  566 +
  567 + // Verify model type
  568 + if (xml.attribute("type") != "LogParabola") {
  569 + throw GException::model_invalid_spectral(G_WRITE, xml.attribute("type"),
  570 + "Spectral model is not of type \"LogParabola\".");
  571 + }
  572 +
  573 + // If XML element has 0 nodes then append 3 parameter nodes
  574 + if (xml.elements() == 0) {
  575 + xml.append(new GXmlElement("parameter name=\"Prefactor\""));
  576 + xml.append(new GXmlElement("parameter name=\"Index\""));
  577 + xml.append(new GXmlElement("parameter name=\"Curvature\""));
  578 + xml.append(new GXmlElement("parameter name=\"Scale\""));
  579 + }
  580 +
  581 + // Verify that XML element has exactly 3 parameters
  582 + if (xml.elements() != 4 || xml.elements("parameter") != 4) {
  583 + throw GException::model_invalid_parnum(G_WRITE, xml,
  584 + "LogParabola law model requires exactly 4 parameters.");
  585 + }
  586 +
  587 + // Set or update model parameter attributes
  588 + int npar[] = {0, 0, 0,0};
  589 + for (int i = 0; i < 4; ++i) {
  590 +
  591 + // Get parameter element
  592 + GXmlElement* par = static_cast<GXmlElement*>(xml.element("parameter", i));
  593 +
  594 + // Handle prefactor
  595 + if (par->attribute("name") == "Prefactor"){
  596 + npar[0]++;
  597 + m_norm.write(*par);
  598 + }
  599 +
  600 + // Handle index
  601 + else if (par->attribute("name") == "Index") {
  602 + npar[1]++;
  603 + m_index.write(*par);
  604 + }
  605 +
  606 + // Handle index
  607 + else if (par->attribute("name") == "Curvature"){
  608 + npar[2]++;
  609 + m_curvature.write(*par);
  610 + }
  611 +
  612 + // Handle pivot energy
  613 + else if (par->attribute("name") == "Scale") {
  614 + m_pivot.write(*par);
  615 + npar[3]++;
  616 + }
  617 +
  618 + } // endfor: looped over all parameters
  619 +
  620 + // Check of all required parameters are present
  621 + if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1 || npar[3] != 1) {
  622 + throw GException::model_invalid_parnames(G_WRITE, xml,
  623 + "Require \"Prefactor\", \"Index\" , \"Curvature\" and \"Scale\" parameters.");
  624 + }
  625 +
  626 + // Return
  627 + return;
  628 +}
  629 +
  630 +
  631 +/***********************************************************************//**
  632 + * @brief Print powerlaw information
  633 + ***************************************************************************/
  634 +std::string GModelSpectralLogParabola::print(void) const
  635 +{
  636 + // Initialise result string
  637 + std::string result;
  638 +
  639 + // Append header
  640 + result.append("=== GModelSpectralLogParabola ===\n");
  641 + result.append(parformat("Number of parameters")+str(size()));
  642 + for (int i = 0; i < size(); ++i) {
  643 + result.append("\n"+m_pars[i]->print());
  644 + }
  645 +
  646 + // Return result
  647 + return result;
  648 +}
  649 +
  650 +
  651 +/*==========================================================================
  652 + = =
  653 + = Private methods =
  654 + = =
  655 + ==========================================================================*/
  656 +
  657 +/***********************************************************************//**
  658 + * @brief Initialise class members
  659 + ***************************************************************************/
  660 +void GModelSpectralLogParabola::init_members(void)
  661 +{
  662 + // Initialise powerlaw normalisation
  663 + m_norm.clear();
  664 + m_norm.name("Prefactor");
  665 + m_norm.unit("ph/cm2/s/MeV");
  666 + m_norm.scale(1.0);
  667 + m_norm.value(1.0); // default: 1.0
  668 + m_norm.min(0.0); // min: 0.0
  669 + m_norm.free();
  670 + m_norm.gradient(0.0);
  671 + m_norm.hasgrad(true);
  672 +
  673 + // Initialise powerlaw index
  674 + m_index.clear();
  675 + m_index.name("Index");
  676 + m_index.scale(1.0);
  677 + m_index.value(-2.0); // default: -2.0
  678 + m_index.range(-10.0,+10.0); // range: [-10,+10]
  679 + m_index.free();
  680 + m_index.gradient(0.0);
  681 + m_index.hasgrad(true);
  682 +
  683 + m_curvature.clear();
  684 + m_curvature.name("Curvature");
  685 + m_curvature.scale(1.0);
  686 + m_curvature.value(-.1); // default: -2.0
  687 + m_curvature.range(-10.0,+10.0); // range: [-10,+10]
  688 + m_curvature.free();
  689 + m_curvature.gradient(0.0);
  690 + m_curvature.hasgrad(true);
  691 +
  692 + // Initialise pivot energy
  693 + m_pivot.clear();
  694 + m_pivot.name("PivotEnergy");
  695 + m_pivot.unit("MeV");
  696 + m_pivot.scale(1.0);
  697 + m_pivot.value(100.0); // default: 100
  698 + m_pivot.fix();
  699 + m_pivot.gradient(0.0);
  700 + m_pivot.hasgrad(true);
  701 +
  702 + // Set parameter pointer(s)
  703 + m_pars.clear();
  704 + m_pars.push_back(&m_norm);
  705 + m_pars.push_back(&m_index);
  706 + m_pars.push_back(&m_curvature);
  707 + m_pars.push_back(&m_pivot);
  708 +
  709 + // Return
  710 + return;
  711 +}
  712 +
  713 +
  714 +/***********************************************************************//**
  715 + * @brief Copy class members
  716 + *
  717 + * @param[in] model GModelSpectralPlaw members which should be copied.
  718 + ***************************************************************************/
  719 +void GModelSpectralLogParabola::copy_members(const GModelSpectralLogParabola& model)
  720 +{
  721 + // Copy members
  722 + m_norm = model.m_norm;
  723 + m_index = model.m_index;
  724 + m_curvature = model.m_curvature;
  725 + m_pivot = model.m_pivot;
  726 +
  727 + // Set parameter pointer(s)
  728 + m_pars.clear();
  729 + m_pars.push_back(&m_norm);
  730 + m_pars.push_back(&m_index);
  731 + m_pars.push_back(&m_curvature);
  732 + m_pars.push_back(&m_pivot);
  733 +
  734 + // Return
  735 + return;
  736 +}
  737 +
  738 +
  739 +/***********************************************************************//**
  740 + * @brief Delete class members
  741 + ***************************************************************************/
  742 +void GModelSpectralLogParabola::free_members(void)
  743 +{
  744 + // Return
  745 + return;
  746 +}
  747 +
  748 +
  749 +/*==========================================================================
  750 + = =
  751 + = Friends =
  752 + = =
  753 + ==========================================================================*/
  754 +
  755 +
... ...
src/model/Makefile.am
  View file @d1cd172
... ... @@ -34,6 +34,7 @@ sources = GModelPar.cpp \
34 34 GModelSpectralNodes.cpp \
35 35 GModelSpectralPlaw.cpp \
36 36 GModelSpectralPlaw2.cpp \
  37 + GModelSpectralLogParabola.cpp \
37 38 GModelTemporal.cpp \
38 39 GModelTemporalRegistry.cpp \
39 40 GModelTemporalConst.cpp \
... ...