src/nqr_blochsimulator/__init__.py

@@ -1,7 +0,0 @@
-"""The nqr_blochsimulator package contains the classes necessary to simulate NQR experiments using the Bloch equations."""
-
-from .classes.sample import Sample
-from .classes.simulation import Simulation
-from .classes.pulse import PulseArray
-
-__all__ = ["Sample", "Simulation", "PulseArray"]

src/nqr_blochsimulator/classes/sample.py

@@ -1,7 +1,5 @@
-import logging
-import numpy as np
+from math import pi, sqrt
 
-logger = logging.getLogger(__name__)
 
 class Sample:
     """
@@ -12,18 +10,18 @@ class Sample:
 
     def __init__(
         self,
-        name : str,
-        density : float,
-        molar_mass : float,
-        resonant_frequency : float,
-        gamma : float,
-        nuclear_spin : float,
-        spin_factor : float,
-        powder_factor : float,
-        filling_factor : float,
-        T1 : float,
-        T2 : float,
-        T2_star : float,
+        name,
+        density,
+        molar_mass,
+        resonant_frequency,
+        gamma,
+        nuclear_spin,
+        spin_factor,
+        powder_factor,
+        filling_factor,
+        T1,
+        T2,
+        T2_star,
         atom_density=None,
         sample_volume=None,
         sample_length=None,
@@ -47,7 +45,7 @@ class Sample:
             nuclear_spin : float
                 The nuclear spin quantum number of the sample.
             spin_factor : float
-                The spin transition factor of the sample.
+                The spin factor of the sample.
             powder_factor : float
                 The powder factor of the sample.
             filling_factor : float
@@ -102,69 +100,3 @@ 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)

tests/simulation.py

@@ -1,13 +1,10 @@
 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):
@@ -16,14 +13,14 @@ class TestSimulation(unittest.TestCase):
             density=1.585e6,  # g/m^3
             molar_mass=440.3,  # g/mol
             resonant_frequency=83.56e6,  # Hz
-            gamma=43.42,  # MHz/T
-            nuclear_spin= 9 / 2,
-            spin_factor=1.94
+            gamma=4.342e7,  # Hz/T
+            nuclear_spin=9 / 2,
+            spin_factor=2,
             powder_factor=0.75,
             filling_factor=0.7,
-            T1=83.5,  # µs
-            T2=396,  # µs
-            T2_star=50,  # µs
+            T1=83.5e-5,  # s
+            T2=396e-6,  # s
+            T2_star=50e-6,  # s
         )
 
         simulation_length = 300e-6
@@ -48,8 +45,8 @@ class TestSimulation(unittest.TestCase):
             initial_magnetization=1,
             gradient=1,
             noise=0.5,
-            length_coil=6,  # mm
-            diameter_coil=3,  # mm
+            length_coil=6e-3,
+            diameter_coil=3e-3,
             number_turns=9,
             q_factor_transmit=100,
             q_factor_receive=100,
@@ -60,7 +57,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):
@@ -72,15 +69,3 @@ 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 = 1
-
-        spin_factor = self.sample.calculate_spin_transition_factor(spin, transition)
-        logger.info("Spin factor: " + str(spin_factor))
-
-
-if __name__ == "__main__":
-    unittest.main()