diff --git a/doc/dev/maths/gammalib_maths.tex b/doc/dev/maths/gammalib_maths.tex index e30ee5f..08a629f 100644 --- a/doc/dev/maths/gammalib_maths.tex +++ b/doc/dev/maths/gammalib_maths.tex @@ -347,9 +347,9 @@ The parameter value gradients for the power law are given by %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\subsection{PowerLaw2} +\subsection{PowerLawPhotonFlux} -This flavour of the power law spectral model uses the integral flux $f$ within the +This flavour of the power law spectral model uses the integral photon flux $f$ within the energy range $E_{\rm min}$ and $E_{\rm max}$ as free parameter instead of the normalization $k$. The integral flux $f$ is given by @@ -380,7 +380,7 @@ is obtained from Eq.~(\ref{eq:fluxint}). Each of the 4 parameters is factorised into a scaling factor and a value, e.g. $f=f_s f_v$ and $\gamma=\gamma_s \gamma_v$. It is assumed that $E_{\rm min}$ and $E_{\rm max}$ are fixed parameters, -and {\tt GModelSpectralPlaw2::eval\_gradients} returns valid gradients only for +and {\tt GModelSpectralPlawPhotonFlux::eval\_gradients} returns valid gradients only for $f_v$ and $\gamma_v$. The flux value gradient $f_v$ is given by @@ -429,4 +429,88 @@ Note that for $\gamma \ne -1$ {\left( E_{\rm max}^{\gamma+1}-E_{\rm min}^{\gamma+1} \right)} \end{equation} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\subsection{PowerLawEnergyFlux} + +This flavour of the power law spectral model uses the integral energy flux $F$ within the +energy range $E_{\rm min}$ and $E_{\rm max}$ as free parameter instead of +the normalization $k$. +The integral energy flux $F$ is given by +\begin{equation} +F = \int_{E_{\rm min}}^{E_{\rm max}} I(E) {\rm d}E +\label{eq:fluxint} +\end{equation} +The power law model is then defined by +\begin{equation} +I(E) = \tilde{k} \, E^{\gamma} +\end{equation} +where +\begin{equation} + \tilde{k} = \left \{ + \begin{array}{l l} + \displaystyle + \frac{F}{\ln E_{\rm max} - \ln E_{\rm min}} + & \mbox{if $\gamma = -2$} \\ + \\ + \displaystyle + \frac{F (2+\gamma)}{E_{\rm max}^{\gamma+2}-E_{\rm min}^{\gamma+2}} + & \mbox{else} + \end{array} + \right . +\end{equation} +is obtained from Eq.~(\ref{eq:fluxenergy}). + +Each of the 4 parameters is factorised into a scaling factor and a value, e.g. +$F=F_s F_v$ and $\gamma=\gamma_s \gamma_v$. +It is assumed that $E_{\rm min}$ and $E_{\rm max}$ are fixed parameters, +and {\tt GModelSpectralPlawEnergyFlux::eval\_gradients} returns valid gradients only for +$F_v$ and $\gamma_v$. + +The flux value gradient $f_v$ is given by +\begin{equation} +\frac{\delta I}{\delta F_v} = + \frac{\delta I}{\delta F} \frac{\delta F}{\delta F_v} = + \frac{\delta I}{\delta F} F_s = + \frac{\delta \tilde{k}}{\delta F} \, E^{\gamma} F_s = + \frac{\tilde{k}}{F} \, E^{\gamma} F_s = + \frac{I(E)}{F} F_s = + \frac{I(E)}{F_v} +\end{equation} + +The index value gradient $\gamma_v$ is given by +\begin{equation} +\frac{\delta I}{\delta \gamma_v} = + \frac{\delta I}{\delta \gamma} \frac{\delta \gamma}{\delta \gamma_v} = + \frac{\delta I}{\delta \gamma} \gamma_s = + \left( \frac{\delta \tilde{k}}{\delta \gamma} \, E^{\gamma} + \tilde{k} \, E^{\gamma} \ln E \right) \gamma_s = + \left( \frac{1}{\tilde{k}} \frac{\delta \tilde{k}}{\delta \gamma} + \ln E \right) \tilde{k} E^{\gamma} \gamma_s = + \left( \frac{1}{\tilde{k}} \frac{\delta \tilde{k}}{\delta \gamma} + \ln E \right) I(E) \gamma_s +\end{equation} +where +\begin{equation} + \frac{\delta \tilde{k}}{\delta \gamma} = \left \{ + \begin{array}{l l} + \displaystyle + 0 & \mbox{if $\gamma = -2$} \\ + \\ + \\ + \displaystyle + \frac{F \left( E_{\rm max}^{\gamma+2}-E_{\rm min}^{\gamma+2} \right) - + F (1+\gamma) \left( E_{\rm max}^{\gamma+2} \ln E_{\rm max} - + E_{\rm min}^{\gamma+2} \ln E_{\rm min} \right)} + {\left( E_{\rm max}^{\gamma+2}-E_{\rm min}^{\gamma+2} \right)^2} + & \mbox{else} + \end{array} + \right . +\end{equation} +Note that for $\gamma \ne -2$ +\begin{equation} +\frac{1}{\tilde{k}} \frac{\delta \tilde{k}}{\delta \gamma} = + \frac{1}{2+\gamma} - + \frac{\left( E_{\rm max}^{\gamma+2} \ln E_{\rm max} - + E_{\rm min}^{\gamma+2} \ln E_{\rm min} \right)} + {\left( E_{\rm max}^{\gamma+2}-E_{\rm min}^{\gamma+2} \right)} +\end{equation} + \end{document} diff --git a/doc/source/develop/coding/design.rst b/doc/source/develop/coding/design.rst index 498bc6f..0451f80 100644 --- a/doc/source/develop/coding/design.rst +++ b/doc/source/develop/coding/design.rst @@ -397,11 +397,11 @@ non-accessible class members need then to be declared ``mutable`` to avoid compiler errors. Alternatively, the ``const_cast`` declaration can be used to allow member modifications within a ``const`` method. -As example we show here part of the definition of :doxy:`GModelSpectralPlaw2`: +As example we show here part of the definition of :doxy:`GModelSpectralPlawPhotonFlux`: .. code-block:: cpp - class GModelSpectralPlaw2 : public GModelSpectral { + class GModelSpectralPlawPhotonFlux : public GModelSpectral { public: virtual double eval(const GEnergy& srcEng) const; virtual void read(const GXmlElement& xml); @@ -430,11 +430,11 @@ As example we show here part of the definition of :doxy:`GModelSpectralPlaw2`: mutable double m_last_g_index; //!< Last spectral index gradient This class has an internal cache for precomputation, which is -potentially updated when :doxy:`GModelSpectralPlaw2::eval` is called. Here the corresponding code: +potentially updated when :doxy:`GModelSpectralPlawPhotonFlux::eval` is called. Here the corresponding code: .. code-block:: cpp - double GModelSpectralPlaw2::eval(const GEnergy& srcEng) const + double GModelSpectralPlawPhotonFlux::eval(const GEnergy& srcEng) const { // Update precomputed values update(srcEng); diff --git a/doc/source/users/user_manual/modules/model.rst b/doc/source/users/user_manual/modules/model.rst index fd58c70..10bda61 100644 --- a/doc/source/users/user_manual/modules/model.rst +++ b/doc/source/users/user_manual/modules/model.rst @@ -445,17 +445,17 @@ format: -An alternative power law function is defined by the :doxy:`GModelSpectralPlaw2` +An alternative power law function is defined by the :doxy:`GModelSpectralPlawPhotonFlux` class that uses the integral photon flux as parameter rather than the Prefactor: .. math:: - \frac{dN}{dE} = \frac{N(\gamma+1)E^{\gamma}} + \frac{dN}{dE} = \frac{F_{\rm Ph}(\gamma+1)E^{\gamma}} {E_{\rm max}^{\gamma+1} - E_{\rm min}^{\gamma+1}} where the parameters in the XML definition have the following mappings: -* :math:`N` = ``PhotonFlux`` +* :math:`F_{\rm Ph}` = ``PhotonFlux`` * :math:`\gamma` = ``Index`` * :math:`E_{\rm min}` = ``LowerLimit`` * :math:`E_{\rm max}` = ``UpperLimit`` @@ -487,12 +487,47 @@ format: .. note:: The UpperLimit and LowerLimit parameters are always treated as fixed and, - as should be apparent from this definition, the flux given by the Integral + as should be apparent from this definition, the flux given by the PhotonFlux parameter is over the range (LowerLimit, UpperLimit). Use of this model - allows the errors on the integrated flux to be evaluated directly by + allows the errors on the integrated photon flux to be evaluated directly by likelihood, obviating the need to propagate the errors if one is using the PowerLaw form. +Another alternative power law function is defined by the :doxy:`GModelSpectralPlawEnergyFlux` +class that uses the integral energy flux as parameter rather than the +Prefactor: + +.. math:: + \frac{dN}{dE} = \frac{F_{\rm E}(\gamma+2)E^{\gamma}} + {E_{\rm max}^{\gamma+2} - E_{\rm min}^{\gamma+2}} + +where the parameters in the XML definition have the following mappings: + +* :math:`F_{\rm E}` = ``EnergyFlux`` +* :math:`\gamma` = ``Index`` +* :math:`E_{\rm min}` = ``LowerLimit`` +* :math:`E_{\rm max}` = ``UpperLimit`` + +The XML format for specifying a power law defined by the integral energy flux +is: + +.. code-block:: xml + + + + + + + + +.. note:: + + The UpperLimit and LowerLimit parameters are always treated as fixed and, + as should be apparent from this definition, the flux given by the EnergyFlux + parameter is over the range (LowerLimit, UpperLimit). Use of this model + allows the errors on the integrated energy flux to be evaluated directly by + likelihood, obviating the need to propagate the errors if one is using + the PowerLaw form. Exponentially cut-off power law =============================== diff --git a/include/GModelSpectralPlaw2.hpp b/include/GModelSpectralPlawEnergyFlux.hpp index 1379d93..36dc293 100644 --- a/include/GModelSpectralPlaw2.hpp +++ b/include/GModelSpectralPlawEnergyFlux.hpp @@ -1,7 +1,7 @@ /*************************************************************************** - * GModelSpectralPlaw2.hpp - Spectral power law model class * + * GModelSpectralPlawEnergyFlux.hpp - Spectral power law model class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2016 by Juergen Knoedlseder * + * copyright (C) 2016 by Michael Mayer * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -19,13 +19,13 @@ * * ***************************************************************************/ /** - * @file GModelSpectralPlaw2.hpp - * @brief Flux normalized power law spectral model class interface definition - * @author Juergen Knoedlseder + * @file GModelSpectralPlawEnergyFlux.hpp + * @brief Energy flux normalized power law spectral model class interface definition + * @author Michael Mayer */ -#ifndef GMODELSPECTRALPLAW2_HPP -#define GMODELSPECTRALPLAW2_HPP +#ifndef GMODELSPECTRALPLAWENERGYFLUX_HPP +#define GMODELSPECTRALPLAWENERGYFLUX_HPP /* __ Includes ___________________________________________________________ */ #include @@ -40,60 +40,60 @@ class GXmlElement; /***********************************************************************//** - * @class GModelSpectralPlaw2 + * @class GModelSpectralPlawEnergyFlux * - * @brief Flux normalized power law spectral model class + * @brief Energy flux normalized power law spectral model class * * This class implements a power law spectrum. The model is defined by * * \f[ - * S_{\rm E}(E | t) = {\tt m\_integral} + * S_{\rm E}(E | t) = {\tt m\_energy\_flux} * \frac{{\tt m\_index}+1} - * {{\tt e\_max}^{{\tt m\_index}+1} - - * {\tt e\_min}^{{\tt m\_index}+1}} + * {{\tt e\_max}^{{\tt m\_index}+2} - + * {\tt e\_min}^{{\tt m\_index}+2}} * E^{\tt m\_index} * \f] * - * for \f${\tt m\_index} \ne -1\f$ and + * for \f${\tt m\_index} \ne -2\f$ and * * \f[ * S_{\rm E}(E | t) = - * \frac{{\tt m\_integral}} + * \frac{{\tt m\_energy\_flux}} * {\log {\tt e\_max} - \log {\tt e\_min}} * E^{\tt m\_index} * \f] * - * for \f${\tt m\_index} = -1\f$, where + * for \f${\tt m\_index} = -2\f$, where * - \f${\tt e\_min}\f$ is the minimum energy of an interval, * - \f${\tt e\_max}\f$ is the maximum energy of an interval, - * - \f${\tt m\_integral}\f$ is the integral flux between + * - \f${\tt m\_energy\_flux}\f$ is the energy flux between * \f${\tt e\_min}\f$ and \f${\tt e\_max}\f$, and * - \f${\tt m\_index}\f$ is the spectral index. ***************************************************************************/ -class GModelSpectralPlaw2 : public GModelSpectral { +class GModelSpectralPlawEnergyFlux : public GModelSpectral { public: // Constructors and destructors - GModelSpectralPlaw2(void); - GModelSpectralPlaw2(const std::string& type, - const std::string& integral, + GModelSpectralPlawEnergyFlux(void); + GModelSpectralPlawEnergyFlux(const std::string& type, + const std::string& energy_flux, const std::string& index, const std::string& emin, const std::string& emax); - GModelSpectralPlaw2(const double& integral, - const double& index, - const GEnergy& emin, - const GEnergy& emax); - explicit GModelSpectralPlaw2(const GXmlElement& xml); - GModelSpectralPlaw2(const GModelSpectralPlaw2& model); - virtual ~GModelSpectralPlaw2(void); + explicit GModelSpectralPlawEnergyFlux(const double& eflux, + const double& index, + const GEnergy& emin, + const GEnergy& emax); + explicit GModelSpectralPlawEnergyFlux(const GXmlElement& xml); + GModelSpectralPlawEnergyFlux(const GModelSpectralPlawEnergyFlux& model); + virtual ~GModelSpectralPlawEnergyFlux(void); // Operators - virtual GModelSpectralPlaw2& operator=(const GModelSpectralPlaw2& model); + virtual GModelSpectralPlawEnergyFlux& operator=(const GModelSpectralPlawEnergyFlux& model); // Implemented pure virtual base class methods virtual void clear(void); - virtual GModelSpectralPlaw2* clone(void) const; + virtual GModelSpectralPlawEnergyFlux* clone(void) const; virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, @@ -113,8 +113,8 @@ public: virtual std::string print(const GChatter& chatter = NORMAL) const; // Other methods - double integral(void) const; - void integral(const double& integral); + double energy_flux(void) const; + void energy_flux(const double& energy_flux); double index(void) const; void index(const double& index); GEnergy emin(void) const; @@ -125,45 +125,45 @@ public: protected: // Protected methods void init_members(void); - void copy_members(const GModelSpectralPlaw2& model); + void copy_members(const GModelSpectralPlawEnergyFlux& model); void free_members(void); void update(const GEnergy& srcEng) const; // Protected members - std::string m_type; //!< Model type - GModelPar m_integral; //!< Integral flux - GModelPar m_index; //!< Spectral index - GModelPar m_emin; //!< Lower energy limit (MeV) - GModelPar m_emax; //!< Upper energy limit (MeV) - - // Cached members used for pre-computations - mutable double m_log_emin; //!< Log(emin) - mutable double m_log_emax; //!< Log(emax) - mutable double m_pow_emin; //!< emin^(index+1) - mutable double m_pow_emax; //!< emax^(index+1) - mutable double m_norm; //!< Power-law normalization (for pivot energy 1 MeV) - mutable double m_g_norm; //!< Power-law normalization gradient - mutable double m_power; //!< Power-law factor - mutable double m_last_integral; //!< Last integral flux - mutable double m_last_index; //!< Last spectral index (MeV) - mutable GEnergy m_last_emin; //!< Last lower energy limit - mutable GEnergy m_last_emax; //!< Last upper energy limit - mutable GEnergy m_last_energy; //!< Last source energy - mutable double m_last_value; //!< Last function value - mutable double m_last_g_integral; //!< Last integral flux gradient - mutable double m_last_g_index; //!< Last spectral index gradient + std::string m_type; //!< Model type + GModelPar m_energy_flux; //!< Energy flux + GModelPar m_index; //!< Spectral index + GModelPar m_emin; //!< Lower energy limit (MeV) + GModelPar m_emax; //!< Upper energy limit (MeV) + + // Cached members used for pre-computations + mutable double m_log_emin; //!< Log(emin) + mutable double m_log_emax; //!< Log(emax) + mutable double m_pow_emin; //!< emin^(index+1) + mutable double m_pow_emax; //!< emax^(index+1) + mutable double m_norm; //!< Power-law normalization (for pivot energy 1 MeV) + mutable double m_g_norm; //!< Power-law normalization gradient + mutable double m_power; //!< Power-law factor + mutable double m_last_phflux; //!< Last energy flux + mutable double m_last_index; //!< Last spectral index (MeV) + mutable GEnergy m_last_emin; //!< Last lower energy limit + mutable GEnergy m_last_emax; //!< Last upper energy limit + mutable GEnergy m_last_energy; //!< Last source energy + mutable double m_last_value; //!< Last function value + mutable double m_last_g_phflux; //!< Last energy flux gradient + mutable double m_last_g_index; //!< Last spectral index gradient }; /***********************************************************************//** * @brief Return class name * - * @return String containing the class name ("GModelSpectralPlaw2"). + * @return String containing the class name ("GModelSpectralPlawEnergyFlux"). ***************************************************************************/ inline -std::string GModelSpectralPlaw2::classname(void) const +std::string GModelSpectralPlawEnergyFlux::classname(void) const { - return ("GModelSpectralPlaw2"); + return ("GModelSpectralPlawEnergyFlux"); } @@ -175,43 +175,41 @@ std::string GModelSpectralPlaw2::classname(void) const * Returns the type of the spectral power law model. ***************************************************************************/ inline -std::string GModelSpectralPlaw2::type(void) const +std::string GModelSpectralPlawEnergyFlux::type(void) const { return (m_type); } /***********************************************************************//** - * @brief Return integral flux + * @brief Return energy flux * - * @return Integral flux (ph/cm2/s). + * @return Energy flux (erg/cm2/s). * - * Returns the integral flux. + * Returns the energy flux. ***************************************************************************/ inline -double GModelSpectralPlaw2::integral(void) const +double GModelSpectralPlawEnergyFlux::energy_flux(void) const { - return (m_integral.value()); + return (m_energy_flux.value()); } /***********************************************************************//** - * @brief Set integral flux + * @brief Set energy flux * - * @param[in] integral Integral flux (ph/cm2/s). + * @param[in] energy_flux Energy flux (erg/cm2/s). * - * Sets the integral flux. + * Sets the energy flux. ***************************************************************************/ inline -void GModelSpectralPlaw2::integral(const double& integral) +void GModelSpectralPlawEnergyFlux::energy_flux(const double& energy_flux) { - m_integral.value(integral); + m_energy_flux.value(energy_flux); return; } - - /***********************************************************************//** * @brief Return power law index * @@ -220,7 +218,7 @@ void GModelSpectralPlaw2::integral(const double& integral) * Returns the power law index. ***************************************************************************/ inline -double GModelSpectralPlaw2::index(void) const +double GModelSpectralPlawEnergyFlux::index(void) const { return (m_index.value()); } @@ -234,7 +232,7 @@ double GModelSpectralPlaw2::index(void) const * Sets the power law index. ***************************************************************************/ inline -void GModelSpectralPlaw2::index(const double& index) +void GModelSpectralPlawEnergyFlux::index(const double& index) { m_index.value(index); return; @@ -249,7 +247,7 @@ void GModelSpectralPlaw2::index(const double& index) * Returns the minimum energy. ***************************************************************************/ inline -GEnergy GModelSpectralPlaw2::emin(void) const +GEnergy GModelSpectralPlawEnergyFlux::emin(void) const { GEnergy energy; energy.MeV(m_emin.value()); @@ -265,7 +263,7 @@ GEnergy GModelSpectralPlaw2::emin(void) const * Sets the minimum energy. ***************************************************************************/ inline -void GModelSpectralPlaw2::emin(const GEnergy& emin) +void GModelSpectralPlawEnergyFlux::emin(const GEnergy& emin) { m_emin.value(emin.MeV()); return; @@ -280,7 +278,7 @@ void GModelSpectralPlaw2::emin(const GEnergy& emin) * Returns the maximum energy. ***************************************************************************/ inline -GEnergy GModelSpectralPlaw2::emax(void) const +GEnergy GModelSpectralPlawEnergyFlux::emax(void) const { GEnergy energy; energy.MeV(m_emax.value()); @@ -296,10 +294,10 @@ GEnergy GModelSpectralPlaw2::emax(void) const * Sets the maximum energy. ***************************************************************************/ inline -void GModelSpectralPlaw2::emax(const GEnergy& emax) +void GModelSpectralPlawEnergyFlux::emax(const GEnergy& emax) { m_emax.value(emax.MeV()); return; } -#endif /* GMODELSPECTRALPLAW2_HPP */ +#endif /* GMODELSPECTRALPLAWENERGYFLUX_HPP */ diff --git a/include/GModelSpectralPlawPhotonFlux.hpp b/include/GModelSpectralPlawPhotonFlux.hpp new file mode 100644 index 0000000..11e55ee --- /dev/null +++ b/include/GModelSpectralPlawPhotonFlux.hpp @@ -0,0 +1,305 @@ +/*************************************************************************** + * GModelSpectralPlawPhotonFlux.hpp - Spectral power law model class * + * ----------------------------------------------------------------------- * + * copyright (C) 2011-2016 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GModelSpectralPlawPhotonFlux.hpp + * @brief Flux normalized power law spectral model class interface definition + * @author Juergen Knoedlseder + */ + +#ifndef GMODELSPECTRALPLAWPHOTONFLUX_HPP +#define GMODELSPECTRALPLAWPHOTONFLUX_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include "GModelSpectral.hpp" +#include "GModelPar.hpp" +#include "GEnergy.hpp" + +/* __ Forward declarations _______________________________________________ */ +class GRan; +class GTime; +class GXmlElement; + + +/***********************************************************************//** + * @class GModelSpectralPlawPhotonFlux + * + * @brief Photon flux normalized power law spectral model class + * + * This class implements a power law spectrum. The model is defined by + * + * \f[ + * S_{\rm E}(E | t) = {\tt m\_photon\_flux} + * \frac{{\tt m\_index}+1} + * {{\tt e\_max}^{{\tt m\_index}+1} - + * {\tt e\_min}^{{\tt m\_index}+1}} + * E^{\tt m\_index} + * \f] + * + * for \f${\tt m\_index} \ne -1\f$ and + * + * \f[ + * S_{\rm E}(E | t) = + * \frac{{\tt m\_photon\_flux}} + * {\log {\tt e\_max} - \log {\tt e\_min}} + * E^{\tt m\_index} + * \f] + * + * for \f${\tt m\_index} = -1\f$, where + * - \f${\tt e\_min}\f$ is the minimum energy of an interval, + * - \f${\tt e\_max}\f$ is the maximum energy of an interval, + * - \f${\tt m\_photon\_flux}\f$ is the photon flux between + * \f${\tt e\_min}\f$ and \f${\tt e\_max}\f$, and + * - \f${\tt m\_index}\f$ is the spectral index. + ***************************************************************************/ +class GModelSpectralPlawPhotonFlux : public GModelSpectral { + +public: + // Constructors and destructors + GModelSpectralPlawPhotonFlux(void); + GModelSpectralPlawPhotonFlux(const std::string& type, + const std::string& integral, + const std::string& index, + const std::string& emin, + const std::string& emax); + GModelSpectralPlawPhotonFlux(const double& integral, + const double& index, + const GEnergy& emin, + const GEnergy& emax); + explicit GModelSpectralPlawPhotonFlux(const GXmlElement& xml); + GModelSpectralPlawPhotonFlux(const GModelSpectralPlawPhotonFlux& model); + virtual ~GModelSpectralPlawPhotonFlux(void); + + // Operators + virtual GModelSpectralPlawPhotonFlux& operator=(const GModelSpectralPlawPhotonFlux& model); + + // Implemented pure virtual base class methods + virtual void clear(void); + virtual GModelSpectralPlawPhotonFlux* 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 photon_flux(void) const; + void photon_flux(const double& photon_flux); + double index(void) const; + void index(const double& index); + GEnergy emin(void) const; + void emin(const GEnergy& emin); + GEnergy emax(void) const; + void emax(const GEnergy& emax); + +protected: + // Protected methods + void init_members(void); + void copy_members(const GModelSpectralPlawPhotonFlux& model); + void free_members(void); + void update(const GEnergy& srcEng) const; + + // Protected members + std::string m_type; //!< Model type + GModelPar m_photon_flux; //!< Photon flux + GModelPar m_index; //!< Spectral index + GModelPar m_emin; //!< Lower energy limit (MeV) + GModelPar m_emax; //!< Upper energy limit (MeV) + + // Cached members used for pre-computations + mutable double m_log_emin; //!< Log(emin) + mutable double m_log_emax; //!< Log(emax) + mutable double m_pow_emin; //!< emin^(index+1) + mutable double m_pow_emax; //!< emax^(index+1) + mutable double m_norm; //!< Power-law normalization (for pivot energy 1 MeV) + mutable double m_g_norm; //!< Power-law normalization gradient + mutable double m_power; //!< Power-law factor + mutable double m_last_phflux; //!< Last photon flux + mutable double m_last_index; //!< Last spectral index (MeV) + mutable GEnergy m_last_emin; //!< Last lower energy limit + mutable GEnergy m_last_emax; //!< Last upper energy limit + mutable GEnergy m_last_energy; //!< Last source energy + mutable double m_last_value; //!< Last function value + mutable double m_last_g_phflux; //!< Last photon flux gradient + mutable double m_last_g_index; //!< Last spectral index gradient +}; + + +/***********************************************************************//** + * @brief Return class name + * + * @return String containing the class name ("GModelSpectralPlawPhotonFlux"). + ***************************************************************************/ +inline +std::string GModelSpectralPlawPhotonFlux::classname(void) const +{ + return ("GModelSpectralPlawPhotonFlux"); +} + + +/***********************************************************************//** + * @brief Return model type + * + * @return Model type. + * + * Returns the type of the spectral power law model. + ***************************************************************************/ +inline +std::string GModelSpectralPlawPhotonFlux::type(void) const +{ + return (m_type); +} + + +/***********************************************************************//** + * @brief Return photon flux + * + * @return Photon flux (ph/cm2/s). + * + * Returns the photon flux. + ***************************************************************************/ +inline +double GModelSpectralPlawPhotonFlux::photon_flux(void) const +{ + return (m_photon_flux.value()); +} + + +/***********************************************************************//** + * @brief Set photon flux + * + * @param[in] photon_flux Photon flux (ph/cm2/s). + * + * Sets the photon flux. + ***************************************************************************/ +inline +void GModelSpectralPlawPhotonFlux::photon_flux(const double& photon_flux) +{ + m_photon_flux.value(photon_flux); + return; +} + + + + +/***********************************************************************//** + * @brief Return power law index + * + * @return Power law index. + * + * Returns the power law index. + ***************************************************************************/ +inline +double GModelSpectralPlawPhotonFlux::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 GModelSpectralPlawPhotonFlux::index(const double& index) +{ + m_index.value(index); + return; +} + + +/***********************************************************************//** + * @brief Return minimum energy + * + * @return Minimum energy. + * + * Returns the minimum energy. + ***************************************************************************/ +inline +GEnergy GModelSpectralPlawPhotonFlux::emin(void) const +{ + GEnergy energy; + energy.MeV(m_emin.value()); + return energy; +} + + +/***********************************************************************//** + * @brief Set minimum energy + * + * @param[in] emin Minimum energy. + * + * Sets the minimum energy. + ***************************************************************************/ +inline +void GModelSpectralPlawPhotonFlux::emin(const GEnergy& emin) +{ + m_emin.value(emin.MeV()); + return; +} + + +/***********************************************************************//** + * @brief Return maximum energy + * + * @return Maximum energy. + * + * Returns the maximum energy. + ***************************************************************************/ +inline +GEnergy GModelSpectralPlawPhotonFlux::emax(void) const +{ + GEnergy energy; + energy.MeV(m_emax.value()); + return energy; +} + + +/***********************************************************************//** + * @brief Set maximum energy + * + * @param[in] emax Maximum energy. + * + * Sets the maximum energy. + ***************************************************************************/ +inline +void GModelSpectralPlawPhotonFlux::emax(const GEnergy& emax) +{ + m_emax.value(emax.MeV()); + return; +} + +#endif /* GMODELSPECTRALPLAWPHOTONFLUX_HPP */ diff --git a/include/GammaLib.hpp b/include/GammaLib.hpp index 5de31eb..fadad96 100644 --- a/include/GammaLib.hpp +++ b/include/GammaLib.hpp @@ -215,7 +215,8 @@ #include "GModelSpectralLogParabola.hpp" #include "GModelSpectralNodes.hpp" #include "GModelSpectralPlaw.hpp" -#include "GModelSpectralPlaw2.hpp" +#include "GModelSpectralPlawPhotonFlux.hpp" +#include "GModelSpectralPlawenergyFlux.hpp" #include "GModelTemporal.hpp" #include "GModelTemporalRegistry.hpp" #include "GModelTemporalConst.hpp" diff --git a/include/Makefile.am b/include/Makefile.am index f34f05d..5f27cd4 100644 --- a/include/Makefile.am +++ b/include/Makefile.am @@ -173,7 +173,8 @@ pkginclude_HEADERS = GammaLib.hpp \ GModelSpectralLogParabola.hpp \ GModelSpectralNodes.hpp \ GModelSpectralPlaw.hpp \ - GModelSpectralPlaw2.hpp \ + GModelSpectralPlawPhotonFlux.hpp \ + GModelSpectralPlawEnergyFlux.hpp \ GModelTemporal.hpp \ GModelTemporalRegistry.hpp \ GModelTemporalConst.hpp diff --git a/pyext/GModelSpectralPlawEnergyFlux.i b/pyext/GModelSpectralPlawEnergyFlux.i new file mode 100644 index 0000000..b5c133b --- /dev/null +++ b/pyext/GModelSpectralPlawEnergyFlux.i @@ -0,0 +1,93 @@ +/*************************************************************************** + * GModelSpectralPlawEnergyFlux.i - Spectral power law model class * + * ----------------------------------------------------------------------- * + * copyright (C) 2016 by Michael Mayer * + * ----------------------------------------------------------------------- * + * * + * 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 GModelSpectralPlawEnergyFlux.i + * @brief Energy flux normalized power law spectral model class Python interface definition + * @author Michael Mayer + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GModelSpectralPlawEnergyFlux.hpp" +%} + + +/***********************************************************************//** + * @class GModelSpectralPlawEnergyFlux + * + * @brief Flux normalized power law spectral model class + ***************************************************************************/ +class GModelSpectralPlawEnergyFlux : public GModelSpectral { +public: + // Constructors and destructors + GModelSpectralPlawEnergyFlux(void); + GModelSpectralPlawEnergyFlux(const std::string& type, + const std::string& energy_flux, + const std::string& index, + const std::string& emin, + const std::string& emax); + explicit GModelSpectralPlawEnergyFlux(const double& energy_flux, + const double& index, + const GEnergy& emin, + const GEnergy& emax); + explicit GModelSpectralPlawEnergyFlux(const GXmlElement& xml); + GModelSpectralPlawEnergyFlux(const GModelSpectralPlawEnergyFlux& model); + virtual ~GModelSpectralPlawEnergyFlux(void); + + // Implemented pure virtual methods + virtual void clear(void); + virtual GModelSpectralPlawEnergyFlux* 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; + + // Other methods + double energy_flux(void) const; + void energy_flux(const double& energy_flux); + double index(void) const; + void index(const double& index); + GEnergy emin(void) const; + void emin(const GEnergy& emin); + GEnergy emax(void) const; + void emax(const GEnergy& emax); +}; + + +/***********************************************************************//** + * @brief GModelSpectralPlawEnergyFlux class extension + ***************************************************************************/ +%extend GModelSpectralPlawEnergyFlux { + GModelSpectralPlawEnergyFlux copy() { + return (*self); + } +}; diff --git a/pyext/GModelSpectralPlaw2.i b/pyext/GModelSpectralPlawPhotonFlux.i index 78893e3..e50e32e 100644 --- a/pyext/GModelSpectralPlaw2.i +++ b/pyext/GModelSpectralPlawPhotonFlux.i @@ -1,5 +1,5 @@ /*************************************************************************** - * GModelSpectralPlaw2.i - Spectral power law model class * + * GModelSpectralPlawPhotonFlux.i - Spectral power law model class * * ----------------------------------------------------------------------- * * copyright (C) 2011-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * @@ -19,41 +19,41 @@ * * ***************************************************************************/ /** - * @file GModelSpectralPlaw2.i + * @file GModelSpectralPlawPhotonFlux.i * @brief Flux normalized power law spectral model class Python interface definition * @author Juergen Knoedlseder */ %{ /* Put headers and other declarations here that are needed for compilation */ -#include "GModelSpectralPlaw2.hpp" +#include "GModelSpectralPlawPhotonFlux.hpp" %} /***********************************************************************//** - * @class GModelSpectralPlaw2 + * @class GModelSpectralPlawPhotonFlux * * @brief Flux normalized power law spectral model class ***************************************************************************/ -class GModelSpectralPlaw2 : public GModelSpectral { +class GModelSpectralPlawPhotonFlux : public GModelSpectral { public: // Constructors and destructors - GModelSpectralPlaw2(void); - GModelSpectralPlaw2(const std::string& type, - const std::string& integral, + GModelSpectralPlawPhotonFlux(void); + GModelSpectralPlawPhotonFlux(const std::string& type, + const std::string& photon_flux, const std::string& index, const std::string& emin, const std::string& emax); - GModelSpectralPlaw2(const double& integral, + GModelSpectralPlawPhotonFlux(const double& photon_flux, const double& index, const GEnergy& emin, const GEnergy& emax); - explicit GModelSpectralPlaw2(const GXmlElement& xml); - GModelSpectralPlaw2(const GModelSpectralPlaw2& model); - virtual ~GModelSpectralPlaw2(void); + explicit GModelSpectralPlawPhotonFlux(const GXmlElement& xml); + GModelSpectralPlawPhotonFlux(const GModelSpectralPlawPhotonFlux& model); + virtual ~GModelSpectralPlawPhotonFlux(void); // Implemented pure virtual methods virtual void clear(void); - virtual GModelSpectralPlaw2* clone(void) const; + virtual GModelSpectralPlawPhotonFlux* clone(void) const; virtual std::string classname(void) const; virtual std::string type(void) const; virtual double eval(const GEnergy& srcEng, @@ -72,8 +72,8 @@ public: virtual void write(GXmlElement& xml) const; // Other methods - double integral(void) const; - void integral(const double& integral); + double photon_flux(void) const; + void photon_flux(const double& photon_flux); double index(void) const; void index(const double& index); GEnergy emin(void) const; @@ -84,10 +84,10 @@ public: /***********************************************************************//** - * @brief GModelSpectralPlaw2 class extension + * @brief GModelSpectralPlawPhotonFlux class extension ***************************************************************************/ -%extend GModelSpectralPlaw2 { - GModelSpectralPlaw2 copy() { +%extend GModelSpectralPlawPhotonFlux { + GModelSpectralPlawPhotonFlux copy() { return (*self); } }; diff --git a/pyext/gammalib/model.i b/pyext/gammalib/model.i index d17a4a2..7b4ed9f 100644 --- a/pyext/gammalib/model.i +++ b/pyext/gammalib/model.i @@ -152,7 +152,8 @@ %include "GModelSpectralLogParabola.i" %include "GModelSpectralNodes.i" %include "GModelSpectralPlaw.i" -%include "GModelSpectralPlaw2.i" +%include "GModelSpectralPlawPhotonFlux.i" +%include "GModelSpectralPlawEnergyFlux.i" %include "GModelTemporal.i" %include "GModelTemporalRegistry.i" %include "GModelTemporalConst.i" diff --git a/src/model/GModelSpectralPlawEnergyFlux.cpp b/src/model/GModelSpectralPlawEnergyFlux.cpp new file mode 100644 index 0000000..aa10979 --- /dev/null +++ b/src/model/GModelSpectralPlawEnergyFlux.cpp @@ -0,0 +1,934 @@ +/*************************************************************************** + * GModelSpectralPlawEnergyFlux.cpp - Spectral power law model class * + * ----------------------------------------------------------------------- * + * copyright (C) 2016 by Michael Mayer * + * ----------------------------------------------------------------------- * + * * + * 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 GModelSpectralPlawEnergyFlux.cpp + * @brief Energy flux normalized power law spectral model class implementation + * @author Michael Mayer + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "GException.hpp" +#include "GTools.hpp" +#include "GMath.hpp" +#include "GRan.hpp" +#include "GModelSpectralPlawEnergyFlux.hpp" +#include "GModelSpectralRegistry.hpp" + +/* __ Constants __________________________________________________________ */ + +/* __ Globals ____________________________________________________________ */ +const GModelSpectralPlawEnergyFlux g_spectral_plaw_eflux_seed("PowerLaw", + "EnergyFlux", + "Index", + "LowerLimit", + "UpperLimit"); +const GModelSpectralRegistry g_spectral_plaw_eflux_registry(&g_spectral_plaw_eflux_seed); + +/* __ Method name definitions ____________________________________________ */ +#define G_FLUX "GModelSpectralPlawEnergyFlux::flux(GEnergy&, GEnergy&)" +#define G_EFLUX "GModelSpectralPlawEnergyFlux::eflux(GEnergy&, GEnergy&)" +#define G_MC "GModelSpectralPlawEnergyFlux::mc(GEnergy&, GEnergy&, GTime&, GRan&)" +#define G_READ "GModelSpectralPlawEnergyFlux::read(GXmlElement&)" +#define G_WRITE "GModelSpectralPlawEnergyFlux::write(GXmlElement&)" + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + +/* __ Debug definitions __________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +GModelSpectralPlawEnergyFlux::GModelSpectralPlawEnergyFlux(void) : GModelSpectral() +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Model type and parameter name constructor + * + * @param[in] type Model type. + * @param[in] energy_flux Name of energy flux parameter. + * @param[in] index Name of index parameter. + * @param[in] emin Name of emin parameter. + * @param[in] emax Name of emax parameter. + ***************************************************************************/ +GModelSpectralPlawEnergyFlux::GModelSpectralPlawEnergyFlux(const std::string& type, + const std::string& energy_flux, + const std::string& index, + const std::string& emin, + const std::string& emax) : + GModelSpectral() +{ + // Initialise members + init_members(); + + // Set model type + m_type = type; + + // Set parameter names + m_energy_flux.name(energy_flux); + m_index.name(index); + m_emin.name(emin); + m_emax.name(emax); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Constructor + * + * @param[in] energy_flux Energy flux (erg/cm2/s). + * @param[in] index Power law index. + * @param[in] emin Minimum energy. + * @param[in] emax Maximum energy. + * + * Construct a spectral power law from the + * - energy flux (in erg/cm2/s), + * - spectral index, + * - minimum energy and + * - maximum energy. + ***************************************************************************/ +GModelSpectralPlawEnergyFlux::GModelSpectralPlawEnergyFlux(const double& energy_flux, + const double& index, + const GEnergy& emin, + const GEnergy& emax) : + GModelSpectral() +{ + // Initialise members + init_members(); + + // Set parameters + m_energy_flux.value(energy_flux); + m_index.value(index); + m_emin.value(emin.MeV()); + m_emax.value(emax.MeV()); + + // Return + return; +} + + +/***********************************************************************//** + * @brief XML constructor + * + * @param[in] xml XML element. + * + * Constructs energy flux normalized 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. + ***************************************************************************/ +GModelSpectralPlawEnergyFlux::GModelSpectralPlawEnergyFlux(const GXmlElement& xml) : + GModelSpectral() +{ + // Initialise members + init_members(); + + // Read information from XML element + read(xml); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] model Spectral power law model. + ***************************************************************************/ +GModelSpectralPlawEnergyFlux::GModelSpectralPlawEnergyFlux(const GModelSpectralPlawEnergyFlux& model) : + GModelSpectral(model) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(model); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +GModelSpectralPlawEnergyFlux::~GModelSpectralPlawEnergyFlux(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] model Spectral power law model. + * @return Spectral power law model. + ***************************************************************************/ +GModelSpectralPlawEnergyFlux& GModelSpectralPlawEnergyFlux::operator=(const GModelSpectralPlawEnergyFlux& 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 power law model + ***************************************************************************/ +void GModelSpectralPlawEnergyFlux::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 power law model + * + * @return Pointer to deep copy of power law model + ***************************************************************************/ +GModelSpectralPlawEnergyFlux* GModelSpectralPlawEnergyFlux::clone(void) const +{ + // Clone power law model + return new GModelSpectralPlawEnergyFlux(*this); +} + + +/***********************************************************************//** + * @brief Evaluate function + * + * @param[in] srcEng True photon energy. + * @param[in] srcTime True photon arrival time. + * @return Model value (ph/cm2/s/MeV). + * + * Computes + * + * \f[ + * S_{\rm E}(E | t) = {\tt m\_energy\_flux} + * \frac{{\tt m\_index}+2} + * {{\tt e\_max}^{{\tt m\_index}+2} - + * {\tt e\_min}^{{\tt m\_index}+2}} + * E^{\tt m\_index} + * \f] + * + * for \f${\tt m\_index} \ne -2\f$ and + * + * \f[ + * S_{\rm E}(E | t) = + * \frac{{\tt m\_energy\_flux}} + * {\log {\tt e\_max} - \log {\tt e\_min}} + * E^{\tt m\_index} + * \f] + * + * for \f${\tt m\_index} = -2\f$, where + * - \f${\tt e\_min}\f$ is the minimum energy of an interval, + * - \f${\tt e\_max}\f$ is the maximum energy of an interval, + * - \f${\tt m\_energy\_flux}\f$ is the energy flux between + * \f${\tt e\_min}\f$ and \f${\tt e\_max}\f$, and + * - \f${\tt m\_index}\f$ is the spectral index. + ***************************************************************************/ +double GModelSpectralPlawEnergyFlux::eval(const GEnergy& srcEng, + const GTime& srcTime) const +{ + // Update precomputed values + update(srcEng); + + // Compute function value + double value = m_energy_flux.value() * gammalib::erg2MeV * m_norm * m_power; + + // Compile option: Check for NaN/Inf + #if defined(G_NAN_CHECK) + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { + std::cout << "*** ERROR: GModelSpectralPlawEnergyFlux::eval"; + std::cout << "(srcEng=" << srcEng; + std::cout << ", srcTime=" << srcTime << "):"; + std::cout << " NaN/Inf encountered"; + std::cout << " (value=" << value; + std::cout << ", energy flux=" << energy_flux(); + std::cout << ", m_norm=" << m_norm; + std::cout << ", m_power=" << m_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 (ph/cm2/s/MeV). + * + * Computes + * + * \f[ + * S_{\rm E}(E | t) = {\tt m\_energy\_flux} + * \frac{{\tt m\_index}+2} + * {{\tt e\_max}^{{\tt m\_index}+2} - + * {\tt e\_min}^{{\tt m\_index}+2}} + * E^{\tt m\_index} + * \f] + * + * for \f${\tt m\_index} \ne -2\f$ and + * + * \f[ + * S_{\rm E}(E | t) = + * \frac{{\tt m\_energy\_flux}} + * {\log {\tt e\_max} - \log {\tt e\_min}} + * E^{\tt m\_index} + * \f] + * + * for \f${\tt m\_index} = -2\f$, where + * - \f${\tt e\_min}\f$ is the minimum energy of an interval, + * - \f${\tt e\_max}\f$ is the maximum energy of an interval, + * - \f${\tt m\_energy\_flux}\f$ is the energy flux between + * \f${\tt e\_min}\f$ and \f${\tt e\_max}\f$, and + * - \f${\tt m\_index}\f$ is the spectral index. + * + * 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\_energy\_flux}} = + * \frac{S_{\rm E}(E | t)}{{\tt m\_energy\_flux}} + * \f] + * + * \f[ + * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_index}} = + * S_{\rm E}(E | t) \, + * \left( \frac{1}{{\tt m\_index}+2} - + * \frac{\log({\tt e\_max}) {\tt e\_max}^{{\tt m\_index}+2} - + * \log({\tt e\_min}) {\tt e\_min}^{{\tt m\_index}+2}} + * {{\tt e\_max}^{{\tt m\_index}+2} - + * {\tt e\_min}^{{\tt m\_index}+2}} + * + \ln(E) \right) + * \f] + * + * for \f${\tt m\_index} \ne -2\f$ and + * + * \f[ + * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_index}} = + * S_{\rm E}(E | t) \, \ln(E) + * \f] + * + * for \f${\tt m\_index} = -2\f$. + * + * No partial derivatives are supported for the energy boundaries. + ***************************************************************************/ +double GModelSpectralPlawEnergyFlux::eval_gradients(const GEnergy& srcEng, + const GTime& srcTime) +{ + // Initialise gradients + double g_energy_flux = 0.0; + double g_index = 0.0; + + // Update precomputed values + update(srcEng); + + // Compute function value + double value = m_energy_flux.value() * gammalib::erg2MeV * m_norm * m_power; + + // Integral flux gradient + if (m_energy_flux.is_free() && m_energy_flux.factor_value() > 0.0) { + g_energy_flux = value / m_energy_flux.factor_value(); + } + + // Index gradient + if (m_index.is_free()) { + g_index = value * (m_g_norm + gammalib::ln10 * srcEng.log10MeV()) * + m_index.scale(); + } + + // Set gradients + m_energy_flux.factor_gradient(g_energy_flux); + m_index.factor_gradient(g_index); + + // Compile option: Check for NaN/Inf + #if defined(G_NAN_CHECK) + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { + std::cout << "*** ERROR: GModelSpectralPlawEnergyFlux::eval_gradients"; + std::cout << "(srcEng=" << srcEng; + std::cout << ", srcTime=" << srcTime << "):"; + std::cout << " NaN/Inf encountered"; + std::cout << " (value=" << value; + std::cout << ", energy flux=" << energy_flux(); + std::cout << ", m_norm=" << m_norm; + std::cout << ", m_power=" << m_power; + std::cout << ", g_energy_flux=" << g_energy_flux; + std::cout << ", g_index=" << g_index; + 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 Minimum 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). + ***************************************************************************/ +double GModelSpectralPlawEnergyFlux::flux(const GEnergy& emin, + const GEnergy& emax) const +{ + // Initialise flux + double flux = 0.0; + + // Compute only if integration range is valid + if (emin < emax) { + + // Compute power law normalization + double norm = 0.0; + if (index() != -2.0) { + double gamma = m_index.value() + 2.0; + double pow_ref_emin = std::pow(this->emin().MeV(), gamma); + double pow_ref_emax = std::pow(this->emax().MeV(), gamma); + norm = m_energy_flux.value() * gammalib::erg2MeV * gamma / (pow_ref_emax - pow_ref_emin); + } + else { + double log_ref_emin = std::log(this->emin().MeV()); + double log_ref_emax = std::log(this->emax().MeV()); + norm = m_energy_flux.value() * gammalib::erg2MeV / (log_ref_emax - log_ref_emin); + } + + // Compute photon flux + if (index() != -1.0) { + double gamma = m_index.value() + 1.0; + double pow_emin = std::pow(emin.MeV(), gamma); + double pow_emax = std::pow(emax.MeV(), gamma); + flux = norm / gamma * (pow_emax - pow_emin); + } + + // Case B: Index is -1 + else { + double log_emin = std::log(emin.MeV()); + double log_emax = std::log(emax.MeV()); + flux = norm * (log_emax - log_emin); + } + + } // endif: integration range was valid + + // Return flux + return flux; +} + + +/***********************************************************************//** + * @brief Returns model energy flux between [emin, emax] (units: erg/cm2/s) + * + * @param[in] emin Minimum photon energy. + * @param[in] emax Minimum 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). + ***************************************************************************/ +double GModelSpectralPlawEnergyFlux::eflux(const GEnergy& emin, + const GEnergy& emax) const +{ + // Initialise flux + double eflux = 0.0; + + // Compute only if integration range is valid + if (emin < emax) { + + // Compute power law normalization + double norm; + if (index() != -2.0) { + double gamma = m_index.value() + 2.0; + double pow_ref_emin = std::pow(this->emin().MeV(), gamma); + double pow_ref_emax = std::pow(this->emax().MeV(), gamma); + double pow_emin = std::pow(emin.MeV(), gamma); + double pow_emax = std::pow(emax.MeV(), gamma); + double factor = (pow_emax - pow_emin) / + (pow_ref_emax - pow_ref_emin); + eflux = m_energy_flux.value() * factor; + } + else { + double log_emin = std::log(emin.MeV()); + double log_emax = std::log(emax.MeV()); + double log_ref_emin = std::log(this->emin().MeV()); + double log_ref_emax = std::log(this->emax().MeV()); + double factor = (log_emax - log_emin) / + (log_ref_emax - log_ref_emin); + eflux = m_energy_flux.value() * factor; + } + + } // endif: integration range was valid + + // Return flux + 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). + * + * Returns Monte Carlo energy by randomly drawing from a power law. + ***************************************************************************/ +GEnergy GModelSpectralPlawEnergyFlux::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; + + // Case A: Index is not -1 + if (index() != -1.0) { + double exponent = index() + 1.0; + double e_max = std::pow(emax.MeV(), exponent); + double e_min = std::pow(emin.MeV(), exponent); + double u = ran.uniform(); + double eng = (u > 0.0) + ? std::exp(std::log(u * (e_max - e_min) + e_min) / exponent) + : 0.0; + energy.MeV(eng); + } + + // Case B: Index is -1 + else { + double e_max = std::log(emax.MeV()); + double e_min = std::log(emin.MeV()); + double u = ran.uniform(); + double eng = std::exp(u * (e_max - e_min) + e_min); + energy.MeV(eng); + } + + // Return energy + return energy; +} + + +/***********************************************************************//** + * @brief Read model from XML element + * + * @param[in] xml XML element. + * + * Reads the spectral information from an XML element. + ***************************************************************************/ +void GModelSpectralPlawEnergyFlux::read(const GXmlElement& xml) +{ + // Get parameter pointers + const GXmlElement* eflux = gammalib::xml_get_par(G_READ, xml, m_energy_flux.name()); + const GXmlElement* index = gammalib::xml_get_par(G_READ, xml, m_index.name()); + const GXmlElement* emin = gammalib::xml_get_par(G_READ, xml, m_emin.name()); + const GXmlElement* emax = gammalib::xml_get_par(G_READ, xml, m_emax.name()); + + // Read parameters + m_energy_flux.read(*eflux); + m_index.read(*index); + m_emin.read(*emin); + m_emax.read(*emax); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Write model into XML element + * + * @param[in] xml XML element. + * + * @exception GException::model_invalid_spectral + * Existing XML element is not of type "PowerLaw" + * + * Writes the spectral information into an XML element. + ***************************************************************************/ +void GModelSpectralPlawEnergyFlux::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* eflux = gammalib::xml_need_par(G_WRITE, xml, m_energy_flux.name()); + GXmlElement* index = gammalib::xml_need_par(G_WRITE, xml, m_index.name()); + GXmlElement* emin = gammalib::xml_need_par(G_WRITE, xml, m_emin.name()); + GXmlElement* emax = gammalib::xml_need_par(G_WRITE, xml, m_emax.name()); + + // Write parameters + m_energy_flux.write(*eflux); + m_index.write(*index); + m_emin.write(*emin); + m_emax.write(*emax); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Print power law information + * + * @param[in] chatter Chattiness (defaults to NORMAL). + * @return String containing power law information. + ***************************************************************************/ +std::string GModelSpectralPlawEnergyFlux::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("=== GModelSpectralPlawEnergyFlux ==="); + + // 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 GModelSpectralPlawEnergyFlux::init_members(void) +{ + // Initialise model type + m_type = "PowerLaw"; + + // Initialise energy flux + m_energy_flux.clear(); + m_energy_flux.name("EnergyFlux"); + m_energy_flux.unit("erg/cm2/s"); + m_energy_flux.scale(1.0); + m_energy_flux.value(1.0); // default: 1.0 + m_energy_flux.range(0.0, 10.0); // range: [0,10] + m_energy_flux.free(); + m_energy_flux.gradient(0.0); + m_energy_flux.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 lower limit + m_emin.clear(); + m_emin.name("LowerLimit"); + m_emin.unit("MeV"); + m_emin.scale(1.0); + m_emin.value(100.0); // default: 100 + m_emin.range(0.001, 1.0e15); // range: [0.001, 1e15] + m_emin.fix(); + m_emin.gradient(0.0); + m_emin.has_grad(false); + + // Initialise upper limit + m_emax.clear(); + m_emax.name("UpperLimit"); + m_emax.unit("MeV"); + m_emax.scale(1.0); + m_emax.value(500000.0); // default: 500000 + m_emin.range(0.001, 1.0e15); // range: [0.001, 1e15] + m_emax.fix(); + m_emax.gradient(0.0); + m_emax.has_grad(false); + + // Set parameter pointer(s) + m_pars.clear(); + m_pars.push_back(&m_energy_flux); + m_pars.push_back(&m_index); + m_pars.push_back(&m_emin); + m_pars.push_back(&m_emax); + + // Initialise last parameters (for fast computation) + m_log_emin = 0.0; + m_log_emax = 0.0; + m_pow_emin = 0.0; + m_pow_emax = 0.0; + m_norm = 0.0; + m_power = 0.0; + m_last_index = 1000.0; + m_last_emin.MeV(0.0); + m_last_emax.MeV(0.0); + m_last_energy.MeV(0.0); + m_last_value = 0.0; + m_last_g_index = 0.0; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] model Spectral power law model. + ***************************************************************************/ +void GModelSpectralPlawEnergyFlux::copy_members(const GModelSpectralPlawEnergyFlux& model) +{ + // Copy members + m_type = model.m_type; + m_energy_flux = model.m_energy_flux; + m_index = model.m_index; + m_emin = model.m_emin; + m_emax = model.m_emax; + + // Set parameter pointer(s) + m_pars.clear(); + m_pars.push_back(&m_energy_flux); + m_pars.push_back(&m_index); + m_pars.push_back(&m_emin); + m_pars.push_back(&m_emax); + + // Copy bookkeeping information + m_log_emin = model.m_log_emin; + m_log_emax = model.m_log_emax; + m_pow_emin = model.m_pow_emin; + m_pow_emax = model.m_pow_emax; + m_norm = model.m_norm; + m_power = model.m_power; + m_last_index = model.m_last_index; + m_last_emin = model.m_last_emin; + m_last_emax = model.m_last_emax; + m_last_energy = model.m_last_energy; + m_last_value = model.m_last_value; + m_last_g_index = model.m_last_g_index; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void GModelSpectralPlawEnergyFlux::free_members(void) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief Update precomputed values + * + * @param[in] srcEng Source energy + ***************************************************************************/ +void GModelSpectralPlawEnergyFlux::update(const GEnergy& srcEng) const +{ + // Compute index+1 + double gamma = index() + 2.0; + + // Change in spectral index? + if (index() != m_last_index) { + + // Save actual spectral index + m_last_index = index(); + + // Change in energy boundaries? + if (emin() != m_last_emin || emax() != m_last_emax) { + m_log_emin = std::log(emin().MeV()); + m_last_emin = emin(); + m_log_emax = std::log(emax().MeV()); + m_last_emax = emax(); + } + + // Compute normalization factors + if (gamma != 0.0) { + m_pow_emin = std::pow(emin().MeV(), gamma); + m_pow_emax = std::pow(emax().MeV(), gamma); + double d = m_pow_emax - m_pow_emin; + m_norm = gamma / d; + m_g_norm = 1.0/gamma - + (m_pow_emax*m_log_emax - m_pow_emin*m_log_emin)/d; + } + else { + m_norm = 1.0 / (m_log_emax - m_log_emin); + m_g_norm = 0.0; + } + + // Update power law factor + m_power = std::pow(srcEng.MeV(), index()); + m_last_energy = srcEng; + + } // endif: change in spectral index + + // ... no change in spectral index + else { + + // Change in energy boundaries? + if (emin() != m_last_emin || emax() != m_last_emax) { + + // Update energy boundaries + m_log_emin = std::log(emin().MeV()); + m_last_emin = emin(); + m_log_emax = std::log(emax().MeV()); + m_last_emax = emax(); + + // Compute power law normalization + if (gamma != 0.0) { + m_pow_emin = std::pow(emin().MeV(), gamma); + m_pow_emax = std::pow(emax().MeV(), gamma); + double d = m_pow_emax - m_pow_emin; + m_norm = gamma / d; + m_g_norm = 1.0/gamma - + (m_pow_emax*m_log_emax - m_pow_emin*m_log_emin)/d; + } + else { + m_norm = 1.0 / (m_log_emax - m_log_emin); + m_g_norm = 0.0; + } + + } // endif: change in energy boundaries + + // Change in energy? + if (srcEng != m_last_energy) { + m_power = std::pow(srcEng.MeV(), index()); + m_last_energy = srcEng; + } + + } // endelse: no change in spectral index + + // Return + return; +} diff --git a/src/model/GModelSpectralPlaw2.cpp b/src/model/GModelSpectralPlawPhotonFlux.cpp index 7dbb393..e485bf1 100644 --- a/src/model/GModelSpectralPlaw2.cpp +++ b/src/model/GModelSpectralPlawPhotonFlux.cpp @@ -1,5 +1,5 @@ /*************************************************************************** - * GModelSpectralPlaw2.cpp - Spectral power law model class * + * GModelSpectralPlawPhotonFlux.cpp - Spectral power law model class * * ----------------------------------------------------------------------- * * copyright (C) 2010-2016 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * @@ -19,7 +19,7 @@ * * ***************************************************************************/ /** - * @file GModelSpectralPlaw2.cpp + * @file GModelSpectralPlawPhotonFlux.cpp * @brief Flux normalized power law spectral model class implementation * @author Juergen Knoedlseder */ @@ -33,33 +33,33 @@ #include "GTools.hpp" #include "GMath.hpp" #include "GRan.hpp" -#include "GModelSpectralPlaw2.hpp" +#include "GModelSpectralPlawPhotonFlux.hpp" #include "GModelSpectralRegistry.hpp" /* __ Constants __________________________________________________________ */ /* __ Globals ____________________________________________________________ */ -const GModelSpectralPlaw2 g_spectral_plaw2_seed1("PowerLaw", +const GModelSpectralPlawPhotonFlux g_spectral_plaw_phflux_seed1("PowerLaw", "PhotonFlux", "Index", "LowerLimit", "UpperLimit"); -const GModelSpectralRegistry g_spectral_plaw2_registry1(&g_spectral_plaw2_seed1); +const GModelSpectralRegistry g_spectral_plaw_phflux_registry1(&g_spectral_plaw_phflux_seed1); #if defined(G_LEGACY_XML_FORMAT) -const GModelSpectralPlaw2 g_spectral_plaw2_seed2("PowerLaw2", +const GModelSpectralPlawPhotonFlux g_spectral_plaw_phflux_seed2("PowerLaw2", "Integral", "Index", "LowerLimit", "UpperLimit"); -const GModelSpectralRegistry g_spectral_plaw2_registry2(&g_spectral_plaw2_seed2); +const GModelSpectralRegistry g_spectral_plaw_phflux_registry2(&g_spectral_plaw_phflux_seed2); #endif /* __ Method name definitions ____________________________________________ */ -#define G_FLUX "GModelSpectralPlaw2::flux(GEnergy&, GEnergy&)" -#define G_EFLUX "GModelSpectralPlaw2::eflux(GEnergy&, GEnergy&)" -#define G_MC "GModelSpectralPlaw2::mc(GEnergy&, GEnergy&, GTime&, GRan&)" -#define G_READ "GModelSpectralPlaw2::read(GXmlElement&)" -#define G_WRITE "GModelSpectralPlaw2::write(GXmlElement&)" +#define G_FLUX "GModelSpectralPlawPhotonFlux::flux(GEnergy&, GEnergy&)" +#define G_EFLUX "GModelSpectralPlawPhotonFlux::eflux(GEnergy&, GEnergy&)" +#define G_MC "GModelSpectralPlawPhotonFlux::mc(GEnergy&, GEnergy&, GTime&, GRan&)" +#define G_READ "GModelSpectralPlawPhotonFlux::read(GXmlElement&)" +#define G_WRITE "GModelSpectralPlawPhotonFlux::write(GXmlElement&)" /* __ Macros _____________________________________________________________ */ @@ -79,7 +79,7 @@ const GModelSpectralRegistry g_spectral_plaw2_registry2(&g_spectral_plaw2_seed2) * * Constructs empty power law photon flux model. ***************************************************************************/ -GModelSpectralPlaw2::GModelSpectralPlaw2(void) : GModelSpectral() +GModelSpectralPlawPhotonFlux::GModelSpectralPlawPhotonFlux(void) : GModelSpectral() { // Initialise members init_members(); @@ -93,13 +93,13 @@ GModelSpectralPlaw2::GModelSpectralPlaw2(void) : GModelSpectral() * @brief Model type and parameter name constructor * * @param[in] type Model type. - * @param[in] integral Name of integral parameter. + * @param[in] photon_flux Name of photon flux parameter. * @param[in] index Name of index parameter. * @param[in] emin Name of emin parameter. * @param[in] emax Name of emax parameter. ***************************************************************************/ -GModelSpectralPlaw2::GModelSpectralPlaw2(const std::string& type, - const std::string& integral, +GModelSpectralPlawPhotonFlux::GModelSpectralPlawPhotonFlux(const std::string& type, + const std::string& photon_flux, const std::string& index, const std::string& emin, const std::string& emax) : @@ -112,7 +112,7 @@ GModelSpectralPlaw2::GModelSpectralPlaw2(const std::string& type, m_type = type; // Set parameter names - m_integral.name(integral); + m_photon_flux.name(photon_flux); m_index.name(index); m_emin.name(emin); m_emax.name(emax); @@ -125,18 +125,18 @@ GModelSpectralPlaw2::GModelSpectralPlaw2(const std::string& type, /***********************************************************************//** * @brief Constructor * - * @param[in] integral Integral flux (ph/cm2/s). + * @param[in] photon_flux Photon flux (ph/cm2/s). * @param[in] index Power law index. * @param[in] emin Minimum energy. * @param[in] emax Maximum energy. * * Construct a spectral power law from the - * - integral flux (in ph/cm2/s), + * - integral photon flux (in ph/cm2/s), * - spectral index, * - minimum energy and * - maximum energy. ***************************************************************************/ -GModelSpectralPlaw2::GModelSpectralPlaw2(const double& integral, +GModelSpectralPlawPhotonFlux::GModelSpectralPlawPhotonFlux(const double& photon_flux, const double& index, const GEnergy& emin, const GEnergy& emax) : @@ -146,7 +146,7 @@ GModelSpectralPlaw2::GModelSpectralPlaw2(const double& integral, init_members(); // Set parameters - m_integral.value(integral); + m_photon_flux.value(photon_flux); m_index.value(index); m_emin.value(emin.MeV()); m_emax.value(emax.MeV()); @@ -165,7 +165,7 @@ GModelSpectralPlaw2::GModelSpectralPlaw2(const double& integral, * information from an XML element. See the read() method for more * information about the expected structure of the XML element. ***************************************************************************/ -GModelSpectralPlaw2::GModelSpectralPlaw2(const GXmlElement& xml) : +GModelSpectralPlawPhotonFlux::GModelSpectralPlawPhotonFlux(const GXmlElement& xml) : GModelSpectral() { // Initialise members @@ -184,7 +184,7 @@ GModelSpectralPlaw2::GModelSpectralPlaw2(const GXmlElement& xml) : * * @param[in] model Spectral power law model. ***************************************************************************/ -GModelSpectralPlaw2::GModelSpectralPlaw2(const GModelSpectralPlaw2& model) : +GModelSpectralPlawPhotonFlux::GModelSpectralPlawPhotonFlux(const GModelSpectralPlawPhotonFlux& model) : GModelSpectral(model) { // Initialise members @@ -201,7 +201,7 @@ GModelSpectralPlaw2::GModelSpectralPlaw2(const GModelSpectralPlaw2& model) : /***********************************************************************//** * @brief Destructor ***************************************************************************/ -GModelSpectralPlaw2::~GModelSpectralPlaw2(void) +GModelSpectralPlawPhotonFlux::~GModelSpectralPlawPhotonFlux(void) { // Free members free_members(); @@ -223,7 +223,7 @@ GModelSpectralPlaw2::~GModelSpectralPlaw2(void) * @param[in] model Spectral power law model. * @return Spectral power law model. ***************************************************************************/ -GModelSpectralPlaw2& GModelSpectralPlaw2::operator=(const GModelSpectralPlaw2& model) +GModelSpectralPlawPhotonFlux& GModelSpectralPlawPhotonFlux::operator=(const GModelSpectralPlawPhotonFlux& model) { // Execute only if object is not identical if (this != &model) { @@ -256,7 +256,7 @@ GModelSpectralPlaw2& GModelSpectralPlaw2::operator=(const GModelSpectralPlaw2& m /***********************************************************************//** * @brief Clear power law model ***************************************************************************/ -void GModelSpectralPlaw2::clear(void) +void GModelSpectralPlawPhotonFlux::clear(void) { // Free class members (base and derived classes, derived class first) free_members(); @@ -276,10 +276,10 @@ void GModelSpectralPlaw2::clear(void) * * @return Pointer to deep copy of power law model ***************************************************************************/ -GModelSpectralPlaw2* GModelSpectralPlaw2::clone(void) const +GModelSpectralPlawPhotonFlux* GModelSpectralPlawPhotonFlux::clone(void) const { // Clone power law model - return new GModelSpectralPlaw2(*this); + return new GModelSpectralPlawPhotonFlux(*this); } @@ -293,7 +293,7 @@ GModelSpectralPlaw2* GModelSpectralPlaw2::clone(void) const * Computes * * \f[ - * S_{\rm E}(E | t) = {\tt m\_integral} + * S_{\rm E}(E | t) = {\tt m\_photon\_flux} * \frac{{\tt m\_index}+1} * {{\tt e\_max}^{{\tt m\_index}+1} - * {\tt e\_min}^{{\tt m\_index}+1}} @@ -304,7 +304,7 @@ GModelSpectralPlaw2* GModelSpectralPlaw2::clone(void) const * * \f[ * S_{\rm E}(E | t) = - * \frac{{\tt m\_integral}} + * \frac{{\tt m\_photon\_flux}} * {\log {\tt e\_max} - \log {\tt e\_min}} * E^{\tt m\_index} * \f] @@ -312,28 +312,28 @@ GModelSpectralPlaw2* GModelSpectralPlaw2::clone(void) const * for \f${\tt m\_index} = -1\f$, where * - \f${\tt e\_min}\f$ is the minimum energy of an interval, * - \f${\tt e\_max}\f$ is the maximum energy of an interval, - * - \f${\tt m\_integral}\f$ is the integral flux between + * - \f${\tt m\_photon\_flux}\f$ is the photon flux flux between * \f${\tt e\_min}\f$ and \f${\tt e\_max}\f$, and * - \f${\tt m\_index}\f$ is the spectral index. ***************************************************************************/ -double GModelSpectralPlaw2::eval(const GEnergy& srcEng, +double GModelSpectralPlawPhotonFlux::eval(const GEnergy& srcEng, const GTime& srcTime) const { // Update precomputed values update(srcEng); // Compute function value - double value = m_integral.value() * m_norm * m_power; + double value = m_photon_flux.value() * m_norm * m_power; // Compile option: Check for NaN/Inf #if defined(G_NAN_CHECK) if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: GModelSpectralPlaw2::eval"; + std::cout << "*** ERROR: GModelSpectralPlawPhotonFlux::eval"; std::cout << "(srcEng=" << srcEng; std::cout << ", srcTime=" << srcTime << "):"; std::cout << " NaN/Inf encountered"; std::cout << " (value=" << value; - std::cout << ", integral=" << integral(); + std::cout << ", photon_flux=" << photon_flux(); std::cout << ", m_norm=" << m_norm; std::cout << ", m_power=" << m_power; std::cout << ")" << std::endl; @@ -355,7 +355,7 @@ double GModelSpectralPlaw2::eval(const GEnergy& srcEng, * Computes * * \f[ - * S_{\rm E}(E | t) = {\tt m\_integral} + * S_{\rm E}(E | t) = {\tt m\_photon\_flux} * \frac{{\tt m\_index}+1} * {{\tt e\_max}^{{\tt m\_index}+1} - * {\tt e\_min}^{{\tt m\_index}+1}} @@ -366,7 +366,7 @@ double GModelSpectralPlaw2::eval(const GEnergy& srcEng, * * \f[ * S_{\rm E}(E | t) = - * \frac{{\tt m\_integral}} + * \frac{{\ttm\_photon\_flux}} * {\log {\tt e\_max} - \log {\tt e\_min}} * E^{\tt m\_index} * \f] @@ -374,7 +374,7 @@ double GModelSpectralPlaw2::eval(const GEnergy& srcEng, * for \f${\tt m\_index} = -1\f$, where * - \f${\tt e\_min}\f$ is the minimum energy of an interval, * - \f${\tt e\_max}\f$ is the maximum energy of an interval, - * - \f${\tt m\_integral}\f$ is the integral flux between + * - \f${\tt m\_photon\_flux}\f$ is the photon flux between * \f${\tt e\_min}\f$ and \f${\tt e\_max}\f$, and * - \f${\tt m\_index}\f$ is the spectral index. * @@ -382,8 +382,8 @@ double GModelSpectralPlaw2::eval(const GEnergy& srcEng, * respect to the parameters using * * \f[ - * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_integral}} = - * \frac{S_{\rm E}(E | t)}{{\tt m\_integral}} + * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_photon\_flux}} = + * \frac{S_{\rm E}(E | t)}{{\tt m\_photon\_flux}} * \f] * * \f[ @@ -408,22 +408,22 @@ double GModelSpectralPlaw2::eval(const GEnergy& srcEng, * * No partial derivatives are supported for the energy boundaries. ***************************************************************************/ -double GModelSpectralPlaw2::eval_gradients(const GEnergy& srcEng, +double GModelSpectralPlawPhotonFlux::eval_gradients(const GEnergy& srcEng, const GTime& srcTime) { // Initialise gradients - double g_integral = 0.0; + double g_photon_flux = 0.0; double g_index = 0.0; // Update precomputed values update(srcEng); // Compute function value - double value = m_integral.value() * m_norm * m_power; + double value = m_photon_flux.value() * m_norm * m_power; - // Integral flux gradient - if (m_integral.is_free() && m_integral.factor_value() > 0.0) { - g_integral = value / m_integral.factor_value(); + // Photon flux gradient + if (m_photon_flux.is_free() && m_photon_flux.factor_value() > 0.0) { + g_photon_flux = value / m_photon_flux.factor_value(); } // Index gradient @@ -433,21 +433,21 @@ double GModelSpectralPlaw2::eval_gradients(const GEnergy& srcEng, } // Set gradients - m_integral.factor_gradient(g_integral); + m_photon_flux.factor_gradient(g_photon_flux); m_index.factor_gradient(g_index); // Compile option: Check for NaN/Inf #if defined(G_NAN_CHECK) if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: GModelSpectralPlaw2::eval_gradients"; + std::cout << "*** ERROR: GModelSpectralPlawPhotonFlux::eval_gradients"; std::cout << "(srcEng=" << srcEng; std::cout << ", srcTime=" << srcTime << "):"; std::cout << " NaN/Inf encountered"; std::cout << " (value=" << value; - std::cout << ", integral=" << integral(); + std::cout << ", photon_flux=" << photon_flux(); std::cout << ", m_norm=" << m_norm; std::cout << ", m_power=" << m_power; - std::cout << ", g_integral=" << g_integral; + std::cout << ", g_photon_flux=" << g_photon_flux; std::cout << ", g_index=" << g_index; std::cout << ")" << std::endl; } @@ -475,7 +475,7 @@ double GModelSpectralPlaw2::eval_gradients(const GEnergy& srcEng, * - [@p emin, @p emax] is an energy interval, and * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV). ***************************************************************************/ -double GModelSpectralPlaw2::flux(const GEnergy& emin, +double GModelSpectralPlawPhotonFlux::flux(const GEnergy& emin, const GEnergy& emax) const { // Initialise flux @@ -493,7 +493,7 @@ double GModelSpectralPlaw2::flux(const GEnergy& emin, double pow_ref_emax = std::pow(this->emax().MeV(), gamma); double factor = (pow_emax - pow_emin) / (pow_ref_emax - pow_ref_emin); - flux = m_integral.value() * factor; + flux = m_photon_flux.value() * factor; } // Case B: Index is -1 @@ -504,7 +504,7 @@ double GModelSpectralPlaw2::flux(const GEnergy& emin, double log_ref_emax = std::log(this->emax().MeV()); double factor = (log_emax - log_emin) / (log_ref_emax - log_ref_emin); - flux = m_integral.value() * factor; + flux = m_photon_flux.value() * factor; } } // endif: integration range was valid @@ -531,7 +531,7 @@ double GModelSpectralPlaw2::flux(const GEnergy& emin, * - [@p emin, @p emax] is an energy interval, and * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV). ***************************************************************************/ -double GModelSpectralPlaw2::eflux(const GEnergy& emin, +double GModelSpectralPlawPhotonFlux::eflux(const GEnergy& emin, const GEnergy& emax) const { // Initialise flux @@ -546,12 +546,12 @@ double GModelSpectralPlaw2::eflux(const GEnergy& emin, double gamma = m_index.value() + 1.0; double pow_ref_emin = std::pow(this->emin().MeV(), gamma); double pow_ref_emax = std::pow(this->emax().MeV(), gamma); - norm = m_integral.value() * gamma / (pow_ref_emax - pow_ref_emin); + norm = m_photon_flux.value() * gamma / (pow_ref_emax - pow_ref_emin); } else { double log_ref_emin = std::log(this->emin().MeV()); double log_ref_emax = std::log(this->emax().MeV()); - norm = m_integral.value() / (log_ref_emax - log_ref_emin); + norm = m_photon_flux.value() / (log_ref_emax - log_ref_emin); } // Compute energy flux @@ -593,7 +593,7 @@ double GModelSpectralPlaw2::eflux(const GEnergy& emin, * * Returns Monte Carlo energy by randomly drawing from a power law. ***************************************************************************/ -GEnergy GModelSpectralPlaw2::mc(const GEnergy& emin, +GEnergy GModelSpectralPlawPhotonFlux::mc(const GEnergy& emin, const GEnergy& emax, const GTime& time, GRan& ran) const @@ -640,16 +640,16 @@ GEnergy GModelSpectralPlaw2::mc(const GEnergy& emin, * * Reads the spectral information from an XML element. ***************************************************************************/ -void GModelSpectralPlaw2::read(const GXmlElement& xml) +void GModelSpectralPlawPhotonFlux::read(const GXmlElement& xml) { // Get parameter pointers - const GXmlElement* integral = gammalib::xml_get_par(G_READ, xml, m_integral.name()); - const GXmlElement* index = gammalib::xml_get_par(G_READ, xml, m_index.name()); - const GXmlElement* emin = gammalib::xml_get_par(G_READ, xml, m_emin.name()); - const GXmlElement* emax = gammalib::xml_get_par(G_READ, xml, m_emax.name()); + const GXmlElement* phflux = gammalib::xml_get_par(G_READ, xml, m_photon_flux.name()); + const GXmlElement* index = gammalib::xml_get_par(G_READ, xml, m_index.name()); + const GXmlElement* emin = gammalib::xml_get_par(G_READ, xml, m_emin.name()); + const GXmlElement* emax = gammalib::xml_get_par(G_READ, xml, m_emax.name()); // Read parameters - m_integral.read(*integral); + m_photon_flux.read(*phflux); m_index.read(*index); m_emin.read(*emin); m_emax.read(*emax); @@ -665,11 +665,11 @@ void GModelSpectralPlaw2::read(const GXmlElement& xml) * @param[in] xml XML element. * * @exception GException::model_invalid_spectral - * Existing XML element is not of type "PowerLaw2" + * Existing XML element is not of type "PowerLaw" * * Writes the spectral information into an XML element. ***************************************************************************/ -void GModelSpectralPlaw2::write(GXmlElement& xml) const +void GModelSpectralPlawPhotonFlux::write(GXmlElement& xml) const { // Set model type if (xml.attribute("type") == "") { @@ -683,13 +683,13 @@ void GModelSpectralPlaw2::write(GXmlElement& xml) const } // Get XML parameters - GXmlElement* integral = gammalib::xml_need_par(G_WRITE, xml, m_integral.name()); - GXmlElement* index = gammalib::xml_need_par(G_WRITE, xml, m_index.name()); - GXmlElement* emin = gammalib::xml_need_par(G_WRITE, xml, m_emin.name()); - GXmlElement* emax = gammalib::xml_need_par(G_WRITE, xml, m_emax.name()); + GXmlElement* phflux = gammalib::xml_need_par(G_WRITE, xml, m_photon_flux.name()); + GXmlElement* index = gammalib::xml_need_par(G_WRITE, xml, m_index.name()); + GXmlElement* emin = gammalib::xml_need_par(G_WRITE, xml, m_emin.name()); + GXmlElement* emax = gammalib::xml_need_par(G_WRITE, xml, m_emax.name()); // Write parameters - m_integral.write(*integral); + m_photon_flux.write(*phflux); m_index.write(*index); m_emin.write(*emin); m_emax.write(*emax); @@ -705,7 +705,7 @@ void GModelSpectralPlaw2::write(GXmlElement& xml) const * @param[in] chatter Chattiness (defaults to NORMAL). * @return String containing power law information. ***************************************************************************/ -std::string GModelSpectralPlaw2::print(const GChatter& chatter) const +std::string GModelSpectralPlawPhotonFlux::print(const GChatter& chatter) const { // Initialise result string std::string result; @@ -714,7 +714,7 @@ std::string GModelSpectralPlaw2::print(const GChatter& chatter) const if (chatter != SILENT) { // Append header - result.append("=== GModelSpectralPlaw2 ==="); + result.append("=== GModelSpectralPlawPhotonFlux ==="); // Append information result.append("\n"+gammalib::parformat("Number of parameters")); @@ -739,21 +739,21 @@ std::string GModelSpectralPlaw2::print(const GChatter& chatter) const /***********************************************************************//** * @brief Initialise class members ***************************************************************************/ -void GModelSpectralPlaw2::init_members(void) +void GModelSpectralPlawPhotonFlux::init_members(void) { // Initialise model type m_type = "PowerLaw"; - // Initialise integral flux - m_integral.clear(); - m_integral.name("PhotonFlux"); - m_integral.unit("ph/cm2/s"); - m_integral.scale(1.0); - m_integral.value(1.0); // default: 1.0 - m_integral.range(0.0, 10.0); // range: [0,10] - m_integral.free(); - m_integral.gradient(0.0); - m_integral.has_grad(true); + // Initialise photon flux + m_photon_flux.clear(); + m_photon_flux.name("PhotonFlux"); + m_photon_flux.unit("ph/cm2/s"); + m_photon_flux.scale(1.0); + m_photon_flux.value(1.0); // default: 1.0 + m_photon_flux.range(0.0, 10.0); // range: [0,10] + m_photon_flux.free(); + m_photon_flux.gradient(0.0); + m_photon_flux.has_grad(true); // Initialise powerlaw index m_index.clear(); @@ -789,7 +789,7 @@ void GModelSpectralPlaw2::init_members(void) // Set parameter pointer(s) m_pars.clear(); - m_pars.push_back(&m_integral); + m_pars.push_back(&m_photon_flux); m_pars.push_back(&m_index); m_pars.push_back(&m_emin); m_pars.push_back(&m_emax); @@ -801,13 +801,13 @@ void GModelSpectralPlaw2::init_members(void) m_pow_emax = 0.0; m_norm = 0.0; m_power = 0.0; - m_last_integral = 0.0; + m_last_phflux = 0.0; m_last_index = 1000.0; m_last_emin.MeV(0.0); m_last_emax.MeV(0.0); m_last_energy.MeV(0.0); m_last_value = 0.0; - m_last_g_integral = 0.0; + m_last_g_phflux = 0.0; m_last_g_index = 0.0; // Return @@ -820,18 +820,18 @@ void GModelSpectralPlaw2::init_members(void) * * @param[in] model Spectral power law model. ***************************************************************************/ -void GModelSpectralPlaw2::copy_members(const GModelSpectralPlaw2& model) +void GModelSpectralPlawPhotonFlux::copy_members(const GModelSpectralPlawPhotonFlux& model) { // Copy members - m_type = model.m_type; - m_integral = model.m_integral; - m_index = model.m_index; - m_emin = model.m_emin; - m_emax = model.m_emax; + m_type = model.m_type; + m_photon_flux = model.m_photon_flux; + m_index = model.m_index; + m_emin = model.m_emin; + m_emax = model.m_emax; // Set parameter pointer(s) m_pars.clear(); - m_pars.push_back(&m_integral); + m_pars.push_back(&m_photon_flux); m_pars.push_back(&m_index); m_pars.push_back(&m_emin); m_pars.push_back(&m_emax); @@ -843,13 +843,13 @@ void GModelSpectralPlaw2::copy_members(const GModelSpectralPlaw2& model) m_pow_emax = model.m_pow_emax; m_norm = model.m_norm; m_power = model.m_power; - m_last_integral = model.m_last_integral; + m_last_phflux = model.m_last_phflux; m_last_index = model.m_last_index; m_last_emin = model.m_last_emin; m_last_emax = model.m_last_emax; m_last_energy = model.m_last_energy; m_last_value = model.m_last_value; - m_last_g_integral = model.m_last_g_integral; + m_last_g_phflux = model.m_last_g_phflux; m_last_g_index = model.m_last_g_index; // Return @@ -860,7 +860,7 @@ void GModelSpectralPlaw2::copy_members(const GModelSpectralPlaw2& model) /***********************************************************************//** * @brief Delete class members ***************************************************************************/ -void GModelSpectralPlaw2::free_members(void) +void GModelSpectralPlawPhotonFlux::free_members(void) { // Return return; @@ -872,7 +872,7 @@ void GModelSpectralPlaw2::free_members(void) * * @param[in] srcEng Source energy ***************************************************************************/ -void GModelSpectralPlaw2::update(const GEnergy& srcEng) const +void GModelSpectralPlawPhotonFlux::update(const GEnergy& srcEng) const { // Compute index+1 double gamma = index() + 1.0; diff --git a/src/model/Makefile.am b/src/model/Makefile.am index 211fc2b..4481ae6 100644 --- a/src/model/Makefile.am +++ b/src/model/Makefile.am @@ -56,7 +56,8 @@ sources = GModelPar.cpp \ GModelSpectralLogParabola.cpp \ GModelSpectralNodes.cpp \ GModelSpectralPlaw.cpp \ - GModelSpectralPlaw2.cpp \ + GModelSpectralPlawPhotonFlux.cpp \ + GModelSpectralPlawEnergyFlux.cpp \ GModelTemporal.cpp \ GModelTemporalRegistry.cpp \ GModelTemporalConst.cpp \ diff --git a/test/data/model_point_plaw_eflux.xml b/test/data/model_point_plaw_eflux.xml new file mode 100644 index 0000000..68df7c2 --- /dev/null +++ b/test/data/model_point_plaw_eflux.xml @@ -0,0 +1,15 @@ + + + + + + + + + + + + + + + diff --git a/test/data/model_point_plaw2.xml b/test/data/model_point_plaw_phflux.xml index c4c5e5c..c4c5e5c 100644 --- a/test/data/model_point_plaw2.xml +++ b/test/data/model_point_plaw_phflux.xml diff --git a/test/test_GModel.cpp b/test/test_GModel.cpp index 1f3ca89..a0faba1 100644 --- a/test/test_GModel.cpp +++ b/test/test_GModel.cpp @@ -58,7 +58,8 @@ void TestGModel::set(void) m_xml_model_point_const = datadir + "/model_point_const.xml"; m_xml_model_point_gauss = datadir + "/model_point_gauss.xml"; m_xml_model_point_plaw = datadir + "/model_point_plaw.xml"; - m_xml_model_point_plaw2 = datadir + "/model_point_plaw2.xml"; + m_xml_model_point_plaw_phflux = datadir + "/model_point_plaw_phflux.xml"; + m_xml_model_point_plaw_eflux = datadir + "/model_point_plaw_eflux.xml"; m_xml_model_point_eplaw = datadir + "/model_point_eplaw.xml"; m_xml_model_point_bplaw = datadir + "/model_point_bplaw.xml"; m_xml_model_point_supeplaw = datadir + "/model_point_supeplaw.xml"; @@ -124,8 +125,10 @@ void TestGModel::set(void) "Test GModelSpectralGauss"); append(static_cast(&TestGModel::test_plaw), "Test GModelSpectralPlaw"); - append(static_cast(&TestGModel::test_plaw2), - "Test GModelSpectralPlaw2"); + append(static_cast(&TestGModel::test_plaw_phflux), + "Test GModelSpectralPlawPhotonFkux"); + append(static_cast(&TestGModel::test_plaw_eflux), + "Test GModelSpectralPlawEnergyFlux"); append(static_cast(&TestGModel::test_eplaw), "Test GModelSpectralExpPlaw"); append(static_cast(&TestGModel::test_supeplaw), @@ -1551,35 +1554,35 @@ void TestGModel::test_plaw(void) /***********************************************************************//** - * @brief Test GModelSpectralPlaw2 class + * @brief Test GModelSpectralPlawPhotonFlux class ***************************************************************************/ -void TestGModel::test_plaw2(void) +void TestGModel::test_plaw_phflux(void) { // Test void constructor - GModelSpectralPlaw2 model1; + GModelSpectralPlawPhotonFlux model1; test_value(model1.type(), "PowerLaw", "Check type of void model"); // Test value constructor - GModelSpectralPlaw2 model2(2.0, -2.1, GEnergy(10.0, "MeV"), GEnergy(100.0, "MeV")); - test_value(model2.integral(), 2.0); + GModelSpectralPlawPhotonFlux model2(2.0, -2.1, GEnergy(10.0, "MeV"), GEnergy(100.0, "MeV")); + test_value(model2.photon_flux(), 2.0); test_value(model2.index(), -2.1); test_value(model2.emin().MeV(), 10.0); test_value(model2.emax().MeV(), 100.0); - + // Test XML constructor and value - GXml xml(m_xml_model_point_plaw2); + GXml xml(m_xml_model_point_plaw_phflux); GXmlElement* element = xml.element(0)->element(0)->element("spectrum", 0); - GModelSpectralPlaw2 model3(*element); + GModelSpectralPlawPhotonFlux model3(*element); test_value(model3.size(), 4); test_value(model3.type(), "PowerLaw", "Check model type"); - test_value(model3.integral(), 1.0e-7); + test_value(model3.photon_flux(), 1.0e-7); test_value(model3.index(), -2.0); test_value(model3.emin().MeV(), 100.0); test_value(model3.emax().MeV(), 500000.0); // Test integral method - model3.integral(2.1e-7); - test_value(model3.integral(), 2.1e-7); + model3.photon_flux(2.1e-7); + test_value(model3.photon_flux(), 2.1e-7); // Test index method model3.index(-2.3); @@ -1612,6 +1615,67 @@ void TestGModel::test_plaw2(void) /***********************************************************************//** + * @brief Test GModelSpectralPlawEnergyFlux class + ***************************************************************************/ +void TestGModel::test_plaw_eflux(void) +{ + // Test void constructor + GModelSpectralPlawEnergyFlux model1; + test_value(model1.type(), "PowerLaw", "Check type of void model"); + + // Test value constructor + GModelSpectralPlawEnergyFlux model2(2.0, -2.1, GEnergy(10.0, "MeV"), GEnergy(100.0, "MeV")); + test_value(model2.energy_flux(), 2.0); + test_value(model2.index(), -2.1); + test_value(model2.emin().MeV(), 10.0); + test_value(model2.emax().MeV(), 100.0); + + // Test XML constructor and value + GXml xml(m_xml_model_point_plaw_eflux); + GXmlElement* element = xml.element(0)->element(0)->element("spectrum", 0); + GModelSpectralPlawEnergyFlux model3(*element); + test_value(model3.size(), 4); + test_value(model3.type(), "PowerLaw", "Check model type"); + test_value(model3.energy_flux(), 1.0e-7); + test_value(model3.index(), -2.0); + test_value(model3.emin().MeV(), 100.0); + test_value(model3.emax().MeV(), 500000.0); + + // Test integral method + model3.energy_flux(2.1e-7); + test_value(model3.energy_flux(), 2.1e-7); + + // Test index method + model3.index(-2.3); + test_value(model3.index(), -2.3); + + // Test emin method + model3.emin(GEnergy(10.0, "MeV")); + test_value(model3.emin().MeV(), 10.0); + + // Test emax method + model3.emax(GEnergy(10.0, "MeV")); + test_value(model3.emax().MeV(), 10.0); + + // Test operator access + const char* strarray[] = {"EnergyFlux", "Index", "LowerLimit", "UpperLimit"}; + for (int i = 0; i < 4; ++i) { + std::string keyname(strarray[i]); + model3[keyname].remove_range(); // To allow setting of any value + model3[keyname].value(2.1); + model3[keyname].error(1.9); + model3[keyname].gradient(0.8); + test_value(model3[keyname].value(), 2.1); + test_value(model3[keyname].error(), 1.9); + test_value(model3[keyname].gradient(), 0.8); + } + + // Exit test + return; +} + + +/***********************************************************************//** * @brief Test GModelSpectralExpPlaw class ***************************************************************************/ void TestGModel::test_eplaw(void) @@ -2094,15 +2158,16 @@ void TestGModel::test_spatial_model(void) void TestGModel::test_spectral_model(void) { // Test spectral models XML interface - test_xml_model("GModelSpectralConst", m_xml_model_point_const); - 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("GModelSpectralBrokenPlaw", m_xml_model_point_bplaw); - test_xml_model("GModelSpectralSuperExpPlaw", m_xml_model_point_supeplaw); - test_xml_model("GModelSpectralLogParabola", m_xml_model_point_logparabola); - test_xml_model("GModelSpectralNodes", m_xml_model_point_nodes); - test_xml_model("GModelSpectralFunc", m_xml_model_point_filefct); + test_xml_model("GModelSpectralConst", m_xml_model_point_const); + test_xml_model("GModelSpectralPlaw", m_xml_model_point_plaw); + test_xml_model("GModelSpectralPlawPhotonFlux", m_xml_model_point_plaw_phflux); + test_xml_model("GModelSpectralPlawEnergyFlux", m_xml_model_point_plaw_eflux); + test_xml_model("GModelSpectralExpPaw", m_xml_model_point_eplaw); + 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); + test_xml_model("GModelSpectralNodes", m_xml_model_point_nodes); + test_xml_model("GModelSpectralFunc", m_xml_model_point_filefct); // Return return; @@ -2823,13 +2888,13 @@ void TestGModel::test_legacy_model_point_plaw2(void) GModels models(m_xml_legacy_point_plaw2); // Extract spectral component - GModelSpectralPlaw2* spectral = static_cast( + GModelSpectralPlawPhotonFlux* spectral = static_cast( static_cast( models[0])->spectral()); // Test model values test_value(spectral->type(), "PowerLaw2", "Check model type"); - test_value(spectral->integral(), 1.0e-7, 1.0e-7, "Check integral photon flux"); + test_value(spectral->photon_flux(), 1.0e-7, 1.0e-7, "Check integral photon flux"); test_value(spectral->index(), -2.0, 1.0e-7, "Check index"); test_value(spectral->emin().MeV(), 100.0, 1.0e-7, "Check minimum energy"); test_value(spectral->emax().MeV(), 500000.0, 1.0e-7, "Check maximum energy"); @@ -2840,12 +2905,12 @@ void TestGModel::test_legacy_model_point_plaw2(void) models.load("test_xml_legacy_point_plaw2.xml"); // Extract spectral component - spectral = static_cast( + spectral = static_cast( static_cast(models[0])->spectral()); // Re-test model values test_value(spectral->type(), "PowerLaw2", "Check model type"); - test_value(spectral->integral(), 1.0e-7, 1.0e-7, "Check integral photon flux"); + test_value(spectral->photon_flux(), 1.0e-7, 1.0e-7, "Check integral photon flux"); test_value(spectral->index(), -2.0, 1.0e-7, "Check index"); test_value(spectral->emin().MeV(), 100.0, 1.0e-7, "Check minimum energy"); test_value(spectral->emax().MeV(), 500000.0, 1.0e-7, "Check maximum energy"); diff --git a/test/test_GModel.hpp b/test/test_GModel.hpp index 1833271..0c45428 100644 --- a/test/test_GModel.hpp +++ b/test/test_GModel.hpp @@ -62,7 +62,8 @@ public: void test_const(void); void test_gauss(void); void test_plaw(void); - void test_plaw2(void); + void test_plaw_phflux(void); + void test_plaw_eflux(void); void test_eplaw(void); void test_bplaw(void); void test_supeplaw(void); @@ -102,7 +103,8 @@ private: std::string m_xml_model_point_const; std::string m_xml_model_point_gauss; std::string m_xml_model_point_plaw; - std::string m_xml_model_point_plaw2; + std::string m_xml_model_point_plaw_phflux; + std::string m_xml_model_point_plaw_eflux; std::string m_xml_model_point_eplaw; std::string m_xml_model_point_bplaw; std::string m_xml_model_point_supeplaw;