diff --git a/src/nmreval/clib/integrate.c b/src/nmreval/clib/integrate.c index 28c591a..74ead29 100644 --- a/src/nmreval/clib/integrate.c +++ b/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); +} + diff --git a/src/nmreval/distributions/gengamma.py b/src/nmreval/distributions/gengamma.py index f463658..f5ab1b3 100644 --- a/src/nmreval/distributions/gengamma.py +++ b/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()