slow but working czjzek central transition (#314)

Co-authored-by: Dominik Demuth <dominik.demuth@physik.tu-darmstadt.de>
Reviewed-on: #314

src/nmreval/models/temperature.py

@@ -79,7 +79,7 @@ class Ecoop:
 class Tanaka:
     name = 'Tanaka two-state model'
     type = 'Temperature'
-    equation = r'\tau_0 exp[(E_{a}^{f} + (E_{a}^{s}-E_{a}^{f})[1/(1+exp[(\DeltaE-T \Delta\sigma)/(k_{}B T)])]/(k_{B} T)]'
+    equation = r'\tau_0 exp[(E_{a}^{f} + (E_{a}^{s}-E_{a}^{f})[1/(1+exp[(\DeltaE-T \Delta\sigma)/(k_{B} T)])]/(k_{B} T)]'
     params = [r'\tau_{0}', 'E_{a}^{f}', 'E_{a}^{s}', r'\DeltaE', r'\Delta\sigma']
     bounds = [(0, None), (0, None), (0, None), (None, 0), (None, 0)]
     choices = [('temperature', 'temp_axis', {'T': 'T', '1000 K/T': '1000/T'})]

src/nmreval/models/wideline.py

@@ -7,7 +7,7 @@ except ImportError:
 from ..math.orientations import zcw_spherical as crystallites
 
 
-__all__ = ['CSA', 'Pake', 'SecCentralLine']
+__all__ = ['CSA', 'Pake', 'SecCentralLine', 'CzjzekCT']
 
 
 def _make_broadening(x: np.ndarray, sigma: float, mode: str):
@@ -39,6 +39,11 @@ def _make_x(x: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
     return _x, bins
 
 
+def _make_quad_prefactor(cq, f_l, spin):
+    omega_q = 2 * np.pi * cq / (2 * spin * (2 * spin - 1))
+    return 1.5 * (omega_q ** 2 / (2 * np.pi * f_l)) * (spin * (spin + 1) - 0.75)
+
+
 class Pake:
     type = 'Spectrum'
     name = 'Pake'
@@ -141,8 +146,7 @@ class SecCentralLine:
         a, b, _ = crystallites(200000)
 
         # coupling constant
-        omega_q = 2 * np.pi * cq / (2*spin*(2*spin-1))
+        coupling = _make_quad_prefactor(cq, f_l, spin)
-        coupling = 1.5 * (omega_q**2 / (2*np.pi*f_l)) * (spin*(spin+1)-0.75)
 
         # orientation
         cos2phi = np.cos(2*a)
@@ -166,3 +170,73 @@ class SecCentralLine:
             ret_val = s
 
         return c * ret_val / simpson(y=ret_val, x=x)
+
+
+class CzjzekCT:
+    type = 'Spectrum'
+    name = 'Czjzek (Central Line)'
+    equation = ''
+    params = ['A', r'\sigma', 'GB', r'\omega_{L}']
+    bounds = [(0, None), (0, None), (0, None), (0, None)]
+    choices = [('Spin', 'spin', {'3/2': 1.5, '5/2': 2.5}),
+               ('Broadening', 'broad', {'Gaussian': 'g', 'Lorentzian': 'l'})]
+
+    @staticmethod
+    def func(
+            x: np.ndarray,
+            c: float,
+            sigma: float,
+            gb: float,
+            f_l: float,
+            spin: float = 1.5,
+            broad: str = 'g',
+    ) -> np.ndarray:
+
+        cq = np.linspace(0, 5 * sigma, num=200)
+        eta = np.linspace(0, 1, num=200)
+
+        a, b, _ = crystallites(10000)
+        a = a.reshape(-1, 1)
+        b = b.reshape(-1, 1)
+
+        # coupling constant
+        dist = CzjzekCT.czjzek(cq, eta[:, None], sigma)
+
+        coupling = _make_quad_prefactor(cq, f_l, spin)
+
+        # orientation
+        e_times_phi = eta * np.cos(2 * a)
+        e_times_phi_sq = e_times_phi * e_times_phi
+
+        cos_theta_square = 0.5 + 0.5 * np.cos(2 * b)  # cos^2(x) = 1/2 (1+cos(2x)
+        prefactor_a = -3.375 + 2.25 * e_times_phi - 0.375 * e_times_phi_sq
+        prefactor_b = 3.75 - 0.5 * eta ** 2 - 2 * e_times_phi + 0.75 * e_times_phi_sq
+        prefactor_c = -0.375 + (eta ** 2) / 3 - 0.25 * e_times_phi - 0.375 * e_times_phi_sq
+
+        orient = np.zeros((a.size, eta.size))
+        orient += prefactor_a * cos_theta_square ** 2
+        orient += prefactor_b * cos_theta_square
+        orient += prefactor_c
+
+        vQ = -orient[..., None] * coupling[None, None, :]
+        weights = np.ones_like(vQ) * dist
+
+        bins = _make_bins(x)
+
+        s, _ = np.histogram(
+            vQ.reshape(-1),
+            weights=weights.reshape(-1),
+            bins=bins,
+        )
+
+        if gb != 0:
+            apd = _make_broadening(x, gb, broad)
+            ret_val = np.convolve(s, apd, mode='same')
+        else:
+            ret_val = s
+
+        return c * ret_val / simpson(y=ret_val, x=x)
+
+    @staticmethod
+    def czjzek(cq, eta, sigma):
+        return np.exp(-cq ** 2 * (1 + eta**2 / 3) / 2 / sigma**2) * cq**4 * eta * (1 - eta**2 / 9) * np.sqrt(2 / np.pi) / sigma**5