diff --git a/src/nmreval/distributions/helper.py b/src/nmreval/distributions/helper.py
index 1138cdd..4927831 100644
--- a/src/nmreval/distributions/helper.py
+++ b/src/nmreval/distributions/helper.py
@@ -55,4 +55,3 @@ try:
 except OSError:
     HAS_C_FUNCS = False
     logger.info('Use python functions')
-
diff --git a/src/nmreval/models/diffusion.py b/src/nmreval/models/diffusion.py
index f43594f..c479d06 100644
--- a/src/nmreval/models/diffusion.py
+++ b/src/nmreval/models/diffusion.py
@@ -98,27 +98,6 @@ class AnisotropicDiffusion(object):
 
     @staticmethod
     def func(x, m0, d_perp, d_par, trel, brel, g, t2, nucleus=2.67522128e8, t_axis='tm'):
-        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_squared = np.power(g * nucleus * tp, 2)
-        t = 2 * tp / 3 + tm
-        z = np.sqrt(q_squared * (d_par - d_perp) * t)
-        # Callaghan eq (6.89)
-        diffs = np.exp(-q_squared*t*d_perp) * special.erf(z) / z
-
-        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
@@ -130,6 +109,18 @@ class AnisotropicDiffusion(object):
 
         q = g * nucleus * tp
         t = 2 * tp / 3 + tm
+        # Callaghan eq (6.89)
+        if HAS_C_FUNCS:
+            diffusion_decay = AnisotropicDiffusion._integrate_c(q, t, d_perp, d_par)
+        else:
+            z = np.sqrt(q**2 * (d_par - d_perp) * t)
+            diffusion_decay = np.exp(-q**2 * t * d_perp) * special.erf(z) / z
+
+        return m0 * diffusion_decay * relax
+
+    @staticmethod
+    def _integrate_c(q, t, d_perp, d_par) -> np.ndarray:
+        diffusion_decay = np.zeros_like(t)
 
         for (i, t_i) in enumerate(t):
             c = (c_double * 4)(q, t_i, d_perp, d_par)
@@ -137,7 +128,7 @@ class AnisotropicDiffusion(object):
 
             diffusion_decay[i] = quad(LowLevelCallable(diffusion_lib.anistropicDiffusion, user_data), 0, np.pi, epsabs=1e-13)[0]
 
-        return m0 * diffusion_decay * relax
+        return diffusion_decay
 
 
 class Peschier: