From 774c383f34968eb836aaac2a7ed7f5b5f1c1ada6 Mon Sep 17 00:00:00 2001
From: SeanvdMeer <18538762+MiniSean@users.noreply.github.com>
Date: Fri, 2 Feb 2024 13:45:00 +0100
Subject: [PATCH] Working distance-3 repetition code conditions (Z4-Z2).
 Includes 0 round QEC-cycles

---
 pycqed/analysis_v2/GBT_analysis.py            |  72 +-
 .../analysis_v2/Parity_benchmark_analysis.py  |   8 +-
 pycqed/analysis_v2/Two_qubit_gate_analysis.py |  39 +-
 pycqed/analysis_v2/cryoscope_v2_analysis.py   |   8 +-
 .../Surface17_dependency_graph.py             | 378 +++++-----
 .../meta_instrument/device_object_CCL.py      | 237 ++++++-
 .../flux_crosstalk_ac/measurement_function.py | 230 ++++++-
 .../measurement/flux_crosstalk_ac/schedule.py |  47 +-
 .../config_cc_s17_direct_iq.json.in           | 644 ++++++++++++++++--
 .../openql_experiments/multi_qubit_oql.py     | 277 +++++++-
 .../connectivity_surface_code.py              | 284 +++++++-
 .../intrf_channel_identifier.py               |  12 +
 12 files changed, 1809 insertions(+), 427 deletions(-)

diff --git a/pycqed/analysis_v2/GBT_analysis.py b/pycqed/analysis_v2/GBT_analysis.py
index e392f1ff1..e08893193 100644
--- a/pycqed/analysis_v2/GBT_analysis.py
+++ b/pycqed/analysis_v2/GBT_analysis.py
@@ -1,4 +1,5 @@
 import os
+import warnings
 import matplotlib.pyplot as plt
 import numpy as np
 import pycqed.analysis_v2.base_analysis as ba
@@ -177,6 +178,17 @@ def extract_data(self):
 			self.raw_data_dict[q]['t_T2'] = t2.sweep_points[:-4]
 			self.raw_data_dict[q]['p_T2'] = t2.normalized_data_points
 			self.raw_data_dict[q]['T2'] = t2.fit_res.params['tau'].value
+			# Bundle all errors in array so it gets saved in hdf5
+			self.raw_data_dict[q]['Errors'] = np.array([
+				self.raw_data_dict[q]['SQG_err'],
+				self.raw_data_dict[q]['L1_err'],
+				self.raw_data_dict[q]['RO_err'],
+				self.raw_data_dict[q]['allxy_err'],
+				self.raw_data_dict[q]['T1'],
+				self.raw_data_dict[q]['T2']])
+
+		self.proc_data_dict = {}
+		self.proc_data_dict['quantities_of_interest'] = self.raw_data_dict
 
 	def process_data(self):
 		pass
@@ -955,30 +967,42 @@ def extract_data(self):
 			Data_qubits = fp_tomo.split(f'tomography_{stab}_')[-1]
 			Data_qubits = Data_qubits.split('_sim')[0]
 			Data_qubits = Data_qubits.split('_')
-			if stab == 'Z1':
-				exc_qubits = ['D1']
-			elif stab == 'Z2':
-				exc_qubits = ['D3']
-			elif stab == 'Z3':
-				exc_qubits = []
-			elif stab == 'Z4':
-				exc_qubits = ['D5']
-			else:
-			    raise ValueError('exception qubits must be given to analysis')
-			a = ma2.pba.Weight_n_parity_tomography(
-				label = label_tomo,
-				sim_measurement=True,
-				n_rounds=2,
-				exception_qubits=exc_qubits,
-				post_selection=False,
-				extract_only=True)
-			self.raw_data_dict[stab]['rho_0'] = a.proc_data_dict['rho_0']
-			self.raw_data_dict[stab]['rho_1'] = a.proc_data_dict['rho_1']
-			self.raw_data_dict[stab]['fidelity_0'] = a.proc_data_dict['Fid_0']
-			self.raw_data_dict[stab]['fidelity_1'] = a.proc_data_dict['Fid_1']
-			self.raw_data_dict[stab]['M1'] = a.proc_data_dict['M1']
-			self.raw_data_dict[stab]['M2'] = a.proc_data_dict['M2']
-			self.raw_data_dict[stab]['repeatability'] = a.proc_data_dict['repeatability']
+			try:
+				if stab == 'Z1':
+					exc_qubits = ['D4', 'D5']  # ['D1']
+				elif stab == 'Z2':
+					exc_qubits = ['D3']
+				elif stab == 'Z3':
+					exc_qubits = ['D7']
+				elif stab == 'Z4':
+					exc_qubits = ['D5']
+				elif stab == 'X1':
+					exc_qubits = ['D2']
+				elif stab == 'X2':
+					exc_qubits = ['D2', 'D3']
+				elif stab == 'X3':
+					exc_qubits = ['D8']
+				elif stab == 'X4':
+					exc_qubits = []
+				else:
+					raise ValueError('exception qubits must be given to analysis')
+
+				a = ma2.pba.Weight_n_parity_tomography(
+					label = label_tomo,
+					sim_measurement=True,
+					n_rounds=2,
+					exception_qubits=exc_qubits,
+					post_selection=False,
+					extract_only=True)
+				self.raw_data_dict[stab]['rho_0'] = a.proc_data_dict['rho_0']
+				self.raw_data_dict[stab]['rho_1'] = a.proc_data_dict['rho_1']
+				self.raw_data_dict[stab]['fidelity_0'] = a.proc_data_dict['Fid_0']
+				self.raw_data_dict[stab]['fidelity_1'] = a.proc_data_dict['Fid_1']
+				self.raw_data_dict[stab]['M1'] = a.proc_data_dict['M1']
+				self.raw_data_dict[stab]['M2'] = a.proc_data_dict['M2']
+				self.raw_data_dict[stab]['repeatability'] = a.proc_data_dict['repeatability']
+			except Exception as e:
+				warnings.warn(e)
 
 	def process_data(self):
 		pass
diff --git a/pycqed/analysis_v2/Parity_benchmark_analysis.py b/pycqed/analysis_v2/Parity_benchmark_analysis.py
index ea0271e64..5e2dc3eec 100644
--- a/pycqed/analysis_v2/Parity_benchmark_analysis.py
+++ b/pycqed/analysis_v2/Parity_benchmark_analysis.py
@@ -1901,7 +1901,7 @@ def defect_rate_plotn(
     axs[0].set_ylabel('error probability')
     axs[0].set_xlabel('rounds')
     axs[0].legend(frameon=False, bbox_to_anchor = (1.01, 1))
-    axs[0].set_title('Deffect rate')
+    axs[0].set_title('defect rate')
 
     axs[1].plot((np.arange(n_rounds)+1), p_0*100, 'C0-', label='Normal')
     axs[1].plot((np.arange(n_rounds)+1), p_1*100, 'C1-', label='LRU data')
@@ -2430,11 +2430,13 @@ def _get_pop_vector(Shots):
             # select experiments for which the Pij matrix will be computed
             k_of_interest = []
             if ('surface_13' in self.experiments) or \
-               ('surface_13_LRU' in self.experiments):
+               ('surface_13_LRU' in self.experiments) or True:
                 if ('surface_13' in self.experiments):
                     k_of_interest.append(self.experiments.index('surface_13'))
                 if ('surface_13_LRU' in self.experiments):
                     k_of_interest.append(self.experiments.index('surface_13_LRU'))
+                else:
+                    k_of_interest = [0]
                 # Calculate defects each stabilizer ancilla qubits
                 _Ancilla_qubits = [q for q in self.Qubits if 'Z' in q]
                 Defects = { q : { k:{} for k in range(n_kernels) } for q in _Ancilla_qubits }
@@ -2598,7 +2600,7 @@ def defect_rate_k_plotfn(
     axs[0].grid(ls='--')
     axs[0].set_ylabel('error probability')
     axs[0].set_xlabel('rounds')
-    axs[0].set_title('Deffect rate')
+    axs[0].set_title('defect rate')
     axs[0].set_ylim(0, .5)
     # axs[0].set_ylim(0, 1)
     axs[0].set_yticks([0, .1, .2, .3, .4, .5])
diff --git a/pycqed/analysis_v2/Two_qubit_gate_analysis.py b/pycqed/analysis_v2/Two_qubit_gate_analysis.py
index f92f21d74..a78a23de8 100644
--- a/pycqed/analysis_v2/Two_qubit_gate_analysis.py
+++ b/pycqed/analysis_v2/Two_qubit_gate_analysis.py
@@ -364,7 +364,7 @@ def get_plot_axis(vals, rang=None):
     # Frequency qubit population
     vmax = min([1, np.max(Pop)])
     vmax = max([vmax, 0.15])
-    im = ax.pcolormesh(Detunings*1e-6, Times*1e9, Pop, vmax=vmax)
+    im = ax.pcolormesh(Detunings*1e-6, Times*1e9, Pop, vmax=1)#vmax)
     fig.colorbar(im, ax=ax, label='Population')
     # plot two-qubit gate frequencies:
     if interaction_freqs:
@@ -3365,7 +3365,8 @@ def CZ_frequency_trajectory_plotfn(
         axR.set_position([pos.x0+pos.width*1.005, pos.y0, pos.width*0.2, pos.height])
         def get_plot_axis(vals, rang=None):
             if len(vals)>1:
-                dx = vals[1]-vals[0]
+                n = len(vals)//2
+                dx = vals[n]-vals[n-1]
                 X = np.concatenate((vals, [vals[-1]+dx])) - dx/2
             else:
                 X = vals
@@ -3375,25 +3376,34 @@ def get_plot_axis(vals, rang=None):
         Pop = TLS_analysis_dict[qH].proc_data_dict['Pop'] 
         # Frequency qubit population
         vmax = min([1, np.max(Pop)])
-        vmax = max([vmax, 0.15])
+        vmax = 1#max([vmax, 0.15])
         im = axR.pcolormesh(Times*1e9, Detunings*1e-9, Pop.transpose(), vmax=vmax)
         axR.text(Times[len(Times)//2]*1e9, Detunings[0]*1e-9-.05, qH, ha='center', va='top', color='w')
-        axR.set_title('Gate qubits', size=7)
+        axR.set_title('High qubit', size=7)
         if qL in TLS_analysis_dict.keys():
+            axL = fig.add_subplot(221)
+            # using previous axis position <pos>
+            axL.set_position([pos.x0+pos.width*(1.21), pos.y0,
+                              pos.width*0.2, pos.height])
             Detunings = data[qL]['frequency'] - get_plot_axis(TLS_analysis_dict[qL].proc_data_dict['Detunings'])
-            Pop = TLS_analysis_dict[qL].proc_data_dict['Pop'] 
+            Pop = TLS_analysis_dict[qL].proc_data_dict['Pop']
             # Frequency qubit population
             vmax = min([1, np.max(Pop)])
-            vmax = max([vmax, 0.15])
-            im = axR.pcolormesh(Times*1e9, Detunings*1e-9, Pop.transpose(), vmax=vmax)
-            axR.text(Times[len(Times)//2]*1e9, Detunings[0]*1e-9-.05, qL, ha='center', va='top', color='w')
-            axR.axhline(Detunings[0]*1e-9, color='w')
+            vmax = 1#max([vmax, 0.15])
+            im = axL.pcolormesh(Times*1e9, Detunings*1e-9, Pop.transpose(), vmax=vmax)
+            axL.text(Times[len(Times)//2]*1e9, max(Detunings)*1e-9-.05, qL, ha='center', va='top', color='w')
+            # axR.axhline(max(Detunings)*1e-9, color='w')
+            axL.set_title('Low qubit', size=7)
+            axL.set_ylim(ax.get_ylim())
+            axL.yaxis.tick_right()
+            axL.set_xticks([])
+            axL.axis('off')
         axR.set_ylim(ax.get_ylim())
         axR.yaxis.tick_right()
         axR.set_xticks([])
         axR.axis('off')
         # Parked qubit plots
-        i = 0
+        i = 1
         for q in parked_qubits:
             if q in TLS_analysis_dict.keys():
                 axP = fig.add_subplot(221+i)
@@ -3405,11 +3415,10 @@ def get_plot_axis(vals, rang=None):
                 Pop = TLS_analysis_dict[q].proc_data_dict['Pop'] 
                 # Frequency qubit population
                 vmax = min([1, np.max(Pop)])
-                vmax = max([vmax, 0.15])
+                vmax = 1#max([vmax, 0.15])
                 im = axP.pcolormesh(Times*1e9, Detunings*1e-9, Pop.transpose(), vmax=vmax)
-                axP.text(Times[len(Times)//2]*1e9, Detunings[0]*1e-9-.05, q, ha='center', va='top', color='w')
-                
-                axP.set_title('Park qubits', size=7)
+                axP.text(Times[len(Times)//2]*1e9, max(Detunings)*1e-9-.05, q, ha='center', va='top', color='w')
+                axP.set_title('Park qubit', size=7)
                 axP.set_ylim(ax.get_ylim())
                 axP.yaxis.tick_right()
                 axP.set_xticks([])
@@ -3552,8 +3561,10 @@ def func(x, phi, A, B):
         # self.qoi['L_1'] = L_1
         self.qoi['P_excited'] = P_excited
         self.qoi['Phases'] = {}
+        self.qoi['Contrast'] = {}
         for q in self.Q_target:
             self.qoi['Phases'][q] = { c:Fit_res[q][c][0] for c in self.control_cases }
+            self.qoi['Contrast'][q] = { c:Fit_res[q][c][1] for c in self.control_cases }
         if self.solve_for_phase_gate_model:
             self.qoi['Phase_model'] = Phase_model
 
diff --git a/pycqed/analysis_v2/cryoscope_v2_analysis.py b/pycqed/analysis_v2/cryoscope_v2_analysis.py
index 4029e7f36..5e3f270e6 100644
--- a/pycqed/analysis_v2/cryoscope_v2_analysis.py
+++ b/pycqed/analysis_v2/cryoscope_v2_analysis.py
@@ -1,8 +1,6 @@
 """
 Created: 2020-07-15
-Author: Victor Negirneac
 """
-
 import matplotlib.pyplot as plt
 from pycqed.analysis.analysis_toolbox import get_datafilepath_from_timestamp
 import pycqed.analysis_v2.cryoscope_v2_tools as cv2_tools
@@ -738,7 +736,7 @@ def __init__(self,
                  extract_only: bool = False,
                  auto=True,
                  poly_params: dict = None,
-                 derivative_window_length: float=3e-9,
+                 derivative_window_length: float=5e-9,
                 ):
         super().__init__(t_start=t_start, t_stop=t_stop,
                          label=label,
@@ -1265,8 +1263,8 @@ def skewed_gauss(x, x0, sigma, alpha, a, b):
         # Fit exponential to trace
         if self.update_IIR:
             try:
-                p0 = [-.01, 100e-9, 1.0085]
-                popt, pcov = curve_fit(filter_func, Time[:]*1e-9, Trace[:], p0=p0)
+                p0 = [+.001, 400e-9, 1.0085]
+                popt, pcov = curve_fit(filter_func, Time[8:]*1e-9, Trace[8:], p0=p0)
             except:
                 print_exception()
                 print('Fit failed. Trying new initial guess')
diff --git a/pycqed/instrument_drivers/meta_instrument/Surface17_dependency_graph.py b/pycqed/instrument_drivers/meta_instrument/Surface17_dependency_graph.py
index ba97fbb0a..e97870074 100644
--- a/pycqed/instrument_drivers/meta_instrument/Surface17_dependency_graph.py
+++ b/pycqed/instrument_drivers/meta_instrument/Surface17_dependency_graph.py
@@ -179,18 +179,18 @@ def create_dep_graph(self):
 			['Z2', 'D3'], 
 			['Z1', 'D1'], 
 			['D5', 'Z4'],
-			# ['X1', 'D2'], 
-			# ['D6', 'X2'], 
-			# ['X3', 'D8'],
-			# ['X1', 'D1'], 
-			# ['D5', 'X2'], 
-			# ['X3', 'D7'],
-			# ['X4', 'D9'], 
-			# ['D5', 'X3'], 
-			# ['X2', 'D3'],
-			# ['X4', 'D8'], 
-			# ['D4', 'X3'], 
-			# ['X2', 'D2'],
+			['X1', 'D2'], 
+			['D6', 'X2'], 
+			['X3', 'D8'],
+			['X1', 'D1'], 
+			['D5', 'X2'], 
+			['X3', 'D7'],
+			['X4', 'D9'], 
+			['D5', 'X3'], 
+			['X2', 'D3'],
+			['X4', 'D8'], 
+			['D4', 'X3'], 
+			['X2', 'D2'],
 			]
 			# Single-qubit nodes
 		Qubits = np.unique(np.array(Qubit_pairs).flatten())
@@ -203,12 +203,15 @@ def create_dep_graph(self):
 		          	})
 		# Two-qubit nodes
 		QL_detunings = {
-			('Z1', 'D2') : 250e6,#400e6,
-			('Z1', 'D1') : 400e6,
+			# After detuning search (Yuejie&Sean)
+			('D5', 'X2') : 60e6,
+			('D5', 'X3') : 110e6,
+			('Z1', 'D1') : 380e6,
+			('Z1', 'D2') : 310e6,
+			('X1', 'D1') : 0e6,
+			('X3', 'D8') : 60e6,
 			('Z4', 'D8') : 100e6,
-			# ('Z4', 'D9') : 100e6,
-			('X3', 'D7') : 100e6,
-			('X3', 'D8') : 100e6,
+			('X1', 'D2') : 165e6,	# After second round search
 		}
 		for pair in Qubit_pairs:
 			self.add_node(f'{pair[0]}, {pair[1]} Chevron',
@@ -697,6 +700,12 @@ def SSRO_wrapper(qubit:str, station):
 	Returns True if successful calibration otherwise
 	returns False.
 	'''
+	file_cfg = gc.generate_config(in_filename=input_file,
+                                  out_filename=config_fn,
+                                  mw_pulse_duration=20,
+                                  ro_duration=420,
+                                  flux_pulse_duration=40,
+                                  init_duration=200000)
 	station.components['MC'].live_plot_enabled(False)
 	station.components['nested_MC'].live_plot_enabled(False)
 	# Set initial parameters for calibration
@@ -884,6 +893,7 @@ def drive_mixer_wrapper(qubit:str, station):
                        'openql_experiments', 'output_cc_s17','config_cc_s17_direct_iq.json')
 
 TWOQ_GATE_DURATION = 60e-9
+TWOQ_GATE_DURATION_NS = 60
 
 OFFSET_QUBITS = ['X2', 'X3', 'X4', 'D7', 'D9']
 
@@ -1131,7 +1141,7 @@ def Flux_arc_wrapper(Qubit, station,
 	for detuning to voltage conversion.
 	'''
 	# Set gate duration
-	TQG_duration_ns = int(TWOQ_GATE_DURATION*1e9) 
+	TQG_duration_ns = TWOQ_GATE_DURATION_NS 
 	file_cfg = gc.generate_config(in_filename=input_file,
 	                              out_filename=config_fn,
 	                              mw_pulse_duration=20,
@@ -1224,7 +1234,7 @@ def Chevron_wrapper(qH, qL, station,
 					   qubit. 
 	'''
 	# Set gate duration
-	TQG_duration_ns = int(TWOQ_GATE_DURATION*1e9)
+	TQG_duration_ns = TWOQ_GATE_DURATION_NS
 	file_cfg = gc.generate_config(in_filename=input_file,
 	                              out_filename=config_fn,
 	                              mw_pulse_duration=20,
@@ -1313,13 +1323,14 @@ def Chevron_wrapper(qH, qL, station,
 		('Z2', 'D3'): 2.1666666666666665e-08+1/2.4e9,
 		('Z1', 'D1'): 2.5416666666666666e-08,
 		('D5', 'Z4'): 1.875e-08,
-		('X1', 'D1'): 2e-08,
+		('X1', 'D1'): 2.0833333333333335e-08,  # (48 sampling points) -> Go to 50 sampling points (2.0833333333333335e-08)
 		('X1', 'D2'): 2.2083333333333333e-08+2/2.4e9,
-		('D5', 'X2'): 1.875e-08,
+		('D5', 'X2'): 1.875e-08-2/2.4e9,
 		('D6', 'X2'): 1.9583333333333333e-08-1/2.4e9,
 		('D4', 'X3'): 2.2083333333333333e-08,
-		# ('X2', 'D2'): 2.4166666666666668e-08+4/2.4e9,
-		# ('X2', 'D3'): 2.0416666666666668e-08+8/2.4e9,
+		('D5', 'X3'): 2.0416666666666668e-08,
+		('X2', 'D2'): 2.0833333333333335e-08-2/2.4e9,
+		('X2', 'D3'): 1.9583333333333333e-08-1/2.4e9,
 		('X3', 'D7'): 2.0416666666666668e-08-1/2.4e9,
 		('X3', 'D8'): 2.1666666666666665e-08-1/2.4e9,
 		('X4', 'D8'): 2.0833333333333335e-08,
@@ -1393,7 +1404,7 @@ def SNZ_tmid_wrapper(qH, qL, station,
 					   parked qubits to low-frequency
 	'''
 	# Set gate duration
-	TQG_duration_ns = int(TWOQ_GATE_DURATION*1e9) - 20
+	TQG_duration_ns = TWOQ_GATE_DURATION_NS - 20
 	file_cfg = gc.generate_config(in_filename=input_file,
 	                              out_filename=config_fn,
 	                              mw_pulse_duration=20,
@@ -1556,7 +1567,7 @@ def SNZ_AB_wrapper(qH, qL, station,
 					   parked qubits to low-frequency
 	'''
 	# Set gate duration
-	TQG_duration_ns = int(TWOQ_GATE_DURATION*1e9) - 20
+	TQG_duration_ns = TWOQ_GATE_DURATION_NS - 20
 	file_cfg = gc.generate_config(in_filename=input_file,
 	                              out_filename=config_fn,
 	                              mw_pulse_duration=20,
@@ -1706,7 +1717,7 @@ def Unipolar_wrapper(qH, qL, station,
 					   parked qubits to low-frequency
 	'''
 	# Set gate duration
-	TQG_duration_ns = int(TWOQ_GATE_DURATION*1e9) - 20
+	TQG_duration_ns = TWOQ_GATE_DURATION_NS - 20
 	file_cfg = gc.generate_config(in_filename=input_file,
 	                              out_filename=config_fn,
 	                              mw_pulse_duration=20,
@@ -1810,7 +1821,7 @@ def Asymmetry_wrapper(qH, qL, station):
 	returns True.
 	'''
 	# Set gate duration
-	TQG_duration_ns = int(TWOQ_GATE_DURATION*1e9) - 20
+	TQG_duration_ns = TWOQ_GATE_DURATION_NS - 20
 	file_cfg = gc.generate_config(in_filename=input_file,
 	                              out_filename=config_fn,
 	                              mw_pulse_duration=20,
@@ -1890,20 +1901,21 @@ def Asymmetry_wrapper(qH, qL, station):
 def Single_qubit_phase_calibration_wrapper(qH, qL, station,
 										   park_distance=700e6,
 										   apply_parking_settings: bool = True,
-										   fine_cphase_calibration: bool = True):
+										   fine_cphase_calibration: bool = False,
+										   pc_repetitions = 1):
 	'''
 	Wrapper function for fine-tunig CP 180 phase, SQ phase updates of 360, and verification. 
 	Returns True if successful calibration otherwise
 	returns False.
 	'''
 	# Set gate duration
-	TQG_duration_ns = int(TWOQ_GATE_DURATION*1e9) - 20
-	file_cfg = gc.generate_config(in_filename=input_file,
-	                              out_filename=config_fn,
-	                              mw_pulse_duration=20,
-	                              ro_duration=1000,
-	                              flux_pulse_duration=TQG_duration_ns,
-	                              init_duration=200000)
+	TQG_duration_ns = TWOQ_GATE_DURATION_NS - 20
+	# file_cfg = gc.generate_config(in_filename=input_file,
+	#                               out_filename=config_fn,
+	#                               mw_pulse_duration=20,
+	#                               ro_duration=1000,
+	#                               flux_pulse_duration=TQG_duration_ns,
+	#                               init_duration=200000)
 	# Setup for measurement
 	dircts = get_gate_directions(qH, qL)
 	station.components['MC'].live_plot_enabled(False)
@@ -1945,6 +1957,7 @@ def Single_qubit_phase_calibration_wrapper(qH, qL, station,
 		Q_inst = station.components[q]
 		flux_lm_p = Q_inst.instr_LutMan_Flux.get_instr()
 		park_det = Q_inst.freq_qubit()-park_freq
+		print('park_det', park_det)
 		if apply_parking_settings:
 			# Only park if the qubit is closer than then 350 MHz
 			if park_det>20e6:
@@ -1956,6 +1969,7 @@ def Single_qubit_phase_calibration_wrapper(qH, qL, station,
 				amp_park_idx = np.argmax(np.abs(_Amps))
 				# Update parking amplitude in lookup table
 				flux_lm_p.park_amp(_Amps[amp_park_idx])
+				print('park_amp', flux_lm_p.park_amp())
 			else:
 				flux_lm_p.park_amp(0)
 		load_single_waveform_on_HDAWG(flux_lm_p, 'park')
@@ -1987,8 +2001,13 @@ def Single_qubit_phase_calibration_wrapper(qH, qL, station,
 		for q in Parked_qubits:
 			flux_lm_p = Q_inst.instr_LutMan_Flux.get_instr()
 			check_flux_wf_duration(flux_lm_p)
+
+
 	# device.prepare_for_timedomain(qubits=[qH, qL])
 	# Check if mw phase pulses are uploaded
+
+
+	# prepare_for_parity_check('X4', station, Data_qubits=['D8', 'D9'])
 	for q in [qH, qL]:
 		Q = station.components[q]
 		Q._prep_td_sources()
@@ -2005,6 +2024,7 @@ def Single_qubit_phase_calibration_wrapper(qH, qL, station,
 		    Q_ancilla = [qH],
 		    Q_control = [qL],
 		    Q_pair_target = [qH, qL],
+		    # flux_cw_list = ['repetition_code_3', 'repetition_code_4'],
 		    flux_cw_list = [flux_cw],
 		    downsample_angle_points = 3,
 		    prepare_for_timedomain = False,
@@ -2016,22 +2036,27 @@ def Single_qubit_phase_calibration_wrapper(qH, qL, station,
 	device.measure_parity_check_ramsey(
 	    Q_target = [qH],
 	    Q_control = [qL],
+	    # flux_cw_list = ['repetition_code_3', 'repetition_code_4'],
 	    flux_cw_list = [flux_cw],
 	    prepare_for_timedomain = False,
 	    downsample_angle_points = 3,
 	    update_mw_phase=True,
 	    mw_phase_param=f'vcz_virtual_q_ph_corr_{dircts[0]}',
-	    disable_metadata=True)
+	    disable_metadata=True,
+	    pc_repetitions=pc_repetitions)
 	# Calibrate low frequency qubit phase
 	device.measure_parity_check_ramsey(
 	    Q_target = [qL],
 	    Q_control = [qH],
+	    # flux_cw_list = ['repetition_code_3', 'repetition_code_4'],
 	    flux_cw_list = [flux_cw],
 	    prepare_for_timedomain = False,
 	    downsample_angle_points = 3,
 	    update_mw_phase=True,
 	    mw_phase_param=f'vcz_virtual_q_ph_corr_{dircts[1]}',
-	    disable_metadata=True)
+	    disable_metadata=True,
+	    pc_repetitions=pc_repetitions)
+
 	mw_lutman_H.upload_single_qubit_phase_corrections()
 	mw_lutman_L.upload_single_qubit_phase_corrections()
 	return True
@@ -2044,7 +2069,7 @@ def TwoQ_Randomized_benchmarking_wrapper(qH, qL, station, **kw):
 	returns False.
 	'''
 	# Set gate duration
-	TQG_duration_ns = int(TWOQ_GATE_DURATION*1e9) - 20
+	TQG_duration_ns = TWOQ_GATE_DURATION_NS - 20
 	# Buffer time after gate
 	if 'buffer_time_ns' in kw.keys():
 		buffer_time_ns = kw['buffer_time_ns']
@@ -2125,7 +2150,8 @@ def TwoQ_Randomized_benchmarking_wrapper(qH, qL, station, **kw):
 
 def TLS_density_wrapper(Qubit, station,
 		                detuning = None,
-		                max_duration = 120e-9):
+		                max_duration = 120e-9,
+						use_second_excited_state: bool = False):
 	'''
 	Wrapper function for measurement of TLS density.
 	Using a dynamical square pulse to flux the qubit
@@ -2158,7 +2184,7 @@ def TLS_density_wrapper(Qubit, station,
 		delta = int(np.round((max_duration-TWOQ_GATE_DURATION)*1e9/20)*20)
 	else:
 		delta = 0
-	TQG_duration_ns = int(TWOQ_GATE_DURATION*1e9) - 20
+	TQG_duration_ns = TWOQ_GATE_DURATION_NS - 20
 	file_cfg = gc.generate_config(in_filename=input_file,
 								out_filename=config_fn,
 								mw_pulse_duration=20,
@@ -2232,6 +2258,7 @@ def TLS_density_wrapper(Qubit, station,
 	    buffer_time=40e-9,
 	    prepare_for_timedomain=False,
 	    disable_metadata=True,
+		second_excited_state=use_second_excited_state,
 	)
 	# Reset waveform durations
 	if max_duration > TWOQ_GATE_DURATION:
@@ -2273,6 +2300,7 @@ def Parking_experiment_wrapper(qH, qL, qP, station):
 								ro_duration=1000,
 								flux_pulse_duration=TQG_duration_ns,
 								init_duration=200000)
+	TQG_duration_ns = TWOQ_GATE_DURATION_NS - 20
 	# Setup for measurement
 	dircts = get_gate_directions(qH, qL)
 	station.components['MC'].live_plot_enabled(False)
@@ -2576,7 +2604,7 @@ def Prepare_for_parity_check_wrapper(station):
 
 
 def Horizontal_calibration_wrapper(stabilizer_qubit, station,
-								   Data_qubits:list = None,
+								   Q_control: list = None,
 								   flux_cw_list = None,
 								   mw_phase_param = None):
 	'''
@@ -2589,9 +2617,11 @@ def Horizontal_calibration_wrapper(stabilizer_qubit, station,
 	device = station.components['device']
 	device.ro_acq_averages(2**8)
 	device.ro_acq_digitized(False)
-	prepare_for_parity_check(stabilizer_qubit, station)
+	prepare_for_parity_check(stabilizer_qubit, station,
+							 Data_qubits = Q_control)
 	Q_ancilla = [stabilizer_qubit]
-	Q_control = list(get_nearest_neighbors(stabilizer_qubit).keys())
+	if not Q_control:
+		Q_control = list(get_nearest_neighbors(stabilizer_qubit).keys())
 	# Parity check settings
 	if flux_cw_list == None:
 		if 'X' in stabilizer_qubit:
@@ -2600,14 +2630,14 @@ def Horizontal_calibration_wrapper(stabilizer_qubit, station,
 			flux_cw_list = [f'flux_dance_{i}' for i in [5, 6, 7, 8]]
 		flux_cw_list = ['cz' for q in Q_control]
 	if mw_phase_param == None:
-		if 'X' in stabilizer_qubit:
-			mw_phase_param = 'vcz_virtual_q_ph_corr_step_4'
-		else:
-			mw_phase_param = 'vcz_virtual_q_ph_corr_step_8'
+		# if 'X' in stabilizer_qubit:
+		# 	mw_phase_param = 'vcz_virtual_q_ph_corr_step_4'
+		# else:
+		# 	mw_phase_param = 'vcz_virtual_q_ph_corr_step_8'
+		dircts = get_gate_directions(q0=stabilizer_qubit, q1=Q_control[0])
+		mw_phase_param = f'vcz_virtual_q_ph_corr_{dircts[0]}'
 	# Calibrate CZ with each data qubit
-	if not Data_qubits:
-		Data_qubits = Q_control
-	for q in Data_qubits:
+	for q in Q_control:
 		# If high frequency qubit
 		if q in ['D4','D5','D6']:
 			# Order of qubits requires high freq. qubit first
@@ -2638,7 +2668,11 @@ def Horizontal_calibration_wrapper(stabilizer_qubit, station,
 
 
 def Measure_parity_check_phase_wrapper(stabilizer_qubit, station, 
-								 flux_cw_list = None, mw_phase_param = None):
+									   Q_control: list = None,
+									   flux_cw_list = None, 
+									   mw_phase_param = None,
+									   pc_repetitions: int = 1,
+									   ):
 	'''
 	Wrapper function to measure pairty checks 
 	returns True
@@ -2649,10 +2683,12 @@ def Measure_parity_check_phase_wrapper(stabilizer_qubit, station,
 	device = station.components['device']
 	device.ro_acq_averages(2**8)
 	device.ro_acq_digitized(False)
-	prepare_for_parity_check(stabilizer_qubit, station)
+	prepare_for_parity_check(stabilizer_qubit, station,
+							 Data_qubits = Q_control)
 	# Parity check settings
 	Q_ancilla = [stabilizer_qubit]
-	Q_control = list(get_nearest_neighbors(stabilizer_qubit).keys())
+	if not Q_control:
+		Q_control = list(get_nearest_neighbors(stabilizer_qubit).keys())
 	if flux_cw_list == None:
 		if 'X' in stabilizer_qubit:
 			flux_cw_list = [f'flux_dance_{i}' for i in [1, 2, 3, 4]]
@@ -2660,18 +2696,19 @@ def Measure_parity_check_phase_wrapper(stabilizer_qubit, station,
 			flux_cw_list = [f'flux_dance_{i}' for i in [5, 6, 7, 8]]
 		flux_cw_list = ['cz' for q in Q_control]
 	if mw_phase_param == None:
-		if 'X' in stabilizer_qubit:
-			mw_phase_param = 'vcz_virtual_q_ph_corr_step_4'
-		else:
-			mw_phase_param = 'vcz_virtual_q_ph_corr_step_8'
-		# dircts = get_gate_directions(q0 = stabilizer_qubit, q1 = Q_control[0])
-		# mw_phase_param = f'vcz_virtual_q_ph_corr_{dircts[0]}'
+		# if 'X' in stabilizer_qubit:
+		# 	mw_phase_param = 'vcz_virtual_q_ph_corr_step_4'
+		# else:
+		# 	mw_phase_param = 'vcz_virtual_q_ph_corr_step_8'
+		dircts = get_gate_directions(q0 = stabilizer_qubit, q1 = Q_control[0])
+		mw_phase_param = f'vcz_virtual_q_ph_corr_{dircts[0]}'
 	# Measure
-	device.measure_parity_check_ramsey(
+	qoi = device.measure_parity_check_ramsey(
 		Q_target = Q_ancilla,
 		Q_control = Q_control,
 		flux_cw_list = flux_cw_list,
 		mw_phase_param=mw_phase_param,
+		pc_repetitions=pc_repetitions,
 		update_mw_phase=True,
 		prepare_for_timedomain = False,
 		disable_metadata=True)
@@ -2679,11 +2716,15 @@ def Measure_parity_check_phase_wrapper(stabilizer_qubit, station,
 	Q_S = station.components[stabilizer_qubit]
 	mw_lm_q = Q_S.instr_LutMan_MW.get_instr()
 	mw_lm_q.upload_single_qubit_phase_corrections()
-	return True
+	return qoi
 
 
 def Data_qubit_phase_calibration_wrapper(stabilizer_qubit, station, 
-								 flux_cw_list = None, mw_phase_param = None):
+										 Q_data: list = None,
+										 flux_cw_list = None, 
+										 mw_phase_param = None,
+									   	 pc_repetitions: int = 1,
+										 ):
 	'''
 	Wrapper function to calibrate single qubit phases of data-qubits in pairty checks 
 	returns True
@@ -2694,10 +2735,12 @@ def Data_qubit_phase_calibration_wrapper(stabilizer_qubit, station,
 	device = station.components['device']
 	device.ro_acq_averages(2**8)
 	device.ro_acq_digitized(False)
-	prepare_for_parity_check(stabilizer_qubit, station)
+	prepare_for_parity_check(stabilizer_qubit, station,
+							 Data_qubits = Q_data)
 	# Parity check settings
 	Q_ancilla = [stabilizer_qubit]
-	Q_data = list(get_nearest_neighbors(stabilizer_qubit).keys())
+	if not Q_data:
+		Q_data = list(get_nearest_neighbors(stabilizer_qubit).keys())
 	if flux_cw_list == None:
 		if 'X' in stabilizer_qubit:
 			flux_cw_list = [f'flux_dance_{i}' for i in [1, 2, 3, 4]]
@@ -2705,22 +2748,23 @@ def Data_qubit_phase_calibration_wrapper(stabilizer_qubit, station,
 			flux_cw_list = [f'flux_dance_{i}' for i in [5, 6, 7, 8]]
 		flux_cw_list = ['cz' for q in Q_data]
 	if mw_phase_param == None:
-		if 'X' in stabilizer_qubit:
-			mw_phase_param = 'vcz_virtual_q_ph_corr_step_4'
-		else:
-			mw_phase_param = 'vcz_virtual_q_ph_corr_step_8'
-		# mw_phase_param = []
-		# for q in Q_data:
-		# 	dircts = get_gate_directions(q0 = q, q1 = stabilizer_qubit)
-		# 	mw_phase_param.append(f'vcz_virtual_q_ph_corr_{dircts[0]}')
+		# if 'X' in stabilizer_qubit:
+		# 	mw_phase_param = 'vcz_virtual_q_ph_corr_step_4'
+		# else:
+		# 	mw_phase_param = 'vcz_virtual_q_ph_corr_step_8'
+		mw_phase_param = []
+		for q in Q_data:
+			dircts = get_gate_directions(q0 = q, q1 = stabilizer_qubit)
+			mw_phase_param.append(f'vcz_virtual_q_ph_corr_{dircts[0]}')
 	# Measure
-	device.measure_parity_check_ramsey(
+	qoi = device.measure_parity_check_ramsey(
 		Q_target = Q_data,
 		Q_control = Q_ancilla,
 		flux_cw_list = flux_cw_list,
 		downsample_angle_points = 1,
 		update_mw_phase=True,
 		mw_phase_param=mw_phase_param,
+		pc_repetitions=pc_repetitions,
 		prepare_for_timedomain = False,
 		disable_metadata=True)
 	# upload phase corrections
@@ -2728,22 +2772,17 @@ def Data_qubit_phase_calibration_wrapper(stabilizer_qubit, station,
 		Q_c = station.components[q]
 		mw_lm_q = Q_c.instr_LutMan_MW.get_instr()
 		mw_lm_q.upload_single_qubit_phase_corrections()
-	return True
+	return qoi
 
 
 def Parity_check_fidelity_wrapper(stabilizer_qubit, station, 
+								  Q_data: list = None,
 								  heralded_init = False,
 								  flux_cw_list = None):
 	'''
 	Wrapper function to measure pairty checks 
 	returns True
 	'''
-	file_cfg = gc.generate_config(in_filename=input_file,
-                              out_filename=config_fn,
-                              mw_pulse_duration=20,
-                              ro_duration=500,
-                              flux_pulse_duration=40,
-                              init_duration=200000)
 	station.components['MC'].live_plot_enabled(False)
 	station.components['nested_MC'].live_plot_enabled(False)
 	# Prepare for timedomain of parity check
@@ -2751,10 +2790,11 @@ def Parity_check_fidelity_wrapper(stabilizer_qubit, station,
 	device.ro_acq_averages(2**10)
 	device.ro_acq_digitized(False)
 	device.ro_acq_integration_length(500e-9)
-	prepare_for_parity_check(stabilizer_qubit, station)
+	prepare_for_parity_check(stabilizer_qubit, station, Data_qubits = Q_data)
 	# Parity check settings
 	Q_ancilla = [stabilizer_qubit]
-	Q_data = list(get_nearest_neighbors(stabilizer_qubit).keys())
+	if not Q_data:
+		Q_data = list(get_nearest_neighbors(stabilizer_qubit).keys())
 	if flux_cw_list == None:
 		if 'X' in stabilizer_qubit:
 			flux_cw_list = [f'flux_dance_{i}' for i in [1, 2, 3, 4]]
@@ -2775,6 +2815,7 @@ def Parity_check_fidelity_wrapper(stabilizer_qubit, station,
 
 
 def Parity_check_repeatability_wrapper(stabilizer_qubit, station, 
+									   Q_data: list = None,
 									   sim_measurement = True,
 									   readout_duration_ns = 420,
 									   repetitions = 5,
@@ -2794,8 +2835,10 @@ def Parity_check_repeatability_wrapper(stabilizer_qubit, station,
 	device.ro_acq_averages(2**10)
 	device.ro_acq_digitized(False)
 	device.ro_acq_integration_length(500e-9)
-	prepare_for_parity_check(stabilizer_qubit, station)
-	Q_data = list(get_nearest_neighbors(stabilizer_qubit).keys())
+	prepare_for_parity_check(stabilizer_qubit, station,
+							 Data_qubits = Q_data)
+	if not Q_data:
+		Q_data = list(get_nearest_neighbors(stabilizer_qubit).keys())
 	if flux_cw_list == None:
 		if 'X' in stabilizer_qubit:
 			flux_cw_list = [f'flux_dance_{i}' for i in [1, 2, 3, 4]]
@@ -2830,9 +2873,9 @@ def Surface_13_wrapper(station, log_zero = None,
 	#######################################################
 	# Preparation
 	#######################################################
-	assert (measurement_time_ns-20)%20 == 0, 'Not a valid measurement time!'
+	assert (measurement_time_ns-20)%40 == 0, 'Not a valid measurement time!'
 	# Set configuration
-	TQG_duration_ns = int(TWOQ_GATE_DURATION*1e9) - 20
+	TQG_duration_ns = TWOQ_GATE_DURATION_NS - 20
 	file_cfg = gc.generate_config(in_filename=input_file,
 	                              out_filename=config_fn,
 	                              mw_pulse_duration=20,
@@ -3047,82 +3090,6 @@ def Surface_13_wrapper(station, log_zero = None,
 				extract_only = False)
 
 
-# def Spectator_data_qubits_wrapper(data_qubits, station, 
-# 						  flux_cw_list = None, mw_phase_param = None):
-# 	'''
-# 	Wrapper function to check spectator effects on data-qubits by toggling 
-# 	ancilla qubits
-# 	returns True
-# 	'''
-# 	# Set Parity check duration
-# 	file_cfg = gc.generate_config(in_filename=input_file,
-# 	                              out_filename=config_fn,
-# 	                              mw_pulse_duration=20,
-# 	                              ro_duration=600,
-# 	                              flux_pulse_duration=40,
-# 	                              init_duration=200000)
-	
-# 	device = station.components['device']
-# 	device.ro_acq_weight_type('optimal')
-# 	device.ro_acq_averages(2**10)
-# 	if data_qubits == None:
-# 		data_qubits =  ['D1', 'D2', 'D3',
-# 						'D4', 'D5', 'D6',
-# 						'D7', 'D8', 'D9',]
-	
-# 	if flux_cw_list == None:
-# 		flux_cw_list = [f'flux_dance_{i}' for i in [5, 6, 7, 8]]
-# 	if mw_phase_param == None:
-# 		mw_phase_param = False
-# 	#hardcode for Z-type ancillas
-# 	Q_control = ['Z1', 'Z2', 'Z3', 'Z4']
-# 	# prepare for timedomain
-# 	device.prepare_for_timedomain(qubits=Q_ancilla+Q_control)
-# 	for q in Q_ancilla:
-# 		mw_lm = device.find_instrument(q).instr_LutMan_MW.get_instr()
-# 		mw_lm.set_default_lutmap()
-# 		mw_lm.load_phase_pulses_to_AWG_lookuptable()
-# 	for q in data_qubits:
-# 		device.measure_parity_check_ramsey(
-# 			Q_target = q,
-# 			Q_control = Q_control,
-# 			flux_cw_list = flux_cw_list,
-# 			prepare_for_timedomain = False,
-# 			update_mw_phase=True,
-# 			mw_phase_param=mw_phase_param,
-# 			disable_metadata=True)
-	
-# 	return True
-
-
-# Dictionary for necessary parking for each interaction
-Park_dict = {
-			 ('Z4', 'D9'): ['D8'],
-			 ('Z4', 'D8'): ['D9'],
-			 ('X4', 'D9'): ['D8'],
-			 ('X4', 'D8'): ['D9'],
-			 ('X3', 'D7'): ['D8'],
-			 ('X3', 'D8'): ['D7'],
-			 ('Z3', 'D7'): [],
-			 ('D4', 'Z3'): ['Z1', 'X3'],
-			 ('D4', 'X3'): ['Z1', 'Z3'],
-			 ('D4', 'Z1'): ['Z3', 'X3'],
-			 ('D5', 'X3'): ['Z1', 'Z4', 'X2'],
-			 ('D5', 'Z1'): ['Z4', 'X2', 'X3'],
-			 ('D5', 'Z4'): ['Z1', 'X2', 'X3'],
-			 ('D5', 'X2'): ['Z1', 'Z4', 'X3'],
-			 ('D6', 'Z4'): ['Z2', 'X2'],
-			 ('D6', 'X2'): ['Z4', 'Z2'],
-			 ('D6', 'Z2'): ['Z4', 'X2'],
-			 ('X1', 'D1'): ['D2'],
-			 ('X1', 'D2'): ['D1'],
-			 ('Z1', 'D1'): ['D2'],
-			 ('Z1', 'D2'): ['D1'],
-			 ('X2', 'D2'): ['D3'],
-			 ('X2', 'D3'): ['D2'],
-			 ('Z2', 'D3'): []
-			 }
-
 ###########################################
 # Helper functions for theory predictions #
 ###########################################
@@ -3262,6 +3229,7 @@ def set_combined_waveform_amplitudes(flux_lutman):
 	parking for a flux lutman.
 	'''
 	# print(f'Setting common amps in {flux_lutman.name}.')
+	qubit = flux_lutman.name.split('_')[-1]
 	# calculate park detuning
 	park_det = get_frequency_waveform('park_amp', flux_lutman)
 	# Remove default values (need to fix this in pycqed instead)
@@ -3283,29 +3251,31 @@ def set_combined_waveform_amplitudes(flux_lutman):
 	# Find waveform with maximum detuning
 	max_wf = max(Detunings, key=Detunings.get)
 	# Set amplitude of DAC to 0.5 and scale gain accordingly
-	flux_lutman.set(max_wf, 0.3)
+	if qubit in OFFSET_QUBITS:
+		flux_lutman.set(max_wf, -0.3)
+	else:
+		flux_lutman.set(max_wf, 0.3)
 	max_wf_gain = get_Ch_amp_frequency(Detunings[max_wf], flux_lutman,
 									   DAC_param = max_wf)
 	flux_lutman.cfg_awg_channel_amplitude(max_wf_gain)
-	# print(max_wf, 0.5)
 	Detunings.pop(max_wf) # remove waveform from detuning dict
 	# Set remaining waveform amplitudes
 	for wf, det in Detunings.items():
 		if det > 20e6:
-			wf_amp = get_DAC_amp_frequency(det, flux_lutman)
+			wf_amp = get_DAC_amp_frequency(det, flux_lutman,
+						negative_amp=True if qubit in OFFSET_QUBITS else False)
 		else:
 			wf_amp = 0
 		flux_lutman.set(wf, wf_amp)
-		# print(wf, wf_amp)
-	# print('')
 
-def prepare_for_parity_check(stabilizer_qubit, station):
+def prepare_for_parity_check(stabilizer_qubit, station,
+							 Data_qubits: list = None):
 	'''
 	Wrapper function to prepare for timedomain of parity check of
 	a stabilizer.
 	'''
 	# Set configuration
-	TQG_duration_ns = int(TWOQ_GATE_DURATION*1e9) - 20
+	TQG_duration_ns = TWOQ_GATE_DURATION_NS - 20
 	file_cfg = gc.generate_config(in_filename=input_file,
 	                              out_filename=config_fn,
 	                              mw_pulse_duration=20,
@@ -3315,7 +3285,8 @@ def prepare_for_parity_check(stabilizer_qubit, station):
 	device = station.components['device']
 	device.ro_acq_weight_type('optimal')
 	# Get qubits directly involved in parity check
-	Data_qubits = list(get_nearest_neighbors(stabilizer_qubit).keys())
+	if not Data_qubits:
+		Data_qubits = list(get_nearest_neighbors(stabilizer_qubit).keys())
 	PC_qubits = [stabilizer_qubit]+Data_qubits
 	# Get spectator qubits of parity check
 	Spec_qubits = []
@@ -3579,7 +3550,7 @@ def check_flux_wf_duration(flux_lutman):
 		flux_lutman.cfg_max_wf_length(TWOQ_GATE_DURATION)
 		awg = flux_lutman.AWG.get_instr()
 		awg.reset_waveforms_zeros()
-		print(f'Loading waveforms to match {int(TWOQ_GATE_DURATION*1e9):.0f} '+\
+		print(f'Loading waveforms to match {TWOQ_GATE_DURATION_NS:.0f} '+\
 			   'ns gate duration')
 		flux_lutman.load_waveforms_onto_AWG_lookuptable()
 
@@ -3747,15 +3718,16 @@ def save_snapshot_metadata(station, Qubits=None, Qubit_pairs = None,
 		ma2.tqg.TwoQubitGate_frequency_trajectory_analysis(Qubit_pairs=Qubit_pairs)
 	# Plot parity-check benchmarks
 	if parity_check:
-		ma2.gbta.ParityCheckGBT_analysis(Stabilizers=['Z1', 'Z2', 'Z3', 'Z4'])
+		ma2.gbta.ParityCheckGBT_analysis(Stabilizers=['Z1', 'Z2', 'Z3', 'Z4', 'X1', 'X2', 'X3', 'X4'])
 	return True
 
-def DIO_calibration(station):
+def DIO_calibration(station, force: bool = False):
 	'''
 	Checks for DIO errors in all instruments and calibrates
 	them if error is found.
 	'''
 	# Get all intruments
+	no_error = True
 	cc = station.components['cc']
 	UHFQC_1 = station.components['UHFQC_1']
 	UHFQC_2 = station.components['UHFQC_2']
@@ -3792,7 +3764,8 @@ def _prep_awg(awg):
 	UHFQC_1.check_errors()
 	_errors = UHFQC_1._errors
 	# if 'AWGDIOTIMING' in _errors.keys():
-	if _errors != {}:
+	if _errors != {} or force:
+		# no_error = False
 		UHFQC_1._errors = {}
 		print(f'Calibrating DIO on UHFQC_1.')
 		try:
@@ -3809,7 +3782,8 @@ def _prep_awg(awg):
 	UHFQC_2.check_errors()
 	_errors = UHFQC_2._errors
 	# if 'AWGDIOTIMING' in _errors.keys():
-	if _errors != {}:
+	if _errors != {} or force:
+		# no_error = False
 		UHFQC_2._errors = {}
 		print(f'Calibrating DIO on UHFQC_2.')
 		try:
@@ -3826,7 +3800,8 @@ def _prep_awg(awg):
 	UHFQC_3.check_errors()
 	_errors = UHFQC_3._errors
 	# if 'AWGDIOTIMING' in _errors.keys():
-	if _errors != {}:
+	if _errors != {} or force:
+		# no_error = False	# this is commented because some holdoff error cannot be fixed
 		UHFQC_3._errors = {}
 		print(f'Calibrating DIO on UHFQC_3.')
 		try:
@@ -3843,7 +3818,8 @@ def _prep_awg(awg):
 	UHFQC_4.check_errors()
 	_errors = UHFQC_4._errors
 	# if 'AWGDIOTIMING' in _errors.keys():
-	if _errors != {}:
+	if _errors != {} or force:
+		# no_error = False
 		UHFQC_4._errors = {}
 		print(f'Calibrating DIO on UHFQC_4.')
 		try:
@@ -3864,7 +3840,8 @@ def _prep_awg(awg):
 	AWG8_8481.check_errors()
 	_errors = AWG8_8481._errors
 	# if 'AWGDIOTIMING' in _errors.keys():
-	if _errors != {}:
+	if _errors != {} or force:
+		no_error = False
 		AWG8_8481._errors = {}
 		print(f'Calibrating DIO on AWG8_8481.')
 		AWG8_8481.set('dios_0_interface', 0)
@@ -3882,7 +3859,8 @@ def _prep_awg(awg):
 	AWG8_8068.check_errors()
 	_errors = AWG8_8068._errors
 	# if 'AWGDIOTIMING' in _errors.keys():
-	if _errors != {}:
+	if _errors != {} or force:
+		no_error = False
 		AWG8_8068._errors = {}
 		print(f'Calibrating DIO on AWG8_8068.')
 		AWG8_8068.set('dios_0_interface', 0)
@@ -3900,7 +3878,8 @@ def _prep_awg(awg):
 	AWG8_8074.check_errors()
 	_errors = AWG8_8074._errors
 	# if 'AWGDIOTIMING' in _errors.keys():
-	if _errors != {}:
+	if _errors != {} or force:
+		no_error = False
 		AWG8_8074._errors = {}
 		print(f'Calibrating DIO on AWG8_8074.')	
 		AWG8_8074.set('dios_0_interface', 0)
@@ -3918,7 +3897,8 @@ def _prep_awg(awg):
 	AWG8_8076.check_errors()
 	_errors = AWG8_8076._errors
 	# if 'AWGDIOTIMING' in _errors.keys():
-	if _errors != {}:
+	if _errors != {} or force:
+		no_error = False
 		AWG8_8076._errors = {}
 		print(f'Calibrating DIO on AWG8_8076.')	
 		AWG8_8076.set('dios_0_interface', 0)
@@ -3936,7 +3916,8 @@ def _prep_awg(awg):
 	AWG8_8499.check_errors()
 	_errors = AWG8_8499._errors
 	# if 'AWGDIOTIMING' in _errors.keys():
-	if _errors != {}:
+	if _errors != {} or force:
+		no_error = False
 		AWG8_8499._errors = {}
 		print(f'Calibrating DIO on AWG8_8499.')	
 		AWG8_8499.set('dios_0_interface', 0)
@@ -3958,7 +3939,8 @@ def _prep_awg(awg):
 	AWG8_8279.check_errors()
 	_errors = AWG8_8279._errors
 	# if 'AWGDIOTIMING' in _errors.keys():
-	if _errors != {}:
+	if _errors != {} or force:
+		no_error = False
 		AWG8_8279._errors = {}
 		print(f'Calibrating DIO on AWG8_8279.')
 		AWG8_8279.set('dios_0_interface', 0)
@@ -3972,24 +3954,25 @@ def _prep_awg(awg):
 		AWG8_8279._set_dio_calibration_delay(6)
 		AWG8_8279_channels = [0, 1, 2, 3, 4, 5, 6, 7]
 		for this_ch in AWG8_8279_channels:
-		    AWG8_8279.setd('sigouts/%d/precompensation/enable'%(int(this_ch)),True)
-		    AWG8_8279.setd('sigouts/%d/precompensation/exponentials/0/enable'%(int(this_ch)),True)
-		    AWG8_8279.setd('sigouts/%d/precompensation/exponentials/1/enable'%(int(this_ch)),True)
-		    AWG8_8279.setd('sigouts/%d/precompensation/exponentials/2/enable'%(int(this_ch)),True)
-		    AWG8_8279.setd('sigouts/%d/precompensation/exponentials/3/enable'%(int(this_ch)),True)
-		    AWG8_8279.setd('sigouts/%d/precompensation/exponentials/4/enable'%(int(this_ch)),True)
-		    AWG8_8279.setd('sigouts/%d/precompensation/exponentials/5/enable'%(int(this_ch)),True)
-		    AWG8_8279.setd('sigouts/%d/precompensation/exponentials/6/enable'%(int(this_ch)),True)
-		    AWG8_8279.setd('sigouts/%d/precompensation/exponentials/7/enable'%(int(this_ch)),True)
-		    AWG8_8279.setd('sigouts/%d/precompensation/fir/enable'%(int(this_ch)),True)
-		    AWG8_8279.set('sigouts_{}_delay'.format(int(this_ch)), 0e-9+4*10/3*1e-9-2*3.33e-9)
+			AWG8_8279.setd('sigouts/%d/precompensation/enable'%(int(this_ch)),True)
+			AWG8_8279.setd('sigouts/%d/precompensation/exponentials/0/enable'%(int(this_ch)),True)
+			AWG8_8279.setd('sigouts/%d/precompensation/exponentials/1/enable'%(int(this_ch)),True)
+			AWG8_8279.setd('sigouts/%d/precompensation/exponentials/2/enable'%(int(this_ch)),True)
+			AWG8_8279.setd('sigouts/%d/precompensation/exponentials/3/enable'%(int(this_ch)),True)
+			AWG8_8279.setd('sigouts/%d/precompensation/exponentials/4/enable'%(int(this_ch)),True)
+			AWG8_8279.setd('sigouts/%d/precompensation/exponentials/5/enable'%(int(this_ch)),True)
+			AWG8_8279.setd('sigouts/%d/precompensation/exponentials/6/enable'%(int(this_ch)),True)
+			AWG8_8279.setd('sigouts/%d/precompensation/exponentials/7/enable'%(int(this_ch)),True)
+			AWG8_8279.setd('sigouts/%d/precompensation/fir/enable'%(int(this_ch)),True)
+			AWG8_8279.set('sigouts_{}_delay'.format(int(this_ch)), 0e-9+4*10/3*1e-9-2*3.33e-9)
 		AWG8_8279.clear_errors()
 		_prep_awg('AWG8_8279')
 	# AWG8_8320
 	AWG8_8320.check_errors()
 	_errors = AWG8_8320._errors
 	# if 'AWGDIOTIMING' in _errors.keys():
-	if _errors != {}:
+	if _errors != {} or force:
+		no_error = False
 		AWG8_8320._errors = {}
 		print(f'Calibrating DIO on AWG8_8320.')
 		AWG8_8320.set('dios_0_interface', 0)
@@ -4020,7 +4003,8 @@ def _prep_awg(awg):
 	AWG8_8071.check_errors()
 	_errors = AWG8_8071._errors
 	# if 'AWGDIOTIMING' in _errors.keys():
-	if _errors != {}:
+	if _errors != {} or force:
+		no_error = False
 		AWG8_8071._errors = {}
 		print(f'Calibrating DIO on AWG8_8071.')
 		AWG8_8071.set('dios_0_interface', 0)
@@ -4053,7 +4037,7 @@ def _prep_awg(awg):
 	device.tim_mw_latency_3(0) # 8068
 	device.tim_mw_latency_4(-10e-9) # 8481
 	device.prepare_timing()
-	return True
+	return no_error
 
 ###############################################################################
 # LRU calibration graph
diff --git a/pycqed/instrument_drivers/meta_instrument/device_object_CCL.py b/pycqed/instrument_drivers/meta_instrument/device_object_CCL.py
index 94f91f348..bc6445c99 100644
--- a/pycqed/instrument_drivers/meta_instrument/device_object_CCL.py
+++ b/pycqed/instrument_drivers/meta_instrument/device_object_CCL.py
@@ -2160,6 +2160,7 @@ def measure_chevron(
         target_qubit_sequence: str = "ramsey",
         waveform_name="square",
         recover_q_spec: bool = False,
+        second_excited_state: bool = False,
         disable_metadata: bool = False,
         ):
         """
@@ -2218,6 +2219,9 @@ def measure_chevron(
             recover_q_spec (bool):
                 applies the first gate of qspec at the end as well if `True`
 
+            second_excited_state (bool):
+                Applies f12 transition pulse before flux pulse.
+
         Circuit:
             q0    -x180-flux-x180-RO-
             qspec --x90-----(x90)-RO- (target_qubit_sequence='ramsey')
@@ -2293,6 +2297,7 @@ def measure_chevron(
             target_qubit_sequence=target_qubit_sequence,
             cc=self.instr_CC.get_instr().name,
             recover_q_spec=recover_q_spec,
+            second_excited_state=second_excited_state,
         )
         self.instr_CC.get_instr().eqasm_program(p.filename)
         self.instr_CC.get_instr().start()
@@ -2309,7 +2314,8 @@ def measure_chevron(
         MC.set_sweep_function_2D(sw)
         MC.set_detector_function(d)
 
-        label = "Chevron {} {} {}".format(q0, q_spec, target_qubit_sequence)
+        prepared_state_label: str = "_prep_state_f" if second_excited_state else "_prep_state_e"
+        label = f"Chevron{prepared_state_label} {q0} {q_spec} {target_qubit_sequence}"
 
         if not adaptive_sampling:
             MC.set_sweep_points(amps)
@@ -3455,7 +3461,7 @@ def measure_multi_qubit_simultaneous_randomized_benchmarking(
         """
         Performs simultaneous single qubit RB on multiple qubits.
         The data of this experiment should be compared to the results of single
-        qubit RB to reveal differences due to crosstalk and residual coupling
+        qubit RB to reveal differences due to MW crosstalk and residual coupling
 
         Args:
             qubits (list):
@@ -3485,8 +3491,6 @@ def measure_multi_qubit_simultaneous_randomized_benchmarking(
         self.ro_acq_weight_type(ro_acq_weight_type)
         self.ro_acq_digitized(False)
 
-        if prepare_for_timedomain:
-            self.prepare_for_timedomain(qubits=qubits, bypass_flux=True)
         if MC is None:
             MC = self.instr_MC.get_instr()
         MC.soft_avg(1)
@@ -3494,13 +3498,18 @@ def measure_multi_qubit_simultaneous_randomized_benchmarking(
         # The detector needs to be defined before setting back parameters
         d = self.get_int_logging_detector(qubits=qubits)
         # set back the settings
-        self.ro_acq_weight_type(old_weight_type)
-        self.ro_acq_digitized(old_digitized)
+        # self.ro_acq_weight_type(old_weight_type)
+        # self.ro_acq_digitized(old_digitized)
 
         for q in qubits:
             q_instr = self.find_instrument(q)
             mw_lutman = q_instr.instr_LutMan_MW.get_instr()
-            mw_lutman.load_ef_rabi_pulses_to_AWG_lookuptable()
+            # mw_lutman.load_ef_rabi_pulses_to_AWG_lookuptable()
+            mw_lutman.set_default_lutmap()
+
+
+        if prepare_for_timedomain:
+            self.prepare_for_timedomain(qubits=qubits, bypass_flux=True)
 
         MC.soft_avg(1)
 
@@ -3606,6 +3615,187 @@ def send_rb_tasks(pool_):
 
         return Analysis
 
+    def measure_two_qubit_simultaneous_randomized_benchmarking(
+        self,
+        qubits,
+        MC= None,
+        nr_cliffords=2 ** np.arange(11),
+        nr_seeds=100,
+        interleaving_cliffords=[None],
+        label="TwoQubit_sim_RB_{}seeds_recompile={}_{}_{}",
+        recompile: bool = "as needed",
+        cal_points: bool = True,
+        ro_acq_weight_type: str = "optimal IQ",
+        compile_only: bool = False,
+        pool=None,  # a multiprocessing.Pool()
+        rb_tasks=None  # used after called with `compile_only=True`
+    ):
+        """
+        Performs simultaneous single qubit RB on two qubits.
+        The data of this experiment should be compared to the results of single
+        qubit RB to reveal differences due to crosstalk and residual coupling
+
+        Args:
+            qubits (list):
+                pair of the qubit names on which to perform RB
+
+            nr_cliffords (array):
+                lengths of the clifford sequences to perform
+
+            nr_seeds (int):
+                number of different clifford sequences of each length
+
+            interleaving_cliffords (list):
+                list of integers (or None) which specifies which cliffords
+                to interleave the sequence with (for interleaved RB)
+                For indices of Clifford group elements go to
+                two_qubit_clifford_group.py
+
+            label (str):
+                string for formatting the measurement name
+
+            recompile (bool, str {'as needed'}):
+                indicate whether to regenerate the sequences of clifford gates.
+                By default it checks whether the needed sequences were already
+                generated since the most recent change of OpenQL file
+                specified in self.cfg_openql_platform_fn
+
+            cal_points (bool):
+                should calibration point (qubits in 0, 1 and 2 states)
+                be included in the measurement
+        """
+
+        # Settings that have to be preserved, change is required for
+        # 2-state readout and postprocessing
+        old_weight_type = self.ro_acq_weight_type()
+        old_digitized = self.ro_acq_digitized()
+        self.ro_acq_weight_type(ro_acq_weight_type)
+        self.ro_acq_digitized(False)
+
+        for q in qubits:
+            q_instr = self.find_instrument(q)
+            mw_lutman = q_instr.instr_LutMan_MW.get_instr()
+            mw_lutman.set_default_lutmap()
+
+        self.prepare_for_timedomain(qubits=qubits)
+        if MC is None:
+            MC = self.instr_MC.get_instr()
+        MC.soft_avg(1)
+
+        # The detector needs to be defined before setting back parameters
+        d = self.get_int_logging_detector(qubits=qubits)
+        # set back the settings
+        self.ro_acq_weight_type(old_weight_type)
+        self.ro_acq_digitized(old_digitized)
+
+        def send_rb_tasks(pool_):
+            tasks_inputs = []
+            for i in range(nr_seeds):
+                task_dict = dict(
+                    qubits=[self.find_instrument(q).cfg_qubit_nr() for q in qubits],
+                    nr_cliffords=nr_cliffords,
+                    nr_seeds=1,
+                    platf_cfg=self.cfg_openql_platform_fn(),
+                    program_name="TwoQ_Sim_RB_int_cl{}_s{}_ncl{}_{}_{}_double".format(
+                        i,
+                        list(map(int, nr_cliffords)),
+                        interleaving_cliffords,
+                        qubits[0],
+                        qubits[1],
+                    ),
+                    interleaving_cliffords=interleaving_cliffords,
+                    simultaneous_single_qubit_RB=True,
+                    cal_points=cal_points,
+                    net_cliffords=[0, 3],  # measures with and without inverting
+                    f_state_cal_pts=True,
+                    recompile=recompile,
+                )
+                tasks_inputs.append(task_dict)
+            # pool.starmap_async can be used for positional arguments
+            # but we are using a wrapper
+            rb_tasks = pool_.map_async(cl_oql.parallel_friendly_rb, tasks_inputs)
+
+            return rb_tasks
+
+        if compile_only:
+            assert pool is not None
+            rb_tasks = send_rb_tasks(pool)
+            return rb_tasks
+
+        if rb_tasks is None:
+            # Using `with ...:` makes sure the other processes will be terminated
+            # avoid starting too mane processes,
+            # nr_processes = None will start as many as the PC can handle
+            nr_processes = None if recompile else 1
+            with multiprocessing.Pool(
+                nr_processes,
+                maxtasksperchild=cl_oql.maxtasksperchild  # avoid RAM issues
+            ) as pool:
+                rb_tasks = send_rb_tasks(pool)
+                cl_oql.wait_for_rb_tasks(rb_tasks)
+
+        programs_filenames = rb_tasks.get()
+
+        # to include calibration points
+        if cal_points:
+            sweep_points = np.append(
+                np.repeat(nr_cliffords, 2),
+                [nr_cliffords[-1] + 0.5] * 2
+                + [nr_cliffords[-1] + 1.5] * 2
+                + [nr_cliffords[-1] + 2.5] * 3,
+            )
+        else:
+            sweep_points = np.repeat(nr_cliffords, 2)
+
+        counter_param = ManualParameter("name_ctr", initial_value=0)
+        prepare_function_kwargs = {
+            "counter_param": counter_param,
+            "programs_filenames": programs_filenames,
+            "CC": self.instr_CC.get_instr(),
+        }
+
+        # Using the first detector of the multi-detector as this is
+        # in charge of controlling the CC (see self.get_int_logging_detector)
+        d.set_prepare_function(
+            oqh.load_range_of_oql_programs_from_filenames,
+            prepare_function_kwargs, detectors="first"
+        )
+        # d.nr_averages = 128
+
+        reps_per_seed = 4094 // len(sweep_points)
+        d.set_child_attr("nr_shots", reps_per_seed * len(sweep_points))
+
+        s = swf.None_Sweep(parameter_name="Number of Cliffords", unit="#")
+
+        MC.set_sweep_function(s)
+        MC.set_sweep_points(np.tile(sweep_points, reps_per_seed * nr_seeds))
+
+        MC.set_detector_function(d)
+        label = label.format(nr_seeds, recompile, qubits[0], qubits[1])
+        MC.run(label, exp_metadata={"bins": sweep_points})
+
+        # N.B. if interleaving cliffords are used, this won't work
+        # [2020-07-11 Victor] not sure if NB still holds
+
+        cal_2Q = ["00", "01", "10", "11", "02", "20", "22"]
+
+        rates_I_quad_ch_idx = 0
+        cal_1Q = [state[rates_I_quad_ch_idx // 2] for state in cal_2Q]
+        a_q0 = ma2.RandomizedBenchmarking_SingleQubit_Analysis(
+            label=label,
+            rates_I_quad_ch_idx=rates_I_quad_ch_idx,
+            cal_pnts_in_dset=cal_1Q
+        )
+        rates_I_quad_ch_idx = 2
+        cal_1Q = [state[rates_I_quad_ch_idx // 2] for state in cal_2Q]
+        a_q1 = ma2.RandomizedBenchmarking_SingleQubit_Analysis(
+            label=label,
+            rates_I_quad_ch_idx=rates_I_quad_ch_idx,
+            cal_pnts_in_dset=cal_1Q
+        )
+
+        return a_q0, a_q1
+
     def measure_gate_process_tomography(
         self,
         meas_qubit: str,
@@ -4875,6 +5065,7 @@ def measure_parity_check_ramsey(
                     # Calculate new virtual phase
                     phi0 = mw_lm.get(param)
                     phi_new = list(a.qoi['Phase_model'][Q.name].values())[0]
+                    phi_new = phi_new / pc_repetitions  # Divide by number of CZ repetitions
                     phi = np.mod(phi0+phi_new, 360)
                     mw_lm.set(param, phi)
                     print(f'{Q.name}.{param} changed to {phi} deg.')
@@ -4928,10 +5119,10 @@ def calibrate_parity_check_phase(
         B0 = fl_lm.get(fl_par)
         if B_amps is None:
             B_amps = np.linspace(-.1, .1, 3)+B0
-            if np.min(B_amps) < 0:
-                B_amps -= np.min(B_amps)
-            if np.max(B_amps) > 1:
-                B_amps -= np.max(B_amps)-1
+        if np.min(B_amps) < 0:
+            B_amps -= np.min(B_amps)
+        if np.max(B_amps) > 1:
+            B_amps -= np.max(B_amps)-1
 
         # Prepare for timedomain
         if prepare_for_timedomain:
@@ -5843,6 +6034,7 @@ def measure_defect_rate(
             initial_state_qubits: list = None,
             measurement_time_ns: int = 500,
             analyze: bool = True,
+            Pij_matrix: bool = True,
             ):
         # assert self.ro_acq_weight_type() == 'optimal IQ'
         assert self.ro_acq_digitized() == False
@@ -5868,12 +6060,14 @@ def measure_defect_rate(
         ######################################################
         # Prepare for timedomain
         ######################################################
+        # prepare mw lutmans
+        # for q in [ancilla_qubit]+data_qubits:
+        for q in Data_qubits + X_ancillas + Z_ancillas:
+            mw_lm = self.find_instrument(f'MW_lutman_{q}')
+            mw_lm.set_default_lutmap()
+            mw_lm.load_waveforms_onto_AWG_lookuptable()
+
         if prepare_for_timedomain:
-            # prepare mw lutmans
-            # for q in [ancilla_qubit]+data_qubits:
-            for q in Data_qubits+X_ancillas+Z_ancillas:
-                mw_lm = self.find_instrument(f'MW_lutman_{q}')
-                mw_lm.set_default_lutmap()
             # Redundancy just to be sure we are uploading every parameter
             # for q_name in data_qubits+[ancilla_qubit]:
             for q_name in Data_qubits+X_ancillas+Z_ancillas:
@@ -5989,6 +6183,8 @@ def measure_defect_rate(
         else:
             _title = f'Repeated_stab_meas_{"_".join([str(r) for r in Rounds])}rounds'+\
                      f'_{ancilla_qubit}_{data_qubits}_data_qubit_measurement'
+        if len(_title) > 96:
+            _title = _title[:96]    # this is to avoid failure in creating hdf5 file.
         try:
             MC.run(_title)
             a = None
@@ -6023,6 +6219,7 @@ def measure_repetition_code_defect_rate(
             analyze: bool = True,
             Pij_matrix: bool = True,
             disable_metadata: bool = False,
+            initial_state: list = None,
             ):
         # assert self.ro_acq_weight_type() == 'optimal IQ'
         assert self.ro_acq_digitized() == False
@@ -6135,7 +6332,7 @@ def internal_prepare_for_timedomain():
         internal_prepare_for_timedomain()
         
         # Generate compiler sequence
-        p = mqo.repetition_code_sequence(
+        p = mqo.repetition_code_sequence_old(
             involved_ancilla_indices=involved_ancilla_indices,
             involved_data_indices=involved_data_indices,
             all_ancilla_indices=all_ancilla_indices,
@@ -6144,13 +6341,15 @@ def internal_prepare_for_timedomain():
             platf_cfg=self.cfg_openql_platform_fn(),
             stabilizer_type=stabilizer_type,
             measurement_time_ns=measurement_time_ns,
+            initial_state=initial_state
         )
         # Set up nr_shots on detector
         d = self.int_log_det
-        uhfqc_max_avg = 2**17  # 2**19
+        uhfqc_max_avg = 2**19  # 2**19
         number_of_kernels: int = 1  # Performing only a single experiment
         for det in d.detectors:
-            readouts_per_round = np.sum(np.array(rounds)+heralded_init) * number_of_kernels + 3*(1+heralded_init)
+            nr_acquisitions = [1 if round == 0 else round for round in rounds]
+            readouts_per_round = np.sum(np.array(nr_acquisitions)+heralded_init) * number_of_kernels + 3*(1+heralded_init)
             det.nr_shots = int(uhfqc_max_avg/readouts_per_round)*readouts_per_round
 
         s = swf.OpenQL_Sweep(openql_program=p, CCL=self.instr_CC.get_instr())
diff --git a/pycqed/measurement/flux_crosstalk_ac/measurement_function.py b/pycqed/measurement/flux_crosstalk_ac/measurement_function.py
index 1f09a28ec..ce8a6a192 100644
--- a/pycqed/measurement/flux_crosstalk_ac/measurement_function.py
+++ b/pycqed/measurement/flux_crosstalk_ac/measurement_function.py
@@ -1,7 +1,12 @@
 # -------------------------------------------
 # Module describing the measurement functionality used for AC flux-crosstalk experiment
 # -------------------------------------------
+from abc import ABC, abstractmethod
+import os
+import functools
 import numpy as np
+import warnings
+from typing import TypeVar, Type, Dict, Optional
 from pycqed.instrument_drivers.meta_instrument.device_object_CCL import DeviceCCL as Device
 from pycqed.instrument_drivers.meta_instrument.qubit_objects.CCL_Transmon import CCLight_Transmon as Transmon
 from pycqed.measurement.measurement_control import MeasurementControl as MeasurementControl
@@ -16,8 +21,185 @@
     schedule_flux_crosstalk,
     schedule_ramsey,
 )
+from pycqed.analysis.analysis_toolbox import get_datafilepath_from_timestamp
+from pycqed.analysis_v2.base_analysis import BaseDataAnalysis
+import pycqed.measurement.hdf5_data as hd5
+import matplotlib.pyplot as plt
+from pycqed.qce_utils.custom_exceptions import InterfaceMethodException
 
 
+class IBaseDataAnalysis(ABC):
+    
+    # region Interface Methods
+    @abstractmethod
+    def process_data(self):
+        """
+        process_data: overloaded in child classes,
+        takes care of mundane tasks such as binning filtering etc
+        """
+        raise InterfaceMethodException
+
+    @abstractmethod
+    def prepare_plots(self):
+        """
+        Defines a default plot by setting up the plotting dictionaries to
+        specify what is to be plotted
+        """
+        raise InterfaceMethodException
+
+    @abstractmethod
+    def analyze_fit_results(self):
+        """
+        Do analysis on the results of the fits to extract quantities of
+        interest.
+        """
+        raise InterfaceMethodException
+    # endregion
+
+
+def plot_example(x_array: np.ndarray, y_array: np.ndarray, ax: Optional[plt.Axes] = None, **kwargs):
+    # Get keyword argument defaults
+    x_label: str = kwargs.pop('x_label', 'Default Label [a.u.]')
+    y_label: str = kwargs.pop('y_label', 'Default Label [a.u.]')
+    # Plot figure
+    print(kwargs)
+    if ax is None:
+        fig, ax = plt.subplots()
+    ax.plot(
+        x_array,
+        y_array,
+        '.-',
+    )
+    ax.grid(True, alpha=0.5, linestyle='dashed')
+    ax.set_axisbelow(True)
+    ax.set_xlabel(x_label)
+    ax.set_ylabel(y_label)
+
+TBaseDataAnalysis = TypeVar('TBaseDataAnalysis', bound=BaseDataAnalysis)
+
+class FluxCrosstalkAnalysis(BaseDataAnalysis, IBaseDataAnalysis):
+    
+    # region Class Constructor
+    def __init__(self, t_start: str = None, t_stop: str = None, label: str = '', data_file_path: str = None, close_figs: bool = True, options_dict: dict = None, extract_only: bool = False, do_fitting: bool = False, save_qois: bool = True):
+        super().__init__(t_start, t_stop, label, data_file_path, close_figs, options_dict, extract_only, do_fitting, save_qois)
+        self.params_dict: Dict = {}
+        self.numeric_params: Dict = {}
+    # endregion
+    
+    # region Interface Methods
+    def extract_data(self):
+        """
+        This is a new style (sept 2019) data extraction.
+        This could at some point move to a higher level class.
+        """
+        self.get_timestamps()
+        self.timestamp = self.timestamps[0]
+        data_filepath = get_datafilepath_from_timestamp(self.timestamp)
+        self._raw_data_key: str = 'data'
+        param_spec = {
+            self._raw_data_key: ('Experimental Data/Data', 'dset'),
+            'value_names': ('Experimental Data', 'attr:value_names'),
+        }
+        self.raw_data_dict = hd5.extract_pars_from_datafile(data_filepath, param_spec)
+        self.raw_data_dict['timestamps'] = self.timestamps
+        self.raw_data_dict['folder'] = os.path.split(data_filepath)[0]
+        
+    def process_data(self):
+        """
+        process_data: overloaded in child classes,
+        takes care of mundane tasks such as binning filtering etc
+        """
+        self.data_array: np.ndarray = self.raw_data_dict[self._raw_data_key]
+
+    def prepare_plots(self):
+        """
+        Defines a default plot by setting up the plotting dictionaries to
+        specify what is to be plotted
+        """
+        self.plot_dicts['Example_plot'] = {
+            'plotfn': plot_example,
+            'x_array': self.data_array[:, 0],
+            'y_array': self.data_array[:, 1],
+            'x_label': 'Times [s]',
+        }
+
+    def analyze_fit_results(self):
+        """
+        Do analysis on the results of the fits to extract quantities of
+        interest.
+        """
+        # raise NotImplemented
+        pass
+    # endregion
+
+
+def decorator_run_analysis(analysis_class: Type[TBaseDataAnalysis], filter_label: str):
+    
+    def decorator(func):
+        """
+        Decorator that constructs analysis and triggers execution.
+        """
+        @functools.wraps(func)
+        def wrapper(*args, **kwargs):
+            result = func(*args, **kwargs)
+            
+            # Construct analysis class
+            instance: TBaseDataAnalysis = analysis_class(label=filter_label)
+            instance.run_analysis()
+            
+            return result
+        return wrapper
+    
+    return decorator
+
+
+def decorator_pyplot_dataset(dataset_key: str, array_dimensions: int, show_figure: bool, *plot_args, **plot_kwargs):
+    
+    def decorator(func):
+        """
+        Decorator to measure and record the execution time of the decorated function or method.
+        """
+        @functools.wraps(func)
+        def wrapper(*args, **kwargs):
+            result = func(*args, **kwargs)
+            # Assumes dictionary structure
+            dataset_key_exists: bool = dataset_key in result and isinstance(result[dataset_key], np.ndarray)
+            # Guard clause, if dataset key does not exists
+            if not dataset_key_exists:
+                warnings.warn(f"Key: {dataset_key} is not present in function output.")
+                return result
+            
+            data_array: np.ndarray = result[dataset_key]
+            # Get keyword argument defaults
+            x_label: str = plot_kwargs.pop('x_label', 'Default Label [a.u.]')
+            y_label: str = plot_kwargs.pop('y_label', 'Default Label [a.u.]')
+            # Plot figure
+            fig, ax = plt.subplots(*plot_args, **plot_kwargs)
+            ax.plot(
+                data_array[:, 0],
+                data_array[:, 1],
+                '.-',
+            )
+            ax.grid(True, alpha=0.5, linestyle='dashed')
+            ax.set_axisbelow(True)
+            ax.set_xlabel(x_label)
+            ax.set_ylabel(y_label)
+            
+            if show_figure:
+                plt.show()
+            else:
+                plt.close(fig)
+            
+            return result
+        return wrapper
+    
+    return decorator
+
+
+@decorator_run_analysis(
+    analysis_class=FluxCrosstalkAnalysis,
+    filter_label='FluxCrosstalk',
+)
 def measure_ac_flux_crosstalk(device: Device, qubit_echo_id: str, prepare_for_timedomain: bool = False, disable_metadata: bool = False):
     """
     Performs an experiment
@@ -25,41 +207,43 @@ def measure_ac_flux_crosstalk(device: Device, qubit_echo_id: str, prepare_for_ti
     # Data allocation
     qubit_echo: Transmon = device.find_instrument(qubit_echo_id)
     flux_lutman: FluxLutMan = qubit_echo.instr_LutMan_Flux.get_instr()
-    times: np.ndarray = np.arange(0, 20e-6, 1e-6)  # Echo time
-    # sweep_function: FluxSweepFunctionObject = FluxSweepFunctionObject(
-    #     flux_lutman,
-    #     flux_lutman.sq_length,  # Square length (qcodes) parameter
-    #     waveform_name="square",  # Meaning we are going to sweep the square-pulse parameters only
-    # )
+    max_duration_ns: int = 10  # [ns]
+    flux_pulse_duration: np.ndarray = np.arange(0e-9, max_duration_ns * 1e-9, 1/2.4e9)
+    # flux_pulse_amplitude: np.ndarray = np.asarray([flux_lutman.sq_amp()])
+    flux_cw = "sf_square"
     meas_control: MeasurementControl = device.instr_MC.get_instr()
     meas_control_nested: MeasurementControl = device.instr_nested_MC.get_instr()
     central_control: CentralControl = device.instr_CC.get_instr()
 
     # Prepare for time-domain if requested
     if prepare_for_timedomain:
-        device.prepare_for_timedomain(qubits=[qubit_echo])
-
-    # schedule: OqlProgram = schedule_flux_crosstalk(
-    #     qubit_echo_index=qubit_echo.cfg_qubit_nr(),
-    #     flux_pulse_cw='fl_cw_06',  # Square pulse?
-    # )
+        device.prepare_for_timedomain(qubits=[qubit_echo_id])
 
-    schedule: OqlProgram = schedule_ramsey(
+    schedule: OqlProgram = schedule_flux_crosstalk(
         qubit_echo_index=qubit_echo.cfg_qubit_nr(),
-        half_echo_delays=times,
-        flux_pulse_cw='fl_cw_06',  # Square pulse?
+        flux_pulse_cw=flux_cw,  # Square pulse?
+        platf_cfg=device.cfg_openql_platform_fn(),
+        half_echo_delay_ns=max_duration_ns,  # [ns]
     )
-
-    sweep_function: OpenQLSweep = OpenQLSweep(
-        openql_program=schedule,
-        CCL=central_control,
-        parameter_name='Time',
-        unit='s',
+    sweep_function: FluxSweepFunctionObject = FluxSweepFunctionObject(
+        flux_lutman,
+        flux_lutman.sq_length,  # Square length (qcodes) parameter
+        amp_for_generation=0.0,
+        waveform_name="square",  # Meaning we are going to sweep the square-pulse parameters only
+    )
+    # flux_pulse_duration = np.concatenate([flux_pulse_duration, np.zeros(shape=4)])  # Include calibration points
+    central_control.eqasm_program(schedule.filename)
+    central_control.start()
+    
+    detector = device.get_int_avg_det(
+        qubits=[qubit_echo_id],
+        single_int_avg=True,
+        always_prepare=True,
     )
 
     meas_control.set_sweep_function(sweep_function)
-    meas_control.set_sweep_points(times)
-    meas_control.set_detector_function(qubit_echo.int_avg_det)
+    meas_control.set_sweep_points(flux_pulse_duration)
+    meas_control.set_detector_function(detector)
     label = f'FluxCrosstalk_{qubit_echo.name}'
     result = meas_control.run(label, disable_snapshot_metadata=disable_metadata)
 
diff --git a/pycqed/measurement/flux_crosstalk_ac/schedule.py b/pycqed/measurement/flux_crosstalk_ac/schedule.py
index c584548e5..1bd2454ff 100644
--- a/pycqed/measurement/flux_crosstalk_ac/schedule.py
+++ b/pycqed/measurement/flux_crosstalk_ac/schedule.py
@@ -5,47 +5,11 @@
 from pycqed.measurement.openql_experiments.openql_helpers import OqlProgram
 
 
-def schedule_ramsey(
-        qubit_echo_index: int,
-        half_echo_delays: np.ndarray,
-        flux_pulse_cw: str,  # = 'fl_cw_06',
-        platf_cfg: str = '',
-    ) -> OqlProgram:
-    """
-    """
-    p = OqlProgram("Ramsey", platf_cfg)
-    flux_pulse_duration_ns: int = 60
-
-    for i, time in enumerate(half_echo_delays[:-4]):
-        half_echo_delay_ns: int = int(round(time/1e-9))
-
-        k = p.create_kernel("Ramsey_{}".format(i))
-        k.prepz(qubit_echo_index)
-        # Echo part 1
-        k.gate('rx90', [qubit_echo_index])
-        k.gate('wait', [], half_echo_delay_ns)
-        k.barrier([])
-        # Echo part 2
-        for i in range(int(half_echo_delay_ns / flux_pulse_duration_ns)):
-            k.gate(flux_pulse_cw, [qubit_echo_index])
-            k.gate("i", [qubit_echo_index])
-        k.gate("wait", [qubit_echo_index], 0)
-        k.gate('rx90', [qubit_echo_index])
-        k.measure(qubit_echo_index)
-        p.add_kernel(k)
-
-    # adding the calibration points
-    p.add_single_qubit_cal_points(qubit_idx=qubit_echo_index)
-
-    p.compile()
-    return p
-
-
 def schedule_flux_crosstalk(
         qubit_echo_index: int,
         flux_pulse_cw: str,  # = 'fl_cw_06',
         platf_cfg: str = '',
-        half_echo_delay: int = 0,  # Delay in ns with respect to end of qubit gate
+        half_echo_delay_ns: int = 0,  # Delay in ns with respect to end of qubit gate
     ) -> OqlProgram:
     """
     """
@@ -58,20 +22,21 @@ def schedule_flux_crosstalk(
     k.barrier([])
     # Echo part 1
     k.gate('rx90', [qubit_echo_index])
-    k.gate('wait', [], half_echo_delay)
+    k.gate('wait', [], half_echo_delay_ns)
     k.barrier([])
     # Echo part 2
+    k.gate('ry180', [qubit_echo_index])
     k.gate(flux_pulse_cw, [qubit_echo_index])
     k.barrier([])  # alignment workaround
-    k.gate('wait', [], half_echo_delay)
-    k.gate('rx90', [qubit_echo_index])
+    k.gate('wait', [], half_echo_delay_ns)
+    k.gate('ry90', [qubit_echo_index])
     k.barrier([])
     # Measure
     k.measure(qubit_echo_index)
     p.add_kernel(k)
 
     # Calibration kernel
-    p.add_single_qubit_cal_points(qubit_idx=qubit_echo_index)
+    # p.add_single_qubit_cal_points(qubit_idx=qubit_echo_index)
 
     p.compile()
     return p
diff --git a/pycqed/measurement/openql_experiments/config_cc_s17_direct_iq.json.in b/pycqed/measurement/openql_experiments/config_cc_s17_direct_iq.json.in
index 9177ef5f6..86feb0bd0 100644
--- a/pycqed/measurement/openql_experiments/config_cc_s17_direct_iq.json.in
+++ b/pycqed/measurement/openql_experiments/config_cc_s17_direct_iq.json.in
@@ -892,88 +892,582 @@
                                      "barrier q9, q5, q15, q3, q8, q2, q4, q12, q7, q0, q1, q10, q16, q13, q14, q11, q6"],
 
 
-        // fluxing steps for parity checks in a distance_9 repetition code
+        // // fluxing steps for parity checks in a distance_9 repetition code
+        // "repetition_code_1 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_nw q11", "sf_cz_se q6",  // X1-D1 2Q gate
+        //                         "sf_cz_nw q7",  "sf_cz_se q15", // Z1-D4 2Q gate
+        //                         "sf_cz_nw q3",  "sf_cz_se q1",  // X3-D7 2Q gate
+        //                         "sf_cz_nw q14", "sf_cz_se q16", // Z2-D6 2Q gate
+        //                         "sf_park q2",  // D2 park
+        //                         "sf_park q8",  // X2 park
+        //                         "sf_park q10", // Z4 park
+        //                         "sf_park q5",  // D8 park
+        //                         "sf_park q12", // Z3 park
+        //                         "barrier q11, q6, q7, q15, q3, q1, q14, q16, q2, q8, q10, q5, q12",
+        //                         "update_ph_nw q11",
+        //                         "update_ph_se q6",
+        //                         "update_ph_nw q7",
+        //                         "update_ph_se q15",
+        //                         "update_ph_nw q3",
+        //                         "update_ph_se q1",
+        //                         "update_ph_nw q14",
+        //                         "update_ph_se q16",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_2 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_ne q11", "sf_cz_sw q2",  // X1-D2 2Q gate
+        //                         "sf_cz_ne q7",  "sf_cz_sw q13", // Z1-D5 2Q gate
+        //                         "sf_cz_ne q3",  "sf_cz_sw q5",  // X3-D8 2Q gate
+        //                         "sf_cz_ne q0",  "sf_cz_sw q14", // Z2-D3 2Q gate
+        //                         "sf_park q6",  // D1 park
+        //                         "sf_park q8",  // X2 park
+        //                         "sf_park q10", // Z4 park
+        //                         "sf_park q1",  // D7 park
+        //                         "barrier q11, q2, q7, q13, q3, q5, q0, q14, q6, q8, q10, q1",
+        //                         "update_ph_ne q11",
+        //                         "update_ph_sw q2",
+        //                         "update_ph_ne q7",
+        //                         "update_ph_sw q13",
+        //                         "update_ph_ne q3",
+        //                         "update_ph_sw q5",
+        //                         "update_ph_ne q0",
+        //                         "update_ph_sw q14",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_3 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_sw q8",  "sf_cz_ne q2",  // X2-D2 2Q gate
+        //                         "sf_cz_sw q10", "sf_cz_ne q13", // Z4-D5 2Q gate
+        //                         "sf_cz_sw q9",  "sf_cz_ne q5",  // X4-D8 2Q gate
+        //                         "sf_cz_sw q1",  "sf_cz_ne q12", // D7-Z3 2Q gate
+        //                         "sf_park q0",  // D3 park
+        //                         "sf_park q7",  // Z1 park
+        //                         "sf_park q3",  // X3 park
+        //                         "sf_park q4",  // D9 park
+        //                         "barrier q8, q2, q10, q13, q9, q5, q1, q12, q0, q7, q3, q4",
+        //                         "update_ph_sw q8",
+        //                         "update_ph_ne q2",
+        //                         "update_ph_sw q10",
+        //                         "update_ph_ne q13",
+        //                         "update_ph_sw q9",
+        //                         "update_ph_ne q5",
+        //                         "update_ph_sw q1",
+        //                         "update_ph_ne q12",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_4 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_se q8",  "sf_cz_nw q0",  // X2-D3 2Q gate
+        //                         "sf_cz_se q10", "sf_cz_nw q16", // Z4-D6 2Q gate
+        //                         "sf_cz_se q9",  "sf_cz_nw q4",  // X4-D9 2Q gate
+        //                         "sf_cz_se q12", "sf_cz_nw q15", // Z3-D4 2Q gate
+        //                         "sf_park q2",  // D2 park
+        //                         "sf_park q7",  // Z1 park
+        //                         "sf_park q14", // Z2 park
+        //                         "sf_park q3",  // X3 park
+        //                         "sf_park q5",  // D8 park
+        //                         "barrier q8, q0, q10, q16, q9, q4, q12, q15, q2, q7, q14, q3, q5",
+        //                         "update_ph_se q8",
+        //                         "update_ph_nw q0",
+        //                         "update_ph_se q10",
+        //                         "update_ph_nw q16",
+        //                         "update_ph_se q9",
+        //                         "update_ph_nw q4",
+        //                         "update_ph_se q12",
+        //                         "update_ph_nw q15",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+
+        // fluxing steps for parity checks in a distance_3 repetition code. Using Z1 and Z3 only
+        // "repetition_code_1 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                          "sf_cz_nw q7",  "sf_cz_se q15", // Z1-D4 2Q gate
+        //                          "sf_park q3",  // X3 park
+        //                          "sf_park q12", // Z3 park
+        //                          "barrier q11, q6, q7, q15, q3, q1, q14, q16, q2, q8, q10, q5, q12",
+        //                          "update_ph_nw q7",
+        //                          "update_ph_se q15",
+        //                          "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_2 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                          "sf_cz_ne q7",  "sf_cz_sw q13", // Z1-D5 2Q gate
+        //                          "sf_park q8",  // X2 park
+        //                          "sf_park q10", // Z4 park
+        //                          "sf_park q3",  // X3 park
+        //                          "barrier q11, q2, q7, q13, q3, q5, q0, q14, q6, q8, q10, q1",
+        //                          "update_ph_ne q7",
+        //                          "update_ph_sw q13",
+        //                          "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_3 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                          "sf_cz_sw q1",  "sf_cz_ne q12", // D7-Z3 2Q gate
+        //                          "sf_park q3",  // X3 park
+        //                          "barrier q8, q2, q10, q13, q9, q5, q1, q12, q0, q7, q3, q4",
+        //                          "update_ph_sw q1",
+        //                          "update_ph_ne q12",
+        //                          "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_4 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                          "sf_cz_se q12", "sf_cz_nw q15", // Z3-D4 2Q gate
+        //                          "sf_park q7",  // Z1 park
+        //                          "sf_park q3",  // X3 park
+        //                          "barrier q8, q0, q10, q16, q9, q4, q12, q15, q2, q7, q14, q3, q5",
+        //                          "update_ph_se q12",
+        //                          "update_ph_nw q15",
+        //                          "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+        
+
+        // fluxing steps for parity checks in a distance_3 repetition code. Using Z4 and Z2 only
         "repetition_code_1 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
-                                 "sf_cz_nw q11", "sf_cz_se q6",  // X1-D1 2Q gate
-                                 "sf_cz_nw q7",  "sf_cz_se q15", // Z1-D4 2Q gate
-                                 "sf_cz_nw q3",  "sf_cz_se q1",  // X3-D7 2Q gate
-                                 "sf_cz_nw q14", "sf_cz_se q16", // Z2-D6 2Q gate
-                                 "sf_park q2",  // D2 park
-                                 "sf_park q8",  // X2 park
-                                 "sf_park q10", // Z4 park
-                                 "sf_park q5",  // D8 park
-                                 "sf_park q12", // Z3 park
-                                 "barrier q11, q6, q7, q15, q3, q1, q14, q16, q2, q8, q10, q5, q12",
-                                 "update_ph_nw q11",
-                                 "update_ph_se q6",
-                                 "update_ph_nw q7",
-                                 "update_ph_se q15",
-                                 "update_ph_nw q3",
-                                 "update_ph_se q1",
-                                 "update_ph_nw q14",
-                                 "update_ph_se q16",
-                                 "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+                                "sf_cz_nw q14", "sf_cz_se q16", // Z2-D6 2Q gate
+                                "sf_park q8",  // X2 park
+                                "sf_park q10", // Z4 park
+                                "barrier q11, q6, q7, q15, q3, q1, q14, q16, q2, q8, q10, q5, q12",
+                                "update_ph_nw q14",
+                                "update_ph_se q16",
+                                "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
 
         "repetition_code_2 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
-                                 "sf_cz_ne q11", "sf_cz_sw q2",  // X1-D2 2Q gate
-                                 "sf_cz_ne q7",  "sf_cz_sw q13", // Z1-D5 2Q gate
-                                 "sf_cz_ne q3",  "sf_cz_sw q5",  // X3-D8 2Q gate
-                                 "sf_cz_ne q0",  "sf_cz_sw q14", // Z2-D3 2Q gate
-                                 "sf_park q6",  // D1 park
-                                 "sf_park q8",  // X2 park
-                                 "sf_park q10", // Z4 park
-                                 "sf_park q1",  // D7 park
-                                 "barrier q11, q2, q7, q13, q3, q5, q0, q14, q6, q8, q10, q1",
-                                 "update_ph_ne q11",
-                                 "update_ph_sw q2",
-                                 "update_ph_ne q7",
-                                 "update_ph_sw q13",
-                                 "update_ph_ne q3",
-                                 "update_ph_sw q5",
-                                 "update_ph_ne q0",
-                                 "update_ph_sw q14",
-                                 "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+                                "sf_cz_ne q0",  "sf_cz_sw q14", // Z2-D3 2Q gate
+                                "barrier q11, q2, q7, q13, q3, q5, q0, q14, q6, q8, q10, q1",
+                                "update_ph_ne q0",
+                                "update_ph_sw q14",
+                                "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
 
         "repetition_code_3 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
-                                 "sf_cz_sw q8",  "sf_cz_ne q2",  // X2-D2 2Q gate
-                                 "sf_cz_sw q10", "sf_cz_ne q13", // Z4-D5 2Q gate
-                                 "sf_cz_sw q9",  "sf_cz_ne q5",  // X4-D8 2Q gate
-                                 "sf_cz_sw q1",  "sf_cz_ne q12", // D7-Z3 2Q gate
-                                 "sf_park q0",  // D3 park
-                                 "sf_park q7",  // Z1 park
-                                 "sf_park q3",  // X3 park
-                                 "sf_park q4",  // D9 park
-                                 "barrier q8, q2, q10, q13, q9, q5, q1, q12, q0, q7, q3, q4",
-                                 "update_ph_sw q8",
-                                 "update_ph_ne q2",
-                                 "update_ph_sw q10",
-                                 "update_ph_ne q13",
-                                 "update_ph_sw q9",
-                                 "update_ph_ne q5",
-                                 "update_ph_sw q1",
-                                 "update_ph_ne q12",
-                                 "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+                                "sf_cz_sw q10", "sf_cz_ne q13", // Z4-D5 2Q gate
+                                "sf_park q7",  // Z1 park
+                                "sf_park q3",  // X3 park
+                                "sf_park q8",  // X2 park
+                                "barrier q8, q2, q10, q13, q9, q5, q1, q12, q0, q7, q3, q4",
+                                "update_ph_sw q10",
+                                "update_ph_ne q13",
+                                "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
 
         "repetition_code_4 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
-                                 "sf_cz_se q8",  "sf_cz_nw q0",  // X2-D3 2Q gate
-                                 "sf_cz_se q10", "sf_cz_nw q16", // Z4-D6 2Q gate
-                                 "sf_cz_se q9",  "sf_cz_nw q4",  // X4-D9 2Q gate
-                                 "sf_cz_se q12", "sf_cz_nw q15", // Z3-D4 2Q gate
-                                 "sf_park q2",  // D2 park
-                                 "sf_park q7",  // Z1 park
-                                 "sf_park q14", // Z2 park
-                                 "sf_park q3",  // X3 park
-                                 "sf_park q5",  // D8 park
-                                 "barrier q8, q0, q10, q16, q9, q4, q12, q15, q2, q7, q14, q3, q5",
-                                 "update_ph_se q8",
-                                 "update_ph_nw q0",
-                                 "update_ph_se q10",
-                                 "update_ph_nw q16",
-                                 "update_ph_se q9",
-                                 "update_ph_nw q4",
-                                 "update_ph_se q12",
-                                 "update_ph_nw q15",
-                                 "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+                                "sf_cz_se q10", "sf_cz_nw q16", // Z4-D6 2Q gate
+                                "sf_park q8",  // X2 park
+                                "sf_park q14", // Z2 park
+                                "barrier q8, q0, q10, q16, q9, q4, q12, q15, q2, q7, q14, q3, q5",
+                                "update_ph_se q10",
+                                "update_ph_nw q16",
+                                "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // // fluxing steps for parity checks in a distance_3 repetition code. Using Z1-Z4 only.
+        // "repetition_code_1 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_nw q7",  "sf_cz_se q15", // Z1-D4 2Q gate
+        //                         "sf_park q3",  // X3 park
+        //                         "sf_park q12", // Z3 park
+        //                         "barrier q11, q6, q7, q15, q3, q1, q14, q16, q2, q8, q10, q5, q12",
+        //                         "update_ph_nw q7",
+        //                         "update_ph_se q15",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_2 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_ne q7",  "sf_cz_sw q13", // Z1-D5 2Q gate
+        //                         "sf_park q3",  // X3 park
+        //                         "sf_park q8",  // X2 park
+        //                         "sf_park q10", // Z4 park
+        //                         "barrier q11, q2, q7, q13, q3, q5, q0, q14, q6, q8, q10, q1",
+        //                         "update_ph_ne q7",
+        //                         "update_ph_sw q13",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_3 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_sw q10", "sf_cz_ne q13", // Z4-D5 2Q gate
+        //                         "sf_park q8",  // X2 park
+        //                         "sf_park q7",  // Z1 park
+        //                         "sf_park q3",  // X3 park
+        //                         "barrier q8, q2, q10, q13, q9, q5, q1, q12, q0, q7, q3, q4",
+        //                         "update_ph_sw q10",
+        //                         "update_ph_ne q13",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_4 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_se q10", "sf_cz_nw q16", // Z4-D6 2Q gate
+        //                         "sf_park q14", // Z2 park
+        //                         "sf_park q8",  // X2 park
+        //                         "barrier q8, q0, q10, q16, q9, q4, q12, q15, q2, q7, q14, q3, q5",
+        //                         "update_ph_se q10",
+        //                         "update_ph_nw q16",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+
+        // fluxing steps for parity checks in a distance_5 repetition code
+        // "repetition_code_1 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_nw q7",  "sf_cz_se q15", // Z1-D4 2Q gate
+        //                         "sf_cz_nw q14", "sf_cz_se q16", // Z2-D6 2Q gate
+        //                         "sf_park q3",  // X3 park
+        //                         "sf_park q12", // Z3 park
+        //                         "sf_park q8",  // X2 park
+        //                         "sf_park q10", // Z4 park
+        //                         "barrier q11, q6, q7, q15, q3, q1, q14, q16, q2, q8, q10, q5, q12",
+        //                         "update_ph_nw q7",
+        //                         "update_ph_se q15",
+        //                         "update_ph_nw q14",
+        //                         "update_ph_se q16",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_2 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_ne q7",  "sf_cz_sw q13", // Z1-D5 2Q gate
+        //                         "sf_cz_ne q0",  "sf_cz_sw q14", // Z2-D3 2Q gate
+        //                         "sf_park q8",  // X2 park
+        //                         "sf_park q3",  // X3 park
+        //                         "sf_park q10", // Z4 park
+        //                         "barrier q11, q2, q7, q13, q3, q5, q0, q14, q6, q8, q10, q1",
+        //                         "update_ph_ne q7",
+        //                         "update_ph_sw q13",
+        //                         "update_ph_ne q0",
+        //                         "update_ph_sw q14",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_3 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_sw q10", "sf_cz_ne q13", // Z4-D5 2Q gate
+        //                         "sf_cz_sw q1",  "sf_cz_ne q12", // D7-Z3 2Q gate
+        //                         "sf_park q7",  // Z1 park
+        //                         "sf_park q8",  // X2 park
+        //                         "sf_park q3",  // X3 park
+        //                         "barrier q8, q2, q10, q13, q9, q5, q1, q12, q0, q7, q3, q4",
+        //                         "update_ph_sw q10",
+        //                         "update_ph_ne q13",
+        //                         "update_ph_sw q1",
+        //                         "update_ph_ne q12",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_4 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_se q10", "sf_cz_nw q16", // Z4-D6 2Q gate
+        //                         "sf_cz_se q12", "sf_cz_nw q15", // Z3-D4 2Q gate
+        //                         "sf_park q7",  // Z1 park
+        //                         "sf_park q14", // Z2 park
+        //                         "sf_park q8",  // X2 park
+        //                         "sf_park q3",  // X3 park
+        //                         "barrier q8, q0, q10, q16, q9, q4, q12, q15, q2, q7, q14, q3, q5",
+        //                         "update_ph_se q10",
+        //                         "update_ph_nw q16",
+        //                         "update_ph_se q12",
+        //                         "update_ph_nw q15",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+
+
+        // // fluxing steps for parity checks in a distance_7 repetition code (X3+X2)
+        // "repetition_code_1 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_nw q7",  "sf_cz_se q15", // Z1-D4 2Q gate
+        //                         "sf_cz_nw q3",  "sf_cz_se q1",  // X3-D7 2Q gate
+        //                         "sf_cz_nw q14", "sf_cz_se q16", // Z2-D6 2Q gate
+        //                         "sf_park q8",  // X2 park
+        //                         "sf_park q10", // Z4 park
+        //                         "sf_park q5",  // D8 park
+        //                         "sf_park q12", // Z3 park
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "update_ph_nw q7",
+        //                         "update_ph_se q15",
+        //                         "update_ph_nw q3",
+        //                         "update_ph_se q1",
+        //                         "update_ph_nw q14",
+        //                         "update_ph_se q16",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_2 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_ne q7",  "sf_cz_sw q13", // Z1-D5 2Q gate
+        //                         "sf_cz_ne q3",  "sf_cz_sw q5",  // X3-D8 2Q gate
+        //                         "sf_cz_ne q0",  "sf_cz_sw q14", // Z2-D3 2Q gate
+        //                         "sf_park q8",  // X2 park
+        //                         "sf_park q10", // Z4 park
+        //                         "sf_park q1",  // D7 park
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "update_ph_ne q7",
+        //                         "update_ph_sw q13",
+        //                         "update_ph_ne q3",
+        //                         "update_ph_sw q5",
+        //                         "update_ph_ne q0",
+        //                         "update_ph_sw q14",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_3 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_sw q8",  "sf_cz_ne q2",  // X2-D2 2Q gate
+        //                         "sf_cz_sw q10", "sf_cz_ne q13", // Z4-D5 2Q gate
+        //                         "sf_cz_sw q1",  "sf_cz_ne q12", // D7-Z3 2Q gate
+        //                         "sf_park q0",  // D3 park
+        //                         "sf_park q7",  // Z1 park
+        //                         "sf_park q3",  // X3 park                               
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "update_ph_sw q8",
+        //                         "update_ph_ne q2",
+        //                         "update_ph_sw q10",
+        //                         "update_ph_ne q13",
+        //                         "update_ph_sw q1",
+        //                         "update_ph_ne q12",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_4 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_se q8",  "sf_cz_nw q0",  // X2-D3 2Q gate
+        //                         "sf_cz_se q10", "sf_cz_nw q16", // Z4-D6 2Q gate
+        //                         "sf_cz_se q12", "sf_cz_nw q15", // Z3-D4 2Q gate
+        //                         "sf_park q2",  // D2 park
+        //                         "sf_park q7",  // Z1 park
+        //                         "sf_park q14", // Z2 park
+        //                         "sf_park q3",  // X3 park                               
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "update_ph_se q8",
+        //                         "update_ph_nw q0",
+        //                         "update_ph_se q10",
+        //                         "update_ph_nw q16",
+        //                         "update_ph_se q12",
+        //                         "update_ph_nw q15",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+
+
+
+        // fluxing steps for parity checks in a distance_7 repetition code (X3+X2) with modified parkings
+        // "repetition_code_1 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_nw q7",  "sf_cz_se q15", // Z1-D4 2Q gate
+        //                         "sf_cz_nw q3",  "sf_cz_se q1",  // X3-D7 2Q gate
+        //                         "sf_cz_nw q14", "sf_cz_se q16", // Z2-D6 2Q gate
+        //                         "sf_park q8",  // X2 park
+        //                         "sf_park q10", // Z4 park
+        //                         "sf_park q5",  // D8 park
+        //                         "sf_park q12", // Z3 park
+        //                         // "sf_park q4",  // D9 park
+        //                         "sf_park q0", // D3 park
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "update_ph_nw q7",
+        //                         "update_ph_se q15",
+        //                         "update_ph_nw q3",
+        //                         "update_ph_se q1",
+        //                         "update_ph_nw q14",
+        //                         "update_ph_se q16",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_2 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_ne q7",  "sf_cz_sw q13", // Z1-D5 2Q gate
+        //                         "sf_cz_ne q3",  "sf_cz_sw q5",  // X3-D8 2Q gate
+        //                         "sf_cz_ne q0",  "sf_cz_sw q14", // Z2-D3 2Q gate
+        //                         "sf_park q8",  // X2 park
+        //                         "sf_park q10", // Z4 park
+        //                         "sf_park q1",  // D7 park
+        //                         // "sf_park q4",  // D9 park
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "update_ph_ne q7",
+        //                         "update_ph_sw q13",
+        //                         "update_ph_ne q3",
+        //                         "update_ph_sw q5",
+        //                         "update_ph_ne q0",
+        //                         "update_ph_sw q14",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_3 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_sw q8",  "sf_cz_ne q2",  // X2-D2 2Q gate
+        //                         "sf_cz_sw q10", "sf_cz_ne q13", // Z4-D5 2Q gate
+        //                         "sf_cz_sw q1",  "sf_cz_ne q12", // D7-Z3 2Q gate
+        //                         "sf_park q0",  // D3 park
+        //                         "sf_park q7",  // Z1 park
+        //                         "sf_park q3",  // X3 park
+        //                         "sf_park q6",  // D1 park
+        //                         // "sf_park q4",  // D9 park                              
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "update_ph_sw q8",
+        //                         "update_ph_ne q2",
+        //                         "update_ph_sw q10",
+        //                         "update_ph_ne q13",
+        //                         "update_ph_sw q1",
+        //                         "update_ph_ne q12",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_4 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_se q8",  "sf_cz_nw q0",  // X2-D3 2Q gate
+        //                         "sf_cz_se q10", "sf_cz_nw q16", // Z4-D6 2Q gate
+        //                         "sf_cz_se q12", "sf_cz_nw q15", // Z3-D4 2Q gate
+        //                         "sf_park q2",  // D2 park
+        //                         "sf_park q7",  // Z1 park
+        //                         "sf_park q14", // Z2 park
+        //                         "sf_park q3",  // X3 park
+        //                         "sf_park q6",  // D1 park
+        //                         // "sf_park q4",  // D9 park                                
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "update_ph_se q8",
+        //                         "update_ph_nw q0",
+        //                         "update_ph_se q10",
+        //                         "update_ph_nw q16",
+        //                         "update_ph_se q12",
+        //                         "update_ph_nw q15",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+
+        // // fluxing steps for parity checks in a distance_7 repetition code (X3+X2) with 6 steps
+        // "repetition_code_1 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_nw q3",  "sf_cz_se q1",  // X3-D7 2Q gate
+        //                         "sf_cz_nw q14", "sf_cz_se q16", // Z2-D6 2Q gate
+        //                         "sf_park q8",  // X2 park
+        //                         "sf_park q10", // Z4 park
+        //                         "sf_park q5",  // D8 park
+        //                         // "sf_park q4",  // D9 park
+        //                         "sf_park q0",  // D3 park
+        //                         "sf_park q2",  // D2 park, check if this is needed?
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "update_ph_nw q3",
+        //                         "update_ph_se q1",
+        //                         "update_ph_nw q14",
+        //                         "update_ph_se q16",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_2 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_ne q3",  "sf_cz_sw q5",  // X3-D8 2Q gate
+        //                         "sf_cz_ne q0",  "sf_cz_sw q14", // Z2-D3 2Q gate
+        //                         "sf_park q1",  // D7 park
+        //                         // "sf_park q4",  // D9 park
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "update_ph_ne q3",
+        //                         "update_ph_sw q5",
+        //                         "update_ph_ne q0",
+        //                         "update_ph_sw q14",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_3 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         // "sf_cz_nw q7",  "sf_cz_se q15", // Z1-D4 2Q gate
+        //                         "sf_cz_sw q8",  "sf_cz_ne q2",  // X2-D2 2Q gate
+        //                         // "sf_cz_sw q10", "sf_cz_ne q13", // Z4-D5 2Q gate
+        //                         // "sf_park q12", // Z3 park
+        //                         // "sf_park q3",  // X3 park
+        //                         // "sf_park q4",  // D9 park
+        //                         "sf_park q0",     // D3 park
+        //                         // "sf_park q7",     // Z1 park
+        //                         // "sf_park q3",     // X3 park
+        //                         // "sf_park q6",     // D1 park
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "update_ph_nw q7",
+        //                         "update_ph_se q15",
+        //                         "update_ph_sw q8",
+        //                         "update_ph_ne q2",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_4 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         // "sf_cz_ne q7",  "sf_cz_sw q13", // Z1-D5 2Q gate
+        //                         "sf_cz_se q8",  "sf_cz_nw q0",  // X2-D3 2Q gate   
+        //                         // "sf_cz_se q10", "sf_cz_nw q16", // Z4-D6 2Q gate                             
+        //                         // "sf_park q10", // Z4 park
+        //                         // "sf_park q3",  // X3 park
+        //                         "sf_park q2",     // D2 park
+        //                         // "sf_park q4",  // D9 park
+        //                         // "sf_park q7",     // Z1 park
+        //                         // "sf_park q3",     // X3 park
+        //                         // "sf_park q6",     // D1 park
+        //                         // "sf_park q14",    // Z2 park
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "update_ph_ne q7",
+        //                         "update_ph_sw q13",
+        //                         "update_ph_se q8",
+        //                         "update_ph_nw q0",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+
+        // "repetition_code_5 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_sw q10", "sf_cz_ne q13", // Z4-D5 2Q gate
+        //                         "sf_cz_sw q1",  "sf_cz_ne q12", // D7-Z3 2Q gate
+        //                         "sf_park q7",  // Z1 park
+        //                         "sf_park q3",  // X3 park
+        //                         "sf_park q8",  // X2 park
+        //                         "sf_park q2",  // D2 park
+        //                         "sf_park q0",  // D3 park
+        //                         // "sf_park q4",  // D9 park                              
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "update_ph_sw q10",
+        //                         "update_ph_ne q13",
+        //                         "update_ph_sw q1",
+        //                         "update_ph_ne q12",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_6 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_se q10", "sf_cz_nw q16", // Z4-D6 2Q gate
+        //                         "sf_cz_se q12", "sf_cz_nw q15", // Z3-D4 2Q gate
+        //                         "sf_park q7",  // Z1 park
+        //                         "sf_park q3",  // X3 park
+        //                         "sf_park q14", // Z2 park
+        //                         "sf_park q8",  // X2 park
+        //                         "sf_park q2",  // D2 park
+        //                         "sf_park q0",  // D3 park
+        //                         // "sf_park q4",  // D9 park                                
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "update_ph_se q10",
+        //                         "update_ph_nw q16",
+        //                         "update_ph_se q12",
+        //                         "update_ph_nw q15",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+
+
+
+
+        // // fluxing steps for parity checks in a distance_7 repetition code (X1+X2)
+        // "repetition_code_1 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_nw q11", "sf_cz_se q6",  // X1-D1 2Q gate
+        //                         "sf_cz_nw q7",  "sf_cz_se q15", // Z1-D4 2Q gate
+        //                         "sf_cz_nw q14", "sf_cz_se q16", // Z2-D6 2Q gate
+        //                         "sf_park q2",  // D2 park
+        //                         "sf_park q8",  // X2 park
+        //                         "sf_park q10", // Z4 park
+        //                         "sf_park q12", // Z3 park
+        //                         "sf_park q3",  // X3 park
+        //                         "barrier q11, q6, q7, q15, q3, q1, q14, q16, q2, q8, q10, q5, q12",
+        //                         "update_ph_nw q11",
+        //                         "update_ph_se q6",
+        //                         "update_ph_nw q7",
+        //                         "update_ph_se q15",
+        //                         "update_ph_nw q14",
+        //                         "update_ph_se q16",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_2 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_ne q11", "sf_cz_sw q2",  // X1-D2 2Q gate
+        //                         "sf_cz_ne q7",  "sf_cz_sw q13", // Z1-D5 2Q gate
+        //                         "sf_cz_ne q0",  "sf_cz_sw q14", // Z2-D3 2Q gate
+        //                         "sf_park q6",  // D1 park
+        //                         "sf_park q8",  // X2 park
+        //                         "sf_park q10", // Z4 park
+        //                         "sf_park q3",  // X3 park
+        //                         "barrier q11, q2, q7, q13, q3, q5, q0, q14, q6, q8, q10, q1",
+        //                         "update_ph_ne q11",
+        //                         "update_ph_sw q2",
+        //                         "update_ph_ne q7",
+        //                         "update_ph_sw q13",
+        //                         "update_ph_ne q0",
+        //                         "update_ph_sw q14",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_3 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_sw q8",  "sf_cz_ne q2",  // X2-D2 2Q gate
+        //                         "sf_cz_sw q10", "sf_cz_ne q13", // Z4-D5 2Q gate
+        //                         "sf_cz_sw q1",  "sf_cz_ne q12", // D7-Z3 2Q gate
+        //                         "sf_park q0",  // D3 park
+        //                         "sf_park q7",  // Z1 park
+        //                         "sf_park q3",  // X3 park
+        //                         "barrier q8, q2, q10, q13, q9, q5, q1, q12, q0, q7, q3, q4",
+        //                         "update_ph_sw q8",
+        //                         "update_ph_ne q2",
+        //                         "update_ph_sw q10",
+        //                         "update_ph_ne q13",
+        //                         "update_ph_sw q1",
+        //                         "update_ph_ne q12",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+        // "repetition_code_4 q0": ["barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25",
+        //                         "sf_cz_se q8",  "sf_cz_nw q0",  // X2-D3 2Q gate
+        //                         "sf_cz_se q10", "sf_cz_nw q16", // Z4-D6 2Q gate
+        //                         "sf_cz_se q12", "sf_cz_nw q15", // Z3-D4 2Q gate
+        //                         "sf_park q2",  // D2 park
+        //                         "sf_park q7",  // Z1 park
+        //                         "sf_park q14", // Z2 park
+        //                         "sf_park q3",  // X3 park
+        //                         "barrier q8, q0, q10, q16, q9, q4, q12, q15, q2, q7, q14, q3, q5",
+        //                         "update_ph_se q8",
+        //                         "update_ph_nw q0",
+        //                         "update_ph_se q10",
+        //                         "update_ph_nw q16",
+        //                         "update_ph_se q12",
+        //                         "update_ph_nw q15",
+        //                         "barrier q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21, q22, q23, q24, q25"],
+
+
+
+
 
         // CC additions
         "cnot_park1 %0 %1 %2": ["ry90 %1", "cz %0 %1", "park_cz %2", "ry90 %1"],
diff --git a/pycqed/measurement/openql_experiments/multi_qubit_oql.py b/pycqed/measurement/openql_experiments/multi_qubit_oql.py
index 14e22ed83..b71e0779a 100644
--- a/pycqed/measurement/openql_experiments/multi_qubit_oql.py
+++ b/pycqed/measurement/openql_experiments/multi_qubit_oql.py
@@ -1926,6 +1926,9 @@ def parity_check_ramsey(
             k.barrier([])
             # Single qubit gates
             for j, state in enumerate(case):
+                if state == '2':
+                    k.gate("rx12", [Q_idxs_control[j]])
+                    k.gate("rxm180", [Q_idxs_control[j]])
                 if state == '1':
                     k.gate("rxm180", [Q_idxs_control[j]])
             # cw_idx corresponds to special hardcoded angles in the lutman
@@ -4321,9 +4324,17 @@ def repeated_stabilizer_data_measurement_sequence(
                 k.prepz(q)
                 k.measure(q)
             for q in _remaining_ancillas+[Q_anc]:
+                k.prepz(q)
                 k.measure(q)
             k.gate('wait', [], 400)  # to avoid UHF trigger holdoff!
             # TODO: Add arbitrary state initialization
+            if initial_state_qubits:
+                if stabilizer_type == 'Z':
+                    for q in Q_D:
+                        k.gate("ry90", [q])
+                if stabilizer_type == 'X':
+                    for q in Q_D:
+                        k.gate("i", [q])
             k.barrier([])
             # QEC Rounds 
             for i in range(n_rounds):
@@ -4348,13 +4359,13 @@ def repeated_stabilizer_data_measurement_sequence(
                              ['X3', 'Z1', 'X1', 'Z2']]:
                         k.gate("rym90", [q])
                     # First Pi/2 pulse
-                    # elif q in [ancilla_qubit_map[qa] for qa in \
-                    #            ['Z3', 'X4', 'Z4', 'X2']]:
-                    #     k.gate("ry90", [q])
+                    elif q in [ancilla_qubit_map[qa] for qa in \
+                               ['Z3', 'X4', 'Z4', 'X2']]:
+                        k.gate("ry90", [q])
                 # Flux dance
-                # k.gate('wait', [])
-                # k.gate(f'repetition_code_3', [0])
-                # k.gate(f'repetition_code_4', [0])
+                k.gate('wait', [])
+                k.gate(f'repetition_code_3', [0])
+                k.gate(f'repetition_code_4', [0])
                 k.gate('wait', [])
                 # Second Pi/2 pulse
                 if stabilizer_type == 'X':
@@ -4446,6 +4457,7 @@ def repetition_code_sequence(
         array_of_round_number: List[int],
         stabilizer_type: str = 'X',
         measurement_time_ns: int = 500,
+        initial_state: List[int] = None,
         ):
     """OpenQL program constructor for repetition code experiement (any distance).
 
@@ -4501,20 +4513,249 @@ def repetition_code_sequence(
             k.measure(q)
         k.gate('wait', [], 400)  # to avoid UHF trigger holdoff!
         # Arbitrary state initialization
-        if stabilizer_type == 'Z':
-            for i, q in enumerate(involved_data_indices):
-                if i % 2 == 0:
-                    k.gate("i", [q])
-                else:
-                    k.gate("rx180", [q])
-                    # k.gate("i", [q])
-        if stabilizer_type == 'X':
-            for i, q in enumerate(involved_data_indices):
-                if i % 2 == 0:
+        if initial_state is None:
+            if stabilizer_type == 'Z':
+                for i, q in enumerate(involved_data_indices):
+                    if i % 2 == 0:
+                        k.gate("i", [q])
+                    else:
+                        k.gate("rx180", [q])
+                        # k.gate("i", [q])
+            if stabilizer_type == 'X':
+                for i, q in enumerate(involved_data_indices):
+                    if i % 2 == 0:
+                        k.gate("ry90", [q])
+                    else:
+                        k.gate("rym90", [q])
+                        # k.gate("i", [q])
+        else:
+            if stabilizer_type == 'Z':
+                for i, q in enumerate(involved_data_indices):
+                    if initial_state[i] == 0:
+                        k.gate("i", [q])
+                    else:
+                        k.gate("rx180", [q])
+                        # k.gate("i", [q])
+            if stabilizer_type == 'X':
+                for i, q in enumerate(involved_data_indices):
+                    if initial_state[i] == 0:
+                        k.gate("ry90", [q])
+                    else:
+                        k.gate("rym90", [q])
+                        # k.gate("i", [q])
+        k.barrier([])
+        # QEC Rounds 
+        for i in range(n_rounds):
+            # First Pi/2 pulse
+            if stabilizer_type == 'X':
+                for q in involved_data_indices:
                     k.gate("ry90", [q])
-                else:
+            # for q in Z_anci_idxs+X_anci_idxs:
+            for q in involved_ancilla_indices:
+                if qubit_id_to_name_lookup[q] in ['X3', 'Z2']:
+                    k.gate("ry90", [q])
+            # Flux dance
+            k.gate('wait', [])
+            k.gate(f'repetition_code_1', [0])
+            k.gate(f'repetition_code_2', [0])
+            k.gate('wait', [])
+            # Second Pi/2 pulse
+            # for q in Z_anci_idxs+X_anci_idxs:
+            for q in involved_ancilla_indices:
+                if qubit_id_to_name_lookup[q] in ['X3', 'Z2']:
+                    k.gate("rym90", [q])
+                # First Pi/2 pulse
+                elif qubit_id_to_name_lookup[q] in ['Z1', 'X2']:
+                    k.gate("ry90", [q])
+            # Flux dance
+            k.gate('wait', [])
+            k.gate(f'repetition_code_3', [0])
+            k.gate(f'repetition_code_4', [0])
+            k.gate('wait', [])
+            # Second Pi/2 pulse
+            # for q in Z_anci_idxs+X_anci_idxs:
+            for q in involved_ancilla_indices:
+                if qubit_id_to_name_lookup[q] in ['Z1', 'X2']:
+                    k.gate("rym90", [q])
+                # First Pi/2 pulse
+                elif qubit_id_to_name_lookup[q] in ['Z3', 'Z4']:
+                    k.gate("ry90", [q])
+            # Flux dance
+            k.gate('wait', [])
+            k.gate(f'repetition_code_5', [0])
+            k.gate(f'repetition_code_6', [0])
+            k.gate('wait', [])
+            # Second Pi/2 pulse
+            if stabilizer_type == 'X':
+                for q in involved_data_indices:
+                    k.gate("rym90", [q])
+            # for q in Z_anci_idxs+X_anci_idxs:
+            for q in involved_ancilla_indices:
+                if qubit_id_to_name_lookup[q] in ['Z3', 'Z4']:
                     k.gate("rym90", [q])
-                    # k.gate("i", [q])
+            k.gate('wait', [])
+            
+            # During any other round, measure only ancilla's and decouple data qubits.
+            at_last_round: bool = i == n_rounds-1
+            if not at_last_round:
+                # Measure (only) all ancilla's, dynamical decoupling on data qubits
+                for q in all_ancilla_indices:
+                    k.measure(q)
+                for q in involved_data_indices:
+                    # Single measurement Echo
+                    idle_time = (measurement_time_ns-20)//2
+                    nr_idles = idle_time//20
+                    for idle in range(nr_idles):
+                        k.gate('i', [q])
+                    k.gate('rx180', [q])
+                    for idle in range(nr_idles):
+                        k.gate('i', [q])
+                k.gate("wait", [], 0)
+
+        # Make sure all qubits are measured in the last round
+        # Before last round apply correction gate to data qubits, depending on the stabilizer type.
+        # Final measurement and data qubit correction
+        if stabilizer_type == 'X':
+            for q in involved_data_indices:
+                k.gate('rym90', [q])
+        for q in all_qubit_indices:
+            k.measure(q)
+        k.gate("wait", [], 0)
+        p.add_kernel(k)
+
+    ######################
+    # Calibration points #
+    ######################
+    # Calibration 000
+    k = p.create_kernel('Cal_zeros')
+    for q in all_qubit_indices:
+        k.prepz(q)
+        k.measure(q)
+    k.gate('wait', [], 400)
+    for q in all_qubit_indices:
+        k.measure(q)
+    p.add_kernel(k)
+    # Calibration 111
+    k = p.create_kernel('Cal_ones')
+    for q in all_qubit_indices:
+        k.prepz(q)
+        k.measure(q)
+    k.gate('wait', [], 400)
+    for q in all_qubit_indices:
+        k.gate('rx180', [q])
+        k.measure(q)
+    p.add_kernel(k)
+    # Calibration 222
+    k = p.create_kernel('Cal_twos')
+    for q in all_qubit_indices:
+        k.prepz(q)
+        k.measure(q)
+    k.gate('wait', [], 400)
+    for q in all_qubit_indices:
+        k.gate('rx180', [q])
+        k.gate('rx12', [q])
+        k.measure(q)
+    p.add_kernel(k)
+
+    p.compile()
+    return p
+
+
+def repetition_code_sequence_old(
+        involved_ancilla_indices: List[int],
+        involved_data_indices: List[int],
+        all_ancilla_indices: List[int],
+        all_data_indices: List[int],
+        platf_cfg: str,
+        array_of_round_number: List[int],
+        stabilizer_type: str = 'X',
+        measurement_time_ns: int = 500,
+        initial_state: List[int] = None,
+        ):
+    """OpenQL program constructor for repetition code experiement (any distance).
+
+    Args:
+        involved_ancilla_indices (List[int]): Ancilla-qubit indices part of the repetition code. Used by the central controller. Also known as 'cfg_qubit_nr'.
+        involved_data_indices (List[int]): Data-qubit indices part of the repetition code. Used by the central controller. Also known as 'cfg_qubit_nr'.
+        all_ancilla_indices (List[int]): All (ancilla) qubit indices on the device. Used for preparation and calibration points.
+        all_data_indices (List[int]): All (ancilla) qubit indices on the device. Used for preparation and calibration points.
+        platf_cfg (str): Config-file path, required for OqlProgram construction.
+        array_of_round_number (list): Array-like of parity-check repetition number. Example [1, 5, 10] will schedule a stabilizer reptition of x1, x5 and x10 respectively.
+        stabilizer_type (str, optional): _description_. Defaults to 'X'.
+        measurement_time_ns (int, optional): _description_. Defaults to 500.
+
+    Returns:
+        OqlProgram: _description_
+    """
+    p = OqlProgram("Repeated_stabilizer_seq", platf_cfg)
+    assert stabilizer_type in ['X', 'Z'], '"stabilizer_type" must be "X" or "Z"'
+    # Data qubit idx dictionary
+    qubit_id_to_name_lookup: Dict[int, str] = {
+        0: "D3",
+        1: "D7",
+        2: "D2",
+        3: "X3",
+        4: "D9",
+        5: "D8",
+        6: "D1",
+        7: "Z1",
+        8: "X2",
+        9: "X4",
+        10: "Z4",
+        11: "X1",
+        12: "Z3",
+        13: "D5",
+        14: "Z2",
+        15: "D4",
+        16: "D6",
+    }
+    
+    # remove X4 
+    qubit_nr_x4: int = 9
+    if qubit_nr_x4 in all_ancilla_indices:
+        all_ancilla_indices.remove(qubit_nr_x4)
+    
+    all_qubit_indices = all_ancilla_indices + all_data_indices
+
+    for n_rounds in array_of_round_number:
+        
+        k = p.create_kernel(f'Repetition_code_seq_{n_rounds}rounds')        
+        # Preparation & heralded_init
+        for q in all_qubit_indices:
+            k.prepz(q)
+            k.measure(q)
+        k.gate('wait', [], 400)  # to avoid UHF trigger holdoff!
+        # Arbitrary state initialization
+        if initial_state is None:
+            if stabilizer_type == 'Z':
+                for i, q in enumerate(involved_data_indices):
+                    if i % 2 == 0:
+                        k.gate("i", [q])
+                    else:
+                        k.gate("rx180", [q])
+                        # k.gate("i", [q])
+            if stabilizer_type == 'X':
+                for i, q in enumerate(involved_data_indices):
+                    if i % 2 == 0:
+                        k.gate("ry90", [q])
+                    else:
+                        k.gate("rym90", [q])
+                        # k.gate("i", [q])
+        else:
+            if stabilizer_type == 'Z':
+                for i, q in enumerate(involved_data_indices):
+                    if initial_state[i] == 0:
+                        k.gate("i", [q])
+                    else:
+                        k.gate("rx180", [q])
+                        # k.gate("i", [q])
+            if stabilizer_type == 'X':
+                for i, q in enumerate(involved_data_indices):
+                    if initial_state[i] == 0:
+                        k.gate("ry90", [q])
+                    else:
+                        k.gate("rym90", [q])
+                        # k.gate("i", [q])
         k.barrier([])
         # QEC Rounds 
         for i in range(n_rounds):
diff --git a/pycqed/qce_utils/control_interfaces/connectivity_surface_code.py b/pycqed/qce_utils/control_interfaces/connectivity_surface_code.py
index 9c5d5c9b1..fecdd8bea 100644
--- a/pycqed/qce_utils/control_interfaces/connectivity_surface_code.py
+++ b/pycqed/qce_utils/control_interfaces/connectivity_surface_code.py
@@ -4,6 +4,7 @@
 from dataclasses import dataclass, field
 import warnings
 from typing import List, Union, Dict, Tuple
+from enum import Enum, unique, auto
 from pycqed.qce_utils.definitions import SingletonABCMeta
 from pycqed.qce_utils.custom_exceptions import ElementNotIncludedException
 from pycqed.qce_utils.control_interfaces.intrf_channel_identifier import (
@@ -20,6 +21,57 @@
     IParityGroup,
     ParityType,
 )
+from pycqed.qce_utils.control_interfaces.intrf_connectivity import (
+    IDeviceLayer
+)
+
+
+@unique
+class FrequencyGroup(Enum):
+    LOW = auto()
+    MID = auto()
+    HIGH = auto()
+
+
+@dataclass(frozen=True)
+class FrequencyGroupIdentifier:
+    """
+    Data class, representing (qubit) frequency group identifier.
+    """
+    _id: FrequencyGroup
+
+    # region Class Properties
+    @property
+    def id(self) -> FrequencyGroup:
+        """:return: Self identifier."""
+        return self._id
+    # endregion
+
+    # region Class Methods
+    def is_equal_to(self, other: 'FrequencyGroupIdentifier') -> bool:
+        """:return: Boolean, whether other frequency group identifier is equal self."""
+        return self.id == other.id
+
+    def is_higher_than(self, other: 'FrequencyGroupIdentifier') -> bool:
+        """:return: Boolean, whether other frequency group identifier is 'lower' than self."""
+        # Guard clause, if frequency groups are equal, return False
+        if self.is_equal_to(other):
+            return False
+        if self.id == FrequencyGroup.MID and other.id == FrequencyGroup.LOW:
+            return True
+        if self.id == FrequencyGroup.HIGH:
+            return True
+        return False
+
+    def is_lower_than(self, other: 'FrequencyGroupIdentifier') -> bool:
+        """:return: Boolean, whether other frequency group identifier is 'higher' than self."""
+        # Guard clause, if frequency groups are equal, return False
+        if self.is_equal_to(other):
+            return False
+        if self.is_higher_than(other):
+            return False
+        return True
+    # endregion
 
 
 @dataclass(frozen=True)
@@ -109,6 +161,69 @@ def __post_init__(self):
     # endregion
 
 
+@dataclass(frozen=True)
+class FluxDanceLayer:
+    """
+    Data class, containing directional gates played during 'flux-dance' layer.
+    """
+    _edge_ids: List[IEdgeID]
+    """Non-directional edges, part of flux-dance layer."""
+
+    # region Class Properties
+    @property
+    def qubit_ids(self) -> List[IQubitID]:
+        """:return: All qubit-ID's."""
+        return list(set([qubit_id for edge in self.edge_ids for qubit_id in edge.qubit_ids]))
+
+    @property
+    def edge_ids(self) -> List[IEdgeID]:
+        """:return: Array-like of directional edge identifiers, specific for this flux dance."""
+        return self._edge_ids
+    # endregion
+
+    # region Class Methods
+    def contains(self, element: Union[IQubitID, IEdgeID]) -> bool:
+        """:return: Boolean, whether element is part of flux-dance layer or not."""
+        if element in self.qubit_ids:
+            return True
+        if element in self.edge_ids:
+            return True
+        return False
+
+    def get_involved_edge(self, qubit_id: IQubitID) -> IEdgeID:
+        """:return: Edge in which qubit-ID is involved. If qubit-ID not part of self, raise error."""
+        for edge in self.edge_ids:
+            if edge.contains(element=qubit_id):
+                return edge
+        raise ElementNotIncludedException(f'Element {qubit_id} is not part of self ({self}) and cannot be part of an edge.')
+
+    def get_spectating_qubit_ids(self, device_layer: IDeviceLayer) -> List[IQubitID]:
+        """:return: Direct spectator (nearest neighbor) to qubit-ID's participating in flux-dance."""
+        participating_qubit_ids: List[IQubitID] = self.qubit_ids
+        nearest_neighbor_ids: List[IQubitID] = [neighbor_id for qubit_id in participating_qubit_ids for neighbor_id in device_layer.get_neighbors(qubit_id, order=1)]
+        filtered_nearest_neighbor_ids: List[IQubitID] = list(set([qubit_id for qubit_id in nearest_neighbor_ids if qubit_id not in participating_qubit_ids]))
+        return filtered_nearest_neighbor_ids
+
+    def requires_parking(self, qubit_id: IQubitID, device_layer: ISurfaceCodeLayer) -> bool:
+        """
+        Determines whether qubit-ID is required to park based on participation in flux dance and frequency group.
+        :return: Boolean, whether qubit-ID requires some form of parking.
+        """
+        spectating_qubit_ids: List[IQubitID] = self.get_spectating_qubit_ids(device_layer=device_layer)
+        # Guard clause, if qubit-ID does not spectate the flux-dance, no need for parking
+        if qubit_id not in spectating_qubit_ids:
+            return False
+        # Check if qubit-ID requires parking based on its frequency group ID and active two-qubit gates.
+        frequency_group: FrequencyGroupIdentifier = device_layer.get_frequency_group_identifier(element=qubit_id)
+        # Parking is required if any neighboring qubit from a higher frequency group is part of an edge.
+        neighboring_qubit_ids: List[IQubitID] = device_layer.get_neighbors(qubit=qubit_id, order=1)
+        involved_neighbors: List[IQubitID] = [qubit_id for qubit_id in neighboring_qubit_ids if self.contains(qubit_id)]
+        involved_frequency_groups: List[FrequencyGroupIdentifier] = [device_layer.get_frequency_group_identifier(element=qubit_id) for qubit_id in involved_neighbors]
+        return any([neighbor_frequency_group.is_higher_than(frequency_group) for neighbor_frequency_group in involved_frequency_groups])
+    # endregion
+
+
+
 @dataclass(frozen=True)
 class VirtualPhaseIdentifier(IChannelIdentifier):
     """
@@ -136,6 +251,33 @@ def __eq__(self, other):
     # endregion
 
 
+@dataclass(frozen=True)
+class FluxOperationIdentifier(IChannelIdentifier):
+    """
+    Data class, describing (code-word) identifier for flux operation.
+    """
+    _id: str
+
+    # region Interface Properties
+    @property
+    def id(self) -> str:
+        """:returns: Reference Identifier."""
+        return self._id
+    # endregion
+
+    # region Interface Methods
+    def __hash__(self):
+        """:returns: Identifiable hash."""
+        return self._id.__hash__()
+
+    def __eq__(self, other):
+        """:returns: Boolean if other shares equal identifier, else InterfaceMethodException."""
+        if isinstance(other, FluxOperationIdentifier):
+            return self.id.__eq__(other.id)
+        return False
+    # endregion
+
+
 class Surface17Layer(ISurfaceCodeLayer, metaclass=SingletonABCMeta):
     """
     Singleton class, implementing ISurfaceCodeLayer interface to describe a surface-17 layout.
@@ -215,6 +357,25 @@ class Surface17Layer(ISurfaceCodeLayer, metaclass=SingletonABCMeta):
             _data_qubits=[QubitIDObj('D5'), QubitIDObj('D6'), QubitIDObj('D8'), QubitIDObj('D9')]
         ),
     ]
+    _frequency_group_lookup: Dict[IQubitID, FrequencyGroupIdentifier] = {
+        QubitIDObj('D1'): FrequencyGroupIdentifier(_id=FrequencyGroup.LOW),
+        QubitIDObj('D2'): FrequencyGroupIdentifier(_id=FrequencyGroup.LOW),
+        QubitIDObj('D3'): FrequencyGroupIdentifier(_id=FrequencyGroup.LOW),
+        QubitIDObj('D4'): FrequencyGroupIdentifier(_id=FrequencyGroup.HIGH),
+        QubitIDObj('D5'): FrequencyGroupIdentifier(_id=FrequencyGroup.HIGH),
+        QubitIDObj('D6'): FrequencyGroupIdentifier(_id=FrequencyGroup.HIGH),
+        QubitIDObj('D7'): FrequencyGroupIdentifier(_id=FrequencyGroup.LOW),
+        QubitIDObj('D8'): FrequencyGroupIdentifier(_id=FrequencyGroup.LOW),
+        QubitIDObj('D9'): FrequencyGroupIdentifier(_id=FrequencyGroup.LOW),
+        QubitIDObj('Z1'): FrequencyGroupIdentifier(_id=FrequencyGroup.MID),
+        QubitIDObj('Z2'): FrequencyGroupIdentifier(_id=FrequencyGroup.MID),
+        QubitIDObj('Z3'): FrequencyGroupIdentifier(_id=FrequencyGroup.MID),
+        QubitIDObj('Z4'): FrequencyGroupIdentifier(_id=FrequencyGroup.MID),
+        QubitIDObj('X1'): FrequencyGroupIdentifier(_id=FrequencyGroup.MID),
+        QubitIDObj('X2'): FrequencyGroupIdentifier(_id=FrequencyGroup.MID),
+        QubitIDObj('X3'): FrequencyGroupIdentifier(_id=FrequencyGroup.MID),
+        QubitIDObj('X4'): FrequencyGroupIdentifier(_id=FrequencyGroup.MID),
+    }
 
     # region ISurfaceCodeLayer Interface Properties
     @property
@@ -293,6 +454,10 @@ def contains(self, element: Union[IFeedlineID, IQubitID, IEdgeID]) -> bool:
         if element in self.edge_ids:
             return True
         return False
+    
+    def get_frequency_group_identifier(self, element: IQubitID) -> FrequencyGroupIdentifier:
+        """:return: Frequency group identifier based on qubit-ID."""
+        return self._frequency_group_lookup[element]
     # endregion
 
 
@@ -394,6 +559,52 @@ class Repetition9Layer(ISurfaceCodeLayer, metaclass=SingletonABCMeta):
         DirectionalEdgeIDObj(QubitIDObj('D9'), QubitIDObj('X4')): VirtualPhaseIdentifier('vcz_virtual_q_ph_corr_NW'),
         DirectionalEdgeIDObj(QubitIDObj('X4'), QubitIDObj('D9')): VirtualPhaseIdentifier('vcz_virtual_q_ph_corr_SE'),
     }
+    _flux_dances: List[Tuple[FluxDanceLayer, FluxOperationIdentifier]] = [
+        (
+            FluxDanceLayer(
+                _edge_ids=[
+                    EdgeIDObj(QubitIDObj('X1'), QubitIDObj('D1')),
+                    EdgeIDObj(QubitIDObj('Z1'), QubitIDObj('D4')),
+                    EdgeIDObj(QubitIDObj('X3'), QubitIDObj('D7')),
+                    EdgeIDObj(QubitIDObj('Z2'), QubitIDObj('D6')),
+                ]
+            ),
+            FluxOperationIdentifier(_id='repetition_code_1')
+        ),
+        (
+            FluxDanceLayer(
+                _edge_ids=[
+                    EdgeIDObj(QubitIDObj('X1'), QubitIDObj('D2')),
+                    EdgeIDObj(QubitIDObj('Z1'), QubitIDObj('D5')),
+                    EdgeIDObj(QubitIDObj('X3'), QubitIDObj('D8')),
+                    EdgeIDObj(QubitIDObj('Z2'), QubitIDObj('D3')),
+                ]
+            ),
+            FluxOperationIdentifier(_id='repetition_code_2')
+        ),
+        (
+            FluxDanceLayer(
+                _edge_ids=[
+                    EdgeIDObj(QubitIDObj('Z3'), QubitIDObj('D7')),
+                    EdgeIDObj(QubitIDObj('X4'), QubitIDObj('D8')),
+                    EdgeIDObj(QubitIDObj('Z4'), QubitIDObj('D5')),
+                    EdgeIDObj(QubitIDObj('X2'), QubitIDObj('D2')),
+                ]
+            ),
+            FluxOperationIdentifier(_id='repetition_code_3')
+        ),
+        (
+            FluxDanceLayer(
+                _edge_ids=[
+                    EdgeIDObj(QubitIDObj('Z3'), QubitIDObj('D4')),
+                    EdgeIDObj(QubitIDObj('X4'), QubitIDObj('D9')),
+                    EdgeIDObj(QubitIDObj('Z4'), QubitIDObj('D6')),
+                    EdgeIDObj(QubitIDObj('X2'), QubitIDObj('D3')),
+                ]
+            ),
+            FluxOperationIdentifier(_id='repetition_code_4')
+        ),
+    ]
 
     # region ISurfaceCodeLayer Interface Properties
     @property
@@ -434,6 +645,16 @@ def get_parity_group(self, element: Union[IQubitID, IEdgeID]) -> IParityGroup:
         raise ElementNotIncludedException(f"Element: {element} is not included in any parity group.")
     # endregion
 
+    # region IGateDanceLayer Interface Methods
+    def get_flux_dance_at_round(self, index: int) -> FluxDanceLayer:
+        """:return: Flux-dance object based on round index."""
+        try:
+            flux_dance_layer: FluxDanceLayer = self._flux_dances[index]
+            return flux_dance_layer
+        except:
+            raise ElementNotIncludedException(f"Index: {index} is out of bounds for flux dance of length: {len(self._flux_dances)}.")
+    # endregion
+
     # region IDeviceLayer Interface Methods
     def get_connected_qubits(self, feedline: IFeedlineID) -> List[IQubitID]:
         """:return: Qubit-ID's connected to feedline-ID."""
@@ -456,6 +677,41 @@ def contains(self, element: Union[IFeedlineID, IQubitID, IEdgeID]) -> bool:
     # endregion
 
     # region Class Methods
+    def _get_flux_dance_layer(self, element: IEdgeID) -> FluxDanceLayer:
+        """:return: Flux-dance layer of which edge element is part of."""
+        # Assumes element is part of only a single flux-dance layer
+        for flux_dance_layer, _ in self._flux_dances:
+            if flux_dance_layer.contains(element=element):
+                return flux_dance_layer
+        raise ElementNotIncludedException(f"Element: {element} is not included in any flux-dance layer.")
+
+    def _get_flux_operation_identifier(self, element: IEdgeID) -> FluxOperationIdentifier:
+        """:return: Identifier describing flux-dance layer."""
+        for flux_dance_layer, flux_operation_identifier in self._flux_dances:
+            if flux_dance_layer.contains(element=element):
+                return flux_operation_identifier
+        raise ElementNotIncludedException(f"Element: {element} is not included in any flux-dance layer.")
+
+
+    def get_flux_operation_identifier(self, qubit_id0: str, qubit_id1: str) -> str:
+        """:return: Identifier describing flux-dance layer."""
+        edge: IEdgeID = EdgeIDObj(
+            qubit_id0=QubitIDObj(_id=qubit_id0),
+            qubit_id1=QubitIDObj(_id=qubit_id1),
+        )
+        return self._get_flux_operation_identifier(element=edge).id
+
+    def get_edge_flux_operation_identifier(self, ancilla_qubit: str) -> List[str]:
+        """:return: Identifier describing flux-dance layer."""
+        qubit_id: IQubitID = QubitIDObj(_id=ancilla_qubit)
+        parity_group: IParityGroup = self.get_parity_group(element=qubit_id)
+        return [
+            self._get_flux_operation_identifier(
+                element=edge_id,
+            ).id
+            for edge_id in parity_group.edge_ids
+        ]
+
     def _get_virtual_phase_identifier(self, directional_edge: DirectionalEdgeIDObj) -> VirtualPhaseIdentifier:
         """:return: Identifier for virtual phase correction. Based on element and parity group."""
         return self._virtual_phase_lookup[directional_edge]
@@ -492,6 +748,19 @@ def get_data_virtual_phase_identifiers(self, ancilla_qubit: str) -> List[str]:
             for data_id in parity_group.data_ids
         ]
 
+    def get_parity_data_identifier(self, ancilla_qubit: str) -> List[str]:
+        """
+        Iterates over provided ancilla qubit ID's.
+        Construct corresponding IQubitID's.
+        Obtain corresponding IParityGroup's.
+        Flatten list of (unique) data qubit ID's part of these parity groups.
+        :return: Array-like of (unique) data qubit ID's part of ancilla qubit parity groups.
+        """
+        ancilla_qubit_id: IQubitID = QubitIDObj(ancilla_qubit)
+        parity_group: IParityGroup = self.get_parity_group(element=ancilla_qubit_id)
+        data_qubit_ids: List[IQubitID] = [qubit_id for qubit_id in parity_group.data_ids]
+        return [qubit_id.id for qubit_id in data_qubit_ids]
+
     def get_parity_data_identifiers(self, ancilla_qubits: List[str]) -> List[str]:
         """
         Iterates over provided ancilla qubit ID's.
@@ -500,16 +769,15 @@ def get_parity_data_identifiers(self, ancilla_qubits: List[str]) -> List[str]:
         Flatten list of (unique) data qubit ID's part of these parity groups.
         :return: Array-like of (unique) data qubit ID's part of ancilla qubit parity groups.
         """
-        ancilla_qubit_ids: List[IQubitID] = [QubitIDObj(ancilla_qubit) for ancilla_qubit in ancilla_qubits]
-        parity_groups: List[IParityGroup] = [self.get_parity_group(element=qubit_id) for qubit_id in ancilla_qubit_ids]
-        data_qubit_ids: List[IQubitID] = [qubit_id for parity_group in parity_groups for qubit_id in parity_group.data_ids]
-        return [unique_qubit_id.id for unique_qubit_id in set(data_qubit_ids)]
+        return [unique_qubit_id for ancilla_qubit in ancilla_qubits for unique_qubit_id in set(self.get_parity_data_identifier(ancilla_qubit=ancilla_qubit))]
+    
+    def get_frequency_group_identifier(self, element: IQubitID) -> FrequencyGroupIdentifier:
+        """:return: Frequency group identifier based on qubit-ID."""
+        return Surface17Layer().get_frequency_group_identifier(element=element)
     # endregion
 
 
 if __name__ == '__main__':
 
-    for parity_group in Repetition9Layer().parity_group_x + Repetition9Layer().parity_group_z:
-        print(parity_group.ancilla_id.id)
-        print(f'(Ancilla) phase cw: {Repetition9Layer().get_ancilla_virtual_phase_identifier(parity_group.ancilla_id.id)}')
-        print(f'(Data) phase cw: ', [phase_id for phase_id in Repetition9Layer().get_data_virtual_phase_identifiers(parity_group.ancilla_id.id)])
+    flux_dance_0 = Repetition9Layer().get_flux_dance_at_round(0)
+    print(flux_dance_0.edge_ids)    
diff --git a/pycqed/qce_utils/control_interfaces/intrf_channel_identifier.py b/pycqed/qce_utils/control_interfaces/intrf_channel_identifier.py
index 24886c4e2..bdf6ffd94 100644
--- a/pycqed/qce_utils/control_interfaces/intrf_channel_identifier.py
+++ b/pycqed/qce_utils/control_interfaces/intrf_channel_identifier.py
@@ -64,6 +64,13 @@ class IEdgeID(IChannelIdentifier, metaclass=ABCMeta):
     Interface class, for qubit-to-qubit edge reference.
     """
 
+    # region Interface Properties
+    @property
+    def qubit_ids(self) -> List[IQubitID]:
+        """:return: All qubit-ID's."""
+        raise InterfaceMethodException
+    # endregion
+
     # region Interface Methods
     @abstractmethod
     def contains(self, element: IQubitID) -> bool:
@@ -246,6 +253,11 @@ class EdgeIDObj(IEdgeID):
     def id(self) -> QID:
         """:returns: Reference ID for edge."""
         return f"{self.qubit_id0.id}-{self.qubit_id1.id}"
+
+    @property
+    def qubit_ids(self) -> List[IQubitID]:
+        """:return: All qubit-ID's."""
+        return [self.qubit_id0, self.qubit_id1]
     # endregion
 
     # region Interface Methods