2023-08-22 14:06:38 +00:00
import numpy as np
2023-08-23 12:51:57 +00:00
import logging
2023-08-23 14:18:22 +00:00
from scipy . constants import h , Boltzmann
2023-08-23 12:51:57 +00:00
from . sample import Sample
from . pulse import PulseArray
2023-08-22 14:06:38 +00:00
2023-08-23 12:51:57 +00:00
logger = logging . getLogger ( __name__ )
logger . setLevel ( logging . DEBUG )
logger . addHandler ( logging . StreamHandler ( ) )
2023-08-23 14:47:41 +00:00
2023-08-22 14:06:38 +00:00
class Simulation :
2023-08-23 12:51:57 +00:00
""" Class for the simulation of the Bloch equations. """
2023-08-23 14:47:41 +00:00
2023-08-23 12:51:57 +00:00
def __init__ (
self ,
2023-08-23 14:47:41 +00:00
sample : Sample ,
number_isochromats : int ,
initial_magnetization : float ,
gradient : float ,
noise : float ,
length_coil : float ,
diameter_coil : float ,
number_turns : float ,
2024-03-14 13:56:06 +00:00
q_factor_transmit : float ,
2023-12-14 11:39:44 +00:00
q_factor_receive : float ,
2023-08-23 14:47:41 +00:00
power_amplifier_power : float ,
pulse : PulseArray ,
2023-08-23 13:48:17 +00:00
averages : int ,
2023-08-23 14:18:22 +00:00
gain : float ,
2023-08-23 14:33:02 +00:00
temperature : float ,
loss_TX : float = 0 ,
loss_RX : float = 0 ,
2023-12-14 11:39:44 +00:00
conversion_factor : float = 1 ,
2023-08-23 13:48:17 +00:00
) - > None :
2023-08-23 12:51:57 +00:00
"""
Constructs all the necessary attributes for the simulation object .
Parameters
- - - - - - - - - -
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
2023-08-23 14:44:39 +00:00
The RMS Noise of the measurement setup in µVolts .
2023-08-23 12:51:57 +00:00
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 .
2024-03-02 19:58:07 +00:00
q_factor_transmit : float
The Q - factor of the transmit path of the probe coil .
2023-12-14 11:39:44 +00:00
q_factor_receive : float
The Q - factor of the receive path of the probe coil .
2023-08-23 12:51:57 +00:00
power_amplifier_power : float
The power of the power amplifier in Watts .
2023-08-23 13:48:17 +00:00
pulse : PulseArray
2023-08-23 12:51:57 +00:00
The pulse that is used for the simulation .
2023-08-23 13:48:17 +00:00
averages :
The number of averages that are used for the simulation .
gain :
The gain of the amplifier .
2023-08-23 14:18:22 +00:00
temperature :
The temperature of the sample in Kelvin .
2023-08-23 14:33:02 +00:00
loss_TX :
The loss of the transmitter in dB .
loss_RX :
The loss of the receiver in dB .
2023-12-14 11:39:44 +00:00
conversion_factor :
The conversion factor of the receiver in spectromter units / Volt .
2023-08-23 14:47:41 +00:00
2023-08-23 12:51:57 +00:00
"""
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
2024-03-02 19:58:07 +00:00
self . q_factor_transmit = q_factor_transmit
2023-12-14 11:39:44 +00:00
self . q_factor_receive = q_factor_receive
2023-08-23 12:51:57 +00:00
self . power_amplifier_power = power_amplifier_power
self . pulse = pulse
2023-08-23 13:48:17 +00:00
self . averages = averages
self . gain = gain
2023-08-23 14:18:22 +00:00
self . temperature = temperature
2023-08-23 14:33:02 +00:00
self . loss_TX = loss_TX
self . loss_RX = loss_RX
2023-12-14 11:39:44 +00:00
self . conversion_factor = conversion_factor
2023-08-23 12:51:57 +00:00
def simulate ( self ) :
2023-08-23 14:18:22 +00:00
reference_voltage = self . calculate_reference_voltage ( )
2023-08-23 14:44:39 +00:00
2023-08-23 14:47:41 +00:00
B1 = (
self . calc_B1 ( ) * 1e3
) # I think this is multiplied by 1e3 because everything is in mT
2023-08-24 11:44:30 +00:00
# B1 = 17.3 # Something might be wrong with the calculation of the B1 field. This has to be checked.
2023-08-23 14:47:41 +00:00
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
2023-08-23 13:48:17 +00:00
2023-08-23 14:33:02 +00:00
# Calculate the x distribution of the isochromats
2023-08-23 12:51:57 +00:00
xdis = self . calc_xdis ( )
real_pulsepower = self . pulse . get_real_pulsepower ( )
imag_pulsepower = self . pulse . get_imag_pulsepower ( )
2023-08-23 14:33:02 +00:00
# Calculate losses on the pulse
real_pulsepower = real_pulsepower * ( 1 - 10 * * ( - self . loss_TX / 20 ) )
imag_pulsepower = imag_pulsepower * ( 1 - 10 * * ( - self . loss_TX / 20 ) )
# Calculate the magnetization
2023-08-23 14:47:41 +00:00
M_sy1 = self . bloch_symmetric_strang_splitting (
B1 , xdis , real_pulsepower , imag_pulsepower
)
2023-08-22 14:06:38 +00:00
# Z-Component
2023-08-23 14:47:41 +00:00
Mlong = np . squeeze ( M_sy1 [ 2 , : , : ] ) # Indices start at 0 in Python
2023-08-23 12:51:57 +00:00
Mlong_avg = np . mean ( Mlong , axis = 0 )
Mlong_avg = np . delete ( Mlong_avg , - 1 ) # Remove the last element
2023-08-22 14:06:38 +00:00
# XY-Component
2023-08-23 14:47:41 +00:00
Mtrans = np . squeeze (
2024-01-10 11:26:39 +00:00
M_sy1 [ 1 , : , : ] + 1 j * M_sy1 [ 0 , : , : ]
2023-08-23 14:47:41 +00:00
) # Indices start at 0 in Python
2023-08-23 12:51:57 +00:00
Mtrans_avg = np . mean ( Mtrans , axis = 0 )
Mtrans_avg = np . delete ( Mtrans_avg , - 1 ) # Remove the last element
2023-08-23 14:47:41 +00:00
2023-12-14 12:56:00 +00:00
2023-08-23 14:33:02 +00:00
# Scale the signal according to the reference voltage, averages and gain
2023-12-14 12:56:00 +00:00
timedomain_signal = Mtrans_avg * reference_voltage
2023-08-23 14:33:02 +00:00
# Add the losses of the receiver - this should probably be done before the scaling
timedomain_signal = timedomain_signal * ( 1 - 10 * * ( - self . loss_RX / 20 ) )
2023-08-23 14:44:39 +00:00
# Add noise to the signal
noise_data = self . calculate_noise ( timedomain_signal )
2023-12-14 12:56:00 +00:00
timedomain_signal = ( timedomain_signal * self . averages * self . gain ) + ( noise_data * self . gain )
# print(abs(timedomain_signal))
2023-08-23 14:44:39 +00:00
2023-12-14 12:56:36 +00:00
timedomain_signal = timedomain_signal
2023-12-14 14:01:05 +00:00
return timedomain_signal * self . conversion_factor
2023-08-22 14:06:38 +00:00
2023-08-23 14:47:41 +00:00
def bloch_symmetric_strang_splitting (
self , B1 , xdis , real_pulsepower , imag_pulsepower , relax = 1
) :
2023-08-23 12:51:57 +00:00
""" This method simulates the Bloch equations using the symmetric strang splitting method.
2023-08-22 14:06:38 +00:00
Parameters
- - - - - - - - - -
2023-08-23 12:51:57 +00:00
B1 : float
The B1 field of the solenoid coil .
xdis : np . array
The x distribution of the isochromats .
2023-08-23 14:18:22 +00:00
real_pulsepower : np . array
The real part of the pulse power .
imag_pulsepower : np . array
The imaginary part of the pulse power .
relax : float
If relax = 1 , the relaxation is taken into account . If relax = 0 , the relaxation is not taken into account .
2023-08-22 14:06:38 +00:00
"""
2023-08-23 12:51:57 +00:00
Nx = self . number_isochromats
Nu = real_pulsepower . shape [ 0 ]
M0 = np . array ( [ np . zeros ( Nx ) , np . zeros ( Nx ) , np . ones ( Nx ) ] )
2023-08-23 14:47:41 +00:00
dt = self . pulse . dwell_time * 1e3 # We need our dwell time in ms
2023-08-22 14:06:38 +00:00
2023-08-23 14:47:41 +00:00
w = np . ones ( ( Nu , 1 ) ) * self . gradient
2023-08-23 12:51:57 +00:00
# Bloch simulation in magnetization domain
2023-08-23 14:47:41 +00:00
gadt = self . sample . gamma * dt / 2
B1 = np . tile (
( gadt * ( real_pulsepower - 1 j * imag_pulsepower ) * B1 ) . reshape ( - 1 , 1 ) , Nx
)
2024-01-10 11:26:39 +00:00
2023-08-23 12:51:57 +00:00
K = gadt * xdis * w * self . gradient
2023-08-23 14:47:41 +00:00
phi = - np . sqrt ( np . abs ( B1 ) * * 2 + K * * 2 )
2023-08-22 14:06:38 +00:00
cs = np . cos ( phi )
si = np . sin ( phi )
n1 = np . real ( B1 ) / np . abs ( phi )
n2 = np . imag ( B1 ) / np . abs ( phi )
n3 = K / np . abs ( phi )
n1 [ np . isnan ( n1 ) ] = 1
n2 [ np . isnan ( n2 ) ] = 0
n3 [ np . isnan ( n3 ) ] = 0
Bd1 = n1 * n1 * ( 1 - cs ) + cs
Bd2 = n1 * n2 * ( 1 - cs ) - n3 * si
Bd3 = n1 * n3 * ( 1 - cs ) + n2 * si
Bd4 = n2 * n1 * ( 1 - cs ) + n3 * si
Bd5 = n2 * n2 * ( 1 - cs ) + cs
Bd6 = n2 * n3 * ( 1 - cs ) - n1 * si
Bd7 = n3 * n1 * ( 1 - cs ) - n2 * si
Bd8 = n3 * n2 * ( 1 - cs ) + n1 * si
Bd9 = n3 * n3 * ( 1 - cs ) + cs
2023-08-23 14:47:41 +00:00
M = np . zeros ( ( 3 , Nx , Nu + 1 ) )
2023-08-22 14:06:38 +00:00
M [ : , : , 0 ] = M0
Mt = M0
2023-08-23 14:47:41 +00:00
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
2023-08-23 12:51:57 +00:00
Mrot = np . zeros ( ( 3 , Nx ) )
2023-08-23 14:47:41 +00:00
Mrot [ 0 , : ] = (
2024-01-10 11:26:39 +00:00
Bd1 . conj ( ) . T [ : , n ] * Mt [ 0 , : ] + Bd2 . conj ( ) . T [ : , n ] * Mt [ 1 , : ] + Bd3 . conj ( ) . T [ : , n ] * Mt [ 2 , : ]
2023-08-23 14:47:41 +00:00
)
Mrot [ 1 , : ] = (
2024-01-10 11:26:39 +00:00
Bd4 . conj ( ) . T [ : , n ] * Mt [ 0 , : ] + Bd5 . conj ( ) . T [ : , n ] * Mt [ 1 , : ] + Bd6 . conj ( ) . T [ : , n ] * Mt [ 2 , : ]
2023-08-23 14:47:41 +00:00
)
Mrot [ 2 , : ] = (
2024-01-10 11:26:39 +00:00
Bd7 . conj ( ) . T [ : , n ] * Mt [ 0 , : ] + Bd8 . conj ( ) . T [ : , n ] * Mt [ 1 , : ] + Bd9 . conj ( ) . T [ : , n ] * Mt [ 2 , : ]
2023-08-23 14:47:41 +00:00
)
2023-08-23 12:51:57 +00:00
Mt = np . dot ( D , Mrot ) + np . tile ( b , ( Nx , 1 ) ) . T
2023-08-23 14:47:41 +00:00
Mrot [ 0 , : ] = (
2024-01-10 11:26:39 +00:00
Bd1 . conj ( ) . T [ : , n ] * Mt [ 0 , : ] + Bd2 . conj ( ) . T [ : , n ] * Mt [ 1 , : ] + Bd3 . conj ( ) . T [ : , n ] * Mt [ 2 , : ]
2023-08-23 14:47:41 +00:00
)
Mrot [ 1 , : ] = (
2024-01-10 11:26:39 +00:00
Bd4 . conj ( ) . T [ : , n ] * Mt [ 0 , : ] + Bd5 . conj ( ) . T [ : , n ] * Mt [ 1 , : ] + Bd6 . conj ( ) . T [ : , n ] * Mt [ 2 , : ]
2023-08-23 14:47:41 +00:00
)
Mrot [ 2 , : ] = (
2024-01-10 11:26:39 +00:00
Bd7 . conj ( ) . T [ : , n ] * Mt [ 0 , : ] + Bd8 . conj ( ) . T [ : , n ] * Mt [ 1 , : ] + Bd9 . conj ( ) . T [ : , n ] * Mt [ 2 , : ]
2023-08-23 14:47:41 +00:00
)
2023-08-22 14:06:38 +00:00
Mt = Mrot
2023-08-23 14:47:41 +00:00
M [ : , : , n + 1 ] = Mrot
2023-08-22 14:06:38 +00:00
return M
2023-08-23 12:51:57 +00:00
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.
2023-08-23 14:47:41 +00:00
2023-08-23 12:51:57 +00:00
Returns
- - - - - - -
B1 : float
2023-08-23 13:48:17 +00:00
The B1 field of the solenoid coil in T . """
2023-08-23 12:51:57 +00:00
2023-08-23 14:47:41 +00:00
B1 = (
np . sqrt ( 2 * self . power_amplifier_power / 50 )
* np . pi
2024-03-14 13:56:06 +00:00
* np . sqrt ( self . q_factor_transmit )
2023-08-23 14:47:41 +00:00
* 4e-7
* self . number_turns
/ self . length_coil
)
2023-08-23 12:51:57 +00:00
return B1
2023-08-23 14:47:41 +00:00
2023-08-23 12:51:57 +00:00
def calc_xdis ( self ) - > np . array :
2023-08-23 14:47:41 +00:00
""" Calculates the x distribution of the isochromats.
2023-08-23 12:51:57 +00:00
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
# 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
2023-08-23 14:47:41 +00:00
xdis = np . linspace ( - 1 , 1 , self . number_isochromats )
xdis = (
( foffr . T * 1e-6 ) / ( self . sample . gamma / 2 / np . pi ) / ( self . gradient * 1e-3 )
)
2023-08-23 13:48:17 +00:00
2023-08-23 12:51:57 +00:00
return xdis
2023-08-23 14:47:41 +00:00
2023-08-23 14:44:39 +00:00
def calculate_reference_voltage ( self ) - > float :
2023-08-23 14:18:22 +00:00
""" This calculates the reference voltage of the measurement setup for the sample at a certain temperature.
Returns
- - - - - - -
reference_voltage : float
The reference voltage of the measurement setup for the sample at a certain temperature in Volts .
"""
2023-08-23 14:47:41 +00:00
u = 4 * np . pi * 1e-7 # permeability of free space
magnetization = (
2023-12-14 12:56:00 +00:00
( ( self . sample . gamma
2023-08-23 14:47:41 +00:00
* 2
2023-12-14 12:56:00 +00:00
* self . sample . atoms )
/ ( 2 * self . sample . nuclear_spin + 1 ) )
* ( h * * 2
2023-12-14 15:42:19 +00:00
* self . sample . resonant_frequency )
2023-12-14 12:56:00 +00:00
/ ( Boltzmann
* self . temperature )
2023-08-23 14:47:41 +00:00
* self . sample . spin_factor
)
2023-08-23 14:18:22 +00:00
coil_crossection = np . pi * ( self . diameter_coil / 2 ) * * 2
2023-12-14 12:56:00 +00:00
2023-08-23 14:47:41 +00:00
reference_voltage = (
self . number_turns
* coil_crossection
* u
2023-12-14 15:42:19 +00:00
* ( self . sample . resonant_frequency )
2023-12-14 12:56:00 +00:00
* magnetization
2023-08-23 14:47:41 +00:00
)
reference_voltage = (
reference_voltage * self . sample . powder_factor * self . sample . filling_factor
)
2023-12-14 11:39:44 +00:00
2023-12-14 15:42:19 +00:00
# This is assumes that our noise is dominated by everything after the resonator - this is not true for low Q probe coils
2023-12-14 11:39:44 +00:00
reference_voltage = reference_voltage * np . sqrt ( self . q_factor_receive )
2023-08-23 14:18:22 +00:00
return reference_voltage
2023-08-23 14:47:41 +00:00
def calculate_noise ( self , timedomain_signal : np . array ) - > np . array :
2023-08-23 14:44:39 +00:00
""" Calculates the noise array that is added to the signal.
2023-08-23 14:47:41 +00:00
2023-08-23 14:44:39 +00:00
Parameters
- - - - - - - - - -
timedomain_signal : np . array
The time domain signal that is used for the simulation .
2023-08-23 14:47:41 +00:00
2023-08-23 14:44:39 +00:00
Returns
- - - - - - -
noise_data : np . array
The noise array that is added to the signal . """
n_timedomain_points = timedomain_signal . shape [ 0 ]
2023-12-14 12:56:00 +00:00
noise_data = self . noise * 1e-6 * np . random . randn (
2023-08-23 14:47:41 +00:00
self . averages , n_timedomain_points
2023-12-14 12:56:00 +00:00
) + 1 j * self . noise * 1e-6 * np . random . randn ( self . averages , n_timedomain_points )
2023-08-23 14:44:39 +00:00
noise_data = np . sum ( noise_data , axis = 0 ) # Sum along the first axis
return noise_data
2023-08-23 12:51:57 +00:00
@property
def sample ( self ) - > Sample :
""" Sample that is used for the simulation. """
return self . _sample
2023-08-23 14:47:41 +00:00
2023-08-23 12:51:57 +00:00
@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
2023-08-23 14:47:41 +00:00
2023-08-23 12:51:57 +00:00
@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
2023-08-23 14:47:41 +00:00
2023-08-23 12:51:57 +00:00
@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
2023-08-23 14:47:41 +00:00
2023-08-23 12:51:57 +00:00
@gradient.setter
def gradient ( self , gradient ) :
self . _gradient = gradient
@property
def noise ( self ) - > float :
2023-08-23 14:47:41 +00:00
""" RMS Noise of the measurement setup in Volts """
2023-08-23 12:51:57 +00:00
return self . _noise
2023-08-23 14:47:41 +00:00
2023-08-23 12:51:57 +00:00
@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
2023-08-23 14:47:41 +00:00
2023-08-23 12:51:57 +00:00
@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
2023-08-23 14:47:41 +00:00
2023-08-23 12:51:57 +00:00
@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
2023-08-23 14:47:41 +00:00
2023-08-23 12:51:57 +00:00
@number_turns.setter
def number_turns ( self , number_turns ) :
self . _number_turns = number_turns
2023-12-14 11:39:44 +00:00
@property
2024-03-14 13:56:06 +00:00
def q_factor_transmit ( self ) - > float :
2024-03-02 19:58:07 +00:00
""" Q-factor of the transmit path of the probe coil. """
return self . _q_factor_transmit
2023-12-14 11:39:44 +00:00
2024-03-02 19:58:07 +00:00
@q_factor_transmit.setter
def q_factor_transmit ( self , q_factor_transmit ) :
self . _q_factor_transmit = q_factor_transmit
2023-12-14 11:39:44 +00:00
@property
def q_factor_receive ( self ) - > float :
""" Q-factor of the receive path of the probe coil. """
return self . _q_factor_receive
@q_factor_receive.setter
def q_factor_receive ( self , q_factor_receive ) :
self . _q_factor_receive = q_factor_receive
2023-08-23 12:51:57 +00:00
@property
def power_amplifier_power ( self ) - > float :
""" Power of the power amplifier in Watts. """
return self . _power_amplifier_power
2023-08-23 14:47:41 +00:00
2023-08-23 12:51:57 +00:00
@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
2023-08-23 14:47:41 +00:00
2023-08-23 12:51:57 +00:00
@pulse.setter
def pulse ( self , pulse ) :
self . _pulse = pulse
2023-08-23 13:48:17 +00:00
@property
def averages ( self ) - > int :
""" Number of averages that are used for the simulation. """
return self . _averages
2023-08-23 14:47:41 +00:00
2023-08-23 13:48:17 +00:00
@averages.setter
def averages ( self , averages ) :
self . _averages = averages
@property
def gain ( self ) - > float :
""" Gain of the amplifier. """
return self . _gain
2023-08-23 14:47:41 +00:00
2023-08-23 13:48:17 +00:00
@gain.setter
def gain ( self , gain ) :
self . _gain = gain
2023-08-23 14:18:22 +00:00
@property
def temperature ( self ) - > float :
""" Temperature of the sample. """
return self . _temperature
2023-08-23 14:47:41 +00:00
2023-08-23 14:18:22 +00:00
@temperature.setter
def temperature ( self , temperature ) :
self . _temperature = temperature