[gnuastro-commits] master d6447d2b: Library (cosmology.h): slower/robust gsl integration when fast fails

[Mohammad Akhlaghi]

branch: master
commit d6447d2b26eec1577faf91d5150a36e77db29f9b
Author: Thorsten Alteholz <thorsten@alteholz.dev>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Library (cosmology.h): slower/robust gsl integration when fast fails
    
    Until now, when CosmicCalculator was given an input redshift that is higher
    than 100, GSL would produce a code-dump. This was due to the choice of
    'gsl_integration_qng' for the numerical intergration (because it is fast).
    
    With this commit, after 'gsl_integration_qng' is finished we check whether
    it was successful or not. If it failed, the slower, but more robust
    'gsl_integration_qag' is used. In order to use the 'error' function to
    report the failure, the variable that was previously named 'error' was
    renamed to 'gslerr'.
---
 lib/cosmology.c | 61 ++++++++++++++++++++++++++++++++++++++++++++++++++-------
 1 file changed, 54 insertions(+), 7 deletions(-)

diff --git a/lib/cosmology.c b/lib/cosmology.c
index a3735062..62bf775f 100644
--- a/lib/cosmology.c
+++ b/lib/cosmology.c
@@ -28,6 +28,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <stdio.h>
 #include <stdlib.h>
 
+#include <gsl/gsl_errno.h>
 #include <gsl/gsl_const_mksa.h>
 #include <gsl/gsl_integration.h>
 
@@ -249,9 +250,12 @@ gal_cosmology_proper_distance(double z, double H0, double 
o_lambda_0,
                               double o_matter_0, double o_radiation_0)
 {
   cosmology_density_check(o_lambda_0, o_matter_0, o_radiation_0);
+  int status;
   size_t neval;
   gsl_function F;
-  double result, error, c=GSL_CONST_MKSA_SPEED_OF_LIGHT;
+  gsl_integration_workspace *w;
+  gsl_error_handler_t *gslerrhandler;
+  double result, gslerr, c=GSL_CONST_MKSA_SPEED_OF_LIGHT;
   double o_curv_0 = 1.0 - ( o_lambda_0 + o_matter_0 + o_radiation_0 );
   double H0s=H0/1000/GSL_CONST_MKSA_PARSEC;  /* H0 in units of seconds. */
   struct cosmology_integrand_t p={o_lambda_0, o_curv_0, o_matter_0,
@@ -261,9 +265,29 @@ gal_cosmology_proper_distance(double z, double H0, double 
o_lambda_0,
   F.params=&p;
   F.function=&cosmology_integrand_proper_dist;
 
-  /* Do the integration. */
-  gsl_integration_qng(&F, 0.0f, z, GSLIEPSABS, GSLIEPSREL, &result,
-                      &error, &neval);
+  /* Temporarily switch off error handling */
+  gslerrhandler=gsl_set_error_handler_off();
+
+  /* Do the integration (first with the fast GSL function). */
+  status=gsl_integration_qng(&F, 0.0f, z, GSLIEPSABS, GSLIEPSREL, &result,
+                             &gslerr, &neval);
+
+  /* If the first integration failed, try a slower, but more robust one. */
+  if (status!=GSL_SUCCESS)
+    {
+      w=gsl_integration_workspace_alloc(GSLILIMIT);
+      status=gsl_integration_qag(&F, 0.0f, z, GSLIEPSABS, GSLIEPSREL,
+                                 GSLILIMIT, GSL_INTEG_GAUSS21, w,
+                                 &result, &gslerr);
+      gsl_integration_workspace_free(w);
+      if (status!=GSL_SUCCESS)
+        error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' to "
+              "fix the problem. The status of the second integration is "
+              "%i.",  __func__, PACKAGE_BUGREPORT, status);
+  }
+
+  /* Reactivate "normal" error handling */
+  gslerrhandler=gsl_set_error_handler(gslerrhandler);
 
   /* Return the result. */
   return result * c / H0s / (1e6 * GSL_CONST_MKSA_PARSEC);
@@ -279,9 +303,12 @@ gal_cosmology_comoving_volume(double z, double H0, double 
o_lambda_0,
                               double o_matter_0, double o_radiation_0)
 {
   cosmology_density_check(o_lambda_0, o_matter_0, o_radiation_0);
+  int status;
   size_t neval;
   gsl_function F;
-  double result, error;
+  double result, gslerr;
+  gsl_integration_workspace *w;
+  gsl_error_handler_t *gslerrhandler;
   double c=GSL_CONST_MKSA_SPEED_OF_LIGHT;
   double H0s=H0/1000/GSL_CONST_MKSA_PARSEC;     /* H0 in units of seconds. */
   double cH = c / H0s / (1e6 * GSL_CONST_MKSA_PARSEC);
@@ -293,9 +320,29 @@ gal_cosmology_comoving_volume(double z, double H0, double 
o_lambda_0,
   F.params=&p;
   F.function=&cosmology_integrand_comoving_volume;
 
+  /* temporarily switch off error handling */
+  gslerrhandler=gsl_set_error_handler_off();
+
   /* Do the integration. */
-  gsl_integration_qng(&F, 0.0f, z, GSLIEPSABS, GSLIEPSREL,
-                      &result, &error, &neval);
+  status=gsl_integration_qng(&F, 0.0f, z, GSLIEPSABS, GSLIEPSREL,
+                      &result, &gslerr, &neval);
+
+  /* If the first integration failed, try a slower, but more robust one. */
+  if (status!=GSL_SUCCESS)
+    {
+      w=gsl_integration_workspace_alloc(GSLILIMIT);
+     status=gsl_integration_qag(&F, 0.0f, z, GSLIEPSABS, GSLIEPSREL,
+                                GSLILIMIT, GSL_INTEG_GAUSS21, w,
+                                &result, &gslerr);
+     gsl_integration_workspace_free(w);
+     if (status!=GSL_SUCCESS)
+       error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' to "
+             "fix the problem. The status of the second integration is "
+             "%i.",  __func__, PACKAGE_BUGREPORT, status);
+    }
+
+  /* Reactivate "normal" error handling */
+  gslerrhandler=gsl_set_error_handler(gslerrhandler);
 
   /* Return the result. */
   return result * 4 * M_PI * cH*cH*cH;