GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpectralPlawEnergyFlux.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpectralPlawEnergyFlux.cpp - Spectral power law model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2016 by Michael Mayer *
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 GModelSpectralPlawEnergyFlux.cpp
23  * @brief Energy flux normalized power law spectral model class implementation
24  * @author Michael Mayer
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 "GMath.hpp"
35 #include "GRan.hpp"
38 
39 /* __ Constants __________________________________________________________ */
40 
41 /* __ Globals ____________________________________________________________ */
43  "EnergyFlux",
44  "Index",
45  "LowerLimit",
46  "UpperLimit");
47 const GModelSpectralRegistry g_spectral_plaw_eflux_registry(&g_spectral_plaw_eflux_seed);
48 
49 /* __ Method name definitions ____________________________________________ */
50 #define G_FLUX "GModelSpectralPlawEnergyFlux::flux(GEnergy&, GEnergy&)"
51 #define G_EFLUX "GModelSpectralPlawEnergyFlux::eflux(GEnergy&, GEnergy&)"
52 #define G_MC "GModelSpectralPlawEnergyFlux::mc(GEnergy&, GEnergy&, GTime&, "\
53  "GRan&)"
54 #define G_READ "GModelSpectralPlawEnergyFlux::read(GXmlElement&)"
55 #define G_WRITE "GModelSpectralPlawEnergyFlux::write(GXmlElement&)"
56 
57 /* __ Macros _____________________________________________________________ */
58 
59 /* __ Coding definitions _________________________________________________ */
60 
61 /* __ Debug definitions __________________________________________________ */
62 
63 
64 /*==========================================================================
65  = =
66  = Constructors/destructors =
67  = =
68  ==========================================================================*/
69 
70 /***********************************************************************//**
71  * @brief Void constructor
72  ***************************************************************************/
75 {
76  // Initialise members
77  init_members();
78 
79  // Return
80  return;
81 }
82 
83 
84 /***********************************************************************//**
85  * @brief Model type and parameter name constructor
86  *
87  * @param[in] type Model type.
88  * @param[in] eflux Name of energy flux parameter.
89  * @param[in] index Name of index parameter.
90  * @param[in] emin Name of emin parameter.
91  * @param[in] emax Name of emax parameter.
92  ***************************************************************************/
94  const std::string& eflux,
95  const std::string& index,
96  const std::string& emin,
97  const std::string& emax) :
99 {
100  // Initialise members
101  init_members();
102 
103  // Set model type
104  m_type = type;
105 
106  // Set parameter names
107  m_eflux.name(eflux);
108  m_index.name(index);
109  m_emin.name(emin);
110  m_emax.name(emax);
111 
112  // Return
113  return;
114 }
115 
116 
117 /***********************************************************************//**
118  * @brief Constructor
119  *
120  * @param[in] eflux Energy flux (erg/cm2/s).
121  * @param[in] index Power law index.
122  * @param[in] emin Minimum energy.
123  * @param[in] emax Maximum energy.
124  *
125  * Construct a spectral power law from the
126  * - energy flux (in erg/cm2/s),
127  * - spectral index,
128  * - minimum energy and
129  * - maximum energy.
130  ***************************************************************************/
132  const double& index,
133  const GEnergy& emin,
134  const GEnergy& emax) :
136 {
137  // Initialise members
138  init_members();
139 
140  // Set parameters
141  m_eflux.value(eflux);
142  m_index.value(index);
143  m_emin.value(emin.MeV());
144  m_emax.value(emax.MeV());
145 
146  // Return
147  return;
148 }
149 
150 
151 /***********************************************************************//**
152  * @brief XML constructor
153  *
154  * @param[in] xml XML element.
155  *
156  * Constructs energy flux normalized power law spectral model by extracting
157  * information from an XML element.
158  ***************************************************************************/
161 {
162  // Initialise members
163  init_members();
164 
165  // Read information from XML element
166  read(xml);
167 
168  // Return
169  return;
170 }
171 
172 
173 /***********************************************************************//**
174  * @brief Copy constructor
175  *
176  * @param[in] model Spectral power law model.
177  ***************************************************************************/
179  GModelSpectral(model)
180 {
181  // Initialise members
182  init_members();
183 
184  // Copy members
185  copy_members(model);
186 
187  // Return
188  return;
189 }
190 
191 
192 /***********************************************************************//**
193  * @brief Destructor
194  ***************************************************************************/
196 {
197  // Free members
198  free_members();
199 
200  // Return
201  return;
202 }
203 
204 
205 /*==========================================================================
206  = =
207  = Operators =
208  = =
209  ==========================================================================*/
210 
211 /***********************************************************************//**
212  * @brief Assignment operator
213  *
214  * @param[in] model Spectral power law model.
215  * @return Spectral power law model.
216  ***************************************************************************/
218 {
219  // Execute only if object is not identical
220  if (this != &model) {
221 
222  // Copy base class members
223  this->GModelSpectral::operator=(model);
224 
225  // Free members
226  free_members();
227 
228  // Initialise members
229  init_members();
230 
231  // Copy members
232  copy_members(model);
233 
234  } // endif: object was not identical
235 
236  // Return
237  return *this;
238 }
239 
240 
241 /*==========================================================================
242  = =
243  = Public methods =
244  = =
245  ==========================================================================*/
246 
247 /***********************************************************************//**
248  * @brief Clear power law model
249  ***************************************************************************/
251 {
252  // Free class members (base and derived classes, derived class first)
253  free_members();
255 
256  // Initialise members
258  init_members();
259 
260  // Return
261  return;
262 }
263 
264 
265 /***********************************************************************//**
266  * @brief Clone power law model
267  *
268  * @return Pointer to deep copy of power law model
269  ***************************************************************************/
271 {
272  // Clone power law model
273  return new GModelSpectralPlawEnergyFlux(*this);
274 }
275 
276 
277 /***********************************************************************//**
278  * @brief Evaluate function
279  *
280  * @param[in] srcEng True photon energy.
281  * @param[in] srcTime True photon arrival time.
282  * @param[in] gradients Compute gradients?
283  * @return Model value (ph/cm2/s/MeV).
284  *
285  * Computes
286  *
287  * \f[
288  * S_{\rm E}(E | t) = {\tt m\_eflux}
289  * \frac{{\tt m\_index}+2}
290  * {{\tt e\_max}^{{\tt m\_index}+2} -
291  * {\tt e\_min}^{{\tt m\_index}+2}}
292  * E^{\tt m\_index}
293  * \f]
294  *
295  * for \f${\tt m\_index} \ne -2\f$ and
296  *
297  * \f[
298  * S_{\rm E}(E | t) =
299  * \frac{{\tt m\_eflux}}
300  * {\log {\tt e\_max} - \log {\tt e\_min}}
301  * E^{\tt m\_index}
302  * \f]
303  *
304  * for \f${\tt m\_index} = -2\f$, where
305  * - \f${\tt e\_min}\f$ is the minimum energy of an interval,
306  * - \f${\tt e\_max}\f$ is the maximum energy of an interval,
307  * - \f${\tt m\_eflux}\f$ is the energy flux between
308  * \f${\tt e\_min}\f$ and \f${\tt e\_max}\f$, and
309  * - \f${\tt m\_index}\f$ is the spectral index.
310  *
311  * If the @p gradients flag is true the method will also compute the
312  * partial derivatives of the model with respect to the parameters using
313  *
314  * \f[
315  * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_eflux}} =
316  * \frac{S_{\rm E}(E | t)}{{\tt m\_eflux}}
317  * \f]
318  *
319  * \f[
320  * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_index}} =
321  * S_{\rm E}(E | t) \,
322  * \left( \frac{1}{{\tt m\_index}+2} -
323  * \frac{\log({\tt e\_max}) {\tt e\_max}^{{\tt m\_index}+2} -
324  * \log({\tt e\_min}) {\tt e\_min}^{{\tt m\_index}+2}}
325  * {{\tt e\_max}^{{\tt m\_index}+2} -
326  * {\tt e\_min}^{{\tt m\_index}+2}}
327  * + \ln(E) \right)
328  * \f]
329  *
330  * for \f${\tt m\_index} \ne -2\f$ and
331  *
332  * \f[
333  * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_index}} =
334  * S_{\rm E}(E | t) \, \ln(E)
335  * \f]
336  *
337  * for \f${\tt m\_index} = -2\f$.
338  *
339  * No partial derivatives are supported for the energy boundaries.
340  ***************************************************************************/
342  const GTime& srcTime,
343  const bool& gradients) const
344 {
345  // Update precomputed values
346  update(srcEng);
347 
348  // Compute function value
349  double value = m_eflux.value() * gammalib::erg2MeV * m_norm * m_power;
350 
351  // Optionally compute gradients
352  if (gradients) {
353 
354  // Compute partial derivatives of the parameter values
355  double g_eflux = (m_eflux.is_free() && m_eflux.factor_value() > 0.0)
356  ? g_eflux = value / m_eflux.factor_value()
357  : 0.0;
358  double g_index = (m_index.is_free())
359  ? value * (m_g_norm + gammalib::ln10 *
360  srcEng.log10MeV()) * m_index.scale()
361  : 0.0;
362 
363  // Set gradients
364  m_eflux.factor_gradient(g_eflux);
365  m_index.factor_gradient(g_index);
366 
367  } // endif: gradient computation was requested
368 
369  // Compile option: Check for NaN/Inf
370  #if defined(G_NAN_CHECK)
371  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
372  std::cout << "*** ERROR: GModelSpectralPlawEnergyFlux::eval";
373  std::cout << "(srcEng=" << srcEng;
374  std::cout << ", srcTime=" << srcTime << "):";
375  std::cout << " NaN/Inf encountered";
376  std::cout << " (value=" << value;
377  std::cout << ", eflux=" << eflux();
378  std::cout << ", m_norm=" << m_norm;
379  std::cout << ", m_power=" << m_power;
380  std::cout << ")" << std::endl;
381  }
382  #endif
383 
384  // Return
385  return value;
386 }
387 
388 
389 /***********************************************************************//**
390  * @brief Returns model photon flux between [emin, emax] (units: ph/cm2/s)
391  *
392  * @param[in] emin Minimum photon energy.
393  * @param[in] emax Minimum photon energy.
394  * @return Photon flux (ph/cm2/s).
395  *
396  * Computes
397  *
398  * \f[
399  * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) dE
400  * \f]
401  *
402  * where
403  * - [@p emin, @p emax] is an energy interval, and
404  * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
405  ***************************************************************************/
407  const GEnergy& emax) const
408 {
409  // Initialise flux
410  double flux = 0.0;
411 
412  // Compute only if integration range is valid
413  if (emin < emax) {
414 
415  // Compute power law normalization
416  double norm = 0.0;
417  if (index() != -2.0) {
418  double gamma = m_index.value() + 2.0;
419  double pow_ref_emin = std::pow(this->emin().MeV(), gamma);
420  double pow_ref_emax = std::pow(this->emax().MeV(), gamma);
421  norm = m_eflux.value() * gammalib::erg2MeV * gamma / (pow_ref_emax - pow_ref_emin);
422  }
423  else {
424  double log_ref_emin = std::log(this->emin().MeV());
425  double log_ref_emax = std::log(this->emax().MeV());
426  norm = m_eflux.value() * gammalib::erg2MeV / (log_ref_emax - log_ref_emin);
427  }
428 
429  // Compute photon flux
430  if (index() != -1.0) {
431  double gamma = m_index.value() + 1.0;
432  double pow_emin = std::pow(emin.MeV(), gamma);
433  double pow_emax = std::pow(emax.MeV(), gamma);
434  flux = norm / gamma * (pow_emax - pow_emin);
435  }
436 
437  // Case B: Index is -1
438  else {
439  double log_emin = std::log(emin.MeV());
440  double log_emax = std::log(emax.MeV());
441  flux = norm * (log_emax - log_emin);
442  }
443 
444  } // endif: integration range was valid
445 
446  // Return flux
447  return flux;
448 }
449 
450 
451 /***********************************************************************//**
452  * @brief Returns model energy flux between [emin, emax] (units: erg/cm2/s)
453  *
454  * @param[in] emin Minimum photon energy.
455  * @param[in] emax Minimum photon energy.
456  * @return Energy flux (erg/cm2/s).
457  *
458  * Computes
459  *
460  * \f[
461  * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) E \, dE
462  * \f]
463  *
464  * where
465  * - [@p emin, @p emax] is an energy interval, and
466  * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
467  ***************************************************************************/
469  const GEnergy& emax) const
470 {
471  // Initialise flux
472  double eflux = 0.0;
473 
474  // Compute only if integration range is valid
475  if (emin < emax) {
476 
477  // Compute power law normalization
478  double norm;
479  if (index() != -2.0) {
480  double gamma = m_index.value() + 2.0;
481  double pow_ref_emin = std::pow(this->emin().MeV(), gamma);
482  double pow_ref_emax = std::pow(this->emax().MeV(), gamma);
483  double pow_emin = std::pow(emin.MeV(), gamma);
484  double pow_emax = std::pow(emax.MeV(), gamma);
485  double factor = (pow_emax - pow_emin) /
486  (pow_ref_emax - pow_ref_emin);
487  eflux = m_eflux.value() * factor;
488  }
489  else {
490  double log_emin = std::log(emin.MeV());
491  double log_emax = std::log(emax.MeV());
492  double log_ref_emin = std::log(this->emin().MeV());
493  double log_ref_emax = std::log(this->emax().MeV());
494  double factor = (log_emax - log_emin) /
495  (log_ref_emax - log_ref_emin);
496  eflux = m_eflux.value() * factor;
497  }
498 
499  } // endif: integration range was valid
500 
501  // Return flux
502  return eflux;
503 }
504 
505 
506 /***********************************************************************//**
507  * @brief Returns MC energy between [emin, emax]
508  *
509  * @param[in] emin Minimum photon energy.
510  * @param[in] emax Maximum photon energy.
511  * @param[in] time True photon arrival time.
512  * @param[in,out] ran Random number generator.
513  * @return Energy.
514  *
515  * @exception GException::erange_invalid
516  * Energy range is invalid (emin < emax required).
517  *
518  * Returns Monte Carlo energy by randomly drawing from a power law.
519  ***************************************************************************/
521  const GEnergy& emax,
522  const GTime& time,
523  GRan& ran) const
524 {
525  // Throw an exception if energy range is invalid
526  if (emin >= emax) {
527  throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(),
528  "Minimum energy < maximum energy required.");
529  }
530 
531  // Allocate energy
532  GEnergy energy;
533 
534  // Case A: Index is not -1
535  if (index() != -1.0) {
536  double exponent = index() + 1.0;
537  double e_max = std::pow(emax.MeV(), exponent);
538  double e_min = std::pow(emin.MeV(), exponent);
539  double u = ran.uniform();
540  double eng = (u > 0.0)
541  ? std::exp(std::log(u * (e_max - e_min) + e_min) / exponent)
542  : 0.0;
543  energy.MeV(eng);
544  }
545 
546  // Case B: Index is -1
547  else {
548  double e_max = std::log(emax.MeV());
549  double e_min = std::log(emin.MeV());
550  double u = ran.uniform();
551  double eng = std::exp(u * (e_max - e_min) + e_min);
552  energy.MeV(eng);
553  }
554 
555  // Return energy
556  return energy;
557 }
558 
559 
560 /***********************************************************************//**
561  * @brief Read model from XML element
562  *
563  * @param[in] xml XML element.
564  *
565  * Reads the spectral information from an XML element.
566  ***************************************************************************/
568 {
569  // Get parameter pointers
574 
575  // Read parameters
576  m_eflux.read(*eflux);
577  m_index.read(*index);
578  m_emin.read(*emin);
579  m_emax.read(*emax);
580 
581  // Return
582  return;
583 }
584 
585 
586 /***********************************************************************//**
587  * @brief Write model into XML element
588  *
589  * @param[in] xml XML element.
590  *
591  * @exception GException::model_invalid_spectral
592  * Existing XML element is not of type "PowerLaw"
593  *
594  * Writes the spectral information into an XML element.
595  ***************************************************************************/
597 {
598  // Set model type
599  if (xml.attribute("type") == "") {
600  xml.attribute("type", type());
601  }
602 
603  // Verify model type
604  if (xml.attribute("type") != type()) {
606  "Spectral model is not of type \""+type()+"\".");
607  }
608 
609  // Get XML parameters
614 
615  // Write parameters
616  m_eflux.write(*eflux);
617  m_index.write(*index);
618  m_emin.write(*emin);
619  m_emax.write(*emax);
620 
621  // Return
622  return;
623 }
624 
625 
626 /***********************************************************************//**
627  * @brief Print power law information
628  *
629  * @param[in] chatter Chattiness.
630  * @return String containing power law information.
631  ***************************************************************************/
632 std::string GModelSpectralPlawEnergyFlux::print(const GChatter& chatter) const
633 {
634  // Initialise result string
635  std::string result;
636 
637  // Continue only if chatter is not silent
638  if (chatter != SILENT) {
639 
640  // Append header
641  result.append("=== GModelSpectralPlawEnergyFlux ===");
642 
643  // Append information
644  result.append("\n"+gammalib::parformat("Number of parameters"));
645  result.append(gammalib::str(size()));
646  for (int i = 0; i < size(); ++i) {
647  result.append("\n"+m_pars[i]->print(chatter));
648  }
649 
650  } // endif: chatter was not silent
651 
652  // Return result
653  return result;
654 }
655 
656 
657 /*==========================================================================
658  = =
659  = Private methods =
660  = =
661  ==========================================================================*/
662 
663 /***********************************************************************//**
664  * @brief Initialise class members
665  ***************************************************************************/
667 {
668  // Initialise model type
669  m_type = "PowerLaw";
670 
671  // Initialise energy flux
672  m_eflux.clear();
673  m_eflux.name("EnergyFlux");
674  m_eflux.unit("erg/cm2/s");
675  m_eflux.scale(1.0);
676  m_eflux.value(1.0); // default: 1.0
677  m_eflux.range(0.0, 10.0); // range: [0,10]
678  m_eflux.free();
679  m_eflux.gradient(0.0);
680  m_eflux.has_grad(true);
681 
682  // Initialise powerlaw index
683  m_index.clear();
684  m_index.name("Index");
685  m_index.scale(1.0);
686  m_index.value(-2.0); // default: -2.0
687  m_index.range(-10.0,+10.0); // range: [-10,+10]
688  m_index.free();
689  m_index.gradient(0.0);
690  m_index.has_grad(true);
691 
692  // Initialise lower limit
693  m_emin.clear();
694  m_emin.name("LowerLimit");
695  m_emin.unit("MeV");
696  m_emin.scale(1.0);
697  m_emin.value(100.0); // default: 100
698  m_emin.range(0.001, 1.0e15); // range: [0.001, 1e15]
699  m_emin.fix();
700  m_emin.gradient(0.0);
701  m_emin.has_grad(false);
702 
703  // Initialise upper limit
704  m_emax.clear();
705  m_emax.name("UpperLimit");
706  m_emax.unit("MeV");
707  m_emax.scale(1.0);
708  m_emax.value(500000.0); // default: 500000
709  m_emin.range(0.001, 1.0e15); // range: [0.001, 1e15]
710  m_emax.fix();
711  m_emax.gradient(0.0);
712  m_emax.has_grad(false);
713 
714  // Set parameter pointer(s)
715  m_pars.clear();
716  m_pars.push_back(&m_eflux);
717  m_pars.push_back(&m_index);
718  m_pars.push_back(&m_emin);
719  m_pars.push_back(&m_emax);
720 
721  // Initialise last parameters (for fast computation)
722  m_log_emin = 0.0;
723  m_log_emax = 0.0;
724  m_pow_emin = 0.0;
725  m_pow_emax = 0.0;
726  m_norm = 0.0;
727  m_power = 0.0;
728  m_last_index = 1000.0;
729  m_last_emin.MeV(0.0);
730  m_last_emax.MeV(0.0);
731  m_last_energy.MeV(0.0);
732 
733  // Return
734  return;
735 }
736 
737 
738 /***********************************************************************//**
739  * @brief Copy class members
740  *
741  * @param[in] model Spectral power law model.
742  ***************************************************************************/
744 {
745  // Copy members
746  m_type = model.m_type;
747  m_eflux = model.m_eflux;
748  m_index = model.m_index;
749  m_emin = model.m_emin;
750  m_emax = model.m_emax;
751 
752  // Set parameter pointer(s)
753  m_pars.clear();
754  m_pars.push_back(&m_eflux);
755  m_pars.push_back(&m_index);
756  m_pars.push_back(&m_emin);
757  m_pars.push_back(&m_emax);
758 
759  // Copy bookkeeping information
760  m_log_emin = model.m_log_emin;
761  m_log_emax = model.m_log_emax;
762  m_pow_emin = model.m_pow_emin;
763  m_pow_emax = model.m_pow_emax;
764  m_norm = model.m_norm;
765  m_power = model.m_power;
766  m_last_index = model.m_last_index;
767  m_last_emin = model.m_last_emin;
768  m_last_emax = model.m_last_emax;
770 
771  // Return
772  return;
773 }
774 
775 
776 /***********************************************************************//**
777  * @brief Delete class members
778  ***************************************************************************/
780 {
781  // Return
782  return;
783 }
784 
785 
786 /***********************************************************************//**
787  * @brief Update precomputed values
788  *
789  * @param[in] srcEng Source energy
790  ***************************************************************************/
792 {
793  // Compute index+1
794  double gamma = index() + 2.0;
795 
796  // Change in spectral index?
797  if (index() != m_last_index) {
798 
799  // Save actual spectral index
800  m_last_index = index();
801 
802  // Change in energy boundaries?
803  if (emin() != m_last_emin || emax() != m_last_emax) {
804  m_log_emin = std::log(emin().MeV());
805  m_last_emin = emin();
806  m_log_emax = std::log(emax().MeV());
807  m_last_emax = emax();
808  }
809 
810  // Compute normalization factors
811  if (gamma != 0.0) {
812  m_pow_emin = std::pow(emin().MeV(), gamma);
813  m_pow_emax = std::pow(emax().MeV(), gamma);
814  double d = m_pow_emax - m_pow_emin;
815  m_norm = gamma / d;
816  m_g_norm = 1.0/gamma -
817  (m_pow_emax*m_log_emax - m_pow_emin*m_log_emin)/d;
818  }
819  else {
820  m_norm = 1.0 / (m_log_emax - m_log_emin);
821  m_g_norm = 0.0;
822  }
823 
824  // Update power law factor
825  m_power = std::pow(srcEng.MeV(), index());
826  m_last_energy = srcEng;
827 
828  } // endif: change in spectral index
829 
830  // ... no change in spectral index
831  else {
832 
833  // Change in energy boundaries?
834  if (emin() != m_last_emin || emax() != m_last_emax) {
835 
836  // Update energy boundaries
837  m_log_emin = std::log(emin().MeV());
838  m_last_emin = emin();
839  m_log_emax = std::log(emax().MeV());
840  m_last_emax = emax();
841 
842  // Compute power law normalization
843  if (gamma != 0.0) {
844  m_pow_emin = std::pow(emin().MeV(), gamma);
845  m_pow_emax = std::pow(emax().MeV(), gamma);
846  double d = m_pow_emax - m_pow_emin;
847  m_norm = gamma / d;
848  m_g_norm = 1.0/gamma -
849  (m_pow_emax*m_log_emax - m_pow_emin*m_log_emin)/d;
850  }
851  else {
852  m_norm = 1.0 / (m_log_emax - m_log_emin);
853  m_g_norm = 0.0;
854  }
855 
856  } // endif: change in energy boundaries
857 
858  // Change in energy?
859  if (srcEng != m_last_energy) {
860  m_power = std::pow(srcEng.MeV(), index());
861  m_last_energy = srcEng;
862  }
863 
864  } // endelse: no change in spectral index
865 
866  // Return
867  return;
868 }
const double erg2MeV
Definition: GTools.hpp:46
virtual void write(GXmlElement &xml) const
Write model into XML element.
double m_last_index
Last spectral index (MeV)
const double & factor_gradient(void) const
Return parameter gradient factor.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:821
const std::string & name(void) const
Return parameter name.
Random number generator class definition.
Abstract spectral model base class.
double gradient(void) const
Return parameter gradient.
int size(void) const
Return number of parameters.
virtual GModelSpectral & operator=(const GModelSpectral &model)
Assignment operator.
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:353
std::vector< GModelPar * > m_pars
Parameter pointers.
XML element node class.
Definition: GXmlElement.hpp:47
Spectral model registry class definition.
GEnergy emin(void) const
Return minimum energy.
Random number generator class.
Definition: GRan.hpp:44
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
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].
Energy flux normalized power law spectral model class interface definition.
virtual ~GModelSpectralPlawEnergyFlux(void)
Destructor.
bool is_free(void) const
Signal if parameter is free.
virtual GModelSpectralPlawEnergyFlux & operator=(const GModelSpectralPlawEnergyFlux &model)
Assignment operator.
virtual GModelSpectralPlawEnergyFlux * clone(void) const
Clone power law model.
const double ln10
Definition: GMath.hpp:46
double log10MeV(void) const
Return log10 of energy in MeV.
Definition: GEnergy.cpp:423
GModelPar m_emin
Lower energy limit (MeV)
void update(const GEnergy &srcEng) const
Update precomputed values.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:184
const double & scale(void) const
Return parameter scale.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
GModelPar m_eflux
Energy flux (erg/cm2/s)
GEnergy m_last_emax
Last upper energy limit.
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
GEnergy m_last_emin
Last lower energy limit.
virtual void read(const GXmlElement &xml)
Read model from XML element.
const GXmlAttribute * attribute(const int &index) const
Return attribute.
virtual std::string type(void) const
Return model type.
double m_g_norm
Power-law normalization gradient.
void free(void)
Free a parameter.
void fix(void)
Fix a parameter.
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
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
void clear(void)
Clear parameter.
GModelPar m_emax
Upper energy limit (MeV)
GModelSpectralPlawEnergyFlux(void)
Void constructor.
GChatter
Definition: GTypemaps.hpp:33
virtual double flux(const GEnergy &emin, const GEnergy &emax) const
Returns model photon flux between emin, emax
Interface definition for the spectral model registry class.
GEnergy m_last_energy
Last source energy.
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const
Evaluate function.
void init_members(void)
Initialise class members.
void copy_members(const GModelSpectralPlawEnergyFlux &model)
Copy class members.
double m_norm
Power-law normalization (for pivot energy 1 MeV)
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
double value(void) const
Return parameter value.
double index(void) const
Return power law index.
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
void free_members(void)
Delete class members.
const GModelSpectralPlawEnergyFlux g_spectral_plaw_eflux_seed("PowerLaw","EnergyFlux","Index","LowerLimit","UpperLimit")
Exception handler interface definition.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
double eflux(void) const
Return energy flux.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
Energy flux normalized power law spectral model class.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print power law information.
GEnergy emax(void) const
Return maximum energy.
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.
virtual void clear(void)
Clear power law model.
const double & factor_value(void) const
Return parameter value factor.
Mathematical function definitions.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413