mirror of
https://github.com/vale981/fibre_walk_project_code
synced 2025-03-04 09:21:38 -05:00
try to tease something out of the 11_07 data
This commit is contained in:
parent
bd9f8e1704
commit
8b4d0ec1bf
6 changed files with 296 additions and 170 deletions
|
|
@ -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}
|
||||
"""
|
||||
)
|
||||
|
|
|
|||
|
|
@ -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()
|
||||
|
|
|
|||
|
|
@ -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]
|
||||
|
|
|
|||
|
|
@ -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()
|
||||
|
|
|
|||
|
|
@ -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:
|
||||
|
|
|
|||
|
|
@ -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()
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue