From ea45cec799ff2491bcfff8eed6320f04255ef815 Mon Sep 17 00:00:00 2001 From: Alexander Ziegler Date: Mon, 25 Jul 2016 15:40:02 +0200 Subject: [PATCH] Code clean up after review. See change list below. --- doc/source/users/user_manual/modules/model.rst | 18 +++++++++--------- include/GModelSpectralExpInvPlaw.hpp | 359 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ include/GModelSpectralExpPlaw2.hpp | 321 --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- include/GammaLib.hpp | 2 +- include/Makefile.am | 2 +- pyext/GModelSpectralExpInvPlaw.i | 100 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ pyext/GModelSpectralExpPlaw2.i | 93 --------------------------------------------------------------------------------------------- pyext/gammalib/model.i | 2 +- src/model/GModelSpectralExpInvPlaw.cpp | 976 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/model/GModelSpectralExpPlaw2.cpp | 1014 ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ src/model/Makefile.am | 2 +- test/data/model_point_einvplaw.xml | 15 +++++++++++++++ test/data/model_point_eplaw2.xml | 15 --------------- test/test_GModel.cpp | 48 ++++++++++++++++++++++++++++++++---------------- test/test_GModel.hpp | 4 ++-- 15 files changed, 1497 insertions(+), 1474 deletions(-) create mode 100644 include/GModelSpectralExpInvPlaw.hpp delete mode 100644 include/GModelSpectralExpPlaw2.hpp create mode 100644 pyext/GModelSpectralExpInvPlaw.i delete mode 100644 pyext/GModelSpectralExpPlaw2.i create mode 100644 src/model/GModelSpectralExpInvPlaw.cpp delete mode 100644 src/model/GModelSpectralExpPlaw2.cpp create mode 100644 test/data/model_point_einvplaw.xml delete mode 100644 test/data/model_point_eplaw2.xml diff --git a/doc/source/users/user_manual/modules/model.rst b/doc/source/users/user_manual/modules/model.rst index 9118e68..bf4c47b 100644 --- a/doc/source/users/user_manual/modules/model.rst +++ b/doc/source/users/user_manual/modules/model.rst @@ -424,30 +424,30 @@ The XML format for specifying an exponentially cut-off power law is: An alternative exponentially cut-off power law function is defined by the -:doxy:`GModelSpectralExpPlaw2` class which makes use of the inverse of the +:doxy:`GModelSpectralExpInvPlaw` class which makes use of the inverse of the cut-off energy for function parametrisation: .. math:: \frac{dN}{dE} = k_0 \left( \frac{E}{E_0} \right)^{\gamma} - \exp \left( -\lambda E \right) + \exp \left( -\alpha E \right) where the parameters in the XML definition have the following mappings: * :math:`k_0` = ``Prefactor`` * :math:`\gamma` = ``Index`` -* :math:`E_0` = ``Scale`` -* :math:`\lambda` = ``lambda`` +* :math:`E_0` = ``PivotEnergy`` +* :math:`\alpha` = ``InverseCutoff`` The XML format for specifying an exponentially cut-off power law using this alternative parametrisation is: .. code-block:: xml - - - - - + + + + + diff --git a/include/GModelSpectralExpInvPlaw.hpp b/include/GModelSpectralExpInvPlaw.hpp new file mode 100644 index 0000000..1cc4181 --- /dev/null +++ b/include/GModelSpectralExpInvPlaw.hpp @@ -0,0 +1,359 @@ +/*************************************************************************** + * GModelSpectralExpInvPlaw.hpp - Exponential cut off power law model * + * ----------------------------------------------------------------------- * + * copyright (C) 2016 by Alexander Ziegler * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GModelSpectralExpInvPlaw.hpp + * @brief Exponential cut off power law spectral class interface definition + * @author Alexander Ziegler + */ + +#ifndef GMODELSPECTRALEXPINVPLAW_HPP +#define GMODELSPECTRALEXPINVPLAW_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include "GModelSpectral.hpp" +#include "GFunction.hpp" +#include "GModelPar.hpp" +#include "GEnergy.hpp" + +/* __ Forward declarations _______________________________________________ */ +class GRan; +class GTime; +class GXmlElement; + + +/***********************************************************************//** + * @class GModelSpectralExpInvPlaw + * + * @brief Exponential cut off power law spectral class + * + * This class implements a power law spectrum with exponential cut off. The + * model is defined by + * + * \f[ + * S_{\rm E}(E | t) = {\tt m\_norm} + * \left( \frac{E}{\tt m\_pivot} \right)^{\tt m\_index} + * \exp \left( -{\tt m\_alpha}*E \right) + * \f] + * + * where + * - \f${\tt m\_norm}\f$ is the normalization or prefactor, + * - \f${\tt m\_pivot}\f$ is the pivot energy, + * - \f${\tt m\_index}\f$ is the spectral index, and + * - \f${\tt m\_alpha}\f$ is the cut-off parameter (1/cut-off energy). + ***************************************************************************/ +class GModelSpectralExpInvPlaw : public GModelSpectral { + +public: + // Constructors and destructors + GModelSpectralExpInvPlaw(void); + explicit GModelSpectralExpInvPlaw(const double& prefactor, + const double& index, + const GEnergy& pivot, + const double& alpha); + explicit GModelSpectralExpInvPlaw(const double& prefactor, + const double& index, + const GEnergy& pivot, + const GEnergy& cutoff); + explicit GModelSpectralExpInvPlaw(const GXmlElement& xml); + GModelSpectralExpInvPlaw(const GModelSpectralExpInvPlaw& model); + virtual ~GModelSpectralExpInvPlaw(void); + + // Operators + virtual GModelSpectralExpInvPlaw& operator=(const GModelSpectralExpInvPlaw& model); + + // Implemented pure virtual methods + virtual void clear(void); + virtual GModelSpectralExpInvPlaw* clone(void) const; + virtual std::string classname(void) const; + virtual std::string type(void) const; + virtual double eval(const GEnergy& srcEng, + const GTime& srcTime = GTime()) const; + virtual double eval_gradients(const GEnergy& srcEng, + const GTime& srcTime = GTime()); + virtual double flux(const GEnergy& emin, + const GEnergy& emax) const; + virtual double eflux(const GEnergy& emin, + const GEnergy& emax) const; + virtual GEnergy mc(const GEnergy& emin, + const GEnergy& emax, + const GTime& time, + GRan& ran) const; + virtual void read(const GXmlElement& xml); + virtual void write(GXmlElement& xml) const; + virtual std::string print(const GChatter& chatter = NORMAL) const; + + // Other methods + double prefactor(void) const; + void prefactor(const double& prefactor); + double index(void) const; + void index(const double& index); + double inverse_cutoff(void) const; + void inverse_cutoff(const double& alpha); + GEnergy cutoff(void) const; + void cutoff(const GEnergy& cutoff); + GEnergy pivot(void) const; + void pivot(const GEnergy& pivot); + +protected: + // Protected methods + void init_members(void); + void copy_members(const GModelSpectralExpInvPlaw& model); + void free_members(void); + void update_eval_cache(const GEnergy& energy) const; + void update_mc_cache(const GEnergy& emin, const GEnergy& emax) const; + + // Photon flux integration kernel + class flux_kernel : public GFunction { + public: + flux_kernel(const double& norm, + const double& index, + const double& pivot, + const double& alpha) : + m_norm(norm), + m_index(index), + m_inv_pivot(1.0/pivot), + m_alpha(alpha) {} + double eval(const double& eng); + protected: + double m_norm; //!< Normalization + double m_index; //!< Index + double m_inv_pivot; //!< 1 / Pivot energy + double m_alpha; //!< Cut-off parameter + }; + + // Energy flux integration kernel + class eflux_kernel : public GFunction { + public: + eflux_kernel(const double& norm, + const double& index, + const double& pivot, + const double& alpha) : + m_norm(norm), + m_index(index), + m_inv_pivot(1.0/pivot), + m_alpha(alpha) {} + double eval(const double& eng); + protected: + double m_norm; //!< Normalization + double m_index; //!< Index + double m_inv_pivot; //!< 1 / Pivot energy + double m_alpha; //!< Cut-off parameter + }; + + // Protected members + GModelPar m_norm; //!< Normalization factor + GModelPar m_index; //!< Spectral index + GModelPar m_alpha; //!< Cut-off parameter + GModelPar m_pivot; //!< Pivot energy + + // Cached members used for pre-computations + mutable GEnergy m_last_energy; //!< Last energy value + mutable double m_last_index; //!< Last index parameter + mutable double m_last_alpha; //!< Last cut-off parameter + mutable double m_last_pivot; //!< Last pivot parameter + mutable double m_last_e_norm; //!< Last E/Epivot value + mutable double m_last_e_alpha; //!< Last E*alpha value + mutable double m_last_power; //!< Last power value + mutable double m_mc_emin; //!< Minimum energy + mutable double m_mc_emax; //!< Maximum energy + mutable double m_mc_exponent; //!< Exponent (index+1) + mutable double m_mc_pow_emin; //!< Power of minimum energy + mutable double m_mc_pow_ewidth; //!< Power of energy width +}; + + +/***********************************************************************//** + * @brief Return class name + * + * @return String containing the class name ("GModelSpectralExpInvPlaw"). + ***************************************************************************/ +inline +std::string GModelSpectralExpInvPlaw::classname(void) const +{ + return ("GModelSpectralExpInvPlaw"); +} + + +/***********************************************************************//** + * @brief Return model type + * + * @return "ExponentialInverseCutoffPowerLaw". + * + * Returns the type of the exponentially cut off power law model. + ***************************************************************************/ +inline +std::string GModelSpectralExpInvPlaw::type(void) const +{ + return "ExponentialInverseCutoffPowerLaw"; +} + + +/****ExponentialInverseCutoffPowerLaw*********************************************************//** + * @brief Return pre factor + * + * @return Pre factor (ph/cm2/s/MeV). + * + * Returns the pre factor. + ***************************************************************************/ +inline +double GModelSpectralExpInvPlaw::prefactor(void) const +{ + return (m_norm.value()); +} + + +/***********************************************************************//** + * @brief Set pre factor + * + * @param[in] prefactor Pre factor (ph/cm2/s/MeV). + * + * Sets the pre factor. + ***************************************************************************/ +inline +void GModelSpectralExpInvPlaw::prefactor(const double& prefactor) +{ + m_norm.value(prefactor); + return; +} + + +/***********************************************************************//** + * @brief Return power law index + * + * @return Power law index. + * + * Returns the power law index. + ***************************************************************************/ +inline +double GModelSpectralExpInvPlaw::index(void) const +{ + return (m_index.value()); +} + + +/***********************************************************************//** + * @brief Set power law index + * + * @param[in] index Power law index. + * + * Sets the power law index. + ***************************************************************************/ +inline +void GModelSpectralExpInvPlaw::index(const double& index) +{ + m_index.value(index); + return; +} + + +/***********************************************************************//** + * @brief Return pivot energy + * + * @return Pivot energy. + * + * Returns the pivot energy. + ***************************************************************************/ +inline +GEnergy GModelSpectralExpInvPlaw::pivot(void) const +{ + GEnergy energy; + energy.MeV(m_pivot.value()); + return energy; +} + + +/***********************************************************************//** + * @brief Set pivot energy + * + * @param[in] pivot Pivot energy. + * + * Sets the pivot energy. + ***************************************************************************/ +inline +void GModelSpectralExpInvPlaw::pivot(const GEnergy& pivot) +{ + m_pivot.value(pivot.MeV()); + return; +} + + +/***********************************************************************//** + * @brief Return exponential cut-off parameter + * + * @return Exponential cut-off parameter. + * + * Returns the exponential cut-off parameter. + ***************************************************************************/ +inline +double GModelSpectralExpInvPlaw::inverse_cutoff(void) const +{ + return (m_alpha.value()); +} + + +/***********************************************************************//** + * @brief Set exponential cut-off parameter + * + * @param[in] alpha Exponential cut-off parameter. + * + * Sets the exponential cut-off parameter. + ***************************************************************************/ +inline +void GModelSpectralExpInvPlaw::inverse_cutoff(const double& alpha) +{ + m_alpha.value(alpha); + return; +} + + +/***********************************************************************//** + * @brief Return cut-off energy (1/cut-off parameter) + * + * @return Exponential cut-off energy. + * + * Returns the exponential cut-off energy. + ***************************************************************************/ +inline +GEnergy GModelSpectralExpInvPlaw::cutoff(void) const +{ + GEnergy cutoff; + cutoff.MeV(1.0/m_alpha.value()); + return cutoff; +} + + +/***********************************************************************//** + * @brief Set exponential cut-off parameter via passing cut-off energy + * + * @param[in] cutoff Exponential cut-off energy. + * + * Sets the exponential cut-off parameter. + ***************************************************************************/ +inline +void GModelSpectralExpInvPlaw::cutoff(const GEnergy& cutoff) +{ + m_alpha.value(1./cutoff.MeV()); + return; +} + + +#endif /* GMODELSPECTRALEXPINVPLAW_HPP */ diff --git a/include/GModelSpectralExpPlaw2.hpp b/include/GModelSpectralExpPlaw2.hpp deleted file mode 100644 index f705f31..0000000 --- a/include/GModelSpectralExpPlaw2.hpp +++ /dev/null @@ -1,321 +0,0 @@ -/*************************************************************************** - * GModelSpectralExpPlaw2.hpp - Exponential cut off power law model * - * ----------------------------------------------------------------------- * - * copyright (C) 2016 by Alexander Ziegler * - * ----------------------------------------------------------------------- * - * * - * This program is free software: you can redistribute it and/or modify * - * it under the terms of the GNU General Public License as published by * - * the Free Software Foundation, either version 3 of the License, or * - * (at your option) any later version. * - * * - * This program is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * - * GNU General Public License for more details. * - * * - * You should have received a copy of the GNU General Public License * - * along with this program. If not, see . * - * * - ***************************************************************************/ -/** - * @file GModelSpectralExpPlaw2.hpp - * @brief Exponential cut off power law spectral class interface definition - * @author Alexander Ziegler - */ - -#ifndef GMODELSPECTRALEXPPLAW2_HPP -#define GMODELSPECTRALEXPPLAW2_HPP - -/* __ Includes ___________________________________________________________ */ -#include -#include "GModelSpectral.hpp" -#include "GFunction.hpp" -#include "GModelPar.hpp" -#include "GEnergy.hpp" - -/* __ Forward declarations _______________________________________________ */ -class GRan; -class GTime; -class GXmlElement; - - -/***********************************************************************//** - * @class GModelSpectralExpPlaw2 - * - * @brief Exponential cut off power law spectral class - * - * This class implements a power law spectrum with exponential cut off. The - * model is defined by - * - * \f[ - * S_{\rm E}(E | t) = {\tt m\_norm} - * \left( \frac{E}{\tt m\_pivot} \right)^{\tt m\_index} - * \exp \left( -{\tt m\_lamda}*E \right) - * \f] - * - * where - * - \f${\tt m\_norm}\f$ is the normalization or prefactor, - * - \f${\tt m\_pivot}\f$ is the pivot energy, - * - \f${\tt m\_index}\f$ is the spectral index, and - * - \f${\tt m\_lambda}\f$ is the cut-off parameter (1/cut-off energy). - ***************************************************************************/ -class GModelSpectralExpPlaw2 : public GModelSpectral { - -public: - // Constructors and destructors - GModelSpectralExpPlaw2(void); - explicit GModelSpectralExpPlaw2(const double& prefactor, - const double& index, - const GEnergy& pivot, - const double& lambda); - explicit GModelSpectralExpPlaw2(const GXmlElement& xml); - GModelSpectralExpPlaw2(const GModelSpectralExpPlaw2& model); - virtual ~GModelSpectralExpPlaw2(void); - - // Operators - virtual GModelSpectralExpPlaw2& operator=(const GModelSpectralExpPlaw2& model); - - // Implemented pure virtual methods - virtual void clear(void); - virtual GModelSpectralExpPlaw2* clone(void) const; - virtual std::string classname(void) const; - virtual std::string type(void) const; - virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); - virtual double flux(const GEnergy& emin, - const GEnergy& emax) const; - virtual double eflux(const GEnergy& emin, - const GEnergy& emax) const; - virtual GEnergy mc(const GEnergy& emin, - const GEnergy& emax, - const GTime& time, - GRan& ran) const; - virtual void read(const GXmlElement& xml); - virtual void write(GXmlElement& xml) const; - virtual std::string print(const GChatter& chatter = NORMAL) const; - - // Other methods - double prefactor(void) const; - void prefactor(const double& prefactor); - double index(void) const; - void index(const double& index); - double lambda(void) const; - void lambda(const double& lambda); - GEnergy pivot(void) const; - void pivot(const GEnergy& pivot); - -protected: - // Protected methods - void init_members(void); - void copy_members(const GModelSpectralExpPlaw2& model); - void free_members(void); - void update_eval_cache(const GEnergy& energy) const; - void update_mc_cache(const GEnergy& emin, const GEnergy& emax) const; - - // Photon flux integration kernel - class flux_kernel : public GFunction { - public: - flux_kernel(const double& norm, - const double& index, - const double& pivot, - const double& lambda) : - m_norm(norm), - m_index(index), - m_inv_pivot(1.0/pivot), - m_lambda(lambda) {} - double eval(const double& eng); - protected: - double m_norm; //!< Normalization - double m_index; //!< Index - double m_inv_pivot; //!< 1 / Pivot energy - double m_lambda; //!< Cut-off parameter - }; - - // Energy flux integration kernel - class eflux_kernel : public GFunction { - public: - eflux_kernel(const double& norm, - const double& index, - const double& pivot, - const double& lambda) : - m_norm(norm), - m_index(index), - m_inv_pivot(1.0/pivot), - m_lambda(lambda) {} - double eval(const double& eng); - protected: - double m_norm; //!< Normalization - double m_index; //!< Index - double m_inv_pivot; //!< 1 / Pivot energy - double m_lambda; //!< Cut-off parameter - }; - - // Protected members - GModelPar m_norm; //!< Normalization factor - GModelPar m_index; //!< Spectral index - GModelPar m_lambda; //!< Cut-off parameter - GModelPar m_pivot; //!< Pivot energy - - // Cached members used for pre-computations - mutable GEnergy m_last_energy; //!< Last energy value - mutable double m_last_index; //!< Last index parameter - mutable double m_last_lambda; //!< Last cut-off parameter - mutable double m_last_pivot; //!< Last pivot parameter - mutable double m_last_e_norm; //!< Last E/Epivot value - mutable double m_last_e_lambda; //!< Last E*lambda value - mutable double m_last_power; //!< Last power value - mutable double m_mc_emin; //!< Minimum energy - mutable double m_mc_emax; //!< Maximum energy - mutable double m_mc_exponent; //!< Exponent (index+1) - mutable double m_mc_pow_emin; //!< Power of minimum energy - mutable double m_mc_pow_ewidth; //!< Power of energy width -}; - - -/***********************************************************************//** - * @brief Return class name - * - * @return String containing the class name ("GModelSpectralExpPlaw2"). - ***************************************************************************/ -inline -std::string GModelSpectralExpPlaw2::classname(void) const -{ - return ("GModelSpectralExpPlaw2"); -} - - -/***********************************************************************//** - * @brief Return model type - * - * @return "ExpCutoff2". - * - * Returns the type of the exponentially cut off power law model. - ***************************************************************************/ -inline -std::string GModelSpectralExpPlaw2::type(void) const -{ - return "ExpCutoff2"; -} - - -/***********************************************************************//** - * @brief Return pre factor - * - * @return Pre factor (ph/cm2/s/MeV). - * - * Returns the pre factor. - ***************************************************************************/ -inline -double GModelSpectralExpPlaw2::prefactor(void) const -{ - return (m_norm.value()); -} - - -/***********************************************************************//** - * @brief Set pre factor - * - * @param[in] prefactor Pre factor (ph/cm2/s/MeV). - * - * Sets the pre factor. - ***************************************************************************/ -inline -void GModelSpectralExpPlaw2::prefactor(const double& prefactor) -{ - m_norm.value(prefactor); - return; -} - - -/***********************************************************************//** - * @brief Return power law index - * - * @return Power law index. - * - * Returns the power law index. - ***************************************************************************/ -inline -double GModelSpectralExpPlaw2::index(void) const -{ - return (m_index.value()); -} - - -/***********************************************************************//** - * @brief Set power law index - * - * @param[in] index Power law index. - * - * Sets the power law index. - ***************************************************************************/ -inline -void GModelSpectralExpPlaw2::index(const double& index) -{ - m_index.value(index); - return; -} - - -/***********************************************************************//** - * @brief Return pivot energy - * - * @return Pivot energy. - * - * Returns the pivot energy. - ***************************************************************************/ -inline -GEnergy GModelSpectralExpPlaw2::pivot(void) const -{ - GEnergy energy; - energy.MeV(m_pivot.value()); - return energy; -} - - -/***********************************************************************//** - * @brief Set pivot energy - * - * @param[in] pivot Pivot energy. - * - * Sets the pivot energy. - ***************************************************************************/ -inline -void GModelSpectralExpPlaw2::pivot(const GEnergy& pivot) -{ - m_pivot.value(pivot.MeV()); - return; -} - - -/***********************************************************************//** - * @brief Return exponential cut-off parameter - * - * @return Exponential cut-off parameter. - * - * Returns the exponential cut-off parameter. - ***************************************************************************/ -inline -double GModelSpectralExpPlaw2::lambda(void) const -{ - return (m_lambda.value()); -} - - -/***********************************************************************//** - * @brief Set exponential cut-off parameter - * - * @param[in] lambda Exponential cut-off parameter. - * - * Sets the exponential cut-off parameter. - ***************************************************************************/ -inline -void GModelSpectralExpPlaw2::lambda(const double& lambda) -{ - m_lambda.value(lambda); - return; -} - -#endif /* GMODELSPECTRALEXPPLAW2_HPP */ diff --git a/include/GammaLib.hpp b/include/GammaLib.hpp index 4d275dc..9783598 100644 --- a/include/GammaLib.hpp +++ b/include/GammaLib.hpp @@ -209,7 +209,7 @@ #include "GModelSpectralBrokenPlaw.hpp" #include "GModelSpectralConst.hpp" #include "GModelSpectralExpPlaw.hpp" -#include "GModelSpectralExpPlaw2.hpp" +#include "GModelSpectralExpInvPlaw.hpp" #include "GModelSpectralSuperExpPlaw.hpp" #include "GModelSpectralFunc.hpp" #include "GModelSpectralGauss.hpp" diff --git a/include/Makefile.am b/include/Makefile.am index fd26f1a..08a59f8 100644 --- a/include/Makefile.am +++ b/include/Makefile.am @@ -145,8 +145,8 @@ pkginclude_HEADERS = GammaLib.hpp \ GModelSpectralRegistry.hpp \ GModelSpectralBrokenPlaw.hpp \ GModelSpectralConst.hpp \ + GModelSpectralExpInvPlaw.hpp \ GModelSpectralExpPlaw.hpp \ - GModelSpectralExpPlaw2.hpp \ GModelSpectralSuperExpPlaw.hpp \ GModelSpectralFunc.hpp \ GModelSpectralGauss.hpp \ diff --git a/pyext/GModelSpectralExpInvPlaw.i b/pyext/GModelSpectralExpInvPlaw.i new file mode 100644 index 0000000..cdd19b4 --- /dev/null +++ b/pyext/GModelSpectralExpInvPlaw.i @@ -0,0 +1,100 @@ +/*************************************************************************** + * GModelSpectralExpInvPlaw.i - Exponential cut off power law model * + * ----------------------------------------------------------------------- * + * copyright (C) 2016 by Alexander Ziegler * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GModelSpectralExpInvPlaw.i + * @brief Exponential cut off power law spectral class interface definition + * @author Alexander Ziegler + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GModelSpectralExpInvPlaw.hpp" +%} + + +/***********************************************************************//** + * @class GModelSpectralExpInvPlaw + * + * @brief Exponential cut off power law spectral class + ***************************************************************************/ +class GModelSpectralExpInvPlaw : public GModelSpectral { + +public: + // Constructors and destructors + GModelSpectralExpInvPlaw(void); + explicit GModelSpectralExpInvPlaw(const double& prefactor, + const double& index, + const GEnergy& pivot, + const double& alpha); + explicit GModelSpectralExpInvPlaw(const double& prefactor, + const double& index, + const GEnergy& pivot, + const GEnergy& cutoff); + explicit GModelSpectralExpInvPlaw(const GXmlElement& xml); + GModelSpectralExpInvPlaw(const GModelSpectralExpInvPlaw& model); + virtual ~GModelSpectralExpInvPlaw(void); + + // Operators + virtual GModelSpectralExpInvPlaw& operator=(const GModelSpectralExpInvPlaw& model); + + // Implemented pure virtual methods + virtual void clear(void); + virtual GModelSpectralExpInvPlaw* clone(void) const; + virtual std::string classname(void) const; + virtual std::string type(void) const; + virtual double eval(const GEnergy& srcEng, + const GTime& srcTime = GTime()) const; + virtual double eval_gradients(const GEnergy& srcEng, + const GTime& srcTime = GTime()); + virtual double flux(const GEnergy& emin, + const GEnergy& emax) const; + virtual double eflux(const GEnergy& emin, + const GEnergy& emax) const; + virtual GEnergy mc(const GEnergy& emin, + const GEnergy& emax, + const GTime& time, + GRan& ran) const; + virtual void read(const GXmlElement& xml); + virtual void write(GXmlElement& xml) const; + virtual std::string print(const GChatter& chatter = NORMAL) const; + + // Other methods + double prefactor(void) const; + void prefactor(const double& prefactor); + double index(void) const; + void index(const double& index); + double inverse_cutoff(void) const; + void inverse_cutoff(const double& alpha); + GEnergy cutoff(void) const; + void cutoff(const GEnergy& cutoff); + GEnergy pivot(void) const; + void pivot(const GEnergy& pivot); + +}; + + +/***********************************************************************//** + * @brief GModelSpectralExpInvPlaw class extension + ***************************************************************************/ +%extend GModelSpectralExpInvPlaw { + GModelSpectralExpInvPlaw copy() { + return (*self); + } +}; diff --git a/pyext/GModelSpectralExpPlaw2.i b/pyext/GModelSpectralExpPlaw2.i deleted file mode 100644 index 8c34bc9..0000000 --- a/pyext/GModelSpectralExpPlaw2.i +++ /dev/null @@ -1,93 +0,0 @@ -/*************************************************************************** - * GModelSpectralExpPlaw2.i - Exponential cut off power law model * - * ----------------------------------------------------------------------- * - * copyright (C) 2016 by Alexander Ziegler * - * ----------------------------------------------------------------------- * - * * - * This program is free software: you can redistribute it and/or modify * - * it under the terms of the GNU General Public License as published by * - * the Free Software Foundation, either version 3 of the License, or * - * (at your option) any later version. * - * * - * This program is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * - * GNU General Public License for more details. * - * * - * You should have received a copy of the GNU General Public License * - * along with this program. If not, see . * - * * - ***************************************************************************/ -/** - * @file GModelSpectralExpPlaw2.i - * @brief Exponential cut off power law spectral class interface definition - * @author Alexander Ziegler - */ -%{ -/* Put headers and other declarations here that are needed for compilation */ -#include "GModelSpectralExpPlaw2.hpp" -%} - - -/***********************************************************************//** - * @class GModelSpectralExpPlaw2 - * - * @brief Exponential cut off power law spectral class - ***************************************************************************/ -class GModelSpectralExpPlaw2 : public GModelSpectral { - -public: - // Constructors and destructors - GModelSpectralExpPlaw2(void); - explicit GModelSpectralExpPlaw2(const double& prefactor, - const double& index, - const GEnergy& pivot, - const double& lambda); - explicit GModelSpectralExpPlaw2(const GXmlElement& xml); - GModelSpectralExpPlaw2(const GModelSpectralExpPlaw2& model); - virtual ~GModelSpectralExpPlaw2(void); - - // Operators - virtual GModelSpectralExpPlaw2& operator=(const GModelSpectralExpPlaw2& model); - - // Implemented pure virtual methods - virtual void clear(void); - virtual GModelSpectralExpPlaw2* clone(void) const; - virtual std::string classname(void) const; - virtual std::string type(void) const; - virtual double eval(const GEnergy& srcEng, - const GTime& srcTime = GTime()) const; - virtual double eval_gradients(const GEnergy& srcEng, - const GTime& srcTime = GTime()); - virtual double flux(const GEnergy& emin, - const GEnergy& emax) const; - virtual double eflux(const GEnergy& emin, - const GEnergy& emax) const; - virtual GEnergy mc(const GEnergy& emin, - const GEnergy& emax, - const GTime& time, - GRan& ran) const; - virtual void read(const GXmlElement& xml); - virtual void write(GXmlElement& xml) const; - virtual std::string print(const GChatter& chatter = NORMAL) const; - - // Other methods - double prefactor(void) const; - void prefactor(const double& prefactor); - double index(void) const; - void index(const double& index); - double lambda(void) const; - void lambda(const double& lambda); - GEnergy pivot(void) const; - void pivot(const GEnergy& pivot); -}; - - -/***********************************************************************//** - * @brief GModelSpectralExpPlaw2 class extension - ***************************************************************************/ -%extend GModelSpectralExpPlaw2 { - GModelSpectralExpPlaw2 copy() { - return (*self); - } -}; diff --git a/pyext/gammalib/model.i b/pyext/gammalib/model.i index 381fe3f..c8f93d5 100644 --- a/pyext/gammalib/model.i +++ b/pyext/gammalib/model.i @@ -145,8 +145,8 @@ %include "GModelSpectralRegistry.i" %include "GModelSpectralBrokenPlaw.i" %include "GModelSpectralConst.i" +%include "GModelSpectralExpInvPlaw.i" %include "GModelSpectralExpPlaw.i" -%include "GModelSpectralExpPlaw2.i" %include "GModelSpectralSuperExpPlaw.i" %include "GModelSpectralFunc.i" %include "GModelSpectralGauss.i" diff --git a/src/model/GModelSpectralExpInvPlaw.cpp b/src/model/GModelSpectralExpInvPlaw.cpp new file mode 100644 index 0000000..d8747ad --- /dev/null +++ b/src/model/GModelSpectralExpInvPlaw.cpp @@ -0,0 +1,976 @@ +/*************************************************************************** + * GModelSpectralExpInvPlaw.cpp - Exponential cut off power law model * + * ----------------------------------------------------------------------- * + * copyright (C) 2016 by Alexander Ziegler * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GModelSpectralExpInvPlaw.cpp + * @brief Exponential cut off power law spectral class implementation + * @author Alexander Ziegler + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "GException.hpp" +#include "GTools.hpp" +#include "GRan.hpp" +#include "GIntegral.hpp" +#include "GModelSpectralExpInvPlaw.hpp" +#include "GModelSpectralRegistry.hpp" + +/* __ Constants __________________________________________________________ */ + +/* __ Globals ____________________________________________________________ */ +const GModelSpectralExpInvPlaw g_spectral_expinvplaw_seed; +const GModelSpectralRegistry g_spectral_expinvplaw_registry(&g_spectral_expinvplaw_seed); + +/* __ Method name definitions ____________________________________________ */ +#define G_MC "GModelSpectralExpInvPlaw::mc(GEnergy&, GEnergy&, GTime&,"\ + " GRan&)" +#define G_READ "GModelSpectralExpInvPlaw::read(GXmlElement&)" +#define G_WRITE "GModelSpectralExpInvPlaw::write(GXmlElement&)" + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + +/* __ Debug definitions __________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +GModelSpectralExpInvPlaw::GModelSpectralExpInvPlaw(void) : GModelSpectral() +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Constructor + * + * @param[in] prefactor Pre factor normalization (ph/cm2/s/MeV). + * @param[in] index Power law index. + * @param[in] pivot Pivot energy. + * @param[in] alpha Cut-off parameter. + * + * Construct an exponentially cut off power law from + * - a prefactor value (in units of ph/cm2/s/MeV) + * - a spectral index, + * - a pivot energy, and + * - a cut-off parameter alpha (in units 1/MeV). + ***************************************************************************/ +GModelSpectralExpInvPlaw::GModelSpectralExpInvPlaw(const double& prefactor, + const double& index, + const GEnergy& pivot, + const double& alpha) : + GModelSpectral() +{ + // Initialise members + init_members(); + + // Set parameters + m_norm.value(prefactor); + m_index.value(index); + m_pivot.value(pivot.MeV()); // Internally stored in MeV + m_alpha.value(alpha); // Internally stored in 1/MeV + + // Autoscale parameters + autoscale(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Constructor (construct model via cut-off energy) + * + * @param[in] prefactor Pre factor normalization (ph/cm2/s/MeV). + * @param[in] index Power law index. + * @param[in] pivot Pivot energy. + * @param[in] cutoff Cut-off energy. + * + * Construct an exponentially cut off power law from + * - a prefactor value (in units of ph/cm2/s/MeV) + * - a spectral index, + * - a pivot energy, and + * - a cut-off energy (in unit MeV). + ***************************************************************************/ +GModelSpectralExpInvPlaw::GModelSpectralExpInvPlaw(const double& prefactor, + const double& index, + const GEnergy& pivot, + const GEnergy& cutoff) : + GModelSpectral() +{ + // Initialise members + init_members(); + + // Set parameters + m_norm.value(prefactor); + m_index.value(index); + m_pivot.value(pivot.MeV()); // Internally stored in MeV + m_alpha.value(1./cutoff.MeV()); // Internally stored in 1/MeV + + // Autoscale parameters + autoscale(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief XML constructor + * + * @param[in] xml XML element. + * + * Constructs an exponentially cut off power law spectral model by extracting + * information from an XML element. See the read() method for more + * information about the expected structure of the XML element. + ***************************************************************************/ +GModelSpectralExpInvPlaw::GModelSpectralExpInvPlaw(const GXmlElement& xml) : + GModelSpectral() +{ + // Initialise members + init_members(); + + // Read information from XML element + read(xml); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] model Exponentially cut off power law model. + ***************************************************************************/ +GModelSpectralExpInvPlaw::GModelSpectralExpInvPlaw(const GModelSpectralExpInvPlaw& model) : + GModelSpectral(model) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(model); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +GModelSpectralExpInvPlaw::~GModelSpectralExpInvPlaw(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] model Exponentially cut off power law model. + * @return Exponentially cut off power law model. + ***************************************************************************/ +GModelSpectralExpInvPlaw& GModelSpectralExpInvPlaw::operator=(const GModelSpectralExpInvPlaw& model) +{ + // Execute only if object is not identical + if (this != &model) { + + // Copy base class members + this->GModelSpectral::operator=(model); + + // Free members + free_members(); + + // Initialise members + init_members(); + + // Copy members + copy_members(model); + + } // endif: object was not identical + + // Return + return *this; +} + + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Clear exponentially cut off power law model + ***************************************************************************/ +void GModelSpectralExpInvPlaw::clear(void) +{ + // Free class members (base and derived classes, derived class first) + free_members(); + this->GModelSpectral::free_members(); + + // Initialise members + this->GModelSpectral::init_members(); + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Clone exponentially cut off power law model + * + * @return Pointer to deep copy of exponentially cut off power law model. + ***************************************************************************/ +GModelSpectralExpInvPlaw* GModelSpectralExpInvPlaw::clone(void) const +{ + // Clone exponentially cut off power law model + return new GModelSpectralExpInvPlaw(*this); +} + + +/***********************************************************************//** + * @brief Evaluate function + * + * @param[in] srcEng True photon energy. + * @param[in] srcTime True photon arrival time. + * @return Model value. + * + * Evaluates + * + * \f[ + * S_{\rm E}(E | t) = {\tt m\_norm} + * \left( \frac{E}{\tt m\_pivot} \right)^{\tt m\_index} + * \exp \left( -{\tt m\_alpha}*E \right) + * \f] + * + * where + * - \f${\tt m\_norm}\f$ is the normalization or prefactor, + * - \f${\tt m\_pivot}\f$ is the pivot energy, + * - \f${\tt m\_index}\f$ is the spectral index, and + * - \f${\tt m\_alpha}\f$ is the cut-off parameter. + * + * @todo The method expects that pivot!=0. Otherwise Inf or NaN + * may result. We should add a test that prevents using invalid + * values. + ***************************************************************************/ +double GModelSpectralExpInvPlaw::eval(const GEnergy& srcEng, + const GTime& srcTime) const +{ + // Update the evaluation cache + update_eval_cache(srcEng); + + // Compute function value + double value = m_norm.value() * m_last_power; + + // Compile option: Check for NaN/Inf + #if defined(G_NAN_CHECK) + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { + std::cout << "*** ERROR: GModelSpectralExpPlaw::eval"; + std::cout << "(srcEng=" << srcEng; + std::cout << ", srcTime=" << srcTime << "):"; + std::cout << " NaN/Inf encountered"; + std::cout << " (value=" << value; + std::cout << ", power=" << m_last_power; + std::cout << ")" << std::endl; + } + #endif + + // Return + return value; +} + + +/***********************************************************************//** + * @brief Evaluate function and gradients + * + * @param[in] srcEng True photon energy. + * @param[in] srcTime True photon arrival time. + * @return Model value. + * + * Evaluates + * + * \f[ + * S_{\rm E}(E | t) = {\tt m\_norm} + * \left( \frac{E}{\tt m\_pivot} \right)^{\tt m\_index} + * \exp \left( -{\tt m\_alpha}*E \right) + * \f] + * + * where + * - \f${\tt m\_norm}\f$ is the normalization or prefactor, + * - \f${\tt m\_pivot}\f$ is the pivot energy, + * - \f${\tt m\_index}\f$ is the spectral index, and + * - \f${\tt m\_alpha}\f$ is the cut-off parameter. + * + * The method also evaluates the partial derivatives of the model with + * respect to the parameters using + * + * \f[ + * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_norm}} = + * \frac{S_{\rm E}(E | t)}{{\tt m\_norm}} + * \f] + * + * \f[ + * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_index}} = + * S_{\rm E}(E | t) \, \ln(E/{\tt m_pivot}) + * \f] + * + * \f[ + * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_alpha}} = + * S_{\rm E}(E | t) \, \left( -E \right) + * \f] + * + * \f[ + * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_pivot}} = + * -S_{\rm E}(E | t) \, + * \left( \frac{{\tt m\_index}}{{\tt m\_pivot}} \right) + * \f] + * + * @todo The method expects that pivot!=0. Otherwise Inf or NaN + * may result. We should add a test that prevents using invalid + * values. + ***************************************************************************/ +double GModelSpectralExpInvPlaw::eval_gradients(const GEnergy& srcEng, + const GTime& srcTime) +{ + // Update the evaluation cache + update_eval_cache(srcEng); + + // Compute function value + double value = m_norm.value() * m_last_power; + + // Compute partial derivatives with respect to the parameter factor + // values. The partial derivatives with respect to the parameter + // values are obtained by division by the scale factor. + double g_norm = (m_norm.is_free()) + ? m_norm.scale() * m_last_power : 0.0; + double g_index = (m_index.is_free()) + ? value * m_index.scale() * std::log(m_last_e_norm) : 0.0; + double g_alpha = (m_alpha.is_free()) + ? -value * m_alpha.scale() * m_last_energy.MeV() : 0.0; + double g_pivot = (m_pivot.is_free()) + ? -value * m_last_index / m_pivot.factor_value() : 0.0; + + // Set gradients + m_norm.factor_gradient(g_norm); + m_index.factor_gradient(g_index); + m_alpha.factor_gradient(g_alpha); + m_pivot.factor_gradient(g_pivot); + + // Compile option: Check for NaN/Inf + #if defined(G_NAN_CHECK) + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { + std::cout << "*** ERROR: GModelSpectralExpPlaw::eval_gradients"; + std::cout << "(srcEng=" << srcEng; + std::cout << ", srcTime=" << srcTime << "):"; + std::cout << " NaN/Inf encountered"; + std::cout << " (value=" << value; + std::cout << ", e_norm=" << m_last_e_norm; + std::cout << ", e_alpha=" << m_last_e_alpha; + std::cout << ", power=" << m_last_power; + std::cout << ")" << std::endl; + } + #endif + + // Return + return value; +} + + +/***********************************************************************//** + * @brief Returns model photon flux between [emin, emax] (units: ph/cm2/s) + * + * @param[in] emin Minimum photon energy. + * @param[in] emax Maximum photon energy. + * @return Photon flux (ph/cm2/s). + * + * Computes + * + * \f[ + * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) dE + * \f] + * + * where + * - [@p emin, @p emax] is an energy interval, and + * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV). + * The integration is done numerically. + ***************************************************************************/ +double GModelSpectralExpInvPlaw::flux(const GEnergy& emin, + const GEnergy& emax) const +{ + // Initialise flux + double flux = 0.0; + + // Compute only if integration range is valid + if (emin < emax) { + + // Setup integration kernel + flux_kernel integrand(m_norm.value(), m_index.value(), + m_pivot.value(), m_alpha.value()); + GIntegral integral(&integrand); + + // Get integration boundaries in MeV + double e_min = emin.MeV(); + double e_max = emax.MeV(); + + // Perform integration + flux = integral.romberg(e_min, e_max); + + } // endif: integration range was valid + + // Return + return flux; +} + + +/***********************************************************************//** + * @brief Returns model energy flux between [emin, emax] (units: erg/cm2/s) + * + * @param[in] emin Minimum photon energy. + * @param[in] emax Maximum photon energy. + * @return Energy flux (erg/cm2/s). + * + * Computes + * + * \f[ + * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) E \, dE + * \f] + * + * where + * - [@p emin, @p emax] is an energy interval, and + * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV). + * The integration is done numerically. + ***************************************************************************/ +double GModelSpectralExpInvPlaw::eflux(const GEnergy& emin, + const GEnergy& emax) const +{ + // Initialise flux + double eflux = 0.0; + + // Compute only if integration range is valid + if (emin < emax) { + + // Setup integration kernel + eflux_kernel integrand(m_norm.value(), m_index.value(), + m_pivot.value(), m_alpha.value()); + GIntegral integral(&integrand); + + // Get integration boundaries in MeV + double e_min = emin.MeV(); + double e_max = emax.MeV(); + + // Perform integration + eflux = integral.romberg(e_min, e_max); + + // Convert from MeV/cm2/s to erg/cm2/s + eflux *= gammalib::MeV2erg; + + } // endif: integration range was valid + + // Return + return eflux; +} + + +/***********************************************************************//** + * @brief Returns MC energy between [emin, emax] + * + * @param[in] emin Minimum photon energy. + * @param[in] emax Maximum photon energy. + * @param[in] time True photon arrival time. + * @param[in,out] ran Random number generator. + * @return Energy. + * + * @exception GException::erange_invalid + * Energy range is invalid (emin < emax required). + * + * Simulates a random energy in the interval [emin, emax] for an + * exponentially cut off power law. The simulation is done using a rejection + * method. First, a random energy within [emin, emax] is drawn from an + * power law distribution. Then the energy is accepted or rejected based + * on an acceptance fraction that is computed from the exponential cut off. + ***************************************************************************/ +GEnergy GModelSpectralExpInvPlaw::mc(const GEnergy& emin, + const GEnergy& emax, + const GTime& time, + GRan& ran) const +{ + // Throw an exception if energy range is invalid + if (emin >= emax) { + throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(), + "Minimum energy < maximum energy required."); + } + + // Allocate energy + GEnergy energy; + + // Update cache + update_mc_cache(emin, emax); + + // Initialise energy + double eng; + + // Initialse acceptance fraction + double acceptance_fraction; + double inv_ecut = m_alpha.value(); + + // Use rejection method to draw a random energy. We first draw + // analytically from a power law, and then compare the power law + // at the drawn energy to the exponentially cut off function. This + // gives an acceptance fraction, and we accept the energy only if + // a uniform random number is <= the acceptance fraction. + do { + + // Get uniform random number + double u = ran.uniform(); + + // Case A: Index is not -1 + if (index() != -1.0) { + if (u > 0.0) { + eng = std::exp(std::log(u * m_mc_pow_ewidth + m_mc_pow_emin) / + m_mc_exponent); + } + else { + eng = 0.0; + } + } + + // Case B: Index is -1 + else { + eng = std::exp(u * m_mc_pow_ewidth + m_mc_pow_emin); + } + + // Compute acceptance fraction + acceptance_fraction = std::exp(-eng * inv_ecut); + + } while (ran.uniform() > acceptance_fraction); + + // Set energy + energy.MeV(eng); + + // Return energy + return energy; +} + + +/***********************************************************************//** + * @brief Read model from XML element + * + * @param[in] xml XML element containing power law model information. + * + * @exception GException::model_invalid_parnum + * Invalid number of model parameters found in XML element. + * + * Reads the spectral information from an XML element. The format of the XML + * elements is + * + * + * + * + * + * + * + ***************************************************************************/ +void GModelSpectralExpInvPlaw::read(const GXmlElement& xml) +{ + // Verify that XML element has exactly 4 parameters + if (xml.elements() != 4 || xml.elements("parameter") != 4) { + throw GException::model_invalid_parnum(G_READ, xml, + "ExponentialInverseCutoffPowerLaw model requires exactly 4 parameters."); + } + + // Get XML parameters + const GXmlElement* prefactor = gammalib::xml_get_par(G_READ, xml, "Prefactor"); + const GXmlElement* index = gammalib::xml_get_par(G_READ, xml, "Index"); + const GXmlElement* alpha = gammalib::xml_get_par(G_READ, xml, "InverseCutoff"); + const GXmlElement* pivot = gammalib::xml_get_par(G_READ, xml, "PivotEnergy"); + + // Read parameters + m_norm.read(*prefactor); + m_index.read(*index); + m_alpha.read(*alpha); + m_pivot.read(*pivot); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Write model into XML element + * + * @param[in] xml XML element into which model information is written. + * + * @exception GException::model_invalid_spectral + * Existing XML element is not of type "ExponentialInverseCutoffPowerLaw" + * + * Writes the spectral information into an XML element. The format of the XML + * element is + * + * + * + * + * + * + * + ***************************************************************************/ +void GModelSpectralExpInvPlaw::write(GXmlElement& xml) const +{ + // Set model type + if (xml.attribute("type") == "") { + xml.attribute("type", type()); + } + + // Verify model type + if (xml.attribute("type") != type()) { + throw GException::model_invalid_spectral(G_WRITE, xml.attribute("type"), + "Spectral model is not of type \""+type()+"\"."); + } + + // Get XML parameters + GXmlElement* norm = gammalib::xml_need_par(G_WRITE, xml, m_norm.name()); + GXmlElement* index = gammalib::xml_need_par(G_WRITE, xml, m_index.name()); + GXmlElement* alpha = gammalib::xml_need_par(G_WRITE, xml, m_alpha.name()); + GXmlElement* pivot = gammalib::xml_need_par(G_WRITE, xml, m_pivot.name()); + + // Write parameters + m_norm.write(*norm); + m_index.write(*index); + m_alpha.write(*alpha); + m_pivot.write(*pivot); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Print model information + * + * @param[in] chatter Chattiness (defaults to NORMAL). + * @return String containing model information. + ***************************************************************************/ +std::string GModelSpectralExpInvPlaw::print(const GChatter& chatter) const +{ + // Initialise result string + std::string result; + + // Continue only if chatter is not silent + if (chatter != SILENT) { + + // Append header + result.append("=== GModelSpectralExpInvPlaw ==="); + + // Append information + result.append("\n"+gammalib::parformat("Number of parameters")); + result.append(gammalib::str(size())); + for (int i = 0; i < size(); ++i) { + result.append("\n"+m_pars[i]->print(chatter)); + } + + } // endif: chatter was not silent + + // Return result + return result; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void GModelSpectralExpInvPlaw::init_members(void) +{ + // Initialise powerlaw normalisation + m_norm.clear(); + m_norm.name("Prefactor"); + m_norm.unit("ph/cm2/s/MeV"); + m_norm.scale(1.0); + m_norm.value(1.0); // default: 1.0 + m_norm.min(0.0); // min: 0.0 + m_norm.free(); + m_norm.gradient(0.0); + m_norm.has_grad(true); + + // Initialise powerlaw index + m_index.clear(); + m_index.name("Index"); + m_index.scale(1.0); + m_index.value(-2.0); // default: -2.0 + m_index.range(-10.0,+10.0); // range: [-10,+10] + m_index.free(); + m_index.gradient(0.0); + m_index.has_grad(true); + + // Initialise cut off parameter + m_alpha.clear(); + m_alpha.name("InverseCutoff"); + m_alpha.unit("1/MeV"); + m_alpha.scale(1.0); + m_alpha.value(0.001); // default: 0.001 + m_alpha.min(0.0); // min: 0.0 + m_alpha.max(10.0); // max: 10.0 + m_alpha.free(); + m_alpha.gradient(0.0); + m_alpha.has_grad(true); + + // Initialise pivot energy + m_pivot.clear(); + m_pivot.name("PivotEnergy"); + m_pivot.unit("MeV"); + m_pivot.scale(1.0); + m_pivot.value(100.0); + m_pivot.fix(); + m_pivot.gradient(0.0); + m_pivot.has_grad(true); + + // Set parameter pointer(s) + m_pars.clear(); + m_pars.push_back(&m_norm); + m_pars.push_back(&m_index); + m_pars.push_back(&m_alpha); + m_pars.push_back(&m_pivot); + + // Initialise eval cache + m_last_energy.clear(); + m_last_index = 1.0e30; + m_last_alpha = 1.0e30; + m_last_pivot = 1.0e30; + m_last_e_norm = 0.0; + m_last_e_alpha = 0.0; + m_last_power = 0.0; + + // Initialise MC cache + m_mc_emin = 0.0; + m_mc_emax = 0.0; + m_mc_exponent = 0.0; + m_mc_pow_emin = 0.0; + m_mc_pow_ewidth = 0.0; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] model Exponential cut off power law model. + ***************************************************************************/ +void GModelSpectralExpInvPlaw::copy_members(const GModelSpectralExpInvPlaw& model) +{ + // Copy members + m_norm = model.m_norm; + m_index = model.m_index; + m_alpha = model.m_alpha; + m_pivot = model.m_pivot; + + // Set parameter pointer(s) + m_pars.clear(); + m_pars.push_back(&m_norm); + m_pars.push_back(&m_index); + m_pars.push_back(&m_alpha); + m_pars.push_back(&m_pivot); + + // Copy eval cache + m_last_energy = model.m_last_energy; + m_last_index = model.m_last_index; + m_last_alpha = model.m_last_alpha; + m_last_pivot = model.m_last_pivot; + m_last_e_norm = model.m_last_e_norm; + m_last_e_alpha = model.m_last_e_alpha; + m_last_power = model.m_last_power; + + // Copy MC cache + m_mc_emin = model.m_mc_emin; + m_mc_emax = model.m_mc_emax; + m_mc_exponent = model.m_mc_exponent; + m_mc_pow_emin = model.m_mc_pow_emin; + m_mc_pow_ewidth = model.m_mc_pow_ewidth; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void GModelSpectralExpInvPlaw::free_members(void) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief Update eval precomputation cache + * + * @param[in] energy Energy. + * + * Updates the precomputation cache for eval() and eval_gradients() methods. + ***************************************************************************/ +void GModelSpectralExpInvPlaw::update_eval_cache(const GEnergy& energy) const +{ + // Get parameter values (takes 3 multiplications which are difficult + // to avoid) + double index = m_index.value(); + double alpha = m_alpha.value(); + double pivot = m_pivot.value(); + + // If the energy or one of the parameters index, alpha or pivot + // energy has changed then recompute the cache + if ((m_last_energy != energy) || + (m_last_index != index) || + (m_last_alpha != alpha) || + (m_last_pivot != pivot)) { + + // Store actual energy and parameter values + m_last_energy = energy; + m_last_index = index; + m_last_alpha = alpha; + m_last_pivot = pivot; + + // Compute and store value + double eng = energy.MeV(); + m_last_e_norm = eng / m_last_pivot; + m_last_e_alpha = eng * m_last_alpha; + m_last_power = std::pow(m_last_e_norm, m_last_index) * + std::exp(-m_last_e_alpha); + } // endif: recomputation was required + + // Return + return; +} + + +/***********************************************************************//** + * @brief Update Monte Carlo pre computation cache + * + * @param[in] emin Minimum photon energy. + * @param[in] emax Maximum photon energy. + * + * Updates the precomputation cache for Monte Carlo simulations. + ***************************************************************************/ +void GModelSpectralExpInvPlaw::update_mc_cache(const GEnergy& emin, + const GEnergy& emax) const + +{ + // Case A: Index is not -1 + if (index() != -1.0) { + + // Change in energy boundaries? + if (emin.MeV() != m_mc_emin || emax.MeV() != m_mc_emax) { + m_mc_emin = emin.MeV(); + m_mc_emax = emax.MeV(); + m_mc_exponent = index() + 1.0; + m_mc_pow_emin = std::pow(m_mc_emin, m_mc_exponent); + m_mc_pow_ewidth = std::pow(m_mc_emax, m_mc_exponent) - m_mc_pow_emin; + } + + } + + // Case B: Index is -1 + else { + + // Change in energy boundaries? + if (emin.MeV() != m_mc_emin || emax.MeV() != m_mc_emax) { + m_mc_emin = emin.MeV(); + m_mc_emax = emax.MeV(); + m_mc_exponent = 0.0; + m_mc_pow_emin = std::log(m_mc_emin); + m_mc_pow_ewidth = std::log(m_mc_emax) - m_mc_pow_emin; + } + + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Kernel for photon flux integration + * + * @param[in] energy Energy (MeV). + ***************************************************************************/ +double GModelSpectralExpInvPlaw::flux_kernel::eval(const double& energy) +{ + // Evaluate function value + double e_norm = energy * m_inv_pivot; + double e_cut = energy * m_alpha; + double power = std::pow(e_norm, m_index) * std::exp(-e_cut); + double value = m_norm * power; + + // Return value + return value; +} + + +/***********************************************************************//** + * @brief Kernel for energy flux integration + * + * @param[in] energy Energy (MeV). + ***************************************************************************/ +double GModelSpectralExpInvPlaw::eflux_kernel::eval(const double& energy) +{ + // Evaluate function value + double e_norm = energy * m_inv_pivot; + double e_cut = energy * m_alpha; + double power = std::pow(e_norm, m_index) * std::exp(-e_cut); + double value = m_norm * power * energy; + + // Return value + return value; +} diff --git a/src/model/GModelSpectralExpPlaw2.cpp b/src/model/GModelSpectralExpPlaw2.cpp deleted file mode 100644 index 01328ee..0000000 --- a/src/model/GModelSpectralExpPlaw2.cpp +++ /dev/null @@ -1,1014 +0,0 @@ -/*************************************************************************** - * GModelSpectralExpPlaw2.cpp - Exponential cut off power law model * - * ----------------------------------------------------------------------- * - * copyright (C) 2016 by Alexander Ziegler * - * ----------------------------------------------------------------------- * - * * - * This program is free software: you can redistribute it and/or modify * - * it under the terms of the GNU General Public License as published by * - * the Free Software Foundation, either version 3 of the License, or * - * (at your option) any later version. * - * * - * This program is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * - * GNU General Public License for more details. * - * * - * You should have received a copy of the GNU General Public License * - * along with this program. If not, see . * - * * - ***************************************************************************/ -/** - * @file GModelSpectralExpPlaw2.cpp - * @brief Exponential cut off power law spectral class implementation - * @author Alexander Ziegler - */ - -/* __ Includes ___________________________________________________________ */ -#ifdef HAVE_CONFIG_H -#include -#endif -#include -#include "GException.hpp" -#include "GTools.hpp" -#include "GRan.hpp" -#include "GIntegral.hpp" -#include "GModelSpectralExpPlaw2.hpp" -#include "GModelSpectralRegistry.hpp" - -/* __ Constants __________________________________________________________ */ - -/* __ Globals ____________________________________________________________ */ -const GModelSpectralExpPlaw2 g_spectral_expplaw2_seed; -const GModelSpectralRegistry g_spectral_expplaw2_registry(&g_spectral_expplaw2_seed); - -/* __ Method name definitions ____________________________________________ */ -#define G_MC "GModelSpectralExpPlaw2::mc(GEnergy&, GEnergy&, GTime&, GRan&)" -#define G_READ "GModelSpectralExpPlaw2::read(GXmlElement&)" -#define G_WRITE "GModelSpectralExpPlaw2::write(GXmlElement&)" - -/* __ Macros _____________________________________________________________ */ - -/* __ Coding definitions _________________________________________________ */ - -/* __ Debug definitions __________________________________________________ */ - - -/*========================================================================== - = = - = Constructors/destructors = - = = - ==========================================================================*/ - -/***********************************************************************//** - * @brief Void constructor - ***************************************************************************/ -GModelSpectralExpPlaw2::GModelSpectralExpPlaw2(void) : GModelSpectral() -{ - // Initialise members - init_members(); - - // Return - return; -} - - -/***********************************************************************//** - * @brief Constructor - * - * @param[in] prefactor Pre factor normalization (ph/cm2/s/MeV). - * @param[in] index Power law index. - * @param[in] pivot Pivot energy. - * @param[in] lambda Cut-off parameter. - * - * Construct an exponentially cut off power law from - * - a prefactor value (in units of ph/cm2/s/MeV) - * - a spectral index, - * - a pivot energy, and - * - a cut-off parameter lambda. - ***************************************************************************/ -GModelSpectralExpPlaw2::GModelSpectralExpPlaw2(const double& prefactor, - const double& index, - const GEnergy& pivot, - const double & lambda) : - GModelSpectral() -{ - // Initialise members - init_members(); - - // Set parameters - m_norm.value(prefactor); - m_index.value(index); - m_pivot.value(pivot.MeV()); // Internally stored in MeV - m_lambda.value(lambda); // Internally stored in 1/MeV - - // Autoscale parameters - autoscale(); - - // Return - return; -} - - -/***********************************************************************//** - * @brief XML constructor - * - * @param[in] xml XML element. - * - * Constructs an exponentially cut off power law spectral model by extracting - * information from an XML element. See the read() method for more - * information about the expected structure of the XML element. - ***************************************************************************/ -GModelSpectralExpPlaw2::GModelSpectralExpPlaw2(const GXmlElement& xml) : - GModelSpectral() -{ - // Initialise members - init_members(); - - // Read information from XML element - read(xml); - - // Return - return; -} - - -/***********************************************************************//** - * @brief Copy constructor - * - * @param[in] model Exponentially cut off power law model. - ***************************************************************************/ -GModelSpectralExpPlaw2::GModelSpectralExpPlaw2(const GModelSpectralExpPlaw2& model) : - GModelSpectral(model) -{ - // Initialise members - init_members(); - - // Copy members - copy_members(model); - - // Return - return; -} - - -/***********************************************************************//** - * @brief Destructor - ***************************************************************************/ -GModelSpectralExpPlaw2::~GModelSpectralExpPlaw2(void) -{ - // Free members - free_members(); - - // Return - return; -} - - -/*========================================================================== - = = - = Operators = - = = - ==========================================================================*/ - -/***********************************************************************//** - * @brief Assignment operator - * - * @param[in] model Exponentially cut off power law model. - * @return Exponentially cut off power law model. - ***************************************************************************/ -GModelSpectralExpPlaw2& GModelSpectralExpPlaw2::operator=(const GModelSpectralExpPlaw2& model) -{ - // Execute only if object is not identical - if (this != &model) { - - // Copy base class members - this->GModelSpectral::operator=(model); - - // Free members - free_members(); - - // Initialise members - init_members(); - - // Copy members - copy_members(model); - - } // endif: object was not identical - - // Return - return *this; -} - - -/*========================================================================== - = = - = Public methods = - = = - ==========================================================================*/ - -/***********************************************************************//** - * @brief Clear exponentially cut off power law model - ***************************************************************************/ -void GModelSpectralExpPlaw2::clear(void) -{ - // Free class members (base and derived classes, derived class first) - free_members(); - this->GModelSpectral::free_members(); - - // Initialise members - this->GModelSpectral::init_members(); - init_members(); - - // Return - return; -} - - -/***********************************************************************//** - * @brief Clone exponentially cut off power law model - * - * @return Pointer to deep copy of exponentially cut off power law model. - ***************************************************************************/ -GModelSpectralExpPlaw2* GModelSpectralExpPlaw2::clone(void) const -{ - // Clone exponentially cut off power law model - return new GModelSpectralExpPlaw2(*this); -} - - -/***********************************************************************//** - * @brief Evaluate function - * - * @param[in] srcEng True photon energy. - * @param[in] srcTime True photon arrival time. - * @return Model value. - * - * Evaluates - * - * \f[ - * S_{\rm E}(E | t) = {\tt m\_norm} - * \left( \frac{E}{\tt m\_pivot} \right)^{\tt m\_index} - * \exp \left( -{\tt m\_lamda}*E \right) - * \f] - * - * where - * - \f${\tt m\_norm}\f$ is the normalization or prefactor, - * - \f${\tt m\_pivot}\f$ is the pivot energy, - * - \f${\tt m\_index}\f$ is the spectral index, and - * - \f${\tt m\_lambda}\f$ is the cut-off parameter. - * - * @todo The method expects that pivot!=0. Otherwise Inf or NaN - * may result. We should add a test that prevents using invalid - * values. - ***************************************************************************/ -double GModelSpectralExpPlaw2::eval(const GEnergy& srcEng, - const GTime& srcTime) const -{ - // Update the evaluation cache - update_eval_cache(srcEng); - - // Compute function value - double value = m_norm.value() * m_last_power; - - // Compile option: Check for NaN/Inf - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: GModelSpectralExpPlaw::eval"; - std::cout << "(srcEng=" << srcEng; - std::cout << ", srcTime=" << srcTime << "):"; - std::cout << " NaN/Inf encountered"; - std::cout << " (value=" << value; - std::cout << ", power=" << m_last_power; - std::cout << ")" << std::endl; - } - #endif - - // Return - return value; -} - - -/***********************************************************************//** - * @brief Evaluate function and gradients - * - * @param[in] srcEng True photon energy. - * @param[in] srcTime True photon arrival time. - * @return Model value. - * - * Evaluates - * - * \f[ - * S_{\rm E}(E | t) = {\tt m\_norm} - * \left( \frac{E}{\tt m\_pivot} \right)^{\tt m\_index} - * \exp \left( -{\tt m\_lamda}*E \right) - * \f] - * - * where - * - \f${\tt m\_norm}\f$ is the normalization or prefactor, - * - \f${\tt m\_pivot}\f$ is the pivot energy, - * - \f${\tt m\_index}\f$ is the spectral index, and - * - \f${\tt m\_lambda}\f$ is the cut-off parameter. - * - * The method also evaluates the partial derivatives of the model with - * respect to the parameters using - * - * \f[ - * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_norm}} = - * \frac{S_{\rm E}(E | t)}{{\tt m\_norm}} - * \f] - * - * \f[ - * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_index}} = - * S_{\rm E}(E | t) \, \ln(E/{\tt m_pivot}) - * \f] - * - * \f[ - * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_lambda}} = - * S_{\rm E}(E | t) \, \left( -E \right) - * \f] - * - * \f[ - * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_pivot}} = - * -S_{\rm E}(E | t) \, - * \left( \frac{{\tt m\_index}}{{\tt m\_pivot}} \right) - * \f] - * - * @todo The method expects that pivot!=0. Otherwise Inf or NaN - * may result. We should add a test that prevents using invalid - * values. - ***************************************************************************/ -double GModelSpectralExpPlaw2::eval_gradients(const GEnergy& srcEng, - const GTime& srcTime) -{ - // Update the evaluation cache - update_eval_cache(srcEng); - - // Compute function value - double value = m_norm.value() * m_last_power; - - // Compute partial derivatives with respect to the parameter factor - // values. The partial derivatives with respect to the parameter - // values are obtained by division by the scale factor. - double g_norm = (m_norm.is_free()) - ? m_norm.scale() * m_last_power : 0.0; - double g_index = (m_index.is_free()) - ? value * m_index.scale() * std::log(m_last_e_norm) : 0.0; - double g_lambda = (m_lambda.is_free()) - ? -value * m_lambda.scale() * m_last_energy.MeV() : 0.0; - double g_pivot = (m_pivot.is_free()) - ? -value * m_last_index / m_pivot.factor_value() : 0.0; - - // Set gradients - m_norm.factor_gradient(g_norm); - m_index.factor_gradient(g_index); - m_lambda.factor_gradient(g_lambda); - m_pivot.factor_gradient(g_pivot); - - // Compile option: Check for NaN/Inf - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: GModelSpectralExpPlaw::eval_gradients"; - std::cout << "(srcEng=" << srcEng; - std::cout << ", srcTime=" << srcTime << "):"; - std::cout << " NaN/Inf encountered"; - std::cout << " (value=" << value; - std::cout << ", e_norm=" << m_last_e_norm; - std::cout << ", e_lambda=" << m_last_e_lambda; - std::cout << ", power=" << m_last_power; - std::cout << ")" << std::endl; - } - #endif - - // Return - return value; -} - - -/***********************************************************************//** - * @brief Returns model photon flux between [emin, emax] (units: ph/cm2/s) - * - * @param[in] emin Minimum photon energy. - * @param[in] emax Maximum photon energy. - * @return Photon flux (ph/cm2/s). - * - * Computes - * - * \f[ - * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) dE - * \f] - * - * where - * - [@p emin, @p emax] is an energy interval, and - * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV). - * The integration is done numerically. - ***************************************************************************/ -double GModelSpectralExpPlaw2::flux(const GEnergy& emin, - const GEnergy& emax) const -{ - // Initialise flux - double flux = 0.0; - - // Compute only if integration range is valid - if (emin < emax) { - - // Setup integration kernel - flux_kernel integrand(m_norm.value(), m_index.value(), - m_pivot.value(), m_lambda.value()); - GIntegral integral(&integrand); - - // Get integration boundaries in MeV - double e_min = emin.MeV(); - double e_max = emax.MeV(); - - // Perform integration - flux = integral.romberg(e_min, e_max); - - } // endif: integration range was valid - - // Return - return flux; -} - - -/***********************************************************************//** - * @brief Returns model energy flux between [emin, emax] (units: erg/cm2/s) - * - * @param[in] emin Minimum photon energy. - * @param[in] emax Maximum photon energy. - * @return Energy flux (erg/cm2/s). - * - * Computes - * - * \f[ - * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) E \, dE - * \f] - * - * where - * - [@p emin, @p emax] is an energy interval, and - * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV). - * The integration is done numerically. - ***************************************************************************/ -double GModelSpectralExpPlaw2::eflux(const GEnergy& emin, - const GEnergy& emax) const -{ - // Initialise flux - double eflux = 0.0; - - // Compute only if integration range is valid - if (emin < emax) { - - // Setup integration kernel - eflux_kernel integrand(m_norm.value(), m_index.value(), - m_pivot.value(), m_lambda.value()); - GIntegral integral(&integrand); - - // Get integration boundaries in MeV - double e_min = emin.MeV(); - double e_max = emax.MeV(); - - // Perform integration - eflux = integral.romberg(e_min, e_max); - - // Convert from MeV/cm2/s to erg/cm2/s - eflux *= gammalib::MeV2erg; - - } // endif: integration range was valid - - // Return - return eflux; -} - - -/***********************************************************************//** - * @brief Returns MC energy between [emin, emax] - * - * @param[in] emin Minimum photon energy. - * @param[in] emax Maximum photon energy. - * @param[in] time True photon arrival time. - * @param[in,out] ran Random number generator. - * @return Energy. - * - * @exception GException::erange_invalid - * Energy range is invalid (emin < emax required). - * - * Simulates a random energy in the interval [emin, emax] for an - * exponentially cut off power law. The simulation is done using a rejection - * method. First, a random energy within [emin, emax] is drawn from an - * power law distribution. Then the energy is accepted or rejected based - * on an acceptance fraction that is computed from the exponential cut off. - ***************************************************************************/ -GEnergy GModelSpectralExpPlaw2::mc(const GEnergy& emin, - const GEnergy& emax, - const GTime& time, - GRan& ran) const -{ - // Throw an exception if energy range is invalid - if (emin >= emax) { - throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(), - "Minimum energy < maximum energy required."); - } - - // Allocate energy - GEnergy energy; - - // Update cache - update_mc_cache(emin, emax); - - // Initialise energy - double eng; - - // Initialse acceptance fraction - double acceptance_fraction; - double inv_ecut = m_lambda.value(); - - // Use rejection method to draw a random energy. We first draw - // analytically from a power law, and then compare the power law - // at the drawn energy to the exponentially cut off function. This - // gives an acceptance fraction, and we accept the energy only if - // a uniform random number is <= the acceptance fraction. - do { - - // Get uniform random number - double u = ran.uniform(); - - // Case A: Index is not -1 - if (index() != -1.0) { - if (u > 0.0) { - eng = std::exp(std::log(u * m_mc_pow_ewidth + m_mc_pow_emin) / - m_mc_exponent); - } - else { - eng = 0.0; - } - } - - // Case B: Index is -1 - else { - eng = std::exp(u * m_mc_pow_ewidth + m_mc_pow_emin); - } - - // Compute acceptance fraction - acceptance_fraction = std::exp(-eng * inv_ecut); - - } while (ran.uniform() > acceptance_fraction); - - // Set energy - energy.MeV(eng); - - // Return energy - return energy; -} - - -/***********************************************************************//** - * @brief Read model from XML element - * - * @param[in] xml XML element containing power law model information. - * - * @exception GException::model_invalid_parnum - * Invalid number of model parameters found in XML element. - * @exception GException::model_invalid_parnames - * Invalid model parameter names found in XML element. - * - * Reads the spectral information from an XML element. The format of the XML - * elements is - * - * - * - * - * - * - * - ***************************************************************************/ -void GModelSpectralExpPlaw2::read(const GXmlElement& xml) -{ - // Verify that XML element has exactly 4 parameters - if (xml.elements() != 4 || xml.elements("parameter") != 4) { - throw GException::model_invalid_parnum(G_READ, xml, - "Power law model requires exactly 4 parameters."); - } - - // Extract model parameters - int npar[] = {0, 0, 0, 0}; - for (int i = 0; i < 4; ++i) { - - // Get parameter element - const GXmlElement* par = xml.element("parameter", i); - - // Handle prefactor - if (par->attribute("name") == "Prefactor") { - m_norm.read(*par); - npar[0]++; - } - - // Handle index - else if (par->attribute("name") == "Index") { - m_index.read(*par); - npar[1]++; - } - - // Handle cutoff parameter - else if (par->attribute("name") == "lambda") { - m_lambda.read(*par); - npar[2]++; - } - - // Handle pivot energy - else if (par->attribute("name") == "Scale") { - m_pivot.read(*par); - npar[3]++; - } - - } // endfor: looped over all parameters - - // Verify that all parameters were found - if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1 || npar[3] != 1) { - throw GException::model_invalid_parnames(G_READ, xml, - "Require \"Prefactor\", \"Index\", \"lambda\" and \"Scale\"" - " parameters."); - } - - // Return - return; -} - - -/***********************************************************************//** - * @brief Write model into XML element - * - * @param[in] xml XML element into which model information is written. - * - * @exception GException::model_invalid_spectral - * Existing XML element is not of type "ExpCutoff" - * @exception GException::model_invalid_parnum - * Invalid number of model parameters or nodes found in XML element. - * @exception GException::model_invalid_parnames - * Invalid model parameter names found in XML element. - * - * Writes the spectral information into an XML element. The format of the XML - * element is - * - * - * - * - * - * - * - ***************************************************************************/ -void GModelSpectralExpPlaw2::write(GXmlElement& xml) const -{ - // Set model type - if (xml.attribute("type") == "") { - xml.attribute("type", "ExpCutoff2"); - } - - // Verify model type - if (xml.attribute("type") != "ExpCutoff2") { - throw GException::model_invalid_spectral(G_WRITE, xml.attribute("type"), - "Spectral model is not of type \"ExpCutoff2\"."); - } - - // If XML element has 0 nodes then append 4 parameter nodes - if (xml.elements() == 0) { - xml.append(GXmlElement("parameter name=\"Prefactor\"")); - xml.append(GXmlElement("parameter name=\"Index\"")); - xml.append(GXmlElement("parameter name=\"lambda\"")); - xml.append(GXmlElement("parameter name=\"Scale\"")); - } - - // Verify that XML element has exactly 4 parameters - if (xml.elements() != 4 || xml.elements("parameter") != 4) { - throw GException::model_invalid_parnum(G_WRITE, xml, - "Power law model requires exactly 4 parameters."); - } - - // Set or update model parameter attributes - int npar[] = {0, 0, 0, 0}; - for (int i = 0; i < 4; ++i) { - - // Get parameter element - GXmlElement* par = xml.element("parameter", i); - - // Handle prefactor - if (par->attribute("name") == "Prefactor") { - npar[0]++; - m_norm.write(*par); - } - - // Handle index - else if (par->attribute("name") == "Index") { - npar[1]++; - m_index.write(*par); - } - - // Handle cutoff parameter - else if (par->attribute("name") == "lambda") { - npar[2]++; - m_lambda.write(*par); - } - - // Handle pivot energy - else if (par->attribute("name") == "Scale") { - m_pivot.write(*par); - npar[3]++; - } - - } // endfor: looped over all parameters - - // Check of all required parameters are present - if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1 || npar[3] != 1) { - throw GException::model_invalid_parnames(G_WRITE, xml, - "Require \"Prefactor\", \"Index\", \"lambda\" and \"Scale\"" - " parameters."); - } - - // Return - return; -} - - -/***********************************************************************//** - * @brief Print model information - * - * @param[in] chatter Chattiness (defaults to NORMAL). - * @return String containing model information. - ***************************************************************************/ -std::string GModelSpectralExpPlaw2::print(const GChatter& chatter) const -{ - // Initialise result string - std::string result; - - // Continue only if chatter is not silent - if (chatter != SILENT) { - - // Append header - result.append("=== GModelSpectralExpPlaw2 ==="); - - // Append information - result.append("\n"+gammalib::parformat("Number of parameters")); - result.append(gammalib::str(size())); - for (int i = 0; i < size(); ++i) { - result.append("\n"+m_pars[i]->print(chatter)); - } - - } // endif: chatter was not silent - - // Return result - return result; -} - - -/*========================================================================== - = = - = Private methods = - = = - ==========================================================================*/ - -/***********************************************************************//** - * @brief Initialise class members - ***************************************************************************/ -void GModelSpectralExpPlaw2::init_members(void) -{ - // Initialise powerlaw normalisation - m_norm.clear(); - m_norm.name("Prefactor"); - m_norm.unit("ph/cm2/s/MeV"); - m_norm.scale(1.0); - m_norm.value(1.0); // default: 1.0 - m_norm.min(0.0); // min: 0.0 - m_norm.free(); - m_norm.gradient(0.0); - m_norm.has_grad(true); - - // Initialise powerlaw index - m_index.clear(); - m_index.name("Index"); - m_index.scale(1.0); - m_index.value(-2.0); // default: -2.0 - m_index.range(-10.0,+10.0); // range: [-10,+10] - m_index.free(); - m_index.gradient(0.0); - m_index.has_grad(true); - - // Initialise cut off parameter - m_lambda.clear(); - m_lambda.name("lambda"); - m_lambda.unit("1/MeV"); - m_lambda.scale(1.0); - m_lambda.value(0.001); // default: 0.001 - m_lambda.min(0.0); // min: 0.0 - m_lambda.max(10.0); // max: 10.0 - m_lambda.free(); - m_lambda.gradient(0.0); - m_lambda.has_grad(true); - - // Initialise pivot energy - m_pivot.clear(); - m_pivot.name("PivotEnergy"); - m_pivot.unit("MeV"); - m_pivot.scale(1.0); - m_pivot.value(100.0); - m_pivot.fix(); - m_pivot.gradient(0.0); - m_pivot.has_grad(true); - - // Set parameter pointer(s) - m_pars.clear(); - m_pars.push_back(&m_norm); - m_pars.push_back(&m_index); - m_pars.push_back(&m_lambda); - m_pars.push_back(&m_pivot); - - // Initialise eval cache - m_last_energy.clear(); - m_last_index = 1.0e30; - m_last_lambda = 1.0e30; - m_last_pivot = 1.0e30; - m_last_e_norm = 0.0; - m_last_e_lambda = 0.0; - m_last_power = 0.0; - - // Initialise MC cache - m_mc_emin = 0.0; - m_mc_emax = 0.0; - m_mc_exponent = 0.0; - m_mc_pow_emin = 0.0; - m_mc_pow_ewidth = 0.0; - - // Return - return; -} - - -/***********************************************************************//** - * @brief Copy class members - * - * @param[in] model Exponential cut off power law model. - ***************************************************************************/ -void GModelSpectralExpPlaw2::copy_members(const GModelSpectralExpPlaw2& model) -{ - // Copy members - m_norm = model.m_norm; - m_index = model.m_index; - m_lambda = model.m_lambda; - m_pivot = model.m_pivot; - - // Set parameter pointer(s) - m_pars.clear(); - m_pars.push_back(&m_norm); - m_pars.push_back(&m_index); - m_pars.push_back(&m_lambda); - m_pars.push_back(&m_pivot); - - // Copy eval cache - m_last_energy = model.m_last_energy; - m_last_index = model.m_last_index; - m_last_lambda = model.m_last_lambda; - m_last_pivot = model.m_last_pivot; - m_last_e_norm = model.m_last_e_norm; - m_last_e_lambda = model.m_last_e_lambda; - m_last_power = model.m_last_power; - - // Copy MC cache - m_mc_emin = model.m_mc_emin; - m_mc_emax = model.m_mc_emax; - m_mc_exponent = model.m_mc_exponent; - m_mc_pow_emin = model.m_mc_pow_emin; - m_mc_pow_ewidth = model.m_mc_pow_ewidth; - - // Return - return; -} - - -/***********************************************************************//** - * @brief Delete class members - ***************************************************************************/ -void GModelSpectralExpPlaw2::free_members(void) -{ - // Return - return; -} - - -/***********************************************************************//** - * @brief Update eval precomputation cache - * - * @param[in] energy Energy. - * - * Updates the precomputation cache for eval() and eval_gradients() methods. - ***************************************************************************/ -void GModelSpectralExpPlaw2::update_eval_cache(const GEnergy& energy) const -{ - // Get parameter values (takes 3 multiplications which are difficult - // to avoid) - double index = m_index.value(); - double lambda = m_lambda.value(); - double pivot = m_pivot.value(); - - // If the energy or one of the parameters index, lambda or pivot - // energy has changed then recompute the cache - if ((m_last_energy != energy) || - (m_last_index != index) || - (m_last_lambda != lambda) || - (m_last_pivot != pivot)) { - - // Store actual energy and parameter values - m_last_energy = energy; - m_last_index = index; - m_last_lambda = lambda; - m_last_pivot = pivot; - - // Compute and store value - double eng = energy.MeV(); - m_last_e_norm = eng / m_last_pivot; - m_last_e_lambda = eng * m_last_lambda; - m_last_power = std::pow(m_last_e_norm, m_last_index) * - std::exp(-m_last_e_lambda); - } // endif: recomputation was required - - // Return - return; -} - - -/***********************************************************************//** - * @brief Update Monte Carlo pre computation cache - * - * @param[in] emin Minimum photon energy. - * @param[in] emax Maximum photon energy. - * - * Updates the precomputation cache for Monte Carlo simulations. - ***************************************************************************/ -void GModelSpectralExpPlaw2::update_mc_cache(const GEnergy& emin, - const GEnergy& emax) const - -{ - // Case A: Index is not -1 - if (index() != -1.0) { - - // Change in energy boundaries? - if (emin.MeV() != m_mc_emin || emax.MeV() != m_mc_emax) { - m_mc_emin = emin.MeV(); - m_mc_emax = emax.MeV(); - m_mc_exponent = index() + 1.0; - m_mc_pow_emin = std::pow(m_mc_emin, m_mc_exponent); - m_mc_pow_ewidth = std::pow(m_mc_emax, m_mc_exponent) - m_mc_pow_emin; - } - - } - - // Case B: Index is -1 - else { - - // Change in energy boundaries? - if (emin.MeV() != m_mc_emin || emax.MeV() != m_mc_emax) { - m_mc_emin = emin.MeV(); - m_mc_emax = emax.MeV(); - m_mc_exponent = 0.0; - m_mc_pow_emin = std::log(m_mc_emin); - m_mc_pow_ewidth = std::log(m_mc_emax) - m_mc_pow_emin; - } - - } - - // Return - return; -} - - -/***********************************************************************//** - * @brief Kernel for photon flux integration - * - * @param[in] energy Energy (MeV). - ***************************************************************************/ -double GModelSpectralExpPlaw2::flux_kernel::eval(const double& energy) -{ - // Evaluate function value - double e_norm = energy * m_inv_pivot; - double e_cut = energy * m_lambda; - double power = std::pow(e_norm, m_index) * std::exp(-e_cut); - double value = m_norm * power; - - // Return value - return value; -} - - -/***********************************************************************//** - * @brief Kernel for energy flux integration - * - * @param[in] energy Energy (MeV). - ***************************************************************************/ -double GModelSpectralExpPlaw2::eflux_kernel::eval(const double& energy) -{ - // Evaluate function value - double e_norm = energy * m_inv_pivot; - double e_cut = energy * m_lambda; - double power = std::pow(e_norm, m_index) * std::exp(-e_cut); - double value = m_norm * power * energy; - - // Return value - return value; -} diff --git a/src/model/Makefile.am b/src/model/Makefile.am index f5bc5f9..c9060ad 100644 --- a/src/model/Makefile.am +++ b/src/model/Makefile.am @@ -30,8 +30,8 @@ sources = GModelPar.cpp \ GModelSpectralRegistry.cpp \ GModelSpectralBrokenPlaw.cpp \ GModelSpectralConst.cpp \ + GModelSpectralExpInvPlaw.cpp \ GModelSpectralExpPlaw.cpp \ - GModelSpectralExpPlaw2.cpp \ GModelSpectralSuperExpPlaw.cpp \ GModelSpectralFunc.cpp \ GModelSpectralGauss.cpp \ diff --git a/test/data/model_point_einvplaw.xml b/test/data/model_point_einvplaw.xml new file mode 100644 index 0000000..15081ee --- /dev/null +++ b/test/data/model_point_einvplaw.xml @@ -0,0 +1,15 @@ + + + + + + + + + + + + + + + diff --git a/test/data/model_point_eplaw2.xml b/test/data/model_point_eplaw2.xml deleted file mode 100644 index 3598800..0000000 --- a/test/data/model_point_eplaw2.xml +++ /dev/null @@ -1,15 +0,0 @@ - - - - - - - - - - - - - - - diff --git a/test/test_GModel.cpp b/test/test_GModel.cpp index 1ff2310..a98dd39 100644 --- a/test/test_GModel.cpp +++ b/test/test_GModel.cpp @@ -53,7 +53,7 @@ void TestGModel::set(void) m_xml_model_point_plaw = "data/model_point_plaw.xml"; m_xml_model_point_plaw2 = "data/model_point_plaw2.xml"; m_xml_model_point_eplaw = "data/model_point_eplaw.xml"; - m_xml_model_point_eplaw2 = "data/model_point_eplaw2.xml"; + m_xml_model_point_einvplaw = "data/model_point_einvplaw.xml"; m_xml_model_point_bplaw = "data/model_point_bplaw.xml"; m_xml_model_point_supeplaw = "data/model_point_supeplaw.xml"; m_xml_model_point_logparabola = "data/model_point_logparabola.xml"; @@ -90,7 +90,7 @@ void TestGModel::set(void) append(static_cast(&TestGModel::test_plaw), "Test GModelSpectralPlaw"); append(static_cast(&TestGModel::test_plaw2), "Test GModelSpectralPlaw2"); append(static_cast(&TestGModel::test_eplaw), "Test GModelSpectralExpPlaw"); - append(static_cast(&TestGModel::test_eplaw2), "Test GModelSpectralExpPlaw2"); + append(static_cast(&TestGModel::test_einvplaw), "Test GModelSpectralExpInvPlaw"); append(static_cast(&TestGModel::test_supeplaw), "Test GModelSpectralSuperExpPlaw"); append(static_cast(&TestGModel::test_bplaw), "Test GModelSpectralBrokenPlaw"); append(static_cast(&TestGModel::test_logparabola), "Test GModelSpectralLogParabola"); @@ -1733,16 +1733,17 @@ void TestGModel::test_eplaw(void) /***********************************************************************//** - * @brief Test GModelSpectralExpPlaw2 class + * @brief Test GModelSpectralExpInvPlaw class ***************************************************************************/ -void TestGModel::test_eplaw2(void) +void TestGModel::test_einvplaw(void) { // Test void constructor test_try("Test void constructor"); try { - GModelSpectralExpPlaw2 model; - test_assert(model.type() == "ExpCutoff2", - "Model type \"ExpCutoff2\" expected."); + GModelSpectralExpInvPlaw model; + + test_assert(model.type() == "ExponentialInverseCutoffPowerLaw", + "Model type \"ExponentialInverseCutoffPowerLaw\" expected."); test_try_success(); } catch (std::exception &e) { @@ -1752,30 +1753,44 @@ void TestGModel::test_eplaw2(void) // Test value constructor test_try("Test value constructor"); try { - GModelSpectralExpPlaw2 model(2.0, -2.1, GEnergy(100.0, "MeV"), 1.0e-3); + GModelSpectralExpInvPlaw model(2.0, -2.1, GEnergy(100.0, "MeV"), 1.0e-3); test_value(model.prefactor(), 2.0); test_value(model.index(), -2.1); test_value(model.pivot().MeV(), 100.0); - test_value(model.lambda(), 1.0e-3); + test_value(model.inverse_cutoff(), 1.0e-3); test_try_success(); } catch (std::exception &e) { test_try_failure(e); } + // Test alternative value constructor + test_try("Test alternative value constructor"); + try { + GModelSpectralExpInvPlaw model(2.0, -2.1, GEnergy(100.0, "MeV"), GEnergy(1.0e+3, "MeV")); + test_value(model.prefactor(), 2.0); + test_value(model.index(), -2.1); + test_value(model.pivot().MeV(), 100.0); + test_value(model.inverse_cutoff(), 1.0e-3); + test_try_success(); + } + catch (std::exception &e) { + test_try_failure(e); + } + // Test XML constructor and value test_try("Test XML constructor, value and gradients"); try { // Test XML constructor - GXml xml(m_xml_model_point_eplaw2); + GXml xml(m_xml_model_point_einvplaw); GXmlElement* element = xml.element(0)->element(0)->element("spectrum", 0); - GModelSpectralExpPlaw2 model(*element); + GModelSpectralExpInvPlaw model(*element); test_value(model.size(), 4); - test_assert(model.type() == "ExpCutoff2", "Expected \"ExpCutoff2\""); + test_assert(model.type() == "ExponentialInverseCutoffPowerLaw", "Expected \"ExponentialInverseCutoffPowerLaw\""); test_value(model.prefactor(), 5.7e-16); test_value(model.index(), -2.48); test_value(model.pivot().TeV(), 0.3); - test_value(model.lambda(), 1.0e-6); + test_value(model.inverse_cutoff(), 1.0e-6); // Test prefactor method model.prefactor(2.3e-16); @@ -1790,11 +1805,11 @@ void TestGModel::test_eplaw2(void) test_value(model.pivot().TeV(), 0.5); // Test cutoff parameter method - model.lambda(4.2e-9); - test_value(model.lambda(), 4.2e-9); + model.inverse_cutoff(4.2e-9); + test_value(model.inverse_cutoff(), 4.2e-9); // Test operator access - const char* strarray[] = {"Prefactor", "Index", "PivotEnergy", "lambda"}; + const char* strarray[] = {"Prefactor", "Index", "PivotEnergy", "InverseCutoff"}; for (int i = 0; i < 4; ++i) { std::string keyname(strarray[i]); model[keyname].remove_range(); // To allow setting of any value @@ -2408,6 +2423,7 @@ void TestGModel::test_spectral_model(void) test_xml_model("GModelSpectralPlaw", m_xml_model_point_plaw); test_xml_model("GModelSpectralPlaw2", m_xml_model_point_plaw2); test_xml_model("GModelSpectralExpPaw", m_xml_model_point_eplaw); + test_xml_model("GModelSpectralExpInvPaw", m_xml_model_point_einvplaw); test_xml_model("GModelSpectralBrokenPlaw", m_xml_model_point_bplaw); test_xml_model("GModelSpectralSuperExpPlaw", m_xml_model_point_supeplaw); test_xml_model("GModelSpectralLogParabola", m_xml_model_point_logparabola); diff --git a/test/test_GModel.hpp b/test/test_GModel.hpp index 272ae12..ebedc6e 100644 --- a/test/test_GModel.hpp +++ b/test/test_GModel.hpp @@ -64,7 +64,7 @@ public: void test_plaw(void); void test_plaw2(void); void test_eplaw(void); - void test_eplaw2(void); + void test_einvplaw(void); void test_bplaw(void); void test_supeplaw(void); void test_logparabola(void); @@ -89,7 +89,7 @@ private: std::string m_xml_model_point_plaw; std::string m_xml_model_point_plaw2; std::string m_xml_model_point_eplaw; - std::string m_xml_model_point_eplaw2; + std::string m_xml_model_point_einvplaw; std::string m_xml_model_point_bplaw; std::string m_xml_model_point_supeplaw; std::string m_xml_model_point_logparabola; -- libgit2 0.21.4