|
|
|
@ -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)
|
|
|
|
|