From 8b4d0ec1bf55240b870ccfc59671ef07ff021b30 Mon Sep 17 00:00:00 2001 From: Valentin Boettcher Date: Tue, 16 Jul 2024 15:26:11 -0400 Subject: [PATCH] try to tease something out of the 11_07 data --- .../007_generate_fake_analysis_data.py | 241 ++++++++++++------ .../008_analyze_11_07_transients.py | 84 +++--- src/rabifun/analysis.py | 63 +++-- src/rabifun/plots.py | 27 +- src/rabifun/system.py | 46 ++-- src/ringfit/plotting.py | 5 +- 6 files changed, 296 insertions(+), 170 deletions(-) diff --git a/scripts/experiments/007_generate_fake_analysis_data.py b/scripts/experiments/007_generate_fake_analysis_data.py index 0f189be..29acfb3 100644 --- a/scripts/experiments/007_generate_fake_analysis_data.py +++ b/scripts/experiments/007_generate_fake_analysis_data.py @@ -8,7 +8,6 @@ import functools # %% interactive -@functools.lru_cache() def make_params_and_solve( total_lifetimes, eom_off_lifetime, @@ -17,6 +16,9 @@ def make_params_and_solve( g_0=1 / 4, small_loop_detuning=0, laser_detuning=0, + η=1.0, + η_hybrid=1.0, + drive_mode="all", ): """ Make a set of parameters for the system with the current @@ -26,7 +28,8 @@ def make_params_and_solve( Ω = 13 params = Params( - η=0.5, + η=η, + η_hybrid=η_hybrid, Ω=Ω, δ=1 / 4, ω_c=ω_c, @@ -34,7 +37,7 @@ def make_params_and_solve( laser_detuning=laser_detuning, N=N, N_couplings=N, - measurement_detuning=20, + measurement_detuning=0, α=0, rwa=False, flat_energies=False, @@ -46,13 +49,34 @@ def make_params_and_solve( params.laser_off_time = params.lifetimes(eom_off_lifetime) params.drive_off_time = params.lifetimes(eom_off_lifetime) - params.drive_override = ( - np.array([params.Ω, params.Ω * (1 - params.δ), params.δ * 2 * params.Ω]), - np.ones(3), - ) - params.drive_override[1][-1] *= 2 + match drive_mode: + case "hybrid": + params.drive_override = ( + np.array([params.Ω * params.δ * 2]), + np.ones(1), + ) - t = time_axis(params, lifetimes=total_lifetimes, resolution=0.01) + case "hybrid_to_one_bath": + params.drive_override = ( + np.array([params.Ω * (1 - params.δ), params.Ω * params.δ * 2]), + np.ones(2), + ) + + case "bath": + params.drive_override = ( + np.array([params.Ω]), + np.ones(1), + ) + + case _: + params.drive_override = ( + np.array( + [params.Ω, params.Ω * (1 - params.δ), params.δ * 2 * params.Ω] + ), + np.ones(3), + ) + + t = time_axis(params, lifetimes=total_lifetimes, resolution=0.001) solution = solve(t, params) return params, t, solution @@ -64,86 +88,143 @@ def generate_phase_one_data(): spectrum. """ - total_lifetimes = 20 - eom_off_lifetime = total_lifetimes * 1 / 2 + total_lifetimes = 30 + eom_off_lifetime = total_lifetimes * 0.210 fluct_size = 0.05 - SNR = 0.5 + noise_amp = 0.1 + η = 0.2 params, t, solution = make_params_and_solve( total_lifetimes, eom_off_lifetime, - N=30, - g_0=1, - small_loop_detuning=0.1, - laser_detuning=0.1, - ) - signal = output_signal(t, solution.y, params) - - rng = np.random.default_rng(seed=0) - signal += ( - np.sqrt(np.mean(abs(signal) ** 2)) / SNR * rng.standard_normal(len(signal)) + N=10, + g_0=0.7, + small_loop_detuning=0, + laser_detuning=0, + η=η, + η_hybrid=η * 3.5, + drive_mode="full", ) - fig = make_figure() + raw_signal = output_signal(t, solution.y, params) + for noise in [True, False]: + rng = np.random.default_rng(seed=0) + signal = raw_signal.copy() + if noise: + signal += noise_amp * rng.standard_normal(len(signal)) - ax_realtime, ax_rotating, ax_spectrum = fig.subplot_mosaic(""" - AB - CC - """).values() + fig = make_figure(f"simulation_noise_{noise}") - 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])) - - ringdown_params = RingdownParams( - fω_shift=params.measurement_detuning, - mode_window=(5, 5), - fΩ_guess=params.Ω * (1 + rng.standard_normal() * fluct_size), - fδ_guess=params.Ω * params.δ * (1 + rng.standard_normal() * fluct_size), - η_guess=params.η * (1 + rng.standard_normal() * fluct_size), - ) - - scan = ScanData(np.ones_like(signal), signal, t) - 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, - start_peaks=6, - bifurcations=5, - ) - - for index, type in ladder: - freq_index = peak_info.peaks[index] - print(type, type is StepType.BATH) - ax_spectrum.plot( - peak_info.freq[freq_index], - peak_info.normalized_power[freq_index], - "o" if type.value == StepType.BATH.value else "*", - color="C4" if type.value == StepType.BATH.value else "C5", - label=type, + ax_realtime, ax_rotating_bath, ax_rotating_system, ax_spectrum = ( + fig.subplot_mosaic(""" + AA + BC + DD + DD + """).values() ) - runtime = RuntimeParams(params) - fig.suptitle( - f"""Calibration Phase One Demonstration\n N={params.N} * 2 modes g_0={params.g_0}, SNR={SNR}, Ω (input) = {params.Ω:.2f}, δ (input) = {(runtime.mode_splitting):.2f} -Ω={Ω:.2f} ± {ΔΩ:.2f}, δ={δ:.2f} ± {Δδ:.2f} -""" - ) + ax_rotating_system.sharey(ax_rotating_bath) + ax_rotating_system.sharex(ax_rotating_bath) - return Ω, ΔΩ, δ, Δδ, ladder + for mode in range(2): + ax_rotating_system.plot( + t[::50], + abs(np.imag(solution.y[mode, ::50])) ** 2 + / 2 + * (params.η / params.η_hybrid) ** 2, + ) + + for mode in range(2, solution.y.shape[0]): + ax_rotating_bath.plot(t[::50], abs(np.imag(solution.y[mode, ::50])) ** 2) + + for ax in [ax_rotating_bath, ax_rotating_system]: + ax.set_xlabel("Time") + ax.set_ylabel("Intensity") + + ax_rotating_bath.set_title("Bath Modes") + ax_rotating_system.set_title( + "Hybridized Modes * 1/sqrt(2) * correction for decay rate" + ) + + ax_realtime.plot(t[::50], signal[::50]) + ax_realtime.axvline(params.drive_off_time, color="black", linestyle="--") + ax_realtime.set_xlabel("Time") + ax_realtime.set_ylabel("Intensity") + ax_realtime.set_title("Photo-diode AC Intensity") + + # now we plot the power spectrum + window = (float(params.laser_off_time or 0), t[-1]) + # window = (0, float(params.laser_off_time or 0)) + + ringdown_params = RingdownParams( + fω_shift=params.measurement_detuning, + mode_window=(5, 5), + fΩ_guess=params.Ω * (1 + rng.standard_normal() * fluct_size), + fδ_guess=params.Ω * params.δ * (1 + rng.standard_normal() * fluct_size), + η_guess=0.5, # params.η * (1 + rng.standard_normal() * fluct_size), + absolute_low_cutoff=2, + ) + + scan = ScanData(np.ones_like(signal), signal, t) + peak_info = find_peaks(scan, ringdown_params, window, prominence=0.01) + peak_info = refine_peaks(peak_info, ringdown_params) + plot_spectrum_and_peak_info(ax_spectrum, peak_info, ringdown_params) + + runtime = RuntimeParams(params) + for i, freq in enumerate(runtime.Ωs.real): + pos = np.abs( + params.measurement_detuning + - freq / (2 * np.pi) + + params.δ * params.Ω + + params.laser_detuning, + ) + + ax_spectrum.axvline( + pos, + color="red", + linestyle="--", + zorder=-100, + ) + ax_spectrum.annotate(mode_name(i), (pos + 0.1, peak_info.power.max())) + + ax_spectrum.set_xlim(ringdown_params.low_cutoff, ringdown_params.high_cutoff) + + try: + Ω, ΔΩ, δ, Δδ, ladder = extract_Ω_δ( + peak_info, + ringdown_params, + Ω_threshold=0.1, + ladder_threshold=0.1, + start_peaks=6, + bifurcations=3, + ) + + if ladder: + for index, type in ladder: + freq_index = peak_info.peaks[index] + print(type, type is StepType.BATH) + ax_spectrum.plot( + peak_info.freq[freq_index], + peak_info.power[freq_index], + "o" if type.value == StepType.BATH.value else "*", + color="C4" if type.value == StepType.BATH.value else "C5", + label=type, + ) + + else: + δ = Δδ = np.nan + + except Exception: + Ω = ΔΩ = δ = Δδ = np.nan + + fig.suptitle( + f"""Calibration Phase One Demonstration {'with noise' if noise else ''} +N={params.N} * 2 + 2 modes g_0={params.g_0}Ω, Noise Amp={noise_amp}, η={params.η}MHz, η_hybrid={params.η_hybrid}MHz + + +Measured: +Ω (input) = {params.Ω:.2f}, δ (input) = {(runtime.mode_splitting):.2f} +Ω={Ω:.2f} ± {ΔΩ:.2f}, δ={δ:.2f} ± {Δδ:.2f} + """ + ) diff --git a/scripts/experiments/008_analyze_11_07_transients.py b/scripts/experiments/008_analyze_11_07_transients.py index f4f4ce7..0e5d517 100644 --- a/scripts/experiments/008_analyze_11_07_transients.py +++ b/scripts/experiments/008_analyze_11_07_transients.py @@ -4,10 +4,10 @@ from rabifun.utilities import * from rabifun.analysis import * from ringfit.data import ScanData from ringfit.plotting import * +from scipy.signal import butter, sosfilt import functools import gc -# %% load data path = "../../data/11_07_24/second_signal" scan = ScanData.from_dir(path, extension="npz") @@ -16,75 +16,67 @@ scan = ScanData.from_dir(path, extension="npz") gc.collect() fig = plt.figure("interactive", constrained_layout=True, figsize=(20, 3 * 5)) fig.clf() -(ax, ax2, ax_signal, ax_window, ax_spectrum) = fig.subplot_mosaic("AB;CC;DE").values() -plot_scan(scan, ax=ax_signal, smoothe_output=1e-8, linewidth=0.5) +(ax_signal, ax_window, ax_spectrum) = fig.subplot_mosaic("AA;BC").values() +plot_scan(scan, ax=ax_signal, linewidth=0.5, every=1000) + +# %% filter +# high_pass_filter = butter( +# 30, 3e6, btype="highpass", output="sos", fs=1 / scan.timestep, analog=False +# ) +# scan._output = sosfilt(high_pass_filter, scan.output) # %% window T_step = 0.0002 -t_peak = 0.032201 -N = 100 - 3 -t_scan_peak = t_peak - T_step * N + 4.07e-4 # T * N_steps +N = 100 +t_scan_peak = 0.0057815 + 0.1e-6 # T * N_steps +t_scan_peak = 0.0057815 - 0.09e-6 + 11 * T_step # T * N_steps +t_peak = t_scan_peak + N * T_step win_length = 5e-05 -window = t_scan_peak - win_length / 2, t_scan_peak + win_length / 2 * 0.8 +window = t_scan_peak, t_scan_peak + win_length * 0.3 ax_signal.axvline(t_peak, color="r", linestyle="--") ax_signal.axvline(t_scan_peak, color="r", linestyle="--") -ax_signal.axvspan( - *window, - color="r", - linestyle="--", -) +ax_signal.axvspan(*window, color="r", linestyle="--") mask = (scan.time > window[0]) & (scan.time < window[1]) ax_window.clear() + + ax_window.plot(scan.time[mask], scan.output[mask], linewidth=0.1) -freq, fft = fourier_transform( - scan.time, scan.output, window=window, low_cutoff=2e6, high_cutoff=90e6 -) -ax.clear() -ax.plot(freq, np.abs(fft) ** 2) ringdown_params = RingdownParams( fω_shift=0, - mode_window=(0, 4), - fΩ_guess=13e6, - fδ_guess=0.2 * 13e6, - η_guess=0.2e6, + mode_window=(0, 5), + fΩ_guess=12.9e6, + fδ_guess=0.2 * 12.9e6, + η_guess=0.5e6, absolute_low_cutoff=2e6, ) - -peak_info = find_peaks(scan, ringdown_params, window, prominence=0.01) +peak_info = find_peaks(scan, ringdown_params, window, prominence=0.08) peak_info = refine_peaks(peak_info, ringdown_params) -plot_spectrum_and_peak_info(ax_spectrum, peak_info, ringdown_params) +plot_spectrum_and_peak_info(ax_spectrum, peak_info, ringdown_params, annotate=True) print(peak_info.peak_freqs * 1e-6) +b1, b2, b3, reflected = peak_info.peak_freqs[[0, 1, 3, 2]] +Ω = b3 - b2 +δ = abs(b3 - 3 * Ω) +δ_2 = reflected - 2 * Ω +detuning = (δ + δ_2) / 2 -Ω, ΔΩ, δ, Δδ, ladder = extract_Ω_δ( - peak_info, - ringdown_params, - Ω_threshold=0.1, - ladder_threshold=0.1, - start_peaks=3, - bifurcations=2, +hyb_amp = peak_info.power[peak_info.peaks[0]] / (np.sqrt(2) * 3.5) ** 2 +hyb_width = peak_info.peak_widths[0] * 3.5 +hyb_freq = 2 * detuning + +ax_spectrum.plot( + peak_info.freq, + lorentzian(peak_info.freq, hyb_amp, hyb_freq, hyb_width), + color="C3", + label="hybridized mode?", ) - -for index, type in ladder: - freq_index = peak_info.peaks[index] - print( - type, - index, - (peak_info.peak_freqs[index] - peak_info.peak_freqs[ladder[0][0]]) / Ω, - ) - ax_spectrum.plot( - peak_info.freq[freq_index], - peak_info.normalized_power[freq_index], - "o" if type.value == StepType.BATH.value else "*", - color="C4" if type.value == StepType.BATH.value else "C5", - label=type, - ) +ax_spectrum.legend() diff --git a/src/rabifun/analysis.py b/src/rabifun/analysis.py index d407e0f..767cc74 100644 --- a/src/rabifun/analysis.py +++ b/src/rabifun/analysis.py @@ -38,21 +38,17 @@ def fourier_transform( return (freq[mask], fft[mask], t) if ret_time else (freq[mask], fft[mask]) -def lorentzian(ω, A, ω0, γ, offset): +def lorentzian(ω, A, ω0, γ): """A Lorentzian function with amplitude ``A``, center frequency ``ω0``, and decay rate ``γ`` and offset ``offset``. :param ω: Frequency array. :param A: Amplitude. :param ω0: Center frequency. - :param γ: Decay rate. The decay of an amplitude is :math:`\frac{γ}{2}`. - :param offset: Vertical offset. + :param γ: Decay rate. The decay of a squared amplitude is :math:`2πγ`. """ - return ( - A * (γ / (4 * np.pi)) ** 2 * (1 / ((ω - ω0) ** 2 + (γ / (4 * np.pi)) ** 2)) - + offset - ) + return A * (γ / 2) ** 2 * (1 / ((ω - ω0) ** 2 + (γ / 2) ** 2)) def complex_lorentzian(ω, A, ω0, γ): @@ -62,10 +58,10 @@ def complex_lorentzian(ω, A, ω0, γ): :param ω: Frequency array. :param A: Amplitude. :param ω0: Center frequency. - :param γ: Decay rate. The decay of an amplitude is :math:`\frac{γ}{2}`. + :param γ: Decay rate. The decay of a squared amplitude is :math:`2πγ`. """ - return A * (γ / 2) * 1 / (-1j * (ω - ω0) - (γ / (4 * np.pi))) + return A * (γ / 2) * 1 / (-1j * (ω - ω0) - (γ / (2))) ############################################################################### @@ -118,8 +114,8 @@ class RingdownPeakData: fft: np.ndarray """The fft amplitudes.""" - normalized_power: np.ndarray - """The normalized power spectrum of the fft.""" + power: np.ndarray + """The power spectrum of the fft.""" peaks: np.ndarray """The indices of the peaks.""" @@ -141,6 +137,9 @@ class RingdownPeakData: Δpeak_widths: np.ndarray | None = None """The uncertainty in the peak widths.""" + lorentz_params: list | None = None + """The lorentzian fit params to be fed into :any:`lorentzian`.""" + @property def is_refined(self) -> bool: """Whether the peaks have been refined with :any:`refine_peaks`.""" @@ -153,7 +152,7 @@ def find_peaks( window: tuple[float, float], prominence: float = 0.005, ) -> RingdownPeakData: - """Determine the peaks of the normalized power spectrum of the + """Determine the peaks of the power spectrum of the ringdown data. :param data: The oscilloscope data. @@ -174,11 +173,13 @@ def find_peaks( freq_step = freq[1] - freq[0] power = np.abs(fft) ** 2 - power /= power.max() distance = params.fδ_guess / 2 / freq_step peaks, peak_info = scipy.signal.find_peaks( - power, distance=distance, wlen=distance // 4, prominence=prominence + power / power.max(), + distance=distance, + wlen=distance // 4, + prominence=prominence, ) peak_freqs = freq[peaks] @@ -189,7 +190,7 @@ def find_peaks( peaks=peaks, peak_freqs=peak_freqs, peak_info=peak_info, - normalized_power=power, + power=power, ) @@ -209,12 +210,13 @@ def refine_peaks( peaks = dataclasses.replace(peaks) freqs = peaks.freq peak_freqs = peaks.peak_freqs - power = peaks.normalized_power + power = peaks.power new_freqs = [] new_widths = [] Δfreqs = [] Δwidths = [] + lorentz_params = [] window = params.η_guess * 3 deleted_peaks = [] @@ -223,10 +225,23 @@ def refine_peaks( windowed_freqs = freqs[mask] windowed_power = power[mask] - p0 = [1, peak_freq, params.η_guess, 0] + scale_freqs = windowed_freqs[-1] - windowed_freqs[0] + root_freq = windowed_freqs[0] + scale_power = windowed_power.max() + + windowed_freqs -= root_freq + windowed_freqs /= scale_freqs + + windowed_power /= scale_power + + p0 = [ + peaks.power[peaks.peaks[i]] / scale_power, + (peak_freq - root_freq) / scale_freqs, + params.η_guess / scale_freqs, + ] bounds = ( - [0, windowed_freqs[0], 0, 0], - [np.inf, windowed_freqs[-1], np.inf, 1], + [0, windowed_freqs[0], 0], + [np.inf, windowed_freqs[-1], np.inf], ) try: @@ -238,6 +253,11 @@ def refine_peaks( bounds=bounds, ) perr = np.sqrt(np.diag(pcov)) + popt[0] = popt[0] * scale_power + popt[1] = popt[1] * scale_freqs + root_freq + popt[2] = popt[2] * scale_freqs + + lorentz_params.append(popt) if perr[1] > uncertainty_threshold * params.fΩ_guess: deleted_peaks.append(i) @@ -260,6 +280,7 @@ def refine_peaks( peaks.Δpeak_freqs = np.array(Δfreqs) peaks.peak_widths = np.array(new_widths) peaks.Δpeak_widths = np.array(Δwidths) + peaks.lorentz_params = lorentz_params return peaks @@ -450,7 +471,9 @@ def extract_Ω_δ( costs = [cost / len(ladder) for cost, ladder in zip(costs, ladders)] if len(costs) == 0: - raise ValueError("No valid ladders/spectra found.") + print("No valid ladders/spectra found.") + + return Ω, ΔΩ, None, None, None best = np.argmin(costs) best_ladder = ladders[best] diff --git a/src/rabifun/plots.py b/src/rabifun/plots.py index 205f2ee..b873831 100644 --- a/src/rabifun/plots.py +++ b/src/rabifun/plots.py @@ -9,7 +9,7 @@ from .system import ( uncoupled_mode_indices, correct_for_decay, ) -from .analysis import fourier_transform, RingdownPeakData, RingdownParams +from .analysis import fourier_transform, RingdownPeakData, RingdownParams, lorentzian import matplotlib.pyplot as plt import numpy as np @@ -213,27 +213,43 @@ def plot_rwa_vs_real_amplitudes(ax, solution_no_rwa, solution_rwa, params, **kwa return no_rwa_lines, rwa_lines -def plot_spectrum_and_peak_info(ax, peaks: RingdownPeakData, params: RingdownParams): +def plot_spectrum_and_peak_info( + ax, peaks: RingdownPeakData, params: RingdownParams, annotate=False +): """Plot the fft spectrum with peaks. :param ax: The axis to plot on. :param peaks: The peak data. :param params: The ringdown parameters. + :param annotate: Whether to annotate the peaks. """ ax.clear() - ax.plot(peaks.freq, peaks.normalized_power, label="FFT Power", color="C0") + ax.plot(peaks.freq, peaks.power, label="FFT Power", color="C0") ax.plot( peaks.peak_freqs, - peaks.normalized_power[peaks.peaks], + peaks.power[peaks.peaks], "x", label="Peaks", color="C2", ) + + if annotate: + for i, (freq, height, lorentz) in enumerate( + zip(peaks.peak_freqs, peaks.power[peaks.peaks], peaks.lorentz_params) + ): + ax.annotate(f"{i} ({freq:.2e})", (freq, height)) + ax.plot( + peaks.freq, + lorentzian(peaks.freq, *lorentz), + "--", + color="C2", + alpha=0.5, + ) + ax.set_title("FFT Spectrum") ax.set_xlabel("ω [linear]") ax.set_ylabel("Power") - ax.set_ylim(0, 1.1) ax.axvline( params.fω_shift, color="gray", @@ -241,5 +257,4 @@ def plot_spectrum_and_peak_info(ax, peaks: RingdownPeakData, params: RingdownPar zorder=-10, label="Frequency Shift", ) - ax.legend() diff --git a/src/rabifun/system.py b/src/rabifun/system.py index 07f14f0..abd05ca 100644 --- a/src/rabifun/system.py +++ b/src/rabifun/system.py @@ -26,7 +26,13 @@ class Params: """Mode splitting in units of :any:`Ω`.""" η: float = 0.5 - """Decay rate :math:`\eta/2` of the system in angular frequency units.""" + """Decay rate :math:`\eta/2` of the system in linear frequency units. + + Squared amplitudes decay like :math:`exp(-2π η t)`. + """ + + η_hybrid: float | None = None + """Decay rate :math:`\eta/2` of the hybridized modes. If ``None``, the rate :any:`η` will be used.""" g_0: float = 0.01 """Drive amplitude in units of :any:`Ω`.""" @@ -100,6 +106,9 @@ class Params: if self.rwa: raise ValueError("Drive override is not compatible with the RWA.") + if self.η_hybrid is None: + self.η_hybrid = self.η + def periods(self, n: float): """ Returns the number of periods of the system that correspond to @@ -112,7 +121,7 @@ class Params: Returns the number of lifetimes of the system that correspond to `n` cycles. """ - return 2 * n / self.η + return 2 * n / (self.η * np.pi * 2) @property def rabi_splitting(self): @@ -151,7 +160,10 @@ class RuntimeParams: ) ) - decay_rates = -1j * np.repeat(params.η / 2, 2 * params.N + 2) + decay_rates = -1j * np.repeat(np.pi * params.η, 2 * params.N + 2) + if params.η_hybrid: + decay_rates[[0, 1]] = -1j * np.pi * np.repeat(params.η_hybrid, 2) + Ωs = freqs + decay_rates self.drive_frequencies, self.detunings, self.g, a_shift = ( @@ -160,10 +172,13 @@ class RuntimeParams: if params.drive_override is not None: self.drive_frequencies = params.drive_override[0] self.g = params.drive_override[1] - a_shift = 0 self.detunings *= 0 - self.g *= params.g_0 / np.sqrt(np.sum(self.g**2)) + + norm = np.sqrt(np.sum(self.g**2)) + + if norm > 0: + self.g *= params.g_0 self.g *= 2 * np.pi self.Ωs = Ωs @@ -183,7 +198,8 @@ class RuntimeParams: + decay_rates ) - self.bath_ε = 2 * np.pi * self.detunings - 1j * params.η / 2 + self.bath_ε = 2 * np.pi * self.detunings - 1j * 2 * np.pi * params.η / 2 + self.a_shift = 2 * np.pi * a_shift self.detuned_Ωs = freqs - self.ε.real @@ -295,7 +311,6 @@ def make_righthand_side(runtime_params: RuntimeParams, params: Params): 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 @@ -361,17 +376,16 @@ def output_signal(t: np.ndarray, amplitudes: np.ndarray, params: Params): """ runtime = RuntimeParams(params) - rotating = amplitudes * np.exp(-1j * runtime.detuned_Ωs[:, None] * t[None, :]) + + laser_times = ( + laser_frequency(params, t) + 2 * np.pi * params.measurement_detuning + ) * t + rotating = amplitudes.copy() * np.exp( + -1j * (runtime.detuned_Ωs[:, None] * t[None, :] - laser_times[None, :]) + ) rotating[0:2, :] *= runtime.a_weights[:, None].conjugate() - return ( - np.sum(rotating, axis=0) - * np.exp( - 1j - * (laser_frequency(params, t) + 2 * np.pi * params.measurement_detuning) - * t - ) - ).imag + return (np.sum(rotating, axis=0)).imag def bath_energies(N_couplings: int, ω_c: float) -> np.ndarray: diff --git a/src/ringfit/plotting.py b/src/ringfit/plotting.py index 54bb75c..7f22436 100644 --- a/src/ringfit/plotting.py +++ b/src/ringfit/plotting.py @@ -33,6 +33,7 @@ def plot_scan( normalize=False, smoothe_output: bool | float = False, max_frequency: float | bool = 0, + every=1, ax=None, **kwargs, ): @@ -60,7 +61,7 @@ def plot_scan( laser_data.max() - laser_data.min() ) - lines = ax.plot(time, laser_data, **kwargs) + lines = ax.plot(time[::every], laser_data[::every], **kwargs) if output: if normalize: @@ -68,7 +69,7 @@ def plot_scan( output_data.max() - output_data.min() ) - lines = ax.plot(time, output_data, **kwargs) + lines = ax.plot(time[::every], output_data[::every], **kwargs) if isinstance(steps, bool) and steps: peaks = data.laser_steps()