Rolf Buehler / gammalib

Commit 938b5ea0397ab19dad3a396bdff22ea5803d3067

Authored by Jürgen Knödlseder
1 parent a6d78bad

Add time dependence to eval(), eval_gradients() and mc() methods. 

The eval_gradients() method is no longer const to reflect the fact that the method explicitly sets the gradient attribute of the model parameters.
Inline Side-by-side
Showing 23 changed files with 609 additions and 235 deletions
include/GModelSpectral.hpp
  View file @938b5ea
... ... @@ -33,6 +33,7 @@
33 33 #include "GBase.hpp"
34 34 #include "GModelPar.hpp"
35 35 #include "GEnergy.hpp"
  36 +#include "GTime.hpp"
36 37 #include "GRan.hpp"
37 38 #include "GXmlElement.hpp"
38 39  
... ... @@ -42,8 +43,25 @@
42 43 *
43 44 * @brief Abstract spectral model base class
44 45 *
45   - * This class implements the spectral component of the factorised source
46   - * model.
  46 + * This class implements the spectral component of the factorized source
  47 + * model
  48 + *
  49 + * \f[
  50 + * S_{\rm E}(E | t)
  51 + * \f]
  52 + *
  53 + * where
  54 + * - \f$E\f$ is the true photon energy, and
  55 + * - \f$t\f$ is the true photon arrival time.
  56 + *
  57 + * The spectral component describes the spatially integrated time dependent
  58 + * spectral distribution of the source. It satisfies
  59 + * \f[
  60 + * \int_{E} S_{\rm E}(E | t) dE = \Phi
  61 + * \f]
  62 + * for all \f$t\f$, where \f$\Phi\f$ is the spatially and spectrally
  63 + * integrated total source flux. The spectral component does not impact
  64 + * the temporal properties of the integrated flux \f$\Phi\f$.
47 65 ***************************************************************************/
48 66 class GModelSpectral : public GBase {
49 67  
... ... @@ -64,11 +82,16 @@ public:
64 82 virtual void clear(void) = 0;
65 83 virtual GModelSpectral* clone(void) const = 0;
66 84 virtual std::string type(void) const = 0;
67   - virtual double eval(const GEnergy& srcEng) const = 0;
68   - virtual double eval_gradients(const GEnergy& srcEng) const = 0;
69   - virtual double flux(const GEnergy& emin, const GEnergy& emax) const = 0;
70   - virtual double eflux(const GEnergy& emin, const GEnergy& emax) const = 0;
71   - virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const = 0;
  85 + virtual double eval(const GEnergy& srcEng,
  86 + const GTime& srcTime) const = 0;
  87 + virtual double eval_gradients(const GEnergy& srcEng,
  88 + const GTime& srcTime) = 0;
  89 + virtual double flux(const GEnergy& emin,
  90 + const GEnergy& emax) const = 0;
  91 + virtual double eflux(const GEnergy& emin,
  92 + const GEnergy& emax) const = 0;
  93 + virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax,
  94 + const GTime& time, GRan& ran) const = 0;
72 95 virtual void read(const GXmlElement& xml) = 0;
73 96 virtual void write(GXmlElement& xml) const = 0;
74 97 virtual std::string print(void) const = 0;
... ...
include/GModelSpectralConst.hpp
  View file @938b5ea
... ... @@ -40,11 +40,15 @@
40 40 *
41 41 * @brief Constant spectral model class
42 42 *
43   - * This class implements a constant as the spectral component of the
44   - * gamma-ray sky model. The function is defined as
45   - * \f[I(E)=norm\f]
  43 + * This class implements a constant spectrum. The model is defined by
  44 + *
  45 + * \f[
  46 + * S_{\rm E}(E | t) = {\tt m\_norm}
  47 + * \f]
  48 + *
46 49 * where
47   - * \f$norm\f$ is the normalization constant.
  50 + * \f${\tt m\_norm}\f$ is the normalization constant in units of
  51 + * ph/cm2/s/MeV.
48 52 ***************************************************************************/
49 53 class GModelSpectralConst : public GModelSpectral {
50 54  
... ... @@ -52,27 +56,36 @@ public:
52 56 // Constructors and destructors
53 57 GModelSpectralConst(void);
54 58 explicit GModelSpectralConst(const GXmlElement& xml);
  59 + explicit GModelSpectralConst(const double& norm);
55 60 GModelSpectralConst(const GModelSpectralConst& model);
56 61 virtual ~GModelSpectralConst(void);
57 62  
58 63 // Operators
59 64 virtual GModelSpectralConst& operator=(const GModelSpectralConst& model);
60 65  
61   - // Implemented pure virtual methods
  66 + // Implemented pure virtual base class methods
62 67 virtual void clear(void);
63 68 virtual GModelSpectralConst* clone(void) const;
64   - virtual std::string type(void) const { return "ConstantValue"; }
65   - virtual double eval(const GEnergy& srcEng) const;
66   - virtual double eval_gradients(const GEnergy& srcEng) const;
67   - virtual double flux(const GEnergy& emin, const GEnergy& emax) const;
68   - virtual double eflux(const GEnergy& emin, const GEnergy& emax) const;
69   - virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
  69 + virtual std::string type(void) const;
  70 + virtual double eval(const GEnergy& srcEng,
  71 + const GTime& srcTime) const;
  72 + virtual double eval_gradients(const GEnergy& srcEng,
  73 + const GTime& srcTime);
  74 + virtual double flux(const GEnergy& emin,
  75 + const GEnergy& emax) const;
  76 + virtual double eflux(const GEnergy& emin,
  77 + const GEnergy& emax) const;
  78 + virtual GEnergy mc(const GEnergy& emin,
  79 + const GEnergy& emax,
  80 + const GTime& time,
  81 + GRan& ran) const;
70 82 virtual void read(const GXmlElement& xml);
71 83 virtual void write(GXmlElement& xml) const;
72 84 virtual std::string print(void) const;
73 85  
74 86 // Other methods
75   - double norm(void) const { return m_norm.value(); }
  87 + double norm(void) const;
  88 + void norm(const double& norm);
76 89  
77 90 protected:
78 91 // Protected methods
... ... @@ -81,7 +94,50 @@ protected:
81 94 void free_members(void);
82 95  
83 96 // Protected members
84   - GModelPar m_norm; //!< Normalization factor
  97 + GModelPar m_norm; //!< Normalization factor
85 98 };
86 99  
  100 +
  101 +/***********************************************************************//**
  102 + * @brief Return model type
  103 + *
  104 + * @return "ConstantValue".
  105 + *
  106 + * Returns the type of the constant spectral model.
  107 + ***************************************************************************/
  108 +inline
  109 +std::string GModelSpectralConst::type(void) const
  110 +{
  111 + return "ConstantValue";
  112 +}
  113 +
  114 +
  115 +/***********************************************************************//**
  116 + * @brief Return normalization factor
  117 + *
  118 + * @return Normalization factor (ph/cm2/s/MeV).
  119 + *
  120 + * Returns the normalization factor.
  121 + ***************************************************************************/
  122 +inline
  123 +double GModelSpectralConst::norm(void) const
  124 +{
  125 + return (m_norm.value());
  126 +}
  127 +
  128 +
  129 +/***********************************************************************//**
  130 + * @brief Set normalization factor
  131 + *
  132 + * @param[in] norm Normalization factor (ph/cm2/s/MeV).
  133 + *
  134 + * Sets the normalization factor.
  135 + ***************************************************************************/
  136 +inline
  137 +void GModelSpectralConst::norm(const double& norm)
  138 +{
  139 + m_norm.value(norm);
  140 + return;
  141 +}
  142 +
87 143 #endif /* GMODELSPECTRALCONST_HPP */
... ...
include/GModelSpectralExpPlaw.hpp
  View file @938b5ea
... ... @@ -41,37 +41,51 @@
41 41 *
42 42 * @brief Exponential cut off power law spectral class
43 43 *
44   - * This class implements a power law as the spectral component of the
45   - * gamma-ray sky model. The power law is defined as
46   - * \f[I(E)=norm (E/pivot)^{index} \exp(-E/ecut)\f]
  44 + * This class implements a power law spectrum with exponential cut off. The
  45 + * model is defined by
  46 + *
  47 + * \f[
  48 + * S_{\rm E}(E | t) = {\tt m\_norm}
  49 + * \left( \frac{E}{\tt m\_pivot} \right)^{\tt m\_index}
  50 + * \exp \left( \frac{-E}{\tt m\_ecut} \right)
  51 + * \f]
  52 + *
47 53 * where
48   - * \f$norm\f$ is the normalization or prefactor,
49   - * \f$pivot\f$ is the pivot energy,
50   - * \f$index\f$ is the spectral index, and
51   - * \f$ecut\f$ is the cut off energy.
  54 + * - \f${\tt m\_norm}\f$ is the normalization or prefactor,
  55 + * - \f${\tt m\_pivot}\f$ is the pivot energy,
  56 + * - \f${\tt m\_index}\f$ is the spectral index, and
  57 + * - \f${\tt m\_ecut}\f$ is the cut off energy.
52 58 ***************************************************************************/
53 59 class GModelSpectralExpPlaw : public GModelSpectral {
54 60  
55 61 public:
56 62 // Constructors and destructors
57 63 GModelSpectralExpPlaw(void);
58   - explicit GModelSpectralExpPlaw(const double& norm, const double& index, const double& ecut);
  64 + explicit GModelSpectralExpPlaw(const double& norm, const double& index,
  65 + const double& ecut);
59 66 explicit GModelSpectralExpPlaw(const GXmlElement& xml);
60 67 GModelSpectralExpPlaw(const GModelSpectralExpPlaw& model);
61 68 virtual ~GModelSpectralExpPlaw(void);
62 69  
63 70 // Operators
64   - virtual GModelSpectralExpPlaw& operator= (const GModelSpectralExpPlaw& model);
  71 + virtual GModelSpectralExpPlaw& operator=(const GModelSpectralExpPlaw& model);
65 72  
66 73 // Implemented pure virtual methods
67 74 virtual void clear(void);
68 75 virtual GModelSpectralExpPlaw* clone(void) const;
69   - virtual std::string type(void) const { return "ExpCutoff"; }
70   - virtual double eval(const GEnergy& srcEng) const;
71   - virtual double eval_gradients(const GEnergy& srcEng) const;
72   - virtual double flux(const GEnergy& emin, const GEnergy& emax) const;
73   - virtual double eflux(const GEnergy& emin, const GEnergy& emax) const;
74   - virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
  76 + virtual std::string type(void) const;
  77 + virtual double eval(const GEnergy& srcEng,
  78 + const GTime& srcTime) const;
  79 + virtual double eval_gradients(const GEnergy& srcEng,
  80 + const GTime& srcTime);
  81 + virtual double flux(const GEnergy& emin,
  82 + const GEnergy& emax) const;
  83 + virtual double eflux(const GEnergy& emin,
  84 + const GEnergy& emax) const;
  85 + virtual GEnergy mc(const GEnergy& emin,
  86 + const GEnergy& emax,
  87 + const GTime& time,
  88 + GRan& ran) const;
75 89 virtual void read(const GXmlElement& xml);
76 90 virtual void write(GXmlElement& xml) const;
77 91 virtual std::string print(void) const;
... ... @@ -141,4 +155,18 @@ protected:
141 155 mutable double m_mc_pow_ewidth; //!< Power of energy width
142 156 };
143 157  
  158 +
  159 +/***********************************************************************//**
  160 + * @brief Return model type
  161 + *
  162 + * @return "ExpCutoff".
  163 + *
  164 + * Returns the type of the exponentially cut off power law model.
  165 + ***************************************************************************/
  166 +inline
  167 +std::string GModelSpectralExpPlaw::type(void) const
  168 +{
  169 + return "ExpCutoff";
  170 +}
  171 +
144 172 #endif /* GMODELSPECTRALEXPPLAW_HPP */
... ...
include/GModelSpectralFunc.hpp
  View file @938b5ea
... ... @@ -21,7 +21,7 @@
21 21 /**
22 22 * @file GModelSpectralFunc.hpp
23 23 * @brief Spectral function model class definition
24   - * @author J. Knoedlseder
  24 + * @author Juergen Knoedlseder
25 25 */
26 26  
27 27 #ifndef GMODELSPECTRALFUNC_HPP
... ... @@ -61,15 +61,20 @@ public:
61 61 // Operators
62 62 virtual GModelSpectralFunc& operator=(const GModelSpectralFunc& model);
63 63  
64   - // Implemented pure virtual methods
  64 + // Implemented pure virtual base class methods
65 65 virtual void clear(void);
66 66 virtual GModelSpectralFunc* clone(void) const;
67   - virtual std::string type(void) const { return "FileFunction"; }
68   - virtual double eval(const GEnergy& srcEng) const;
69   - virtual double eval_gradients(const GEnergy& srcEng) const;
70   - virtual double flux(const GEnergy& emin, const GEnergy& emax) const;
71   - virtual double eflux(const GEnergy& emin, const GEnergy& emax) const;
72   - virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
  67 + virtual std::string type(void) const;
  68 + virtual double eval(const GEnergy& srcEng,
  69 + const GTime& srcTime) const;
  70 + virtual double eval_gradients(const GEnergy& srcEng,
  71 + const GTime& srcTime);
  72 + virtual double flux(const GEnergy& emin,
  73 + const GEnergy& emax) const;
  74 + virtual double eflux(const GEnergy& emin,
  75 + const GEnergy& emax) const;
  76 + virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax,
  77 + const GTime& time, GRan& ran) const;
73 78 virtual void read(const GXmlElement& xml);
74 79 virtual void write(GXmlElement& xml) const;
75 80 virtual std::string print(void) const;
... ... @@ -110,4 +115,18 @@ protected:
110 115 mutable std::vector<double> m_mc_exp; //!< Exponent for MC
111 116 };
112 117  
  118 +
  119 +/***********************************************************************//**
  120 + * @brief Return model type
  121 + *
  122 + * @return "FileFunction".
  123 + *
  124 + * Returns the type of the spectral function model.
  125 + ***************************************************************************/
  126 +inline
  127 +std::string GModelSpectralFunc::type(void) const
  128 +{
  129 + return "FileFunction";
  130 +}
  131 +
113 132 #endif /* GMODELSPECTRALFUNC_HPP */
... ...
include/GModelSpectralLogParabola.hpp
  View file @938b5ea
... ... @@ -60,7 +60,7 @@ public:
60 60 GModelSpectralLogParabola(void);
61 61 explicit GModelSpectralLogParabola(const double& norm,
62 62 const double& index,
63   - const double &curvature);
  63 + const double& curvature);
64 64 explicit GModelSpectralLogParabola(const GXmlElement& xml);
65 65 GModelSpectralLogParabola(const GModelSpectralLogParabola& model);
66 66 virtual ~GModelSpectralLogParabola(void);
... ... @@ -71,15 +71,19 @@ public:
71 71 // Implemented pure virtual methods
72 72 virtual void clear(void);
73 73 virtual GModelSpectralLogParabola* clone(void) const;
74   - virtual std::string type(void) const { return "LogParabola"; }
75   - virtual double eval(const GEnergy& srcEng) const;
76   - virtual double eval_gradients(const GEnergy& srcEng) const;
  74 + virtual std::string type(void) const;
  75 + virtual double eval(const GEnergy& srcEng,
  76 + const GTime& srcTime) const;
  77 + virtual double eval_gradients(const GEnergy& srcEng,
  78 + const GTime& srcTime);
77 79 virtual double flux(const GEnergy& emin,
78 80 const GEnergy& emax) const;
79 81 virtual double eflux(const GEnergy& emin,
80 82 const GEnergy& emax) const;
81 83 virtual GEnergy mc(const GEnergy& emin,
82   - const GEnergy& emax, GRan& ran) const;
  84 + const GEnergy& emax,
  85 + const GTime& time,
  86 + GRan& ran) const;
83 87 virtual void read(const GXmlElement& xml);
84 88 virtual void write(GXmlElement& xml) const;
85 89 virtual std::string print(void) const;
... ... @@ -95,7 +99,8 @@ protected:
95 99 void init_members(void);
96 100 void copy_members(const GModelSpectralLogParabola& model);
97 101 void free_members(void);
98   - void update_mc_cache(const GEnergy& emin, const GEnergy& emax) const;
  102 + void update_mc_cache(const GEnergy& emin, const GEnergy& emax,
  103 + const GTime& time) const;
99 104  
100 105 // Class to determine to the integral photon flux
101 106 class flux_kern : public GFunction {
... ... @@ -146,4 +151,18 @@ protected:
146 151 mutable double m_mc_norm; //!< Norm of powerlaw model at logparabola pivot energy
147 152 };
148 153  
  154 +
  155 +/***********************************************************************//**
  156 + * @brief Return model type
  157 + *
  158 + * @return "LogParabola".
  159 + *
  160 + * Returns the type of the log parabola spectral model.
  161 + ***************************************************************************/
  162 +inline
  163 +std::string GModelSpectralLogParabola::type(void) const
  164 +{
  165 + return "LogParabola";
  166 +}
  167 +
149 168 #endif /* GMODELSPECTRALLOGPARABOLA_HPP */
... ...
include/GModelSpectralNodes.hpp
  View file @938b5ea
1 1 /***************************************************************************
2   - * GModelSpectralNodes.hpp - Spectral nodes model class *
  2 + * GModelSpectralNodes.hpp - Spectral nodes model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2013 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -21,7 +21,7 @@
21 21 /**
22 22 * @file GModelSpectralNodes.hpp
23 23 * @brief Spectral nodes model class definition
24   - * @author J. Knoedlseder
  24 + * @author Juergen Knoedlseder
25 25 */
26 26  
27 27 #ifndef GMODELSPECTRALNODES_HPP
... ... @@ -76,12 +76,19 @@ public:
76 76 // Implemented pure virtual methods
77 77 virtual void clear(void);
78 78 virtual GModelSpectralNodes* clone(void) const;
79   - virtual std::string type(void) const { return "NodeFunction"; }
80   - virtual double eval(const GEnergy& srcEng) const;
81   - virtual double eval_gradients(const GEnergy& srcEng) const;
82   - virtual double flux(const GEnergy& emin, const GEnergy& emax) const;
83   - virtual double eflux(const GEnergy& emin, const GEnergy& emax) const;
84   - virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
  79 + virtual std::string type(void) const;
  80 + virtual double eval(const GEnergy& srcEng,
  81 + const GTime& srcTime) const;
  82 + virtual double eval_gradients(const GEnergy& srcEng,
  83 + const GTime& srcTime);
  84 + virtual double flux(const GEnergy& emin,
  85 + const GEnergy& emax) const;
  86 + virtual double eflux(const GEnergy& emin,
  87 + const GEnergy& emax) const;
  88 + virtual GEnergy mc(const GEnergy& emin,
  89 + const GEnergy& emax,
  90 + const GTime& time,
  91 + GRan& ran) const;
85 92 virtual void read(const GXmlElement& xml);
86 93 virtual void write(GXmlElement& xml) const;
87 94 virtual std::string print(void) const;
... ... @@ -128,4 +135,18 @@ protected:
128 135 mutable std::vector<double> m_mc_exp; //!< Exponent for MC
129 136 };
130 137  
  138 +
  139 +/***********************************************************************//**
  140 + * @brief Return model type
  141 + *
  142 + * @return "NodeFunction".
  143 + *
  144 + * Returns the type of the spectral node function model.
  145 + ***************************************************************************/
  146 +inline
  147 +std::string GModelSpectralNodes::type(void) const
  148 +{
  149 + return "NodeFunction";
  150 +}
  151 +
131 152 #endif /* GMODELSPECTRALNODES_HPP */
... ...
include/GModelSpectralPlaw.hpp
  View file @938b5ea
1 1 /***************************************************************************
2 2 * GModelSpectralPlaw.hpp - Spectral power law model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2013 by Jurgen Knodlseder *
  4 + * copyright (C) 2009-2013 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -65,12 +65,19 @@ public:
65 65 // Implemented pure virtual methods
66 66 virtual void clear(void);
67 67 virtual GModelSpectralPlaw* clone(void) const;
68   - virtual std::string type(void) const { return "PowerLaw"; }
69   - virtual double eval(const GEnergy& srcEng) const;
70   - virtual double eval_gradients(const GEnergy& srcEng) const;
71   - virtual double flux(const GEnergy& emin, const GEnergy& emax) const;
72   - virtual double eflux(const GEnergy& emin, const GEnergy& emax) const;
73   - virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
  68 + virtual std::string type(void) const;
  69 + virtual double eval(const GEnergy& srcEng,
  70 + const GTime& srcTime) const;
  71 + virtual double eval_gradients(const GEnergy& srcEng,
  72 + const GTime& srcTime);
  73 + virtual double flux(const GEnergy& emin,
  74 + const GEnergy& emax) const;
  75 + virtual double eflux(const GEnergy& emin,
  76 + const GEnergy& emax) const;
  77 + virtual GEnergy mc(const GEnergy& emin,
  78 + const GEnergy& emax,
  79 + const GTime& time,
  80 + GRan& ran) const;
74 81 virtual void read(const GXmlElement& xml);
75 82 virtual void write(GXmlElement& xml) const;
76 83 virtual std::string print(void) const;
... ... @@ -92,4 +99,18 @@ protected:
92 99 GModelPar m_pivot; //!< Pivot energy
93 100 };
94 101  
  102 +
  103 +/***********************************************************************//**
  104 + * @brief Return model type
  105 + *
  106 + * @return "PowerLaw".
  107 + *
  108 + * Returns the type of the spectral power law model.
  109 + ***************************************************************************/
  110 +inline
  111 +std::string GModelSpectralPlaw::type(void) const
  112 +{
  113 + return "PowerLaw";
  114 +}
  115 +
95 116 #endif /* GMODELSPECTRALPLAW_HPP */
... ...
include/GModelSpectralPlaw2.hpp
  View file @938b5ea
... ... @@ -65,15 +65,22 @@ public:
65 65 // Operators
66 66 virtual GModelSpectralPlaw2& operator= (const GModelSpectralPlaw2& model);
67 67  
68   - // Implemented pure virtual methods
  68 + // Implemented pure virtual base class methods
69 69 virtual void clear(void);
70 70 virtual GModelSpectralPlaw2* clone(void) const;
71   - virtual std::string type(void) const { return "PowerLaw2"; }
72   - virtual double eval(const GEnergy& srcEng) const;
73   - virtual double eval_gradients(const GEnergy& srcEng) const;
74   - virtual double flux(const GEnergy& emin, const GEnergy& emax) const;
75   - virtual double eflux(const GEnergy& emin, const GEnergy& emax) const;
76   - virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
  71 + virtual std::string type(void) const;
  72 + virtual double eval(const GEnergy& srcEng,
  73 + const GTime& srcTime) const;
  74 + virtual double eval_gradients(const GEnergy& srcEng,
  75 + const GTime& srcTime);
  76 + virtual double flux(const GEnergy& emin,
  77 + const GEnergy& emax) const;
  78 + virtual double eflux(const GEnergy& emin,
  79 + const GEnergy& emax) const;
  80 + virtual GEnergy mc(const GEnergy& emin,
  81 + const GEnergy& emax,
  82 + const GTime& time,
  83 + GRan& ran) const;
77 84 virtual void read(const GXmlElement& xml);
78 85 virtual void write(GXmlElement& xml) const;
79 86 virtual std::string print(void) const;
... ... @@ -115,4 +122,18 @@ protected:
115 122 mutable double m_last_g_index; //!< Last spectral index gradient
116 123 };
117 124  
  125 +
  126 +/***********************************************************************//**
  127 + * @brief Return model type
  128 + *
  129 + * @return "PowerLaw2".
  130 + *
  131 + * Returns the type of the spectral power law model.
  132 + ***************************************************************************/
  133 +inline
  134 +std::string GModelSpectralPlaw2::type(void) const
  135 +{
  136 + return "PowerLaw2";
  137 +}
  138 +
118 139 #endif /* GMODELSPECTRALPLAW2_HPP */
... ...
pyext/GModelSpectral.i
  View file @938b5ea
... ... @@ -47,11 +47,16 @@ public:
47 47 virtual void clear(void) = 0;
48 48 virtual GModelSpectral* clone(void) const = 0;
49 49 virtual std::string type(void) const = 0;
50   - virtual double eval(const GEnergy& srcEng) const = 0;
51   - virtual double eval_gradients(const GEnergy& srcEng) const = 0;
52   - virtual double flux(const GEnergy& emin, const GEnergy& emax) const = 0;
53   - virtual double eflux(const GEnergy& emin, const GEnergy& emax) const = 0;
54   - virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const = 0;
  50 + virtual double eval(const GEnergy& srcEng,
  51 + const GTime& srcTime) const = 0;
  52 + virtual double eval_gradients(const GEnergy& srcEng,
  53 + const GTime& srcTime) = 0;
  54 + virtual double flux(const GEnergy& emin,
  55 + const GEnergy& emax) const = 0;
  56 + virtual double eflux(const GEnergy& emin,
  57 + const GEnergy& emax) const = 0;
  58 + virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax,
  59 + const GTime& time, GRan& ran) const = 0;
55 60 virtual void read(const GXmlElement& xml) = 0;
56 61 virtual void write(GXmlElement& xml) const = 0;
57 62  
... ... @@ -69,10 +74,12 @@ public:
69 74 return tochar(self->print());
70 75 }
71 76 GModelPar& __getitem__(const int& index) {
72   - if (index >= 0 && index < self->size())
  77 + if (index >= 0 && index < self->size()) {
73 78 return (*self)[index];
74   - else
  79 + }
  80 + else {
75 81 throw GException::out_of_range("__getitem__(int)", index, self->size());
  82 + }
76 83 }
77 84 GModelPar& __getitem__(const std::string& name) {
78 85 return (*self)[name];
... ... @@ -82,8 +89,9 @@ public:
82 89 (*self)[index] = val;
83 90 return;
84 91 }
85   - else
  92 + else {
86 93 throw GException::out_of_range("__setitem__(int)", index, self->size());
  94 + }
87 95 }
88 96 void __setitem__(const std::string& name, const GModelPar& val) {
89 97 (*self)[name] = val;
... ...
pyext/GModelSpectralConst.i
  View file @938b5ea
1 1 /***************************************************************************
2   - * GModelSpectralConst.i - Spectral constant model class *
  2 + * GModelSpectralConst.i - Spectral constant model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011 by Jurgen Knodlseder *
  4 + * copyright (C) 2011-2013 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -20,8 +20,8 @@
20 20 ***************************************************************************/
21 21 /**
22 22 * @file GModelSpectralConst.i
23   - * @brief Spectral contant model class Python interface definition
24   - * @author J. Knodlseder
  23 + * @brief Spectral constant model class interface definition
  24 + * @author Juergen Knoedlseder
25 25 */
26 26 %{
27 27 /* Put headers and other declarations here that are needed for compilation */
... ... @@ -40,23 +40,30 @@ public:
40 40 // Constructors and destructors
41 41 GModelSpectralConst(void);
42 42 explicit GModelSpectralConst(const GXmlElement& xml);
  43 + explicit GModelSpectralConst(const double& norm);
43 44 GModelSpectralConst(const GModelSpectralConst& model);
44 45 virtual ~GModelSpectralConst(void);
45 46  
46   - // Implemented pure virtual methods
  47 + // Implemented pure virtual base class methods
47 48 virtual void clear(void);
48 49 virtual GModelSpectralConst* clone(void) const;
49 50 virtual std::string type(void) const;
50   - virtual double eval(const GEnergy& srcEng) const;
51   - virtual double eval_gradients(const GEnergy& srcEng) const;
52   - virtual double flux(const GEnergy& emin, const GEnergy& emax) const;
53   - virtual double eflux(const GEnergy& emin, const GEnergy& emax) const;
54   - virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
  51 + virtual double eval(const GEnergy& srcEng,
  52 + const GTime& srcTime) const;
  53 + virtual double eval_gradients(const GEnergy& srcEng,
  54 + const GTime& srcTime);
  55 + virtual double flux(const GEnergy& emin,
  56 + const GEnergy& emax) const;
  57 + virtual double eflux(const GEnergy& emin,
  58 + const GEnergy& emax) const;
  59 + virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax,
  60 + const GTime& time, GRan& ran) const;
55 61 virtual void read(const GXmlElement& xml);
56 62 virtual void write(GXmlElement& xml) const;
57 63  
58 64 // Other methods
59 65 double norm(void) const;
  66 + void norm(const double& norm);
60 67 };
61 68  
62 69  
... ...
pyext/GModelSpectralExpPlaw.i
  View file @938b5ea
... ... @@ -49,11 +49,18 @@ public:
49 49 virtual void clear(void);
50 50 virtual GModelSpectralExpPlaw* clone(void) const;
51 51 virtual std::string type(void) const;
52   - virtual double eval(const GEnergy& srcEng) const;
53   - virtual double eval_gradients(const GEnergy& srcEng) const;
54   - virtual double flux(const GEnergy& emin, const GEnergy& emax) const;
55   - virtual double eflux(const GEnergy& emin, const GEnergy& emax) const;
56   - virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
  52 + virtual double eval(const GEnergy& srcEng,
  53 + const GTime& srcTime) const;
  54 + virtual double eval_gradients(const GEnergy& srcEng,
  55 + const GTime& srcTime);
  56 + virtual double flux(const GEnergy& emin,
  57 + const GEnergy& emax) const;
  58 + virtual double eflux(const GEnergy& emin,
  59 + const GEnergy& emax) const;
  60 + virtual GEnergy mc(const GEnergy& emin,
  61 + const GEnergy& emax,
  62 + const GTime& time,
  63 + GRan& ran) const;
57 64 virtual void read(const GXmlElement& xml);
58 65 virtual void write(GXmlElement& xml) const;
59 66  
... ...
pyext/GModelSpectralFunc.i
  View file @938b5ea
... ... @@ -48,11 +48,16 @@ public:
48 48 virtual void clear(void);
49 49 virtual GModelSpectralFunc* clone(void) const;
50 50 virtual std::string type(void) const;
51   - virtual double eval(const GEnergy& srcEng) const;
52   - virtual double eval_gradients(const GEnergy& srcEng) const;
53   - virtual double flux(const GEnergy& emin, const GEnergy& emax) const;
54   - virtual double eflux(const GEnergy& emin, const GEnergy& emax) const;
55   - virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
  51 + virtual double eval(const GEnergy& srcEng,
  52 + const GTime& srcTime) const;
  53 + virtual double eval_gradients(const GEnergy& srcEng,
  54 + const GTime& srcTime);
  55 + virtual double flux(const GEnergy& emin,
  56 + const GEnergy& emax) const;
  57 + virtual double eflux(const GEnergy& emin,
  58 + const GEnergy& emax) const;
  59 + virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax,
  60 + const GTime& time, GRan& ran) const;
56 61 virtual void read(const GXmlElement& xml);
57 62 virtual void write(GXmlElement& xml) const;
58 63  
... ...
pyext/GModelSpectralLogParabola.i
  View file @938b5ea
... ... @@ -50,14 +50,18 @@ public:
50 50 virtual void clear(void);
51 51 virtual GModelSpectralLogParabola* clone(void) const;
52 52 virtual std::string type(void) const;
53   - virtual double eval(const GEnergy& srcEng) const;
54   - virtual double eval_gradients(const GEnergy& srcEng) const;
  53 + virtual double eval(const GEnergy& srcEng,
  54 + const GTime& srcTime) const;
  55 + virtual double eval_gradients(const GEnergy& srcEng,
  56 + const GTime& srcTime);
55 57 virtual double flux(const GEnergy& emin,
56 58 const GEnergy& emax) const;
57 59 virtual double eflux(const GEnergy& emin,
58 60 const GEnergy& emax) const;
59 61 virtual GEnergy mc(const GEnergy& emin,
60   - const GEnergy& emax, GRan& ran) const;
  62 + const GEnergy& emax,
  63 + const GTime& time,
  64 + GRan& ran) const;
61 65 virtual void read(const GXmlElement& xml);
62 66 virtual void write(GXmlElement& xml) const;
63 67  
... ...
pyext/GModelSpectralNodes.i
  View file @938b5ea
1 1 /***************************************************************************
2   - * GModelSpectralNodes.i - Spectral nodes model class *
  2 + * GModelSpectralNodes.i - Spectral nodes model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2013 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -21,7 +21,7 @@
21 21 /**
22 22 * @file GModelSpectralNodes.i
23 23 * @brief Spectral nodes model class Python interface definition
24   - * @author J. Knoedlseder
  24 + * @author Juergen Knoedlseder
25 25 */
26 26 %{
27 27 /* Put headers and other declarations here that are needed for compilation */
... ... @@ -47,11 +47,18 @@ public:
47 47 virtual void clear(void);
48 48 virtual GModelSpectralNodes* clone(void) const;
49 49 virtual std::string type(void) const;
50   - virtual double eval(const GEnergy& srcEng) const;
51   - virtual double eval_gradients(const GEnergy& srcEng) const;
52   - virtual double flux(const GEnergy& emin, const GEnergy& emax) const;
53   - virtual double eflux(const GEnergy& emin, const GEnergy& emax) const;
54   - virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
  50 + virtual double eval(const GEnergy& srcEng,
  51 + const GTime& srcTime) const;
  52 + virtual double eval_gradients(const GEnergy& srcEng,
  53 + const GTime& srcTime);
  54 + virtual double flux(const GEnergy& emin,
  55 + const GEnergy& emax) const;
  56 + virtual double eflux(const GEnergy& emin,
  57 + const GEnergy& emax) const;
  58 + virtual GEnergy mc(const GEnergy& emin,
  59 + const GEnergy& emax,
  60 + const GTime& time,
  61 + GRan& ran) const;
55 62 virtual void read(const GXmlElement& xml);
56 63 virtual void write(GXmlElement& xml) const;
57 64 };
... ...
pyext/GModelSpectralPlaw.i
  View file @938b5ea
... ... @@ -49,11 +49,18 @@ public:
49 49 virtual void clear(void);
50 50 virtual GModelSpectralPlaw* clone(void) const;
51 51 virtual std::string type(void) const;
52   - virtual double eval(const GEnergy& srcEng) const;
53   - virtual double eval_gradients(const GEnergy& srcEng) const;
54   - virtual double flux(const GEnergy& emin, const GEnergy& emax) const;
55   - virtual double eflux(const GEnergy& emin, const GEnergy& emax) const;
56   - virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
  52 + virtual double eval(const GEnergy& srcEng,
  53 + const GTime& srcTime) const;
  54 + virtual double eval_gradients(const GEnergy& srcEng,
  55 + const GTime& srcTime);
  56 + virtual double flux(const GEnergy& emin,
  57 + const GEnergy& emax) const;
  58 + virtual double eflux(const GEnergy& emin,
  59 + const GEnergy& emax) const;
  60 + virtual GEnergy mc(const GEnergy& emin,
  61 + const GEnergy& emax,
  62 + const GTime& time,
  63 + GRan& ran) const;
57 64 virtual void read(const GXmlElement& xml);
58 65 virtual void write(GXmlElement& xml) const;
59 66  
... ...
pyext/GModelSpectralPlaw2.i
  View file @938b5ea
1 1 /***************************************************************************
2   - * GModelSpectralPlaw2.i - Spectral power law model class *
  2 + * GModelSpectralPlaw2.i - Spectral power law model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011 by Jurgen Knodlseder *
  4 + * copyright (C) 2011-2013 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -21,7 +21,7 @@
21 21 /**
22 22 * @file GModelSpectralPlaw2.i
23 23 * @brief Flux normalized power law spectral model class Python interface definition
24   - * @author J. Knodlseder
  24 + * @author Juergen Knoedlseder
25 25 */
26 26 %{
27 27 /* Put headers and other declarations here that are needed for compilation */
... ... @@ -48,11 +48,18 @@ public:
48 48 virtual void clear(void);
49 49 virtual GModelSpectralPlaw2* clone(void) const;
50 50 virtual std::string type(void) const;
51   - virtual double eval(const GEnergy& srcEng) const;
52   - virtual double eval_gradients(const GEnergy& srcEng) const;
53   - virtual double flux(const GEnergy& emin, const GEnergy& emax) const;
54   - virtual double eflux(const GEnergy& emin, const GEnergy& emax) const;
55   - virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax, GRan& ran) const;
  51 + virtual double eval(const GEnergy& srcEng,
  52 + const GTime& srcTime) const;
  53 + virtual double eval_gradients(const GEnergy& srcEng,
  54 + const GTime& srcTime);
  55 + virtual double flux(const GEnergy& emin,
  56 + const GEnergy& emax) const;
  57 + virtual double eflux(const GEnergy& emin,
  58 + const GEnergy& emax) const;
  59 + virtual GEnergy mc(const GEnergy& emin,
  60 + const GEnergy& emax,
  61 + const GTime& time,
  62 + GRan& ran) const;
56 63 virtual void read(const GXmlElement& xml);
57 64 virtual void write(GXmlElement& xml) const;
58 65  
... ...
src/model/GModelSpectralConst.cpp
  View file @938b5ea
... ... @@ -42,7 +42,7 @@ const GModelSpectralRegistry g_spectral_const_registry(&amp;g_spectral_const_seed);
42 42 /* __ Method name definitions ____________________________________________ */
43 43 #define G_FLUX "GModelSpectralConst::flux(GEnergy&, GEnergy&)"
44 44 #define G_EFLUX "GModelSpectralConst::eflux(GEnergy&, GEnergy&)"
45   -#define G_MC "GModelSpectralConst::mc(GEnergy&, GEnergy&, GRan&)"
  45 +#define G_MC "GModelSpectralConst::mc(GEnergy&, GEnergy&, GTime&, GRan&)"
46 46 #define G_READ "GModelSpectralConst::read(GXmlElement&)"
47 47 #define G_WRITE "GModelSpectralConst::write(GXmlElement&)"
48 48  
... ... @@ -77,9 +77,9 @@ GModelSpectralConst::GModelSpectralConst(void) : GModelSpectral()
77 77 *
78 78 * @param[in] xml XML element.
79 79 *
80   - * Creates instance of a constant spectral model by extracting information
81   - * from an XML element. See GModelSpectralConst::read() for more information
82   - * about the expected structure of the XML element.
  80 + * Constructs constant spectral model by extracting information from an XML
  81 + * element. See the read() method for more information about the expected
  82 + * structure of the XML element.
83 83 ***************************************************************************/
84 84 GModelSpectralConst::GModelSpectralConst(const GXmlElement& xml) :
85 85 GModelSpectral()
... ... @@ -96,6 +96,27 @@ GModelSpectralConst::GModelSpectralConst(const GXmlElement&amp; xml) :
96 96  
97 97  
98 98 /***********************************************************************//**
  99 + * @brief Value constructor
  100 + *
  101 + * @param[in] norm Normalization factor (ph/cm2/s).
  102 + *
  103 + * Constructs constant spectral model by setting the normalization factor.
  104 + ***************************************************************************/
  105 +GModelSpectralConst::GModelSpectralConst(const double& norm) :
  106 + GModelSpectral()
  107 +{
  108 + // Initialise members
  109 + init_members();
  110 +
  111 + // Set normalization factor
  112 + m_norm.value(norm);
  113 +
  114 + // Return
  115 + return;
  116 +}
  117 +
  118 +
  119 +/***********************************************************************//**
99 120 * @brief Copy constructor
100 121 *
101 122 * @param[in] model Spectral constant model.
... ... @@ -136,8 +157,8 @@ GModelSpectralConst::~GModelSpectralConst(void)
136 157 /***********************************************************************//**
137 158 * @brief Assignment operator
138 159 *
139   - * @param[in] model Spectral constant model.
140   - * @return Spectral model.
  160 + * @param[in] model Constant spectral model.
  161 + * @return Constant spectral model.
141 162 ***************************************************************************/
142 163 GModelSpectralConst& GModelSpectralConst::operator=(const GModelSpectralConst& model)
143 164 {
... ... @@ -170,7 +191,7 @@ GModelSpectralConst&amp; GModelSpectralConst::operator=(const GModelSpectralConst&amp; m
170 191 ==========================================================================*/
171 192  
172 193 /***********************************************************************//**
173   - * @brief Clear instance
  194 + * @brief Clear constant spectral model
174 195 ***************************************************************************/
175 196 void GModelSpectralConst::clear(void)
176 197 {
... ... @@ -188,28 +209,36 @@ void GModelSpectralConst::clear(void)
188 209  
189 210  
190 211 /***********************************************************************//**
191   - * @brief Clone instance
  212 + * @brief Clone constant spectral model
192 213 *
193   - * @return Pointer to deep copy of spectral model.
  214 + * @return Pointer to deep copy of constant spectral model.
194 215 ***************************************************************************/
195 216 GModelSpectralConst* GModelSpectralConst::clone(void) const
196 217 {
  218 + // Clone constant spectral model
197 219 return new GModelSpectralConst(*this);
198 220 }
199 221  
200 222  
201 223 /***********************************************************************//**
202   - * @brief Evaluate function
  224 + * @brief Evaluate model
  225 + *
  226 + * @param[in] srcEng True photon energy.
  227 + * @param[in] srcTime True photon arrival time.
  228 + * @return Model value (ph/cm2/s/MeV).
203 229 *
204   - * @param[in] srcEng True energy of photon.
205   - * @return Model value.
  230 + * The spectral model is defined by
  231 + *
  232 + * \f[
  233 + * S_{\rm E}(E | t) = {\tt m\_norm}
  234 + * \f]
206 235 *
207   - * The spectral model is defined as
208   - * \f[I(E)=norm\f]
209 236 * where
210   - * \f$norm\f$ is the normalization of the function.
  237 + * \f${\tt m\_norm}\f$ is the normalization constant in units of
  238 + * ph/cm2/s/MeV.
211 239 ***************************************************************************/
212   -double GModelSpectralConst::eval(const GEnergy& srcEng) const
  240 +double GModelSpectralConst::eval(const GEnergy& srcEng,
  241 + const GTime& srcTime) const
213 242 {
214 243 // Compute function value
215 244 double value = norm();
... ... @@ -222,21 +251,32 @@ double GModelSpectralConst::eval(const GEnergy&amp; srcEng) const
222 251 /***********************************************************************//**
223 252 * @brief Evaluate function and gradients
224 253 *
225   - * @param[in] srcEng True energy of photon.
226   - * @return Model value.
  254 + * @param[in] srcEng True photon energy.
  255 + * @param[in] srcTime True photon arrival time.
  256 + * @return Model value (ph/cm2/s/MeV).
  257 + *
  258 + * The spectral model is defined by
  259 + *
  260 + * \f[
  261 + * S_{\rm E}(E | t) = {\tt m\_norm}
  262 + * \f]
227 263 *
228   - * The spectral model is defined as
229   - * \f[I(E)=norm\f]
230 264 * where
231   - * \f$norm=n_s n_v\f$ is the normalization of the function.
232   - * Note that the normalization is factorised into a scaling factor and a
233   - * value and that the method is expected to return the gradient with respect
234   - * to the parameter value \f$n_v\f$.
  265 + * \f${\tt m\_norm}\f$ is the normalization constant in units of
  266 + * ph/cm2/s/MeV.
  267 + *
  268 + * The normalization constant is factorized into a scale factor and a value
  269 + * factor. The partial derivative with respect to the normalization constant
  270 + * value factor is given by
235 271 *
236   - * The partial derivative of the normalization value is given by
237   - * \f[dI/dn_v=n_s\f]
  272 + * \f[
  273 + * \frac{\delta S_{\rm E}(E | t)}
  274 + * {\delta {\tt m\_norm.factor\_value()}} =
  275 + * {\tt m\_norm.scale()}
  276 + * \f]
238 277 ***************************************************************************/
239   -double GModelSpectralConst::eval_gradients(const GEnergy& srcEng) const
  278 +double GModelSpectralConst::eval_gradients(const GEnergy& srcEng,
  279 + const GTime& srcTime)
240 280 {
241 281 // Compute function value
242 282 double value = norm();
... ... @@ -244,8 +284,8 @@ double GModelSpectralConst::eval_gradients(const GEnergy&amp; srcEng) const
244 284 // Compute partial derivatives of the parameter values
245 285 double g_norm = (m_norm.isfree()) ? m_norm.scale() : 0.0;
246 286  
247   - // Set gradients (circumvent const correctness)
248   - const_cast<GModelSpectralConst*>(this)->m_norm.factor_gradient(g_norm);
  287 + // Set gradient
  288 + m_norm.factor_gradient(g_norm);
249 289  
250 290 // Return
251 291 return value;
... ... @@ -253,7 +293,7 @@ double GModelSpectralConst::eval_gradients(const GEnergy&amp; srcEng) const
253 293  
254 294  
255 295 /***********************************************************************//**
256   - * @brief Returns model photon flux between [emin, emax] (units: ph/cm2/s)
  296 + * @brief Returns model photon flux between [emin, emax] (ph/cm2/s)
257 297 *
258 298 * @param[in] emin Minimum photon energy.
259 299 * @param[in] emax Maximum photon energy.
... ... @@ -263,14 +303,21 @@ double GModelSpectralConst::eval_gradients(const GEnergy&amp; srcEng) const
263 303 * Energy range is invalid (emin < emax required).
264 304 *
265 305 * Computes
266   - * \f[\int_{E_{\rm min}}^{E_{\rm max}} I(E) dE\f]
  306 + *
  307 + * \f[
  308 + * \int_{E_{\rm min}}^{E_{\rm max}} S_{\rm E}(E | t) dE
  309 + * \f]
  310 + *
267 311 * where
268   - * \f$E_{\rm min}\f$ and \f$E_{\rm max}\f$ are the minimum and maximum
269   - * energy, respectively, and
270   - * \f$I(E)\f$ is the spectral model (units: ph/cm2/s/MeV).
  312 + * - \f$E_{\rm min}\f$ and \f$E_{\rm max}\f$ are the minimum and maximum
  313 + * energy, respectively, and
  314 + * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
271 315 * The integration is done analytically.
  316 + *
  317 + * @todo What about returning a zero flux if the energy range is invalid?
272 318 ***************************************************************************/
273   -double GModelSpectralConst::flux(const GEnergy& emin, const GEnergy& emax) const
  319 +double GModelSpectralConst::flux(const GEnergy& emin,
  320 + const GEnergy& emax) const
274 321 {
275 322 // Throw an exception if energy range is invalid
276 323 if (emin >= emax) {
... ... @@ -288,7 +335,7 @@ double GModelSpectralConst::flux(const GEnergy&amp; emin, const GEnergy&amp; emax) const
288 335  
289 336  
290 337 /***********************************************************************//**
291   - * @brief Returns model energy flux between [emin, emax] (units: erg/cm2/s)
  338 + * @brief Returns model energy flux between [emin, emax] (erg/cm2/s)
292 339 *
293 340 * @param[in] emin Minimum photon energy.
294 341 * @param[in] emax Maximum photon energy.
... ... @@ -298,24 +345,31 @@ double GModelSpectralConst::flux(const GEnergy&amp; emin, const GEnergy&amp; emax) const
298 345 * Energy range is invalid (emin < emax required).
299 346 *
300 347 * Computes
301   - * \f[\int_{E_{\rm min}}^{E_{\rm max}} I(E) E dE\f]
  348 + *
  349 + * \f[
  350 + * \int_{E_{\rm min}}^{E_{\rm max}} S_{\rm E}(E | t) E dE
  351 + * \f]
  352 + *
302 353 * where
303   - * \f$E_{\rm min}\f$ and \f$E_{\rm max}\f$ are the minimum and maximum
304   - * energy, respectively, and
305   - * \f$I(E)\f$ is the spectral model (units: ph/cm2/s/MeV).
  354 + * - \f$E_{\rm min}\f$ and \f$E_{\rm max}\f$ are the minimum and maximum
  355 + * energy, respectively, and
  356 + * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
306 357 * The integration is done analytically.
  358 + *
  359 + * @todo What about returning a zero flux if the energy range is invalid?
307 360 ***************************************************************************/
308   -double GModelSpectralConst::eflux(const GEnergy& emin, const GEnergy& emax) const
  361 +double GModelSpectralConst::eflux(const GEnergy& emin,
  362 + const GEnergy& emax) const
309 363 {
310 364 // Throw an exception if energy range is invalid
311 365 if (emin >= emax) {
312 366 throw GException::erange_invalid(G_EFLUX, emin.MeV(), emax.MeV(),
313 367 "Minimum energy < maximum energy required.");
314   -
315 368 }
316 369  
317 370 // Compute flux for a constant model
318   - double flux = norm() * 0.5 * (emax.MeV()*emax.MeV() - emin.MeV()*emin.MeV());
  371 + double flux = norm() * 0.5 * (emax.MeV()*emax.MeV() -
  372 + emin.MeV()*emin.MeV());
319 373  
320 374 // Convert from MeV/cm2/s to erg/cm2/s
321 375 flux *= MeV2erg;
... ... @@ -330,7 +384,8 @@ double GModelSpectralConst::eflux(const GEnergy&amp; emin, const GEnergy&amp; emax) cons
330 384 *
331 385 * @param[in] emin Minimum photon energy.
332 386 * @param[in] emax Maximum photon energy.
333   - * @param[in] ran Random number generator.
  387 + * @param[in] time True photon arrival time.
  388 + * @param[in,out] ran Random number generator.
334 389 * @return Energy.
335 390 *
336 391 * @exception GException::erange_invalid
... ... @@ -339,14 +394,15 @@ double GModelSpectralConst::eflux(const GEnergy&amp; emin, const GEnergy&amp; emax) cons
339 394 * Returns Monte Carlo energy by randomly drawing from a constant between
340 395 * the minimum and maximum photon energy.
341 396 ***************************************************************************/
342   -GEnergy GModelSpectralConst::mc(const GEnergy& emin, const GEnergy& emax,
343   - GRan& ran) const
  397 +GEnergy GModelSpectralConst::mc(const GEnergy& emin,
  398 + const GEnergy& emax,
  399 + const GTime& time,
  400 + GRan& ran) const
344 401 {
345 402 // Throw an exception if energy range is invalid
346 403 if (emin >= emax) {
347 404 throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(),
348 405 "Minimum energy < maximum energy required.");
349   -
350 406 }
351 407  
352 408 // Allocate energy
... ... @@ -373,9 +429,18 @@ GEnergy GModelSpectralConst::mc(const GEnergy&amp; emin, const GEnergy&amp; emax,
373 429 * @exception GException::model_invalid_parnames
374 430 * Invalid model parameter name found in XML element.
375 431 *
376   - * Read the function information from an XML element and load the nodes
377   - * from the associated file. The XML element is required to have a parameter
378   - * named "Normalization" or "Value".
  432 + * Reads the spectral information from an XML element having either the
  433 + * format
  434 + *
  435 + * <spectrum type="ConstantValue">
  436 + * <parameter name="Value" scale="1" min="0" max="1000" value="1" free="1"/>
  437 + * </spectrum>
  438 + *
  439 + * or the format
  440 + *
  441 + * <spectrum type="ConstantValue">
  442 + * <parameter name="Normalization" scale="1" min="0" max="1000" value="1" free="1"/>
  443 + * </spectrum>
379 444 ***************************************************************************/
380 445 void GModelSpectralConst::read(const GXmlElement& xml)
381 446 {
... ... @@ -416,9 +481,11 @@ void GModelSpectralConst::read(const GXmlElement&amp; xml)
416 481 * @exception GException::model_invalid_parnames
417 482 * Invalid model parameter names found in XML element.
418 483 *
419   - * Write the spectral constant information into an XML element. The XML
420   - * element has to be of type "ConstantValue" and will have 1 parameter leaf
421   - * named "Normalization" or "Value" (default).
  484 + * Writes the spectral information into an XML element with the format
  485 + *
  486 + * <spectrum type="ConstantValue">
  487 + * <parameter name="Value" scale="1" min="0" max="1000" value="1" free="1"/>
  488 + * </spectrum>
422 489 ***************************************************************************/
423 490 void GModelSpectralConst::write(GXmlElement& xml) const
424 491 {
... ...
src/model/GModelSpectralExpPlaw.cpp
  View file @938b5ea
... ... @@ -235,7 +235,8 @@ GModelSpectralExpPlaw* GModelSpectralExpPlaw::clone(void) const
235 235 /***********************************************************************//**
236 236 * @brief Evaluate function
237 237 *
238   - * @param[in] srcEng True energy of photon.
  238 + * @param[in] srcEng True photon energy.
  239 + * @param[in] srcTime True photon arrival time.
239 240 * @return Model value.
240 241 *
241 242 * The power law function is defined as
... ... @@ -250,7 +251,8 @@ GModelSpectralExpPlaw* GModelSpectralExpPlaw::clone(void) const
250 251 * may result. We should add a test that prevents using invalid
251 252 * values.
252 253 ***************************************************************************/
253   -double GModelSpectralExpPlaw::eval(const GEnergy& srcEng) const
  254 +double GModelSpectralExpPlaw::eval(const GEnergy& srcEng,
  255 + const GTime& srcTime) const
254 256 {
255 257 // Compute function value
256 258 double energy = srcEng.MeV();
... ... @@ -285,7 +287,8 @@ double GModelSpectralExpPlaw::eval(const GEnergy&amp; srcEng) const
285 287 /***********************************************************************//**
286 288 * @brief Evaluate function and gradients
287 289 *
288   - * @param[in] srcEng True energy of photon.
  290 + * @param[in] srcEng True photon energy.
  291 + * @param[in] srcTime True photon arrival time.
289 292 * @return Model value.
290 293 *
291 294 * The power law function is defined as
... ... @@ -309,7 +312,8 @@ double GModelSpectralExpPlaw::eval(const GEnergy&amp; srcEng) const
309 312 * may result. We should add a test that prevents using invalid
310 313 * values.
311 314 ***************************************************************************/
312   -double GModelSpectralExpPlaw::eval_gradients(const GEnergy& srcEng) const
  315 +double GModelSpectralExpPlaw::eval_gradients(const GEnergy& srcEng,
  316 + const GTime& srcTime)
313 317 {
314 318 // Compute function value
315 319 double energy = srcEng.MeV();
... ... @@ -449,7 +453,8 @@ double GModelSpectralExpPlaw::eflux(const GEnergy&amp; emin, const GEnergy&amp; emax) co
449 453 *
450 454 * @param[in] emin Minimum photon energy.
451 455 * @param[in] emax Maximum photon energy.
452   - * @param[in] ran Random number generator.
  456 + * @param[in] time True photon arrival time.
  457 + * @param[in,out] ran Random number generator.
453 458 * @return Energy.
454 459 *
455 460 * @exception GException::erange_invalid
... ... @@ -460,14 +465,15 @@ double GModelSpectralExpPlaw::eflux(const GEnergy&amp; emin, const GEnergy&amp; emax) co
460 465 * method. First, a random energy within [emin, emax] is drawn from an
461 466 * exponential distribution.
462 467 ***************************************************************************/
463   -GEnergy GModelSpectralExpPlaw::mc(const GEnergy& emin, const GEnergy& emax,
464   - GRan& ran) const
  468 +GEnergy GModelSpectralExpPlaw::mc(const GEnergy& emin,
  469 + const GEnergy& emax,
  470 + const GTime& time,
  471 + GRan& ran) const
465 472 {
466 473 // Throw an exception if energy range is invalid
467 474 if (emin >= emax) {
468 475 throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(),
469 476 "Minimum energy < maximum energy required.");
470   -
471 477 }
472 478  
473 479 // Allocate energy
... ...
src/model/GModelSpectralFunc.cpp
  View file @938b5ea
... ... @@ -107,8 +107,8 @@ GModelSpectralFunc::GModelSpectralFunc(const std::string&amp; filename)
107 107 * from an XML element. See GModelSpectralFunc::read() for more information
108 108 * about the expected structure of the XML element.
109 109 ***************************************************************************/
110   -GModelSpectralFunc::GModelSpectralFunc(const GXmlElement& xml)
111   - : GModelSpectral()
  110 +GModelSpectralFunc::GModelSpectralFunc(const GXmlElement& xml) :
  111 + GModelSpectral()
112 112 {
113 113 // Initialise members
114 114 init_members();
... ... @@ -126,8 +126,8 @@ GModelSpectralFunc::GModelSpectralFunc(const GXmlElement&amp; xml)
126 126 *
127 127 * @param[in] model Spectral function model.
128 128 ***************************************************************************/
129   -GModelSpectralFunc::GModelSpectralFunc(const GModelSpectralFunc& model)
130   - : GModelSpectral(model)
  129 +GModelSpectralFunc::GModelSpectralFunc(const GModelSpectralFunc& model) :
  130 + GModelSpectral(model)
131 131 {
132 132 // Initialise members
133 133 init_members();
... ... @@ -222,9 +222,11 @@ GModelSpectralFunc* GModelSpectralFunc::clone(void) const
222 222  
223 223  
224 224 /***********************************************************************//**
225   - * @brief Evaluate function
  225 + * @brief Evaluate model
226 226 *
227   - * @param[in] srcEng True energy of photon.
  227 + * @param[in] srcEng True photon energy.
  228 + * @param[in] srcTime True photon arrival time.
  229 + * @return Model value (ph/cm2/s/MeV).
228 230 *
229 231 * The spectral model is defined as
230 232 * \f[I(E)=norm f(E)\f]
... ... @@ -233,7 +235,8 @@ GModelSpectralFunc* GModelSpectralFunc::clone(void) const
233 235 * Note that the node energies are stored as log10 of energy in units of
234 236 * MeV.
235 237 ***************************************************************************/
236   -double GModelSpectralFunc::eval(const GEnergy& srcEng) const
  238 +double GModelSpectralFunc::eval(const GEnergy& srcEng,
  239 + const GTime& srcTime) const
237 240 {
238 241 // Interpolate function. This is done in log10-log10 space, but the
239 242 // linear value is returned.
... ... @@ -264,7 +267,9 @@ double GModelSpectralFunc::eval(const GEnergy&amp; srcEng) const
264 267 /***********************************************************************//**
265 268 * @brief Evaluate function and gradients
266 269 *
267   - * @param[in] srcEng True energy of photon.
  270 + * @param[in] srcEng True photon energy.
  271 + * @param[in] srcTime True photon arrival time.
  272 + * @return Model value (ph/cm2/s/MeV).
268 273 *
269 274 * The spectral model is defined as
270 275 * \f[I(E)=norm f(E)\f]
... ... @@ -277,7 +282,8 @@ double GModelSpectralFunc::eval(const GEnergy&amp; srcEng) const
277 282 * The partial derivative of the normalization value is given by
278 283 * \f[dI/dn_v=n_s f(E)\f]
279 284 ***************************************************************************/
280   -double GModelSpectralFunc::eval_gradients(const GEnergy& srcEng) const
  285 +double GModelSpectralFunc::eval_gradients(const GEnergy& srcEng,
  286 + const GTime& srcTime)
281 287 {
282 288 // Interpolate function. This is done in log10-log10 space, but the
283 289 // linear value is returned.
... ... @@ -507,7 +513,8 @@ double GModelSpectralFunc::eflux(const GEnergy&amp; emin, const GEnergy&amp; emax) const
507 513 *
508 514 * @param[in] emin Minimum photon energy.
509 515 * @param[in] emax Maximum photon energy.
510   - * @param[in] ran Random number generator.
  516 + * @param[in] time True photon arrival time.
  517 + * @param[in,out] ran Random number generator.
511 518 * @return Energy.
512 519 *
513 520 * @exception GException::erange_invalid
... ... @@ -516,14 +523,15 @@ double GModelSpectralFunc::eflux(const GEnergy&amp; emin, const GEnergy&amp; emax) const
516 523 * Simulates a random energy in the interval [emin, emax] for a spectral
517 524 * function.
518 525 ***************************************************************************/
519   -GEnergy GModelSpectralFunc::mc(const GEnergy& emin, const GEnergy& emax,
520   - GRan& ran) const
  526 +GEnergy GModelSpectralFunc::mc(const GEnergy& emin,
  527 + const GEnergy& emax,
  528 + const GTime& time,
  529 + GRan& ran) const
521 530 {
522 531 // Throw an exception if energy range is invalid
523 532 if (emin >= emax) {
524 533 throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(),
525 534 "Minimum energy < maximum energy required.");
526   -
527 535 }
528 536  
529 537 // Allocate energy
... ...
src/model/GModelSpectralLogParabola.cpp
  View file @938b5ea
... ... @@ -232,10 +232,11 @@ GModelSpectralLogParabola* GModelSpectralLogParabola::clone(void) const
232 232  
233 233  
234 234 /***********************************************************************//**
235   - * @brief Evaluate function
  235 + * @brief Evaluate model
236 236 *
237   - * @param[in] srcEng True energy of photon.
238   - * @return Function value.
  237 + * @param[in] srcEng True photon energy.
  238 + * @param[in] srcTime True photon arrival time.
  239 + * @return Model value (ph/cm2/s/MeV).
239 240 *
240 241 * The logparabola function is defined as
241 242 * \f[I(E)=norm (E/pivot)^{index-curvature*log(E/pivot)}\f]
... ... @@ -250,7 +251,8 @@ GModelSpectralLogParabola* GModelSpectralLogParabola::clone(void) const
250 251 * Furthermore, the method expects that energy!=0. Otherwise Inf or NaN
251 252 * may result.
252 253 ***************************************************************************/
253   -double GModelSpectralLogParabola::eval(const GEnergy& srcEng) const
  254 +double GModelSpectralLogParabola::eval(const GEnergy& srcEng,
  255 + const GTime& srcTime) const
254 256 {
255 257 // Compute function value
256 258 double energy = srcEng.MeV() / pivot();
... ... @@ -282,8 +284,9 @@ double GModelSpectralLogParabola::eval(const GEnergy&amp; srcEng) const
282 284 /***********************************************************************//**
283 285 * @brief Evaluate function and gradients
284 286 *
285   - * @param[in] srcEng True energy of photon.
286   - * @return Function value.
  287 + * @param[in] srcEng True photon energy.
  288 + * @param[in] srcTime True photon arrival time.
  289 + * @return Model value (ph/cm2/s/MeV).
287 290 *
288 291 * The logparabol function is defined as
289 292 * \f[I(E)=norm (E/pivot)^{index-curvature\ln(E/pivot)}\f]
... ... @@ -304,7 +307,8 @@ double GModelSpectralLogParabola::eval(const GEnergy&amp; srcEng) const
304 307 * Furthermore, the method expects that energy!=0. Otherwise Inf or NaN
305 308 * may result.
306 309 ***************************************************************************/
307   -double GModelSpectralLogParabola::eval_gradients(const GEnergy& srcEng) const
  310 +double GModelSpectralLogParabola::eval_gradients(const GEnergy& srcEng,
  311 + const GTime& srcTime)
308 312 {
309 313 // Compute function value
310 314 double energy = srcEng.MeV() / pivot();
... ... @@ -433,7 +437,8 @@ double GModelSpectralLogParabola::eflux(const GEnergy&amp; emin,
433 437 *
434 438 * @param[in] emin Minimum photon energy.
435 439 * @param[in] emax Maximum photon energy.
436   - * @param[in] ran Random number generator.
  440 + * @param[in] time True photon arrival time.
  441 + * @param[in,out] ran Random number generator.
437 442 * @return Energy.
438 443 *
439 444 * @exception GException::feature_not_implemented
... ... @@ -444,7 +449,8 @@ double GModelSpectralLogParabola::eflux(const GEnergy&amp; emin,
444 449 ***************************************************************************/
445 450 GEnergy GModelSpectralLogParabola::mc(const GEnergy& emin,
446 451 const GEnergy& emax,
447   - GRan& ran) const
  452 + const GTime& time,
  453 + GRan& ran) const
448 454 {
449 455 if (emin >= emax) {
450 456 throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(),
... ... @@ -455,7 +461,7 @@ GEnergy GModelSpectralLogParabola::mc(const GEnergy&amp; emin,
455 461 GEnergy energy;
456 462  
457 463 // Update cache
458   - update_mc_cache(emin, emax);
  464 + update_mc_cache(emin, emax, time);
459 465  
460 466 // Initialise energy
461 467 double eng;
... ... @@ -845,11 +851,13 @@ void GModelSpectralLogParabola::free_members(void)
845 851 *
846 852 * @param[in] emin Minimum photon energy.
847 853 * @param[in] emax Maximum photon energy.
  854 + * @param[in] time True photon arrival time.
848 855 *
849 856 * Updates the precomputation cache for Monte Carlo simulations.
850 857 ***************************************************************************/
851 858 void GModelSpectralLogParabola::update_mc_cache(const GEnergy& emin,
852   - const GEnergy& emax) const
  859 + const GEnergy& emax,
  860 + const GTime& time) const
853 861  
854 862 {
855 863 // Only update if boundaries have changed
... ... @@ -876,10 +884,11 @@ void GModelSpectralLogParabola::update_mc_cache(const GEnergy&amp; emin,
876 884  
877 885 // Plaw index defined by the slope of a straight line in the
878 886 // log-log-plane
879   - index_pl = std::log(eval(emin)/eval(emax))/std::log(emin.MeV()/emax.MeV());
  887 + index_pl = std::log(eval(emin,time)/eval(emax,time)) /
  888 + std::log(emin.MeV()/emax.MeV());
880 889  
881 890 // Plaw norm defined such that Plaw = LogParabola at emin
882   - m_mc_norm = eval(emin)/std::pow(emin.MeV()/pivot(),index_pl);
  891 + m_mc_norm = eval(emin,time) / std::pow(emin.MeV()/pivot(),index_pl);
883 892  
884 893 }
885 894  
... ...
src/model/GModelSpectralNodes.cpp
  View file @938b5ea
... ... @@ -199,7 +199,9 @@ GModelSpectralNodes* GModelSpectralNodes::clone(void) const
199 199 /***********************************************************************//**
200 200 * @brief Evaluate function
201 201 *
202   - * @param[in] srcEng True energy of photon.
  202 + * @param[in] srcEng True photon energy.
  203 + * @param[in] srcTime True photon arrival time.
  204 + * @return Model value (ph/cm2/s/MeV).
203 205 *
204 206 * The spectral model is defined as a piecewise power law between spectral
205 207 * nodes \f$(E_i, I_i)\f$,
... ... @@ -220,7 +222,8 @@ GModelSpectralNodes* GModelSpectralNodes::clone(void) const
220 222 * \f$(\log E_2, \log I_2)\f$
221 223 * to the requested energy \f$\log E\f$.
222 224 ***************************************************************************/
223   -double GModelSpectralNodes::eval(const GEnergy& srcEng) const
  225 +double GModelSpectralNodes::eval(const GEnergy& srcEng,
  226 + const GTime& srcTime) const
224 227 {
225 228 // Update evaluation cache
226 229 update_eval_cache();
... ... @@ -250,7 +253,9 @@ double GModelSpectralNodes::eval(const GEnergy&amp; srcEng) const
250 253 /***********************************************************************//**
251 254 * @brief Evaluate function and gradients
252 255 *
253   - * @param[in] srcEng True energy of photon.
  256 + * @param[in] srcEng True photon energy.
  257 + * @param[in] srcTime True photon arrival time.
  258 + * @return Model value (ph/cm2/s/MeV).
254 259 *
255 260 * The spectral model is defined as a piecewise power law between spectral
256 261 * nodes \f$(E_i, I_i)\f$,
... ... @@ -275,7 +280,8 @@ double GModelSpectralNodes::eval(const GEnergy&amp; srcEng) const
275 280 * \f[dI/dv_i=I(E) w_i / v_i\f]
276 281 * for \f$i=1,2\f$ and \f$0\f$ for all other nodes.
277 282 ***************************************************************************/
278   -double GModelSpectralNodes::eval_gradients(const GEnergy& srcEng) const
  283 +double GModelSpectralNodes::eval_gradients(const GEnergy& srcEng,
  284 + const GTime& srcTime)
279 285 {
280 286 // Update evaluation cache
281 287 update_eval_cache();
... ... @@ -525,20 +531,22 @@ double GModelSpectralNodes::eflux(const GEnergy&amp; emin, const GEnergy&amp; emax) cons
525 531 *
526 532 * @param[in] emin Minimum photon energy.
527 533 * @param[in] emax Maximum photon energy.
528   - * @param[in] ran Random number generator.
  534 + * @param[in] time True photon arrival time.
  535 + * @param[in,out] ran Random number generator.
529 536 * @return Energy.
530 537 *
531 538 * @exception GException::erange_invalid
532 539 * Energy range is invalid (emin < emax required).
533 540 ***************************************************************************/
534   -GEnergy GModelSpectralNodes::mc(const GEnergy& emin, const GEnergy& emax,
535   - GRan& ran) const
  541 +GEnergy GModelSpectralNodes::mc(const GEnergy& emin,
  542 + const GEnergy& emax,
  543 + const GTime& time,
  544 + GRan& ran) const
536 545 {
537 546 // Throw an exception if energy range is invalid
538 547 if (emin >= emax) {
539 548 throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(),
540 549 "Minimum energy < maximum energy required.");
541   -
542 550 }
543 551  
544 552 // Allocate energy
... ...
src/model/GModelSpectralPlaw.cpp
  View file @938b5ea
... ... @@ -230,7 +230,9 @@ GModelSpectralPlaw* GModelSpectralPlaw::clone(void) const
230 230 /***********************************************************************//**
231 231 * @brief Evaluate function
232 232 *
233   - * @param[in] srcEng True energy of photon.
  233 + * @param[in] srcEng True photon energy.
  234 + * @param[in] srcTime True photon arrival time.
  235 + * @return Model value (ph/cm2/s/MeV).
234 236 *
235 237 * The power law function is defined as
236 238 * \f[I(E)=norm (E/pivot)^{index}\f]
... ... @@ -244,7 +246,8 @@ GModelSpectralPlaw* GModelSpectralPlaw::clone(void) const
244 246 * Furthermore, the method expects that energy!=0. Otherwise Inf or NaN
245 247 * may result.
246 248 ***************************************************************************/
247   -double GModelSpectralPlaw::eval(const GEnergy& srcEng) const
  249 +double GModelSpectralPlaw::eval(const GEnergy& srcEng,
  250 + const GTime& srcTime) const
248 251 {
249 252 // Compute function value
250 253 double energy = srcEng.MeV() / pivot();
... ... @@ -274,7 +277,9 @@ double GModelSpectralPlaw::eval(const GEnergy&amp; srcEng) const
274 277 /***********************************************************************//**
275 278 * @brief Evaluate function and gradients
276 279 *
277   - * @param[in] srcEng True energy of photon.
  280 + * @param[in] srcEng True photon energy.
  281 + * @param[in] srcTime True photon arrival time.
  282 + * @return Model value (ph/cm2/s/MeV).
278 283 *
279 284 * The power law function is defined as
280 285 * \f[I(E)=norm (E/pivot)^{index}\f]
... ... @@ -296,7 +301,8 @@ double GModelSpectralPlaw::eval(const GEnergy&amp; srcEng) const
296 301 * Furthermore, the method expects that energy!=0. Otherwise Inf or NaN
297 302 * may result.
298 303 ***************************************************************************/
299   -double GModelSpectralPlaw::eval_gradients(const GEnergy& srcEng) const
  304 +double GModelSpectralPlaw::eval_gradients(const GEnergy& srcEng,
  305 + const GTime& srcTime)
300 306 {
301 307 // Compute function value
302 308 double energy = srcEng.MeV() / pivot();
... ... @@ -416,7 +422,8 @@ double GModelSpectralPlaw::eflux(const GEnergy&amp; emin, const GEnergy&amp; emax) const
416 422 *
417 423 * @param[in] emin Minimum photon energy.
418 424 * @param[in] emax Maximum photon energy.
419   - * @param[in] ran Random number generator.
  425 + * @param[in] time True photon arrival time.
  426 + * @param[in,out] ran Random number generator.
420 427 * @return Energy.
421 428 *
422 429 * @exception GException::erange_invalid
... ... @@ -424,14 +431,15 @@ double GModelSpectralPlaw::eflux(const GEnergy&amp; emin, const GEnergy&amp; emax) const
424 431 *
425 432 * Returns Monte Carlo energy by randomly drawing from a power law.
426 433 ***************************************************************************/
427   -GEnergy GModelSpectralPlaw::mc(const GEnergy& emin, const GEnergy& emax,
428   - GRan& ran) const
  434 +GEnergy GModelSpectralPlaw::mc(const GEnergy& emin,
  435 + const GEnergy& emax,
  436 + const GTime& time,
  437 + GRan& ran) const
429 438 {
430 439 // Throw an exception if energy range is invalid
431 440 if (emin >= emax) {
432 441 throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(),
433 442 "Minimum energy < maximum energy required.");
434   -
435 443 }
436 444  
437 445 // Allocate energy
... ...
src/model/GModelSpectralPlaw2.cpp
  View file @938b5ea
... ... @@ -223,7 +223,9 @@ GModelSpectralPlaw2* GModelSpectralPlaw2::clone(void) const
223 223 /***********************************************************************//**
224 224 * @brief Evaluate function
225 225 *
226   - * @param[in] srcEng True energy of photon.
  226 + * @param[in] srcEng True photon energy.
  227 + * @param[in] srcTime True photon arrival time.
  228 + * @return Model value (ph/cm2/s/MeV).
227 229 *
228 230 * The power law function is defined as
229 231 * \f[I(E)=integral (index+1)/(emax^{index+1}-emin^{index+1}) E^{index}\f]
... ... @@ -235,7 +237,8 @@ GModelSpectralPlaw2* GModelSpectralPlaw2::clone(void) const
235 237 * \f$emin\f$ is the lower energy limit, and
236 238 * \f$emax\f$ is the upper energy limit.
237 239 ***************************************************************************/
238   -double GModelSpectralPlaw2::eval(const GEnergy& srcEng) const
  240 +double GModelSpectralPlaw2::eval(const GEnergy& srcEng,
  241 + const GTime& srcTime) const
239 242 {
240 243 // Update precomputed values
241 244 update(srcEng);
... ... @@ -265,7 +268,9 @@ double GModelSpectralPlaw2::eval(const GEnergy&amp; srcEng) const
265 268 /***********************************************************************//**
266 269 * @brief Evaluate function and gradients
267 270 *
268   - * @param[in] srcEng True energy of photon.
  271 + * @param[in] srcEng True photon energy.
  272 + * @param[in] srcTime True photon arrival time.
  273 + * @return Model value (ph/cm2/s/MeV).
269 274 *
270 275 * The power law function is defined as
271 276 * \f[I(E)=integral (index+1)/(emax^{index+1}-emin^{index+1}) E^{index}\f]
... ... @@ -294,7 +299,8 @@ double GModelSpectralPlaw2::eval(const GEnergy&amp; srcEng) const
294 299 *
295 300 * No partial derivatives are supported for the energy boundaries.
296 301 ***************************************************************************/
297   -double GModelSpectralPlaw2::eval_gradients(const GEnergy& srcEng) const
  302 +double GModelSpectralPlaw2::eval_gradients(const GEnergy& srcEng,
  303 + const GTime& srcTime)
298 304 {
299 305 // Initialise gradients
300 306 double g_integral = 0.0;
... ... @@ -471,7 +477,8 @@ double GModelSpectralPlaw2::eflux(const GEnergy&amp; emin, const GEnergy&amp; emax) cons
471 477 *
472 478 * @param[in] emin Minimum photon energy.
473 479 * @param[in] emax Maximum photon energy.
474   - * @param[in] ran Random number generator.
  480 + * @param[in] time True photon arrival time.
  481 + * @param[in,out] ran Random number generator.
475 482 * @return Energy.
476 483 *
477 484 * @exception GException::erange_invalid
... ... @@ -479,14 +486,15 @@ double GModelSpectralPlaw2::eflux(const GEnergy&amp; emin, const GEnergy&amp; emax) cons
479 486 *
480 487 * Returns Monte Carlo energy by randomly drawing from a power law.
481 488 ***************************************************************************/
482   -GEnergy GModelSpectralPlaw2::mc(const GEnergy& emin, const GEnergy& emax,
483   - GRan& ran) const
  489 +GEnergy GModelSpectralPlaw2::mc(const GEnergy& emin,
  490 + const GEnergy& emax,
  491 + const GTime& time,
  492 + GRan& ran) const
484 493 {
485 494 // Throw an exception if energy range is invalid
486 495 if (emin >= emax) {
487 496 throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(),
488 497 "Minimum energy < maximum energy required.");
489   -
490 498 }
491 499  
492 500 // Allocate energy
... ...