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