GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpectralNodes.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpectralNodes.cpp - Spectral nodes model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2012-2018 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GModelSpectralNodes.cpp
23  * @brief Spectral nodes model class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cmath>
32 #include "GException.hpp"
33 #include "GTools.hpp"
34 #include "GRan.hpp"
35 #include "GEnergies.hpp"
36 #include "GModelSpectralNodes.hpp"
38 
39 /* __ Constants __________________________________________________________ */
40 
41 /* __ Globals ____________________________________________________________ */
43 const GModelSpectralRegistry g_spectral_nodes_registry(&g_spectral_nodes_seed);
44 
45 /* __ Method name definitions ____________________________________________ */
46 #define G_FLUX "GModelSpectralNodes::flux(GEnergy&, GEnergy&)"
47 #define G_EFLUX "GModelSpectralNodes::eflux(GEnergy&, GEnergy&)"
48 #define G_MC "GModelSpectralNodes::mc(GEnergy&, GEnergy&, GTime&, GRan&)"
49 #define G_READ "GModelSpectralNodes::read(GXmlElement&)"
50 #define G_WRITE "GModelSpectralNodes::write(GXmlElement&)"
51 #define G_APPEND "GModelSpectralNodes::append(GEnergy&, double&)"
52 #define G_INSERT "GModelSpectralNodes::insert(int&, GEnergy&, double&)"
53 #define G_REMOVE "GModelSpectralNodes::remove(int&)"
54 #define G_ENERGY_GET "GModelSpectralNodes::energy(int&)"
55 #define G_ENERGY_SET "GModelSpectralNodes::energy(int&, GEnergy&)"
56 #define G_INTENSITY_GET "GModelSpectralNodes::intensity(int&)"
57 #define G_INTENSITY_SET "GModelSpectralNodes::intensity(int&, double&)"
58 
59 /* __ Macros _____________________________________________________________ */
60 
61 /* __ Coding definitions _________________________________________________ */
62 
63 /* __ Debug definitions __________________________________________________ */
64 
65 
66 /*==========================================================================
67  = =
68  = Constructors/destructors =
69  = =
70  ==========================================================================*/
71 
72 /***********************************************************************//**
73  * @brief Void constructor
74  *
75  * Constructs an empty spectral node model.
76  ***************************************************************************/
78 {
79  // Initialise members
80  init_members();
81 
82  // Return
83  return;
84 }
85 
86 
87 /***********************************************************************//**
88  * @brief Spectral model constructor
89  *
90  * @param[in] model Spectral model.
91  * @param[in] energies Node energies.
92  *
93  * Constructs a spectral node model from any spectral model.
94  ***************************************************************************/
96  const GEnergies& energies) :
98 {
99  // Initialise members
100  init_members();
101 
102  // Append nodes for all energies
103  for (int i = 0; i < energies.size(); ++i) {
104  append(energies[i], model.eval(energies[i]));
105  }
106 
107  // Return
108  return;
109 }
110 
111 
112 /***********************************************************************//**
113  * @brief XML constructor
114  *
115  * @param[in] xml XML element.
116  *
117  * Construct spectral nodes model by extracting information from an XML
118  * element. See the read() method for more information about the expected
119  * structure of the XML element.
120  ***************************************************************************/
123 {
124  // Initialise members
125  init_members();
126 
127  // Read information from XML element
128  read(xml);
129 
130  // Return
131  return;
132 }
133 
134 
135 /***********************************************************************//**
136  * @brief Copy constructor
137  *
138  * @param[in] model Spectral nodes model.
139  ***************************************************************************/
141  GModelSpectral(model)
142 {
143  // Initialise members
144  init_members();
145 
146  // Copy members
147  copy_members(model);
148 
149  // Return
150  return;
151 }
152 
153 
154 /***********************************************************************//**
155  * @brief Destructor
156  ***************************************************************************/
158 {
159  // Free members
160  free_members();
161 
162  // Return
163  return;
164 }
165 
166 
167 /*==========================================================================
168  = =
169  = Operators =
170  = =
171  ==========================================================================*/
172 
173 /***********************************************************************//**
174  * @brief Assignment operator
175  *
176  * @param[in] model Spectral nodes model.
177  * @return Spectral nodes model.
178  ***************************************************************************/
180 {
181  // Execute only if object is not identical
182  if (this != &model) {
183 
184  // Copy base class members
185  this->GModelSpectral::operator=(model);
186 
187  // Free members
188  free_members();
189 
190  // Initialise members
191  init_members();
192 
193  // Copy members
194  copy_members(model);
195 
196  } // endif: object was not identical
197 
198  // Return
199  return *this;
200 }
201 
202 
203 /*==========================================================================
204  = =
205  = Public methods =
206  = =
207  ==========================================================================*/
208 
209 /***********************************************************************//**
210  * @brief Clear spectral nodes model
211  ***************************************************************************/
213 {
214  // Free class members (base and derived classes, derived class first)
215  free_members();
217 
218  // Initialise members
220  init_members();
221 
222  // Return
223  return;
224 }
225 
226 
227 /***********************************************************************//**
228  * @brief Clone spectral nodes model
229 ***************************************************************************/
231 {
232  // Clone spectral nodes model
233  return new GModelSpectralNodes(*this);
234 }
235 
236 
237 /***********************************************************************//**
238  * @brief Evaluate function
239  *
240  * @param[in] srcEng True photon energy.
241  * @param[in] srcTime True photon arrival time.
242  * @param[in] gradients Compute gradients?
243  * @return Model value (ph/cm2/s/MeV).
244  *
245  * Computes
246  *
247  * \f[
248  * S_{\rm E}(E | t) =
249  * 10^{(\log {\tt m\_values[i]}) w_{i} +
250  * (\log {\tt m\_values[i+1]}) w_{i+1}}
251  * \f]
252  *
253  * where
254  * - \f${\tt m\_values[i]}\f$ is the intensity of node \f$i\f$,
255  * - \f${\tt m\_values[i+1]}\f$ is the intensity of node \f$i+1\f$,
256  * - \f$w_{i}\f$ is the weighting of node \f$i\f$,
257  * - \f$w_{i+1}\f$ is the weighting of node \f$i+1\f$, and
258  * - \f${\tt m\_energies[i]} <= E <= {\tt m\_energies[i+1]}\f$.
259  *
260  * The weightings \f$w_{i}\f$ and \f$w_{i+1}\f$ are computed by linear
261  * interpolation (in the log-log plane) between the nodes
262  * \f$(\log {\tt m\_energies[i]}, \log{\tt m\_values[i]})\f$
263  * and
264  * \f$(\log {\tt m\_energies[i+1]}, \log{\tt m\_values[i+1]})\f$
265  * to the requested energy \f$\log E\f$.
266  *
267  * If the @p gradients flag is true the method will also compute the partial
268  * derivatives of the model with respect to the parameters using
269  *
270  * \f[
271  * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_values[i]}} =
272  * \frac{S_{\rm E}(E | t) \, w_i}{{\tt m\_values[i]}}
273  * \f]
274  ***************************************************************************/
275 double GModelSpectralNodes::eval(const GEnergy& srcEng,
276  const GTime& srcTime,
277  const bool& gradients) const
278 {
279  // Update evaluation cache
281 
282  // Set interpolation value
284 
285  // Get indices and weights for interpolation
286  int inx_left = m_log_energies.inx_left();
287  int inx_right = m_log_energies.inx_right();
288  double wgt_left = m_log_energies.wgt_left();
289  double wgt_right = m_log_energies.wgt_right();
290 
291  // Interpolate function
292  double exponent = m_log_values[inx_left] * wgt_left +
293  m_log_values[inx_right] * wgt_right;
294 
295  // Compute linear value
296  double value = std::pow(10.0, exponent);
297 
298  // Optionally compute partial derivatives
299  if (gradients) {
300 
301  // Initialise gradients
302  for (int i = 0; i < m_values.size(); ++i) {
303  m_values[i].factor_gradient(0.0);
304  }
305 
306  // Gradient for left node
307  if (m_values[inx_left].is_free()) {
308  double grad = value * wgt_left / m_values[inx_left].factor_value();
309  m_values[inx_left].factor_gradient(grad);
310  }
311 
312  // Gradient for right node
313  if (m_values[inx_right].is_free()) {
314  double grad = value * wgt_right / m_values[inx_right].factor_value();
315  m_values[inx_right].factor_gradient(grad);
316  }
317 
318  } // endif: computed partial derivatives
319 
320  // Compile option: Check for NaN/Inf
321  #if defined(G_NAN_CHECK)
322  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
323  std::cout << "*** ERROR: GModelSpectralNodes::eval";
324  std::cout << "(srcEng=" << srcEng;
325  std::cout << ", srcTime=" << srcTime << "):";
326  std::cout << " NaN/Inf encountered";
327  std::cout << " (value=" << value;
328  std::cout << ", exponent=" << exponent;
329  std::cout << ")" << std::endl;
330  }
331  #endif
332 
333  // Return
334  return value;
335 }
336 
337 
338 /***********************************************************************//**
339  * @brief Returns model photon flux between [emin, emax] (units: ph/cm2/s)
340  *
341  * @param[in] emin Minimum photon energy.
342  * @param[in] emax Maximum photon energy.
343  * @return Photon flux (ph/cm2/s).
344  *
345  * Computes
346  *
347  * \f[
348  * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) dE
349  * \f]
350  *
351  * where
352  * - [@p emin, @p emax] is an energy interval, and
353  * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
354  ***************************************************************************/
355 double GModelSpectralNodes::flux(const GEnergy& emin, const GEnergy& emax) const
356 {
357  // Initialise flux
358  double flux = 0.0;
359 
360  // Compute only if integration range is valid
361  if (emin < emax) {
362 
363  // Update flux cache
365 
366  // Get energy range in MeV
367  double e_min = emin.MeV();
368  double e_max = emax.MeV();
369 
370  // Determine left node index for minimum energy
371  m_lin_energies.set_value(e_min);
372  int inx_emin = m_lin_energies.inx_left();
373 
374  // Determine left node index for maximum energy
375  m_lin_energies.set_value(e_max);
376  int inx_emax = m_lin_energies.inx_left();
377 
378  // If both energies are within the same nodes then simply
379  // integrate over the energy interval using the appropriate power
380  // law parameters
381  if (inx_emin == inx_emax) {
382  flux = m_prefactor[inx_emin] *
384  e_max,
385  m_epivot[inx_emin],
386  m_gamma[inx_emin]);
387  }
388 
389  // ... otherwise integrate over the nodes where emin and emax
390  // resides and all the remaining nodes
391  else {
392 
393  // If we are requested to extrapolate beyond first node,
394  // the use the first nodes lower energy and upper energy
395  // boundary for the initial integration.
396  int i_start = (e_min < m_lin_energies[0]) ? inx_emin : inx_emin+1;
397 
398  // Integrate from emin to the node boundary
399  flux = m_prefactor[inx_emin] *
401  m_lin_energies[i_start],
402  m_epivot[inx_emin],
403  m_gamma[inx_emin]);
404 
405  // Integrate over all nodes between
406  for (int i = i_start; i < inx_emax; ++i) {
407  flux += m_flux[i];
408  }
409 
410  // Integrate from node boundary to emax
411  flux += m_prefactor[inx_emax] *
413  e_max,
414  m_epivot[inx_emax],
415  m_gamma[inx_emax]);
416  }
417 
418  } // endif: integration range was valid
419 
420  // Return
421  return flux;
422 }
423 
424 
425 /***********************************************************************//**
426  * @brief Returns model energy flux between [emin, emax] (units: erg/cm2/s)
427  *
428  * @param[in] emin Minimum photon energy.
429  * @param[in] emax Maximum photon energy.
430  * @return Energy flux (erg/cm2/s).
431  *
432  * Computes
433  *
434  * \f[
435  * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) E \, dE
436  * \f]
437  *
438  * where
439  * - [@p emin, @p emax] is an energy interval, and
440  * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
441  ***************************************************************************/
442 double GModelSpectralNodes::eflux(const GEnergy& emin, const GEnergy& emax) const
443 {
444  // Initialise flux
445  double eflux = 0.0;
446 
447  // Compute only if integration range is valid
448  if (emin < emax) {
449 
450  // Update flux cache
452 
453  // Get energy range in MeV
454  double e_min = emin.MeV();
455  double e_max = emax.MeV();
456 
457  // Determine left node index for minimum energy
458  m_lin_energies.set_value(e_min);
459  int inx_emin = m_lin_energies.inx_left();
460 
461  // Determine left node index for maximum energy
462  m_lin_energies.set_value(e_max);
463  int inx_emax = m_lin_energies.inx_left();
464 
465  // If both energies are within the same nodes then simply
466  // integrate over the energy interval using the appropriate power
467  // law parameters
468  if (inx_emin == inx_emax) {
469  eflux = m_prefactor[inx_emin] *
471  e_max,
472  m_epivot[inx_emin],
473  m_gamma[inx_emin]) *
475  }
476 
477  // ... otherwise integrate over the nodes where emin and emax
478  // resides and all the remaining nodes
479  else {
480 
481  // If we are requested to extrapolate beyond first node,
482  // the use the first nodes lower energy and upper energy
483  // boundary for the initial integration.
484  int i_start = (e_min < m_lin_energies[0]) ? inx_emin : inx_emin+1;
485 
486  // Integrate from emin to the node boundary
487  eflux = m_prefactor[inx_emin] *
489  m_lin_energies[i_start],
490  m_epivot[inx_emin],
491  m_gamma[inx_emin]) *
493 
494  // Integrate over all nodes between
495  for (int i = i_start; i < inx_emax; ++i) {
496  eflux += m_eflux[i];
497  }
498 
499  // Integrate from node boundary to emax
500  eflux += m_prefactor[inx_emax] *
502  e_max,
503  m_epivot[inx_emax],
504  m_gamma[inx_emax]) *
506 
507  } // endelse: emin and emax not between same nodes
508 
509  } // endif: integration range was valid
510 
511  // Return flux
512  return eflux;
513 }
514 
515 
516 /***********************************************************************//**
517  * @brief Returns MC energy between [emin, emax]
518  *
519  * @param[in] emin Minimum photon energy.
520  * @param[in] emax Maximum photon energy.
521  * @param[in] time True photon arrival time.
522  * @param[in,out] ran Random number generator.
523  * @return Energy.
524  *
525  * @exception GException::erange_invalid
526  * Energy range is invalid (emin < emax required).
527  * @exception GException::invalid_return_value
528  * No valid Monte Carlo cache
529  *
530  * Returns Monte Carlo energy by randomly drawing from node function.
531  ***************************************************************************/
533  const GEnergy& emax,
534  const GTime& time,
535  GRan& ran) const
536 {
537  // Throw an exception if energy range is invalid
538  if (emin >= emax) {
539  throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(),
540  "Minimum energy < maximum energy required.");
541  }
542 
543  // Allocate energy
544  GEnergy energy;
545 
546  // Update cache
547  mc_update(emin, emax);
548 
549  // Determine in which bin we reside. If this fails then thrown an
550  // exception
551  int inx = 0;
552  if (m_mc_cum.size() > 1) {
553  double u = ran.uniform();
554  for (inx = m_mc_cum.size()-1; inx > 0; --inx) {
555  if (m_mc_cum[inx-1] <= u) {
556  break;
557  }
558  }
559  }
560  else if (m_mc_cum.size() == 0) {
561  std::string msg = "No valid nodes found for energy interval ["+
562  emin.print()+","+emax.print()+"]. Either restrict "
563  "the energy range to the one covered by the "
564  "spectral nodes or extend the spectral nodes "
565  "in energy.";
567  }
568 
569  // Get random energy for specific bin
570  if (m_mc_exp[inx] != 0.0) {
571  double e_min = m_mc_min[inx];
572  double e_max = m_mc_max[inx];
573  double u = ran.uniform();
574  double eng = (u > 0.0)
575  ? std::exp(std::log(u * (e_max - e_min) + e_min) / m_mc_exp[inx])
576  : 0.0;
577  energy.MeV(eng);
578  }
579  else {
580  double e_min = m_mc_min[inx];
581  double e_max = m_mc_max[inx];
582  double u = ran.uniform();
583  double eng = std::exp(u * (e_max - e_min) + e_min);
584  energy.MeV(eng);
585  }
586 
587  // Return energy
588  return energy;
589 }
590 
591 
592 /***********************************************************************//**
593  * @brief Read model from XML element
594  *
595  * @param[in] xml XML element.
596  *
597  * @exception GException::model_invalid_parnum
598  * Invalid number of model parameters found in XML element.
599  * @exception GException::model_invalid_parnames
600  * Invalid model parameter name found in XML element.
601  * @exception GException::invalid_value
602  * Energy or intensity are not positive.
603  *
604  * Reads the spectral information from an XML element. The format of the XML
605  * elements is
606  *
607  * <spectrum type="NodeFunction">
608  * <node>
609  * <parameter name="Energy" ../>
610  * <parameter name="Intensity" ../>
611  * </node>
612  * ...
613  * <node>
614  * <parameter name="Energy" ../>
615  * <parameter name="Intensity" ../>
616  * </node>
617  * </spectrum>
618  *
619  * @todo Check that nodes are ordered
620  * @todo Check that energy boundaries are not overlapping
621  ***************************************************************************/
623 {
624  // Free space for nodes
625  m_energies.clear();
626  m_values.clear();
627 
628  // Get number of nodes from XML file
629  int nodes = xml.elements("node");
630 
631  // Throw an error if there not at least two nodes
632  if (nodes < 1) {
633  std::string message = "Node function model requires at least two"
634  " nodes.";
635  throw GException::model_invalid_parnum(G_READ, xml, message);
636  }
637 
638  // Loop over all nodes
639  for (int i = 0; i < nodes; ++i) {
640 
641  // Allocate node parameters
644 
645  // Get node
646  const GXmlElement* node = xml.element("node", i);
647 
648  // Get parameters
649  const GXmlElement* epar = gammalib::xml_get_par(G_READ, *node, "Energy");
650  const GXmlElement* ipar = gammalib::xml_get_par(G_READ, *node, "Intensity");
651 
652  // Read parameters
653  energy.read(*epar);
654  intensity.read(*ipar);
655 
656  // Throw an exception if either energy or intensity is not positive
657  if (energy.value() <= 0.0) {
658  std::string msg = "Non positive energy "+
659  gammalib::str(energy.value())+" MeV encountered "
660  "in model definition XML file. Please specify "
661  "only nodes with positive energy.";
662  throw GException::invalid_value(G_READ, msg);
663  }
664  if (intensity.value() <= 0.0) {
665  std::string msg = "Non positive intensity "+
666  gammalib::str(intensity.value())+" ph/cm2/s/MeV "
667  "encountered in model definition XML file. "
668  "Please specify only nodes with positive "
669  "intensity.";
670  throw GException::invalid_value(G_READ, msg);
671  }
672 
673  // Set parameter names
674  std::string energy_name = "Energy"+gammalib::str(i);
675  std::string intensity_name = "Intensity"+gammalib::str(i);
676 
677  // Set energy attributes
678  energy.name(energy_name);
679  energy.unit("MeV");
680  energy.has_grad(false);
681 
682  // Set intensity attributes
683  intensity.name(intensity_name);
684  intensity.unit("ph/cm2/s/MeV");
685  intensity.has_grad(true);
686 
687  // Push node parameters on list
688  m_energies.push_back(energy);
689  m_values.push_back(intensity);
690 
691  } // endfor: looped over nodes
692 
693  // Update parameter mapping
694  update_pars();
695 
696  // Set pre-computation cache
697  set_cache();
698 
699  // Return
700  return;
701 }
702 
703 
704 /***********************************************************************//**
705  * @brief Write model into XML element
706  *
707  * @param[in] xml XML element into which model information is written.
708  *
709  * @exception GException::model_invalid_spectral
710  * Existing XML element is not of required type
711  * @exception GException::model_invalid_parnum
712  * Invalid number of model parameters or nodes found in XML element.
713  *
714  * Writes the spectral information into an XML element. The format of the XML
715  * element is
716  *
717  * <spectrum type="NodeFunction">
718  * <node>
719  * <parameter name="Energy" ../>
720  * <parameter name="Intensity" ../>
721  * </node>
722  * ...
723  * <node>
724  * <parameter name="Energy" ../>
725  * <parameter name="Intensity" ../>
726  * </node>
727  * </spectrum>
728  ***************************************************************************/
730 {
731  // Determine number of nodes
732  int nodes = m_energies.size();
733 
734  // Set model type
735  if (xml.attribute("type") == "") {
736  xml.attribute("type", type());
737  }
738 
739  // Verify model type
740  if (xml.attribute("type") != type()) {
742  "Spectral model is not of type \""+type()+"\".");
743  }
744 
745  // If XML element has 0 nodes then append nodes
746  if (xml.elements() == 0) {
747  for (int i = 0; i < nodes; ++i) {
748  xml.append(GXmlElement("node"));
749  }
750  }
751 
752  // Verify that XML element has the required number of nodes
753  if (xml.elements() != nodes || xml.elements("node") != nodes) {
754  std::string message = "Spectral model requires exactly " +
755  gammalib::str(nodes) + " nodes.";
756  throw GException::model_invalid_parnum(G_WRITE, xml, message);
757  }
758 
759  // Loop over all nodes
760  for (int i = 0; i < nodes; ++i) {
761 
762  // Get node
763  GXmlElement* node = xml.element("node", i);
764 
765  // Get copy of parameters so that we can change their names
768 
769  // Set names since we appended for internal usage the indices to the
770  // parameter names, and we want to get rid of them for writing the
771  // model into the XML element
772  energy.name("Energy");
773  intensity.name("Intensity");
774 
775  // Get XML parameters
776  GXmlElement* eng = gammalib::xml_need_par(G_WRITE, *node, energy.name());
777  GXmlElement* val = gammalib::xml_need_par(G_WRITE, *node, intensity.name());
778 
779  // Write parameters
780  energy.write(*eng);
781  intensity.write(*val);
782 
783  } // endfor: looped over nodes
784 
785  // Return
786  return;
787 }
788 
789 
790 /***********************************************************************//**
791  * @brief Append node
792  *
793  * @param[in] energy Node energy.
794  * @param[in] intensity Node intensity.
795  *
796  * @exception GException::invalid_argument
797  * Non-positive energy or intensity specified.
798  *
799  * Appends one node to the node function. By default the energy of the node
800  * is fixed while the intensity of the node is free.
801  ***************************************************************************/
803  const double& intensity)
804 {
805  // Throw an exception if either energy or intensity is not positive
806  if (energy.MeV() <= 0.0) {
807  std::string msg = "Non-positive energy "+energy.print()+" not allowed. "
808  "Please specify only positive energies.";
810  }
811  if (intensity <= 0.0) {
812  std::string msg = "Non-positive intensity "+
813  gammalib::str(intensity)+" ph/cm2/s/MeV "
814  "not allowed. Please specify only positive "
815  "intensities.";
817  }
818 
819  // Allocate node parameters
820  GModelPar e_par;
821  GModelPar i_par;
822 
823  // Set energy attributes
824  e_par.name("Energy");
825  e_par.value(energy.MeV());
826  e_par.unit("MeV");
827  e_par.has_grad(false);
828  e_par.fix();
829 
830  // Set intensity attributes
831  i_par.name("Intensity");
832  i_par.value(intensity);
833  i_par.unit("ph/cm2/s/MeV");
834  i_par.has_grad(true);
835  i_par.free();
836 
837  // Append to nodes
838  m_energies.push_back(e_par);
839  m_values.push_back(i_par);
840 
841  // Update parameter mapping
842  update_pars();
843 
844  // Set pre-computation cache
845  set_cache();
846 
847  // Return
848  return;
849 }
850 
851 
852 /***********************************************************************//**
853  * @brief Insert node
854  *
855  * @param[in] index Node index [0,...,nodes()-1].
856  * @param[in] energy Node energy.
857  * @param[in] intensity Node intensity.
858  *
859  * @exception GException::out_of_range
860  * Node index is out of range.
861  * @exception GException::invalid_argument
862  * Non-positive energy or intensity specified.
863  *
864  * Inserts a node into the node function before the node with the specified
865  * @p index. By default the energy of the node is fixed while the intensity
866  * of the node is free.
867  ***************************************************************************/
868 void GModelSpectralNodes::insert(const int& index,
869  const GEnergy& energy,
870  const double& intensity)
871 {
872  // Raise exception if index is outside boundary
873  #if defined(G_RANGE_CHECK)
874  if (m_energies.empty()) {
875  if (index > 0) {
876  throw GException::out_of_range(G_INSERT, index, 0, nodes()-1);
877  }
878  }
879  else {
880  if (index < 0 || index >= nodes()) {
881  throw GException::out_of_range(G_INSERT, index, 0, nodes()-1);
882  }
883  }
884  #endif
885 
886  // Throw an exception if either energy or intensity is not positive
887  if (energy.MeV() <= 0.0) {
888  std::string msg = "Non-positive energy "+energy.print()+" not allowed. "
889  "Please specify only positive energies.";
891  }
892  if (intensity <= 0.0) {
893  std::string msg = "Non-positive intensity "+
894  gammalib::str(intensity)+" ph/cm2/s/MeV "
895  "not allowed. Please specify only positive "
896  "intensities.";
898  }
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.has_grad(false);
909  e_par.fix();
910 
911  // Set intensity attributes
912  i_par.name("Intensity");
913  i_par.value(intensity);
914  i_par.unit("ph/cm2/s/MeV");
915  i_par.has_grad(true);
916  i_par.free();
917 
918  // Insert node
919  m_energies.insert(m_energies.begin()+index, e_par);
920  m_values.insert(m_values.begin()+index, i_par);
921 
922  // Update parameter mapping
923  update_pars();
924 
925  // Set pre-computation cache
926  set_cache();
927 
928  // Return
929  return;
930 }
931 
932 
933 /***********************************************************************//**
934  * @brief Remove node
935  *
936  * @param[in] index Node index [0,...,nodes()-1].
937  *
938  * @exception GException::out_of_range
939  * Node index is out of range.
940  *
941  * Removes node of specified @p index from the node function.
942  ***************************************************************************/
943 void GModelSpectralNodes::remove(const int& index)
944 {
945  // Raise exception if index is outside boundary
946  #if defined(G_RANGE_CHECK)
947  if (index < 0 || index >= nodes()) {
948  throw GException::out_of_range(G_REMOVE, index, 0, nodes()-1);
949  }
950  #endif
951 
952  // Erase energy and intensity
953  m_energies.erase(m_energies.begin() + index);
954  m_values.erase(m_values.begin() + index);
955 
956  // Update parameter mapping
957  update_pars();
958 
959  // Set pre-computation cache
960  set_cache();
961 
962  // Return
963  return;
964 }
965 
966 
967 /***********************************************************************//**
968  * @brief Reserve space for nodes
969  *
970  * @param[in] num Number of reseved nodes.
971  *
972  * Reserves space for @p num nodes.
973  ***************************************************************************/
974 void GModelSpectralNodes::reserve(const int& num)
975 {
976  // Reserve space
977  m_energies.reserve(num);
978  m_values.reserve(num);
979 
980  // Return
981  return;
982 }
983 
984 
985 /***********************************************************************//**
986  * @brief Append nodes from node function
987  *
988  * @param[in] nodes Node function.
989  *
990  * Appends all nodes from a node function to current object.
991  ***************************************************************************/
993 {
994  // Get number of nodes in node function. Note that we extract the size
995  // first to avoid an endless loop that arises when a container is
996  // appended to itself.
997  int num = nodes.nodes();
998 
999  // Continue only if node function is not empty
1000  if (num > 0) {
1001 
1002  // Reserve enough space
1003  reserve(this->nodes() + num);
1004 
1005  // Loop over all nodes and append them to the node function
1006  for (int i = 0; i < num; ++i) {
1007  m_energies.push_back(nodes.m_energies[i]);
1008  m_values.push_back(nodes.m_values[i]);
1009  }
1010 
1011  // Update parameter mapping
1012  update_pars();
1013 
1014  // Set pre-computation cache
1015  set_cache();
1016 
1017  } // endif: node function was not empty
1018 
1019  // Return
1020  return;
1021 }
1022 
1023 
1024 /***********************************************************************//**
1025  * @brief Return node energy
1026  *
1027  * @param[in] index Node index [0,...,nodes()-1].
1028  * @return Energy of node @p index.
1029  *
1030  * @exception GException::out_of_range
1031  * Index is out of range.
1032  *
1033  * Returns the energy of node @p index.
1034  ***************************************************************************/
1035 GEnergy GModelSpectralNodes::energy(const int& index) const
1036 {
1037  // Raise an exception if index is out of range
1038  #if defined(G_RANGE_CHECK)
1039  if (index < 0 || index >= nodes()) {
1040  throw GException::out_of_range(G_ENERGY_GET, index, nodes()-1);
1041  }
1042  #endif
1043 
1044  // Retrieve energy
1045  GEnergy energy;
1046  energy.MeV(m_energies[index].value());
1047 
1048  // Return energy
1049  return energy;
1050 }
1051 
1052 
1053 /***********************************************************************//**
1054  * @brief Set node energy
1055  *
1056  * @param[in] index Node index [0,...,nodes()-1].
1057  * @param[in] energy Node energy.
1058  *
1059  * @exception GException::out_of_range
1060  * Index is out of range.
1061  * @exception GException::invalid_argument
1062  * Non-positive energy specified.
1063  *
1064  * Sets the energy of node @p index.
1065  ***************************************************************************/
1066 void GModelSpectralNodes::energy(const int& index, const GEnergy& energy)
1067 {
1068  // Raise an exception if index is out of range
1069  #if defined(G_RANGE_CHECK)
1070  if (index < 0 || index >= nodes()) {
1071  throw GException::out_of_range(G_ENERGY_SET, index, nodes()-1);
1072  }
1073  #endif
1074 
1075  // Throw an exception if energy is not positive
1076  if (energy.MeV() <= 0.0) {
1077  std::string msg = "Non-positive energy "+energy.print()+" not allowed. "
1078  "Please specify only positive energies.";
1080  }
1081 
1082  // Set energy
1083  m_energies[index].value(energy.MeV());
1084 
1085  // Set pre-computation cache
1086  set_cache();
1087 
1088  // Return
1089  return;
1090 }
1091 
1092 
1093 /***********************************************************************//**
1094  * @brief Return node intensity
1095  *
1096  * @param[in] index Node index [0,...,nodes()-1].
1097  * @return Intensity of node @p index.
1098  *
1099  * @exception GException::out_of_range
1100  * Index is out of range.
1101  *
1102  * Returns the intensity of node @p index.
1103  ***************************************************************************/
1104 double GModelSpectralNodes::intensity(const int& index) const
1105 {
1106  // Raise an exception if index is out of range
1107  #if defined(G_RANGE_CHECK)
1108  if (index < 0 || index >= nodes()) {
1109  throw GException::out_of_range(G_INTENSITY_GET, index, nodes()-1);
1110  }
1111  #endif
1112 
1113  // Return intensity
1114  return (m_values[index].value());
1115 }
1116 
1117 
1118 /***********************************************************************//**
1119  * @brief Set node intensity
1120  *
1121  * @param[in] index Node index [0,...,nodes()-1].
1122  * @param[in] intensity Node Intensity.
1123  *
1124  * @exception GException::out_of_range
1125  * Index is out of range.
1126  * @exception GException::invalid_argument
1127  * Non-positive intensity specified.
1128  *
1129  * Set the intensity of node @p index.
1130  ***************************************************************************/
1131 void GModelSpectralNodes::intensity(const int& index, const double& intensity)
1132 {
1133  // Raise an exception if index is out of range
1134  #if defined(G_RANGE_CHECK)
1135  if (index < 0 || index >= nodes()) {
1136  throw GException::out_of_range(G_INTENSITY_SET, index, nodes()-1);
1137  }
1138  #endif
1139 
1140  // Throw an exception if either energy or intensity is not positive
1141  if (intensity <= 0.0) {
1142  std::string msg = "Non-positive intensity "+
1143  gammalib::str(intensity)+" ph/cm2/s/MeV "
1144  "not allowed. Please specify only positive "
1145  "intensities.";
1147  }
1148 
1149  // Set intensity
1150  m_values[index].value(intensity);
1151 
1152  // Set pre-computation cache
1153  set_cache();
1154 
1155  // Return
1156  return;
1157 }
1158 
1159 
1160 /***********************************************************************//**
1161  * @brief Print node function information
1162  *
1163  * @param[in] chatter Chattiness.
1164  * @return String containing node function information.
1165  ***************************************************************************/
1166 std::string GModelSpectralNodes::print(const GChatter& chatter) const
1167 {
1168  // Initialise result string
1169  std::string result;
1170 
1171  // Continue only if chatter is not silent
1172  if (chatter != SILENT) {
1173 
1174  // Append header
1175  result.append("=== GModelSpectralNodes ===");
1176 
1177  // Append information
1178  result.append("\n"+gammalib::parformat("Number of nodes"));
1179  result.append(gammalib::str(m_energies.size()));
1180  result.append("\n"+gammalib::parformat("Number of parameters"));
1181  result.append(gammalib::str(size()));
1182  for (int i = 0; i < size(); ++i) {
1183  result.append("\n"+m_pars[i]->print(chatter));
1184  }
1185 
1186  // VERBOSE: Append information about power laws between node
1187  if (chatter == VERBOSE) {
1188  for (int i = 0; i < m_prefactor.size(); ++i) {
1189  result.append("\n"+gammalib::parformat("Interval "+gammalib::str(i+1)));
1190  result.append("Epivot="+gammalib::str(m_epivot[i]));
1191  result.append(" Prefactor="+gammalib::str(m_prefactor[i]));
1192  result.append(" Gamma="+gammalib::str(m_gamma[i]));
1193  result.append(" Flux="+gammalib::str(m_flux[i]));
1194  result.append(" EFlux="+gammalib::str(m_eflux[i]));
1195  }
1196  }
1197 
1198  } // endif: chatter was not silent
1199 
1200  // Return result
1201  return result;
1202 }
1203 
1204 
1205 /*==========================================================================
1206  = =
1207  = Private methods =
1208  = =
1209  ==========================================================================*/
1210 
1211 /***********************************************************************//**
1212  * @brief Initialise class members
1213  ***************************************************************************/
1215 {
1216  // Initialise node energies and values
1217  m_energies.clear();
1218  m_values.clear();
1219 
1220  // Initialise evaluation cache
1221  m_old_energies.clear();
1222  m_old_values.clear();
1224  m_log_values.clear();
1225 
1226  // Initialise flux computation cache
1228  m_lin_values.clear();
1229  m_prefactor.clear();
1230  m_gamma.clear();
1231  m_epivot.clear();
1232  m_flux.clear();
1233  m_eflux.clear();
1234 
1235  // Initialise MC cache
1236  m_mc_emin.clear();
1237  m_mc_emax.clear();
1238  m_mc_cum.clear();
1239  m_mc_min.clear();
1240  m_mc_max.clear();
1241  m_mc_exp.clear();
1242 
1243  // Update parameter mapping
1244  update_pars();
1245 
1246  // Return
1247  return;
1248 }
1249 
1250 
1251 /***********************************************************************//**
1252  * @brief Copy class members
1253  *
1254  * @param[in] model Spectral function model.
1255  ***************************************************************************/
1257 {
1258  // Copy node energies and values
1259  m_energies = model.m_energies;
1260  m_values = model.m_values;
1261 
1262  // Copy evaluation cache
1264  m_old_values = model.m_old_values;
1266  m_log_values = model.m_log_values;
1267 
1268  // Copy flux computation cache
1270  m_lin_values = model.m_lin_values;
1271  m_prefactor = model.m_prefactor;
1272  m_gamma = model.m_gamma;
1273  m_epivot = model.m_epivot;
1274  m_flux = model.m_flux;
1275  m_eflux = model.m_eflux;
1276 
1277  // Copy MC cache
1278  m_mc_emin = model.m_mc_emin;
1279  m_mc_emax = model.m_mc_emax;
1280  m_mc_cum = model.m_mc_cum;
1281  m_mc_min = model.m_mc_min;
1282  m_mc_max = model.m_mc_max;
1283  m_mc_exp = model.m_mc_exp;
1284 
1285  // Update parameter mapping
1286  update_pars();
1287 
1288  // Return
1289  return;
1290 }
1291 
1292 
1293 /***********************************************************************//**
1294  * @brief Delete class members
1295  ***************************************************************************/
1297 {
1298  // Return
1299  return;
1300 }
1301 
1302 
1303 /***********************************************************************//**
1304  * @brief Update parameter mapping
1305  *
1306  * Sets the parameter pointers in the m_pars array, enabling iterating over
1307  * all model parameters. This method needs to be called after changing the
1308  * number of nodes in the spectral model. The method needs not to be called
1309  * after value update.
1310  ***************************************************************************/
1312 {
1313  // Clear parameter pointer(s)
1314  m_pars.clear();
1315 
1316  // Get number of nodes
1317  int nodes = m_energies.size();
1318 
1319  // Set parameter pointers for all nodes
1320  for (int i = 0; i < nodes; ++i) {
1321 
1322  // Set parameter names
1323  std::string energy_name = "Energy"+gammalib::str(i);
1324  std::string intensity_name = "Intensity"+gammalib::str(i);
1325 
1326  // Set energy attributes
1327  m_energies[i].name(energy_name);
1328  m_energies[i].unit("MeV");
1329  m_energies[i].has_grad(false);
1330 
1331  // Set intensity attributes
1332  m_values[i].name(intensity_name);
1333  m_values[i].unit("ph/cm2/s/MeV");
1334  m_values[i].has_grad(true);
1335 
1336  // Set pointer
1337  m_pars.push_back(&(m_energies[i]));
1338  m_pars.push_back(&(m_values[i]));
1339 
1340  } // endfor: looped over nodes
1341 
1342  // Return
1343  return;
1344 }
1345 
1346 
1347 /***********************************************************************//**
1348  * @brief Set pre-computation cache
1349  ***************************************************************************/
1351 {
1352  // Set evaluation cache
1353  set_eval_cache();
1354 
1355  // Set flux computation cache
1356  set_flux_cache();
1357 
1358  // Return
1359  return;
1360 }
1361 
1362 
1363 /***********************************************************************//**
1364  * @brief Set evaluation cache
1365  *
1366  * The evaluation cache contains pre-computed results that are needed for
1367  * fast function evaluation.
1368  *
1369  * @todo Check that all energies and intensities are > 0
1370  ***************************************************************************/
1372 {
1373  // Clear any existing values
1374  m_old_energies.clear();
1375  m_old_values.clear();
1377  m_log_values.clear();
1378 
1379  // Compute log10 of energies and intensities
1380  for (int i = 0; i < m_energies.size(); ++i) {
1381 
1382  // Set log10(energy)
1383  double log10energy = std::log10(m_energies[i].value());
1384  m_old_energies.push_back(log10energy);
1385  m_log_energies.append(log10energy);
1386 
1387  // Set log10(intensity)
1388  double log10value = std::log10(m_values[i].value());
1389  m_old_values.push_back(log10value);
1390  m_log_values.push_back(log10value);
1391 
1392  }
1393 
1394  // Return
1395  return;
1396 }
1397 
1398 
1399 /***********************************************************************//**
1400  * @brief Set flux computation cache
1401  *
1402  * Pre-computes some values that are needed for flux computation.
1403  *
1404  * @todo Handle special case emin=emax and fmin=fmax
1405  ***************************************************************************/
1407 {
1408  // Clear any existing values
1410  m_lin_values.clear();
1411  m_prefactor.clear();
1412  m_gamma.clear();
1413  m_epivot.clear();
1414  m_flux.clear();
1415  m_eflux.clear();
1416 
1417  // Store linear energies and values
1418  for (int i = 0; i < m_energies.size(); ++i) {
1419  m_lin_energies.append(m_energies[i].value());
1420  m_lin_values.push_back(m_values[i].value());
1421  }
1422 
1423  // Loop over all nodes-1
1424  for (int i = 0; i < m_energies.size()-1; ++i) {
1425 
1426  // Get energies and function values
1427  double emin = m_lin_energies[i];
1428  double emax = m_lin_energies[i+1];
1429  double fmin = m_values[i].value();
1430  double fmax = m_values[i+1].value();
1431 
1432  // Compute pivot energy (MeV). We use here the geometric mean of
1433  // the node boundaries.
1434  double epivot = std::sqrt(emin*emax);
1435 
1436  // Compute spectral index
1437  double gamma = std::log(fmin/fmax) / std::log(emin/emax);
1438 
1439  // Compute power law normalisation
1440  double prefactor = fmin / std::pow(emin/epivot, gamma);
1441 
1442  // Compute photon flux between nodes
1443  double flux = prefactor *
1444  gammalib::plaw_photon_flux(emin, emax, epivot, gamma);
1445 
1446  // Compute energy flux between nodes
1447  double eflux = prefactor *
1448  gammalib::plaw_energy_flux(emin, emax, epivot, gamma);
1449 
1450  // Convert energy flux from MeV/cm2/s to erg/cm2/s
1451  eflux *= gammalib::MeV2erg;
1452 
1453  // Push back values on pre-computation cache
1454  m_prefactor.push_back(prefactor);
1455  m_gamma.push_back(gamma);
1456  m_epivot.push_back(epivot);
1457  m_flux.push_back(flux);
1458  m_eflux.push_back(eflux);
1459 
1460  } // endfor: looped over all nodes
1461 
1462  // Return
1463  return;
1464 }
1465 
1466 
1467 /***********************************************************************//**
1468  * @brief Update evaluation cache
1469  *
1470  * Updates the evaluation cache by computing only results for values that
1471  * changed.
1472  *
1473  * @todo Check that all energies and intensities are > 0
1474  ***************************************************************************/
1476 {
1477  // Update energies
1478  for (int i = 0; i < m_energies.size(); ++i) {
1479  double energy = m_energies[i].value();
1480  if (energy != m_old_energies[i]) {
1481  m_log_energies[i] = std::log10(energy);
1482  m_old_energies[i] = energy;
1483  }
1484  }
1485 
1486  // Update intensities
1487  for (int i = 0; i < m_values.size(); ++i) {
1488  double value = m_values[i].value();
1489  if (value != m_old_values[i]) {
1490  m_log_values[i] = std::log10(value);
1491  m_old_values[i] = value;
1492  }
1493  }
1494 
1495  // Return
1496  return;
1497 }
1498 
1499 
1500 /***********************************************************************//**
1501  * @brief Update flux computation cache
1502  *
1503  * Updates the flux computation cache if either the energy boundaries or the
1504  * intensity values have changed.
1505  ***************************************************************************/
1507 {
1508  // Loop over all nodes. If energy or value of a node has changed then
1509  // set cache
1510  for (int i = 0; i < m_energies.size(); ++i) {
1511  if ((m_lin_energies[i] != m_energies[i].value()) ||
1512  (m_lin_values[i] != m_values[i].value())) {
1513  set_flux_cache();
1514  break;
1515  }
1516  }
1517 
1518  // Return
1519  return;
1520 }
1521 
1522 
1523 /***********************************************************************//**
1524  * @brief Set MC pre-computation cache
1525  *
1526  * @param[in] emin Minimum energy.
1527  * @param[in] emax Maximum energy.
1528  *
1529  * This method sets up an array of indices and the cumulative distribution
1530  * function needed for MC simulations.
1531  ***************************************************************************/
1532 void GModelSpectralNodes::mc_update(const GEnergy& emin, const GEnergy& emax) const
1533 {
1534  // Check if we need to update the cache
1535  if (emin != m_mc_emin || emax != m_mc_emax) {
1536 
1537  // Store new energy interval
1538  m_mc_emin = emin;
1539  m_mc_emax = emax;
1540 
1541  // Initialise cache
1542  m_mc_cum.clear();
1543  m_mc_min.clear();
1544  m_mc_max.clear();
1545  m_mc_exp.clear();
1546 
1547  // Get energy range in MeV
1548  double e_min = emin.MeV();
1549  double e_max = emax.MeV();
1550 
1551  // Continue only if e_max > e_min
1552  if (e_max > e_min) {
1553 
1554  // Allocate flux
1555  double flux;
1556 
1557  // Determine left node index for minimum energy
1558  m_lin_energies.set_value(e_min);
1559  int inx_emin = m_lin_energies.inx_left();
1560 
1561  // Determine left node index for maximum energy
1562  m_lin_energies.set_value(e_max);
1563  int inx_emax = m_lin_energies.inx_left();
1564 
1565  // If both energies are within the same node then just
1566  // add this one node on the stack
1567  if (inx_emin == inx_emax) {
1568  flux = m_prefactor[inx_emin] *
1570  e_max,
1571  m_epivot[inx_emin],
1572  m_gamma[inx_emin]);
1573  m_mc_cum.push_back(flux);
1574  m_mc_min.push_back(e_min);
1575  m_mc_max.push_back(e_max);
1576  m_mc_exp.push_back(m_gamma[inx_emin]);
1577  }
1578 
1579  // ... otherwise integrate over the nodes where emin and emax
1580  // resides and all the remaining nodes
1581  else {
1582 
1583  // If we are requested to extrapolate beyond first node,
1584  // the use the first nodes lower energy and upper energy
1585  // boundary for the initial integration.
1586  int i_start = (e_min < m_lin_energies[0]) ? inx_emin : inx_emin+1;
1587 
1588  // Add emin to the node boundary
1589  flux = m_prefactor[inx_emin] *
1591  m_lin_energies[i_start],
1592  m_epivot[inx_emin],
1593  m_gamma[inx_emin]);
1594  m_mc_cum.push_back(flux);
1595  m_mc_min.push_back(e_min);
1596  m_mc_max.push_back(m_lin_energies[i_start]);
1597  m_mc_exp.push_back(m_gamma[inx_emin]);
1598 
1599  // Add all nodes between
1600  for (int i = i_start; i < inx_emax; ++i) {
1601  flux = m_flux[i];
1602  m_mc_cum.push_back(flux);
1603  m_mc_min.push_back(m_lin_energies[i]);
1604  m_mc_max.push_back(m_lin_energies[i+1]);
1605  m_mc_exp.push_back(m_gamma[i]);
1606  }
1607 
1608  // Add node boundary to emax
1609  flux = m_prefactor[inx_emax] *
1611  e_max,
1612  m_epivot[inx_emax],
1613  m_gamma[inx_emax]);
1614  m_mc_cum.push_back(flux);
1615  m_mc_min.push_back(m_lin_energies[inx_emax]);
1616  m_mc_max.push_back(e_max);
1617  m_mc_exp.push_back(m_gamma[inx_emax]);
1618 
1619  } // endelse: emin and emax not between same nodes
1620 
1621  // Build cumulative distribution
1622  for (int i = 1; i < m_mc_cum.size(); ++i) {
1623  m_mc_cum[i] += m_mc_cum[i-1];
1624  }
1625  double norm = m_mc_cum[m_mc_cum.size()-1];
1626  if (norm > 0.0) {
1627  for (int i = 0; i < m_mc_cum.size(); ++i) {
1628  m_mc_cum[i] /= norm;
1629  }
1630  }
1631 
1632  // Set MC values
1633  for (int i = 0; i < m_mc_cum.size(); ++i) {
1634 
1635  // Compute exponent
1636  double exponent = m_mc_exp[i] + 1.0;
1637 
1638  // Exponent dependend computation
1639  if (std::abs(exponent) > 1.0e-11) {
1640 
1641  // If the exponent is too small then use lower energy
1642  // boundary
1643  if (exponent < -50.0) {
1644  m_mc_exp[i] = 0.0;
1645  m_mc_min[i] = std::log(m_mc_min[i]);
1646  m_mc_max[i] = m_mc_min[i];
1647  }
1648 
1649  // ... otherwise if exponent is too large then use
1650  // upper energy boundary
1651  else if (exponent > +50.0) {
1652  m_mc_exp[i] = 0.0;
1653  m_mc_min[i] = std::log(m_mc_max[i]);
1654  m_mc_max[i] = m_mc_min[i];
1655  }
1656 
1657  // ... otherwise use transformation formula
1658  else {
1659  m_mc_exp[i] = exponent;
1660  m_mc_min[i] = std::pow(m_mc_min[i], exponent);
1661  m_mc_max[i] = std::pow(m_mc_max[i], exponent);
1662  }
1663  }
1664  else {
1665  m_mc_exp[i] = 0.0;
1666  m_mc_min[i] = std::log(m_mc_min[i]);
1667  m_mc_max[i] = std::log(m_mc_max[i]);
1668  }
1669 
1670  } // endfor: set MC values
1671 
1672  } // endif: e_max > e_min
1673 
1674  } // endif: Update was required
1675 
1676  // Return
1677  return;
1678 }
std::vector< double > m_flux
Photon fluxes.
void update_flux_cache(void) const
Update flux computation cache.
virtual GEnergy mc(const GEnergy &emin, const GEnergy &emax, const GTime &time, GRan &ran) const
Returns MC energy between [emin, emax].
std::vector< double > m_eflux
Energy fluxes.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:821
virtual GModelSpectralNodes & operator=(const GModelSpectralNodes &model)
Assignment operator.
const std::string & name(void) const
Return parameter name.
Random number generator class definition.
#define G_INTENSITY_GET
Abstract spectral model base class.
virtual std::string type(void) const
Return model type.
int size(void) const
Return number of parameters.
void insert(const int &index, const GEnergy &energy, const double &intensity)
Insert node.
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1163
virtual GModelSpectral & operator=(const GModelSpectral &model)
Assignment operator.
void extend(const GModelSpectralNodes &nodes)
Append nodes from node function.
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:353
std::vector< double > m_mc_max
Upper boundary for MC.
std::vector< double > m_old_energies
Old energies.
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const =0
std::vector< GModelPar * > m_pars
Parameter pointers.
XML element node class.
Definition: GXmlElement.hpp:47
const double & wgt_left(void) const
Returns left node weight.
Definition: GNodeArray.hpp:260
Spectral model registry class definition.
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
Random number generator class.
Definition: GRan.hpp:44
virtual double flux(const GEnergy &emin, const GEnergy &emax) const
Returns model photon flux between emin, emax
virtual double eflux(const GEnergy &emin, const GEnergy &emax) const
Returns model energy flux between emin, emax
Spectral nodes model class definition.
std::vector< double > m_epivot
Power-law pivot energies.
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
virtual int elements(void) const
Return number of GXMLElement children of node.
Definition: GXmlNode.cpp:580
Time class.
Definition: GTime.hpp:54
std::vector< double > m_mc_cum
Cumulative distribution.
Gammalib tools definition.
void reserve(const int &num)
Reserve space for nodes.
void set_value(const double &value) const
Set indices and weighting factors for interpolation.
Definition: GNodeArray.cpp:580
GNodeArray m_lin_energies
Energy of nodes.
double plaw_energy_flux(const double &emin, const double &emax, const double &epivot, const double &gamma)
Compute energy flux between two energies for a power law.
Definition: GTools.cpp:1127
Energy container class.
Definition: GEnergies.hpp:60
double log10MeV(void) const
Return log10 of energy in MeV.
Definition: GEnergy.cpp:423
Spectral nodes model class.
std::vector< double > m_log_values
log10(value) of nodes
std::vector< GModelPar > m_values
Node values.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:184
void set_flux_cache(void) const
Set flux computation cache.
#define G_WRITE
virtual std::string print(const GChatter &chatter=NORMAL) const
Print node function information.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
void append(const GEnergy &energy, const double &intensity)
Append node.
Model parameter class.
Definition: GModelPar.hpp:87
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
const double & wgt_right(void) const
Returns right node weight.
Definition: GNodeArray.hpp:275
virtual GModelSpectralNodes * clone(void) const
Clone spectral nodes model.
const GXmlAttribute * attribute(const int &index) const
Return attribute.
void copy_members(const GModelSpectralNodes &model)
Copy class members.
std::string print(const GChatter &chatter=NORMAL) const
Print energy.
Definition: GEnergy.cpp:748
void free(void)
Free a parameter.
virtual void write(GXmlElement &xml) const
Write model into XML element.
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const
Evaluate function.
Energy container class definition.
void fix(void)
Fix a parameter.
int size(void) const
Return number of energies in container.
Definition: GEnergies.hpp:162
GXmlElement * xml_need_par(const std::string &origin, GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1513
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1184
GModelSpectralNodes(void)
Void constructor.
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
virtual ~GModelSpectralNodes(void)
Destructor.
std::vector< double > m_mc_exp
Exponent for MC.
void update_pars(void)
Update parameter mapping.
void set_eval_cache(void) const
Set evaluation cache.
GChatter
Definition: GTypemaps.hpp:33
GEnergy m_mc_emax
Maximum energy.
const int & inx_left(void) const
Returns left node index.
Definition: GNodeArray.hpp:231
GEnergy energy(const int &index) const
Return node energy.
std::vector< double > m_lin_values
Values of nodes.
#define G_INTENSITY_SET
const GModelSpectralNodes g_spectral_nodes_seed
Interface definition for the spectral model registry class.
const int & inx_right(void) const
Returns right node index.
Definition: GNodeArray.hpp:245
#define G_MC
#define G_ENERGY_GET
double intensity(const int &index) const
Return node intensity.
std::vector< double > m_mc_min
Lower boundary for MC.
void free_members(void)
Delete class members.
virtual GXmlElement * element(const int &index)
Return pointer to GXMLElement child.
Definition: GXmlNode.cpp:634
void update_eval_cache(void) const
Update evaluation cache.
virtual void clear(void)
Clear spectral nodes model.
void init_members(void)
Initialise class members.
std::vector< GModelPar > m_energies
Node energies.
void init_members(void)
Initialise class members.
const double MeV2erg
Definition: GTools.hpp:45
int nodes(void) const
Return number of nodes.
void read(const GXmlElement &xml)
Extract parameter attributes from XML element.
Definition: GModelPar.cpp:229
double value(void) const
Return parameter value.
const std::string & unit(void) const
Return parameter unit.
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1332
std::vector< double > m_prefactor
Power-law normalisations.
GEnergy m_mc_emin
Minimum energy.
Exception handler interface definition.
virtual void read(const GXmlElement &xml)
Read model from XML element.
void mc_update(const GEnergy &emin, const GEnergy &emax) const
Set MC pre-computation cache.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
#define G_REMOVE
#define G_INSERT
virtual GXmlNode * append(const GXmlNode &node)
Append XML child node.
Definition: GXmlNode.cpp:284
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
#define G_READ
void append(const double &node)
Append one node to array.
Definition: GNodeArray.cpp:351
#define G_APPEND
void remove(const int &index)
Remove node.
const GXmlElement * xml_get_par(const std::string &origin, const GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1562
void free_members(void)
Delete class members.
void clear(void)
Clear instance.
Definition: GEnergy.cpp:261
void set_cache(void) const
Set pre-computation cache.
#define G_ENERGY_SET
std::vector< double > m_old_values
Old values.
GVector log10(const GVector &vector)
Computes base10 logarithm of vector elements.
Definition: GVector.cpp:1205
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
double plaw_photon_flux(const double &emin, const double &emax, const double &epivot, const double &gamma)
Compute photon flux between two energies for a power law.
Definition: GTools.cpp:1079
std::vector< double > m_gamma
Power-law indices.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413
GNodeArray m_log_energies
log10(energy) of nodes