From e19795b51cba08e299297f500cf117a04bbcc917 Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Mon, 9 Dec 2024 12:28:38 +0100 Subject: [PATCH] wide-line spectra handle missing x values better --- src/gui_qt/main/management.py | 4 +- src/nmreval/models/wideline.py | 120 +++++++++++++++++++++------------ 2 files changed, 81 insertions(+), 43 deletions(-) diff --git a/src/gui_qt/main/management.py b/src/gui_qt/main/management.py index 8c923ae..186c9f8 100644 --- a/src/gui_qt/main/management.py +++ b/src/gui_qt/main/management.py @@ -542,7 +542,9 @@ class UpperManagement(QtCore.QObject): elif fit_limits[0] == 'in': inside = np.where((_x >= fit_limits[1][0]) & (_x <= fit_limits[1][1])) else: - inside = np.where((_x < fit_limits[1][0]) | (_x > fit_limits[1][1])) + x_lim, _ = self.graphs[self.current_graph].ranges + inside_graph = (_x >= x_lim[0]) & (_x <= x_lim[1]) + inside = np.where(((_x < fit_limits[1][0]) | (_x > fit_limits[1][1])) & inside_graph) try: if isinstance(we, str): diff --git a/src/nmreval/models/wideline.py b/src/nmreval/models/wideline.py index ac5176b..74bee7b 100644 --- a/src/nmreval/models/wideline.py +++ b/src/nmreval/models/wideline.py @@ -3,11 +3,42 @@ try: from scipy.integrate import simpson except ImportError: from scipy.integrate import simps as simpson -from numpy import pi from ..math.orientations import zcw_spherical as crystallites +__all__ = ['CSA', 'Pake', 'SecCentralLine'] + + +def _make_broadening(x: np.ndarray, sigma: float, mode: str): + dx = x[1] - x[0] + _x = np.arange(len(x)) * dx + _x -= 0.5 * _x[-1] + if mode == 'l': + apd = 2 * sigma / (4*_x**2 + sigma**2) / np.pi + else: + ln2 = np.log(2) + apd = np.exp(-4*ln2 * (_x/sigma)**2) * 2 * np.sqrt(ln2/np.pi) / sigma + return apd + + +def _make_bins(x: np.ndarray) -> np.ndarray: + bins = 0.5 * (x[1:] + x[:-1]) + return np.r_[0.5 * (-x[1] + 3 * x[0]), bins, 0.5 * (3 * x[-1] - x[-2])] + + +def _make_x(x: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + _x = x + dx = x[1:] - x[:-1] + dx = np.min(dx) + width = x[-1] - x[0] + _x = np.arange(width/dx - 1) * dx + x[0] + + bins = (_x[1:] + _x[:-1]) / 2 + bins = np.r_[_x[0]-dx/2, bins, _x[-1] + dx/2] + return _x, bins + + class Pake: type = 'Spectrum' name = 'Pake' @@ -17,38 +48,39 @@ class Pake: choices = [('Broadening', 'broad', {'Gaussian': 'g', 'Lorentzian': 'l'})] @staticmethod - def func(x, c, delta, eta, sigma, t_pulse, broad='g'): + def func( + x: np.ndarray, + c: float, + delta: float, + eta: float, + sigma: float, + t_pulse: float, + broad: str = 'g', + ) -> np.ndarray: a, b, _ = crystallites(100000) - bins = 0.5 * (x[1:] + x[:-1]) - bins = np.r_[0.5*(3*x[0]-x[1]), bins, 0.5*(3*x[-1]-x[-2])] - omega = delta * 0.5 * (3*np.cos(b)**2 - 1 - eta * np.sin(b)**2 * np.cos(2*a)) + x_used, bins = _make_x(x) s_left = np.histogram(omega, bins=bins)[0] s_right = np.histogram(-omega, bins=bins)[0] s = s_left + s_right if sigma != 0: - _x = np.arange(len(x))*(x[1]-x[0]) - _x -= 0.5*_x[-1] - - if broad == 'l': - apd = 2 * sigma / (4 * _x**2 + sigma**2) / pi - else: - apd = np.exp(-4 * np.log(2) * (_x/sigma)**2) * 2 * np.sqrt(np.log(2) / pi) / sigma - + apd = _make_broadening(x_used, sigma, broad) ret_val = np.convolve(s, apd, mode='same') - else: ret_val = s - - omega_1 = pi/2/t_pulse - attn = omega_1 * np.sin(t_pulse*np.sqrt(omega_1**2+0.5*(2*pi*x)**2)) / \ - np.sqrt(omega_1**2+(np.pi*x)**2) + + omega_1 = np.pi/2/t_pulse + attn = omega_1 * np.sin(t_pulse*np.sqrt(omega_1**2 + 0.5*(2*np.pi*x_used)**2)) / np.sqrt(omega_1**2 + (np.pi*x_used)**2) ret_val *= attn + ret_val /= simpson(y=ret_val, x=x_used) - return c * ret_val / simpson(ret_val, x) + if x_used.size == x.size: + return c * ret_val + else: + return c * np.interp(x=x, xp=x_used, fp=ret_val) class CSA: @@ -60,28 +92,29 @@ class CSA: choices = [('Broadening', 'broad', {'Gaussian': 'g', 'Lorentzian': 'l'})] @staticmethod - def func(x, c, delta, eta, w_iso, sigma, broad='g'): + def func( + x: np.ndarray, + c: float, + delta: float, + eta: float, + w_iso: float, + sigma: float, + broad: str = 'g', + ) -> np.ndarray: a, b, _ = crystallites(100000) - bins = 0.5 * (x[1:] + x[:-1]) - bins = np.r_[0.5*(-x[1] + 3*x[0]), bins, 0.5*(3*x[-1] - x[-2])] omega = w_iso + delta * 0.5 * (3*np.cos(b)**2 - 1 - eta * np.sin(b)**2 * np.cos(2*a)) - s_left = np.histogram(omega, bins=bins)[0] - s = s_left + s = np.histogram(omega, bins=_make_bins(x))[0] if sigma != 0: - _x = np.arange(len(x)) * (x[1] - x[0]) - _x -= 0.5 * _x[-1] - if broad == 'l': - apd = 2 * sigma / (4*_x**2 + sigma**2) / pi - else: - apd = np.exp(-4 * np.log(2) * (_x / sigma) ** 2) * 2 * np.sqrt(np.log(2) / pi) / sigma + print(len(s)) + apd = _make_broadening(x, sigma, broad) ret_val = np.convolve(s, apd, mode='same') else: ret_val = s - return c * ret_val / simpson(ret_val, x) + return c * ret_val / simpson(y=ret_val, x=x) class SecCentralLine: @@ -94,10 +127,18 @@ class SecCentralLine: ('Broadening', 'broad', {'Gaussian': 'g', 'Lorentzian': 'l'})] @staticmethod - def func(x, c, cq, eta, f_iso, gb, f_l, spin=2.5, broad='g'): + def func( + x: np.ndarray, + c: float, + cq: float, + eta: float, + f_iso: float, + gb: float, + f_l: float, + spin: float = 2.5, + broad: str = 'g', + ) -> np.ndarray: a, b, _ = crystallites(200000) - bins = 0.5 * (x[1:] + x[:-1]) - bins = np.r_[0.5*(-x[1] + 3*x[0]), bins, 0.5*(3*x[-1] - x[-2])] # coupling constant omega_q = 2 * np.pi * cq / (2*spin*(2*spin-1)) @@ -116,17 +157,12 @@ class SecCentralLine: orient += prefactor_c omega = 2*np.pi*f_iso + coupling * orient - s = np.histogram(omega / (2*np.pi), bins=bins)[0] + s = np.histogram(omega / (2*np.pi), bins=_make_bins(x))[0] if gb != 0: - _x = np.arange(len(x)) * (x[1]-x[0]) - _x -= 0.5*_x[-1] - if broad == 'l': - apd = 2*gb / (4*_x**2 + gb**2) / np.pi - else: - apd = np.exp(-4*np.log(2) * (_x/gb)**2) * 2 * np.sqrt(np.log(2)/np.pi) / gb + apd = _make_broadening(x, gb, broad) ret_val = np.convolve(s, apd, mode='same') else: ret_val = s - return c * ret_val / simpson(ret_val, x) + return c * ret_val / simpson(y=ret_val, x=x)