Commit 94fedb394affe23b77c4c2f8bd40bb8581b0814e
1 parent
c63649fd
Add norm argument to value constructor, implement member access methods, update documentation.
Showing
3 changed files
with
237 additions
and
164 deletions
include/GModelSpectralFunc.hpp
| @@ -42,18 +42,23 @@ | @@ -42,18 +42,23 @@ | ||
| 42 | * | 42 | * |
| 43 | * @brief Spectral function model class | 43 | * @brief Spectral function model class |
| 44 | * | 44 | * |
| 45 | - * This class implements an arbitrary function as the spectral component of | ||
| 46 | - * the model. The function is defined as | ||
| 47 | - * \f[I(E)=norm f(E)\f] | 45 | + * This class implements an arbitrary function as the spectral model |
| 46 | + * component. The model is defined by | ||
| 47 | + * | ||
| 48 | + * \f[ | ||
| 49 | + * S_{\rm E}(E | t) = {\tt m\_norm} | ||
| 50 | + * \f] | ||
| 51 | + * | ||
| 48 | * where | 52 | * where |
| 49 | - * \f$norm\f$ is the normalization of the function. | 53 | + * - \f${\tt m\_norm}\f$ is the normalization of the function. |
| 50 | ***************************************************************************/ | 54 | ***************************************************************************/ |
| 51 | class GModelSpectralFunc : public GModelSpectral { | 55 | class GModelSpectralFunc : public GModelSpectral { |
| 52 | 56 | ||
| 53 | public: | 57 | public: |
| 54 | // Constructors and destructors | 58 | // Constructors and destructors |
| 55 | GModelSpectralFunc(void); | 59 | GModelSpectralFunc(void); |
| 56 | - explicit GModelSpectralFunc(const std::string& filename); | 60 | + explicit GModelSpectralFunc(const std::string& filename, |
| 61 | + const double& norm); | ||
| 57 | explicit GModelSpectralFunc(const GXmlElement& xml); | 62 | explicit GModelSpectralFunc(const GXmlElement& xml); |
| 58 | GModelSpectralFunc(const GModelSpectralFunc& model); | 63 | GModelSpectralFunc(const GModelSpectralFunc& model); |
| 59 | virtual ~GModelSpectralFunc(void); | 64 | virtual ~GModelSpectralFunc(void); |
| @@ -73,14 +78,19 @@ public: | @@ -73,14 +78,19 @@ public: | ||
| 73 | const GEnergy& emax) const; | 78 | const GEnergy& emax) const; |
| 74 | virtual double eflux(const GEnergy& emin, | 79 | virtual double eflux(const GEnergy& emin, |
| 75 | const GEnergy& emax) const; | 80 | const GEnergy& emax) const; |
| 76 | - virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax, | ||
| 77 | - const GTime& time, GRan& ran) const; | 81 | + virtual GEnergy mc(const GEnergy& emin, |
| 82 | + const GEnergy& emax, | ||
| 83 | + const GTime& time, | ||
| 84 | + GRan& ran) const; | ||
| 78 | virtual void read(const GXmlElement& xml); | 85 | virtual void read(const GXmlElement& xml); |
| 79 | virtual void write(GXmlElement& xml) const; | 86 | virtual void write(GXmlElement& xml) const; |
| 80 | virtual std::string print(void) const; | 87 | virtual std::string print(void) const; |
| 81 | 88 | ||
| 82 | // Other methods | 89 | // Other methods |
| 83 | - double norm(void) const { return m_norm.value(); } | 90 | + const std::string& filename(void) const; |
| 91 | + void filename(const std::string& filename); | ||
| 92 | + double norm(void) const; | ||
| 93 | + void norm(const double& norm); | ||
| 84 | 94 | ||
| 85 | protected: | 95 | protected: |
| 86 | // Protected methods | 96 | // Protected methods |
| @@ -129,4 +139,62 @@ std::string GModelSpectralFunc::type(void) const | @@ -129,4 +139,62 @@ std::string GModelSpectralFunc::type(void) const | ||
| 129 | return "FileFunction"; | 139 | return "FileFunction"; |
| 130 | } | 140 | } |
| 131 | 141 | ||
| 142 | + | ||
| 143 | +/***********************************************************************//** | ||
| 144 | + * @brief Return normalization factor | ||
| 145 | + * | ||
| 146 | + * @return Normalization factor. | ||
| 147 | + * | ||
| 148 | + * Returns the normalization factor. | ||
| 149 | + ***************************************************************************/ | ||
| 150 | +inline | ||
| 151 | +double GModelSpectralFunc::norm(void) const | ||
| 152 | +{ | ||
| 153 | + return (m_norm.value()); | ||
| 154 | +} | ||
| 155 | + | ||
| 156 | + | ||
| 157 | +/***********************************************************************//** | ||
| 158 | + * @brief Set normalization factor | ||
| 159 | + * | ||
| 160 | + * @param[in] norm Normalization factor. | ||
| 161 | + * | ||
| 162 | + * Sets the normalization factor. | ||
| 163 | + ***************************************************************************/ | ||
| 164 | +inline | ||
| 165 | +void GModelSpectralFunc::norm(const double& norm) | ||
| 166 | +{ | ||
| 167 | + m_norm.value(norm); | ||
| 168 | + return; | ||
| 169 | +} | ||
| 170 | + | ||
| 171 | + | ||
| 172 | +/***********************************************************************//** | ||
| 173 | + * @brief Return file name | ||
| 174 | + * | ||
| 175 | + * @return Name of node file. | ||
| 176 | + * | ||
| 177 | + * Returns the name of the file function node file. | ||
| 178 | + ***************************************************************************/ | ||
| 179 | +inline | ||
| 180 | +const std::string& GModelSpectralFunc::filename(void) const | ||
| 181 | +{ | ||
| 182 | + return (m_filename); | ||
| 183 | +} | ||
| 184 | + | ||
| 185 | + | ||
| 186 | +/***********************************************************************//** | ||
| 187 | + * @brief Loads nodes from node file and set filename | ||
| 188 | + * | ||
| 189 | + * @param[in] filename Node file name. | ||
| 190 | + * | ||
| 191 | + * Loads the nodes from a file function node file and sets the filename. | ||
| 192 | + ***************************************************************************/ | ||
| 193 | +inline | ||
| 194 | +void GModelSpectralFunc::filename(const std::string& filename) | ||
| 195 | +{ | ||
| 196 | + load_nodes(filename); | ||
| 197 | + return; | ||
| 198 | +} | ||
| 199 | + | ||
| 132 | #endif /* GMODELSPECTRALFUNC_HPP */ | 200 | #endif /* GMODELSPECTRALFUNC_HPP */ |
pyext/GModelSpectralFunc.i
| 1 | /*************************************************************************** | 1 | /*************************************************************************** |
| 2 | - * GModelSpectralFunc.i - Spectral function model class * | 2 | + * GModelSpectralFunc.i - Spectral function model class * |
| 3 | * ----------------------------------------------------------------------- * | 3 | * ----------------------------------------------------------------------- * |
| 4 | - * copyright (C) 2011 by Jurgen Knodlseder * | 4 | + * copyright (C) 2011-2013 by Juergen Knoedlseder * |
| 5 | * ----------------------------------------------------------------------- * | 5 | * ----------------------------------------------------------------------- * |
| 6 | * * | 6 | * * |
| 7 | * This program is free software: you can redistribute it and/or modify * | 7 | * This program is free software: you can redistribute it and/or modify * |
| @@ -21,7 +21,7 @@ | @@ -21,7 +21,7 @@ | ||
| 21 | /** | 21 | /** |
| 22 | * @file GModelSpectralFunc.i | 22 | * @file GModelSpectralFunc.i |
| 23 | * @brief Spectral function model class Python interface definition | 23 | * @brief Spectral function model class Python interface definition |
| 24 | - * @author J. Knodlseder | 24 | + * @author Juergen Knoedlseder |
| 25 | */ | 25 | */ |
| 26 | %{ | 26 | %{ |
| 27 | /* Put headers and other declarations here that are needed for compilation */ | 27 | /* Put headers and other declarations here that are needed for compilation */ |
| @@ -39,7 +39,8 @@ class GModelSpectralFunc : public GModelSpectral { | @@ -39,7 +39,8 @@ class GModelSpectralFunc : public GModelSpectral { | ||
| 39 | public: | 39 | public: |
| 40 | // Constructors and destructors | 40 | // Constructors and destructors |
| 41 | GModelSpectralFunc(void); | 41 | GModelSpectralFunc(void); |
| 42 | - explicit GModelSpectralFunc(const std::string& filename); | 42 | + explicit GModelSpectralFunc(const std::string& filename, |
| 43 | + const double& norm); | ||
| 43 | explicit GModelSpectralFunc(const GXmlElement& xml); | 44 | explicit GModelSpectralFunc(const GXmlElement& xml); |
| 44 | GModelSpectralFunc(const GModelSpectralFunc& model); | 45 | GModelSpectralFunc(const GModelSpectralFunc& model); |
| 45 | virtual ~GModelSpectralFunc(void); | 46 | virtual ~GModelSpectralFunc(void); |
| @@ -56,13 +57,18 @@ public: | @@ -56,13 +57,18 @@ public: | ||
| 56 | const GEnergy& emax) const; | 57 | const GEnergy& emax) const; |
| 57 | virtual double eflux(const GEnergy& emin, | 58 | virtual double eflux(const GEnergy& emin, |
| 58 | const GEnergy& emax) const; | 59 | const GEnergy& emax) const; |
| 59 | - virtual GEnergy mc(const GEnergy& emin, const GEnergy& emax, | ||
| 60 | - const GTime& time, GRan& ran) const; | 60 | + virtual GEnergy mc(const GEnergy& emin, |
| 61 | + const GEnergy& emax, | ||
| 62 | + const GTime& time, | ||
| 63 | + GRan& ran) const; | ||
| 61 | virtual void read(const GXmlElement& xml); | 64 | virtual void read(const GXmlElement& xml); |
| 62 | virtual void write(GXmlElement& xml) const; | 65 | virtual void write(GXmlElement& xml) const; |
| 63 | 66 | ||
| 64 | // Other methods | 67 | // Other methods |
| 65 | - double norm(void) const; | 68 | + const std::string& filename(void) const; |
| 69 | + void filename(const std::string& filename); | ||
| 70 | + double norm(void) const; | ||
| 71 | + void norm(const double& norm); | ||
| 66 | }; | 72 | }; |
| 67 | 73 | ||
| 68 | 74 |
src/model/GModelSpectralFunc.cpp
| @@ -44,7 +44,7 @@ const GModelSpectralRegistry g_spectral_func_registry(&g_spectral_func_seed); | @@ -44,7 +44,7 @@ const GModelSpectralRegistry g_spectral_func_registry(&g_spectral_func_seed); | ||
| 44 | /* __ Method name definitions ____________________________________________ */ | 44 | /* __ Method name definitions ____________________________________________ */ |
| 45 | #define G_FLUX "GModelSpectralFunc::flux(GEnergy&, GEnergy&)" | 45 | #define G_FLUX "GModelSpectralFunc::flux(GEnergy&, GEnergy&)" |
| 46 | #define G_EFLUX "GModelSpectralFunc::eflux(GEnergy&, GEnergy&)" | 46 | #define G_EFLUX "GModelSpectralFunc::eflux(GEnergy&, GEnergy&)" |
| 47 | -#define G_MC "GModelSpectralFunc::mc(GEnergy&, GEnergy&, GRan&)" | 47 | +#define G_MC "GModelSpectralFunc::mc(GEnergy&, GEnergy&, GTime&, GRan&)" |
| 48 | #define G_READ "GModelSpectralFunc::read(GXmlElement&)" | 48 | #define G_READ "GModelSpectralFunc::read(GXmlElement&)" |
| 49 | #define G_WRITE "GModelSpectralFunc::write(GXmlElement&)" | 49 | #define G_WRITE "GModelSpectralFunc::write(GXmlElement&)" |
| 50 | #define G_LOAD_NODES "GModelSpectralFunc::load_nodes(std::string&)" | 50 | #define G_LOAD_NODES "GModelSpectralFunc::load_nodes(std::string&)" |
| @@ -79,13 +79,15 @@ GModelSpectralFunc::GModelSpectralFunc(void) : GModelSpectral() | @@ -79,13 +79,15 @@ GModelSpectralFunc::GModelSpectralFunc(void) : GModelSpectral() | ||
| 79 | * @brief File constructor | 79 | * @brief File constructor |
| 80 | * | 80 | * |
| 81 | * @param[in] filename File name of nodes. | 81 | * @param[in] filename File name of nodes. |
| 82 | + * @param[in] norm Normalization factor. | ||
| 82 | * | 83 | * |
| 83 | - * Creates instance of a spectral function model from a list of nodes that is | ||
| 84 | - * found in the specified file. See GModelSpectralFunc::load_nodes() for more | 84 | + * Constructs spectral file function model from a list of nodes that is found |
| 85 | + * in the specified ASCII file. See the load_nodes() method for more | ||
| 85 | * information about the expected structure of the file. | 86 | * information about the expected structure of the file. |
| 86 | ***************************************************************************/ | 87 | ***************************************************************************/ |
| 87 | -GModelSpectralFunc::GModelSpectralFunc(const std::string& filename) | ||
| 88 | - : GModelSpectral() | 88 | +GModelSpectralFunc::GModelSpectralFunc(const std::string& filename, |
| 89 | + const double& norm) : | ||
| 90 | + GModelSpectral() | ||
| 89 | { | 91 | { |
| 90 | // Initialise members | 92 | // Initialise members |
| 91 | init_members(); | 93 | init_members(); |
| @@ -103,9 +105,9 @@ GModelSpectralFunc::GModelSpectralFunc(const std::string& filename) | @@ -103,9 +105,9 @@ GModelSpectralFunc::GModelSpectralFunc(const std::string& filename) | ||
| 103 | * | 105 | * |
| 104 | * @param[in] xml XML element. | 106 | * @param[in] xml XML element. |
| 105 | * | 107 | * |
| 106 | - * Creates instance of a spectral function model by extracting information | ||
| 107 | - * from an XML element. See GModelSpectralFunc::read() for more information | ||
| 108 | - * about the expected structure of the XML element. | 108 | + * Constructs spectral file function model by extracting information from an |
| 109 | + * XML element. See the read() method for more information about the expected | ||
| 110 | + * structure of the XML element. | ||
| 109 | ***************************************************************************/ | 111 | ***************************************************************************/ |
| 110 | GModelSpectralFunc::GModelSpectralFunc(const GXmlElement& xml) : | 112 | GModelSpectralFunc::GModelSpectralFunc(const GXmlElement& xml) : |
| 111 | GModelSpectral() | 113 | GModelSpectral() |
| @@ -124,7 +126,7 @@ GModelSpectralFunc::GModelSpectralFunc(const GXmlElement& xml) : | @@ -124,7 +126,7 @@ GModelSpectralFunc::GModelSpectralFunc(const GXmlElement& xml) : | ||
| 124 | /***********************************************************************//** | 126 | /***********************************************************************//** |
| 125 | * @brief Copy constructor | 127 | * @brief Copy constructor |
| 126 | * | 128 | * |
| 127 | - * @param[in] model Spectral function model. | 129 | + * @param[in] model File function model. |
| 128 | ***************************************************************************/ | 130 | ***************************************************************************/ |
| 129 | GModelSpectralFunc::GModelSpectralFunc(const GModelSpectralFunc& model) : | 131 | GModelSpectralFunc::GModelSpectralFunc(const GModelSpectralFunc& model) : |
| 130 | GModelSpectral(model) | 132 | GModelSpectral(model) |
| @@ -162,7 +164,8 @@ GModelSpectralFunc::~GModelSpectralFunc(void) | @@ -162,7 +164,8 @@ GModelSpectralFunc::~GModelSpectralFunc(void) | ||
| 162 | /***********************************************************************//** | 164 | /***********************************************************************//** |
| 163 | * @brief Assignment operator | 165 | * @brief Assignment operator |
| 164 | * | 166 | * |
| 165 | - * @param[in] model Spectral function model. | 167 | + * @param[in] model File function model. |
| 168 | + * @return File function model. | ||
| 166 | ***************************************************************************/ | 169 | ***************************************************************************/ |
| 167 | GModelSpectralFunc& GModelSpectralFunc::operator= (const GModelSpectralFunc& model) | 170 | GModelSpectralFunc& GModelSpectralFunc::operator= (const GModelSpectralFunc& model) |
| 168 | { | 171 | { |
| @@ -195,7 +198,7 @@ GModelSpectralFunc& GModelSpectralFunc::operator= (const GModelSpectralFunc& mod | @@ -195,7 +198,7 @@ GModelSpectralFunc& GModelSpectralFunc::operator= (const GModelSpectralFunc& mod | ||
| 195 | ==========================================================================*/ | 198 | ==========================================================================*/ |
| 196 | 199 | ||
| 197 | /***********************************************************************//** | 200 | /***********************************************************************//** |
| 198 | - * @brief Clear instance | 201 | + * @brief Clear file function |
| 199 | ***************************************************************************/ | 202 | ***************************************************************************/ |
| 200 | void GModelSpectralFunc::clear(void) | 203 | void GModelSpectralFunc::clear(void) |
| 201 | { | 204 | { |
| @@ -213,27 +216,30 @@ void GModelSpectralFunc::clear(void) | @@ -213,27 +216,30 @@ void GModelSpectralFunc::clear(void) | ||
| 213 | 216 | ||
| 214 | 217 | ||
| 215 | /***********************************************************************//** | 218 | /***********************************************************************//** |
| 216 | - * @brief Clone instance | 219 | + * @brief Clone file function |
| 217 | ***************************************************************************/ | 220 | ***************************************************************************/ |
| 218 | GModelSpectralFunc* GModelSpectralFunc::clone(void) const | 221 | GModelSpectralFunc* GModelSpectralFunc::clone(void) const |
| 219 | { | 222 | { |
| 223 | + // Clone file function | ||
| 220 | return new GModelSpectralFunc(*this); | 224 | return new GModelSpectralFunc(*this); |
| 221 | } | 225 | } |
| 222 | 226 | ||
| 223 | 227 | ||
| 224 | /***********************************************************************//** | 228 | /***********************************************************************//** |
| 225 | - * @brief Evaluate model | 229 | + * @brief Evaluate function |
| 226 | * | 230 | * |
| 227 | * @param[in] srcEng True photon energy. | 231 | * @param[in] srcEng True photon energy. |
| 228 | * @param[in] srcTime True photon arrival time. | 232 | * @param[in] srcTime True photon arrival time. |
| 229 | * @return Model value (ph/cm2/s/MeV). | 233 | * @return Model value (ph/cm2/s/MeV). |
| 230 | * | 234 | * |
| 231 | - * The spectral model is defined as | ||
| 232 | - * \f[I(E)=norm f(E)\f] | 235 | + * Evaluates |
| 236 | + * | ||
| 237 | + * \f[ | ||
| 238 | + * S_{\rm E}(E | t) = {\tt m\_norm} | ||
| 239 | + * \f] | ||
| 240 | + * | ||
| 233 | * where | 241 | * where |
| 234 | - * \f$norm\f$ is the normalization of the function. | ||
| 235 | - * Note that the node energies are stored as log10 of energy in units of | ||
| 236 | - * MeV. | 242 | + * - \f${\tt m\_norm}\f$ is the normalization factor. |
| 237 | ***************************************************************************/ | 243 | ***************************************************************************/ |
| 238 | double GModelSpectralFunc::eval(const GEnergy& srcEng, | 244 | double GModelSpectralFunc::eval(const GEnergy& srcEng, |
| 239 | const GTime& srcTime) const | 245 | const GTime& srcTime) const |
| @@ -244,13 +250,14 @@ double GModelSpectralFunc::eval(const GEnergy& srcEng, | @@ -244,13 +250,14 @@ double GModelSpectralFunc::eval(const GEnergy& srcEng, | ||
| 244 | double func = std::pow(10.0, arg); | 250 | double func = std::pow(10.0, arg); |
| 245 | 251 | ||
| 246 | // Compute function value | 252 | // Compute function value |
| 247 | - double value = norm() * func; | 253 | + double value = m_norm.value() * func; |
| 248 | 254 | ||
| 249 | // Compile option: Check for NaN/Inf | 255 | // Compile option: Check for NaN/Inf |
| 250 | #if defined(G_NAN_CHECK) | 256 | #if defined(G_NAN_CHECK) |
| 251 | if (isnotanumber(value) || isinfinite(value)) { | 257 | if (isnotanumber(value) || isinfinite(value)) { |
| 252 | std::cout << "*** ERROR: GModelSpectralFunc::eval"; | 258 | std::cout << "*** ERROR: GModelSpectralFunc::eval"; |
| 253 | - std::cout << "(srcEng=" << srcEng << "):"; | 259 | + std::cout << "(srcEng=" << srcEng; |
| 260 | + std::cout << ", srcTime=" << srcTime << "):"; | ||
| 254 | std::cout << " NaN/Inf encountered"; | 261 | std::cout << " NaN/Inf encountered"; |
| 255 | std::cout << " (value=" << value; | 262 | std::cout << " (value=" << value; |
| 256 | std::cout << ", norm=" << norm(); | 263 | std::cout << ", norm=" << norm(); |
| @@ -271,16 +278,22 @@ double GModelSpectralFunc::eval(const GEnergy& srcEng, | @@ -271,16 +278,22 @@ double GModelSpectralFunc::eval(const GEnergy& srcEng, | ||
| 271 | * @param[in] srcTime True photon arrival time. | 278 | * @param[in] srcTime True photon arrival time. |
| 272 | * @return Model value (ph/cm2/s/MeV). | 279 | * @return Model value (ph/cm2/s/MeV). |
| 273 | * | 280 | * |
| 274 | - * The spectral model is defined as | ||
| 275 | - * \f[I(E)=norm f(E)\f] | 281 | + * Evaluates |
| 282 | + * | ||
| 283 | + * \f[ | ||
| 284 | + * S_{\rm E}(E | t) = {\tt m\_norm} | ||
| 285 | + * \f] | ||
| 286 | + * | ||
| 276 | * where | 287 | * where |
| 277 | - * \f$norm=n_s n_v\f$ is the normalization of the function. | ||
| 278 | - * Note that the normalization is factorised into a scaling factor and a | ||
| 279 | - * value and that the method is expected to return the gradient with respect | ||
| 280 | - * to the parameter value \f$n_v\f$. | 288 | + * - \f${\tt m\_norm}\f$ is the normalization factor. |
| 289 | + * | ||
| 290 | + * The method also evaluates the partial derivatives of the model with | ||
| 291 | + * respect to the normalization parameter using | ||
| 281 | * | 292 | * |
| 282 | - * The partial derivative of the normalization value is given by | ||
| 283 | - * \f[dI/dn_v=n_s f(E)\f] | 293 | + * \f[ |
| 294 | + * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_norm}} = | ||
| 295 | + * \frac{S_{\rm E}(E | t)}{{\tt m\_norm}} | ||
| 296 | + * \f] | ||
| 284 | ***************************************************************************/ | 297 | ***************************************************************************/ |
| 285 | double GModelSpectralFunc::eval_gradients(const GEnergy& srcEng, | 298 | double GModelSpectralFunc::eval_gradients(const GEnergy& srcEng, |
| 286 | const GTime& srcTime) | 299 | const GTime& srcTime) |
| @@ -291,19 +304,20 @@ double GModelSpectralFunc::eval_gradients(const GEnergy& srcEng, | @@ -291,19 +304,20 @@ double GModelSpectralFunc::eval_gradients(const GEnergy& srcEng, | ||
| 291 | double func = std::pow(10.0, arg); | 304 | double func = std::pow(10.0, arg); |
| 292 | 305 | ||
| 293 | // Compute function value | 306 | // Compute function value |
| 294 | - double value = norm() * func; | 307 | + double value = m_norm.value() * func; |
| 295 | 308 | ||
| 296 | // Compute partial derivatives of the parameter values | 309 | // Compute partial derivatives of the parameter values |
| 297 | double g_norm = (m_norm.isfree()) ? m_norm.scale() * func : 0.0; | 310 | double g_norm = (m_norm.isfree()) ? m_norm.scale() * func : 0.0; |
| 298 | 311 | ||
| 299 | - // Set gradients (circumvent const correctness) | ||
| 300 | - const_cast<GModelSpectralFunc*>(this)->m_norm.factor_gradient(g_norm); | 312 | + // Set gradients |
| 313 | + m_norm.factor_gradient(g_norm); | ||
| 301 | 314 | ||
| 302 | // Compile option: Check for NaN/Inf | 315 | // Compile option: Check for NaN/Inf |
| 303 | #if defined(G_NAN_CHECK) | 316 | #if defined(G_NAN_CHECK) |
| 304 | if (isnotanumber(value) || isinfinite(value)) { | 317 | if (isnotanumber(value) || isinfinite(value)) { |
| 305 | std::cout << "*** ERROR: GModelSpectralFunc::eval_gradients"; | 318 | std::cout << "*** ERROR: GModelSpectralFunc::eval_gradients"; |
| 306 | - std::cout << "(srcEng=" << srcEng << "):"; | 319 | + std::cout << "(srcEng=" << srcEng; |
| 320 | + std::cout << ", srcTime=" << srcTime << "):"; | ||
| 307 | std::cout << " NaN/Inf encountered"; | 321 | std::cout << " NaN/Inf encountered"; |
| 308 | std::cout << " (value=" << value; | 322 | std::cout << " (value=" << value; |
| 309 | std::cout << ", norm=" << norm(); | 323 | std::cout << ", norm=" << norm(); |
| @@ -325,34 +339,27 @@ double GModelSpectralFunc::eval_gradients(const GEnergy& srcEng, | @@ -325,34 +339,27 @@ double GModelSpectralFunc::eval_gradients(const GEnergy& srcEng, | ||
| 325 | * @param[in] emax Maximum photon energy. | 339 | * @param[in] emax Maximum photon energy. |
| 326 | * @return Photon flux (ph/cm2/s). | 340 | * @return Photon flux (ph/cm2/s). |
| 327 | * | 341 | * |
| 328 | - * @exception GException::erange_invalid | ||
| 329 | - * Energy range is invalid (emin < emax required). | ||
| 330 | - * | ||
| 331 | * Computes | 342 | * Computes |
| 332 | - * \f[\int_{E_{\rm min}}^{E_{\rm max}} I(E) dE\f] | 343 | + * |
| 344 | + * \f[ | ||
| 345 | + * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) dE | ||
| 346 | + * \f] | ||
| 347 | + * | ||
| 333 | * where | 348 | * where |
| 334 | - * \f$E_{\rm min}\f$ and \f$E_{\rm max}\f$ are the minimum and maximum | ||
| 335 | - * energy, respectively, and | ||
| 336 | - * \f$I(E)\f$ is the spectral model (units: ph/cm2/s/MeV). | 349 | + * - [@p emin, @p emax] is an energy interval, and |
| 350 | + * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV). | ||
| 337 | ***************************************************************************/ | 351 | ***************************************************************************/ |
| 338 | double GModelSpectralFunc::flux(const GEnergy& emin, const GEnergy& emax) const | 352 | double GModelSpectralFunc::flux(const GEnergy& emin, const GEnergy& emax) const |
| 339 | { | 353 | { |
| 340 | - // Throw an exception if energy range is invalid | ||
| 341 | - if (emin >= emax) { | ||
| 342 | - throw GException::erange_invalid(G_FLUX, emin.MeV(), emax.MeV(), | ||
| 343 | - "Minimum energy < maximum energy required."); | ||
| 344 | - | ||
| 345 | - } | ||
| 346 | - | ||
| 347 | // Initialise flux | 354 | // Initialise flux |
| 348 | double flux = 0.0; | 355 | double flux = 0.0; |
| 349 | - | ||
| 350 | - // Get energy range in MeV | ||
| 351 | - double e_min = emin.MeV(); | ||
| 352 | - double e_max = emax.MeV(); | ||
| 353 | - | ||
| 354 | - // Continue only if e_max > e_min | ||
| 355 | - if (e_max > e_min) { | 356 | + |
| 357 | + // Compute only if integration range is valid | ||
| 358 | + if (emin < emax) { | ||
| 359 | + | ||
| 360 | + // Get energy range in MeV | ||
| 361 | + double e_min = emin.MeV(); | ||
| 362 | + double e_max = emax.MeV(); | ||
| 356 | 363 | ||
| 357 | // Determine left node index for minimum energy | 364 | // Determine left node index for minimum energy |
| 358 | m_lin_nodes.set_value(e_min); | 365 | m_lin_nodes.set_value(e_min); |
| @@ -404,9 +411,9 @@ double GModelSpectralFunc::flux(const GEnergy& emin, const GEnergy& emax) const | @@ -404,9 +411,9 @@ double GModelSpectralFunc::flux(const GEnergy& emin, const GEnergy& emax) const | ||
| 404 | } // endelse: emin and emax not between same nodes | 411 | } // endelse: emin and emax not between same nodes |
| 405 | 412 | ||
| 406 | // Multiply flux by normalisation factor | 413 | // Multiply flux by normalisation factor |
| 407 | - flux *= norm(); | 414 | + flux *= m_norm.value(); |
| 408 | 415 | ||
| 409 | - } // endif: e_max > e_min | 416 | + } // endif: integration range was valid |
| 410 | 417 | ||
| 411 | // Return | 418 | // Return |
| 412 | return flux; | 419 | return flux; |
| @@ -420,34 +427,27 @@ double GModelSpectralFunc::flux(const GEnergy& emin, const GEnergy& emax) const | @@ -420,34 +427,27 @@ double GModelSpectralFunc::flux(const GEnergy& emin, const GEnergy& emax) const | ||
| 420 | * @param[in] emax Maximum photon energy. | 427 | * @param[in] emax Maximum photon energy. |
| 421 | * @return Energy flux (erg/cm2/s). | 428 | * @return Energy flux (erg/cm2/s). |
| 422 | * | 429 | * |
| 423 | - * @exception GException::erange_invalid | ||
| 424 | - * Energy range is invalid (emin < emax required). | ||
| 425 | - * | ||
| 426 | * Computes | 430 | * Computes |
| 427 | - * \f[\int_{E_{\rm min}}^{E_{\rm max}} I(E) E dE\f] | 431 | + * |
| 432 | + * \f[ | ||
| 433 | + * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) E \, dE | ||
| 434 | + * \f] | ||
| 435 | + * | ||
| 428 | * where | 436 | * where |
| 429 | - * \f$E_{\rm min}\f$ and \f$E_{\rm max}\f$ are the minimum and maximum | ||
| 430 | - * energy, respectively, and | ||
| 431 | - * \f$I(E)\f$ is the spectral model (units: ph/cm2/s/MeV). | 437 | + * - [@p emin, @p emax] is an energy interval, and |
| 438 | + * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV). | ||
| 432 | ***************************************************************************/ | 439 | ***************************************************************************/ |
| 433 | double GModelSpectralFunc::eflux(const GEnergy& emin, const GEnergy& emax) const | 440 | double GModelSpectralFunc::eflux(const GEnergy& emin, const GEnergy& emax) const |
| 434 | { | 441 | { |
| 435 | - // Throw an exception if energy range is invalid | ||
| 436 | - if (emin >= emax) { | ||
| 437 | - throw GException::erange_invalid(G_EFLUX, emin.MeV(), emax.MeV(), | ||
| 438 | - "Minimum energy < maximum energy required."); | ||
| 439 | - | ||
| 440 | - } | ||
| 441 | - | ||
| 442 | // Initialise flux | 442 | // Initialise flux |
| 443 | - double flux = 0.0; | ||
| 444 | - | ||
| 445 | - // Get energy range in MeV | ||
| 446 | - double e_min = emin.MeV(); | ||
| 447 | - double e_max = emax.MeV(); | 443 | + double eflux = 0.0; |
| 444 | + | ||
| 445 | + // Compute only if integration range is valid | ||
| 446 | + if (emin < emax) { | ||
| 448 | 447 | ||
| 449 | - // Continue only if e_max > e_min | ||
| 450 | - if (e_max > e_min) { | 448 | + // Get energy range in MeV |
| 449 | + double e_min = emin.MeV(); | ||
| 450 | + double e_max = emax.MeV(); | ||
| 451 | 451 | ||
| 452 | // Determine left node index for minimum energy | 452 | // Determine left node index for minimum energy |
| 453 | m_lin_nodes.set_value(e_min); | 453 | m_lin_nodes.set_value(e_min); |
| @@ -461,11 +461,11 @@ double GModelSpectralFunc::eflux(const GEnergy& emin, const GEnergy& emax) const | @@ -461,11 +461,11 @@ double GModelSpectralFunc::eflux(const GEnergy& emin, const GEnergy& emax) const | ||
| 461 | // integrate over the energy interval using the appropriate power | 461 | // integrate over the energy interval using the appropriate power |
| 462 | // law parameters | 462 | // law parameters |
| 463 | if (inx_emin == inx_emax) { | 463 | if (inx_emin == inx_emax) { |
| 464 | - flux = m_prefactor[inx_emin] * | ||
| 465 | - plaw_energy_flux(e_min, | ||
| 466 | - e_max, | ||
| 467 | - m_epivot[inx_emin], | ||
| 468 | - m_gamma[inx_emin]) * MeV2erg; | 464 | + eflux = m_prefactor[inx_emin] * |
| 465 | + plaw_energy_flux(e_min, | ||
| 466 | + e_max, | ||
| 467 | + m_epivot[inx_emin], | ||
| 468 | + m_gamma[inx_emin]) * MeV2erg; | ||
| 469 | } | 469 | } |
| 470 | 470 | ||
| 471 | // ... otherwise integrate over the nodes where emin and emax | 471 | // ... otherwise integrate over the nodes where emin and emax |
| @@ -478,33 +478,33 @@ double GModelSpectralFunc::eflux(const GEnergy& emin, const GEnergy& emax) const | @@ -478,33 +478,33 @@ double GModelSpectralFunc::eflux(const GEnergy& emin, const GEnergy& emax) const | ||
| 478 | int i_start = (e_min < m_lin_nodes[0]) ? inx_emin : inx_emin+1; | 478 | int i_start = (e_min < m_lin_nodes[0]) ? inx_emin : inx_emin+1; |
| 479 | 479 | ||
| 480 | // Integrate from emin to the node boundary | 480 | // Integrate from emin to the node boundary |
| 481 | - flux = m_prefactor[inx_emin] * | ||
| 482 | - plaw_energy_flux(e_min, | ||
| 483 | - m_lin_nodes[i_start], | ||
| 484 | - m_epivot[inx_emin], | ||
| 485 | - m_gamma[inx_emin]) * MeV2erg; | 481 | + eflux = m_prefactor[inx_emin] * |
| 482 | + plaw_energy_flux(e_min, | ||
| 483 | + m_lin_nodes[i_start], | ||
| 484 | + m_epivot[inx_emin], | ||
| 485 | + m_gamma[inx_emin]) * MeV2erg; | ||
| 486 | 486 | ||
| 487 | // Integrate over all nodes between | 487 | // Integrate over all nodes between |
| 488 | for (int i = i_start; i < inx_emax; ++i) { | 488 | for (int i = i_start; i < inx_emax; ++i) { |
| 489 | - flux += m_eflux[i]; | 489 | + eflux += m_eflux[i]; |
| 490 | } | 490 | } |
| 491 | 491 | ||
| 492 | // Integrate from node boundary to emax | 492 | // Integrate from node boundary to emax |
| 493 | - flux += m_prefactor[inx_emax] * | ||
| 494 | - plaw_energy_flux(m_lin_nodes[inx_emax], | ||
| 495 | - e_max, | ||
| 496 | - m_epivot[inx_emax], | ||
| 497 | - m_gamma[inx_emax]) * MeV2erg; | 493 | + eflux += m_prefactor[inx_emax] * |
| 494 | + plaw_energy_flux(m_lin_nodes[inx_emax], | ||
| 495 | + e_max, | ||
| 496 | + m_epivot[inx_emax], | ||
| 497 | + m_gamma[inx_emax]) * MeV2erg; | ||
| 498 | 498 | ||
| 499 | } // endelse: emin and emax not between same nodes | 499 | } // endelse: emin and emax not between same nodes |
| 500 | 500 | ||
| 501 | // Multiply flux by normalisation factor | 501 | // Multiply flux by normalisation factor |
| 502 | - flux *= norm(); | 502 | + eflux *= norm(); |
| 503 | 503 | ||
| 504 | - } // endif: e_max > e_min | 504 | + } // endif: integration range was valid |
| 505 | 505 | ||
| 506 | - // Return | ||
| 507 | - return flux; | 506 | + // Return flux |
| 507 | + return eflux; | ||
| 508 | } | 508 | } |
| 509 | 509 | ||
| 510 | 510 | ||
| @@ -520,8 +520,8 @@ double GModelSpectralFunc::eflux(const GEnergy& emin, const GEnergy& emax) const | @@ -520,8 +520,8 @@ double GModelSpectralFunc::eflux(const GEnergy& emin, const GEnergy& emax) const | ||
| 520 | * @exception GException::erange_invalid | 520 | * @exception GException::erange_invalid |
| 521 | * Energy range is invalid (emin < emax required). | 521 | * Energy range is invalid (emin < emax required). |
| 522 | * | 522 | * |
| 523 | - * Simulates a random energy in the interval [emin, emax] for a spectral | ||
| 524 | - * function. | 523 | + * Returns Monte Carlo energy by randomly drawing from a broken power law |
| 524 | + * defined by the file function. | ||
| 525 | ***************************************************************************/ | 525 | ***************************************************************************/ |
| 526 | GEnergy GModelSpectralFunc::mc(const GEnergy& emin, | 526 | GEnergy GModelSpectralFunc::mc(const GEnergy& emin, |
| 527 | const GEnergy& emax, | 527 | const GEnergy& emax, |
| @@ -537,42 +537,37 @@ GEnergy GModelSpectralFunc::mc(const GEnergy& emin, | @@ -537,42 +537,37 @@ GEnergy GModelSpectralFunc::mc(const GEnergy& emin, | ||
| 537 | // Allocate energy | 537 | // Allocate energy |
| 538 | GEnergy energy; | 538 | GEnergy energy; |
| 539 | 539 | ||
| 540 | - // Continue only if emax > emin | ||
| 541 | - if (emax > emin) { | ||
| 542 | - | ||
| 543 | - // Update cache | ||
| 544 | - mc_update(emin, emax); | ||
| 545 | - | ||
| 546 | - // Determine in which bin we reside | ||
| 547 | - int inx = 0; | ||
| 548 | - if (m_mc_cum.size() > 1) { | ||
| 549 | - double u = ran.uniform(); | ||
| 550 | - for (inx = m_mc_cum.size()-1; inx > 0; --inx) { | ||
| 551 | - if (m_mc_cum[inx-1] <= u) { | ||
| 552 | - break; | ||
| 553 | - } | 540 | + // Update cache |
| 541 | + mc_update(emin, emax); | ||
| 542 | + | ||
| 543 | + // Determine in which bin we reside | ||
| 544 | + int inx = 0; | ||
| 545 | + if (m_mc_cum.size() > 1) { | ||
| 546 | + double u = ran.uniform(); | ||
| 547 | + for (inx = m_mc_cum.size()-1; inx > 0; --inx) { | ||
| 548 | + if (m_mc_cum[inx-1] <= u) { | ||
| 549 | + break; | ||
| 554 | } | 550 | } |
| 555 | } | 551 | } |
| 552 | + } | ||
| 556 | 553 | ||
| 557 | - // Get random energy for specific bin | ||
| 558 | - if (m_mc_exp[inx] != 0.0) { | ||
| 559 | - double e_min = m_mc_min[inx]; | ||
| 560 | - double e_max = m_mc_max[inx]; | ||
| 561 | - double u = ran.uniform(); | ||
| 562 | - double eng = (u > 0.0) | ||
| 563 | - ? std::exp(std::log(u * (e_max - e_min) + e_min) / m_mc_exp[inx]) | ||
| 564 | - : 0.0; | ||
| 565 | - energy.MeV(eng); | ||
| 566 | - } | ||
| 567 | - else { | ||
| 568 | - double e_min = m_mc_min[inx]; | ||
| 569 | - double e_max = m_mc_max[inx]; | ||
| 570 | - double u = ran.uniform(); | ||
| 571 | - double eng = std::exp(u * (e_max - e_min) + e_min); | ||
| 572 | - energy.MeV(eng); | ||
| 573 | - } | ||
| 574 | - | ||
| 575 | - } // endif: emax > emin | 554 | + // Get random energy for specific bin |
| 555 | + if (m_mc_exp[inx] != 0.0) { | ||
| 556 | + double e_min = m_mc_min[inx]; | ||
| 557 | + double e_max = m_mc_max[inx]; | ||
| 558 | + double u = ran.uniform(); | ||
| 559 | + double eng = (u > 0.0) | ||
| 560 | + ? std::exp(std::log(u * (e_max - e_min) + e_min) / m_mc_exp[inx]) | ||
| 561 | + : 0.0; | ||
| 562 | + energy.MeV(eng); | ||
| 563 | + } | ||
| 564 | + else { | ||
| 565 | + double e_min = m_mc_min[inx]; | ||
| 566 | + double e_max = m_mc_max[inx]; | ||
| 567 | + double u = ran.uniform(); | ||
| 568 | + double eng = std::exp(u * (e_max - e_min) + e_min); | ||
| 569 | + energy.MeV(eng); | ||
| 570 | + } | ||
| 576 | 571 | ||
| 577 | // Return energy | 572 | // Return energy |
| 578 | return energy; | 573 | return energy; |
| @@ -589,10 +584,12 @@ GEnergy GModelSpectralFunc::mc(const GEnergy& emin, | @@ -589,10 +584,12 @@ GEnergy GModelSpectralFunc::mc(const GEnergy& emin, | ||
| 589 | * @exception GException::model_invalid_parnames | 584 | * @exception GException::model_invalid_parnames |
| 590 | * Invalid model parameter name found in XML element. | 585 | * Invalid model parameter name found in XML element. |
| 591 | * | 586 | * |
| 592 | - * Read the function information from an XML element and load the nodes | ||
| 593 | - * from the associated file. The XML element is required to have an | ||
| 594 | - * attribute "file" that specifies the function nodes and a parameter | ||
| 595 | - * named "Normalization". | 587 | + * Reads the spectral information from an XML element. The format of the XML |
| 588 | + * elements is | ||
| 589 | + * | ||
| 590 | + * <spectrum type="FileFunction" file=".."> | ||
| 591 | + * <parameter name="Normalization" scale=".." value=".." min=".." max=".." free=".."/> | ||
| 592 | + * </spectrum> | ||
| 596 | ***************************************************************************/ | 593 | ***************************************************************************/ |
| 597 | void GModelSpectralFunc::read(const GXmlElement& xml) | 594 | void GModelSpectralFunc::read(const GXmlElement& xml) |
| 598 | { | 595 | { |
| @@ -634,10 +631,15 @@ void GModelSpectralFunc::read(const GXmlElement& xml) | @@ -634,10 +631,15 @@ void GModelSpectralFunc::read(const GXmlElement& xml) | ||
| 634 | * @exception GException::model_invalid_parnames | 631 | * @exception GException::model_invalid_parnames |
| 635 | * Invalid model parameter names found in XML element. | 632 | * Invalid model parameter names found in XML element. |
| 636 | * | 633 | * |
| 637 | - * Write the spectral function information into an XML element. The XML | ||
| 638 | - * element has to be of type "FileFunction" and will have 1 parameter leaf | ||
| 639 | - * named "Normalization". Note that the function nodes will not be written | ||
| 640 | - * since they will not be altered by any method. | 634 | + * Writes the spectral information into an XML element. The format of the XML |
| 635 | + * element is | ||
| 636 | + * | ||
| 637 | + * <spectrum type="FileFunction" file=".."> | ||
| 638 | + * <parameter name="Normalization" scale=".." value=".." min=".." max=".." free=".."/> | ||
| 639 | + * </spectrum> | ||
| 640 | + * | ||
| 641 | + * Note that the function nodes will not be written since they will not be | ||
| 642 | + * altered by any method. | ||
| 641 | ***************************************************************************/ | 643 | ***************************************************************************/ |
| 642 | void GModelSpectralFunc::write(GXmlElement& xml) const | 644 | void GModelSpectralFunc::write(GXmlElement& xml) const |
| 643 | { | 645 | { |
| @@ -681,7 +683,9 @@ void GModelSpectralFunc::write(GXmlElement& xml) const | @@ -681,7 +683,9 @@ void GModelSpectralFunc::write(GXmlElement& xml) const | ||
| 681 | 683 | ||
| 682 | 684 | ||
| 683 | /***********************************************************************//** | 685 | /***********************************************************************//** |
| 684 | - * @brief Print powerlaw information | 686 | + * @brief Print file function information |
| 687 | + * | ||
| 688 | + * @return String containing file function information. | ||
| 685 | ***************************************************************************/ | 689 | ***************************************************************************/ |
| 686 | std::string GModelSpectralFunc::print(void) const | 690 | std::string GModelSpectralFunc::print(void) const |
| 687 | { | 691 | { |
| @@ -690,6 +694,8 @@ std::string GModelSpectralFunc::print(void) const | @@ -690,6 +694,8 @@ std::string GModelSpectralFunc::print(void) const | ||
| 690 | 694 | ||
| 691 | // Append header | 695 | // Append header |
| 692 | result.append("=== GModelSpectralFunc ==="); | 696 | result.append("=== GModelSpectralFunc ==="); |
| 697 | + | ||
| 698 | + // Append information | ||
| 693 | result.append("\n"+parformat("Function file")+m_filename); | 699 | result.append("\n"+parformat("Function file")+m_filename); |
| 694 | result.append("\n"+parformat("Number of nodes")+str(m_lin_nodes.size())); | 700 | result.append("\n"+parformat("Number of nodes")+str(m_lin_nodes.size())); |
| 695 | result.append("\n"+parformat("Number of parameters")+str(size())); | 701 | result.append("\n"+parformat("Number of parameters")+str(size())); |
| @@ -962,9 +968,9 @@ void GModelSpectralFunc::set_cache(void) const | @@ -962,9 +968,9 @@ void GModelSpectralFunc::set_cache(void) const | ||
| 962 | * | 968 | * |
| 963 | * This method sets up an array of indices and the cumulative distribution | 969 | * This method sets up an array of indices and the cumulative distribution |
| 964 | * function needed for MC simulations. | 970 | * function needed for MC simulations. |
| 965 | - * | ||
| 966 | ***************************************************************************/ | 971 | ***************************************************************************/ |
| 967 | -void GModelSpectralFunc::mc_update(const GEnergy& emin, const GEnergy& emax) const | 972 | +void GModelSpectralFunc::mc_update(const GEnergy& emin, |
| 973 | + const GEnergy& emax) const | ||
| 968 | { | 974 | { |
| 969 | // Check if we need to update the cache | 975 | // Check if we need to update the cache |
| 970 | if (emin != m_mc_emin || emax != m_mc_emax) { | 976 | if (emin != m_mc_emin || emax != m_mc_emax) { |
| @@ -1090,10 +1096,3 @@ void GModelSpectralFunc::mc_update(const GEnergy& emin, const GEnergy& emax) con | @@ -1090,10 +1096,3 @@ void GModelSpectralFunc::mc_update(const GEnergy& emin, const GEnergy& emax) con | ||
| 1090 | // Return | 1096 | // Return |
| 1091 | return; | 1097 | return; |
| 1092 | } | 1098 | } |
| 1093 | - | ||
| 1094 | - | ||
| 1095 | -/*========================================================================== | ||
| 1096 | - = = | ||
| 1097 | - = Friends = | ||
| 1098 | - = = | ||
| 1099 | - ==========================================================================*/ |