src/nmreval/clib/diffusion.c

@@ -0,0 +1,16 @@
+/* integrands used in quadrature integration with scipy's LowLevelCallables */
+#include <math.h>
+
+double anistropicDiffusion(double x, void *user_data) {
+    double *c = (double *)user_data;
+
+    double q = c[0];
+    double t = c[1];
+    double d_perp = c[2];
+    double d_par = c[3];
+
+    double cos_theta = cos(x);
+    double sin_theta = sin(x);
+
+    return exp(-q * q * t * (d_par * cos_theta * cos_theta + d_perp * sin_theta * sin_theta)) * sin_theta;
+}

src/nmreval/distributions/helper.py

@@ -5,6 +5,17 @@ from ctypes import CDLL, c_double, c_void_p
 from ..lib.logger import logger
 
 
+diffusion_lib = None
+try:
+    diffusion_lib = CDLL(str(Path(__file__).parents[1] / 'clib' / 'diffusion.so'))
+
+    diffusion_lib.anistropicDiffusion.restype = c_double
+    diffusion_lib.anistropicDiffusion.argtypes = (c_double, c_void_p)
+
+    HAS_C_FUNCS = True
+except OSError:
+    HAS_C_FUNCS = False
+
 lib = None
 try:
     lib = CDLL(str(Path(__file__).parents[1] / 'clib' / 'integrate.so'))
@@ -39,7 +50,6 @@ try:
     lib.energyDistSuscImag.restype = c_double
     lib.energyDistSuscImag.argtypes = (c_double, c_void_p)
 
-
     HAS_C_FUNCS = True
     logger.info('Use C functions')
 except OSError:

src/nmreval/models/diffusion.py

@@ -1,7 +1,11 @@
+from ctypes import c_double, cast, c_void_p, pointer
+
 import numpy as np
-from scipy import special as special
+from scipy import special as special, LowLevelCallable
+from scipy.integrate import quad
 
 from ..utils import gamma
+from nmreval.distributions.helper import HAS_C_FUNCS, diffusion_lib
 
 
 class Diffusion:
@@ -111,6 +115,30 @@ class AnisotropicDiffusion(object):
 
         return m0 * diffs * relax
 
+    @staticmethod
+    def _integrate_c(x, m0, d_perp, d_par, trel, brel, g, t2, nucleus=2.67522128e8, t_axis='tm') -> np.ndarray:
+        diffusion_decay = np.zeros_like(x)
+
+        if t_axis == 'tm':
+            tm = x
+            tp = t2
+            relax = np.exp(-(tm/trel)**brel)
+        else:
+            tm = t2
+            tp = x
+            relax = np.exp(-(tp/trel)**brel)*np.exp(-(tp/trel)**brel)
+
+        q = g * nucleus * tp
+        t = 2 * tp / 3 + tm
+
+        for (i, t_i) in enumerate(t):
+            c = (c_double * 4)(q, t_i, d_perp, d_par)
+            user_data = cast(pointer(c), c_void_p)
+
+            diffusion_decay[i] = quad(LowLevelCallable(diffusion_lib.anistropicDiffusion, user_data), 0, np.pi, epsabs=1e-13)[0]
+
+        return m0 * diffusion_decay * relax
+
 
 class Peschier:
     name = 'Diffusion + Cross-Relaxation'