From 37a6d789bbd12afb748c5d009ab5d4cd30d8d9c5 Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Mon, 2 Oct 2023 12:31:40 +0200 Subject: [PATCH] increase precision in log-gaussian integration --- src/nmreval/distributions/loggaussian.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/nmreval/distributions/loggaussian.py b/src/nmreval/distributions/loggaussian.py index ffe867f..7ff9f45 100644 --- a/src/nmreval/distributions/loggaussian.py +++ b/src/nmreval/distributions/loggaussian.py @@ -105,8 +105,16 @@ def _integrate_susc_c(lowfunc, highfunc, omega, tau, sigma): c = (ctypes.c_double * 3)(o, t, sigma) user_data = ctypes.cast(ctypes.pointer(c), ctypes.c_void_p) - area = quad(LowLevelCallable(highfunc, user_data), 0, np.infty, epsabs=1e-12, epsrel=1e-12)[0] - area += quad(LowLevelCallable(lowfunc, user_data), -np.infty, 0, epsabs=1e-12, epsrel=1e-12)[0] + area = 0 + for (func, limits) in [(highfunc, (0, np.inf)), (lowfunc, (-np.infty, 0))]: + epsabs = 1e-12 + while epsabs > 1e-25: + a = quad(LowLevelCallable(func, user_data), *limits, epsabs=epsabs, epsrel=1e-12, full_output=1) + if a[2]['last'] > 2 or a[0] < 1e-48: + break + epsabs /= 10. + + area += a[0] res.append(area)