mirror of
https://github.com/nqrduck/nqr-blochsimulator.git
synced 2024-09-29 08:30:35 +00:00
Compare commits
3 commits
f845135dc2
...
83a4889748
Author | SHA1 | Date | |
---|---|---|---|
|
83a4889748 | ||
|
8182de30fb | ||
|
6643d8833c |
3 changed files with 113 additions and 23 deletions
|
@ -0,0 +1,7 @@
|
||||||
|
"""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"]
|
|
@ -1,5 +1,7 @@
|
||||||
from math import pi, sqrt
|
import logging
|
||||||
|
import numpy as np
|
||||||
|
|
||||||
|
logger = logging.getLogger(__name__)
|
||||||
|
|
||||||
class Sample:
|
class Sample:
|
||||||
"""
|
"""
|
||||||
|
@ -10,18 +12,18 @@ class Sample:
|
||||||
|
|
||||||
def __init__(
|
def __init__(
|
||||||
self,
|
self,
|
||||||
name,
|
name : str,
|
||||||
density,
|
density : float,
|
||||||
molar_mass,
|
molar_mass : float,
|
||||||
resonant_frequency,
|
resonant_frequency : float,
|
||||||
gamma,
|
gamma : float,
|
||||||
nuclear_spin,
|
nuclear_spin : float,
|
||||||
spin_factor,
|
spin_factor : float,
|
||||||
powder_factor,
|
powder_factor : float,
|
||||||
filling_factor,
|
filling_factor : float,
|
||||||
T1,
|
T1 : float,
|
||||||
T2,
|
T2 : float,
|
||||||
T2_star,
|
T2_star : float,
|
||||||
atom_density=None,
|
atom_density=None,
|
||||||
sample_volume=None,
|
sample_volume=None,
|
||||||
sample_length=None,
|
sample_length=None,
|
||||||
|
@ -45,7 +47,7 @@ class Sample:
|
||||||
nuclear_spin : float
|
nuclear_spin : float
|
||||||
The nuclear spin quantum number of the sample.
|
The nuclear spin quantum number of the sample.
|
||||||
spin_factor : float
|
spin_factor : float
|
||||||
The spin factor of the sample.
|
The spin transition factor of the sample.
|
||||||
powder_factor : float
|
powder_factor : float
|
||||||
The powder factor of the sample.
|
The powder factor of the sample.
|
||||||
filling_factor : float
|
filling_factor : float
|
||||||
|
@ -100,3 +102,69 @@ class Sample:
|
||||||
)
|
)
|
||||||
else:
|
else:
|
||||||
self.atoms = self.avogadro * self.density / self.molar_mass
|
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)
|
||||||
|
|
|
@ -1,10 +1,13 @@
|
||||||
import unittest
|
import unittest
|
||||||
import numpy as np
|
import numpy as np
|
||||||
|
import logging
|
||||||
import matplotlib.pyplot as plt
|
import matplotlib.pyplot as plt
|
||||||
from nqr_blochsimulator.classes.sample import Sample
|
from nqr_blochsimulator.classes.sample import Sample
|
||||||
from nqr_blochsimulator.classes.simulation import Simulation
|
from nqr_blochsimulator.classes.simulation import Simulation
|
||||||
from nqr_blochsimulator.classes.pulse import PulseArray
|
from nqr_blochsimulator.classes.pulse import PulseArray
|
||||||
|
|
||||||
|
logger = logging.getLogger(__name__)
|
||||||
|
logging.basicConfig(level=logging.INFO)
|
||||||
|
|
||||||
class TestSimulation(unittest.TestCase):
|
class TestSimulation(unittest.TestCase):
|
||||||
def setUp(self):
|
def setUp(self):
|
||||||
|
@ -13,14 +16,14 @@ class TestSimulation(unittest.TestCase):
|
||||||
density=1.585e6, # g/m^3
|
density=1.585e6, # g/m^3
|
||||||
molar_mass=440.3, # g/mol
|
molar_mass=440.3, # g/mol
|
||||||
resonant_frequency=83.56e6, # Hz
|
resonant_frequency=83.56e6, # Hz
|
||||||
gamma=4.342e7, # Hz/T
|
gamma=43.42, # MHz/T
|
||||||
nuclear_spin=9 / 2,
|
nuclear_spin= 9 / 2,
|
||||||
spin_factor=2,
|
spin_factor=1.94
|
||||||
powder_factor=0.75,
|
powder_factor=0.75,
|
||||||
filling_factor=0.7,
|
filling_factor=0.7,
|
||||||
T1=83.5e-5, # s
|
T1=83.5, # µs
|
||||||
T2=396e-6, # s
|
T2=396, # µs
|
||||||
T2_star=50e-6, # s
|
T2_star=50, # µs
|
||||||
)
|
)
|
||||||
|
|
||||||
simulation_length = 300e-6
|
simulation_length = 300e-6
|
||||||
|
@ -45,8 +48,8 @@ class TestSimulation(unittest.TestCase):
|
||||||
initial_magnetization=1,
|
initial_magnetization=1,
|
||||||
gradient=1,
|
gradient=1,
|
||||||
noise=0.5,
|
noise=0.5,
|
||||||
length_coil=6e-3,
|
length_coil=6, # mm
|
||||||
diameter_coil=3e-3,
|
diameter_coil=3, # mm
|
||||||
number_turns=9,
|
number_turns=9,
|
||||||
q_factor_transmit=100,
|
q_factor_transmit=100,
|
||||||
q_factor_receive=100,
|
q_factor_receive=100,
|
||||||
|
@ -57,7 +60,7 @@ class TestSimulation(unittest.TestCase):
|
||||||
temperature=300,
|
temperature=300,
|
||||||
loss_TX=12,
|
loss_TX=12,
|
||||||
loss_RX=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):
|
def test_simulation(self):
|
||||||
|
@ -69,3 +72,15 @@ class TestSimulation(unittest.TestCase):
|
||||||
plt.ylabel("Magnetization (a.u.)")
|
plt.ylabel("Magnetization (a.u.)")
|
||||||
plt.title("FID of BiPh3")
|
plt.title("FID of BiPh3")
|
||||||
plt.show()
|
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()
|
||||||
|
|
Loading…
Reference in a new issue