From 2d472bd44e51a171420212896543980a7f8cee7f Mon Sep 17 00:00:00 2001
From: Dominik Demuth <dominik.demuth@physik.tu-darmstadt.de>
Date: Wed, 19 Apr 2023 18:24:29 +0200
Subject: [PATCH] lower epsabs for integration of energy barrier spec dens

---
 Makefile                            |  4 ++--
 src/nmreval/distributions/energy.py |  8 ++++----
 src/nmreval/distributions/helper.py | 20 +++++++++++---------
 src/nmreval/models/fieldcycling.py  |  2 +-
 4 files changed, 18 insertions(+), 16 deletions(-)

diff --git a/Makefile b/Makefile
index f0b77fa..24001f8 100755
--- a/Makefile
+++ b/Makefile
@@ -23,8 +23,8 @@ LDFLAGS = -shared
 
 C_DIR = src/nmreval/clib
 
-all : ui
- 
+all : ui compile
+
 ui : $(PYQT_UI)
 
 rcc : $(PNG_FILES)
diff --git a/src/nmreval/distributions/energy.py b/src/nmreval/distributions/energy.py
index 3ef11df..1120b78 100644
--- a/src/nmreval/distributions/energy.py
+++ b/src/nmreval/distributions/energy.py
@@ -1,5 +1,5 @@
-import ctypes
 from itertools import product
+from ctypes import c_double, cast, pointer, c_void_p
 
 import numpy as np
 from scipy import LowLevelCallable
@@ -82,9 +82,9 @@ class EnergyBarriers(Distribution):
     def spec_dens_c(omega: np.ndarray, temperature: np.ndarray, tau0: float, e_m: float, e_b: float) -> np.ndarray:
         res = []
         for o, t in product(omega, temperature):
-            c = (ctypes.c_double * 5)(o, tau0, e_m, e_b, t)
-            user_data = ctypes.cast(ctypes.pointer(c), ctypes.c_void_p)
-            area = quad(LowLevelCallable(lib.energyDist_SD, user_data), 0, np.infty, epsabs=1e-10)[0]
+            c = (c_double * 5)(o, tau0, e_m, e_b, t)
+            user_data = cast(pointer(c), c_void_p)
+            area = quad(LowLevelCallable(lib.energyDist_SD, user_data), 0, np.infty, epsabs=1e-13)[0]
             res.append(area)
 
         return np.array(res)
diff --git a/src/nmreval/distributions/helper.py b/src/nmreval/distributions/helper.py
index 27c4db2..11c0ba9 100644
--- a/src/nmreval/distributions/helper.py
+++ b/src/nmreval/distributions/helper.py
@@ -1,24 +1,26 @@
 
 from pathlib import Path
-import ctypes
+from ctypes import CDLL, c_double, c_void_p
 
 from ..lib.logger import logger
 
 
-
 lib = None
 try:
-    lib = ctypes.CDLL(str(Path(__file__).parents[1] / 'clib' / 'integrate.so'))
+    lib = CDLL(str(Path(__file__).parents[1] / 'clib' / 'integrate.so'))
 
     # FFHS integrand
-    lib.ffhsSD.restype = ctypes.c_double
-    lib.ffhsSD.argtypes = (ctypes.c_double, ctypes.c_void_p)
+    lib.ffhsSD.restype = c_double
+    lib.ffhsSD.argtypes = (c_double, c_void_p)
 
     # Log-Gaussian integrands
-    lib.logGaussianSD_high.restype = ctypes.c_double
-    lib.logGaussianSD_high.argtypes = (ctypes.c_double, ctypes.c_void_p)
-    lib.logGaussianSD_low.restype = ctypes.c_double
-    lib.logGaussianSD_low.argtypes = (ctypes.c_double, ctypes.c_void_p)
+    lib.logGaussianSD_high.restype = c_double
+    lib.logGaussianSD_high.argtypes = (c_double, c_void_p)
+    lib.logGaussianSD_low.restype = c_double
+    lib.logGaussianSD_low.argtypes = (c_double, c_void_p)
+
+    lib.energyDist_SD.restype = c_double
+    lib.energyDist_SD.argtypes = (c_double, c_void_p)
 
     HAS_C_FUNCS = True
     logger.info('Use C functions')
diff --git a/src/nmreval/models/fieldcycling.py b/src/nmreval/models/fieldcycling.py
index 1da10c5..4d5221e 100644
--- a/src/nmreval/models/fieldcycling.py
+++ b/src/nmreval/models/fieldcycling.py
@@ -86,7 +86,7 @@ class EnergyFC(_AbstractFC):
     name = 'Energy distribution'
     params = ['C', 'T'] + EnergyBarriers.parameter
     bounds = [(0, None), (0, None), (0, None), (0, None)]
-    ralax = Relaxation(distribution=EnergyBarriers)
+    relax = Relaxation(distribution=EnergyBarriers)
 
 
 class _AbstractFCDipolar(_AbstractFC):