generalize to handle imperfect hybridization

scripts/experiments/007_generate_fake_analysis_data.py

@@ -10,7 +10,13 @@ import functools
 
 @functools.lru_cache()
 def make_params_and_solve(
-    total_lifetimes, eom_off_lifetime, ω_c=0.1 / 2, N=10, g_0=1 / 4
+    total_lifetimes,
+    eom_off_lifetime,
+    ω_c=0.1 / 2,
+    N=10,
+    g_0=1 / 4,
+    small_loop_detuning=0,
+    laser_detuning=0,
 ):
     """
     Make a set of parameters for the system with the current
@@ -25,8 +31,8 @@ def make_params_and_solve(
         δ=1 / 4,
         ω_c=ω_c,
         g_0=g_0,
-        laser_detuning=0,
+        laser_detuning=laser_detuning,
-        N=2 * N + 2,
+        N=N,
         N_couplings=N,
         measurement_detuning=Ω * (3),
         α=0,
@@ -34,6 +40,7 @@ def make_params_and_solve(
         flat_energies=False,
         correct_lamb_shift=0,
         laser_off_time=0,
+        small_loop_detuning=small_loop_detuning,
     )
 
     params.laser_off_time = params.lifetimes(eom_off_lifetime)
@@ -43,6 +50,8 @@ def make_params_and_solve(
         np.array([params.Ω, params.Ω * (1 - params.δ), params.δ * 2 * params.Ω]),
         np.repeat(params.Ω, 3),
     )
+    params.drive_override[1][-1] *= 4
+    params.drive_override[1][-1] *= 4
 
     t = time_axis(params, lifetimes=total_lifetimes, resolution=0.01)
     solution = solve(t, params)
@@ -56,25 +65,44 @@ def generate_phase_one_data():
     spectrum.
     """
 
-    total_lifetimes = 50
+    total_lifetimes = 20
-    eom_off_lifetime = total_lifetimes * 2 / 3
+    eom_off_lifetime = total_lifetimes * 1 / 2
     fluct_size = 0.05
+    SNR = 0.8
 
     params, t, solution = make_params_and_solve(
-        total_lifetimes, eom_off_lifetime, N=3, g_0=10
+        total_lifetimes,
+        eom_off_lifetime,
+        N=20,
+        g_0=1,
+        small_loop_detuning=0,  # 0.05,
+        laser_detuning=-0.01,
     )
     signal = output_signal(t, solution.y, params)
 
     rng = np.random.default_rng(seed=0)
-    signal += np.mean(abs(signal)) * rng.standard_normal(len(signal)) * fluct_size * 50
+    signal += (
+        np.sqrt(np.mean(abs(signal) ** 2)) / SNR * rng.standard_normal(len(signal))
+    )
 
     fig = make_figure()
-    ax_realtime, ax_spectrum = fig.subplots(2, 1)
 
-    ax_realtime.plot(t[::500], signal[::500])
+    ax_realtime, ax_rotating, ax_spectrum = fig.subplot_mosaic("""
+    AB
+    CC
+    """).values()
+
+    for mode in range(solution.y.shape[0]):
+        ax_rotating.plot(t[::50], np.abs(solution.y[mode, ::50]) ** 2)
+    ax_rotating.set_xlabel("Time")
+    ax_rotating.set_ylabel("Intensity")
+    ax_rotating.set_title("Mode Intensities in Rotating Frame")
+
+    ax_realtime.plot(t[::50], signal[::50])
     ax_realtime.axvline(params.laser_off_time, color="black", linestyle="--")
     ax_realtime.set_xlabel("Time")
     ax_realtime.set_ylabel("Intensity")
+    ax_realtime.set_title("Measures Intensity")
 
     # now we plot the power spectrum
     window = (float(params.laser_off_time or 0), float(t[-1]))
@@ -88,12 +116,12 @@ def generate_phase_one_data():
     )
 
     scan = ScanData(np.ones_like(signal), signal, t)
-    peak_info = find_peaks(scan, ringdown_params, window, prominence=0.1)
+    peak_info = find_peaks(scan, ringdown_params, window, prominence=0.1 / 4)
     peak_info = refine_peaks(peak_info, ringdown_params)
     plot_spectrum_and_peak_info(ax_spectrum, peak_info, ringdown_params)
 
     Ω, ΔΩ, δ, Δδ, ladder = extract_Ω_δ(
-        peak_info, ringdown_params, Ω_threshold=0.1, ladder_threshold=0.1
+        peak_info, ringdown_params, Ω_threshold=0.1, ladder_threshold=0.1, start_peaks=4
     )
 
     for index, type in ladder:
@@ -106,4 +134,11 @@ def generate_phase_one_data():
             color="C4" if type.value == StepType.BATH.value else "C5",
             label=type,
         )
+
+    fig.suptitle(
+        f"""Calibration Phase One Demonstration\n N={params.N} * 2 modes g_0={params.g_0}, SNR={SNR}, Ω (input) = {params.Ω:.2f}, δ (input) = {params.Ω*params.δ:.2f}
+Ω={Ω:.2f} ± {ΔΩ:.2f}, δ={δ:.2f} ± {Δδ:.2f}
+"""
+    )
+
     return Ω, ΔΩ, δ, Δδ, ladder

src/rabifun/analysis.py

@@ -314,7 +314,7 @@ def extract_Ω_δ(
     ΔΩ = max(np.sqrt(np.sum(Δcandidates**2)) / len(candidates), np.std(candidates))
 
     if np.isnan(Ω):
-        raise ValueError("No FSR found")
+        raise ValueError("No bath modes found!")
 
     # second step: we walk through the peaks and label them as for the
     total_peaks = len(peak_indices)
@@ -333,7 +333,7 @@ def extract_Ω_δ(
         current_peak,
         last_type=StepType.BATH,
         second_last_type=StepType.BATH,
-        bifurcations=3,
+        bifurcations=bifurcations,
     ):
         if current_peak == total_peaks - 1:
             return [], []
@@ -423,11 +423,15 @@ def extract_Ω_δ(
 
     costs = [cost / len(ladder) for cost, ladder in zip(costs, ladders)]
 
+    if len(costs) == 0:
+        raise ValueError("No valid ladders/spectra found.")
+
     best = np.argmin(costs)
     best_ladder = ladders[best]
 
     Ωs = []
     δs = []
+    Δδs = []
     Ω_m_δs = []
 
     for (i, (begin_index, begin_type)), (end_index, _) in zip(
@@ -438,17 +442,15 @@ def extract_Ω_δ(
                 Ωs.append(peak_freqs[end_index] - peak_freqs[begin_index])
             case StepType.BATH_TO_A:
                 Ω_m_δs.append(peak_freqs[end_index] - peak_freqs[begin_index])
-                if i < len(best_ladder) - 2:
-                    Ωs.append(
-                        (peak_freqs[best_ladder[i + 2][0]] - peak_freqs[begin_index])
-                        / 2
-                    )
             case StepType.A_TO_A:
                 δs.append((peak_freqs[end_index] - peak_freqs[begin_index]) / 2)
+                Δδs.append(
+                    np.sqrt(Δpeak_freqs[end_index] ** 2 + Δpeak_freqs[begin_index] ** 2)
+                )
 
     Ω = np.mean(Ωs)
     ΔΩ = np.std(Ωs)
     δ = np.mean(δs)
-    Δδ = np.std(δs)
+    Δδ = np.mean(Δδs)
 
     return Ω, ΔΩ, δ, Δδ, best_ladder

src/rabifun/system.py

@@ -78,6 +78,9 @@ class Params:
     The drive strength is normalized to :any:`g_0`.
     """
 
+    small_loop_detuning: float = 0
+    """The detuning (in units of :any:`Ω`) of the small loop mode relative to the ``A`` mode."""
+
     def __post_init__(self):
         if self.N_couplings > self.N:
             raise ValueError("N_couplings must be less than or equal to N.")
@@ -124,12 +127,28 @@ class RuntimeParams:
     """Secondary Parameters that are required to run the simulation."""
 
     def __init__(self, params: Params):
+        H_A = np.array(
+            [
+                [0, params.δ],
+                [params.δ, params.small_loop_detuning],
+            ]
+        )
+
+        eig = np.linalg.eigh(H_A)
+        idx = np.argsort(eig.eigenvalues)
+        anti_a_frequency, a_frequency = eig.eigenvalues[idx]
+        self.a_weights = np.abs(eig.eigenvectors[:, idx][0, :])
+
         bath = np.arange(1, params.N + 1)
         freqs = (
             2
             * np.pi
             * params.Ω
-            * np.concatenate([[-1 * params.δ, params.δ], bath, -bath])
+            * (
+                np.concatenate([[anti_a_frequency, a_frequency], bath, -bath])
+                - a_frequency
+                + params.δ
+            )
         )
 
         decay_rates = -1j * np.repeat(params.η / 2, 2 * params.N + 2)
@@ -138,16 +157,28 @@ class RuntimeParams:
         self.drive_frequencies, self.detunings, self.g, a_shift = (
             drive_frequencies_and_amplitudes(params)
         )  # linear frequencies!
-
         if params.drive_override is not None:
             self.drive_frequencies = params.drive_override[0]
             self.g = 2 * np.pi * params.drive_override[1]
             a_shift = 0
             self.detunings *= 0
-            self.g /= params.g_0 * np.sqrt(np.sum(self.g**2))
+            self.g *= params.g_0 / np.sqrt(np.sum(self.g**2))
+
+        # print(params.δ * 2 - (a_frequency - anti_a_frequency))
+        # print((params.Ω - params.δ) - (params.Ω - a_frequency))
+        # print(
+        #     (freqs[2] - freqs[1]) / (2 * np.pi),
+        #     (params.Ω * (1 - a_frequency)),
+        #     (params.Ω * (1 - params.δ)),
+        # )
+        # print(np.diff(freqs) / (2 * np.pi))
+        # import ipdb
+
+        # ipdb.set_trace()
 
         self.g *= 2 * np.pi
         self.Ωs = Ωs
+
         self.ε = (
             2
             * np.pi
@@ -206,13 +237,15 @@ def time_axis(
     return np.arange(0, tmax, resolution * np.pi / (params.Ω * params.N))
 
 
-def eom_drive(t, x, ds, ωs, det_matrix):
+def eom_drive(t, x, ds, ω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 det_matrix: detuning matrix
+    :param a_weights: weights of the A modes
     """
 
     # test = abs(det_matrix.copy())
@@ -221,6 +254,14 @@ def eom_drive(t, x, ds, ωs, det_matrix):
     # print(np.argmin(test, keepdims=True))
 
     det_matrix = np.exp(-1j * det_matrix * t)
+    for i, weight in enumerate(a_weights):
+        det_matrix[i, 2:] *= weight
+        det_matrix[2:, i] *= weight.conjugate()
+
+    # FIXME: that's not strictly right for the non symmetric damping
+    prod = a_weights[0] * a_weights[1].conj()
+    det_matrix[0, 1] *= prod
+    det_matrix[1, 0] *= prod.conjugate()
 
     driven_x = np.sum(2 * ds * np.sin(2 * np.pi * ωs * t)) * (det_matrix @ x)
 
@@ -261,21 +302,20 @@ def make_righthand_side(runtime_params: RuntimeParams, params: Params):
                     runtime_params.g,
                     runtime_params.drive_frequencies,
                     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
 
-            laser = np.exp(
+            laser = np.exp(-1j * freqs * t)
-                -1j * (laser_frequency(params, t) - runtime_params.detuned_Ωs.real) * t
-            )
 
             if params.rwa:
                 index = np.argmin(abs(freqs))
                 laser[0:index] = 0
                 laser[index + 1 :] = 0
 
-            differential[0:2] += laser[:2] / np.sqrt(2)
+            differential[0:2] += laser[:2] * runtime_params.a_weights
             differential[2:] += laser[2:]
 
         if params.rwa:
@@ -323,23 +363,6 @@ def solve(t: np.ndarray, params: Params, **kwargs):
     )
 
 
-def in_rotating_frame(
-    t: np.ndarray, amplitudes: np.ndarray, params: Params
-) -> np.ndarray:
-    """Transform the amplitudes to the rotating frame."""
-    Ωs = RuntimeParams(params).Ωs
-
-    detunings = np.concatenate(
-        [[0, 0], drive_detunings(params), np.zeros(params.N - params.N_couplings)]
-    )
-
-    return amplitudes * np.exp(
-        1j
-        * (Ωs[:, None].real - detunings[:, None] + laser_frequency(params, t)[None, :])
-        * t[None, :]
-    )
-
-
 def output_signal(t: np.ndarray, amplitudes: np.ndarray, params: Params):
     """
     Calculate the output signal when mixing with laser light of
@@ -347,7 +370,8 @@ def output_signal(t: np.ndarray, amplitudes: np.ndarray, params: Params):
     """
 
     runtime = RuntimeParams(params)
-    rotating = amplitudes * np.exp(-1j * runtime.detuned_Ωs[:, None] * t)
+    rotating = amplitudes * np.exp(-1j * runtime.detuned_Ωs[:, None] * t[None, :])
+    rotating[0:2, :] *= runtime.a_weights[:, None].conjugate()
 
     return (
         np.sum(rotating, axis=0)