implement drive phase

scripts/experiments/013_step_through.py

@@ -74,9 +74,7 @@ def solve_shot(t, params, t_before, t_after):
     return t_after, amps
 
 
-def make_shots(
+def make_shots(params, total_lifetimes, eom_range, eom_steps, num_freq):
-    params, total_lifetimes, eom_range, eom_steps, σ_modulation_time, num_freq
-):
     solutions = []
 
     analyze_time = params.lifetimes(total_lifetimes) - params.laser_off_time
@@ -91,15 +89,17 @@ def make_shots(
         base = params.Ω * (
             eom_range[0] + (eom_range[1] - eom_range[0]) * step / eom_steps
         )
+
         current_params = copy.deepcopy(params)
         current_params.drive_override = (
             base + params.Ω * np.arange(num_freq),
             np.ones(num_freq),
         )
+        current_params.drive_phases = rng.uniform(0, 2 * np.pi, size=num_freq)
 
-        off_time = rng.normal(
+        off_time = rng.normal(params.laser_off_time, 0.1 * params.laser_off_time)
-            params.laser_off_time, σ_modulation_time * params.lifetimes(1)
+
-        )
+        params.laser_off_time
         current_params.laser_off_time = off_time
         current_params.drive_off_time = off_time
         current_params.total_lifetimes = (off_time + analyze_time) / params.lifetimes(1)
@@ -114,8 +114,6 @@ def make_shots(
 
 
 def process_shots(solutions, noise_amplitude, params):
-    average_power_spectrum = None
-
     rng = np.random.default_rng(seed=0)
 
     # let us get a measure calibrate the noise strength
@@ -167,23 +165,26 @@ def plot_power_spectrum(
     )
 
     peak_info = find_peaks(
-        freq, average_power_spectrum, ringdown_params, prominence=0.1, height=0.5
+        freq, average_power_spectrum, ringdown_params, prominence=0.05, height=0.1
     )
 
     peak_info, lm_result = refine_peaks(
         peak_info, ringdown_params, height_cutoff=0.05, σ=σ_power_spectrum
     )
-    print(lm_result.fit_report())
+
     peak_info.power = average_power_spectrum
     plot_spectrum_and_peak_info(
         ax_spectrum, peak_info, ringdown_params, annotate=annotate
     )
     if lm_result is not None:
+        # print(lm_result.fit_report())
         fine_freq = np.linspace(freq.min(), freq.max(), 5000)
         fine_fit = lm_result.eval(ω=fine_freq)
         ax_spectrum.plot(fine_freq, fine_fit, color="red")
         ax_spectrum.set_ylim(-0.1, max(1, fine_fit.max() * 1.1))
 
+    print(runtime.Ωs.real / (2 * np.pi))
+
     for i, peak_freq in enumerate(runtime.Ωs):
         pos = np.abs(
             params.measurement_detuning
@@ -221,7 +222,6 @@ def generate_data(
     small_loop_detuning=0,
     excitation_lifetimes=2,
     measurement_lifetimes=4,
-    σ_modulation_time=0.01,
     num_freq=3,
 ):
     η = 0.2
@@ -234,7 +234,7 @@ def generate_data(
         δ=1 / 4,
         ω_c=0.1,
         g_0=g_0,
-        laser_detuning=laser_detuning,  # 13 * (-1 - 1 / 4) + laser_detuning,
+        laser_detuning=laser_detuning,
         N=N,
         N_couplings=N,
         measurement_detuning=0,
@@ -244,7 +244,7 @@ def generate_data(
         correct_lamb_shift=0,
         laser_off_time=0,
         small_loop_detuning=small_loop_detuning,
-        drive_override=(np.array([1.0]), np.array([1.0])),
+        drive_override=(np.array([]), np.array([])),
     )
 
     params.laser_off_time = params.lifetimes(excitation_lifetimes)
@@ -255,25 +255,22 @@ def generate_data(
         excitation_lifetimes + measurement_lifetimes,
         eom_ranges,
         eom_steps,
-        σ_modulation_time,
         num_freq,
     )
 
     (sol_on_res) = make_shots(
         params,
         excitation_lifetimes + measurement_lifetimes,
-        ((1 - params.δ), (1 - params.δ)),
+        ((1 + params.δ), (1 + params.δ)),
         1,
-        0,
         num_freq,
     )
 
     (sol_on_res_bath) = make_shots(
         params,
         excitation_lifetimes + measurement_lifetimes,
-        ((1 - params.δ * 1.1), (1 - params.δ * 1.1)),
+        ((1 + params.δ * 1.1), (1 + params.δ * 1.1)),
         1,
-        0,
         num_freq,
     )
 
@@ -295,7 +292,7 @@ def generate_data(
 
     Ω/2π = {params.Ω}MHz, η/2π = {params.η}MHz, g_0 = {params.g_0}Ω, N = {params.N}
     noise amplitude = {noise_amplitude} * 2/η, η_A = {η_factor} x η, EOM stepped from {eom_ranges[0]:.2f}Ω to {eom_ranges[1]:.2f}Ω in {eom_steps} steps
-    total time = {(excitation_lifetimes + measurement_lifetimes) * eom_steps} / η
+    total time = {(excitation_lifetimes + measurement_lifetimes) * eom_steps / (params.η * 1e6)}s
     """)
     ax_multi, ax_single, ax_single_bath = fig.subplot_mosaic("AA\nBC").values()
 
@@ -347,16 +344,15 @@ def generate_data(
 # %% save
 if __name__ == "__main__":
     fig = generate_data(
-        g_0=0.6,
+        g_0=0.5,
         η_factor=5,
-        noise_amplitude=8e-3,
+        noise_amplitude=5e-3,
-        N=4,
+        N=5,
-        eom_ranges=(0.7, 0.9),
+        eom_ranges=(0.7, 0.9),  # (1.9, 2.1),
         eom_steps=100,
         small_loop_detuning=0,
         laser_detuning=0,
         excitation_lifetimes=1,
-        measurement_lifetimes=4,
+        measurement_lifetimes=3,
-        σ_modulation_time=0.2,
         num_freq=4,
     )

src/rabifun/analysis.py

@@ -218,7 +218,7 @@ def filter_peaks(
             i in to_be_deleted
             or Δω0 > uncertainty_threshold * params.fΩ_guess
             or A < height_cutoff
-            or A > 10
+            or A > 5
             or Δγ > uncertainty_threshold * params.fΩ_guess
         ):
             np.delete(peaks.peaks, i)
@@ -255,6 +255,9 @@ def refine_peaks(
         :any:`ringdown_params.fΩ_guess`.
     """
 
+    if len(peaks.peaks) == 0:
+        return peaks, None
+
     peaks = dataclasses.replace(peaks)
     freqs = peaks.freq
     peak_freqs = peaks.peak_freqs
@@ -288,7 +291,6 @@ def refine_peaks(
 
     initial_params = model.make_params(
         A=dict(value=1, min=0, max=np.inf),
-        ω0=dict(value=0, min=0, max=np.inf),
         γ=dict(value=params.η_guess, min=0, max=np.inf),
     )
 

src/rabifun/system.py

@@ -84,6 +84,9 @@ class Params:
     The drive strength is normalized to :any:`g_0`.
     """
 
+    drive_phases: np.ndarray | None = None
+    """The phase for each drive tone. If not specified, the phases are zero."""
+
     small_loop_detuning: float = 0
     """The detuning (in units of :any:`Ω`) of the small loop mode relative to the ``A`` mode."""
 
@@ -106,6 +109,12 @@ class Params:
             if self.rwa:
                 raise ValueError("Drive override is not compatible with the RWA.")
 
+        if self.drive_phases is not None:
+            if len(self.drive_phases) != self.N_couplings:
+                raise ValueError("Need as many drive phases as couplings.")
+        else:
+            self.drive_phases = np.zeros(self.N_couplings)
+
         if self.η_hybrid is None:
             self.η_hybrid = self.η
 
@@ -265,13 +274,14 @@ def time_axis(
     return np.arange(0, tmax, resolution * np.pi / (params.Ω * params.N))
 
 
-def eom_drive(t, x, ds, ωs, det_matrix, a_weights):
+def eom_drive(t, x, ds, ωs, φs, det_matrix, a_weights):
     """The electrooptical modulation drive.
 
     :param t: time
     :param x: amplitudes
     :param ds: drive amplitudes
     :param ωs: linear drive frequencies
+    :param φs: drive phases
     :param det_matrix: detuning matrix
     :param a_weights: weights of the A modes
     """
@@ -291,7 +301,7 @@ def eom_drive(t, x, ds, ωs, det_matrix, a_weights):
     rot_matrix[0, 1] *= prod
     rot_matrix[1, 0] *= prod.conjugate()
 
-    driven_x = np.sum(2 * ds * np.sin(2 * np.pi * ωs * t)) * (rot_matrix @ x)
+    driven_x = np.sum(2 * ds * np.sin(2 * np.pi * ωs * t + φs)) * (rot_matrix @ x)
 
     return driven_x
 
@@ -329,9 +339,11 @@ def make_righthand_side(runtime_params: RuntimeParams, params: Params):
                     x,
                     runtime_params.g,
                     runtime_params.drive_frequencies,
+                    params.drive_phases,
                     runtime_params.detuning_matrix,
                     runtime_params.a_weights,
                 )
+
         if (params.laser_off_time is None) or (t < params.laser_off_time):
             freqs = laser_frequency(params, t) - runtime_params.detuned_Ωs.real