diff --git a/Data-Analyzer/compareSpectralWeights.m b/Data-Analyzer/compareSpectralWeights.m index fd1e852..b8206c8 100644 --- a/Data-Analyzer/compareSpectralWeights.m +++ b/Data-Analyzer/compareSpectralWeights.m @@ -8,7 +8,7 @@ format long font = 'Bahnschrift'; % Load data -Data = load('D:/Results - Experiment/B2.45G/DropletsToStripes.mat', 'unique_scan_parameter_values', 'mean_sc', 'stderr_sc', 'mean_sw', 'stderr_sw'); +Data = load('D:/Results - Experiment/B2.42G/DropletsToStripes.mat', 'unique_scan_parameter_values', 'mean_sc', 'stderr_sc', 'mean_sw', 'stderr_sw'); dts_scan_parameter_values = Data.unique_scan_parameter_values; dts_mean_sc = Data.mean_sc; @@ -16,7 +16,7 @@ dts_stderr_sc = Data.stderr_sc; dts_mean_sw = Data.mean_sw; dts_stderr_sw = Data.stderr_sw; -Data = load('D:/Results - Experiment/B2.45G/StripesToDroplets.mat', 'unique_scan_parameter_values', 'mean_sc', 'stderr_sc', 'mean_sw', 'stderr_sw'); +Data = load('D:/Results - Experiment/B2.42G/StripesToDroplets.mat', 'unique_scan_parameter_values', 'mean_sc', 'stderr_sc', 'mean_sw', 'stderr_sw'); std_scan_parameter_values = Data.unique_scan_parameter_values; std_mean_sw = Data.mean_sw; @@ -48,7 +48,7 @@ errorbar(std_scan_parameter_values, std_mean_sc, std_stderr_sc, 'o--', ... set(gca, 'FontSize', 14); % For tick labels only hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex'); hYLabel = ylabel('Spectral Contrast', 'Interpreter', 'tex'); -hTitle = title('B = 2.45 G', 'Interpreter', 'tex'); +hTitle = title('B = 2.42 G', 'Interpreter', 'tex'); legend set([hXLabel, hYLabel], 'FontName', font) set([hXLabel, hYLabel], 'FontSize', 14) @@ -65,7 +65,7 @@ errorbar(std_scan_parameter_values, std_mean_sw, std_stderr_sw, 'o--', ... set(gca, 'FontSize', 14); % For tick labels only hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex'); hYLabel = ylabel('Spectral Weight', 'Interpreter', 'tex'); -hTitle = title('B = 2.45 G', 'Interpreter', 'tex'); +hTitle = title('B = 2.42 G', 'Interpreter', 'tex'); legend set([hXLabel, hYLabel], 'FontName', font) set([hXLabel, hYLabel], 'FontSize', 14) diff --git a/Data-Analyzer/extractAutocorrelation.m b/Data-Analyzer/extractAutocorrelation.m index 073ab11..903319c 100644 --- a/Data-Analyzer/extractAutocorrelation.m +++ b/Data-Analyzer/extractAutocorrelation.m @@ -4,28 +4,30 @@ groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1_Axi "/images/ODT_2_Axis_Camera/in_situ_absorption", "/images/Horizontal_Axis_Camera/in_situ_absorption", ... "/images/Vertical_Axis_Camera/in_situ_absorption"]; -folderPath = "E:/Data - Experiment/2025/05/22/"; +folderPath = "//DyLabNAS/Data/TwoDGas/2025/06/23/"; -run = '0078'; +run = '0300'; folderPath = strcat(folderPath, run); cam = 5; angle = 0; -center = [1375, 2020]; +center = [1410, 2030]; span = [200, 200]; fraction = [0.1, 0.1]; pixel_size = 5.86e-6; removeFringes = false; -scan_parameter = 'rot_mag_fin_pol_angle'; +scan_parameter = 'ps_rot_mag_fin_pol_angle'; % scan_parameter = 'rot_mag_field'; scan_parameter_text = 'Angle = '; % scan_parameter_text = 'BField = '; -font = 'Bahnschrift'; +savefolderPath = 'D:/Results - Experiment/B2.42G/'; +savefileName = 'DropletsToStripes.mat'; +font = 'Bahnschrift'; skipPreprocessing = true; skipMasking = true; @@ -191,7 +193,7 @@ ylim([-1.5 3.0]); % Set y-axis limits here set(gca, 'FontSize', 14); hXLabel = xlabel('$\delta\theta / \pi$', 'Interpreter', 'latex'); hYLabel = ylabel('$g^{(2)}(\delta\theta)$', 'Interpreter', 'latex'); -hTitle = title('B = 2.45 G - Droplets to Stripes', 'Interpreter', 'tex'); +hTitle = title('Change across transition', 'Interpreter', 'tex'); legend(legend_entries, 'Interpreter', 'latex', 'Location', 'bestoutside'); set([hXLabel, hYLabel], 'FontName', font) set([hXLabel, hYLabel], 'FontSize', 14) diff --git a/Data-Analyzer/extractSpectralWeight.m b/Data-Analyzer/extractSpectralWeight.m index accaf48..5842236 100644 --- a/Data-Analyzer/extractSpectralWeight.m +++ b/Data-Analyzer/extractSpectralWeight.m @@ -4,31 +4,38 @@ groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1_Axi "/images/ODT_2_Axis_Camera/in_situ_absorption", "/images/Horizontal_Axis_Camera/in_situ_absorption", ... "/images/Vertical_Axis_Camera/in_situ_absorption"]; -folderPath = "D:/Data - Experiment/2025/05/22/"; +folderPath = "//DyLabNAS/Data/TwoDGas/2025/06/24/"; -run = '0079'; +run = '0001'; folderPath = strcat(folderPath, run); cam = 5; angle = 0; -center = [1375, 2020]; +center = [1410, 2030]; span = [200, 200]; fraction = [0.1, 0.1]; pixel_size = 5.86e-6; removeFringes = false; -scan_parameter = 'rot_mag_fin_pol_angle'; +scan_parameter = 'ps_rot_mag_fin_pol_angle'; % scan_parameter = 'rot_mag_field'; scan_parameter_text = 'Angle = '; % scan_parameter_text = 'BField = '; -savefodlerPath = 'D:/Results - Experiment/B2.45G/'; -savefileName = 'StripesToDroplets.mat'; +savefolderPath = 'D:/Results - Experiment/B2.42G/'; +savefileName = 'StripesToDroplets'; font = 'Bahnschrift'; +skipUnshuffling = false; +if strcmp(savefileName, 'DropletsToStripes') + scan_groups = 0:5:45; +else + scan_groups = 45:-5:0; +end + skipPreprocessing = true; skipMasking = true; skipIntensityThresholding = true; @@ -86,7 +93,7 @@ for k = 1 : length(files) info = h5info(fullFileName, '/globals'); for i = 1:length(info.Attributes) if strcmp(info.Attributes(i).Name, scan_parameter) - if strcmp(scan_parameter, 'rot_mag_fin_pol_angle') + if strcmp(scan_parameter, 'ps_rot_mag_fin_pol_angle') scan_parameter_values(k) = 180 - info.Attributes(i).Value; else scan_parameter_values(k) = info.Attributes(i).Value; @@ -95,6 +102,43 @@ for k = 1 : length(files) end end +%% Unshuffle if necessary to do so + +if ~skipUnshuffling + n_values = length(scan_groups); + n_total = length(scan_parameter_values); + + % Infer number of repetitions + n_reps = n_total / n_values; + + % Preallocate ordered arrays + ordered_scan_values = zeros(1, n_total); + ordered_od_imgs = cell(1, n_total); + + counter = 1; + + for rep = 1:n_reps + for val = scan_groups + % Find the next unused match for this val + idx = find(scan_parameter_values == val, 1, 'first'); + + % Assign and remove from list to avoid duplicates + ordered_scan_values(counter) = scan_parameter_values(idx); + ordered_od_imgs{counter} = od_imgs{idx}; + + % Mark as used by removing + scan_parameter_values(idx) = NaN; % NaN is safe since original values are 0:5:45 + od_imgs{idx} = []; % empty cell so it won't be matched again + + counter = counter + 1; + end + end + + % Now assign back + scan_parameter_values = ordered_scan_values; + od_imgs = ordered_od_imgs; +end + %% Run Fourier analysis over images fft_imgs = cell(1, nimgs); @@ -107,7 +151,7 @@ Sigma = 2; N_shots = length(od_imgs); % Create VideoWriter object for movie -videoFile = VideoWriter('Single_Shot_FFT.mp4', 'MPEG-4'); +videoFile = VideoWriter([savefileName '.mp4'], 'MPEG-4'); videoFile.Quality = 100; % Set quality to maximum (0–100) videoFile.FrameRate = 2; % Set the frame rate (frames per second) open(videoFile); % Open the video file to write @@ -280,7 +324,7 @@ set([hXLabel, hYLabel], 'FontSize', 14) set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title grid on -save([savefolderPath savefileName], 'unique_scan_parameter_values', 'mean_sc', 'stderr_sc', 'mean_sw', 'stderr_sw'); +save([savefolderPath savefileName '.mat'], 'unique_scan_parameter_values', 'mean_sc', 'stderr_sc', 'mean_sw', 'stderr_sw'); %% k-means Clustering diff --git a/Data-Analyzer/matchTFRadiiAtomNumber.m b/Data-Analyzer/matchTFRadiiAtomNumber.m index 871d1f4..9d17f2a 100644 --- a/Data-Analyzer/matchTFRadiiAtomNumber.m +++ b/Data-Analyzer/matchTFRadiiAtomNumber.m @@ -159,5 +159,5 @@ ylabel('TF Radius - Y ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16); legend('FontSize', 12, 'Interpreter', 'tex', 'Location', 'bestoutside'); axis square; grid on; -sgtitle('[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz', ... +sgtitle('[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 36.87, 17.4, 93.16 ] Hz', ... 'Interpreter', 'tex', 'FontSize', 18); diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m index c41a1aa..78003e8 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -614,11 +614,6 @@ JobNumber = 0; SuppressPlotFlag = false; UCD = Scripts.extractAverageUnitCellDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag); -NUM_ATOMS_LIST = [100000 304167 508333 712500 916667 1120833 1325000 ... - 1529167 1733333 1937500 2141667 2345833 2550000 2754167 ... - 2958333 3162500 3366667 3570833 3775000 3979167 4183333 ... - 4387500 4591667 4795833 5000000]; - %% Plot number of droplets % Parameters Radius = 2; @@ -780,7 +775,7 @@ TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, SCATTERING_LENGTH_RANGE = [95.62]; -NUM_ATOMS_LIST = [712500 916667 1120833 1325000 1529167 1733333 1937500 2141667 2345833 2550000 2754167 2958333 3162500 3366667 3570833]; +NUM_ATOMS_LIST = [712500 916667 1120833 1325000 1529167 1937500 2141667 2345833 2550000 2754167 2958333 3162500 3366667 3570833]; UCD_values = zeros(length(SCATTERING_LENGTH_RANGE), length(NUM_ATOMS_LIST)); @@ -834,7 +829,7 @@ PlanckConstant = 6.62607015E-34; PlanckConstantReduced = 6.62607015E-34/(2*pi); AtomicMassUnit = 1.660539066E-27; BohrMagneton = 9.274009994E-24; - +VacuumPermeability = 1.25663706212E-6; % Dy specific constants Dy164Mass = 163.929174751*AtomicMassUnit; Dy164IsotopicAbundance = 0.2826; diff --git a/Estimations/20260605_ALSData.mat b/Estimations/20260605_ALSData.mat new file mode 100644 index 0000000..21d85f7 Binary files /dev/null and b/Estimations/20260605_ALSData.mat differ diff --git a/Estimations/PhaseCoherenceMap.m b/Estimations/PhaseCoherenceMap.m new file mode 100644 index 0000000..376de5b --- /dev/null +++ b/Estimations/PhaseCoherenceMap.m @@ -0,0 +1,430 @@ +FR_choice = 1; +ABKG_choice = 1; + +% Get the full spectrum +[B, a_s] = getFullFeschbachSpectrum(FR_choice, ABKG_choice); + +% Define the plotting range +x_limits = [1.15, 2.7]; +y_limits = [0, 150]; + +% Find indices within our x-range +mask = (B >= x_limits(1)) & (B <= x_limits(2)); +B_plot = B(mask); +a_s_plot = a_s(mask); + +% Identify resonances to mask (looking at the spectrum, we'll mask around 2.174G and 2.336G) +resonances_to_mask = [2.174, 2.336]; +mask_width = 0.1; % Width to mask around each resonance (in G) + +% Create a mask for the regions to keep (1 = keep, 0 = mask) +keep_mask = true(size(B_plot)); +for res = resonances_to_mask + keep_mask = keep_mask & (B_plot < (res - mask_width) | B_plot > (res + mask_width)); +end + +% Create masked version +B_masked = B_plot(keep_mask); +a_s_masked = a_s_plot(keep_mask); + +% Interpolate over the masked regions +a_s_interp = interp1(B_masked, a_s_masked, B_plot, 'pchip'); + +% --- Find the remaining resonances (peaks) in the interpolated data --- +% Use findpeaks to detect peaks (adjust MinPeakProminence as needed) +[peaks, locs] = findpeaks(abs(a_s_interp), 'MinPeakProminence', 20); % Tune this! +remaining_resonances = B_plot(locs); % Positions of remaining resonances + +% --- Select the two resonances of interest (1.3G and 2.59G) --- +% If automatic detection fails, manually define them: +target_resonances = [1.3, 2.59]; % Manually specified (adjust as needed) +% OR use the closest detected resonances: +% target_resonances = remaining_resonances(abs(remaining_resonances - 1.3) < 0.2 & abs(remaining_resonances - 2.59) < 0.2); + +% --- Extract the curve BETWEEN these two resonances --- +res1 = min(target_resonances); % Lower resonance (1.3G) +res2 = max(target_resonances); % Upper resonance (2.59G) +between_mask = (B_plot > res1) & (B_plot < res2); % No masking here, since resonances are already removed +B_between = B_plot(between_mask); +a_s_between = a_s_interp(between_mask); + +% --- Plotting --- +figure(1); +set(gcf,'Position',[100 100 950 750]) +set(gca,'FontSize',16,'Box','On','Linewidth',2); +hold on; + +% Plot the interpolated spectrum (without 2.174G and 2.336G) +plot(B_plot, a_s_interp, 'k-', 'LineWidth', 1, 'DisplayName', 'Interpolated spectrum'); + +% Highlight the region between 1.3G and 2.59G +plot(B_between, a_s_between, 'm-', 'LineWidth', 2, 'DisplayName', 'Between 1.3G and 2.59G'); + +% Formatting +xlim(x_limits); +ylim([0, 150]); +xlabel('Magnetic Field (G)'); +ylabel('Scattering Length (a_0)'); +title('Interpolated Spectrum: Region Between 1.3G and 2.59G'); +legend('Location', 'southeast'); +grid on; +hold off; +%% +% Load data +FSData = load('20260605_ALSData.mat'); +x = FSData.data(:,1); +y = FSData.data(:,2) * 166 * 1E-5; +y_err = FSData.data(:,3) * 166 * 1E-5; + +% --- Plotting --- +figure(2); clf; +set(gcf,'Position',[100 100 950 750]) +set(gca,'FontSize',16,'Box','On','Linewidth',2); +hold on; + +% Plot data with error bars, no connecting lines (only points) +errorbar(x, y, y_err, 'ko', ... + 'MarkerFaceColor', 'r', ... % Filled black circles + 'MarkerSize', 6, ... % Size of data points + 'LineStyle', 'none', ... % No connecting line + 'LineWidth', 1.5, ... % Width of error bars + 'DisplayName', 'As measured with ALS on a cold thermal cloud'); + +% Formatting +xlabel('Magnetic Field (G)'); +ylabel('Atom Number'); +legend('Location', 'southeast'); +grid on; +box on; +hold off; +%% Combined plot + +figure(3); clf; +set(gcf, 'Position', [100, 100, 950, 750]); +set(gca, 'FontSize', 16, 'Box', 'On', 'LineWidth', 2); +hold on; + +% ========== Define Shaded Regions ========== +% Format: [x_start, x_end, label] +regions = [ + 2.45, 2.52, "BEC"; + 2.35, 2.45, "SSD/SS"; + 2.1, 2.35, "ID/IS"; + 1.4, 2.1, "TBD" +]; + +% Colors for each region (semi-transparent) +region_colors = [ + 0.85, 0.9, 1.0; % light blue + 0.9, 1.0, 0.85; % light green + 1.0, 0.9, 0.85; % light red + 0.95, 0.95, 0.95 % light gray +]; + +ylims = [0, 150]; % Adjust as needed + +for i = 1:size(regions,1) + x1 = str2double(regions(i,1)); + x2 = str2double(regions(i,2)); + label = string(regions(i,3)); + fill([x1 x2 x2 x1], [ylims(1) ylims(1) ylims(2) ylims(2)], ... + region_colors(i,:), ... + 'FaceAlpha', 0.9, 'EdgeColor', 'none', 'HandleVisibility', 'off'); + text((x1 + x2)/2, mean(ylims), label, ... + 'HorizontalAlignment', 'center', ... + 'VerticalAlignment', 'middle', ... + 'FontSize', 14, 'FontWeight', 'bold', ... + 'Rotation', 90); % Vertical orientation + +end + +% ========== LEFT Y-AXIS ========== +yyaxis left +plot(B_plot, a_s_interp, 'k-', 'LineWidth', 1.5, 'DisplayName', 'Interpolated Dy-164 Feshbach spectrum'); +plot(B_between, a_s_between, 'm-', 'LineWidth', 2.5, 'DisplayName', 'Region between 1.3G and 2.59G'); + +ylabel('Scattering Length (a_0)'); +ylim(ylims); + +% ========== RIGHT Y-AXIS ========== +yyaxis right +ax = gca; +ax.YColor = [0.75 0 0]; % Make right y-axis red to match data points +FSData = load('20260605_ALSData.mat'); +x = FSData.data(:,1); +y = FSData.data(:,2) * 166 * 1E-5; +y_err = FSData.data(:,3) * 166 * 1E-5; + +errorbar(x, y, y_err, 'ko', ... + 'MarkerFaceColor', 'r', ... + 'MarkerSize', 6, ... + 'LineStyle', 'none', ... + 'LineWidth', 1.5, ... + 'DisplayName', 'As measured with ALS on a cold thermal cloud'); + +ylabel('Atom Number'); + +% ========== Shared Formatting ========== +xlabel('Magnetic Field (G)'); +xlim([1.15, 2.7]); +grid on; + +legend('Location', 'southeast'); + +title('BEC-SSD/SS-ID/IS Regimes'); + +hold off; + + + +%% Helper functions +function [t, B_ramp, a_check] = generateSmoothBRamp(FR_choice, ABKG_choice, a_start, a_end, selectedResRange, T, Nt, opts) + % Time array + t = linspace(0, T, Nt); + + if strcmp(opts.rampShape, 'linear') + % Directly call helper for linear LUT ramp generation + [t, B_ramp, a_check] = generateLinearBRampUsingLUT(FR_choice, ABKG_choice, a_start, a_end, selectedResRange, T, Nt); + return + end + + % Get target a_s(t) based on rampShape choice + a_target = getTargetScatteringLength(t, T, a_start, a_end, opts.rampShape); + + % --- a(B) interpolation --- + [B_between, a_between] = extractBetweenResonances(FR_choice, ABKG_choice, selectedResRange); + valid_idx = a_between > 0 & a_between < 150; + [a_sorted, sort_idx] = sort(a_between(valid_idx)); + B_sorted = B_between(valid_idx); + B_sorted = B_sorted(sort_idx); + B_of_a = @(a) interp1(a_sorted, B_sorted, a, 'linear', 'extrap'); + B_raw = B_of_a(a_target); + + % --- Smoothing --- + switch opts.smoothingMethod + case 'sgolay' + B_smooth = sgolayfilt(B_raw, opts.sgolayOrder, opts.sgolayFrameLength); + case 'lowpass' + dt = T / (Nt - 1); Fs = 1 / dt; + B_smooth = lowpass(B_raw, Fs / 20, Fs); + case 'none' + B_smooth = B_raw; + otherwise + error('Unknown smoothing method'); + end + + % --- Bound the ramp --- + B_smooth = min(max(B_smooth, opts.Bmin), opts.Bmax); + + % --- Enforce max dB/dt --- + dt = T / (Nt - 1); + for i = 2:Nt + delta = B_smooth(i) - B_smooth(i-1); + if abs(delta/dt) > opts.maxRampRate + delta = sign(delta) * opts.maxRampRate * dt; + B_smooth(i) = B_smooth(i-1) + delta; + end + end + B_ramp = B_smooth; + + % --- Verify a_s(t) from B_ramp --- + [a_bkg, resonanceB, resonancewB] = getResonanceParams(FR_choice, ABKG_choice, selectedResRange); + a_of_B = @(B) arrayfun(@(b) ... + a_bkg * prod(1 - resonancewB ./ (b - resonanceB)), B); + a_check = a_of_B(B_ramp); +end + +function a_target = getTargetScatteringLength(t, T, a_start, a_end, rampShape) + % If rampShape is a function handle, use it directly (for flexibility) + if isa(rampShape, 'function_handle') + a_target = rampShape(t); + return + end + + switch lower(rampShape) + case 'exponential' + tau = T / 3; + base = (1 - exp(-t / tau)) / (1 - exp(-T / tau)); + a_target = a_start + (a_end - a_start) * base; + + case 'sigmoid' + s = 10 / T; center = T / 2; + sigmoid = @(x) 1 ./ (1 + exp(-s * (x - center))); + base = (sigmoid(t) - sigmoid(t(1))) / (sigmoid(t(end)) - sigmoid(t(1))); + a_target = a_start + (a_end - a_start) * base; + + otherwise + error('Unknown ramp shape: %s', rampShape); + end +end + +function [B_between, a_between] = extractBetweenResonances(FR_choice, ABKG_choice, selectedRange) + [a_bkg, resonanceB, resonancewB] = getResonanceParams(FR_choice, ABKG_choice, selectedRange); + [~, idx] = sort(resonancewB, 'descend'); + B1 = resonanceB(idx(1)); B2 = resonanceB(idx(2)); + w1 = resonancewB(idx(1)); w2 = resonancewB(idx(2)); + Bvis = linspace(min(B1, B2) - 20*min(w1,w2), max(B1, B2) + 20*min(w1,w2), 2000); + + a_of_B = @(B) arrayfun(@(b) ... + a_bkg * prod(1 - resonancewB ./ (b - resonanceB)), B); + avis = a_of_B(Bvis); + between_idx = Bvis >= min(B1,B2) & Bvis <= max(B1,B2); + B_between = Bvis(between_idx); + a_between = avis(between_idx); +end + +function [B_range, a_values] = fullResonanceCurve(FR_choice, ABKG_choice, selectedRange) + [a_bkg, resonanceB, resonancewB] = getResonanceParams(FR_choice, ABKG_choice, selectedRange); + B_range = linspace(min(resonanceB)-0.2, max(resonanceB)+0.2, 3000); + a_of_B = @(B) arrayfun(@(b) ... + a_bkg * prod(1 - resonancewB ./ (b - resonanceB)), B); + a_values = a_of_B(B_range); +end + +function [a_bkg, resonanceB, resonancewB] = getResonanceParams(FR_choice, ABKG_choice, selectedRange) + if FR_choice == 1 + a_bkg_list = [85.5, 93.5, 77.5]; + resonanceB = [1.295, 1.306, 2.174, 2.336, 2.591, 2.740, 2.803, ... + 2.780, 3.357, 4.949, 5.083, 7.172, 7.204, 7.134, 76.9]; + resonancewB = [0.009, 0.010, 0.0005, 0.0005, 0.001, 0.0005, 0.021, ... + 0.015, 0.043, 0.0005, 0.130, 0.024, 0.0005, 0.036, 3.1]; + else + a_bkg_list = [87.2, 95.2, 79.2]; + resonanceB = [1.298, 2.802, 3.370, 5.092, 7.154, 2.592, 2.338, 2.177]; + resonancewB = [0.018, 0.047, 0.048, 0.145, 0.020, 0.008, 0.001, 0.001]; + end + a_bkg = a_bkg_list(ABKG_choice); + + % --- Filter resonanceB and resonancewB if selectedRange is provided --- + if nargin >= 3 && ~isempty(selectedRange) + minB = min(selectedRange); maxB = max(selectedRange); + keep_idx = (resonanceB >= minB) & (resonanceB <= maxB); + + % Keep only the lowest and highest resonance in the selected range + if sum(keep_idx) >= 2 + B_sub = resonanceB(keep_idx); + w_sub = resonancewB(keep_idx); + [~, idx_lo] = min(B_sub); + [~, idx_hi] = max(B_sub); + + resonanceB = [B_sub(idx_lo), B_sub(idx_hi)]; + resonancewB = [w_sub(idx_lo), w_sub(idx_hi)]; + else + error('Selected resonance range does not include at least two resonances.'); + end + end +end + +function [t, B_ramp, a_check] = generateLinearBRampUsingLUT(FR_choice, ABKG_choice, a_start, a_end, selectedResRange, T, Nt) + + % Time vector + t = linspace(0, T, Nt); + + % 1) Generate LUT of B and a_s(B) + [B_between, a_between] = extractBetweenResonances(FR_choice, ABKG_choice, selectedResRange); + + % Restrict to physically meaningful range (optional) + valid_idx = a_between > 0 & a_between < 150; + B_lut = B_between(valid_idx); + a_lut = a_between(valid_idx); + + % 2) Generate linear a_s ramp in time + a_target = linspace(a_start, a_end, Nt); + + % 3) Interpolate B(t) from LUT a_s -> B + % Make sure a_lut is sorted ascending + [a_lut_sorted, idx_sort] = sort(a_lut); + B_lut_sorted = B_lut(idx_sort); + + B_ramp = interp1(a_lut_sorted, B_lut_sorted, a_target, 'linear', 'extrap'); + + % 4) Compute resulting a_s(t) for verification + [a_bkg, resonanceB, resonancewB] = getResonanceParams(FR_choice, ABKG_choice, selectedResRange); + a_of_B = @(B) arrayfun(@(b) ... + a_bkg * prod(1 - resonancewB ./ (b - resonanceB)), B); + a_check = a_of_B(B_ramp); + +end + +function [B, a_s] = getFullFeschbachSpectrum(FR_choice, ABKG_choice) + % plotScatteringLength(FR_choice, ABKG_choice) + % Plots the Dy-164 scattering length vs B field based on PhysRevX.9.021012 Suppl. Fig. S5 + % + % Inputs: + % FR_choice: 1 = new resonance parameters, 2 = old + % ABKG_choice: 1 = lower a_bkg, 2 = mid, 3 = upper + % + % Example: + % plotScatteringLength(1, 1); + + if nargin < 1, FR_choice = 1; end + if nargin < 2, ABKG_choice = 1; end + + % Choose background scattering length + if FR_choice == 1 % New values + switch ABKG_choice + case 1, a_bkg = 85.5; + case 2, a_bkg = 93.5; + case 3, a_bkg = 77.5; + end + + % Resonance positions (G) and widths (G) + resonanceB = [1.295, 1.306, 2.174, 2.336, 2.591, 2.740, 2.803, 2.780, ... + 3.357, 4.949, 5.083, 7.172, 7.204, 7.134, 76.9]; + resonancewB = [0.009, 0.010, 0.0005, 0.0005, 0.0010, 0.0005, 0.021, ... + 0.015, 0.043, 0.0005, 0.130, 0.024, 0.0005, 0.036, 3.1]; + + else % Old values + switch ABKG_choice + case 1, a_bkg = 87.2; + case 2, a_bkg = 95.2; + case 3, a_bkg = 79.2; + end + + % Resonance positions (G) and widths (G) + resonanceB = [1.298, 2.802, 3.370, 5.092, 7.154, 2.592, 2.338, 2.177]; + resonancewB = [0.018, 0.047, 0.048, 0.145, 0.020, 0.008, 0.001, 0.001]; + end + + % Magnetic field range for plotting (G) + B = linspace(0.5, 8, 2000); + + % Compute scattering length using product formula + a_s = a_bkg * ones(size(B)); + for j = 1:length(resonanceB) + a_s = a_s .* (1 - resonancewB(j) ./ (B - resonanceB(j))); + end +end + +function [t, B_ramp, a_check] = generateLinearBRamp(B_between, a_s_between, a_start, a_end, T, Nt) + % Generates a B-field ramp (B_ramp) to produce a linear a_s ramp from a_start to a_end. + % Uses precomputed LUT: B_between (magnetic field) and a_s_between (scattering length). + % + % Inputs: + % B_between - Magnetic field values (G) between resonances [vector] + % a_s_between - Corresponding scattering lengths (a_0) [vector] + % a_start, a_end - Target scattering length range (a_0) + % T - Total ramp time (s) + % Nt - Number of time steps + % + % Outputs: + % t - Time vector (s) + % B_ramp - Generated B-field ramp (G) + % a_check - Verified a_s(t) using interpolation (a_0) + + % --- 1. Time vector --- + t = linspace(0, T, Nt); + + % --- 2. Ensure LUT is sorted and unique --- + [a_s_sorted, sort_idx] = unique(a_s_between); % Remove duplicates and sort + B_sorted = B_between(sort_idx); + + % --- 3. Generate target linear a_s ramp --- + a_target = linspace(a_start, a_end, Nt); + + % --- 4. Interpolate B(t) from a_s -> B --- + B_ramp = interp1(a_s_sorted, B_sorted, a_target, 'pchip', 'extrap'); + + % --- 5. Verify a_s(t) by re-interpolating --- + a_check = interp1(B_sorted, a_s_sorted, B_ramp, 'pchip'); +end