New script to plot SS regimes and mods to analyze new data
This commit is contained in:
parent
dce4a0f067
commit
d874270a62
@ -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)
|
||||
|
||||
@ -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)
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
@ -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);
|
||||
|
||||
@ -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;
|
||||
|
||||
BIN
Estimations/20260605_ALSData.mat
Normal file
BIN
Estimations/20260605_ALSData.mat
Normal file
Binary file not shown.
430
Estimations/PhaseCoherenceMap.m
Normal file
430
Estimations/PhaseCoherenceMap.m
Normal file
@ -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
|
||||
Loading…
Reference in New Issue
Block a user