Justine Devin / gammalib

Commit 5cbe4da0b2559b6d8c6a4b4bb0899608b689bca0

Authored by Jürgen Knödlseder
1 parent 0b7ec5a0

Add node manipulation methods. 

A set of methods have been added for node manipulation, such as retrieval of node energy and intensity, but also appending, inserting and removing of nodes or extending the node function.
Inline Side-by-side
Showing 3 changed files with 563 additions and 176 deletions
include/GModelSpectralNodes.hpp
  View file @5cbe4da
... ... @@ -43,22 +43,37 @@
43 43 * @brief Spectral nodes model class
44 44 *
45 45 * This class implements a piecewise power law between spectral nodes
46   - * \f$(E_i, I_i)\f$,
  46 + *
  47 + * \f[
  48 + * ({\tt m\_energies[i]}, {\tt m\_values[i]})
  49 + * \f]
  50 + *
47 51 * where
48   - * \f$E_i\f$ is the energy and \f$I_i\f$ is the intensity of node \f$i\f$.
49   - * For a given energy \f$E\f$, the piecewise powerlaw is computing by
50   - * finding the bracketing energies \f$E_1 <= E <= E_2\f$ and computing
51   - * \f[I(E)=10^{(\log v_1 + \log s_1) w_1 + (\log v_2 + \log s_2) w_2}\f]
  52 + * - \f${\tt m\_energies[i]}\f$ is the energy, and
  53 + * - \f${\tt m\_values[i]}\f$ is the intensity (in photons/cm2/s/MeV)
  54 + * of node \f$i\f$.
  55 + *
  56 + * For a given energy \f$E\f$, the piecewise power law is computing by
  57 + * finding the bracketing energies
  58 + * \f${\tt m\_energies[i]} <= E <= {\tt m\_energies[i+1]}\f$ and computing
  59 + *
  60 + * \f[
  61 + * S_{\rm E}(E | t) =
  62 + * 10^{(\log {\tt m\_values[i]}) w_{i} +
  63 + * (\log {\tt m\_values[i+1]}) w_{i+1}}
  64 + * \f]
  65 + *
52 66 * where
53   - * \f$I_1 = v_1 s_1\f$ is the intensity of node 1 (\f$E_1\f$),
54   - * \f$I_2 = v_2 s_2\f$ is the intensity of node 2 (\f$E_2\f$),
55   - * \f$w_1\f$ is the weighting of node 1, and
56   - * \f$w_2\f$ is the weighting of node 2.
57   - * The weightings \f$w_1\f$ and \f$w_2\f$ are computed by linear
  67 + * - \f${\tt m\_values[i]}\f$ is the intensity of node \f$i\f$,
  68 + * - \f${\tt m\_values[i+1]}\f$ is the intensity of node \f$i+1\f$,
  69 + * - \f$w_{i}\f$ is the weighting of node \f$i\f$, and
  70 + * - \f$w_{i+1}\f$ is the weighting of node \f$i+1\f$.
  71 + *
  72 + * The weightings \f$w_{i}\f$ and \f$w_{i+1}\f$ are computed by linear
58 73 * interpolation (in the log-log plane) between the nodes
59   - * \f$(\log E_1, \log I_1)\f$
  74 + * \f$(\log {\tt m\_energies[i]}, \log{\tt m\_values[i]})\f$
60 75 * and
61   - * \f$(\log E_2, \log I_2)\f$
  76 + * \f$(\log {\tt m\_energies[i+1]}, \log{\tt m\_values[i+1]})\f$
62 77 * to the requested energy \f$\log E\f$.
63 78 ***************************************************************************/
64 79 class GModelSpectralNodes : public GModelSpectral {
... ... @@ -73,7 +88,7 @@ public:
73 88 // Operators
74 89 virtual GModelSpectralNodes& operator=(const GModelSpectralNodes& model);
75 90  
76   - // Implemented pure virtual methods
  91 + // Implemented pure virtual base class methods
77 92 virtual void clear(void);
78 93 virtual GModelSpectralNodes* clone(void) const;
79 94 virtual std::string type(void) const;
... ... @@ -93,6 +108,19 @@ public:
93 108 virtual void write(GXmlElement& xml) const;
94 109 virtual std::string print(void) const;
95 110  
  111 + // Other methods
  112 + int nodes(void) const;
  113 + void append(const GEnergy& energy, const double& intensity);
  114 + void insert(const int& index, const GEnergy& energy,
  115 + const double& intensity);
  116 + void remove(const int& index);
  117 + void reserve(const int& num);
  118 + void extend(const GModelSpectralNodes& nodes);
  119 + GEnergy energy(const int& index) const;
  120 + void energy(const int& index, const GEnergy& energy);
  121 + double intensity(const int& index) const;
  122 + void intensity(const int& index, const double& intensity);
  123 +
96 124 protected:
97 125 // Protected methods
98 126 void init_members(void);
... ... @@ -149,4 +177,18 @@ std::string GModelSpectralNodes::type(void) const
149 177 return "NodeFunction";
150 178 }
151 179  
  180 +
  181 +/***********************************************************************//**
  182 + * @brief Return number of nodes
  183 + *
  184 + * @return Number of nodes.
  185 + *
  186 + * Returns the number of nodes in the node function model.
  187 + ***************************************************************************/
  188 +inline
  189 +int GModelSpectralNodes::nodes(void) const
  190 +{
  191 + return (m_energies.size());
  192 +}
  193 +
152 194 #endif /* GMODELSPECTRALNODES_HPP */
... ...
pyext/GModelSpectralNodes.i
  View file @5cbe4da
... ... @@ -61,6 +61,19 @@ public:
61 61 GRan& ran) const;
62 62 virtual void read(const GXmlElement& xml);
63 63 virtual void write(GXmlElement& xml) const;
  64 +
  65 + // Other methods
  66 + int nodes(void) const;
  67 + void append(const GEnergy& energy, const double& intensity);
  68 + void insert(const int& index, const GEnergy& energy,
  69 + const double& intensity);
  70 + void remove(const int& index);
  71 + void reserve(const int& num);
  72 + void extend(const GModelSpectralNodes& nodes);
  73 + GEnergy energy(const int& index) const;
  74 + void energy(const int& index, const GEnergy& energy);
  75 + double intensity(const int& index) const;
  76 + void intensity(const int& index, const double& intensity);
64 77 };
65 78  
66 79  
... ...
src/model/GModelSpectralNodes.cpp
  View file @5cbe4da
... ... @@ -43,9 +43,15 @@ const GModelSpectralRegistry g_spectral_nodes_registry(&amp;g_spectral_nodes_seed);
43 43 /* __ Method name definitions ____________________________________________ */
44 44 #define G_FLUX "GModelSpectralNodes::flux(GEnergy&, GEnergy&)"
45 45 #define G_EFLUX "GModelSpectralNodes::eflux(GEnergy&, GEnergy&)"
46   -#define G_MC "GModelSpectralNodes::mc(GEnergy&, GEnergy&, GRan&)"
  46 +#define G_MC "GModelSpectralNodes::mc(GEnergy&, GEnergy&, GTime&, GRan&)"
47 47 #define G_READ "GModelSpectralNodes::read(GXmlElement&)"
48 48 #define G_WRITE "GModelSpectralNodes::write(GXmlElement&)"
  49 +#define G_INSERT "GModelSpectralNodes::insert(int&, GEnergy&, double&)"
  50 +#define G_REMOVE "GModelSpectralNodes::remove(int&)"
  51 +#define G_ENERGY_GET "GModelSpectralNodes::energy(int&)"
  52 +#define G_ENERGY_SET "GModelSpectralNodes::energy(int&, GEnergy&)"
  53 +#define G_INTENSITY_GET "GModelSpectralNodes::intensity(int&)"
  54 +#define G_INTENSITY_SET "GModelSpectralNodes::intensity(int&, double&)"
49 55  
50 56 /* __ Macros _____________________________________________________________ */
51 57  
... ... @@ -78,12 +84,12 @@ GModelSpectralNodes::GModelSpectralNodes(void) : GModelSpectral()
78 84 *
79 85 * @param[in] xml XML element.
80 86 *
81   - * Creates instance of a spectral nodes model by extracting information from
82   - * an XML element. See GModelSpectralNodes::read() for more information about
83   - * the expected structure of the XML element.
  87 + * Construct spectral nodes model by extracting information from an XML
  88 + * element. See the read() method for more information about the expected
  89 + * structure of the XML element.
84 90 ***************************************************************************/
85   -GModelSpectralNodes::GModelSpectralNodes(const GXmlElement& xml)
86   - : GModelSpectral()
  91 +GModelSpectralNodes::GModelSpectralNodes(const GXmlElement& xml) :
  92 + GModelSpectral()
87 93 {
88 94 // Initialise members
89 95 init_members();
... ... @@ -101,8 +107,8 @@ GModelSpectralNodes::GModelSpectralNodes(const GXmlElement&amp; xml)
101 107 *
102 108 * @param[in] model Spectral nodes model.
103 109 ***************************************************************************/
104   -GModelSpectralNodes::GModelSpectralNodes(const GModelSpectralNodes& model)
105   - : GModelSpectral(model)
  110 +GModelSpectralNodes::GModelSpectralNodes(const GModelSpectralNodes& model) :
  111 + GModelSpectral(model)
106 112 {
107 113 // Initialise members
108 114 init_members();
... ... @@ -138,8 +144,9 @@ GModelSpectralNodes::~GModelSpectralNodes(void)
138 144 * @brief Assignment operator
139 145 *
140 146 * @param[in] model Spectral nodes model.
  147 + * @return Spectral nodes model.
141 148 ***************************************************************************/
142   -GModelSpectralNodes& GModelSpectralNodes::operator= (const GModelSpectralNodes& model)
  149 +GModelSpectralNodes& GModelSpectralNodes::operator=(const GModelSpectralNodes& model)
143 150 {
144 151 // Execute only if object is not identical
145 152 if (this != &model) {
... ... @@ -170,7 +177,7 @@ GModelSpectralNodes&amp; GModelSpectralNodes::operator= (const GModelSpectralNodes&amp;
170 177 ==========================================================================*/
171 178  
172 179 /***********************************************************************//**
173   - * @brief Clear instance
  180 + * @brief Clear spectral nodes model
174 181 ***************************************************************************/
175 182 void GModelSpectralNodes::clear(void)
176 183 {
... ... @@ -188,10 +195,11 @@ void GModelSpectralNodes::clear(void)
188 195  
189 196  
190 197 /***********************************************************************//**
191   - * @brief Clone instance
  198 + * @brief Clone spectral nodes model
192 199 ***************************************************************************/
193 200 GModelSpectralNodes* GModelSpectralNodes::clone(void) const
194 201 {
  202 + // Clone spectral nodes model
195 203 return new GModelSpectralNodes(*this);
196 204 }
197 205  
... ... @@ -203,23 +211,26 @@ GModelSpectralNodes* GModelSpectralNodes::clone(void) const
203 211 * @param[in] srcTime True photon arrival time.
204 212 * @return Model value (ph/cm2/s/MeV).
205 213 *
206   - * The spectral model is defined as a piecewise power law between spectral
207   - * nodes \f$(E_i, I_i)\f$,
208   - * where
209   - * \f$E_i\f$ is the energy and \f$I_i\f$ is the intensity of node \f$i\f$.
210   - * For a given energy \f$E\f$, the piecewise powerlaw is computing by
211   - * finding the bracketing energies \f$E_1 <= E <= E_2\f$ and computing
212   - * \f[I(E)=10^{(\log v_1 + \log s_1) w_1 + (\log v_2 + \log s_2) w_2}\f]
  214 + * Computes
  215 + *
  216 + * \f[
  217 + * S_{\rm E}(E | t) =
  218 + * 10^{(\log {\tt m\_values[i]}) w_{i} +
  219 + * (\log {\tt m\_values[i+1]}) w_{i+1}}
  220 + * \f]
  221 + *
213 222 * where
214   - * \f$I_1 = v_1 s_1\f$ is the intensity of node 1 (\f$E_1\f$),
215   - * \f$I_2 = v_2 s_2\f$ is the intensity of node 2 (\f$E_2\f$),
216   - * \f$w_1\f$ is the weighting of node 1, and
217   - * \f$w_2\f$ is the weighting of node 2.
218   - * The weightings \f$w_1\f$ and \f$w_2\f$ are computed by linear
  223 + * - \f${\tt m\_values[i]}\f$ is the intensity of node \f$i\f$,
  224 + * - \f${\tt m\_values[i+1]}\f$ is the intensity of node \f$i+1\f$,
  225 + * - \f$w_{i}\f$ is the weighting of node \f$i\f$,
  226 + * - \f$w_{i+1}\f$ is the weighting of node \f$i+1\f$, and
  227 + * - \f${\tt m\_energies[i]} <= E <= {\tt m\_energies[i+1]}\f$.
  228 + *
  229 + * The weightings \f$w_{i}\f$ and \f$w_{i+1}\f$ are computed by linear
219 230 * interpolation (in the log-log plane) between the nodes
220   - * \f$(\log E_1, \log I_1)\f$
  231 + * \f$(\log {\tt m\_energies[i]}, \log{\tt m\_values[i]})\f$
221 232 * and
222   - * \f$(\log E_2, \log I_2)\f$
  233 + * \f$(\log {\tt m\_energies[i+1]}, \log{\tt m\_values[i+1]})\f$
223 234 * to the requested energy \f$\log E\f$.
224 235 ***************************************************************************/
225 236 double GModelSpectralNodes::eval(const GEnergy& srcEng,
... ... @@ -237,7 +248,8 @@ double GModelSpectralNodes::eval(const GEnergy&amp; srcEng,
237 248 #if defined(G_NAN_CHECK)
238 249 if (isnotanumber(value) || isinfinite(value)) {
239 250 std::cout << "*** ERROR: GModelSpectralNodes::eval";
240   - std::cout << "(srcEng=" << srcEng << "):";
  251 + std::cout << "(srcEng=" << srcEng;
  252 + std::cout << ", srcTime=" << srcTime << "):";
241 253 std::cout << " NaN/Inf encountered";
242 254 std::cout << " (value=" << value;
243 255 std::cout << ", exponent=" << exponent;
... ... @@ -257,28 +269,35 @@ double GModelSpectralNodes::eval(const GEnergy&amp; srcEng,
257 269 * @param[in] srcTime True photon arrival time.
258 270 * @return Model value (ph/cm2/s/MeV).
259 271 *
260   - * The spectral model is defined as a piecewise power law between spectral
261   - * nodes \f$(E_i, I_i)\f$,
262   - * where
263   - * \f$E_i\f$ is the energy and \f$I_i\f$ is the intensity of node \f$i\f$.
264   - * For a given energy \f$E\f$, the piecewise powerlaw is computing by
265   - * finding the bracketing energies \f$E_1 <= E <= E_2\f$ and computing
266   - * \f[I(E)=10^{(\log v_1 + \log s_1) w_1 + (\log v_2 + \log s_2) w_2}\f]
  272 + * Computes
  273 + *
  274 + * \f[
  275 + * S_{\rm E}(E | t) =
  276 + * 10^{(\log {\tt m\_values[i]}) w_{i} +
  277 + * (\log {\tt m\_values[i+1]}) w_{i+1}}
  278 + * \f]
  279 + *
267 280 * where
268   - * \f$I_1 = v_1 s_1\f$ is the intensity of node 1 (\f$E_1\f$),
269   - * \f$I_2 = v_2 s_2\f$ is the intensity of node 2 (\f$E_2\f$),
270   - * \f$w_1\f$ is the weighting of node 1, and
271   - * \f$w_2\f$ is the weighting of node 2.
272   - * The weightings \f$w_1\f$ and \f$w_2\f$ are computed by linear
  281 + * - \f${\tt m\_values[i]}\f$ is the intensity of node \f$i\f$,
  282 + * - \f${\tt m\_values[i+1]}\f$ is the intensity of node \f$i+1\f$,
  283 + * - \f$w_{i}\f$ is the weighting of node \f$i\f$,
  284 + * - \f$w_{i+1}\f$ is the weighting of node \f$i+1\f$, and
  285 + * - \f${\tt m\_energies[i]} <= E <= {\tt m\_energies[i+1]}\f$.
  286 + *
  287 + * The weightings \f$w_{i}\f$ and \f$w_{i+1}\f$ are computed by linear
273 288 * interpolation (in the log-log plane) between the nodes
274   - * \f$(\log E_1, \log I_1)\f$
  289 + * \f$(\log {\tt m\_energies[i]}, \log{\tt m\_values[i]})\f$
275 290 * and
276   - * \f$(\log E_2, \log I_2)\f$
  291 + * \f$(\log {\tt m\_energies[i+1]}, \log{\tt m\_values[i+1]})\f$
277 292 * to the requested energy \f$\log E\f$.
278 293 *
279   - * The partial derivatives are given by
280   - * \f[dI/dv_i=I(E) w_i / v_i\f]
281   - * for \f$i=1,2\f$ and \f$0\f$ for all other nodes.
  294 + * The method also evaluates the partial derivatives of the model with
  295 + * respect to the parameters using
  296 + *
  297 + * \f[
  298 + * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_values[i]}} =
  299 + * \frac{S_{\rm E}(E | t) \, w_i}{{\tt m\_values[i]}}
  300 + * \f]
282 301 ***************************************************************************/
283 302 double GModelSpectralNodes::eval_gradients(const GEnergy& srcEng,
284 303 const GTime& srcTime)
... ... @@ -304,26 +323,27 @@ double GModelSpectralNodes::eval_gradients(const GEnergy&amp; srcEng,
304 323  
305 324 // Initialise gradients
306 325 for (int i = 0; i < m_values.size(); ++i) {
307   - const_cast<GModelSpectralNodes*>(this)->m_values[i].factor_gradient(0.0);
  326 + m_values[i].factor_gradient(0.0);
308 327 }
309 328  
310 329 // Gradient for left node
311 330 if (m_values[inx_left].isfree()) {
312 331 double grad = value * wgt_left / m_values[inx_left].factor_value();
313   - const_cast<GModelSpectralNodes*>(this)->m_values[inx_left].factor_gradient(grad);
  332 + m_values[inx_left].factor_gradient(grad);
314 333 }
315 334  
316 335 // Gradient for right node
317 336 if (m_values[inx_right].isfree()) {
318 337 double grad = value * wgt_right / m_values[inx_right].factor_value();
319   - const_cast<GModelSpectralNodes*>(this)->m_values[inx_right].factor_gradient(grad);
  338 + m_values[inx_right].factor_gradient(grad);
320 339 }
321 340  
322 341 // Compile option: Check for NaN/Inf
323 342 #if defined(G_NAN_CHECK)
324 343 if (isnotanumber(value) || isinfinite(value)) {
325 344 std::cout << "*** ERROR: GModelSpectralNodes::eval_gradients";
326   - std::cout << "(srcEng=" << srcEng << "):";
  345 + std::cout << "(srcEng=" << srcEng;
  346 + std::cout << ", srcTime=" << srcTime << "):";
327 347 std::cout << " NaN/Inf encountered";
328 348 std::cout << " (value=" << value;
329 349 std::cout << ", exponent=" << exponent;
... ... @@ -343,37 +363,30 @@ double GModelSpectralNodes::eval_gradients(const GEnergy&amp; srcEng,
343 363 * @param[in] emax Maximum photon energy.
344 364 * @return Photon flux (ph/cm2/s).
345 365 *
346   - * @exception GException::erange_invalid
347   - * Energy range is invalid (emin < emax required).
348   - *
349 366 * Computes
350   - * \f[\int_{E_{\rm min}}^{E_{\rm max}} I(E) dE\f]
  367 + *
  368 + * \f[
  369 + * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) dE
  370 + * \f]
  371 + *
351 372 * where
352   - * \f$E_{\rm min}\f$ and \f$E_{\rm max}\f$ are the minimum and maximum
353   - * energy, respectively, and
354   - * \f$I(E)\f$ is the spectral model (units: ph/cm2/s/MeV).
  373 + * - [@p emin, @p emax] is an energy interval, and
  374 + * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
355 375 ***************************************************************************/
356 376 double GModelSpectralNodes::flux(const GEnergy& emin, const GEnergy& emax) const
357 377 {
358   - // Throw an exception if energy range is invalid
359   - if (emin >= emax) {
360   - throw GException::erange_invalid(G_FLUX, emin.MeV(), emax.MeV(),
361   - "Minimum energy < maximum energy required.");
362   -
363   - }
364   -
365   - // Update flux cache
366   - update_flux_cache();
367   -
368 378 // Initialise flux
369 379 double flux = 0.0;
370 380  
371   - // Get energy range in MeV
372   - double e_min = emin.MeV();
373   - double e_max = emax.MeV();
374   -
375   - // Continue only if e_max > e_min
376   - if (e_max > e_min) {
  381 + // Compute only if integration range is valid
  382 + if (emin < emax) {
  383 +
  384 + // Update flux cache
  385 + update_flux_cache();
  386 +
  387 + // Get energy range in MeV
  388 + double e_min = emin.MeV();
  389 + double e_max = emax.MeV();
377 390  
378 391 // Determine left node index for minimum energy
379 392 m_lin_energies.set_value(e_min);
... ... @@ -421,10 +434,9 @@ double GModelSpectralNodes::flux(const GEnergy&amp; emin, const GEnergy&amp; emax) const
421 434 e_max,
422 435 m_epivot[inx_emax],
423 436 m_gamma[inx_emax]);
424   -
425   - } // endelse: emin and emax not between same nodes
426   -
427   - } // endif: e_max > e_min
  437 + }
  438 +
  439 + } // endif: integration range was valid
428 440  
429 441 // Return
430 442 return flux;
... ... @@ -438,37 +450,30 @@ double GModelSpectralNodes::flux(const GEnergy&amp; emin, const GEnergy&amp; emax) const
438 450 * @param[in] emax Maximum photon energy.
439 451 * @return Energy flux (erg/cm2/s).
440 452 *
441   - * @exception GException::erange_invalid
442   - * Energy range is invalid (emin < emax required).
443   - *
444 453 * Computes
445   - * \f[\int_{E_{\rm min}}^{E_{\rm max}} I(E) E dE\f]
  454 + *
  455 + * \f[
  456 + * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) E \, dE
  457 + * \f]
  458 + *
446 459 * where
447   - * \f$E_{\rm min}\f$ and \f$E_{\rm max}\f$ are the minimum and maximum
448   - * energy, respectively, and
449   - * \f$I(E)\f$ is the spectral model (units: ph/cm2/s/MeV).
  460 + * - [@p emin, @p emax] is an energy interval, and
  461 + * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
450 462 ***************************************************************************/
451 463 double GModelSpectralNodes::eflux(const GEnergy& emin, const GEnergy& emax) const
452 464 {
453   - // Throw an exception if energy range is invalid
454   - if (emin >= emax) {
455   - throw GException::erange_invalid(G_EFLUX, emin.MeV(), emax.MeV(),
456   - "Minimum energy < maximum energy required.");
457   -
458   - }
459   -
460   - // Update flux cache
461   - update_flux_cache();
462   -
463 465 // Initialise flux
464   - double flux = 0.0;
  466 + double eflux = 0.0;
465 467  
466   - // Get energy range in MeV
467   - double e_min = emin.MeV();
468   - double e_max = emax.MeV();
  468 + // Compute only if integration range is valid
  469 + if (emin < emax) {
  470 +
  471 + // Update flux cache
  472 + update_flux_cache();
469 473  
470   - // Continue only if e_max > e_min
471   - if (e_max > e_min) {
  474 + // Get energy range in MeV
  475 + double e_min = emin.MeV();
  476 + double e_max = emax.MeV();
472 477  
473 478 // Determine left node index for minimum energy
474 479 m_lin_energies.set_value(e_min);
... ... @@ -482,11 +487,11 @@ double GModelSpectralNodes::eflux(const GEnergy&amp; emin, const GEnergy&amp; emax) cons
482 487 // integrate over the energy interval using the appropriate power
483 488 // law parameters
484 489 if (inx_emin == inx_emax) {
485   - flux = m_prefactor[inx_emin] *
486   - plaw_energy_flux(e_min,
487   - e_max,
488   - m_epivot[inx_emin],
489   - m_gamma[inx_emin]) * MeV2erg;
  490 + eflux = m_prefactor[inx_emin] *
  491 + plaw_energy_flux(e_min,
  492 + e_max,
  493 + m_epivot[inx_emin],
  494 + m_gamma[inx_emin]) * MeV2erg;
490 495 }
491 496  
492 497 // ... otherwise integrate over the nodes where emin and emax
... ... @@ -499,30 +504,30 @@ double GModelSpectralNodes::eflux(const GEnergy&amp; emin, const GEnergy&amp; emax) cons
499 504 int i_start = (e_min < m_lin_energies[0]) ? inx_emin : inx_emin+1;
500 505  
501 506 // Integrate from emin to the node boundary
502   - flux = m_prefactor[inx_emin] *
503   - plaw_energy_flux(e_min,
504   - m_lin_energies[i_start],
505   - m_epivot[inx_emin],
506   - m_gamma[inx_emin]) * MeV2erg;
  507 + eflux = m_prefactor[inx_emin] *
  508 + plaw_energy_flux(e_min,
  509 + m_lin_energies[i_start],
  510 + m_epivot[inx_emin],
  511 + m_gamma[inx_emin]) * MeV2erg;
507 512  
508 513 // Integrate over all nodes between
509 514 for (int i = i_start; i < inx_emax; ++i) {
510   - flux += m_eflux[i];
  515 + eflux += m_eflux[i];
511 516 }
512 517  
513 518 // Integrate from node boundary to emax
514   - flux += m_prefactor[inx_emax] *
515   - plaw_energy_flux(m_lin_energies[inx_emax],
516   - e_max,
517   - m_epivot[inx_emax],
518   - m_gamma[inx_emax]) * MeV2erg;
  519 + eflux += m_prefactor[inx_emax] *
  520 + plaw_energy_flux(m_lin_energies[inx_emax],
  521 + e_max,
  522 + m_epivot[inx_emax],
  523 + m_gamma[inx_emax]) * MeV2erg;
519 524  
520 525 } // endelse: emin and emax not between same nodes
521 526  
522   - } // endif: e_max > e_min
  527 + } // endif: integration range was valid
523 528  
524   - // Return
525   - return flux;
  529 + // Return flux
  530 + return eflux;
526 531 }
527 532  
528 533  
... ... @@ -537,6 +542,8 @@ double GModelSpectralNodes::eflux(const GEnergy&amp; emin, const GEnergy&amp; emax) cons
537 542 *
538 543 * @exception GException::erange_invalid
539 544 * Energy range is invalid (emin < emax required).
  545 + *
  546 + * Returns Monte Carlo energy by randomly drawing from node function.
540 547 ***************************************************************************/
541 548 GEnergy GModelSpectralNodes::mc(const GEnergy& emin,
542 549 const GEnergy& emax,
... ... @@ -552,41 +559,36 @@ GEnergy GModelSpectralNodes::mc(const GEnergy&amp; emin,
552 559 // Allocate energy
553 560 GEnergy energy;
554 561  
555   - // Continue only if emax > emin
556   - if (emax > emin) {
557   -
558   - // Update cache
559   - mc_update(emin, emax);
560   -
561   - // Determine in which bin we reside
562   - int inx = 0;
563   - if (m_mc_cum.size() > 1) {
564   - double u = ran.uniform();
565   - for (inx = m_mc_cum.size()-1; inx > 0; --inx) {
566   - if (m_mc_cum[inx-1] <= u)
567   - break;
568   - }
569   - }
570   -
571   - // Get random energy for specific bin
572   - if (m_mc_exp[inx] != 0.0) {
573   - double e_min = m_mc_min[inx];
574   - double e_max = m_mc_max[inx];
575   - double u = ran.uniform();
576   - double eng = (u > 0.0)
577   - ? std::exp(std::log(u * (e_max - e_min) + e_min) / m_mc_exp[inx])
578   - : 0.0;
579   - energy.MeV(eng);
580   - }
581   - else {
582   - double e_min = m_mc_min[inx];
583   - double e_max = m_mc_max[inx];
584   - double u = ran.uniform();
585   - double eng = std::exp(u * (e_max - e_min) + e_min);
586   - energy.MeV(eng);
  562 + // Update cache
  563 + mc_update(emin, emax);
  564 +
  565 + // Determine in which bin we reside
  566 + int inx = 0;
  567 + if (m_mc_cum.size() > 1) {
  568 + double u = ran.uniform();
  569 + for (inx = m_mc_cum.size()-1; inx > 0; --inx) {
  570 + if (m_mc_cum[inx-1] <= u)
  571 + break;
587 572 }
  573 + }
588 574  
589   - } // endif: emax > emin
  575 + // Get random energy for specific bin
  576 + if (m_mc_exp[inx] != 0.0) {
  577 + double e_min = m_mc_min[inx];
  578 + double e_max = m_mc_max[inx];
  579 + double u = ran.uniform();
  580 + double eng = (u > 0.0)
  581 + ? std::exp(std::log(u * (e_max - e_min) + e_min) / m_mc_exp[inx])
  582 + : 0.0;
  583 + energy.MeV(eng);
  584 + }
  585 + else {
  586 + double e_min = m_mc_min[inx];
  587 + double e_max = m_mc_max[inx];
  588 + double u = ran.uniform();
  589 + double eng = std::exp(u * (e_max - e_min) + e_min);
  590 + energy.MeV(eng);
  591 + }
590 592  
591 593 // Return energy
592 594 return energy;
... ... @@ -603,9 +605,20 @@ GEnergy GModelSpectralNodes::mc(const GEnergy&amp; emin,
603 605 * @exception GException::model_invalid_parnames
604 606 * Invalid model parameter name found in XML element.
605 607 *
606   - * Read the node function information from an XML element. The XML
607   - * element is required to have at least two nodes. Each node has two
608   - * parameters named "Energy" and "Intensity".
  608 + * Reads the spectral information from an XML element. The format of the XML
  609 + * elements is
  610 + *
  611 + * <spectrum type="NodeFunction">
  612 + * <node>
  613 + * <parameter name="Energy" scale=".." value=".." min=".." max=".." free=".."/>
  614 + * <parameter name="Intensity" scale=".." value=".." min=".." max=".." free=".."/>
  615 + * </node>
  616 + * ...
  617 + * <node>
  618 + * <parameter name="Energy" scale=".." value=".." min=".." max=".." free=".."/>
  619 + * <parameter name="Intensity" scale=".." value=".." min=".." max=".." free=".."/>
  620 + * </node>
  621 + * </spectrum>
609 622 *
610 623 * @todo Check that all energies and intensities are positive
611 624 * @todo Check that nodes are ordered
... ... @@ -713,9 +726,20 @@ void GModelSpectralNodes::read(const GXmlElement&amp; xml)
713 726 * @exception GException::model_invalid_parnames
714 727 * Invalid model parameter names found in XML element.
715 728 *
716   - * Write the node function information into an XML element. The XML element
717   - * will have node leaves, where each node leave has two parameters named
718   - * "Energy" and "Intensity".
  729 + * Writes the spectral information into an XML element. The format of the XML
  730 + * element is
  731 + *
  732 + * <spectrum type="NodeFunction">
  733 + * <node>
  734 + * <parameter name="Energy" scale=".." value=".." min=".." max=".." free=".."/>
  735 + * <parameter name="Intensity" scale=".." value=".." min=".." max=".." free=".."/>
  736 + * </node>
  737 + * ...
  738 + * <node>
  739 + * <parameter name="Energy" scale=".." value=".." min=".." max=".." free=".."/>
  740 + * <parameter name="Intensity" scale=".." value=".." min=".." max=".." free=".."/>
  741 + * </node>
  742 + * </spectrum>
719 743 ***************************************************************************/
720 744 void GModelSpectralNodes::write(GXmlElement& xml) const
721 745 {
... ... @@ -801,7 +825,320 @@ void GModelSpectralNodes::write(GXmlElement&amp; xml) const
801 825  
802 826  
803 827 /***********************************************************************//**
804   - * @brief Print powerlaw information
  828 + * @brief Append node
  829 + *
  830 + * @param[in] energy Node energy.
  831 + * @param[in] intensity Node intensity.
  832 + *
  833 + * Appends one node to the node function.
  834 + ***************************************************************************/
  835 +void GModelSpectralNodes::append(const GEnergy& energy,
  836 + const double& intensity)
  837 +{
  838 + // Allocate node parameters
  839 + GModelPar e_par;
  840 + GModelPar i_par;
  841 +
  842 + // Set energy attributes
  843 + e_par.name("Energy");
  844 + e_par.value(energy.MeV());
  845 + e_par.unit("MeV");
  846 + e_par.hasgrad(false);
  847 +
  848 + // Set intensity attributes
  849 + i_par.name("Intensity");
  850 + i_par.value(intensity);
  851 + i_par.unit("ph/cm2/s/MeV");
  852 + i_par.hasgrad(true);
  853 +
  854 + // Append to nodes
  855 + m_energies.push_back(e_par);
  856 + m_values.push_back(i_par);
  857 +
  858 + // Update parameter mapping
  859 + update_pars();
  860 +
  861 + // Set pre-computation cache
  862 + set_cache();
  863 +
  864 + // Return
  865 + return;
  866 +}
  867 +
  868 +
  869 +/***********************************************************************//**
  870 + * @brief Insert node
  871 + *
  872 + * @param[in] index Node index [0,...,nodes()-1].
  873 + * @param[in] energy Node energy.
  874 + * @param[in] intensity Node intensity.
  875 + *
  876 + * @exception GException::out_of_range
  877 + * Node index is out of range.
  878 + *
  879 + * Inserts a node into the node function before the node with the specified
  880 + * @p index.
  881 + ***************************************************************************/
  882 +void GModelSpectralNodes::insert(const int& index,
  883 + const GEnergy& energy,
  884 + const double& intensity)
  885 +{
  886 + // Raise exception if index is outside boundary
  887 + #if defined(G_RANGE_CHECK)
  888 + if (m_energies.empty()) {
  889 + if (index > 0) {
  890 + throw GException::out_of_range(G_INSERT, index, 0, nodes()-1);
  891 + }
  892 + }
  893 + else {
  894 + if (index < 0 || index >= nodes()) {
  895 + throw GException::out_of_range(G_INSERT, index, 0, nodes()-1);
  896 + }
  897 + }
  898 + #endif
  899 +
  900 + // Allocate node parameters
  901 + GModelPar e_par;
  902 + GModelPar i_par;
  903 +
  904 + // Set energy attributes
  905 + e_par.name("Energy");
  906 + e_par.value(energy.MeV());
  907 + e_par.unit("MeV");
  908 + e_par.hasgrad(false);
  909 +
  910 + // Set intensity attributes
  911 + i_par.name("Intensity");
  912 + i_par.value(intensity);
  913 + i_par.unit("ph/cm2/s/MeV");
  914 + i_par.hasgrad(true);
  915 +
  916 + // Insert node
  917 + m_energies.insert(m_energies.begin()+index, e_par);
  918 + m_values.insert(m_values.begin()+index, i_par);
  919 +
  920 + // Update parameter mapping
  921 + update_pars();
  922 +
  923 + // Set pre-computation cache
  924 + set_cache();
  925 +
  926 + // Return
  927 + return;
  928 +}
  929 +
  930 +
  931 +/***********************************************************************//**
  932 + * @brief Remove node
  933 + *
  934 + * @param[in] index Node index [0,...,nodes()-1].
  935 + *
  936 + * @exception GException::out_of_range
  937 + * Node index is out of range.
  938 + *
  939 + * Removes node of specified @p index from the node function.
  940 + ***************************************************************************/
  941 +void GModelSpectralNodes::remove(const int& index)
  942 +{
  943 + // Raise exception if index is outside boundary
  944 + #if defined(G_RANGE_CHECK)
  945 + if (index < 0 || index >= nodes()) {
  946 + throw GException::out_of_range(G_REMOVE, index, 0, nodes()-1);
  947 + }
  948 + #endif
  949 +
  950 + // Erase energy and intensity
  951 + m_energies.erase(m_energies.begin() + index);
  952 + m_values.erase(m_values.begin() + index);
  953 +
  954 + // Update parameter mapping
  955 + update_pars();
  956 +
  957 + // Set pre-computation cache
  958 + set_cache();
  959 +
  960 + // Return
  961 + return;
  962 +}
  963 +
  964 +
  965 +/***********************************************************************//**
  966 + * @brief Reserve space for nodes
  967 + *
  968 + * @param[in] num Number of reseved nodes.
  969 + *
  970 + * Reserves space for @p num nodes.
  971 + ***************************************************************************/
  972 +void GModelSpectralNodes::reserve(const int& num)
  973 +{
  974 + // Reserve space
  975 + m_energies.reserve(num);
  976 + m_values.reserve(num);
  977 +
  978 + // Return
  979 + return;
  980 +}
  981 +
  982 +
  983 +/***********************************************************************//**
  984 + * @brief Append nodes from node function
  985 + *
  986 + * @param[in] nodes Node function.
  987 + *
  988 + * Appends all nodes from a node function to current object.
  989 + ***************************************************************************/
  990 +void GModelSpectralNodes::extend(const GModelSpectralNodes& nodes)
  991 +{
  992 + // Get number of nodes in node function. Note that we extract the size
  993 + // first to avoid an endless loop that arises when a container is
  994 + // appended to itself.
  995 + int num = nodes.nodes();
  996 +
  997 + // Continue only if node function is not empty
  998 + if (num > 0) {
  999 +
  1000 + // Reserve enough space
  1001 + reserve(this->nodes() + num);
  1002 +
  1003 + // Loop over all nodes and append them to the node function
  1004 + for (int i = 0; i < num; ++i) {
  1005 + m_energies.push_back(nodes[2*i]);
  1006 + m_values.push_back(nodes[2*i+1]);
  1007 + }
  1008 +
  1009 + // Update parameter mapping
  1010 + update_pars();
  1011 +
  1012 + // Set pre-computation cache
  1013 + set_cache();
  1014 +
  1015 + } // endif: node function was not empty
  1016 +
  1017 + // Return
  1018 + return;
  1019 +}
  1020 +
  1021 +
  1022 +/***********************************************************************//**
  1023 + * @brief Return node energy
  1024 + *
  1025 + * @param[in] index Node index [0,...,nodes()-1].
  1026 + * @return Energy of node @p index.
  1027 + *
  1028 + * @exception GException::out_of_range
  1029 + * Index is out of range.
  1030 + *
  1031 + * Returns the energy of node @p index.
  1032 + ***************************************************************************/
  1033 +GEnergy GModelSpectralNodes::energy(const int& index) const
  1034 +{
  1035 + // Raise an exception if index is out of range
  1036 + #if defined(G_RANGE_CHECK)
  1037 + if (index < 0 || index >= nodes()) {
  1038 + throw GException::out_of_range(G_ENERGY_GET, index, nodes()-1);
  1039 + }
  1040 + #endif
  1041 +
  1042 + // Retrieve energy
  1043 + GEnergy energy;
  1044 + energy.MeV(m_energies[index].value());
  1045 +
  1046 + // Return energy
  1047 + return energy;
  1048 +}
  1049 +
  1050 +
  1051 +/***********************************************************************//**
  1052 + * @brief Set node energy
  1053 + *
  1054 + * @param[in] index Node index [0,...,nodes()-1].
  1055 + * @param[in] energy Node energy.
  1056 + *
  1057 + * @exception GException::out_of_range
  1058 + * Index is out of range.
  1059 + *
  1060 + * Returns the energy of node @p index.
  1061 + ***************************************************************************/
  1062 +void GModelSpectralNodes::energy(const int& index, const GEnergy& energy)
  1063 +{
  1064 + // Raise an exception if index is out of range
  1065 + #if defined(G_RANGE_CHECK)
  1066 + if (index < 0 || index >= nodes()) {
  1067 + throw GException::out_of_range(G_ENERGY_SET, index, nodes()-1);
  1068 + }
  1069 + #endif
  1070 +
  1071 + // Set energy
  1072 + m_energies[index].value(energy.MeV());
  1073 +
  1074 + // Set pre-computation cache
  1075 + set_cache();
  1076 +
  1077 + // Return
  1078 + return;
  1079 +}
  1080 +
  1081 +
  1082 +/***********************************************************************//**
  1083 + * @brief Return node intensity
  1084 + *
  1085 + * @param[in] index Node index [0,...,nodes()-1].
  1086 + * @return Energy of node @p index.
  1087 + *
  1088 + * @exception GException::out_of_range
  1089 + * Index is out of range.
  1090 + *
  1091 + * Returns the energy of node @p index.
  1092 + ***************************************************************************/
  1093 +double GModelSpectralNodes::intensity(const int& index) const
  1094 +{
  1095 + // Raise an exception if index is out of range
  1096 + #if defined(G_RANGE_CHECK)
  1097 + if (index < 0 || index >= nodes()) {
  1098 + throw GException::out_of_range(G_INTENSITY_GET, index, nodes()-1);
  1099 + }
  1100 + #endif
  1101 +
  1102 + // Return intensity
  1103 + return (m_values[index].value());
  1104 +}
  1105 +
  1106 +
  1107 +/***********************************************************************//**
  1108 + * @brief Set node intensity
  1109 + *
  1110 + * @param[in] index Node index [0,...,nodes()-1].
  1111 + * @param[in] energy Node energy.
  1112 + *
  1113 + * @exception GException::out_of_range
  1114 + * Index is out of range.
  1115 + *
  1116 + * Returns the energy of node @p index.
  1117 + ***************************************************************************/
  1118 +void GModelSpectralNodes::intensity(const int& index, const double& intensity)
  1119 +{
  1120 + // Raise an exception if index is out of range
  1121 + #if defined(G_RANGE_CHECK)
  1122 + if (index < 0 || index >= nodes()) {
  1123 + throw GException::out_of_range(G_INTENSITY_SET, index, nodes()-1);
  1124 + }
  1125 + #endif
  1126 +
  1127 + // Set intensity
  1128 + m_values[index].value(intensity);
  1129 +
  1130 + // Set pre-computation cache
  1131 + set_cache();
  1132 +
  1133 + // Return
  1134 + return;
  1135 +}
  1136 +
  1137 +
  1138 +/***********************************************************************//**
  1139 + * @brief Print node function information
  1140 + *
  1141 + * @return String containing node function information.
805 1142 ***************************************************************************/
806 1143 std::string GModelSpectralNodes::print(void) const
807 1144 {
... ... @@ -810,6 +1147,8 @@ std::string GModelSpectralNodes::print(void) const
810 1147  
811 1148 // Append header
812 1149 result.append("=== GModelSpectralNodes ===");
  1150 +
  1151 + // Append information
813 1152 result.append("\n"+parformat("Number of nodes")+str(m_energies.size()));
814 1153 result.append("\n"+parformat("Number of parameters")+str(size()));
815 1154 for (int i = 0; i < size(); ++i) {
... ... @@ -1310,10 +1649,3 @@ void GModelSpectralNodes::mc_update(const GEnergy&amp; emin, const GEnergy&amp; emax) co
1310 1649 // Return
1311 1650 return;
1312 1651 }
1313   -
1314   -
1315   -/*==========================================================================
1316   - = =
1317   - = Friends =
1318   - = =
1319   - ==========================================================================*/
... ...