diff --git a/scripts/experiments/007_generate_fake_analysis_data.py b/scripts/experiments/007_generate_fake_analysis_data.py index c8f18d8..e4a6a97 100644 --- a/scripts/experiments/007_generate_fake_analysis_data.py +++ b/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, - N=2 * N + 2, + laser_detuning=laser_detuning, + 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 - eom_off_lifetime = total_lifetimes * 2 / 3 + total_lifetimes = 20 + 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 diff --git a/src/rabifun/analysis.py b/src/rabifun/analysis.py index bc5ef51..41f1aa6 100644 --- a/src/rabifun/analysis.py +++ b/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 diff --git a/src/rabifun/system.py b/src/rabifun/system.py index ca4b20a..7ec6918 100644 --- a/src/rabifun/system.py +++ b/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( - -1j * (laser_frequency(params, t) - runtime_params.detuned_Ωs.real) * t - ) + laser = np.exp(-1j * freqs * 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)