From 34797214559c637ea660430962513ceb8b18a37d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Wed, 2 Sep 2020 10:51:40 +0300
Subject: [PATCH] Workaround for crashing spherical Bessel funcs at large
 distance.

This leads to loss of precision nevertheless, but it should
not be critical, especially if the main use is in far field patterns.
---
 qpms/bessel.c     | 17 +++++++++++++----
 qpms/qpms_error.h | 13 +++++++++++++
 2 files changed, 26 insertions(+), 4 deletions(-)

diff --git a/qpms/bessel.c b/qpms/bessel.c
index 000ee49..2b047d4 100644
--- a/qpms/bessel.c
+++ b/qpms/bessel.c
@@ -84,7 +84,8 @@ qpms_errno_t qpms_sph_bessel_realx_fill(qpms_bessel_t typ, qpms_l_t lmax, double
 
 // TODO DOC
 qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, complex double x, complex double *res) {
-  if(!cimag(x))
+  if(!cimag(x) && 
+        creal(x) < 10000) // For large arguments (around 30000 or more), gsl_sf_bessel_jl_steed_array fails
     return qpms_sph_bessel_realx_fill(typ, lmax, creal(x), res);
   else if (isnan(creal(x)) || isnan(cimag(x))) 
     for(qpms_l_t l = 0; l <= lmax; ++l) res[l] = NAN + I*NAN;
@@ -128,9 +129,17 @@ try_again: ;
       for (qpms_l_t l = 0; l <= lmax; ++l)
         res[l] = prefac * (cyr[l] + I * cyi[l]);
       if (ierr == 3) 
-        QPMS_WARN("Amos's zbes%c computation done but losses of significance "
-            "by argument reduction produce less than half of machine accuracy.", 
-            kindchar);
+        if (creal(x) >= 10000) {
+          QPMS_WARN_ONCE("Due to large argument, "
+                "Amos's zbes%c computation done but losses of significance "
+                "by argument reduction produce less than half of machine accuracy. "
+                "This warning is shown only once.", 
+                kindchar);
+        } else {
+          QPMS_WARN("Amos's zbes%c computation done but losses of significance "
+                "by argument reduction produce less than half of machine accuracy.", 
+                kindchar);
+        }
       return QPMS_SUCCESS; //TODO maybe something else if ierr == 3
     }
     else {
diff --git a/qpms/qpms_error.h b/qpms/qpms_error.h
index 0cf9f22..f554035 100644
--- a/qpms/qpms_error.h
+++ b/qpms/qpms_error.h
@@ -67,6 +67,19 @@ qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types);
  */
 #define QPMS_WARN(msg, ...) qpms_warn_at_flf(__FILE__,__LINE__,__func__,msg, ##__VA_ARGS__)
 
+/// Print a warning to stderr and flush the buffer – only once per run (or current thread).
+/**
+ * The arguments are the same as in standard printf().
+ */
+#define QPMS_WARN_ONCE(msg, ...) {\
+	static _Bool already_bitched = 0; \
+	if (QPMS_UNLIKELY(!already_bitched)) {\
+		qpms_warn_at_flf(__FILE__,__LINE__,__func__,msg, ##__VA_ARGS__);\
+		already_bitched = 1;\
+	}\
+}
+
+
 /// Print a debugging message to stderr and flush the buffer.
 /**
  * The arguments after \a type are the same as in standard printf().