Alexander Ziegler / gammalib

Commit 95c5fce05607a77bba5d5fbb43038ea35046c619

Authored by Jürgen Knödlseder
1 parent 4dc6a2e5

Avoid overflow in powerlaw flux computation. 

The original powerlaw integration functions were prone to overflows for large gammas as the function computed the pow(energy, -gamma), where -gamma can be a large positive number.
Inline Side-by-side
Showing 1 changed file with 11 additions and 8 deletions
src/support/GTools.cpp
  View file @95c5fce
... ... @@ -791,14 +791,15 @@ double gammalib::plaw_photon_flux(const double& emin, const double& emax,
791 791  
792 792 // Compute photon flux. Computations dependend on the exponent. We
793 793 // add here a kluge to assure numerical accuracy.
794   - flux = std::pow(epivot, -gamma);
  794 + double xmin = emin/epivot;
  795 + double xmax = emax/epivot;
795 796 double exponent = gamma + 1.0;
796 797 if (std::abs(exponent) > 1.0e-11) {
797   - flux *= (std::pow(emax, exponent) -
798   - std::pow(emin, exponent)) / exponent;
  798 + flux = epivot / exponent * (std::pow(xmax, exponent) -
  799 + std::pow(xmin, exponent));
799 800 }
800 801 else {
801   - flux *= (std::log(emax) - std::log(emin));
  802 + flux = epivot * (std::log(xmax) - std::log(xmin));
802 803 }
803 804  
804 805 } // endif: emax > emin
... ... @@ -836,14 +837,16 @@ double gammalib::plaw_energy_flux(const double& emin, const double& emax,
836 837  
837 838 // Compute energy flux. Computations dependend on the exponent. We
838 839 // add here a kluge to assure numerical accuracy.
839   - flux = std::pow(epivot, -gamma);
  840 + double xmin = emin/epivot;
  841 + double xmax = emax/epivot;
840 842 double exponent = gamma + 2.0;
  843 + double epivot2 = epivot * epivot;
841 844 if (std::abs(exponent) > 1.0e-11) {
842   - flux *= (std::pow(emax, exponent) -
843   - std::pow(emin, exponent)) / exponent;
  845 + flux = epivot2 / exponent * (std::pow(xmax, exponent) -
  846 + std::pow(xmin, exponent));
844 847 }
845 848 else {
846   - flux *= (std::log(emax) - std::log(emin));
  849 + flux = epivot2 * (std::log(xmax) - std::log(xmin));
847 850 }
848 851  
849 852 } // endif: emax > emin
... ...