From b68184b554c840de3cf890cfd73441dc2edfed4a Mon Sep 17 00:00:00 2001
From: Emil VATAI <emil.vatai@gmail.com>
Date: Sun, 21 Jan 2024 18:39:54 +0900
Subject: [PATCH] It runs (but wrong output)

---
 examples/sctep.py        | 25 +++++++++++++++++++++++--
 radicalpy/experiments.py | 29 ++++++++++-------------------
 2 files changed, 33 insertions(+), 21 deletions(-)

diff --git a/examples/sctep.py b/examples/sctep.py
index 43c0e10..0cd5178 100644
--- a/examples/sctep.py
+++ b/examples/sctep.py
@@ -1,5 +1,6 @@
 #! /usr/bin/env python
 
+import matplotlib.pyplot as plt
 import numpy as np
 
 from radicalpy.data import Isotope, Molecule, Nucleus, Triplet
@@ -25,16 +26,36 @@ def main(
         E=E,
     )
     print(H)
-    steady_state_mary(
+    rhos, Phi_s = steady_state_mary(
         sim,
-        obs_state=State.TP_SINGLET,
+        obs=State.TP_SINGLET,
         Bs=Bs,
         D=D,
         E=E,
         J=J,
+        theta=np.pi / 4,
+        phi=0,
         kinetics=[Haberkorn(ks, State.TP_SINGLET), HaberkornFree(kd)],
     )
+    rhos *= k0
+    Phi_s *= k0 * ks
+
+    MFE = ((np.abs(Phi_s) - np.abs(Phi_s[0])) / np.abs(Phi_s[0])) * 100
     # print(H)
+    fig = plt.figure(1)
+    ax = fig.add_axes([0, 0, 1, 1])
+    ax.set_facecolor("none")
+    ax.grid(False)
+    plt.axis("on")
+    plt.rc("axes", edgecolor="k")
+    plt.plot(Bs / J, MFE, linewidth=3, color="tab:red")
+    ax.set_xlabel("g$μ_B$$B_0$ / J", size=18)
+    ax.set_ylabel("MFE (%)", size=18)
+    ax.axvline(x=1.5, color="k", linestyle="--")
+    ax.axvline(x=3, color="k", linestyle="--")
+    plt.tick_params(labelsize=14)
+    fig.set_size_inches(8, 5)
+    plt.show()
 
 
 if __name__ == "__main__":
diff --git a/radicalpy/experiments.py b/radicalpy/experiments.py
index 1ba7dc4..e44f93b 100644
--- a/radicalpy/experiments.py
+++ b/radicalpy/experiments.py
@@ -9,34 +9,25 @@
 
 def steady_state_mary(
     sim: HilbertSimulation,
-    obs_state: State,
+    obs: State,
     Bs: np.ndarray,
     D: float,
     E: float,
     J: float,
+    theta: float,
+    phi: float,
     kinetics: list[HilbertIncoherentProcessBase] = [],
-    relaxations: list[HilbertIncoherentProcessBase] = [],
+    # relaxations: list[HilbertIncoherentProcessBase] = [],
 ) -> np.ndarray:
-    theta = np.pi / 4
-    phi = 0
     HZFS = sim.zero_field_splitting_hamiltonian(D, E)
     HJ = sim.exchange_hamiltonian(J)
-    Phi_s = np.zeros((len(Bs)), dtype=complex)
-    Ps = sim.projection_operator(State.TP_SINGLET)
+    rhos = np.zeros(shape=(len(Bs), *HJ.shape))
+    Q = sim.projection_operator(obs)
     for i, B in enumerate(tqdm(Bs)):
         HZ = sim.zeeman_hamiltonian_3d(B, theta, phi)
-        print(f"{HZ.shape=}")
-        print(f"{HZFS.shape=}")
-        print(f"{HJ.shape=}")
         H = HZ + HZFS + HJ
-        sim.apply_liouville_hamiltonian_modifiers(H, kinetics + relaxations)
-        # HL = ry.Hilbert2Liouville(H)
-        # L = 1j * H + sum(kinetics)
-        rho = np.linalg.solve(H, Ps)  # Density operator
+        sim.apply_liouville_hamiltonian_modifiers(H, kinetics)  # + relaxations)
+        rhos[i] = np.linalg.solve(H, Q)
+    Phi_s = sim.product_probability(obs, rhos)
 
-        print(f"{type(rho)=} {rho.shape=}")
-        print(f"{type(Ps)=} {Ps.shape=}")
-
-        Phi_s[i] = np.matmul(Ps.T, rho)
-
-    return rho, Phi_s
+    return rhos, Phi_s