QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

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.
This commit is contained in:
Marek Nečada 2020-09-02 10:51:40 +03:00
parent 6189918348
commit 3479721455
2 changed files with 26 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
qpms/bessel.c
View File

@ -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)
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 {

13
qpms/qpms_error.h
View File

@ -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().