From 52db21e3370295735f9568d3105db83b65a0b44c Mon Sep 17 00:00:00 2001 From: jupfi Date: Wed, 23 Aug 2023 14:51:57 +0200 Subject: [PATCH] Blochsimulator works and produces same results as matlab script. --- .gitignore | 3 +- src/nqr_blochsimulator/classes/pulse.py | 47 +++ src/nqr_blochsimulator/classes/sample.py | 13 +- src/nqr_blochsimulator/classes/simulation.py | 400 +++++++++++-------- tests/simulation.py | 70 +++- 5 files changed, 349 insertions(+), 184 deletions(-) create mode 100644 src/nqr_blochsimulator/classes/pulse.py diff --git a/.gitignore b/.gitignore index eba74f4..4ea05a1 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,2 @@ -venv/ \ No newline at end of file +venv/ +__pycache__/ \ No newline at end of file diff --git a/src/nqr_blochsimulator/classes/pulse.py b/src/nqr_blochsimulator/classes/pulse.py new file mode 100644 index 0000000..e6fc479 --- /dev/null +++ b/src/nqr_blochsimulator/classes/pulse.py @@ -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 diff --git a/src/nqr_blochsimulator/classes/sample.py b/src/nqr_blochsimulator/classes/sample.py index 57844ad..c3e61eb 100644 --- a/src/nqr_blochsimulator/classes/sample.py +++ b/src/nqr_blochsimulator/classes/sample.py @@ -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 diff --git a/src/nqr_blochsimulator/classes/simulation.py b/src/nqr_blochsimulator/classes/simulation.py index 00fda97..5a7df26 100644 --- a/src/nqr_blochsimulator/classes/simulation.py +++ b/src/nqr_blochsimulator/classes/simulation.py @@ -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 diff --git a/tests/simulation.py b/tests/simulation.py index 8db9ad2..ad25b2a 100644 --- a/tests/simulation.py +++ b/tests/simulation.py @@ -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) - \ No newline at end of file + 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() \ No newline at end of file