From 5b2b3645d135a97eee8d0dd484e0faf16c418d99 Mon Sep 17 00:00:00 2001 From: jupfi Date: Wed, 23 Aug 2023 16:47:41 +0200 Subject: [PATCH] Formatting. --- src/nqr_blochsimulator/classes/pulse.py | 7 +- src/nqr_blochsimulator/classes/sample.py | 2 +- src/nqr_blochsimulator/classes/simulation.py | 211 ++++++++++++------- tests/simulation.py | 37 ++-- 4 files changed, 163 insertions(+), 94 deletions(-) diff --git a/src/nqr_blochsimulator/classes/pulse.py b/src/nqr_blochsimulator/classes/pulse.py index e6fc479..48a3d5a 100644 --- a/src/nqr_blochsimulator/classes/pulse.py +++ b/src/nqr_blochsimulator/classes/pulse.py @@ -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 diff --git a/src/nqr_blochsimulator/classes/sample.py b/src/nqr_blochsimulator/classes/sample.py index c3e61eb..91041ab 100644 --- a/src/nqr_blochsimulator/classes/sample.py +++ b/src/nqr_blochsimulator/classes/sample.py @@ -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. diff --git a/src/nqr_blochsimulator/classes/simulation.py b/src/nqr_blochsimulator/classes/simulation.py index a685499..770cdec 100644 --- a/src/nqr_blochsimulator/classes/simulation.py +++ b/src/nqr_blochsimulator/classes/simulation.py @@ -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 diff --git a/tests/simulation.py b/tests/simulation.py index 446c141..6ca9223 100644 --- a/tests/simulation.py +++ b/tests/simulation.py @@ -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() \ No newline at end of file + plt.show()