Merge pull request #145 from Spin-Chemistry-Labs/144-add-mod_mary

144 add mod mary

examples/mod_mary.py

@@ -0,0 +1,103 @@
+#! /usr/bin/env python
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from radicalpy.experiments import modulated_mary_brute_force
+
+
+def main(
+    Bs = np.linspace(-5, 5, 500),
+    modulation_depths = [2, 1.5, 1, 0.5, 0.1],
+    modulation_frequency = 3,
+    time_constant = 0.3,
+    harmonics = [1, 2, 3, 4, 5],
+    lfe_magnitude = 0.02,
+):
+    S = modulated_mary_brute_force(
+        Bs,
+        modulation_depths,
+        modulation_frequency,
+        time_constant,
+        harmonics,
+        lfe_magnitude)
+
+
+    harmonic = 0
+
+    plt.clf()
+    plt.grid(False)
+    plt.axis("on")
+    plt.rc("axes", edgecolor="k")
+    for i, md in enumerate(modulation_depths):
+        plt.plot(Bs, S[harmonic, i, :], label=f"{md} G", linewidth=3)
+    plt.legend()
+    plt.xlabel(r"$B_0$ / G", size=14)
+    plt.ylabel(r"ModMARY signal / au", size=14)
+    plt.tick_params(labelsize=14)
+    path = __file__[:-3] + f"_{0}.png"
+    plt.savefig(path)
+
+    harmonic = 1
+
+    plt.clf()
+    plt.grid(False)
+    plt.axis("on")
+    plt.rc("axes", edgecolor="k")
+    for i, md in enumerate(modulation_depths):
+        plt.plot(Bs, S[harmonic, i, :], label=f"{md} G", linewidth=3)
+    plt.legend()
+    plt.xlabel(r"$B_0$ / G", size=14)
+    plt.ylabel(r"ModMARY signal / au", size=14)
+    plt.tick_params(labelsize=14)
+    path = __file__[:-3] + f"_{1}.png"
+    plt.savefig(path)
+
+    harmonic = 2
+
+    plt.clf()
+    plt.grid(False)
+    plt.axis("on")
+    plt.rc("axes", edgecolor="k")
+    for i, md in enumerate(modulation_depths):
+        plt.plot(Bs, S[harmonic, i, :], label=f"{md} G", linewidth=3)
+    plt.legend()
+    plt.xlabel(r"$B_0$ / G", size=14)
+    plt.ylabel(r"ModMARY signal / au", size=14)
+    plt.tick_params(labelsize=14)
+    path = __file__[:-3] + f"_{2}.png"
+    plt.savefig(path)
+
+    harmonic = 3
+
+    plt.clf()
+    plt.grid(False)
+    plt.axis("on")
+    plt.rc("axes", edgecolor="k")
+    for i, md in enumerate(modulation_depths):
+        plt.plot(Bs, S[harmonic, i, :], label=f"{md} G", linewidth=3)
+    plt.legend()
+    plt.xlabel(r"$B_0$ / G", size=14)
+    plt.ylabel(r"ModMARY signal / au", size=14)
+    plt.tick_params(labelsize=14)
+    path = __file__[:-3] + f"_{3}.png"
+    plt.savefig(path)
+
+    harmonic = 4
+
+    plt.clf()
+    plt.grid(False)
+    plt.axis("on")
+    plt.rc("axes", edgecolor="k")
+    for i, md in enumerate(modulation_depths):
+        plt.plot(Bs, S[harmonic, i, :], label=f"{md} G", linewidth=3)
+    plt.legend()
+    plt.xlabel(r"$B_0$ / G", size=14)
+    plt.ylabel(r"ModMARY signal / au", size=14)
+    plt.tick_params(labelsize=14)
+    path = __file__[:-3] + f"_{4}.png"
+    plt.savefig(path)
+
+
+if __name__ == "__main__":
+    main()

examples/sctep.py

@@ -35,21 +35,19 @@ def main(
     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.clf()
+    plt.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.xlabel("g$μ_B$$B_0$ / J", size=14)
+    plt.ylabel("MFE (%)", size=14)
+    plt.axvline(x=1.5, color="k", linestyle="--")
+    plt.axvline(x=3, color="k", linestyle="--")
     plt.tick_params(labelsize=14)
-    fig.set_size_inches(8, 5)
-    plt.show()
+    path = __file__[:-3] + f"_{0}.png"
+    plt.savefig(path)
 
 
 if __name__ == "__main__":

radicalpy/experiments.py

@@ -4,6 +4,43 @@
 from tqdm import tqdm
 
 from .simulation import HilbertIncoherentProcessBase, LiouvilleSimulation, State
+from .utils import mary_lorentzian, modulated_signal, reference_signal
+
+
+def modulated_mary_brute_force(
+    Bs: np.ndarray,
+    modulation_depths: list,
+    modulation_frequency: float,
+    time_constant: float,
+    harmonics: list,
+    lfe_magnitude: float,
+) -> np.ndarray:
+    t = np.linspace(-3 * time_constant, 0, 100)  # Simulate over 10 time constants
+    S = np.zeros([len(harmonics), len(modulation_depths), len(Bs)])
+    sa = 0
+
+    for i, h in enumerate(tqdm(harmonics)):
+        for j, md in enumerate(modulation_depths):
+            for k, B in enumerate(Bs):
+                for l in range(0, 20):  # Randomise phase
+                    theta = 2 * np.pi * np.random.rand()
+
+                    # Calculate the modulated signal
+                    ms = B + md * modulated_signal(t, theta, modulation_frequency)
+                    s = mary_lorentzian(ms, lfe_magnitude)
+                    s = s - np.mean(s)  # Signal (AC coupled)
+
+                    # Calculate the reference signal
+                    envelope = np.exp(t / time_constant) / time_constant  # Envelope
+                    r = (
+                        reference_signal(t, h, theta, modulation_frequency) * envelope
+                    )  # Reference * envelope
+
+                    # Calculate the MARY spectra
+                    sa = sa + np.trapz(t, s * r)  # Integrate
+                sa = sa
+                S[i, j, k] = sa * np.sqrt(2)  # RMS signal
+    return S
 
 
 def steady_state_mary(

radicalpy/kinetics.py

@@ -2,8 +2,12 @@
 
 import numpy as np
 
-from .simulation import (HilbertIncoherentProcessBase, LiouvilleIncoherentProcessBase,
-                         LiouvilleSimulation, State)
+from .simulation import (
+    HilbertIncoherentProcessBase,
+    LiouvilleIncoherentProcessBase,
+    LiouvilleSimulation,
+    State,
+)
 
 
 class HilbertKineticsBase(HilbertIncoherentProcessBase):

radicalpy/utils.py

@@ -110,6 +110,21 @@ def Lorentzian(B: np.ndarray, amplitude: float, Bhalf: float) -> np.ndarray:
     return (amplitude / Bhalf**2) - (amplitude / (B**2 + Bhalf**2))
 
 
+def mary_lorentzian(mod_signal: np.ndarray, lfe_magnitude: float):
+    """Lorentzian MARY spectral shape.
+
+    Used for brute force modulated MARY simulations.
+
+    Args:
+            mod_signal (np.ndarray): The modulated signal.
+            lfe_magnitude (float): The magnitude of your low field effect (G).
+
+    Returns:
+            np.ndarray: The modulated MARY signal.
+    """
+    return 1 / (1 + mod_signal**2) - lfe_magnitude / (0.1 + mod_signal**2)
+
+
 def MHz_in_angular_frequency(MHz: float) -> float:
     """Convert MHz into angular frequency.
 
@@ -146,6 +161,22 @@ def MHz_to_mT(MHz: float) -> float:
     return MHz / (1e-9 * -C.g_e * C.mu_B / C.h)
 
 
+def modulated_signal(timeconstant: np.ndarray, theta: float, frequency: float):
+    """Modulated MARY signal.
+
+    Used for brute force modulated MARY simulations.
+
+    Args:
+            timeconstant (np.ndarray): The modulation time constant (s).
+            theta (float): The modulation phase (rad).
+            frequency (float): The modulation frequency (Hz).
+
+    Returns:
+            np.ndarray: The modulated signal.
+    """
+    return np.cos(frequency * timeconstant * (2 * np.pi) + theta)
+
+
 def angular_frequency_in_MHz(ang_freq: float) -> float:
     """Convert angular frequency into MHz.
 
@@ -264,6 +295,25 @@ def mT_to_angular_frequency(mT: float) -> float:
     return mT * (C.mu_B / C.hbar * -C.g_e / 1e9)
 
 
+def reference_signal(
+    timeconstant: np.ndarray, harmonic: float, theta: float, frequency: float
+):
+    """Modulated MARY reference signal.
+
+    Used for brute force modulated MARY simulations.
+
+    Args:
+            timeconstant (np.ndarray): The modulation time constant (s).
+            harmonic (float): The harmonic of the modulation.
+            theta (float): The modulation phase (rad).
+            frequency (float): The modulation frequency (Hz).
+
+    Returns:
+            np.ndarray: The modulated reference signal.
+    """
+    return np.cos(harmonic * frequency * timeconstant * (2 * np.pi) + harmonic * theta)
+
+
 def spectral_density(omega: float, tau_c: float) -> float:
     """Frequency at which the motion of the particle exists.