commit fe458ed59ec6d65c6acb6a59a259b641cbd93753 Author: jupfi Date: Wed May 29 08:36:34 2024 +0200 Initial commit. diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..f9af7dd --- /dev/null +++ b/.gitignore @@ -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 \ No newline at end of file diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..7198416 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,5 @@ +# Changelog + +## Version 0.0.1 (15-04-2024) + +- Initial release diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..7617f05 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/README.md b/README.md new file mode 100644 index 0000000..0fb51ba --- /dev/null +++ b/README.md @@ -0,0 +1 @@ +# quackseq-simulator \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..1ed1970 --- /dev/null +++ b/pyproject.toml @@ -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"] \ No newline at end of file diff --git a/src/quackseq_simulator/__init__.py b/src/quackseq_simulator/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/quackseq_simulator/simulator.py b/src/quackseq_simulator/simulator.py new file mode 100644 index 0000000..9462b4e --- /dev/null +++ b/src/quackseq_simulator/simulator.py @@ -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 diff --git a/src/quackseq_simulator/simulator_controller.py b/src/quackseq_simulator/simulator_controller.py new file mode 100644 index 0000000..1fccb40 --- /dev/null +++ b/src/quackseq_simulator/simulator_controller.py @@ -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) diff --git a/src/quackseq_simulator/simulator_model.py b/src/quackseq_simulator/simulator_model.py new file mode 100644 index 0000000..7489bb3 --- /dev/null +++ b/src/quackseq_simulator/simulator_model.py @@ -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 diff --git a/tests/simulator.py b/tests/simulator.py new file mode 100644 index 0000000..b8cc4b3 --- /dev/null +++ b/tests/simulator.py @@ -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()