src/nmreval/clib/integrate.c

@@ -156,3 +156,20 @@ double energyDistCorrelation(double x, void *user_data) {
     return normalDist(x, e_m, e_b) * exp(-t * r);
 }
 
+// Generalised Gamma Function
+double genGammaAlphaDist(double x, void *user_data) {
+    double *c = (double *)user_data;
+
+    double omega = c[0];
+    double tau0 = c[1];
+    double alpha = c[2];
+    double beta = c[3];
+
+    double b_to_a = beta / alpha;
+    double tau_to_tau0 = tau / tau0;
+
+    double norm = exp(-lgamma(b_to_a) + b_to_a * log(b_to_a)) * alpha;
+
+    return norm * exp(-b_to_a * pow(tau_to_tau0, alpha)) * pow(tau_to_tau0, beta);
+}
+

src/nmreval/distributions/gengamma.py

@@ -36,7 +36,10 @@ class AbstractGG(Distribution, ABC):
         taus, ln_tau = AbstractGG._prepare_integration(tau0)
         g_tau = cls.distribution(taus, tau0, *args)
 
-        ret_val = np.array([simpson(g_tau / (1 - 1j*w_i*taus), ln_tau) for w_i in w]).squeeze()
+        ret_val = np.zeros_like(omega, dtype=np.complex128)
+
+        ret_val.real = np.array([simpson(g_tau * taus / (1 + (w_i*taus)**2), ln_tau) for w_i in w]).squeeze()
+        ret_val.imag = np.array([simpson(g_tau * w_i * taus / (1 + (w_i*taus)**2), ln_tau) for w_i in w]).squeeze()
 
         return ret_val
 
@@ -50,7 +53,7 @@ class AbstractGG(Distribution, ABC):
         ret_val = np.zeros((w.size, _t.size))
 
         for i, tau_i in enumerate(_t):
-            taus, ln_tau = AbstractGG._prepare_integration(tau_i)
+            taus, ln_tau = AbstractGG._prepare_integration(tau_i, limits=limits, num_steps=num_steps)
             g_tau = cls.distribution(taus, tau_i, *args)
 
             ret_val[:, i] = np.array([simpson(g_tau * taus / (1 + (w_i*taus)**2), ln_tau) for w_i in w]).squeeze()