nqrduck/quackseq-simulator
1
0
Fork 0
mirror of https://github.com/nqrduck/quackseq-simulator.git synced 2024-12-30 08:38:05 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

Initial commit.

Browse source
This commit is contained in:
jupfi 2024-05-29 08:36:34 +02:00
commit fe458ed59e
10 changed files with 887 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

24
.gitignore vendored Normal file
View file

@ -0,0 +1,24 @@
# Byte-compiled / optimized / DLL files
__pycache__/
*.pyc
*$py.class
# Distribution / packaging
dist/
build/
*.egg-info/
# IDE-specific files
.idea/
.vscode/
# Logs
*.log
# Virtual environments
venv/
# Other
*.DS_Store
*.pos
*.quack

5
CHANGELOG.md Normal file
View file

@ -0,0 +1,5 @@
# Changelog
## Version 0.0.1 (15-04-2024)
- Initial release

20
LICENSE Normal file
View file

@ -0,0 +1,20 @@
MIT License
Copyright (c) 2023-2024 jupfi
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

1
README.md Normal file
View file

@ -0,0 +1 @@
# quackseq-simulator

61
pyproject.toml Normal file
View file

@ -0,0 +1,61 @@
[build-system]
requires = ["hatchling"]
build-backend = "hatchling.build"
[tool.hatch.metadata]
allow-direct-references = true
[project]
name = "quackseq-simulator"
version = "0.0.1"
authors = [
{ name="jupfi", email="support@nqrduck.cool" },
]
description = "Simple Python script to perform magnetic resonance spectroscopy experiments."
readme = "README.md"
license = { file="LICENSE" }
requires-python = ">=3.10"
classifiers = [
"Programming Language :: Python :: 3",
"License :: OSI Approved :: MIT License",
"Operating System :: OS Independent",
]
dependencies = [
"numpy",
"scipy",
"nqr-blochsimulator",
"quackseq",
]
[project.optional-dependencies]
dev = [
"black",
"pydocstyle",
"pyupgrade",
"ruff",
]
[tool.ruff]
[tool.ruff.lint]
extend-select = [
"UP", # pyupgrade
"D", # pydocstyle
]
[tool.ruff.lint.per-file-ignores]
"__init__.py" = ["F401"]
[tool.ruff.lint.pydocstyle]
convention = "google"
[project.urls]
"Homepage" = "https://nqrduck.cool"
"Bug Tracker" = "https://github.com/nqrduck/quackseq/issues"
"Source Code" = "https://github.com/nqrduck/quackseq"
[tool.hatch.build.targets.wheel]
packages = ["src/quackseq_simulator"]

0
src/quackseq_simulator/__init__.py Normal file
View file

17
src/quackseq_simulator/simulator.py Normal file
View file

@ -0,0 +1,17 @@
from quackseq.spectrometer.spectrometer import Spectrometer
from .simulator_model import SimulatorModel
from .simulator_controller import SimulatorController
class Simulator(Spectrometer):
def __init__(self):
self.model = SimulatorModel()
self.controller = SimulatorController(self.model)
def run_sequence(self, sequence):
result = self.controller.run_sequence(sequence)
return result
def set_averages(self, value: int):
self.model.average = value

354
src/quackseq_simulator/simulator_controller.py Normal file
View file

@ -0,0 +1,354 @@
"""The controller module for the simulator spectrometer."""
import logging
from datetime import datetime
import numpy as np
from quackseq.spectrometer.spectrometer_controller import SpectrometerController
from quackseq.measurement import Measurement
from quackseq.pulseparameters import TXPulse, RXReadout
from quackseq.pulsesequence import QuackSequence
from nqr_blochsimulator.classes.pulse import PulseArray
from nqr_blochsimulator.classes.sample import Sample
from nqr_blochsimulator.classes.simulation import Simulation
logger = logging.getLogger(__name__)
class SimulatorController(SpectrometerController):
"""The controller class for the nqrduck simulator module."""
def __init__(self, model):
"""Initializes the SimulatorController."""
super().__init__()
self.model = model
def run_sequence(self, sequence: QuackSequence) -> None:
"""This method is called when the start_measurement signal is received from the core.
It will becalled if the simulator is the active spectrometer.
This will start the simulation based on the settings and the pulse sequence.
"""
logger.debug("Starting simulation")
sample = self.get_sample_from_settings()
logger.debug("Sample: %s", sample.name)
dwell_time = self.calculate_dwelltime(sequence)
logger.debug("Dwell time: %s", dwell_time)
try:
pulse_array = self.translate_pulse_sequence(sequence, dwell_time)
except AttributeError:
logger.warning("Could not translate pulse sequence")
return
simulation = self.get_simulation(sample, pulse_array)
result = simulation.simulate()
tdx = (
np.linspace(0, float(self.calculate_simulation_length(sequence)), len(result)) * 1e6
)
rx_begin, rx_stop = self.translate_rx_event(sequence)
# If we have a RX event, we need to cut the result to the RX event
if rx_begin and rx_stop:
evidx = np.where((tdx > rx_begin) & (tdx < rx_stop))[0]
tdx = tdx[evidx]
result = result[evidx]
# Measurement name date + module + target frequency + averages + sequence name
name = f"{datetime.now().strftime('%Y-%m-%d %H:%M:%S')} - Simulator - {self.model.target_frequency / 1e6} MHz - {self.model.averages} averages - {sequence.name}"
logger.debug(f"Measurement name: {name}")
measurement_data = Measurement(
name,
tdx,
result / simulation.averages,
sample.resonant_frequency,
# frequency_shift=self.module.model.if_frequency,
)
return measurement_data
def get_sample_from_settings(self) -> Sample:
"""This method creates a sample object based on the settings in the model.
Returns:
Sample: The sample object created from the settings.
"""
model = self.model
atom_density = None
sample_volume = None
sample_length = None
sample_diameter = None
for samplesetting in model.settings[self.model.SAMPLE]:
logger.debug("Sample setting: %s", samplesetting.name)
if samplesetting.name == model.NAME:
name = samplesetting.value
elif samplesetting.name == model.DENSITY:
density = float(samplesetting.value)
elif samplesetting.name == model.MOLAR_MASS:
molar_mass = float(samplesetting.value)
elif samplesetting.name == model.RESONANT_FREQUENCY:
resonant_frequency = float(samplesetting.value)
elif samplesetting.name == model.GAMMA:
gamma = float(samplesetting.value)
elif samplesetting.name == model.NUCLEAR_SPIN:
nuclear_spin = float(samplesetting.value)
elif samplesetting.name == model.SPIN_FACTOR:
spin_factor = float(samplesetting.value)
elif samplesetting.name == model.POWDER_FACTOR:
powder_factor = float(samplesetting.value)
elif samplesetting.name == model.FILLING_FACTOR:
filling_factor = float(samplesetting.value)
elif samplesetting.name == model.T1:
T1 = float(samplesetting.value)
elif samplesetting.name == model.T2:
T2 = float(samplesetting.value)
elif samplesetting.name == model.T2_STAR:
T2_star = float(samplesetting.value)
elif samplesetting.name == model.ATOM_DENSITY:
atom_density = float(samplesetting.value)
elif samplesetting.name == model.SAMPLE_VOLUME:
sample_volume = float(samplesetting.value)
elif samplesetting.name == model.SAMPLE_LENGTH:
sample_length = float(samplesetting.value)
elif samplesetting.name == model.SAMPLE_DIAMETER:
sample_diameter = float(samplesetting.value)
else:
logger.warning("Unknown sample setting: %s", samplesetting.name)
self.module.nqrduck_signal.emit(
"notification",
["Error", "Unknown sample setting: " + samplesetting.name],
)
return None
sample = Sample(
name=name,
density=density,
molar_mass=molar_mass,
resonant_frequency=resonant_frequency,
gamma=gamma,
nuclear_spin=nuclear_spin,
spin_factor=spin_factor,
powder_factor=powder_factor,
filling_factor=filling_factor,
T1=T1,
T2=T2,
T2_star=T2_star,
atom_density=atom_density,
sample_volume=sample_volume,
sample_length=sample_length,
sample_diameter=sample_diameter,
)
return sample
def translate_pulse_sequence(self, sequence : QuackSequence, dwell_time: float) -> PulseArray:
"""This method translates the pulse sequence from the core to a PulseArray object needed for the simulation.
Args:
sequence (QuackSequence): The pulse sequence from the core.
dwell_time (float): The dwell time in seconds.
Returns:
PulseArray: The pulse sequence translated to a PulseArray object.
"""
events = sequence.events
amplitude_array = list()
for event in events:
logger.debug("Event %s has parameters: %s", event.name, event.parameters)
for parameter in event.parameters.values():
logger.debug(
"Parameter %s has options: %s", parameter.name, parameter.options
)
if (
parameter.name == sequence.TX_PULSE
and parameter.get_option_by_name(TXPulse.RELATIVE_AMPLITUDE).value
> 0
):
logger.debug(f"Adding pulse: {event.duration} s")
# If we have a pulse, we need to add it to the pulse array
pulse_shape = parameter.get_option_by_name(
TXPulse.TX_PULSE_SHAPE
).value
pulse_amplitude = abs(
pulse_shape.get_pulse_amplitude(
event.duration, resolution=dwell_time
)
)
amplitude_array.append(pulse_amplitude)
elif (
parameter.name == sequence.TX_PULSE
and parameter.get_option_by_name(TXPulse.RELATIVE_AMPLITUDE).value
== 0
):
# If we have a wait, we need to add it to the pulse array
amplitude_array.append(np.zeros(int(event.duration / dwell_time)))
amplitude_array = np.concatenate(amplitude_array)
# This has not yet been implemented right now the phase is always 0
phase_array = np.zeros(len(amplitude_array))
pulse_array = PulseArray(
pulseamplitude=amplitude_array,
pulsephase=phase_array,
dwell_time=float(dwell_time),
)
return pulse_array
def get_simulation(self, sample: Sample, pulse_array: PulseArray) -> Simulation:
"""This method creates a simulation object based on the settings and the pulse sequence.
Args:
sample (Sample): The sample object created from the settings.
pulse_array (PulseArray): The pulse sequence translated to a PulseArray object.
Returns:
Simulation: The simulation object created from the settings and the pulse sequence.
"""
model = self.model
# noise = float(model.get_setting_by_name(model.NOISE).value)
simulation = Simulation(
sample=sample,
pulse=pulse_array,
number_isochromats=int(
model.get_setting_by_name(model.NUMBER_ISOCHROMATS).value
),
initial_magnetization=float(
model.get_setting_by_name(model.INITIAL_MAGNETIZATION).value
),
gradient=float(model.get_setting_by_name(model.GRADIENT).value),
noise=float(model.get_setting_by_name(model.NOISE).value),
length_coil=float(model.get_setting_by_name(model.LENGTH_COIL).value),
diameter_coil=float(model.get_setting_by_name(model.DIAMETER_COIL).value),
number_turns=float(model.get_setting_by_name(model.NUMBER_TURNS).value),
q_factor_transmit=float(
model.get_setting_by_name(model.Q_FACTOR_TRANSMIT).value
),
q_factor_receive=float(
model.get_setting_by_name(model.Q_FACTOR_RECEIVE).value
),
power_amplifier_power=float(
model.get_setting_by_name(model.POWER_AMPLIFIER_POWER).value
),
gain=float(model.get_setting_by_name(model.GAIN).value),
temperature=float(model.get_setting_by_name(model.TEMPERATURE).value),
averages=int(model.averages),
loss_TX=float(model.get_setting_by_name(model.LOSS_TX).value),
loss_RX=float(model.get_setting_by_name(model.LOSS_RX).value),
conversion_factor=float(
model.get_setting_by_name(model.CONVERSION_FACTOR).value
),
)
return simulation
def calculate_dwelltime(self, sequence : QuackSequence) -> float:
"""This method calculates the dwell time based on the settings and the pulse sequence.
Returns:
float: The dwell time in seconds.
"""
n_points = int(
self.model.get_setting_by_name(self.model.NUMBER_POINTS).value
)
simulation_length = self.calculate_simulation_length(sequence)
dwell_time = simulation_length / n_points
return dwell_time
def calculate_simulation_length(self, sequence : QuackSequence) -> float:
"""This method calculates the simulation length based on the settings and the pulse sequence.
Returns:
float: The simulation length in seconds.
"""
events = sequence.events
simulation_length = 0
for event in events:
simulation_length += event.duration
return simulation_length
def translate_rx_event(self, sequence : QuackSequence) -> tuple:
"""This method translates the RX event of the pulse sequence to the limr object.
Returns:
tuple: A tuple containing the start and stop time of the RX event in µs
"""
# This is a correction factor for the RX event. The offset of the first pulse is 2.2µs longer than from the specified samples.
events = sequence.events
previous_events_duration = 0
# offset = 0
rx_duration = 0
for event in events:
logger.debug("Event %s has parameters: %s", event.name, event.parameters)
for parameter in event.parameters.values():
logger.debug(
"Parameter %s has options: %s", parameter.name, parameter.options
)
if (
parameter.name == sequence.RX_READOUT
and parameter.get_option_by_name(RXReadout.RX).value
):
# Get the length of all previous events
previous_events = events[: events.index(event)]
previous_events_duration = sum(
[event.duration for event in previous_events]
)
rx_duration = event.duration
rx_begin = float(previous_events_duration)
if rx_duration:
rx_stop = rx_begin + float(rx_duration)
return rx_begin * 1e6, rx_stop * 1e6
else:
return None, None
def set_frequency(self, value: str) -> None:
"""This method is called when the set_frequency signal is received from the core.
For the simulator this just prints a warning that the simulator is selected.
Args:
value (str) : The new frequency in MHz.
"""
logger.debug("Setting frequency to: %s", value)
try:
self.module.model.target_frequency = float(value)
logger.debug("Successfully set frequency to: %s", value)
except ValueError:
logger.warning("Could not set frequency to: %s", value)
self.module.nqrduck_signal.emit(
"notification", ["Error", "Could not set frequency to: " + value]
)
self.module.nqrduck_signal.emit("failure_set_frequency", value)
def set_averages(self, value: str) -> None:
"""This method is called when the set_averages signal is received from the core.
It sets the averages in the model used for the simulation.
Args:
value (str): The value to set the averages to.
"""
logger.debug("Setting averages to: %s", value)
try:
self.module.model.averages = int(value)
logger.debug("Successfully set averages to: %s", value)
except ValueError:
logger.warning("Could not set averages to: %s", value)
self.module.nqrduck_signal.emit(
"notification", ["Error", "Could not set averages to: " + value]
)
self.module.nqrduck_signal.emit("failure_set_averages", value)

327
src/quackseq_simulator/simulator_model.py Normal file
View file

@ -0,0 +1,327 @@
"""The model module for the simulator spectrometer."""
import logging
from quackseq.spectrometer.spectrometer_model import SpectrometerModel
from quackseq.spectrometer.spectrometer_settings import IntSetting, FloatSetting, StringSetting
from quackseq.pulseparameters import TXPulse, RXReadout
logger = logging.getLogger(__name__)
class SimulatorModel(SpectrometerModel):
"""Model class for the simulator spectrometer."""
# Simulation settings
NUMBER_POINTS = "N. simulation points"
NUMBER_ISOCHROMATS = "N. of isochromats"
INITIAL_MAGNETIZATION = "Initial magnetization"
GRADIENT = "Gradient (mT/m))"
NOISE = "Noise (uV)"
# Hardware settings
LENGTH_COIL = "Length coil (m)"
DIAMETER_COIL = "Diameter coil (m)"
NUMBER_TURNS = "Number turns"
Q_FACTOR_TRANSMIT = "Q factor Transmit"
Q_FACTOR_RECEIVE = "Q factor Receive"
POWER_AMPLIFIER_POWER = "PA power (W)"
GAIN = "Gain"
TEMPERATURE = "Temperature (K)"
AVERAGES = "Averages"
LOSS_TX = "Loss TX (dB)"
LOSS_RX = "Loss RX (dB)"
CONVERSION_FACTOR = "Conversion factor"
# Sample settings, this will be done in a separate module later on
NAME = "Name"
DENSITY = "Density (g/cm^3)"
MOLAR_MASS = "Molar mass (g/mol)"
RESONANT_FREQUENCY = "Resonant freq. (Hz)"
GAMMA = "Gamma (Hz/T)"
NUCLEAR_SPIN = "Nuclear spin"
SPIN_FACTOR = "Spin factor"
POWDER_FACTOR = "Powder factor"
FILLING_FACTOR = "Filling factor"
T1 = "T1 (s)"
T2 = "T2 (s)"
T2_STAR = "T2* (s)"
ATOM_DENSITY = "Atom density (1/cm^3)"
SAMPLE_VOLUME = "Sample volume (m^3)"
SAMPLE_LENGTH = "Sample length (m)"
SAMPLE_DIAMETER = "Sample diameter (m)"
# Categories of the settings
SIMULATION = "Simulation"
HARDWARE = "Hardware"
EXPERIMENTAL_Setup = "Experimental Setup"
SAMPLE = "Sample"
def __init__(self):
"""Initializes the SimulatorModel."""
super().__init__()
# Simulation settings
number_of_points_setting = IntSetting(
self.NUMBER_POINTS,
8192,
"Number of points used for the simulation. This influences the dwell time in combination with the total event simulation given by the pulse sequence.",
min_value=0,
)
self.add_setting(
number_of_points_setting,
self.SIMULATION,
)
number_of_isochromats_setting = IntSetting(
self.NUMBER_ISOCHROMATS,
1000,
"Number of isochromats used for the simulation. This influences the computation time.",
min_value=0,
max_value=10000,
)
self.add_setting(number_of_isochromats_setting, self.SIMULATION)
initial_magnetization_setting = FloatSetting(
self.INITIAL_MAGNETIZATION,
1,
"Initial magnetization",
min_value=0,
)
self.add_setting(initial_magnetization_setting, self.SIMULATION)
# This doesn't really do anything yet
gradient_setting = FloatSetting(
self.GRADIENT,
1,
"Gradient",
)
self.add_setting(gradient_setting, self.SIMULATION)
noise_setting = FloatSetting(
self.NOISE,
2,
"Adds a specified level of random noise to the simulation to mimic real-world signal variations.",
min_value=0,
max_value=100,
)
self.add_setting(noise_setting, self.SIMULATION)
# Hardware settings
coil_length_setting = FloatSetting(
self.LENGTH_COIL,
30e-3,
"The length of the sample coil within the hardware setup.",
min_value=1e-3,
)
self.add_setting(coil_length_setting, self.HARDWARE)
coil_diameter_setting = FloatSetting(
self.DIAMETER_COIL,
8e-3,
"The diameter of the sample coil.",
min_value=1e-3,
)
self.add_setting(coil_diameter_setting, self.HARDWARE)
number_turns_setting = FloatSetting(
self.NUMBER_TURNS,
8,
"The total number of turns of the sample coil.",
min_value=1,
)
self.add_setting(number_turns_setting, self.HARDWARE)
q_factor_transmit_setting = FloatSetting(
self.Q_FACTOR_TRANSMIT,
80,
"The quality factor of the transmit path, which has an effect on the field strength for excitation.",
min_value=1,
)
self.add_setting(q_factor_transmit_setting, self.HARDWARE)
q_factor_receive_setting = FloatSetting(
self.Q_FACTOR_RECEIVE,
80,
"The quality factor of the receive path, which has an effect on the final SNR.",
min_value=1,
)
self.add_setting(q_factor_receive_setting, self.HARDWARE)
power_amplifier_power_setting = FloatSetting(
self.POWER_AMPLIFIER_POWER,
110,
"The power output capability of the power amplifier, determines the strength of pulses that can be generated.",
min_value=0.1,
)
self.add_setting(power_amplifier_power_setting, self.HARDWARE)
gain_setting = FloatSetting(
self.GAIN,
6000,
"The amplification factor of the receiver chain, impacting the final measured signal amplitude.",
min_value=0.1,
)
self.add_setting(gain_setting, self.HARDWARE)
temperature_setting = FloatSetting(
self.TEMPERATURE,
300,
"The absolute temperature during the experiment. This influences the SNR of the measurement.",
min_value=0.1,
max_value=400,
)
self.add_setting(temperature_setting, self.EXPERIMENTAL_Setup)
loss_tx_setting = FloatSetting(
self.LOSS_TX,
25,
"The signal loss occurring in the transmission path, affecting the effective RF pulse power.",
min_value=0.1,
max_value=60,
)
self.add_setting(loss_tx_setting, self.EXPERIMENTAL_Setup)
loss_rx_setting = FloatSetting(
self.LOSS_RX,
25,
"The signal loss in the reception path, which can reduce the signal that is ultimately detected.",
min_value=0.1,
max_value=60,
)
self.add_setting(loss_rx_setting, self.EXPERIMENTAL_Setup)
conversion_factor_setting = FloatSetting(
self.CONVERSION_FACTOR,
2884,
"Conversion factor (spectrometer units / V)",
)
self.add_setting(
conversion_factor_setting,
self.EXPERIMENTAL_Setup,
) # Conversion factor for the LimeSDR based spectrometer
# Sample settings
sample_name_setting = StringSetting(
self.NAME,
"BiPh3",
"The name of the sample.",
)
self.add_setting(sample_name_setting, self.SAMPLE)
density_setting = FloatSetting(
self.DENSITY,
1.585e6,
"The density of the sample. This is used to calculate the number of spins in the sample volume.",
min_value=0.1,
)
self.add_setting(density_setting, self.SAMPLE)
molar_mass_setting = FloatSetting(
self.MOLAR_MASS,
440.3,
"The molar mass of the sample. This is used to calculate the number of spins in the sample volume.",
min_value=0.1,
)
self.add_setting(molar_mass_setting, self.SAMPLE)
resonant_frequency_setting = FloatSetting(
self.RESONANT_FREQUENCY,
83.56e6,
"The resonant frequency of the observed transition.",
min_value=1e5,
)
self.add_setting(resonant_frequency_setting, self.SAMPLE)
gamma_setting = FloatSetting(
self.GAMMA,
4.342e7,
"The gyromagnetic ratio of the sample’s nuclei.",
min_value=0,
)
self.add_setting(gamma_setting, self.SAMPLE)
# This could be updated to a selection setting
nuclear_spin_setting = FloatSetting(
self.NUCLEAR_SPIN,
9 / 2,
"The nuclear spin of the sample’s nuclei.",
min_value=0,
)
self.add_setting(nuclear_spin_setting, self.SAMPLE)
spin_factor_setting = FloatSetting(
self.SPIN_FACTOR,
2,
"The spin factor represents the scaling coefficient for observable nuclear spin transitions along the x-axis, derived from the Pauli I x 0 -matrix elements.",
min_value=0,
)
self.add_setting(spin_factor_setting, self.SAMPLE)
powder_factor_setting = FloatSetting(
self.POWDER_FACTOR,
0.75,
"A factor representing the crystallinity of the solid sample. A value of 0.75 corresponds to a powder sample.",
min_value=0,
max_value=1,
)
self.add_setting(powder_factor_setting, self.SAMPLE)
filling_factor_setting = FloatSetting(
self.FILLING_FACTOR,
0.7,
"The ratio of the sample volume that occupies the coil’s sensitive volume.",
min_value=0,
max_value=1,
)
self.add_setting(filling_factor_setting, self.SAMPLE)
t1_setting = FloatSetting(
self.T1,
83.5e-5,
"The longitudinal or spin-lattice relaxation time of the sample, influencing signal recovery between pulses.",
min_value=1e-6,
)
self.add_setting(t1_setting, self.SAMPLE)
t2_setting = FloatSetting(
self.T2,
396e-6,
"The transverse or spin-spin relaxation time, determining the rate at which spins dephase and the signal decays in the xy plane",
min_value=1e-6,
)
self.add_setting(t2_setting, self.SAMPLE)
t2_star_setting = FloatSetting(
self.T2_STAR,
50e-6,
"The effective transverse relaxation time, incorporating effects of EFG inhomogeneities and other dephasing factors.",
min_value=1e-6,
)
self.add_setting(t2_star_setting, self.SAMPLE)
self.averages = 1
self.target_frequency = 100e6
@property
def averages(self):
"""The number of averages used for the simulation.
More averages improve the signal-to-noise ratio of the simulated signal.
"""
return self._averages
@averages.setter
def averages(self, value):
self._averages = value
@property
def target_frequency(self):
"""The target frequency for the simulation.
Doesn't do anything at the moment.
"""
return self._target_frequency
@target_frequency.setter
def target_frequency(self, value):
self._target_frequency = value

78
tests/simulator.py Normal file
View file

@ -0,0 +1,78 @@
import unittest
import logging
import matplotlib.pyplot as plt
from quackseq.pulsesequence import QuackSequence
from quackseq.event import Event
from quackseq.functions import RectFunction
from quackseq_simulator.simulator import Simulator
# logging.basicConfig(level=logging.DEBUG)
class TestQuackSequence(unittest.TestCase):
def test_event_creation(self):
seq = QuackSequence("test - event creation")
seq.add_pulse_event("tx", "10u", 1, 0, RectFunction())
seq.add_blank_event("blank", "3u")
seq.add_readout_event("rx", "100u")
seq.add_blank_event("TR", "1m")
json = seq.to_json()
print(json)
sim = Simulator()
sim.set_averages(100)
result = sim.run_sequence(seq)
self.assertIsNotNone(result)
self.assertTrue(hasattr(result, "tdx"))
self.assertTrue(hasattr(result, "tdy"))
self.assertGreater(len(result.tdx), 0)
self.assertGreater(len(result.tdy), 0)
# Plotting the result can be useful for visual inspection during development
plt.plot(result.tdx, abs(result.tdy))
plt.show()
def test_simulation_run_sequence(self):
seq = QuackSequence("test - simulation run sequence")
tx = Event("tx", "10u", seq)
seq.add_event(tx)
seq.set_tx_amplitude(tx, 1)
seq.set_tx_phase(tx, 0)
json = seq.to_json()
print(json)
rect = RectFunction()
seq.set_tx_shape(tx, rect)
blank = Event("blank", "3u", seq)
seq.add_event(blank)
rx = Event("rx", "100u", seq)
seq.set_rx(rx, True)
seq.add_event(rx)
TR = Event("TR", "1m", seq)
seq.add_event(TR)
sim = Simulator()
sim.set_averages(100)
result = sim.run_sequence(seq)
self.assertIsNotNone(result)
self.assertTrue(hasattr(result, "tdx"))
self.assertTrue(hasattr(result, "tdy"))
self.assertGreater(len(result.tdx), 0)
self.assertGreater(len(result.tdy), 0)
# Plotting the result can be useful for visual inspection during development
plt.plot(result.tdx, abs(result.tdy))
plt.show()
if __name__ == "__main__":
unittest.main()