GammaLib  1.7.0.dev
 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-2016 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  * @exception GException::erange_invalid
527  * Energy range is invalid (emin < emax required).
528  *
529  * Returns Monte Carlo energy by randomly drawing from a power law.
530  ***************************************************************************/
532  const GEnergy& emax,
533  const GTime& time,
534  GRan& ran) const
535 {
536  // Throw an exception if energy range is invalid
537  if (emin >= emax) {
538  throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(),
539  "Minimum energy < maximum energy required.");
540  }
541 
542  // Allocate energy
543  GEnergy energy;
544 
545  // Case A: Index is not -1
546  if (index() != -1.0) {
547  double exponent = index() + 1.0;
548  double e_max = std::pow(emax.MeV(), exponent);
549  double e_min = std::pow(emin.MeV(), exponent);
550  double u = ran.uniform();
551  double eng = (u > 0.0)
552  ? std::exp(std::log(u * (e_max - e_min) + e_min) / exponent)
553  : 0.0;
554  energy.MeV(eng);
555  }
556 
557  // Case B: Index is -1
558  else {
559  double e_max = std::log(emax.MeV());
560  double e_min = std::log(emin.MeV());
561  double u = ran.uniform();
562  double eng = std::exp(u * (e_max - e_min) + e_min);
563  energy.MeV(eng);
564  }
565 
566  // Return energy
567  return energy;
568 }
569 
570 
571 /***********************************************************************//**
572  * @brief Read model from XML element
573  *
574  * @param[in] xml XML element.
575  *
576  * Reads the spectral information from an XML element.
577  ***************************************************************************/
579 {
580  // Get parameter pointers
585 
586  // Read parameters
587  m_flux.read(*flux);
588  m_index.read(*index);
589  m_emin.read(*emin);
590  m_emax.read(*emax);
591 
592  // Return
593  return;
594 }
595 
596 
597 /***********************************************************************//**
598  * @brief Write model into XML element
599  *
600  * @param[in] xml XML element.
601  *
602  * @exception GException::model_invalid_spectral
603  * Existing XML element is not of type "PowerLaw"
604  *
605  * Writes the spectral information into an XML element.
606  ***************************************************************************/
608 {
609  // Set model type
610  if (xml.attribute("type") == "") {
611  xml.attribute("type", type());
612  }
613 
614  // Verify model type
615  if (xml.attribute("type") != type()) {
617  "Spectral model is not of type \""+type()+"\".");
618  }
619 
620  // Get XML parameters
625 
626  // Write parameters
627  m_flux.write(*flux);
628  m_index.write(*index);
629  m_emin.write(*emin);
630  m_emax.write(*emax);
631 
632  // Return
633  return;
634 }
635 
636 
637 /***********************************************************************//**
638  * @brief Print power law information
639  *
640  * @param[in] chatter Chattiness.
641  * @return String containing power law information.
642  ***************************************************************************/
643 std::string GModelSpectralPlawPhotonFlux::print(const GChatter& chatter) const
644 {
645  // Initialise result string
646  std::string result;
647 
648  // Continue only if chatter is not silent
649  if (chatter != SILENT) {
650 
651  // Append header
652  result.append("=== GModelSpectralPlawPhotonFlux ===");
653 
654  // Append information
655  result.append("\n"+gammalib::parformat("Number of parameters"));
656  result.append(gammalib::str(size()));
657  for (int i = 0; i < size(); ++i) {
658  result.append("\n"+m_pars[i]->print(chatter));
659  }
660 
661  } // endif: chatter was not silent
662 
663  // Return result
664  return result;
665 }
666 
667 
668 /*==========================================================================
669  = =
670  = Private methods =
671  = =
672  ==========================================================================*/
673 
674 /***********************************************************************//**
675  * @brief Initialise class members
676  ***************************************************************************/
678 {
679  // Initialise model type
680  m_type = "PowerLaw";
681 
682  // Initialise photon flux
683  m_flux.clear();
684  m_flux.name("PhotonFlux");
685  m_flux.unit("ph/cm2/s");
686  m_flux.scale(1.0);
687  m_flux.value(1.0); // default: 1.0
688  m_flux.range(0.0, 10.0); // range: [0,10]
689  m_flux.free();
690  m_flux.gradient(0.0);
691  m_flux.has_grad(true);
692 
693  // Initialise powerlaw index
694  m_index.clear();
695  m_index.name("Index");
696  m_index.scale(1.0);
697  m_index.value(-2.0); // default: -2.0
698  m_index.range(-10.0,+10.0); // range: [-10,+10]
699  m_index.free();
700  m_index.gradient(0.0);
701  m_index.has_grad(true);
702 
703  // Initialise lower limit
704  m_emin.clear();
705  m_emin.name("LowerLimit");
706  m_emin.unit("MeV");
707  m_emin.scale(1.0);
708  m_emin.value(100.0); // default: 100
709  m_emin.range(0.001, 1.0e15); // range: [0.001, 1e15]
710  m_emin.fix();
711  m_emin.gradient(0.0);
712  m_emin.has_grad(false);
713 
714  // Initialise upper limit
715  m_emax.clear();
716  m_emax.name("UpperLimit");
717  m_emax.unit("MeV");
718  m_emax.scale(1.0);
719  m_emax.value(500000.0); // default: 500000
720  m_emin.range(0.001, 1.0e15); // range: [0.001, 1e15]
721  m_emax.fix();
722  m_emax.gradient(0.0);
723  m_emax.has_grad(false);
724 
725  // Set parameter pointer(s)
726  m_pars.clear();
727  m_pars.push_back(&m_flux);
728  m_pars.push_back(&m_index);
729  m_pars.push_back(&m_emin);
730  m_pars.push_back(&m_emax);
731 
732  // Initialise last parameters (for fast computation)
733  m_log_emin = 0.0;
734  m_log_emax = 0.0;
735  m_pow_emin = 0.0;
736  m_pow_emax = 0.0;
737  m_norm = 0.0;
738  m_power = 0.0;
739  m_last_index = 1000.0;
740  m_last_emin.MeV(0.0);
741  m_last_emax.MeV(0.0);
742  m_last_energy.MeV(0.0);
743 
744  // Return
745  return;
746 }
747 
748 
749 /***********************************************************************//**
750  * @brief Copy class members
751  *
752  * @param[in] model Spectral power law model.
753  ***************************************************************************/
755 {
756  // Copy members
757  m_type = model.m_type;
758  m_flux = model.m_flux;
759  m_index = model.m_index;
760  m_emin = model.m_emin;
761  m_emax = model.m_emax;
762 
763  // Set parameter pointer(s)
764  m_pars.clear();
765  m_pars.push_back(&m_flux);
766  m_pars.push_back(&m_index);
767  m_pars.push_back(&m_emin);
768  m_pars.push_back(&m_emax);
769 
770  // Copy bookkeeping information
771  m_log_emin = model.m_log_emin;
772  m_log_emax = model.m_log_emax;
773  m_pow_emin = model.m_pow_emin;
774  m_pow_emax = model.m_pow_emax;
775  m_norm = model.m_norm;
776  m_power = model.m_power;
777  m_last_index = model.m_last_index;
778  m_last_emin = model.m_last_emin;
779  m_last_emax = model.m_last_emax;
781 
782  // Return
783  return;
784 }
785 
786 
787 /***********************************************************************//**
788  * @brief Delete class members
789  ***************************************************************************/
791 {
792  // Return
793  return;
794 }
795 
796 
797 /***********************************************************************//**
798  * @brief Update precomputed values
799  *
800  * @param[in] srcEng Source energy
801  ***************************************************************************/
803 {
804  // Compute index+1
805  double gamma = index() + 1.0;
806 
807  // Change in spectral index?
808  if (index() != m_last_index) {
809 
810  // Save actual spectral index
811  m_last_index = index();
812 
813  // Change in energy boundaries?
814  if (emin() != m_last_emin || emax() != m_last_emax) {
815  m_log_emin = std::log(emin().MeV());
816  m_last_emin = emin();
817  m_log_emax = std::log(emax().MeV());
818  m_last_emax = emax();
819  }
820 
821  // Compute normalization factors
822  if (gamma != 0.0) {
823  m_pow_emin = std::pow(emin().MeV(), gamma);
824  m_pow_emax = std::pow(emax().MeV(), gamma);
825  double d = m_pow_emax - m_pow_emin;
826  m_norm = gamma / d;
827  m_g_norm = 1.0/gamma -
828  (m_pow_emax*m_log_emax - m_pow_emin*m_log_emin)/d;
829  }
830  else {
831  m_norm = 1.0 / (m_log_emax - m_log_emin);
832  m_g_norm = 0.0;
833  }
834 
835  // Update power law factor
836  m_power = std::pow(srcEng.MeV(), index());
837  m_last_energy = srcEng;
838 
839  } // endif: change in spectral index
840 
841  // ... no change in spectral index
842  else {
843 
844  // Change in energy boundaries?
845  if (emin() != m_last_emin || emax() != m_last_emax) {
846 
847  // Update energy boundaries
848  m_log_emin = std::log(emin().MeV());
849  m_last_emin = emin();
850  m_log_emax = std::log(emax().MeV());
851  m_last_emax = emax();
852 
853  // Compute power law normalization
854  if (gamma != 0.0) {
855  m_pow_emin = std::pow(emin().MeV(), gamma);
856  m_pow_emax = std::pow(emax().MeV(), gamma);
857  double d = m_pow_emax - m_pow_emin;
858  m_norm = gamma / d;
859  m_g_norm = 1.0/gamma -
860  (m_pow_emax*m_log_emax - m_pow_emin*m_log_emin)/d;
861  }
862  else {
863  m_norm = 1.0 / (m_log_emax - m_log_emin);
864  m_g_norm = 0.0;
865  }
866 
867  } // endif: change in energy boundaries
868 
869  // Change in energy?
870  if (srcEng != m_last_energy) {
871  m_power = std::pow(srcEng.MeV(), index());
872  m_last_energy = srcEng;
873  }
874 
875  } // endelse: no change in spectral index
876 
877  // Return
878  return;
879 }
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 gradient factor.
GEnergy m_last_emax
Last upper energy limit.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:821
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:353
double index(void) const
Return power law index.
std::vector< GModelPar * > m_pars
Parameter pointers.
XML element node class.
Definition: GXmlElement.hpp:47
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:54
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:184
const double & scale(void) const
Return parameter scale.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
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.
const GXmlAttribute * attribute(const int &index) const
Return attribute.
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: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.
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.
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:1332
Exception handler interface definition.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
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:1022
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: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.
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:413