From b8bab2af7b4cf0f4b07f5db1dbe866eb3e02f538 Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Fri, 8 Sep 2023 17:54:27 +0200 Subject: [PATCH] fixed LG+CC spectral density calculation for omega=0; closes #118 --- src/nmreval/distributions/colecole.py | 3 +++ src/nmreval/distributions/loggaussian.py | 25 ++++++++++++------------ 2 files changed, 15 insertions(+), 13 deletions(-) diff --git a/src/nmreval/distributions/colecole.py b/src/nmreval/distributions/colecole.py index b8551a8..b1a5a66 100644 --- a/src/nmreval/distributions/colecole.py +++ b/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 diff --git a/src/nmreval/distributions/loggaussian.py b/src/nmreval/distributions/loggaussian.py index 10c3615..8979750 100644 --- a/src/nmreval/distributions/loggaussian.py +++ b/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):