fixed LG+CC spectral density calculation for omega=0; closes #118

src/nmreval/distributions/colecole.py

@@ -54,6 +54,9 @@ class ColeCole(Distribution):
             tau (array_like):
             alpha (float):
         """
+        if alpha == 1:
+            return tau / (1 + omega**2 * tau**2)
+
         omtau = (omega*tau)**alpha
         return np.sin(alpha*np.pi/2) * omtau / (1 + omtau**2 + 2*np.cos(alpha*np.pi/2)*omtau) / omega
 

src/nmreval/distributions/loggaussian.py

@@ -5,6 +5,7 @@ from typing import Callable
 
 import numpy as np
 from scipy import LowLevelCallable
+from scipy.special import erf
 
 from nmreval.lib.utils import ArrayLike
 
@@ -32,7 +33,7 @@ class LogGaussian(Distribution):
         return np.exp(-0.5*(np.log(tau/tau0)/sigma)**2)/np.sqrt(2*np.pi)/sigma
 
     @staticmethod
-    def correlation(t, tau0, sigma: float):
+    def correlation(t: ArrayLike, tau0: ArrayLike, sigma: float):
         _t = np.atleast_1d(t)
         _tau = np.atleast_1d(tau0)
 
@@ -44,7 +45,7 @@ class LogGaussian(Distribution):
         return res.squeeze()
 
     @staticmethod
-    def susceptibility(omega, tau0, sigma: float):
+    def susceptibility(omega: ArrayLike, tau0: ArrayLike, sigma: float):
         _omega = np.atleast_1d(omega)
         _tau = np.atleast_1d(tau0)
 
@@ -68,6 +69,7 @@ class LogGaussian(Distribution):
             ret_val = _integration_parallel(_omega, _tau, sigma, _integrate_process_imag)
 
         ret_val /= _omega[:, None]
+        ret_val[_omega == 0, :] = tau[None, :] * np.exp(sigma**2 / 2)
 
         return ret_val.squeeze()
 
@@ -113,18 +115,16 @@ def _integrate_susc_c(lowfunc, highfunc, omega, tau, sigma):
     return res
 
 
-def _integrate_process_imag(args):
-    omega_i, tau_j, sigma = args
-    area = quad(_integrand_freq_imag_high, 0, 50, args=(omega_i, tau_j, sigma), epsabs=1e-12, epsrel=1e-12)[0]
-    area += quad(_integrand_freq_imag_low, -50, 0, args=(omega_i, tau_j, sigma), epsabs=1e-12, epsrel=1e-12)[0]
+def _integrate_process_imag(omega, tau, sigma):
+    area = quad(_integrand_freq_imag_high, 0, 50, args=(omega, tau, sigma), epsabs=1e-12, epsrel=1e-12)[0]
+    area += quad(_integrand_freq_imag_low, -50, 0, args=(omega, tau, sigma), epsabs=1e-12, epsrel=1e-12)[0]
 
     return area
 
 
-def _integrate_process_real(args):
-    omega_i, tau_j, sigma = args
-    area = quad(_integrand_freq_real_high, 0, 50, args=(omega_i, tau_j, sigma))[0]
-    area += quad(_integrand_freq_real_low, -50, 0, args=(omega_i, tau_j, sigma))[0]
+def _integrate_process_real(omega: float, tau: float, sigma: float):
+    area = quad(_integrand_freq_real_high, 0, 50, args=(omega, tau, sigma))[0]
+    area += quad(_integrand_freq_real_low, -50, 0, args=(omega, tau, sigma))[0]
 
     return area
 
@@ -145,9 +145,8 @@ def _integrate_correlation_c(t, tau, sigma):
     return res
 
 
-def _integrate_process_time(args):
-    omega_i, tau_j, sigma = args
-    return quad(_integrand_time, -50, 50, args=(omega_i, tau_j, sigma), epsabs=1e-12, epsrel=1e-12)[0]
+def _integrate_process_time(omega, tau, sigma):
+    return quad(_integrand_time, -50, 50, args=(omega, tau, sigma), epsabs=1e-12, epsrel=1e-12)[0]
 
 
 def _integrand_time(u, t, tau, sigma):