nqrduck/quackseq-simulator
1
0
Fork 0
mirror of https://github.com/nqrduck/quackseq-simulator.git synced 2025-01-02 18:08:07 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

Working PhaseCycling.

Browse source
This commit is contained in:
jupfi 2024-06-07 16:48:45 +02:00
parent 471c9e52be
commit 0720f19d87
2 changed files with 68 additions and 47 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

99
src/quackseq_simulator/simulator_controller.py
View file

@ -29,53 +29,64 @@ class SimulatorController(SpectrometerController):
This will start the simulation based on the settings and the pulse sequence. This will start the simulation based on the settings and the pulse sequence.
""" """
logger.debug("Starting simulation") logger.debug("Starting simulation")
sample = self.get_sample_from_settings()
logger.debug("Sample: %s", sample.name)
# This needs to be called to update the phase array # This needs to be called to update the phase array
sequence.phase_table.update_phase_array() sequence.phase_table.update_phase_array()
dwell_time = self.calculate_dwelltime(sequence) # Get the number of phasecycles
logger.debug("Dwell time: %s", dwell_time) number_phasecycles = sequence.phase_table.n_phase_cycles
try: # Empty measurement object
pulse_array = self.translate_pulse_sequence(sequence, dwell_time) measurement_data = None
except AttributeError:
logger.warning("Could not translate pulse sequence")
return
simulation = self.get_simulation(sample, pulse_array) for cycle in range(number_phasecycles):
result = simulation.simulate() sample = self.get_sample_from_settings()
logger.debug("Sample: %s", sample.name)
tdx = ( dwell_time = self.calculate_dwelltime(sequence)
np.linspace( logger.debug("Dwell time: %s", dwell_time)
0, float(self.calculate_simulation_length(sequence)), len(result)
try:
pulse_array = self.translate_pulse_sequence(sequence, dwell_time, cycle)
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
) )
* 1e6
)
rx_begin, rx_stop, phase = self.translate_rx_event(sequence) rx_begin, rx_stop, phase = self.translate_rx_event(sequence)
# If we have a RX event, we need to cut the result to the RX event # If we have a RX event, we need to cut the result to the RX event
if rx_begin and rx_stop: if rx_begin and rx_stop:
evidx = np.where((tdx > rx_begin) & (tdx < rx_stop))[0] evidx = np.where((tdx > rx_begin) & (tdx < rx_stop))[0]
tdx = tdx[evidx] tdx = tdx[evidx]
result = result[evidx] result = result[evidx]
# Measurement name date + module + target frequency + averages + sequence name # Measurement name date + module + target frequency + averages + sequence name
name = f"{datetime.now().strftime('%Y-%m-%d %H:%M:%S')} - Simulator - {self.simulator.model.target_frequency / 1e6} MHz - {self.simulator.model.averages} averages - {sequence.name}" name = f"{datetime.now().strftime('%Y-%m-%d %H:%M:%S')} - Simulator - {self.simulator.model.target_frequency / 1e6} MHz - {self.simulator.model.averages} averages - {sequence.name}"
logger.debug(f"Measurement name: {name}") logger.debug(f"Measurement name: {name}")
measurement_data = Measurement( if not measurement_data:
name, measurement_data = Measurement(
tdx, name,
result / simulation.averages, tdx,
sample.resonant_frequency, result / simulation.averages,
# frequency_shift=self.module.model.if_frequency, sample.resonant_frequency,
) )
else:
measurement_data.add_dataset(tdx, result / simulation.averages)
if phase: if phase:
measurement_data.phase_shift(phase, 0) measurement_data.phase_shift(phase, cycle)
return measurement_data return measurement_data
@ -127,13 +138,14 @@ class SimulatorController(SpectrometerController):
return sample return sample
def translate_pulse_sequence( def translate_pulse_sequence(
self, sequence: QuackSequence, dwell_time: float self, sequence: QuackSequence, dwell_time: float, cycle: int
) -> PulseArray: ) -> PulseArray:
"""This method translates the pulse sequence from the core to a PulseArray object needed for the simulation. """This method translates the pulse sequence from the core to a PulseArray object needed for the simulation.
Args: Args:
sequence (QuackSequence): The pulse sequence from the core. sequence (QuackSequence): The pulse sequence from the core.
dwell_time (float): The dwell time in seconds. dwell_time (float): The dwell time in seconds.
cycle (int): The current phase cycle.
Returns: Returns:
PulseArray: The pulse sequence translated to a PulseArray object. PulseArray: The pulse sequence translated to a PulseArray object.
@ -143,6 +155,12 @@ class SimulatorController(SpectrometerController):
amplitude_array = list() amplitude_array = list()
phase_array = list() phase_array = list()
# Get the phase_table
phase_table = sequence.phase_table.phase_array[cycle]
# Count the number of TX pulses with relative amplitude > 0
n_tx_pulses = 0
for event in events: for event in events:
logger.debug("Event %s has parameters: %s", event.name, event.parameters) logger.debug("Event %s has parameters: %s", event.name, event.parameters)
for parameter in event.parameters.values(): for parameter in event.parameters.values():
@ -156,6 +174,7 @@ class SimulatorController(SpectrometerController):
> 0 > 0
): ):
logger.debug(f"Adding pulse: {event.duration} s") logger.debug(f"Adding pulse: {event.duration} s")
# If we have a pulse, we need to add it to the pulse array # If we have a pulse, we need to add it to the pulse array
pulse_shape = parameter.get_option_by_name( pulse_shape = parameter.get_option_by_name(
TXPulse.TX_PULSE_SHAPE TXPulse.TX_PULSE_SHAPE
@ -167,12 +186,13 @@ class SimulatorController(SpectrometerController):
) )
amplitude_array.append(pulse_amplitude) amplitude_array.append(pulse_amplitude)
# Phase stuff # Phase from the phase table - column is the number of the pulse
phase = parameter.get_option_by_name(TXPulse.TX_PHASE).value phase = phase_table[n_tx_pulses]
# Degrees to radians # Degrees to radians
phase = np.radians(phase) phase = np.radians(phase)
phase_array.append([phase for _ in range(len(pulse_amplitude))]) phase_array.append([phase for _ in range(len(pulse_amplitude))])
n_tx_pulses += 1
elif ( elif (
parameter.name == sequence.TX_PULSE parameter.name == sequence.TX_PULSE
@ -182,10 +202,9 @@ class SimulatorController(SpectrometerController):
# If we have a wait, we need to add it to the pulse array # 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.append(np.zeros(int(event.duration / dwell_time)))
phase_array.append(np.zeros(int(event.duration / dwell_time))) phase_array.append(np.zeros(int(event.duration / dwell_time)))
amplitude_array = np.concatenate(amplitude_array) amplitude_array = np.concatenate(amplitude_array)
phase_array = np.concatenate(phase_array) phase_array = np.concatenate(phase_array)
pulse_array = PulseArray( pulse_array = PulseArray(
pulseamplitude=amplitude_array, pulseamplitude=amplitude_array,

16
tests/simulator.py
View file

@ -93,14 +93,14 @@ class TestQuackSequence(unittest.TestCase):
seq.add_event(tx2) seq.add_event(tx2)
seq.set_tx_amplitude(tx2, 100) seq.set_tx_amplitude(tx2, 100)
seq.set_tx_phase(tx2, 1) seq.set_tx_phase(tx2, 1)
seq.set_tx_n_phase_cycles(tx2, 4) seq.set_tx_n_phase_cycles(tx2, 2)
seq.set_tx_phase_cycle_group(tx2, 1) seq.set_tx_phase_cycle_group(tx2, 1)
tx3 = Event("tx3", "10u", seq) tx3 = Event("tx3", "10u", seq)
seq.add_event(tx3) seq.add_event(tx3)
seq.set_tx_amplitude(tx3, 100) seq.set_tx_amplitude(tx3, 100)
seq.set_tx_phase(tx3, 2) seq.set_tx_phase(tx3, 2)
seq.set_tx_n_phase_cycles(tx3, 3) seq.set_tx_n_phase_cycles(tx3, 2)
seq.set_tx_phase_cycle_group(tx3, 3) seq.set_tx_phase_cycle_group(tx3, 3)
sim = Simulator() sim = Simulator()
@ -109,12 +109,13 @@ class TestQuackSequence(unittest.TestCase):
result = sim.run_sequence(seq) result = sim.run_sequence(seq)
plt.plot(result.tdx[0], abs(result.tdy[0])) plt.plot(result.tdx[0], abs(result.tdy[0]))
plt.plot(result.tdx[1], abs(result.tdy[1]))
plt.show() plt.show()
phase_table = PhaseTable(seq) phase_table = PhaseTable(seq)
logger.info(phase_table.n_phase_cycles) logger.info(phase_table.n_phase_cycles)
self.assertEqual(phase_table.n_phase_cycles, 24) self.assertEqual(phase_table.n_phase_cycles, 8)
self.assertEqual(phase_table.n_parameters, 3) self.assertEqual(phase_table.n_parameters, 3)
def test_phase_cycling(self): def test_phase_cycling(self):
@ -152,12 +153,13 @@ class TestQuackSequence(unittest.TestCase):
phase_table.rx_phase_sign = readout_scheme phase_table.rx_phase_sign = readout_scheme
result = Simulator().run_sequence(seq) result = Simulator().run_sequence(seq)
plt.title("Phase cycling") plt.title("Phase cycling")
plt.plot(result.tdx[0], abs(result.tdy[0]), label="abs") logger.info(f"Number of data sets {len(result.tdy)}")
plt.plot(result.tdx[0], result.tdy[0].real, label="real") plt.plot(result.tdx[0], result.tdy[0].real, label="pc 1")
plt.plot(result.tdx[0], result.tdy[0].imag, label="imag") plt.plot(result.tdx[1], result.tdy[1].real, label="pc 2")
plt.plot(result.tdx[2], result.tdy[2].real, label="pc 3")
plt.plot(result.tdx[3], result.tdy[3].real, label="pc 4")
plt.legend() plt.legend()
plt.show() plt.show()
# rx = Event("rx", "100u", seq) # rx = Event("rx", "100u", seq)