From dbb35cdba4fd2704b31fac723716bd1203725e01 Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Thu, 9 Nov 2023 17:26:57 +0100 Subject: [PATCH] added f-omega option to bds fit functions; closes #133 --- src/nmreval/models/bds.py | 122 ++++++++++++++++++++++++++++---------- 1 file changed, 91 insertions(+), 31 deletions(-) diff --git a/src/nmreval/models/bds.py b/src/nmreval/models/bds.py index 7e770e7..9ccb922 100644 --- a/src/nmreval/models/bds.py +++ b/src/nmreval/models/bds.py @@ -10,14 +10,20 @@ class _AbstractBDS: equation = '' params = [r'\Delta\epsilon', r'\tau_{0}'] bounds = [(0, None), (0, None)] + choices = [ + ('x axis', 'xaxis', {'Frequency': 'freq', 'Omega': 'omega'}) + ] susceptibility = None iscomplex = True @classmethod - def func(cls, x, *args, complex_mode: int = 0, **kwargs): + def func(cls, x, *args, complex_mode: int = 0, xaxis: str = 'freq', **kwargs): # args[0] : Delta epsilon # args[1:] : every other parameter - chi = args[0] * cls.susceptibility(2*np.pi*x, *args[1:], **kwargs) + + _w = _convert_x_to_omega(x, xaxis=xaxis) + + chi = args[0] * cls.susceptibility(_w, *args[1:], **kwargs) if complex_mode == 0: return chi elif complex_mode == 1: @@ -64,11 +70,16 @@ class HavriliakNegamiAlphaGammaBDS: equation = r'\Delta\epsilon / [1-(i\omega\tau)^{\gamma}]^{\alpha}' params = [r'\Delta\epsilon', r'\tau_{0}', r'\alpha', r'\alpha\gamma'] bounds = [(0, None), (0, None), (0, 1), (0, 1)] + choices = [ + ('x axis', 'xaxis', {'Frequency': 'freq', 'Omega': 'omega'}) + ] iscomplex = True @staticmethod - def func(x, deps, tau, alpha, alphagamma, complex_mode: int = 0, **kwargs): - chi = deps * HavriliakNegami.susceptibility(2*np.pi*x, tau, alpha, alphagamma/alpha, **kwargs) + def func(x, deps, tau, alpha, alphagamma, complex_mode: int = 0, xaxis: str = 'freq', **kwargs): + _w = _convert_x_to_omega(x, xaxis=xaxis) + + chi = deps * HavriliakNegami.susceptibility(_w, tau, alpha, alphagamma/alpha, **kwargs) if complex_mode == 0: return chi elif complex_mode == 1: @@ -115,13 +126,17 @@ class HNWithHF: equation = r'\Delta\epsilon HN(\omega, \tau, \alpha, \gamma) / CD(\omega, \tau_{c}, \alpha\gamma-\delta)' params = [r'\Delta\epsilon', r'\tau', r'\alpha', r'\gamma', r'\tau_{c}', '\delta'] bounds = [(0, None), (0, None), (0, 1), (0, 1), (0, None), (0, 1)] + choices = [ + ('x axis', 'xaxis', {'Frequency': 'freq', 'Omega': 'omega'}) + ] iscomplex = True @staticmethod - def func(x, deps, tau, alpha, gamma, tauc, delta, complex_mode: int = 0): - w = 2 * np.pi * x - numerator = (1 + 1j*w*tauc) ** (gamma-delta) - denominator = (1 + (1j*w*tau)**alpha) ** gamma + def func(x, deps, tau, alpha, gamma, tauc, delta, complex_mode: int = 0, xaxis: str = 'freq'): + _w = _convert_x_to_omega(x, xaxis=xaxis) + + numerator = (1 + 1j*_w*tauc) ** (gamma-delta) + denominator = (1 + (1j*_w*tau)**alpha) ** gamma epsilon = deps * np.conjugate(numerator / denominator) @@ -144,8 +159,8 @@ class _CCWithHF: iscomplex = True @staticmethod - def func(x, deps, tau, alpha, tauc, delta, complex_mode: int = 0): - return HNWithHF.func(x, deps, tau, alpha, 1, tauc, delta, complex_mode=complex_mode) + def func(x, deps, tau, alpha, tauc, delta, complex_mode: int = 0, xaxis: str = 'freq'): + return HNWithHF.func(x, deps, tau, alpha, 1, tauc, delta, complex_mode=complex_mode, xaxis=xaxis) class _CDWithHF: @@ -157,8 +172,8 @@ class _CDWithHF: iscomplex = True @staticmethod - def func(x, deps, tau, gamma, tauc, delta, complex_mode: int = 0): - return HNWithHF.func(x, deps, tau, 1, gamma, tauc, delta, complex_mode=complex_mode) + def func(x, deps, tau, gamma, tauc, delta, complex_mode: int = 0, xaxis: str = 'freq'): + return HNWithHF.func(x, deps, tau, 1, gamma, tauc, delta, complex_mode=complex_mode, xaxis=xaxis) class PowerLawBDS: @@ -167,16 +182,21 @@ class PowerLawBDS: equation = r'A / (i\omega)^{n}' params = ['A', 'n'] bounds = [(None, None), (None, None)] + choices = [ + ('x axis', 'xaxis', {'Frequency': 'freq', 'Omega': 'omega'}) + ] iscomplex = True @staticmethod - def func(x, a, n, complex_mode: int = 0): + def func(x, a, n, complex_mode: int = 0, xaxis: str = 'freq'): + _w = _convert_x_to_omega(x, xaxis=xaxis) + if complex_mode == 0: - ret_val = np.exp(1j*n*np.pi/2) * a / x**n + ret_val = np.exp(1j*n*np.pi/2) * a / _w**n elif complex_mode == 1: - ret_val = np.cos(n*np.pi/2) * a / x**n + ret_val = np.cos(n*np.pi/2) * a / _w**n elif complex_mode == 2: - ret_val = np.sin(n*np.pi/2) * a / x**n + ret_val = np.sin(n*np.pi/2) * a / _w**n else: raise ValueError(f'{complex_mode!r} is not 0, 1, 2') @@ -189,17 +209,22 @@ class DCCondBDS: equation = r'i\sigma_{dc}/\epsilon_{0}\omega' params = [r'\sigma_{dc}'] bounds = [(0, None)] + choices = [ + ('x axis', 'xaxis', {'Frequency': 'freq', 'Omega': 'omega'}) + ] iscomplex = True @staticmethod - def func(x, sigma, complex_mode: int = 0): + def func(x, sigma, complex_mode: int = 0, xaxis: str = 'freq'): + _w = _convert_x_to_omega(x, xaxis=xaxis) + if complex_mode == 0: ret_val = np.zeros(x.shape, dtype=complex) - ret_val += 1j * sigma / x / epsilon0 + ret_val += 1j * sigma / _w / epsilon0 elif complex_mode == 1: ret_val = np.zeros(x.shape) elif complex_mode == 2: - ret_val = sigma / x / epsilon0 + ret_val = sigma / _w / epsilon0 else: raise ValueError(f'{complex_mode!r} is not 0, 1, 2') @@ -210,11 +235,16 @@ class DerivativeHavriliakNegami: name = 'Derivative HN' type = 'Dielectric Spectroscopy' params = [r'\Delta\epsilon', r'\tau', r'\alpha', r'\gamma'] + choices = [ + ('x axis', 'xaxis', {'Frequency': 'freq', 'Omega': 'omega'}) + ] bounds = [(0, None), (0, None), (0, 1), (0, 1)] @staticmethod - def func(x, eps, tau, a, g): - omtau = 2*np.pi*x * tau + def func(x, eps, tau, a, g, xaxis: str = 'freq'): + w = _convert_x_to_omega(x, xaxis=xaxis) + + omtau = w * tau theta = np.arctan2(np.sin(np.pi*a/2.), omtau**(-a) + np.cos(np.pi*a/2.)) numer = a*g * omtau**a * np.cos(np.pi*a/2. - (1+g)*theta) @@ -228,10 +258,15 @@ class DerivativeColeCole: type = 'Dielectric Spectroscopy' params = [r'\Delta\epsilon', r'\tau', r'\alpha'] bounds = [(0, None), (0, None), (0, 1)] + choices = [ + ('x axis', 'xaxis', {'Frequency': 'freq', 'Omega': 'omega'}) + ] @staticmethod - def func(x, eps, tau, alpha): - omtau = 2*np.pi*x * tau + def func(x, eps, tau, alpha, xaxis: str = 'freq'): + w = _convert_x_to_omega(x, xaxis=xaxis) + + omtau = w * tau theta = np.arctan2(np.sin(np.pi*alpha/2.), omtau**(-alpha) + np.cos(np.pi*alpha/2.)) numer = alpha * omtau**alpha * np.cos(np.pi*alpha/2. - 2*theta) @@ -245,10 +280,14 @@ class DerivativeColeDavidson: type = 'Dielectric Spectroscopy' params = [r'\Delta\epsilon', r'\tau', r'\gamma'] bounds = [(0, None), (0, None), (0, 1)] + choices = [ + ('x axis', 'xaxis', {'Frequency': 'freq', 'Omega': 'omega'}) + ] @staticmethod - def func(x, eps, tau, g): - omtau = 2*np.pi*x * tau + def func(x, eps, tau, g, xaxis: str = 'freq'): + w = _convert_x_to_omega(x, xaxis=xaxis) + omtau = w * tau numer = g * omtau * np.sin((1+g)*np.arctan(omtau)) denom = (1 + omtau**2)**((1+g)/2.) @@ -260,10 +299,13 @@ class _DerivativeHNWithHF: type = 'Dielectric Spectroscopy' params = [r'\Delta\epsilon', r'\tau', r'\alpha', r'\gamma', r'\tau_{c}', r'\delta'] bounds = [(0, None), (0, None), (0, 1), (0, 1), (0, None), (0, 1)] + choices = [ + ('x axis', 'xaxis', {'Frequency': 'freq', 'Omega': 'omega'}) + ] @staticmethod - def func(x, deps, tau, alpha, gamma, tauc, delta): - w = 2*np.pi*x + def func(x, deps, tau, alpha, gamma, tauc, delta, xaxis: str = 'freq'): + w = _convert_x_to_omega(x, xaxis=xaxis) a_pi2 = alpha*np.pi/2 w_lamb = (w*tau)**alpha @@ -287,10 +329,13 @@ class DerivativeCCWithHF: type = 'Dielectric Spectroscopy' params = [r'\Delta\epsilon', r'\tau', r'\alpha', r'\tau_{c}', r'\delta'] bounds = [(0, None), (0, None), (0, 1), (0, None), (0, 1)] + choices = [ + ('x axis', 'xaxis', {'Frequency': 'freq', 'Omega': 'omega'}) + ] @staticmethod - def func(x, deps, tau, alpha, tauc, delta): - return _DerivativeHNWithHF.func(x, deps, tau, alpha, 1, tauc, delta) + def func(x, deps, tau, alpha, tauc, delta, xaxis: str = 'freq'): + return _DerivativeHNWithHF.func(x, deps, tau, alpha, 1, tauc, delta, xaxis=xaxis) class DerivativeCDWithHF: @@ -298,7 +343,22 @@ class DerivativeCDWithHF: type = 'Dielectric Spectroscopy' params = [r'\Delta\epsilon', r'\tau', r'\gamma', r'\tau_{c}', r'\delta'] bounds = [(0, None), (0, None), (0, 1), (0, None), (0, 1)] + choices = [ + ('x axis', 'xaxis', {'Frequency': 'freq', 'Omega': 'omega'}) + ] @staticmethod - def func(x, deps, tau, gamma, tauc, delta): - return _DerivativeHNWithHF.func(x, deps, tau, 1, gamma, tauc, delta) + def func(x, deps, tau, gamma, tauc, delta, xaxis: str = 'freq'): + return _DerivativeHNWithHF.func(x, deps, tau, 1, gamma, tauc, delta, xaxis=xaxis) + + +def _convert_x_to_omega(x, xaxis: str = 'freq'): + if xaxis not in ['freq', 'omega']: + raise ValueError(f'Argument `xaxis` is `freq` or `omega`, given is {xaxis!r}') + + if xaxis == 'freq': + _w = 2 * np.pi * x + else: + _w = x + + return _w