nqrduck/nqr-blochsimulator
1
0
Fork 0
mirror of https://github.com/nqrduck/nqr-blochsimulator.git synced 2024-06-21 05:03:20 +00:00
Code Issues Packages Projects Releases Wiki Activity

Formatting.

Browse source
This commit is contained in:
jupfi 2023-08-23 16:47:41 +02:00
parent 8ae00f8b98
commit 5b2b3645d1
4 changed files with 163 additions and 94 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

7
src/nqr_blochsimulator/classes/pulse.py
View file

@ -1,5 +1,6 @@
import numpy as np
class PulseArray:
"""A class to represent a pulsearray for a NQR sequence."""
@ -23,7 +24,7 @@ class PulseArray:
def get_real_pulsepower(self) -> np.array:
"""Returns the real part of the pulse power."""
return self.pulseamplitude * np.cos(self.pulsephase)
def get_imag_pulsepower(self) -> np.array:
"""Returns the imaginary part of the pulse power."""
return self.pulseamplitude * np.sin(self.pulsephase)
@ -32,7 +33,7 @@ class PulseArray:
def pulseamplitude(self) -> np.array:
"""Amplitude of the pulse."""
return self._pulseamplitude
@pulseamplitude.setter
def pulseamplitude(self, pulseamplitude):
self._pulseamplitude = pulseamplitude
@ -41,7 +42,7 @@ class PulseArray:
def pulsephase(self) -> np.array:
"""Phase of the pulse."""
return self._pulsephase
@pulsephase.setter
def pulsephase(self, pulsephase):
self._pulsephase = pulsephase

2
src/nqr_blochsimulator/classes/sample.py
View file

@ -25,7 +25,7 @@ class Sample:
atom_density=None,
sample_volume=None,
sample_length=None,
sample_diameter=None
sample_diameter=None,
):
"""
Constructs all the necessary attributes for the sample object.

211
src/nqr_blochsimulator/classes/simulation.py
View file

@ -9,27 +9,27 @@ logger = logging.getLogger(__name__)
logger.setLevel(logging.DEBUG)
logger.addHandler(logging.StreamHandler())
class Simulation:
"""Class for the simulation of the Bloch equations."""
def __init__(
self,
sample : Sample,
number_isochromats : int,
initial_magnetization : float,
gradient : float,
noise : float,
length_coil : float,
diameter_coil : float,
number_turns : float,
power_amplifier_power : float,
pulse : PulseArray,
sample: Sample,
number_isochromats: int,
initial_magnetization: float,
gradient: float,
noise: float,
length_coil: float,
diameter_coil: float,
number_turns: float,
power_amplifier_power: float,
pulse: PulseArray,
averages: int,
gain: float,
temperature: float,
loss_TX: float = 0,
loss_RX: float = 0,
) -> None:
"""
Constructs all the necessary attributes for the simulation object.
@ -66,7 +66,7 @@ class Simulation:
The loss of the transmitter in dB.
loss_RX:
The loss of the receiver in dB.
"""
self.sample = sample
self.number_isochromats = number_isochromats
@ -87,11 +87,13 @@ class Simulation:
def simulate(self):
reference_voltage = self.calculate_reference_voltage()
B1 = self.calc_B1() * 1e3 # I think this is multiplied by 1e3 because everything is in mT
B1 = 17.3 # Something might be wrong with the calculation of the B1 field. This has to be checked.
self.sample.gamma = self.sample.gamma * 1e-6 # We need our gamma in MHz / T
self.sample.T1 = self.sample.T1 * 1e3 # We need our T1 in ms
self.sample.T2 = self.sample.T2 * 1e3 # We need our T2 in ms
B1 = (
self.calc_B1() * 1e3
) # I think this is multiplied by 1e3 because everything is in mT
B1 = 17.3 # Something might be wrong with the calculation of the B1 field. This has to be checked.
self.sample.gamma = self.sample.gamma * 1e-6 # We need our gamma in MHz / T
self.sample.T1 = self.sample.T1 * 1e3 # We need our T1 in ms
self.sample.T2 = self.sample.T2 * 1e3 # We need our T2 in ms
# Calculate the x distribution of the isochromats
xdis = self.calc_xdis()
@ -104,18 +106,22 @@ class Simulation:
imag_pulsepower = imag_pulsepower * (1 - 10 ** (-self.loss_TX / 20))
# Calculate the magnetization
M_sy1 = self.bloch_symmetric_strang_splitting(B1, xdis, real_pulsepower, imag_pulsepower)
M_sy1 = self.bloch_symmetric_strang_splitting(
B1, xdis, real_pulsepower, imag_pulsepower
)
# Z-Component
Mlong = np.squeeze(M_sy1[2,:,:]) # Indices start at 0 in Python
Mlong = np.squeeze(M_sy1[2, :, :]) # Indices start at 0 in Python
Mlong_avg = np.mean(Mlong, axis=0)
Mlong_avg = np.delete(Mlong_avg, -1) # Remove the last element
# XY-Component
Mtrans = np.squeeze(M_sy1[0,:,:] + 1j*M_sy1[1,:,:]) # Indices start at 0 in Python
Mtrans = np.squeeze(
M_sy1[0, :, :] + 1j * M_sy1[1, :, :]
) # Indices start at 0 in Python
Mtrans_avg = np.mean(Mtrans, axis=0)
Mtrans_avg = np.delete(Mtrans_avg, -1) # Remove the last element
# Scale the signal according to the reference voltage, averages and gain
timedomain_signal = Mtrans_avg * reference_voltage * 1e6
@ -129,8 +135,9 @@ class Simulation:
return timedomain_signal
def bloch_symmetric_strang_splitting(self, B1, xdis, real_pulsepower, imag_pulsepower, relax = 1):
def bloch_symmetric_strang_splitting(
self, B1, xdis, real_pulsepower, imag_pulsepower, relax=1
):
"""This method simulates the Bloch equations using the symmetric strang splitting method.
Parameters
@ -149,15 +156,17 @@ class Simulation:
Nx = self.number_isochromats
Nu = real_pulsepower.shape[0]
M0 = np.array([np.zeros(Nx), np.zeros(Nx), np.ones(Nx)])
dt = self.pulse.dwell_time * 1e3 # We need our dwell time in ms
dt = self.pulse.dwell_time * 1e3 # We need our dwell time in ms
w = np.ones((Nu, 1)) * self.gradient
w = np.ones((Nu, 1)) * self.gradient
# Bloch simulation in magnetization domain
gadt = self.sample.gamma * dt /2
B1 = np.tile((gadt * (real_pulsepower - 1j * imag_pulsepower) * B1).reshape(-1, 1), Nx)
gadt = self.sample.gamma * dt / 2
B1 = np.tile(
(gadt * (real_pulsepower - 1j * imag_pulsepower) * B1).reshape(-1, 1), Nx
)
K = gadt * xdis * w * self.gradient
phi = -np.sqrt(np.abs(B1) ** 2 + K ** 2)
phi = -np.sqrt(np.abs(B1) ** 2 + K**2)
cs = np.cos(phi)
si = np.sin(phi)
@ -177,44 +186,73 @@ class Simulation:
Bd8 = n3 * n2 * (1 - cs) + n1 * si
Bd9 = n3 * n3 * (1 - cs) + cs
M = np.zeros((3, Nx, Nu+1))
M = np.zeros((3, Nx, Nu + 1))
M[:, :, 0] = M0
Mt = M0
D = np.diag([np.exp(-1 / self.sample.T2 * relax * dt), np.exp(-1 / self.sample.T2 * relax * dt), np.exp(-1 / self.sample.T1 * relax * dt)])
b = np.array([0, 0, self.initial_magnetization]) - np.array([0, 0, self.initial_magnetization * np.exp(-1 / self.sample.T1 * relax * dt)])
for n in range(Nu): # time loop
D = np.diag(
[
np.exp(-1 / self.sample.T2 * relax * dt),
np.exp(-1 / self.sample.T2 * relax * dt),
np.exp(-1 / self.sample.T1 * relax * dt),
]
)
b = np.array([0, 0, self.initial_magnetization]) - np.array(
[
0,
0,
self.initial_magnetization * np.exp(-1 / self.sample.T1 * relax * dt),
]
)
for n in range(Nu): # time loop
Mrot = np.zeros((3, Nx))
Mrot[0,:] = Bd1.T[:,n]*Mt[0,:] + Bd2.T[:,n]*Mt[1,:] + Bd3.T[:,n]*Mt[2,:]
Mrot[1,:] = Bd4.T[:,n]*Mt[0,:] + Bd5.T[:,n]*Mt[1,:] + Bd6.T[:,n]*Mt[2,:]
Mrot[2,:] = Bd7.T[:,n]*Mt[0,:] + Bd8.T[:,n]*Mt[1,:] + Bd9.T[:,n]*Mt[2,:]
Mrot[0, :] = (
Bd1.T[:, n] * Mt[0, :] + Bd2.T[:, n] * Mt[1, :] + Bd3.T[:, n] * Mt[2, :]
)
Mrot[1, :] = (
Bd4.T[:, n] * Mt[0, :] + Bd5.T[:, n] * Mt[1, :] + Bd6.T[:, n] * Mt[2, :]
)
Mrot[2, :] = (
Bd7.T[:, n] * Mt[0, :] + Bd8.T[:, n] * Mt[1, :] + Bd9.T[:, n] * Mt[2, :]
)
Mt = np.dot(D, Mrot) + np.tile(b, (Nx, 1)).T
Mrot[0,:] = Bd1.T[:,n]*Mt[0,:] + Bd2.T[:,n]*Mt[1,:] + Bd3.T[:,n]*Mt[2,:]
Mrot[1,:] = Bd4.T[:,n]*Mt[0,:] + Bd5.T[:,n]*Mt[1,:] + Bd6.T[:,n]*Mt[2,:]
Mrot[2,:] = Bd7.T[:,n]*Mt[0,:] + Bd8.T[:,n]*Mt[1,:] + Bd9.T[:,n]*Mt[2,:]
Mrot[0, :] = (
Bd1.T[:, n] * Mt[0, :] + Bd2.T[:, n] * Mt[1, :] + Bd3.T[:, n] * Mt[2, :]
)
Mrot[1, :] = (
Bd4.T[:, n] * Mt[0, :] + Bd5.T[:, n] * Mt[1, :] + Bd6.T[:, n] * Mt[2, :]
)
Mrot[2, :] = (
Bd7.T[:, n] * Mt[0, :] + Bd8.T[:, n] * Mt[1, :] + Bd9.T[:, n] * Mt[2, :]
)
Mt = Mrot
M[:, :,n+1] = Mrot
M[:, :, n + 1] = Mrot
return M
def calc_B1(self) -> float:
"""This method calculates the B1 field of our solenoid coil based on the coil parameters and the power amplifier power.
Returns
-------
B1 : float
The B1 field of the solenoid coil in T."""
B1 = np.sqrt(2 * self.power_amplifier_power / 50) * np.pi * 4e-7 * self.number_turns / self.length_coil
B1 = (
np.sqrt(2 * self.power_amplifier_power / 50)
* np.pi
* 4e-7
* self.number_turns
/ self.length_coil
)
return B1
def calc_xdis(self) -> np.array:
""" Calculates the x distribution of the isochromats.
"""Calculates the x distribution of the isochromats.
Returns
-------
xdis : np.array
@ -228,11 +266,13 @@ class Simulation:
foffr = Df / 2 * np.tan(np.pi * uu)
# xdis is a spatial function, but it is being repurposed here to convert through the gradient to a phase difference per time -> T2 dispersion of the isochromats
xdis = np.linspace(-1, 1, self.number_isochromats)
xdis = (foffr.T * 1e-6) / (self.sample.gamma / 2 / np.pi) / (self.gradient * 1e-3)
xdis = np.linspace(-1, 1, self.number_isochromats)
xdis = (
(foffr.T * 1e-6) / (self.sample.gamma / 2 / np.pi) / (self.gradient * 1e-3)
)
return xdis
def calculate_reference_voltage(self) -> float:
"""This calculates the reference voltage of the measurement setup for the sample at a certain temperature.
@ -241,29 +281,56 @@ class Simulation:
reference_voltage : float
The reference voltage of the measurement setup for the sample at a certain temperature in Volts.
"""
u = 4 * np.pi * 1e-7 # permeability of free space
u = 4 * np.pi * 1e-7 # permeability of free space
magnetization = self.sample.gamma * 2 * self.sample.atoms / (2 * self.sample.nuclear_spin +1) * h**2 * self.sample.resonant_frequency/ Boltzmann / self.temperature * self.sample.spin_factor
magnetization = (
self.sample.gamma
* 2
* self.sample.atoms
/ (2 * self.sample.nuclear_spin + 1)
* h**2
* self.sample.resonant_frequency
/ Boltzmann
/ self.temperature
* self.sample.spin_factor
)
coil_crossection = np.pi * (self.diameter_coil / 2) ** 2
reference_voltage = self.number_turns * coil_crossection * u * self.sample.gamma * 2 * self.sample.atoms / (2 * self.sample.nuclear_spin +1) * h**2 * self.sample.resonant_frequency **2 / Boltzmann / self.temperature * self.sample.spin_factor
reference_voltage = reference_voltage * self.sample.powder_factor * self.sample.filling_factor
reference_voltage = (
self.number_turns
* coil_crossection
* u
* self.sample.gamma
* 2
* self.sample.atoms
/ (2 * self.sample.nuclear_spin + 1)
* h**2
* self.sample.resonant_frequency**2
/ Boltzmann
/ self.temperature
* self.sample.spin_factor
)
reference_voltage = (
reference_voltage * self.sample.powder_factor * self.sample.filling_factor
)
return reference_voltage
def calculate_noise(self, timedomain_signal : np.array) -> np.array:
def calculate_noise(self, timedomain_signal: np.array) -> np.array:
"""Calculates the noise array that is added to the signal.
Parameters
----------
timedomain_signal : np.array
The time domain signal that is used for the simulation.
Returns
-------
noise_data : np.array
The noise array that is added to the signal."""
n_timedomain_points = timedomain_signal.shape[0]
noise_data = (self.noise * np.random.randn(self.averages, n_timedomain_points) + 1j * self.noise * np.random.randn(self.averages, n_timedomain_points))
noise_data = self.noise * np.random.randn(
self.averages, n_timedomain_points
) + 1j * self.noise * np.random.randn(self.averages, n_timedomain_points)
noise_data = np.sum(noise_data, axis=0) # Sum along the first axis
return noise_data
@ -271,7 +338,7 @@ class Simulation:
def sample(self) -> Sample:
"""Sample that is used for the simulation."""
return self._sample
@sample.setter
def sample(self, sample):
self._sample = sample
@ -280,7 +347,7 @@ class Simulation:
def number_isochromats(self) -> int:
"""Number of isochromats used for the simulation."""
return self._number_isochromats
@number_isochromats.setter
def number_isochromats(self, number_isochromats):
self._number_isochromats = number_isochromats
@ -289,7 +356,7 @@ class Simulation:
def initial_magnetization(self) -> float:
"""Initial magnetization of the sample."""
return self._initial_magnetization
@initial_magnetization.setter
def initial_magnetization(self, initial_magnetization):
self._initial_magnetization = initial_magnetization
@ -298,16 +365,16 @@ class Simulation:
def gradient(self) -> float:
"""Gradient of the magnetic field in mt/M."""
return self._gradient
@gradient.setter
def gradient(self, gradient):
self._gradient = gradient
@property
def noise(self) -> float:
""" RMS Noise of the measurement setup in Volts"""
"""RMS Noise of the measurement setup in Volts"""
return self._noise
@noise.setter
def noise(self, noise):
self._noise = noise
@ -316,7 +383,7 @@ class Simulation:
def length_coil(self) -> float:
"""Length of the coil in meters."""
return self._length_coil
@length_coil.setter
def length_coil(self, length_coil):
self._length_coil = length_coil
@ -325,7 +392,7 @@ class Simulation:
def diameter_coil(self) -> float:
"""Diameter of the coil in meters."""
return self._diameter_coil
@diameter_coil.setter
def diameter_coil(self, diameter_coil):
self._diameter_coil = diameter_coil
@ -334,7 +401,7 @@ class Simulation:
def number_turns(self) -> float:
"""Number of turns of the coil."""
return self._number_turns
@number_turns.setter
def number_turns(self, number_turns):
self._number_turns = number_turns
@ -343,7 +410,7 @@ class Simulation:
def power_amplifier_power(self) -> float:
"""Power of the power amplifier in Watts."""
return self._power_amplifier_power
@power_amplifier_power.setter
def power_amplifier_power(self, power_amplifier_power):
self._power_amplifier_power = power_amplifier_power
@ -352,7 +419,7 @@ class Simulation:
def pulse(self) -> PulseArray:
"""Pulse that is used for the simulation."""
return self._pulse
@pulse.setter
def pulse(self, pulse):
self._pulse = pulse
@ -361,7 +428,7 @@ class Simulation:
def averages(self) -> int:
"""Number of averages that are used for the simulation."""
return self._averages
@averages.setter
def averages(self, averages):
self._averages = averages
@ -370,7 +437,7 @@ class Simulation:
def gain(self) -> float:
"""Gain of the amplifier."""
return self._gain
@gain.setter
def gain(self, gain):
self._gain = gain
@ -379,7 +446,7 @@ class Simulation:
def temperature(self) -> float:
"""Temperature of the sample."""
return self._temperature
@temperature.setter
def temperature(self, temperature):
self._temperature = temperature

37
tests/simulation.py
View file

@ -5,37 +5,38 @@ from nqr_blochsimulator.classes.sample import Sample
from nqr_blochsimulator.classes.simulation import Simulation
from nqr_blochsimulator.classes.pulse import PulseArray
class TestSimulation(unittest.TestCase):
def setUp(self):
self.sample = Sample(
"BiPh3",
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,
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,
powder_factor=0.75,
filling_factor=0.7,
T1=83.5e-5, #s
T2=396e-6, #s
T2_star=50e-6, #s
T1=83.5e-5, # s
T2=396e-6, # s
T2_star=50e-6, # s
)
simulation_length = 300e-6
dwell_time = 1e-6
self.time_array = np.arange(0, simulation_length, dwell_time)
pulse_length = 3e-6
# Simple FID sequence with pulse length of 3µs
pulse_amplitude_array = np.zeros(int(simulation_length/dwell_time))
pulse_amplitude_array[:int(pulse_length/dwell_time)] = 1
pulse_phase_array = np.zeros(int(simulation_length/dwell_time))
pulse_amplitude_array = np.zeros(int(simulation_length / dwell_time))
pulse_amplitude_array[: int(pulse_length / dwell_time)] = 1
pulse_phase_array = np.zeros(int(simulation_length / dwell_time))
self.pulse = PulseArray(
pulseamplitude=pulse_amplitude_array,
pulsephase=pulse_phase_array,
dwell_time=dwell_time
dwell_time=dwell_time,
)
self.simulation = Simulation(
@ -48,12 +49,12 @@ class TestSimulation(unittest.TestCase):
diameter_coil=3e-3,
number_turns=9,
power_amplifier_power=500,
pulse = self.pulse,
averages = 1,
gain = 6000,
pulse=self.pulse,
averages=1,
gain=6000,
temperature=77,
loss_TX=12,
loss_RX=12
loss_RX=12,
)
def test_simulation(self):
@ -64,4 +65,4 @@ class TestSimulation(unittest.TestCase):
plt.xlabel("Time (µs)")
plt.ylabel("Magnetization (a.u.)")
plt.title("FID of BiPh3")
plt.show()
plt.show()