From f3b48099746518d2396cfa1dafb988610c678625 Mon Sep 17 00:00:00 2001 From: lmantill <72053628+lmantill@users.noreply.github.com> Date: Thu, 25 Jan 2024 08:40:11 +0000 Subject: [PATCH 1/2] Finish ModMARY experiment and example --- examples/mod_mary.py | 103 +++++++++++++++++++++++++++++++++++++++ examples/sctep.py | 20 ++++---- radicalpy/experiments.py | 37 ++++++++++++++ radicalpy/utils.py | 48 ++++++++++++++++++ 4 files changed, 197 insertions(+), 11 deletions(-) create mode 100644 examples/mod_mary.py diff --git a/examples/mod_mary.py b/examples/mod_mary.py new file mode 100644 index 0000000..baf43d7 --- /dev/null +++ b/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() diff --git a/examples/sctep.py b/examples/sctep.py index 2eb7bdb..e5c036d 100644 --- a/examples/sctep.py +++ b/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__": diff --git a/radicalpy/experiments.py b/radicalpy/experiments.py index a3a56cc..7cf40f2 100644 --- a/radicalpy/experiments.py +++ b/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( diff --git a/radicalpy/utils.py b/radicalpy/utils.py index 0692399..ce244b1 100644 --- a/radicalpy/utils.py +++ b/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,23 @@ 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. From 836da4a6ba82ff0e0f2a3782d78dbbd57a98bc1b Mon Sep 17 00:00:00 2001 From: lmantill <72053628+lmantill@users.noreply.github.com> Date: Thu, 25 Jan 2024 08:41:06 +0000 Subject: [PATCH 2/2] Add formatting --- radicalpy/kinetics.py | 8 ++++++-- radicalpy/utils.py | 4 +++- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/radicalpy/kinetics.py b/radicalpy/kinetics.py index ea02f8f..4f8d03d 100644 --- a/radicalpy/kinetics.py +++ b/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): diff --git a/radicalpy/utils.py b/radicalpy/utils.py index ce244b1..e1f0260 100644 --- a/radicalpy/utils.py +++ b/radicalpy/utils.py @@ -295,7 +295,9 @@ 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): +def reference_signal( + timeconstant: np.ndarray, harmonic: float, theta: float, frequency: float +): """Modulated MARY reference signal. Used for brute force modulated MARY simulations.