dominik/nmreval
1
0
Fork 0
forked from IPKM/nmreval
Code Pull Requests Activity

increase precision in log-gaussian integration

Browse Source
This commit is contained in:
Dominik Demuth 2023-10-02 12:31:40 +02:00
parent cde794cb0d
commit 37a6d789bb
1 changed files with 10 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
src/nmreval/distributions/loggaussian.py
View File

@ -105,8 +105,16 @@ def _integrate_susc_c(lowfunc, highfunc, omega, tau, sigma):
c = (ctypes.c_double * 3)(o, t, sigma) c = (ctypes.c_double * 3)(o, t, sigma)
user_data = ctypes.cast(ctypes.pointer(c), ctypes.c_void_p) 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 = 0
area += quad(LowLevelCallable(lowfunc, user_data), -np.infty, 0, epsabs=1e-12, epsrel=1e-12)[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) res.append(area)