Assure numerical accuracy in gammalib::plaw_integral() method

Jürgen Knödlseder
1 parent dc59816a
Showing 2 changed files with 11 additions and 7 deletions
src/numerics/GMath.cpp
test/test_GNumerics.cpp
@@ -546,14 +546,17 @@ double gammalib::plaw_integral(const double&amp; x1,
                                const double& f2)
 {
     // Compute power law slope
-    double fratio = std::log(f2/f1);
-    double xratio = std::log(x2/x1);
-    double slope  = fratio / xratio;
-
-    // Compute integral
+    double x2x1    = x2/x1;
+    double fratio  = std::log(f2/f1);
+    double xratio  = std::log(x2x1);
+    double slope   = fratio / xratio;
+    double slopep1 = slope + 1.0;
+
+    // Compute integral. Computations dependend on the slope. We add here a
+    // kluge to assure numerical accuracy.
     double integral;
-    if (slope != -1.0) {
-        integral = f1 / (slope + 1.0) * (x2 * std::pow(x2/x1, slope) - x1);
+    if (std::abs(slopep1) > 1.0e-11) {
+        integral = f1 / slopep1 * (x2 * std::pow(x2x1, slope) - x1);
     }
     else {
         integral = f1 * x1 * xratio;
@@ -90,6 +90,7 @@ void TestGNumerics::test_math(void)
     // Test power law integration
     test_value(gammalib::plaw_integral(2.0, 3.0, 9.0, 5.0), 29.112577);
     test_value(gammalib::plaw_integral(2.0, 9.0, 9.0, 2.0), 27.073393);
+    test_value(gammalib::plaw_integral(2.0, 9.0, 5.0, 2.0), 12.471065);
  
     // Exit test
     return;