From cf2aca7f035ae01ac06b02335efb109666f067d3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Tue, 17 Nov 2020 10:30:56 +0200
Subject: [PATCH] Silent NAN (with warning) instead of crash for singular
 Bessel at zero

---
 qpms/bessel.c | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/qpms/bessel.c b/qpms/bessel.c
index 2b047d4..93080d7 100644
--- a/qpms/bessel.c
+++ b/qpms/bessel.c
@@ -52,6 +52,11 @@ static inline complex double spherical_yn(qpms_l_t l, complex double z) {
 qpms_errno_t qpms_sph_bessel_realx_fill(qpms_bessel_t typ, qpms_l_t lmax, double x, complex double *result_array) {
   int retval;
   double tmparr[lmax+1];
+  if (typ != QPMS_BESSEL_REGULAR && x == 0) {
+    for (int l = 0; l <= lmax; ++l) result_array[l] = NAN + I*NAN;
+    QPMS_WARN("Evaluating singular Bessel functions at zero!");
+    return QPMS_SUCCESS; // Really?
+  }
   switch(typ) {
     case QPMS_BESSEL_REGULAR:
       retval = gsl_sf_bessel_jl_steed_array(lmax, x, tmparr);