diff --git a/src/nqr_blochsimulator/classes/sample.py b/src/nqr_blochsimulator/classes/sample.py index c869366..d4a2632 100644 --- a/src/nqr_blochsimulator/classes/sample.py +++ b/src/nqr_blochsimulator/classes/sample.py @@ -1,5 +1,7 @@ -from math import pi, sqrt +import logging +import numpy as np +logger = logging.getLogger(__name__) class Sample: """ @@ -10,18 +12,18 @@ class Sample: def __init__( self, - name, - density, - molar_mass, - resonant_frequency, - gamma, - nuclear_spin, - spin_factor, - powder_factor, - filling_factor, - T1, - T2, - T2_star, + name : str, + density : float, + molar_mass : float, + resonant_frequency : float, + gamma : float, + nuclear_spin : float, + spin_transition : int, + powder_factor : float, + filling_factor : float, + T1 : float, + T2 : float, + T2_star : float, atom_density=None, sample_volume=None, sample_length=None, @@ -44,8 +46,13 @@ class Sample: The gamma value of the sample in MHz/T. nuclear_spin : float The nuclear spin quantum number of the sample. - spin_factor : float - The spin factor of the sample. + spin_transition: int + The spin transition of the sample. + 0 is -1/2 -> 1/2 + 1 is 1/2 -> 3/2 + 2 is 3/2 -> 5/2 + 3 is 5/2 -> 7/2 + 4 is 7/2 -> 9/2 powder_factor : float The powder factor of the sample. filling_factor : float @@ -71,7 +78,8 @@ class Sample: self.resonant_frequency = resonant_frequency * 1e6 self.gamma = gamma * 1e6 self.nuclear_spin = nuclear_spin - self.spin_factor = spin_factor + self.spin_transition = spin_transition + self.spin_factor = self.calculate_spin_transition_factor(nuclear_spin, self.spin_transition) self.powder_factor = powder_factor self.filling_factor = filling_factor self.T1 = T1 * 1e-6 @@ -100,3 +108,69 @@ class Sample: ) else: self.atoms = self.avogadro * self.density / self.molar_mass + + def pauli_spin_matrices(self, spin): + """ + Generate the spin matrices for a given spin value. + + Parameters: + spin (float): The spin value, which can be a half-integer or integer. + + Returns: + tuple: A tuple containing the following elements: + Jx (np.ndarray): The x-component of the spin matrix. + Jy (np.ndarray): The y-component of the spin matrix. + Jz (np.ndarray): The z-component of the spin matrix. + J_minus (np.ndarray): The lowering operator matrix. + J_plus (np.ndarray): The raising operator matrix. + m (np.ndarray): The array of magnetic quantum numbers. + """ + + m = np.arange(spin, -spin-1, -1) + paulirowlength = int(spin * 2 + 1) + + pauli_z = np.diag(m) + pauli_plus = np.zeros((paulirowlength, paulirowlength)) + pauli_minus = np.zeros((paulirowlength, paulirowlength)) + + for row_index in range(paulirowlength - 1): + col_index = row_index + 1 + pauli_plus[row_index, col_index] = np.sqrt(spin * (spin + 1) - m[col_index] * (m[col_index] + 1)) + + for row_index in range(1, paulirowlength): + col_index = row_index - 1 + pauli_minus[row_index, col_index] = np.sqrt(spin * (spin + 1) - m[col_index] * (m[col_index] - 1)) + + Jx = 0.5 * (pauli_plus + pauli_minus) + Jy = -0.5j * (pauli_plus - pauli_minus) + Jz = pauli_z + + return Jx, Jy, Jz, pauli_minus, pauli_plus, m + + def calculate_spin_transition_factor(self, I, transition): + """ + Calculate the prefactor for the envisaged spin transition for a given nuclear spin. + + Parameters: + I (float): The nuclear spin value, which can be a half-integer or integer. + transition (int): The index of the transition. + The transition indices represent the shifts between magnetic quantum numbers m: + - 0 represents -1/2 --> 1/2 + - 1 represents 1/2 --> 3/2 + - 2 represents 3/2 --> 5/2 + (only valid transitions based on spin value I are allowed) + + Returns: + float: The prefactor for the envisaged spin transition. + """ + m_values = np.arange(I, -I-1, -1) + if transition < 0 or transition >= len(m_values) - 1: + raise ValueError(f"Invalid transition for spin {I}. Valid range is 0 to {len(m_values) - 2}") + + Jx, Jy, Jz, J_minus, J_plus, m = self.pauli_spin_matrices(I) + trindex = int(len(Jx) / 2 - transition) + spinfactor = Jx[trindex - 1, trindex] + + logger.debug(f"Spin transition factor for I={I}, transition={transition}: {np.real(spinfactor)}") + logger.info(f"Jx is {Jx}") + return np.real(spinfactor) diff --git a/tests/simulation.py b/tests/simulation.py index 499482f..bb5f873 100644 --- a/tests/simulation.py +++ b/tests/simulation.py @@ -1,10 +1,13 @@ import unittest import numpy as np +import logging import matplotlib.pyplot as plt from nqr_blochsimulator.classes.sample import Sample from nqr_blochsimulator.classes.simulation import Simulation from nqr_blochsimulator.classes.pulse import PulseArray +logger = logging.getLogger(__name__) +logging.basicConfig(level=logging.INFO) class TestSimulation(unittest.TestCase): def setUp(self): @@ -13,14 +16,14 @@ class TestSimulation(unittest.TestCase): density=1.585e6, # g/m^3 molar_mass=440.3, # g/mol resonant_frequency=83.56e6, # Hz - gamma=4.342e7, # Hz/T - nuclear_spin=9 / 2, - spin_factor=2, + gamma=43.42, # MHz/T + nuclear_spin= 7 / 2, + spin_transition=2, powder_factor=0.75, filling_factor=0.7, - T1=83.5e-5, # s - T2=396e-6, # s - T2_star=50e-6, # s + T1=83.5, # µs + T2=396, # µs + T2_star=50, # µs ) simulation_length = 300e-6 @@ -45,8 +48,8 @@ class TestSimulation(unittest.TestCase): initial_magnetization=1, gradient=1, noise=0.5, - length_coil=6e-3, - diameter_coil=3e-3, + length_coil=6, # mm + diameter_coil=3, # mm number_turns=9, q_factor_transmit=100, q_factor_receive=100, @@ -57,7 +60,7 @@ class TestSimulation(unittest.TestCase): temperature=300, loss_TX=12, loss_RX=12, - conversion_factor=2884 # This is for the LimeSDR based spectrometer + conversion_factor=2884, # This is for the LimeSDR based spectrometer ) def test_simulation(self): @@ -69,3 +72,15 @@ class TestSimulation(unittest.TestCase): plt.ylabel("Magnetization (a.u.)") plt.title("FID of BiPh3") plt.show() + + def test_spin_factor_calculation(self): + + spin = self.sample.nuclear_spin + transition = self.sample.spin_transition + + spin_factor = self.sample.calculate_spin_transition_factor(spin, transition) + logger.info("Spin factor: " + str(spin_factor)) + + +if __name__ == "__main__": + unittest.main()