GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpectralFunc.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpectralFunc.cpp - Spectral function model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2010-2022 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 GModelSpectralFunc.cpp
23  * @brief Spectral function 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 "GCsv.hpp"
35 #include "GRan.hpp"
36 #include "GEnergies.hpp"
37 #include "GModelSpectralFunc.hpp"
39 
40 /* __ Constants __________________________________________________________ */
41 
42 /* __ Globals ____________________________________________________________ */
44 const GModelSpectralRegistry g_spectral_func_registry(&g_spectral_func_seed);
45 
46 /* __ Method name definitions ____________________________________________ */
47 #define G_FLUX "GModelSpectralFunc::flux(GEnergy&, GEnergy&)"
48 #define G_EFLUX "GModelSpectralFunc::eflux(GEnergy&, GEnergy&)"
49 #define G_MC "GModelSpectralFunc::mc(GEnergy&, GEnergy&, GTime&, GRan&)"
50 #define G_READ "GModelSpectralFunc::read(GXmlElement&)"
51 #define G_WRITE "GModelSpectralFunc::write(GXmlElement&)"
52 #define G_APPEND "GModelSpectralFunc::append(GEnergy&, double&)"
53 #define G_INSERT "GModelSpectralFunc::insert(GEnergy&, double&)"
54 #define G_REMOVE "GModelSpectralFunc::remove(int&)"
55 #define G_ENERGY1 "GModelSpectralFunc::energy(int&)"
56 #define G_ENERGY2 "GModelSpectralFunc::energy(int&, GEnergy&)"
57 #define G_INTENSITY1 "GModelSpectralFunc::intensity(int&)"
58 #define G_INTENSITY2 "GModelSpectralFunc::intensity(int&, double&)"
59 #define G_LOAD_NODES "GModelSpectralFunc::load_nodes(GFilename&)"
60 
61 /* __ Macros _____________________________________________________________ */
62 
63 /* __ Coding definitions _________________________________________________ */
64 
65 /* __ Debug definitions __________________________________________________ */
66 
67 
68 /*==========================================================================
69  = =
70  = Constructors/destructors =
71  = =
72  ==========================================================================*/
73 
74 /***********************************************************************//**
75  * @brief Void constructor
76  ***************************************************************************/
78 {
79  // Initialise members
80  init_members();
81 
82  // Return
83  return;
84 }
85 
86 
87 /***********************************************************************//**
88  * @brief File constructor
89  *
90  * @param[in] filename File name of nodes.
91  * @param[in] norm Normalization factor.
92  *
93  * Constructs spectral file function model from a list of nodes that is found
94  * in the specified ASCII file. See the load_nodes() method for more
95  * information about the expected structure of the file.
96  ***************************************************************************/
98  const double& norm) :
100 {
101  // Initialise members
102  init_members();
103 
104  // Load nodes
105  load_nodes(filename);
106 
107  // Set normalization
108  m_norm.value(norm);
109 
110  // Return
111  return;
112 }
113 
114 
115 /***********************************************************************//**
116  * @brief XML constructor
117  *
118  * @param[in] xml XML element.
119  *
120  * Constructs spectral file function model by extracting information from an
121  * XML element. See the read() method for more information about the expected
122  * structure of the XML element.
123  ***************************************************************************/
126 {
127  // Initialise members
128  init_members();
129 
130  // Read information from XML element
131  read(xml);
132 
133  // Return
134  return;
135 }
136 
137 
138 /***********************************************************************//**
139  * @brief Spectral model constructor
140  *
141  * @param[in] model Spectral model.
142  * @param[in] energies File function node energies.
143  *
144  * Constructs a spectral file function model from any spectral model.
145  * The file function normalisation will be set to unity.
146  ***************************************************************************/
148  const GEnergies& energies) :
150 {
151  // Initialise members
152  init_members();
153 
154  // Reserve space for file function nodes
155  reserve(energies.size());
156 
157  // Append nodes for all energies
158  for (int i = 0; i < energies.size(); ++i) {
159  append(energies[i], model.eval(energies[i]));
160  }
161 
162  // Return
163  return;
164 }
165 
166 
167 /***********************************************************************//**
168  * @brief Copy constructor
169  *
170  * @param[in] model File function model.
171  ***************************************************************************/
173  GModelSpectral(model)
174 {
175  // Initialise members
176  init_members();
177 
178  // Copy members
179  copy_members(model);
180 
181  // Return
182  return;
183 }
184 
185 
186 /***********************************************************************//**
187  * @brief Destructor
188  ***************************************************************************/
190 {
191  // Free members
192  free_members();
193 
194  // Return
195  return;
196 }
197 
198 
199 /*==========================================================================
200  = =
201  = Operators =
202  = =
203  ==========================================================================*/
204 
205 /***********************************************************************//**
206  * @brief Assignment operator
207  *
208  * @param[in] model File function model.
209  * @return File function model.
210  ***************************************************************************/
212 {
213  // Execute only if object is not identical
214  if (this != &model) {
215 
216  // Copy base class members
217  this->GModelSpectral::operator=(model);
218 
219  // Free members
220  free_members();
221 
222  // Initialise members
223  init_members();
224 
225  // Copy members
226  copy_members(model);
227 
228  } // endif: object was not identical
229 
230  // Return
231  return *this;
232 }
233 
234 
235 /*==========================================================================
236  = =
237  = Public methods =
238  = =
239  ==========================================================================*/
240 
241 /***********************************************************************//**
242  * @brief Clear file function
243 ***************************************************************************/
245 {
246  // Free class members (base and derived classes, derived class first)
247  free_members();
249 
250  // Initialise members
252  init_members();
253 
254  // Return
255  return;
256 }
257 
258 
259 /***********************************************************************//**
260  * @brief Clone file function
261 ***************************************************************************/
263 {
264  // Clone file function
265  return new GModelSpectralFunc(*this);
266 }
267 
268 
269 /***********************************************************************//**
270  * @brief Evaluate function
271  *
272  * @param[in] srcEng True photon energy.
273  * @param[in] srcTime True photon arrival time.
274  * @param[in] gradients Compute gradients?
275  * @return Model value (ph/cm2/s/MeV).
276  *
277  * Evaluates
278  *
279  * \f[
280  * S_{\rm E}(E | t) = {\tt m\_norm} F(E)
281  * \f]
282  *
283  * where
284  * - \f${\tt m\_norm}\f$ is the normalization factor and
285  * - \f${F(E)}\f$ is the spectral function.
286  *
287  * If the @p gradients flag is true the method will also compute the
288  * partial derivatives of the model with respect to the parameters using
289  *
290  * \f[
291  * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_norm}} =
292  * \frac{S_{\rm E}(E | t)}{{\tt m\_norm}}
293  * \f]
294  ***************************************************************************/
295 double GModelSpectralFunc::eval(const GEnergy& srcEng,
296  const GTime& srcTime,
297  const bool& gradients) const
298 {
299  // Interpolate function. This is done in log10-log10 space, but the
300  // linear value is returned.
301  double arg = m_log_nodes.interpolate(srcEng.log10MeV(), m_log_values);
302  double func = std::pow(10.0, arg);
303 
304  // Compute function value
305  double value = m_norm.value() * func;
306 
307  // Optionally compute gradients
308  if (gradients) {
309 
310  // Compute partial derivatives of the parameter values
311  double g_norm = (m_norm.is_free()) ? m_norm.scale() * func : 0.0;
312 
313  // Set gradients
314  m_norm.factor_gradient(g_norm);
315 
316  } // endif: gradient computation was requested
317 
318  // Compile option: Check for NaN/Inf
319  #if defined(G_NAN_CHECK)
320  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
321  std::cout << "*** ERROR: GModelSpectralFunc::eval";
322  std::cout << "(srcEng=" << srcEng;
323  std::cout << ", srcTime=" << srcTime << "):";
324  std::cout << " NaN/Inf encountered";
325  std::cout << " (value=" << value;
326  std::cout << ", norm=" << norm();
327  std::cout << ", func=" << func;
328  std::cout << ")" << std::endl;
329  }
330  #endif
331 
332  // Return
333  return value;
334 }
335 
336 
337 /***********************************************************************//**
338  * @brief Returns model photon flux between [emin, emax] (units: ph/cm2/s)
339  *
340  * @param[in] emin Minimum photon energy.
341  * @param[in] emax Maximum photon energy.
342  * @return Photon flux (ph/cm2/s).
343  *
344  * Computes
345  *
346  * \f[
347  * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) dE
348  * \f]
349  *
350  * where
351  * - [@p emin, @p emax] is an energy interval, and
352  * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
353  ***************************************************************************/
354 double GModelSpectralFunc::flux(const GEnergy& emin, const GEnergy& emax) const
355 {
356  // Initialise flux
357  double flux = 0.0;
358 
359  // Compute only if integration range is valid
360  if (emin < emax) {
361 
362  // Get energy range in MeV
363  double e_min = emin.MeV();
364  double e_max = emax.MeV();
365 
366  // Determine left node index for minimum energy
367  m_lin_nodes.set_value(e_min);
368  int inx_emin = m_lin_nodes.inx_left();
369 
370  // Determine left node index for maximum energy
371  m_lin_nodes.set_value(e_max);
372  int inx_emax = m_lin_nodes.inx_left();
373 
374  // If both energies are within the same nodes then simply
375  // integrate over the energy interval using the appropriate power
376  // law parameters
377  if (inx_emin == inx_emax) {
378  flux = m_prefactor[inx_emin] *
380  e_max,
381  m_epivot[inx_emin],
382  m_gamma[inx_emin]);
383  }
384 
385  // ... otherwise integrate over the nodes where emin and emax
386  // resides and all the remaining nodes
387  else {
388 
389  // If we are requested to extrapolate beyond first node,
390  // the use the first nodes lower energy and upper energy
391  // boundary for the initial integration.
392  int i_start = (e_min < m_lin_nodes[0]) ? inx_emin : inx_emin+1;
393 
394  // Integrate from emin to the node boundary
395  flux = m_prefactor[inx_emin] *
397  m_lin_nodes[i_start],
398  m_epivot[inx_emin],
399  m_gamma[inx_emin]);
400 
401  // Integrate over all nodes between
402  for (int i = i_start; i < inx_emax; ++i) {
403  flux += m_flux[i];
404  }
405 
406  // Integrate from node boundary to emax
407  flux += m_prefactor[inx_emax] *
409  e_max,
410  m_epivot[inx_emax],
411  m_gamma[inx_emax]);
412 
413  } // endelse: emin and emax not between same nodes
414 
415  // Multiply flux by normalisation factor
416  flux *= norm();
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 GModelSpectralFunc::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  // Get energy range in MeV
451  double e_min = emin.MeV();
452  double e_max = emax.MeV();
453 
454  // Determine left node index for minimum energy
455  m_lin_nodes.set_value(e_min);
456  int inx_emin = m_lin_nodes.inx_left();
457 
458  // Determine left node index for maximum energy
459  m_lin_nodes.set_value(e_max);
460  int inx_emax = m_lin_nodes.inx_left();
461 
462  // If both energies are within the same nodes then simply
463  // integrate over the energy interval using the appropriate power
464  // law parameters
465  if (inx_emin == inx_emax) {
466  eflux = m_prefactor[inx_emin] *
468  e_max,
469  m_epivot[inx_emin],
470  m_gamma[inx_emin]) *
472  }
473 
474  // ... otherwise integrate over the nodes where emin and emax
475  // resides and all the remaining nodes
476  else {
477 
478  // If we are requested to extrapolate beyond first node,
479  // the use the first nodes lower energy and upper energy
480  // boundary for the initial integration.
481  int i_start = (e_min < m_lin_nodes[0]) ? inx_emin : inx_emin+1;
482 
483  // Integrate from emin to the node boundary
484  eflux = m_prefactor[inx_emin] *
486  m_lin_nodes[i_start],
487  m_epivot[inx_emin],
488  m_gamma[inx_emin]) *
490 
491  // Integrate over all nodes between
492  for (int i = i_start; i < inx_emax; ++i) {
493  eflux += m_eflux[i];
494  }
495 
496  // Integrate from node boundary to emax
497  eflux += m_prefactor[inx_emax] *
499  e_max,
500  m_epivot[inx_emax],
501  m_gamma[inx_emax]) *
503 
504  } // endelse: emin and emax not between same nodes
505 
506  // Multiply flux by normalisation factor
507  eflux *= norm();
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  *
528  * Returns Monte Carlo energy by randomly drawing from a broken power law
529  * defined by the file function.
530  ***************************************************************************/
532  const GEnergy& emax,
533  const GTime& time,
534  GRan& ran) const
535 {
536  // Check energy interval
538 
539  // Allocate energy
540  GEnergy energy;
541 
542  // Update cache
543  mc_update(emin, emax);
544 
545  // Determine in which bin we reside
546  int inx = 0;
547  if (m_mc_cum.size() > 1) {
548  double u = ran.uniform();
549  for (inx = m_mc_cum.size()-1; inx > 0; --inx) {
550  if (m_mc_cum[inx-1] <= u) {
551  break;
552  }
553  }
554  }
555 
556  // Get random energy for specific bin
557  if (m_mc_exp[inx] != 0.0) {
558  double e_min = m_mc_min[inx];
559  double e_max = m_mc_max[inx];
560  double u = ran.uniform();
561  double eng = (u > 0.0)
562  ? std::exp(std::log(u * (e_max - e_min) + e_min) / m_mc_exp[inx])
563  : 0.0;
564  energy.MeV(eng);
565  }
566  else {
567  double e_min = m_mc_min[inx];
568  double e_max = m_mc_max[inx];
569  double u = ran.uniform();
570  double eng = std::exp(u * (e_max - e_min) + e_min);
571  energy.MeV(eng);
572  }
573 
574  // Return energy
575  return energy;
576 }
577 
578 
579 /***********************************************************************//**
580  * @brief Read model from XML element
581  *
582  * @param[in] xml XML element containing power law model information.
583  *
584  * Reads the spectral information from an XML element. The format of the XML
585  * elements is
586  *
587  * <spectrum type="FileFunction" file="..">
588  * <parameter name="Normalization" scale=".." value=".." min=".." max=".." free=".."/>
589  * </spectrum>
590  ***************************************************************************/
592 {
593  // Verify number of model parameters
595 
596  // Get parameter
598 
599  // Read parameter
600  m_norm.read(*norm);
601 
602  // Load nodes from file
604 
605  // Return
606  return;
607 }
608 
609 
610 /***********************************************************************//**
611  * @brief Write model into XML element
612  *
613  * @param[in] xml XML element into which model information is written.
614  *
615  * Writes the spectral information into an XML element. The format of the XML
616  * element is
617  *
618  * <spectrum type="FileFunction" file="..">
619  * <parameter name="Normalization" scale=".." value=".." min=".." max=".." free=".."/>
620  * </spectrum>
621  *
622  * Note that the function nodes will not be written since they will not be
623  * altered by any method.
624  ***************************************************************************/
626 {
627  // Verify model type
629 
630  // Get or create parameter
632 
633  // Write parameter
634  m_norm.write(*norm);
635 
636  // Set file attribute
638 
639  // Return
640  return;
641 }
642 
643 
644 /***********************************************************************//**
645  * @brief Append node to file function
646  *
647  * @param[in] energy Energy.
648  * @param[in] intensity Intensity.
649  *
650  * @exception GException::invalid_argument
651  * Energy not larger than energy of last node or energy or
652  * intensity are not positive.
653  *
654  * Appends one node to the file function.
655  ***************************************************************************/
656 void GModelSpectralFunc::append(const GEnergy& energy, const double& intensity)
657 {
658  // Check that energy is larger than energy of last node
659  if (nodes() > 0) {
660  if (energy.MeV() <= m_lin_nodes[nodes()-1]) {
661  GEnergy last(m_lin_nodes[nodes()-1], "MeV");
662  std::string msg = "Specified energy "+energy.print()+" is not "
663  "larger than the energy "+last.print()+" of the "
664  "last node of the file function. Please append "
665  "nodes in increasing order of energy.";
667  }
668  }
669 
670  // Check that energy is positive
671  if (energy.MeV() <= 0.0) {
672  std::string msg = "Specified energy "+energy.print()+" is not "
673  "positive. Please append only nodes with positive "
674  "energies.";
676  }
677 
678  // Check that intensity is positive
679  if (intensity <= 0.0) {
680  std::string msg = "Specified intensity "+gammalib::str(intensity)+" is "
681  "not positive. Please append only nodes with positive "
682  "intensities.";
684  }
685 
686  // Append node
687  m_lin_nodes.append(energy.MeV());
688  m_log_nodes.append(energy.log10MeV());
689  m_lin_values.push_back(intensity);
690  m_log_values.push_back(std::log10(intensity));
691 
692  // Set pre-computation cache
693  set_cache();
694 
695  // Return
696  return;
697 }
698 
699 
700 /***********************************************************************//**
701  * @brief Insert node into file function
702  *
703  * @param[in] energy Energy.
704  * @param[in] intensity Intensity.
705  *
706  * @exception GException::invalid_argument
707  * Energy or intensity are not positive, or energy does already
708  * exist in file function.
709  *
710  * Inserts one node at the appropriate location into the file function.
711  ***************************************************************************/
712 void GModelSpectralFunc::insert(const GEnergy& energy, const double& intensity)
713 {
714  // Get energy in MeV
715  double energy_MeV = energy.MeV();
716 
717  // Check that energy is positive
718  if (energy_MeV <= 0.0) {
719  std::string msg = "Specified energy "+energy.print()+" is not "
720  "positive. Please append only nodes with positive "
721  "energies.";
723  }
724 
725  // Check that intensity is positive
726  if (intensity <= 0.0) {
727  std::string msg = "Specified intensity "+gammalib::str(intensity)+" is "
728  "not positive. Please append only nodes with positive "
729  "intensities.";
731  }
732 
733  // Find index before which the node should be inserted
734  int index = 0;
735  for (int i = 0; i < nodes(); ++i) {
736  if (m_lin_nodes[i] == energy_MeV) {
737  std::string msg = "A node with the specified energy "+
738  energy.print()+" exists already in the file "
739  "function. Please insert only nodes with "
740  "energies that do not yet exist.";
742  }
743  else if (m_lin_nodes[i] > energy_MeV) {
744  break;
745  }
746  index++;
747  }
748 
749  // Insert node
750  m_lin_nodes.insert(index, energy_MeV);
751  m_log_nodes.insert(index, energy.log10MeV());
752  m_lin_values.insert(m_lin_values.begin()+index, intensity);
753  m_log_values.insert(m_log_values.begin()+index, std::log10(intensity));
754 
755  // Set pre-computation cache
756  set_cache();
757 
758  // Return
759  return;
760 }
761 
762 
763 /***********************************************************************//**
764  * @brief Remove node from file function
765  *
766  * @param[in] index Node index [0,...,nodes()-1].
767  *
768  * @exception GException::out_of_range
769  * Node index is out of range.
770  *
771  * Removes the node of specified @p index from the file function.
772  ***************************************************************************/
773 void GModelSpectralFunc::remove(const int& index)
774 {
775  // Compile option: raise exception if index is out of range
776  #if defined(G_RANGE_CHECK)
777  if (index < 0 || index >= nodes()) {
778  throw GException::out_of_range(G_REMOVE, "Node index", index, nodes());
779  }
780  #endif
781 
782  // Delete node
783  m_lin_nodes.remove(index);
784  m_log_nodes.remove(index);
785  m_lin_values.erase(m_lin_values.begin() + index);
786  m_log_values.erase(m_log_values.begin() + index);
787 
788  // Set pre-computation cache
789  set_cache();
790 
791  // Return
792  return;
793 }
794 
795 
796 /***********************************************************************//**
797  * @brief Append file function
798  *
799  * @param[in] filefct File function.
800  *
801  * Appends file function to the existing file function.
802  ***************************************************************************/
804 {
805  // Do nothing if file function is empty
806  if (!filefct.is_empty()) {
807 
808  // Get file function size
809  int num = filefct.nodes();
810 
811  // Reserve enough space
812  reserve(nodes() + num);
813 
814  // Append all nodes
815  for (int i = 0; i < num; ++i) {
816  append(filefct.energy(i), filefct.intensity(i));
817  }
818 
819  } // endif: file function was not empty
820 
821  // Set pre-computation cache
822  set_cache();
823 
824  // Return
825  return;
826 }
827 
828 
829 /***********************************************************************//**
830  * @brief Return energy of node
831  *
832  * @param[in] index Node index [0,...,nodes()-1].
833  *
834  * @exception GException::out_of_range
835  * Node index is out of range.
836  *
837  * Returns the energy of node with specified @p index in the file function.
838  ***************************************************************************/
839 GEnergy GModelSpectralFunc::energy(const int& index) const
840 {
841  // Compile option: raise exception if index is out of range
842  #if defined(G_RANGE_CHECK)
843  if (index < 0 || index >= nodes()) {
844  throw GException::out_of_range(G_ENERGY1, "Node index", index, nodes());
845  }
846  #endif
847 
848  // Set energy
849  GEnergy energy;
850  energy.MeV(m_lin_nodes[index]);
851 
852  // Return energy
853  return energy;
854 }
855 
856 
857 /***********************************************************************//**
858  * @brief Set energy of node
859  *
860  * @param[in] index Node index [0,...,nodes()-1].
861  * @param[in] energy Energy.
862  *
863  * @exception GException::out_of_range
864  * Node index is out of range.
865  * @exception GException::invalid_argument
866  * Specified energy is not larger than energy of preceeding node
867  * or not smaller than energy of following node.
868  *
869  * Sets the energy of node with specified @p index in the file function.
870  ***************************************************************************/
871 void GModelSpectralFunc::energy(const int& index, const GEnergy& energy)
872 {
873  // Compile option: raise exception if index is out of range
874  #if defined(G_RANGE_CHECK)
875  if (index < 0 || index >= nodes()) {
876  throw GException::out_of_range(G_ENERGY2, "Node index", index, nodes());
877  }
878  #endif
879 
880  // Check that energy is larger than energy of precedent node
881  if (index > 0) {
882  if (energy.MeV() <= m_lin_nodes[index-1]) {
883  GEnergy precedent(m_lin_nodes[index-1], "MeV");
884  std::string msg = "Specified energy "+energy.print()+" is not "
885  "larger than energy "+precedent.print()+" of "
886  "precedent node. Please specify an energy that "
887  "is larger than "+precedent.print()+".";
889  }
890  }
891 
892  // Check that energy is smaller than energy of following node
893  if (index < nodes()-1) {
894  if (energy.MeV() >= m_lin_nodes[index+1]) {
895  GEnergy following(m_lin_nodes[index+1], "MeV");
896  std::string msg = "Specified energy "+energy.print()+" is not "
897  "smaller than energy "+following.print()+" of "
898  "following node. Please specify an energy that "
899  "is smaller than "+following.print()+".";
901  }
902  }
903 
904  // Set node
905  m_lin_nodes[index] = energy.MeV();
906  m_log_nodes[index] = energy.log10MeV();
907 
908  // Set pre-computation cache
909  set_cache();
910 
911  // Return
912  return;
913 }
914 
915 
916 /***********************************************************************//**
917  * @brief Return intensity of node
918  *
919  * @param[in] index Node index [0,...,nodes()-1].
920  * @return Intensity (ph/cm2/s/MeV).
921  *
922  * @exception GException::out_of_range
923  * Node index is out of range.
924  *
925  * Returns the intensity of node with specified @p index in the file function.
926  ***************************************************************************/
927 double GModelSpectralFunc::intensity(const int& index) const
928 {
929  // Compile option: raise exception if index is out of range
930  #if defined(G_RANGE_CHECK)
931  if (index < 0 || index >= nodes()) {
932  throw GException::out_of_range(G_INTENSITY1, "Node index", index, nodes());
933  }
934  #endif
935 
936  // Return intensity
937  return (m_lin_values[index]);
938 }
939 
940 
941 /***********************************************************************//**
942  * @brief Set intensity of node
943  *
944  * @param[in] index Node index [0,...,nodes()-1].
945  * @param[in] intensity Intensity (ph/cm2/s/MeV).
946  *
947  * @exception GException::out_of_range
948  * Node index is out of range.
949  * @exception GException::invalid_argument
950  * Intensity is not positive.
951  *
952  * Sets the intensity of node with specified @p index in the file function.
953  ***************************************************************************/
954 void GModelSpectralFunc::intensity(const int& index, const double& intensity)
955 {
956  // Compile option: raise exception if index is out of range
957  #if defined(G_RANGE_CHECK)
958  if (index < 0 || index >= nodes()) {
959  throw GException::out_of_range(G_INTENSITY2, "Node index", index, nodes());
960  }
961  #endif
962 
963  // Check that intensity is positive
964  if (intensity <= 0.0) {
965  std::string msg = "Specified intensity "+gammalib::str(intensity)+" is "
966  "not positive. Please set only positive intensities.";
968  }
969 
970  // Set intensity
971  m_lin_values[index] = intensity;
972  m_log_values[index] = std::log10(intensity);
973 
974  // Set pre-computation cache
975  set_cache();
976 
977  // Return
978  return;
979 }
980 
981 
982 /***********************************************************************//**
983  * @brief Save file function in ASCII file
984  *
985  * @param[in] filename ASCII filename.
986  * @param[in] clobber Overwrite existing file?
987  *
988  * Saves file function in ASCII file.
989  ***************************************************************************/
990 void GModelSpectralFunc::save(const GFilename& filename, const bool& clobber) const
991 {
992  // Allocate CSV file
993  GCsv csv(nodes(), 2);
994 
995  // Fill CSV file
996  for (int i = 0; i < nodes(); ++i) {
997  csv(i,0) = gammalib::str(m_lin_nodes[i]);
998  csv(i,1) = gammalib::str(m_lin_values[i]);
999  }
1000 
1001  // Save CSV file
1002  csv.save(filename, " ", clobber);
1003 
1004  // Store filename
1005  m_filename = filename;
1006 
1007  // Return
1008  return;
1009 }
1010 
1011 
1012 /***********************************************************************//**
1013  * @brief Print file function information
1014  *
1015  * @param[in] chatter Chattiness.
1016  * @return String containing file function information.
1017  ***************************************************************************/
1018 std::string GModelSpectralFunc::print(const GChatter& chatter) const
1019 {
1020  // Initialise result string
1021  std::string result;
1022 
1023  // Continue only if chatter is not silent
1024  if (chatter != SILENT) {
1025 
1026  // Append header
1027  result.append("=== GModelSpectralFunc ===");
1028 
1029  // Append information
1030  result.append("\n"+gammalib::parformat("Function file"));
1031  result.append(m_filename.url());
1032  result.append("\n"+gammalib::parformat("Number of nodes"));
1033  result.append(gammalib::str(nodes()));
1034  result.append("\n"+gammalib::parformat("Number of parameters"));
1035  result.append(gammalib::str(size()));
1036  for (int i = 0; i < size(); ++i) {
1037  result.append("\n"+m_pars[i]->print(chatter));
1038  }
1039 
1040  // Append node information
1041  for (int i = 0; i < m_prefactor.size(); ++i) {
1042  result.append("\n"+gammalib::parformat("Node "+gammalib::str(i+1)));
1043  result.append("Epivot="+gammalib::str(m_epivot[i]));
1044  result.append(" Prefactor="+gammalib::str(m_prefactor[i]));
1045  result.append(" Gamma="+gammalib::str(m_gamma[i]));
1046  result.append(" Flux="+gammalib::str(m_flux[i]));
1047  result.append(" EFlux="+gammalib::str(m_eflux[i]));
1048  }
1049 
1050  } // endif: chatter was not silent
1051 
1052  // Return result
1053  return result;
1054 }
1055 
1056 
1057 /*==========================================================================
1058  = =
1059  = Private methods =
1060  = =
1061  ==========================================================================*/
1062 
1063 /***********************************************************************//**
1064  * @brief Initialise class members
1065  ***************************************************************************/
1067 {
1068  // Initialise powerlaw normalisation
1069  m_norm.clear();
1070  m_norm.name("Normalization");
1071  m_norm.scale(1.0);
1072  m_norm.value(1.0);
1073  m_norm.range(0.0,1000.0);
1074  m_norm.free();
1075  m_norm.gradient(0.0);
1076  m_norm.has_grad(true);
1077 
1078  // Set parameter pointer(s)
1079  m_pars.clear();
1080  m_pars.push_back(&m_norm);
1081 
1082  // Initialise members
1083  m_lin_nodes.clear();
1084  m_log_nodes.clear();
1085  m_lin_values.clear();
1086  m_log_values.clear();
1087  m_filename.clear();
1088 
1089  // Initialise flux cache
1090  m_prefactor.clear();
1091  m_gamma.clear();
1092  m_epivot.clear();
1093  m_flux.clear();
1094  m_eflux.clear();
1095 
1096  // Initialise MC cache
1097  m_mc_emin.clear();
1098  m_mc_emax.clear();
1099  m_mc_cum.clear();
1100  m_mc_min.clear();
1101  m_mc_max.clear();
1102  m_mc_exp.clear();
1103 
1104  // Return
1105  return;
1106 }
1107 
1108 
1109 /***********************************************************************//**
1110  * @brief Copy class members
1111  *
1112  * @param[in] model Spectral function model.
1113  ***************************************************************************/
1115 {
1116  // Copy members
1117  m_norm = model.m_norm;
1118  m_lin_nodes = model.m_lin_nodes;
1119  m_log_nodes = model.m_log_nodes;
1120  m_lin_values = model.m_lin_values;
1121  m_log_values = model.m_log_values;
1122  m_filename = model.m_filename;
1123 
1124  // Copy flux cache
1125  m_prefactor = model.m_prefactor;
1126  m_gamma = model.m_gamma;
1127  m_epivot = model.m_epivot;
1128  m_flux = model.m_flux;
1129  m_eflux = model.m_eflux;
1130 
1131  // Copy MC cache
1132  m_mc_emin = model.m_mc_emin;
1133  m_mc_emax = model.m_mc_emax;
1134  m_mc_cum = model.m_mc_cum;
1135  m_mc_min = model.m_mc_min;
1136  m_mc_max = model.m_mc_max;
1137  m_mc_exp = model.m_mc_exp;
1138 
1139  // Set parameter pointer(s)
1140  m_pars.clear();
1141  m_pars.push_back(&m_norm);
1142 
1143  // Return
1144  return;
1145 }
1146 
1147 
1148 /***********************************************************************//**
1149  * @brief Delete class members
1150  ***************************************************************************/
1152 {
1153  // Return
1154  return;
1155 }
1156 
1157 
1158 /***********************************************************************//**
1159  * @brief Load nodes from file
1160  *
1161  * @param[in] filename File name.
1162  *
1163  * @exception GException::invalid_argument
1164  * Empty @p filename specified.
1165  * GException::invalid_value
1166  * Invalid file function ASCII file.
1167  *
1168  * The file function is stored as a column separated value table (CSV) in an
1169  * ASCII file with (at least) 2 columns. The first column specifies the
1170  * energy in MeV while the second column specifies the intensity at this
1171  * energy in units of ph/cm2/s/MeV.
1172  *
1173  * The node energies and values will be stored both linearly and as log10.
1174  * The log10 storing requires that node energies and node values are
1175  * positive. Also, at least 2 nodes and 2 columns are required in the file
1176  * function.
1177  ***************************************************************************/
1179 {
1180  // Clear nodes and values
1181  m_lin_nodes.clear();
1182  m_log_nodes.clear();
1183  m_lin_values.clear();
1184  m_log_values.clear();
1185 
1186  // Throw an exception if the filename is empty
1187  if (filename.is_empty()) {
1188  std::string msg = "Empty file function ASCII file name specified. "
1189  "Please specify a valid file function ASCII file "
1190  "name.";
1192  }
1193 
1194  // Set filename
1195  m_filename = filename;
1196 
1197  // Load file
1198  GCsv csv = GCsv(filename.url());
1199 
1200  // Check if there are at least 2 nodes
1201  if (csv.nrows() < 2) {
1202  std::string msg = "File function ASCII file \""+filename.url()+
1203  "\" contains "+gammalib::str(csv.nrows())+
1204  " rows but at least 2 rows are required. Please "
1205  "specify a file function ASCII file with at "
1206  "least two rows.";
1208  }
1209 
1210  // Check if there are at least 2 columns
1211  if (csv.ncols() < 2) {
1212  std::string msg = "File function ASCII file \""+filename.url()+
1213  "\" contains "+gammalib::str(csv.ncols())+
1214  " columns but at least 2 columns are required. "
1215  "Please specify a file function ASCII file with "
1216  "at least two columns.";
1218  }
1219 
1220  // Setup nodes
1221  double last_energy = 0.0;
1222  for (int i = 0; i < csv.nrows(); ++i) {
1223 
1224  // Get log10 of node energy and value. Make sure they are valid.
1225  double log10energy;
1226  double log10value;
1227  if (csv.real(i,0) > 0.0) {
1228  log10energy = std::log10(csv.real(i,0));
1229  }
1230  else {
1231  std::string msg = "Non-positive energy "+
1232  gammalib::str(csv.real(i,0))+" encountered "
1233  "in row "+gammalib::str(i)+" of file "
1234  "function ASCII file. Please specify only "
1235  "positive energies in file function.";
1237  }
1238  if (csv.real(i,1) > 0.0) {
1239  log10value = std::log10(csv.real(i,1));
1240  }
1241  else {
1242  std::string msg = "Non-positive intensity "+
1243  gammalib::str(csv.real(i,1))+" encountered "
1244  "in row "+gammalib::str(i)+" of file "
1245  "function ASCII file. Please specify only "
1246  "positive intensities in file function.";
1248  }
1249 
1250  // Make sure that energies are increasing
1251  if (csv.real(i,0) <= last_energy) {
1252  std::string msg = "Energy "+gammalib::str(csv.real(i,0))+
1253  "in row "+gammalib::str(i)+" of file "
1254  "function ASCII file is equal or smaller "
1255  "than preceding energy "+
1256  gammalib::str(last_energy)+". Please specify "
1257  "monotonically increasing energies in the "
1258  "file function ASCII file.";
1260  }
1261 
1262  // Append log10 of node energy and value
1263  m_lin_nodes.append(csv.real(i,0));
1264  m_log_nodes.append(log10energy);
1265  m_lin_values.push_back(csv.real(i,1));
1266  m_log_values.push_back(log10value);
1267 
1268  // Store last energy for monotonically increasing check
1269  last_energy = csv.real(i,0);
1270 
1271  } // endfor: looped over nodes
1272 
1273  // Set pre-computation cache
1274  set_cache();
1275 
1276  // Return
1277  return;
1278 }
1279 
1280 
1281 /***********************************************************************//**
1282  * @brief Set pre-computation cache
1283  ***************************************************************************/
1285 {
1286  // Clear any existing values
1287  m_prefactor.clear();
1288  m_gamma.clear();
1289  m_epivot.clear();
1290  m_flux.clear();
1291  m_eflux.clear();
1292 
1293  // Loop over all nodes-1
1294  for (int i = 0; i < nodes()-1; ++i) {
1295 
1296  // Get energies and function values
1297  double emin = m_lin_nodes[i];
1298  double emax = m_lin_nodes[i+1];
1299  double fmin = m_lin_values[i];
1300  double fmax = m_lin_values[i+1];
1301 
1302  // Compute pivot energy (MeV). We use here the geometric mean of
1303  // the node boundaries.
1304  double epivot = std::sqrt(emin*emax);
1305 
1306  // Compute spectral index
1307  double gamma = std::log(fmin/fmax) / std::log(emin/emax);
1308 
1309  // Compute power law normalisation
1310  double prefactor = fmin / std::pow(emin/epivot, gamma);
1311 
1312  // Compute photon flux between nodes
1313  double flux = prefactor *
1314  gammalib::plaw_photon_flux(emin, emax, epivot, gamma);
1315 
1316  // Compute energy flux between nodes
1317  double eflux = prefactor *
1318  gammalib::plaw_energy_flux(emin, emax, epivot, gamma);
1319 
1320  // Convert energy flux from MeV/cm2/s to erg/cm2/s
1321  eflux *= gammalib::MeV2erg;
1322 
1323  // Push back values on pre-computation cache
1324  m_prefactor.push_back(prefactor);
1325  m_gamma.push_back(gamma);
1326  m_epivot.push_back(epivot);
1327  m_flux.push_back(flux);
1328  m_eflux.push_back(eflux);
1329 
1330  } // endfor: looped over all nodes
1331 
1332  // Return
1333  return;
1334 }
1335 
1336 
1337 /***********************************************************************//**
1338  * @brief Set MC pre-computation cache
1339  *
1340  * @param[in] emin Minimum energy.
1341  * @param[in] emax Maximum energy.
1342  *
1343  * This method sets up an array of indices and the cumulative distribution
1344  * function needed for MC simulations.
1345  ***************************************************************************/
1347  const GEnergy& emax) const
1348 {
1349  // Check if we need to update the cache
1350  if (emin != m_mc_emin || emax != m_mc_emax) {
1351 
1352  // Store new energy interval
1353  m_mc_emin = emin;
1354  m_mc_emax = emax;
1355 
1356  // Initialise cache
1357  m_mc_cum.clear();
1358  m_mc_min.clear();
1359  m_mc_max.clear();
1360  m_mc_exp.clear();
1361 
1362  // Get energy range in MeV
1363  double e_min = emin.MeV();
1364  double e_max = emax.MeV();
1365 
1366  // Continue only if e_max > e_min
1367  if (e_max > e_min) {
1368 
1369  // Allocate flux
1370  double flux;
1371 
1372  // Determine left node index for minimum energy
1373  m_lin_nodes.set_value(e_min);
1374  int inx_emin = m_lin_nodes.inx_left();
1375 
1376  // Determine left node index for maximum energy
1377  m_lin_nodes.set_value(e_max);
1378  int inx_emax = m_lin_nodes.inx_left();
1379 
1380  // If both energies are within the same node then just
1381  // add this one node on the stack
1382  if (inx_emin == inx_emax) {
1383  flux = m_prefactor[inx_emin] *
1385  e_max,
1386  m_epivot[inx_emin],
1387  m_gamma[inx_emin]);
1388  m_mc_cum.push_back(flux);
1389  m_mc_min.push_back(e_min);
1390  m_mc_max.push_back(e_max);
1391  m_mc_exp.push_back(m_gamma[inx_emin]);
1392  }
1393 
1394  // ... otherwise integrate over the nodes where emin and emax
1395  // resides and all the remaining nodes
1396  else {
1397 
1398  // If we are requested to extrapolate beyond first node,
1399  // the use the first nodes lower energy and upper energy
1400  // boundary for the initial integration.
1401  int i_start = (e_min < m_lin_nodes[0]) ? inx_emin : inx_emin+1;
1402 
1403  // Add emin to the node boundary
1404  flux = m_prefactor[inx_emin] *
1406  m_lin_nodes[i_start],
1407  m_epivot[inx_emin],
1408  m_gamma[inx_emin]);
1409  m_mc_cum.push_back(flux);
1410  m_mc_min.push_back(e_min);
1411  m_mc_max.push_back(m_lin_nodes[i_start]);
1412  m_mc_exp.push_back(m_gamma[inx_emin]);
1413 
1414  // Add all nodes between
1415  for (int i = i_start; i < inx_emax; ++i) {
1416  flux = m_flux[i];
1417  m_mc_cum.push_back(flux);
1418  m_mc_min.push_back(m_lin_nodes[i]);
1419  m_mc_max.push_back(m_lin_nodes[i+1]);
1420  m_mc_exp.push_back(m_gamma[i]);
1421  }
1422 
1423  // Add node boundary to emax
1424  flux = m_prefactor[inx_emax] *
1426  e_max,
1427  m_epivot[inx_emax],
1428  m_gamma[inx_emax]);
1429  m_mc_cum.push_back(flux);
1430  m_mc_min.push_back(m_lin_nodes[inx_emax]);
1431  m_mc_max.push_back(e_max);
1432  m_mc_exp.push_back(m_gamma[inx_emax]);
1433 
1434  } // endelse: emin and emax not between same nodes
1435 
1436  // Build cumulative distribution
1437  for (int i = 1; i < m_mc_cum.size(); ++i) {
1438  m_mc_cum[i] += m_mc_cum[i-1];
1439  }
1440  double norm = m_mc_cum[m_mc_cum.size()-1];
1441  if (norm > 0.0) {
1442  for (int i = 0; i < m_mc_cum.size(); ++i) {
1443  m_mc_cum[i] /= norm;
1444  }
1445  }
1446 
1447  // Set MC values
1448  for (int i = 0; i < m_mc_cum.size(); ++i) {
1449 
1450  // Compute exponent
1451  double exponent = m_mc_exp[i] + 1.0;
1452 
1453  // Exponent dependend computation
1454  if (std::abs(exponent) > 1.0e-11) {
1455 
1456  // If the exponent is too small then use lower energy
1457  // boundary
1458  if (exponent < -50.0) {
1459  m_mc_exp[i] = 0.0;
1460  m_mc_min[i] = std::log(m_mc_min[i]);
1461  m_mc_max[i] = m_mc_min[i];
1462  }
1463 
1464  // ... otherwise if exponent is too large then use
1465  // upper energy boundary
1466  else if (exponent > +50.0) {
1467  m_mc_exp[i] = 0.0;
1468  m_mc_min[i] = std::log(m_mc_max[i]);
1469  m_mc_max[i] = m_mc_min[i];
1470  }
1471 
1472  // ... otherwise use transformation formula
1473  else {
1474  m_mc_exp[i] = exponent;
1475  m_mc_min[i] = std::pow(m_mc_min[i], exponent);
1476  m_mc_max[i] = std::pow(m_mc_max[i], exponent);
1477  }
1478  }
1479  else {
1480  m_mc_exp[i] = 0.0;
1481  m_mc_min[i] = std::log(m_mc_min[i]);
1482  m_mc_max[i] = std::log(m_mc_max[i]);
1483  }
1484 
1485  } // endfor: set MC values
1486 
1487  } // endif: e_max > e_min
1488 
1489  } // endif: Update was required
1490 
1491  // Return
1492  return;
1493 }
std::vector< double > m_gamma
Power-law indices.
void mc_update(const GEnergy &emin, const GEnergy &emax) const
Set MC pre-computation cache.
int nodes(void) const
Return number of nodes in file function.
GEnergy m_mc_emax
Maximum energy.
void remove(const int &index)
Remove node from file function.
void save(const GFilename &filename, const std::string &sep=" ", const bool &clobber=false) const
Save CSV table.
Definition: GCsv.cpp:519
virtual double flux(const GEnergy &emin, const GEnergy &emax) const
Returns model photon flux between emin, emax
std::vector< double > m_lin_values
Function values at nodes.
const double & factor_gradient(void) const
Return parameter factor gradient.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:864
#define G_REMOVE
const std::string & name(void) const
Return parameter name.
Random number generator class definition.
Abstract spectral model base class.
bool is_empty(void) const
Signals if there are nodes in file function.
double gradient(void) const
Return parameter gradient.
#define G_WRITE
std::vector< double > m_mc_max
Upper boundary for MC.
int size(void) const
Return number of parameters.
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1253
virtual GModelSpectral & operator=(const GModelSpectral &model)
Assignment operator.
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:351
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const =0
#define G_ENERGY2
std::vector< GModelPar * > m_pars
Parameter pointers.
bool is_empty(void) const
Signal if filename is empty.
Definition: GFilename.hpp:160
std::vector< double > m_mc_exp
Exponent for MC.
XML element node class.
Definition: GXmlElement.hpp:48
GEnergy energy(const int &index) const
Return energy of node.
virtual void read(const GXmlElement &xml)
Read model from XML element.
Spectral model registry class definition.
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
Random number generator class.
Definition: GRan.hpp:44
Comma-separated values table class.
Definition: GCsv.hpp:57
#define G_MC
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
Time class.
Definition: GTime.hpp:55
void load_nodes(const GFilename &filename)
Load nodes from file.
Gammalib tools definition.
void free_members(void)
Delete class members.
void set_value(const double &value) const
Set indices and weighting factors for interpolation.
Definition: GNodeArray.cpp:587
GFilename xml_file_reduce(const GXmlElement &xml, const std::string &filename)
Reduce file name provided for writing as XML attribute.
Definition: GTools.cpp:1946
const GFilename & filename(void) const
Return file name.
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:1248
bool is_free(void) const
Signal if parameter is free.
Energy container class.
Definition: GEnergies.hpp:60
std::vector< double > m_log_values
log10(Function) values at nodes
virtual void write(GXmlElement &xml) const
Write model into XML element.
double log10MeV(void) const
Return log10 of energy in MeV.
Definition: GEnergy.cpp:423
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:201
const double & scale(void) const
Return parameter scale.
GFilename m_filename
Name of file function.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:184
virtual ~GModelSpectralFunc(void)
Destructor.
const int & ncols(void) const
Return number of columns.
Definition: GCsv.hpp:137
std::vector< double > m_eflux
Energy fluxes.
void insert(const GEnergy &energy, const double &intensity)
Insert node into file function.
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1358
void remove(const int &index)
Remove one node into array.
Definition: GNodeArray.cpp:413
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
#define G_INSERT
const GXmlAttribute * attribute(const int &index) const
Return attribute.
void save(const GFilename &filename, const bool &clobber=false) const
Save file function in ASCII file.
void reserve(const int &num)
Reserves space for nodes in file function.
std::string print(const GChatter &chatter=NORMAL) const
Print energy.
Definition: GEnergy.cpp:748
std::vector< double > m_mc_min
Lower boundary for MC.
GEnergy m_mc_emin
Minimum energy.
void free(void)
Free a parameter.
virtual void clear(void)
Clear file function.
Filename class.
Definition: GFilename.hpp:62
Energy container class definition.
void set_cache(void) const
Set pre-computation cache.
GNodeArray m_log_nodes
log10(Energy) nodes of function
Spectral function model class definition.
int size(void) const
Return number of energies in container.
Definition: GEnergies.hpp:170
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:1637
GModelPar m_norm
Normalization factor.
#define G_INTENSITY2
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1274
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
void check_energy_interval(const std::string &origin, const GEnergy &emin, const GEnergy &emax)
Checks energy interval.
Definition: GException.cpp:364
void clear(void)
Clear parameter.
virtual double eflux(const GEnergy &emin, const GEnergy &emax) const
Returns model energy flux between emin, emax
double interpolate(const double &value, const std::vector< double > &vector) const
Interpolate value.
Definition: GNodeArray.cpp:536
#define G_LOAD_NODES
GChatter
Definition: GTypemaps.hpp:33
const int & inx_left(void) const
Returns left node index.
Definition: GNodeArray.hpp:235
#define G_ENERGY1
const int & nrows(void) const
Return number of rows.
Definition: GCsv.hpp:149
void copy_members(const GModelSpectralFunc &model)
Copy class members.
void xml_check_parnum(const std::string &origin, const GXmlElement &xml, const int &number)
Checks number of parameters.
Definition: GTools.cpp:1777
Interface definition for the spectral model registry class.
virtual std::string type(void) const
Return model type.
std::vector< double > m_epivot
Power-law pivot energies.
double intensity(const int &index) const
Return intensity of node.
double real(const int &row, const int &col) const
Get double precision value.
Definition: GCsv.cpp:324
Spectral function model class.
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition: GFilename.hpp:189
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
void init_members(void)
Initialise class members.
virtual GEnergy mc(const GEnergy &emin, const GEnergy &emax, const GTime &time, GRan &ran) const
Returns MC energy between [emin, emax].
void insert(const int &index, const double &node)
Insert one node into array.
Definition: GNodeArray.cpp:376
const double MeV2erg
Definition: GTools.hpp:45
std::vector< double > m_mc_cum
Cumulative distribution.
void range(const double &min, const double &max)
Set minimum and maximum parameter boundaries.
void read(const GXmlElement &xml)
Extract parameter attributes from XML element.
Definition: GModelPar.cpp:229
virtual GModelSpectralFunc * clone(void) const
Clone file function.
double value(void) const
Return parameter value.
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1422
std::vector< double > m_prefactor
Power-law normalisations.
Exception handler interface definition.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print file function information.
void init_members(void)
Initialise class members.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1232
#define G_READ
GModelSpectralFunc(void)
Void constructor.
GNodeArray m_lin_nodes
Energy nodes of function.
#define G_APPEND
void extend(const GModelSpectralFunc &filefct)
Append file function.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
void append(const double &node)
Append one node to array.
Definition: GNodeArray.cpp:351
const GModelSpectralFunc g_spectral_func_seed
#define G_INTENSITY1
void append(const GEnergy &energy, const double &intensity)
Append node to file function.
virtual GModelSpectralFunc & operator=(const GModelSpectralFunc &model)
Assignment operator.
double norm(void) const
Return normalization factor.
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:1689
void free_members(void)
Delete class members.
void clear(void)
Clear instance.
Definition: GEnergy.cpp:261
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const
Evaluate function.
Comma-separated values table class definition.
GVector log10(const GVector &vector)
Computes base10 logarithm of vector elements.
Definition: GVector.cpp:1295
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:1200
GFilename xml_file_expand(const GXmlElement &xml, const std::string &filename)
Expand file name provided as XML attribute for loading.
Definition: GTools.cpp:1889
std::vector< double > m_flux
Photon fluxes.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489
void xml_check_type(const std::string &origin, GXmlElement &xml, const std::string &type)
Checks the model type.
Definition: GTools.cpp:1819