WIP good results (needs cleanup)

radicalpy/data.py

@@ -582,6 +582,9 @@ def semiclassical_std(self) -> float:
 
         Calculate the standard deviation :math:`\sigma` where
 
+        .. todo::
+           Do the math properly.
+
         .. math::
            \sigma = \sqrt{\frac{2}{\tau^2}}
 
@@ -603,8 +606,7 @@ def semiclassical_std(self) -> float:
             n.spin_quantum_number * (n.spin_quantum_number + 1) * n.hfc.isotropic**2
             for n in self.nuclei
         )
-        tau2 = 6.0 / tmp
-        return np.sqrt(2 / tau2)
+        return tmp**2 / 6
 
 
 class Triplet(Molecule):

radicalpy/experiments.py

@@ -98,7 +98,7 @@ def semiclassical_mary(
     M = 16  # number of spin states
     trace = np.zeros((num_samples, len(ts)))
     mary = np.zeros((len(ts), len(Bs)))
-    gen = sim.semiclassical_gen(num_samples)
+    gen = sim.semiclassical_gen2(num_samples)
 
     for i, B0 in enumerate(tqdm(Bs)):
         Hz = sim.zeeman_hamiltonian(B0)

radicalpy/simulation.py

@@ -1254,17 +1254,15 @@ def semiclassical_gen2(
         self,
         num_samples: int,
     ) -> np.ndarray:
-        num_particles = len(self.radicals)
-        assert num_particles == 2
-        assert num_particles == len(self.molecules)
+        assert len(self.radicals) == 2
         assert self.radicals[0].multiplicity == 2
         assert self.radicals[1].multiplicity == 2
 
         spinops = np.array([self.spin_operator(0, ax) for ax in "xyz"])
-        result = np.zeros((num_samples, 4, 4), dtype=complex)
         cov = np.diag([m.semiclassical_std for m in self.molecules])
         samples = np.random.multivariate_normal([0, 0], cov, size=(num_samples, 3))
-        return np.einsum("nam,axy->nxy", samples, spinops)
+        result = np.einsum("nam,axy->nxy", samples, spinops)
+        return result * self.radicals[0].gamma_mT
 
     @property
     def nuclei(self):

utils/wip_experiments.py

@@ -113,18 +113,20 @@ def algebra():
     from sympy.abc import T, h, l, n, sigma, tau, x
 
     pexpr = (tau**2 / (4 * pi)) ** (3 / 2) * exp(-1 / 4 * x**2 * tau**2)
-    bexpr = (3 / (2 * pi * n * l**2)) ** (3 / 2) * exp(-(3 * h**2) / (2 * n * l**2))
-    nexpr = 1 / ((sigma**2 * 2 * pi) ** (1 / 2)) * exp(-1 / 2 * x**2 / sigma**2)
-
-    pexpr = pexpr.subs(tau**2, 1 / T)
-    pexpr = pexpr.subs(x, h)
+    pexpr = pexpr.subs(tau**2, 1 / T).subs(x, h)
     print(pexpr)
+
+    bexpr = (3 / (2 * pi * n * l**2)) ** (3 / 2) * exp(-(3 * h**2) / (2 * n * l**2))
     bexpr = bexpr.subs(n * l**2, T * 6)
     print(bexpr)
-    nexpr = nexpr.subs(sigma**2, 2 * T).subs(x, h)
-    print(nexpr**3)
     pb_delta = (pexpr - bexpr).simplify()
     print(f"{pb_delta=}")
+    # sigma = 6T = 6 tau^-2
+    # 1/T = tau^2
+
+    nexpr = 1 / ((sigma**2 * 2 * pi) ** (1 / 2)) * exp(-1 / 2 * x**2 / sigma**2)
+    nexpr = nexpr.subs(sigma**2, 2 * T).subs(x, h)
+    print(nexpr**3)
     pn_delta = (pexpr / nexpr**3).simplify().expand()
     print(f"{pn_delta=}")
     # 1/tau2 == n*l2 / 6