nqrduck/nqr-blochsimulator
1
0
Fork 0
mirror of https://github.com/nqrduck/nqr-blochsimulator.git synced 2024-11-13 05:41:57 +00:00
Code Issues Projects Releases Packages Wiki Activity

Blochsimulator works and produces same results as matlab script.

Browse source
This commit is contained in:
jupfi 2023-08-23 14:51:57 +02:00
parent a0470a05e3
commit 52db21e337
5 changed files with 349 additions and 184 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
.gitignore vendored
View file

@ -1 +1,2 @@
venv/
venv/
__pycache__/

47
src/nqr_blochsimulator/classes/pulse.py Normal file
View file

@ -0,0 +1,47 @@
import numpy as np
class PulseArray:
"""A class to represent a pulsearray for a NQR sequence."""
def __init__(self, pulseamplitude, pulsephase, dwell_time) -> None:
"""
Constructs all the necessary attributes for the pulsearray object.
Parameters
----------
pulseamplitude : float
The amplitude of the pulse.
pulsephase : float
The phase of the pulse.
dwell_time : float
The dwell time of the pulse.
"""
self.pulseamplitude = pulseamplitude
self.pulsephase = pulsephase
self.dwell_time = dwell_time
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)
@property
def pulseamplitude(self) -> np.array:
"""Amplitude of the pulse."""
return self._pulseamplitude
@pulseamplitude.setter
def pulseamplitude(self, pulseamplitude):
self._pulseamplitude = pulseamplitude
@property
def pulsephase(self) -> np.array:
"""Phase of the pulse."""
return self._pulsephase
@pulsephase.setter
def pulsephase(self, pulsephase):
self._pulsephase = pulsephase

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

@ -24,6 +24,8 @@ class Sample:
T2_star,
atom_density=None,
sample_volume=None,
sample_length=None,
sample_diameter=None
):
"""
Constructs all the necessary attributes for the sample object.
@ -58,6 +60,10 @@ class Sample:
The atom density of the sample (atoms per cm^3). By default None.
sample_volume : float, optional
The volume of the sample (m^3). By default None.
sample_length : float, optional
The length of the sample (m). By default None.
sample_diameter : float, optional
The diameter of the sample (m). By default None.
"""
self.name = name
self.density = density
@ -73,11 +79,14 @@ class Sample:
self.T2_star = T2_star
self.atom_density = atom_density
self.sample_volume = sample_volume
self.sample_length = sample_length
self.sample_diameter = sample_diameter
self.calculate_atoms()
def calculate_atoms(self):
"""
Calculate the number of atoms in the sample per volume unit.
Calculate the number of atoms in the sample per volume unit. This only works if the sample volume and atom density are provided.
Also the sample should be cylindrical.
If atom density and sample volume are provided, use these to calculate the number of atoms.
If not, use Avogadro's number, density, and molar mass to calculate the number of atoms.
@ -87,7 +96,7 @@ class Sample:
self.atom_density
* self.sample_volume
/ 1e-6
/ (self.sample_volume * 6 / 3)
/ (self.sample_volume * self.sample_length / self.sample_diameter)
)
else:
self.atoms = self.avogadro * self.density / self.molar_mass

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

@ -1,148 +1,114 @@
import numpy as np
from numpy import pi
import logging
from scipy import signal
from .sample import Sample
from .pulse import PulseArray
logger = logging.getLogger(__name__)
logger.setLevel(logging.DEBUG)
logger.addHandler(logging.StreamHandler())
class Simulation:
def __init__(self) -> None:
pass
"""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
):
"""
Constructs all the necessary attributes for the simulation object.
def blochsim(sim_points, sim_time, reference, isochrom, sample, pulse):
# PRE-SETTINGS
d = {"M0c": 1} # initial mag
NISO = 100
if isochrom > 0:
NISO = isochrom # number of isochromates
nsamples = (
sim_points # number of sample/rasterization points for the calculation
)
sim_length = sim_time # in s; Not larger than the repetition time!
modulation = "OFF" # select a optional modulation of the pulse ['OFF','SIN']
# Replace by the NWA power later
B1c_calc = np.sqrt(2 * 500 / 50) * pi * 4e-7 * 9 / 6e-3 # 17 od 8.5
d["B1c"] = 17.3e-3 # for Peak B1 T %12.5%14.3
# SAMPLE SETTINGS
d["T1"] = sample["T1"] # in s; T1, T2
d["T2"] = sample["T2"] #
T2STAR = sample["T2s"] # only used for some calculations
d["gamma"] = (
sample["gamma"] / (2 * pi) / 1e6
) # gamma in MHz/T eg 5e6 % sample.gamma in rad/(T s) eg 0.8
d["relax"] = 1 # Flag 1: with Relax, 0 without Relax
# Parameter preparation
# clear up some unit problems
# DO NOT CHANGE
d["B1c"] = d["B1c"] * 1e3
d["T1"] = d["T1"] * 1e3
d["T2"] = d["T2"] * 1e3
d["gamma"] = d["gamma"] * 2 * pi
d["Nx"] = NISO
d["M0"] = np.array(
[np.zeros(NISO), np.zeros(NISO), np.ones(NISO)]
) # initial magnetization
d["dt"] = sim_length / nsamples # time step width
d["dt"] = (
d["dt"] * 1e3
) # again unit correction. could be changed if necessary, but other time factors
# in later calculations would have to be changed too
# Pulse Designer
u = np.zeros((nsamples, 1))
v = np.zeros((nsamples, 1))
w = np.ones((nsamples, 1))
tt = (np.array(range(1, nsamples + 1)) * d["dt"] - d["dt"]).reshape(
-1, 1
) # time axis in ms
# PULSE TEMPLATES
pulse_dur_pow_pha = pulse
num_pulses, _ = pulse_dur_pow_pha.shape
# loop through every pulse
for count in range(num_pulses):
pulse_begin = pulse_dur_pow_pha[count, 0]
pulse_end = pulse_dur_pow_pha[count, 1]
pha = pulse_dur_pow_pha[count, 3] * (2 * pi / 360) # phase in rad
ind_begin = np.argmin(np.abs(tt * 1e-3 - pulse_begin)) # minValue is unused
ind_end = np.argmin(np.abs(tt * 1e-3 - pulse_end))
ind_end = ind_end - 1
u_pow, v_pow = np.pol2cart(
pha, pulse_dur_pow_pha[count, 2]
) # theta angle; rho abs
if modulation == "OFF":
u[ind_begin:ind_end, 0] = u_pow # set real pulse power
v[ind_begin:ind_end, 0] = v_pow # set imag pulse power
elif modulation == "SIN":
u[ind_begin:ind_end, 0] = u_pow * np.sin(
(pi * 1e-3 / (pulse_end - pulse_begin)) * tt[ind_begin:ind_end, 0]
)
v[ind_begin:ind_end, 0] = v_pow * np.sin(
(pi * 1e-3 / (pulse_end - pulse_begin)) * tt[ind_begin:ind_end, 0]
)
# Some sidenotes that can be ignored
for count in range(1):
d["G3"] = 1 # mT/m, fhwm of 2mm Gradient scaling
w = w * d["G3"] # Gradient in mT/m
# Isochromatic simulaten by modeling with Lorentz distribution
Df = 1 / pi / T2STAR # FWHF of Lorentzian in Hz
foffr = []
uu = np.random.rand(NISO, 1) - 0.5
foffr = Df / 2 * np.tan(pi * uu) # cauchy distributed frequency offset
d["xdis"] = np.linspace(-1, 1, NISO) # in m spatial resolution
d["xdis"] = (
np.array(foffr) * 1e-6 / (d["gamma"] / 2 / pi) / (d["G3"] * 1e-3)
) # Conversion factors: foffr from Hz/T to MHz/T as required for d.gamma/2/pi, conversion from Hz-Gamma to radian gamma, and gradient must be scaled from mT/m to T/m
# USE BLOCH EQUATIONS
# M_sy1 = bloch_symmetric_strang_splitting_vectorised(u, v, w, d) # This function would need to be defined or imported
# Z-Component
# Mlong = np.squeeze(M_sy1[2, :, :]) # Coordinates M: space components - location(isochromat) - time
# Mlong_avg = np.mean(Mlong, 1)
# Mlong_avg = Mlong_avg[:-1]
# siglong = np.abs(Mlong_avg)
# XY-Component
# Mtrans = np.squeeze(M_sy1[0, :, :] + 1j*M_sy1[1, :, :]) # Coordinates M: space components - location(isochromat) - time
# Mtrans_avg = np.mean(Mtrans, 1)
# Mtrans_avg = Mtrans_avg[:-1]
# sigtrans = Mtrans_avg * reference
# return sigtrans
def bloch_symmetric_strang_splitting_vectorised(u, v, w, d):
"""Vectorised version of bloch_symmetric_strang_splitting
Parameters
----------
u : array_like
Real part of pulse
v : array_like
Imaginary part of pulse
w : array_like
Gradient
d : dict
sample : Sample
The sample that is used for the simulation.
number_isochromats : int
The number of isochromats used for the simulation.
initial_magnetization : float
The initial magnetization of the sample.
gradient : float
The gradient of the magnetic field in mt/M.
noise : float
The RMS Noise of the measurement setup in Volts.
length_coil : float
The length of the coil in meters.
diameter_coil : float
The diameter of the coil in meters.
number_turns : float
The number of turns of the coil.
power_amplifier_power : float
The power of the power amplifier in Watts.
puslse: PulseArray
The pulse that is used for the simulation.
"""
xdis = d["xdis"]
Nx = d["Nx"]
Nu = len(u)
M0 = d["M0"]
dt = d["dt"]
self.sample = sample
self.number_isochromats = number_isochromats
self.initial_magnetization = initial_magnetization
self.gradient = gradient
self.noise = noise
self.length_coil = length_coil
self.diameter_coil = diameter_coil
self.number_turns = number_turns
self.power_amplifier_power = power_amplifier_power
self.pulse = pulse
gadt = d["gamma"] * dt / 2
B1 = np.tile(gadt * np.transpose(u - 1j * v) * d["B1c"], (Nx, 1))
K = gadt * xdis * np.transpose(w) * d["G3"]
phi = -np.sqrt(np.abs(B1) ** 2 + K**2)
def simulate(self):
B1 = self.calc_B1()
xdis = self.calc_xdis()
real_pulsepower = self.pulse.get_real_pulsepower()
imag_pulsepower = self.pulse.get_imag_pulsepower()
M_sy1 = self.bloch_symmetric_strang_splitting(B1, xdis, real_pulsepower, imag_pulsepower)
logger.debug("Shape of Msy1: %s", M_sy1.shape)
# Z-Component
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
siglong = np.abs(Mlong_avg)
# XY-Component
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
reference = 4.5502
sigtrans = Mtrans_avg * reference
return sigtrans
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
----------
B1 : float
The B1 field of the solenoid coil.
xdis : np.array
The x distribution of the isochromats.
"""
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
w = np.ones((Nu, 1)) * self.gradient
# Bloch simulation in magnetization domain
gadt = self.sample.gamma * dt /2 * 1e-3
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)
cs = np.cos(phi)
si = np.sin(phi)
@ -162,40 +128,148 @@ class Simulation:
Bd8 = n3 * n2 * (1 - cs) + n1 * si
Bd9 = n3 * n3 * (1 - cs) + cs
M = np.zeros((3, Nx, Nu))
M = np.zeros((3, Nx, Nu+1))
M[:, :, 0] = M0
Mt = M0
D = np.diag(
[
np.exp(-1 / d["T2"] * d["relax"] * dt),
np.exp(-1 / d["T2"] * d["relax"] * dt),
np.exp(-1 / d["T1"] * d["relax"] * dt),
]
)
b = np.array([0, 0, d["M0c"]]) - np.array(
[0, 0, d["M0c"] * np.exp(-1 / d["T1"] * d["relax"] * dt)]
)
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.array(
[
Bd1[:, n] * Mt[0, :] + Bd2[:, n] * Mt[1, :] + Bd3[:, n] * Mt[2, :],
Bd4[:, n] * Mt[0, :] + Bd5[:, n] * Mt[1, :] + Bd6[:, n] * Mt[2, :],
Bd7[:, n] * Mt[0, :] + Bd8[:, n] * Mt[1, :] + Bd9[:, n] * Mt[2, :],
]
)
for n in range(Nu): # time loop
Mt = np.dot(D, Mrot) + np.tile(b, (Nx, 1)).transpose()
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 = np.array(
[
Bd1[:, n] * Mt[0, :] + Bd2[:, n] * Mt[1, :] + Bd3[:, n] * Mt[2, :],
Bd4[:, n] * Mt[0, :] + Bd5[:, n] * Mt[1, :] + Bd6[:, n] * Mt[2, :],
Bd7[:, n] * Mt[0, :] + Bd8[:, n] * Mt[1, :] + Bd9[:, 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,:]
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."""
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.
Returns
-------
xdis : np.array
The x distribution of the isochromats.
"""
# Df is the Full Width at Half Maximum (FWHM) of Lorentzian in Hz
Df = 1 / np.pi / self.sample.T2_star
logger.debug("Df: %s", Df)
# Randomly generating frequency offset using Cauchy distribution
uu = np.random.rand(self.number_isochromats, 1) - 0.5
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) / (self.sample.gamma / (2 * np.pi)) / (self.gradient * 1e-3)
return xdis
@property
def sample(self) -> Sample:
"""Sample that is used for the simulation."""
return self._sample
@sample.setter
def sample(self, sample):
self._sample = sample
@property
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
@property
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
@property
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"""
return self._noise
@noise.setter
def noise(self, noise):
self._noise = noise
@property
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
@property
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
@property
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
@property
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
@property
def pulse(self) -> PulseArray:
"""Pulse that is used for the simulation."""
return self._pulse
@pulse.setter
def pulse(self, pulse):
self._pulse = pulse

70
tests/simulation.py
View file

@ -1,25 +1,59 @@
import unittest
import numpy as np
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
class TestSimulation(unittest.TestCase):
def setUp(self):
self.sample = Sample(
"Ammonium nitrate",
1720,
80.0433
* 1e-3
/ Simulation.avogadro, # molar mass in kg/mol
1.945e6,
2 * 3.436e8,
1.5,
0.5,
1,
0.1,
0.1,
0.1,
0.1,
)
self.simulation = Simulation(self.sample, 1e-6, 1e-6, 1e-6, 1e-6, 1e-6)
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,
spin_factor=2,
powder_factor=1,
filling_factor=0.7,
T1=82.6e-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))
self.pulse = PulseArray(
pulseamplitude=pulse_amplitude_array,
pulsephase=pulse_phase_array,
dwell_time=dwell_time
)
self.simulation = Simulation(
sample=self.sample,
number_isochromats=1000,
initial_magnetization=1,
gradient=1,
noise=0,
length_coil=6e-3,
diameter_coil=3e-3,
number_turns=9,
power_amplifier_power=500,
pulse = self.pulse
)
def test_simulation(self):
M = self.simulation.simulate()
# Plotting the results
plt.plot(self.time_array, abs(M))
plt.show()