Latest script to choose resonance range, ignore resonances between two chosen ones and generate a B Field ramp that results in a linear ramp in scattering length.
This commit is contained in:
parent
0ddc8bda00
commit
658ebb0496
@ -1,143 +1,188 @@
|
||||
% --- User setup ---
|
||||
resIndex = 1; % resonance index (1-based)
|
||||
duration = 0.5; % total ramp time [s]
|
||||
N_points = 300; % number of time points
|
||||
FR_choice = 1; % resonance data choice (1=new, 0=old)
|
||||
ABKG_choice = 1; % background scattering length choice (1,2,3)
|
||||
side = 'right'; % 'left' or 'right' side of resonance
|
||||
a_range = [111, 100]; % desired ramp of scattering length (a0), can be descending
|
||||
doPlot = true; % show resonance and ramp region
|
||||
% --- Main script ---
|
||||
FR_choice = 1;
|
||||
ABKG_choice = 1;
|
||||
|
||||
% --- Generate ramp ---
|
||||
[t, B_t, a_t] = makeLinearScatteringRamp(resIndex, duration, N_points, ...
|
||||
FR_choice, ABKG_choice, side, a_range, doPlot);
|
||||
% Ramp parameters
|
||||
a_start = 91; % initial scattering length (a_0)
|
||||
a_end = 89; % final scattering length (a_0)
|
||||
T = 0.030; % ramp duration (s)
|
||||
Nt = 1000; % number of time points
|
||||
ResonanceRange = [2.3, 3.0];
|
||||
% options.rampShape = @(t) a_start + (a_end - a_start) * sin(pi*t/T/2).^2;
|
||||
|
||||
% --- Plot results ---
|
||||
figure;
|
||||
% Options for ramp generation
|
||||
options = struct(...
|
||||
'smoothingMethod', 'sgolay', ...
|
||||
'sgolayOrder', 3, ...
|
||||
'sgolayFrameLength', 51, ...
|
||||
'maxRampRate', 5, ...
|
||||
'Bmin', 0.5, ...
|
||||
'Bmax', 3.5, ...
|
||||
'rampShape', 'linear' ... % 'linear', 'exponential', 'sigmoid', or function handle
|
||||
);
|
||||
|
||||
[t, B_ramp, a_check] = generateSmoothBRamp(FR_choice, ABKG_choice, a_start, a_end, ResonanceRange, T, Nt, options);
|
||||
|
||||
% Plot
|
||||
figure(1);
|
||||
subplot(2,1,1);
|
||||
plot(t, B_t, 'b-', 'LineWidth', 1.5);
|
||||
ylabel('B (G)', 'Interpreter', 'tex');
|
||||
title('Magnetic field B(t)', 'Interpreter', 'tex');
|
||||
plot(t, B_ramp, 'b-', 'LineWidth', 2);
|
||||
xlabel('Time (s)'); ylabel('B(t) (G)');
|
||||
title('Generated magnetic field ramp'); grid on;
|
||||
|
||||
subplot(2,1,2);
|
||||
plot(t, a_t, 'r-', 'LineWidth', 1.5);
|
||||
ylabel('Scattering length a ($a_0$)', 'Interpreter', 'latex');
|
||||
xlabel('Time (s)', 'Interpreter', 'tex');
|
||||
title('Linear ramp of a(t)', 'Interpreter', 'tex');
|
||||
plot(t, a_check, 'r--', 'LineWidth', 2);
|
||||
xlabel('Time (s)'); ylabel('a_s(B(t)) (a_0)');
|
||||
title('Resulting ramp in scattering length a_s(t)'); grid on;
|
||||
|
||||
% Visualize full resonance curve and selected window
|
||||
[B_full, a_full] = extractBetweenResonances(FR_choice, ABKG_choice, ResonanceRange);
|
||||
[B_curve, a_curve] = fullResonanceCurve(FR_choice, ABKG_choice, ResonanceRange);
|
||||
|
||||
figure(2);
|
||||
plot(B_curve, a_curve, 'k-', 'LineWidth', 1); hold on;
|
||||
plot(B_full, a_full, 'b-', 'LineWidth', 2); % zoomed-in region
|
||||
plot(B_ramp, a_check, 'r-', 'LineWidth', 2); % actual ramp
|
||||
plot(B_ramp([1 end]), a_check([1 end]), 'ro', 'MarkerSize', 8); % endpoints
|
||||
xline(min(B_ramp), '--b', 'B_{min}');
|
||||
xline(max(B_ramp), '--b', 'B_{max}');
|
||||
yline(min(a_check), '--r', 'a_{min}');
|
||||
yline(max(a_check), '--r', 'a_{max}');
|
||||
xlabel('B field (G)');
|
||||
ylabel('Scattering length a_s (a_0)');
|
||||
title('Resonance curve with selected B and a_s range', 'Interpreter','tex');
|
||||
legend('Full a(B)', 'Zoomed Region', 'Ramp a_s(t)', 'Ramp Endpoints', 'Location', 'NorthWest');
|
||||
ylim([0 150])
|
||||
grid on;
|
||||
|
||||
|
||||
%% Function: makeLinearScatteringRamp
|
||||
function [timeGrid, B_ramp, a_ramp] = makeLinearScatteringRamp(resIndex, duration, N_points, ...
|
||||
FR_choice, ABKG_choice, side, a_range, doPlot)
|
||||
% Generates B(t) that linearly ramps scattering length a(t) over time avoiding divergence.
|
||||
% Works with increasing or decreasing a_range.
|
||||
%% --- generateSmoothBRamp ---
|
||||
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);
|
||||
|
||||
% --- Resonance parameters ---
|
||||
if FR_choice == 1
|
||||
switch ABKG_choice
|
||||
case 1, a_bkg = 85.5;
|
||||
case 2, a_bkg = 93.5;
|
||||
case 3, a_bkg = 77.5;
|
||||
otherwise, error('Invalid ABKG_choice');
|
||||
% --- Generate base ramp shape ---
|
||||
ascending = (a_end > a_start);
|
||||
switch lower(class(opts.rampShape))
|
||||
case 'char'
|
||||
switch lower(opts.rampShape)
|
||||
case 'linear'
|
||||
base = t / T;
|
||||
case 'exponential'
|
||||
tau = T / 3;
|
||||
base = (1 - exp(-t / tau)) / (1 - exp(-T / tau));
|
||||
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)));
|
||||
otherwise
|
||||
error('Unknown ramp shape string: %s', opts.rampShape);
|
||||
end
|
||||
case 'function_handle'
|
||||
base = opts.rampShape(t);
|
||||
otherwise
|
||||
error('Invalid type for rampShape. Use string or function handle.');
|
||||
end
|
||||
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
|
||||
switch ABKG_choice
|
||||
case 1, a_bkg = 87.2;
|
||||
case 2, a_bkg = 95.2;
|
||||
case 3, a_bkg = 79.2;
|
||||
otherwise, error('Invalid ABKG_choice');
|
||||
|
||||
a_target = a_start + (a_end - a_start) * base;
|
||||
|
||||
% --- 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
|
||||
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];
|
||||
|
||||
% --- 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
|
||||
|
||||
% --- Scattering length formula ---
|
||||
a_of_B = @(B) arrayfun(@(b) ...
|
||||
a_bkg * prod(1 - resonancewB ./ (b - resonanceB)), ...
|
||||
B);
|
||||
%% --- extractBetweenResonances ---
|
||||
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);
|
||||
|
||||
% --- Select resonance ---
|
||||
B0 = resonanceB(resIndex);
|
||||
wB = resonancewB(resIndex);
|
||||
|
||||
% --- Define B scan range ---
|
||||
Bpad = 15 * wB;
|
||||
switch lower(side)
|
||||
case 'left'
|
||||
Bscan = linspace(B0 - Bpad, B0 - wB, 2000)';
|
||||
case 'right'
|
||||
Bscan = linspace(B0 + wB, B0 + Bpad, 2000)';
|
||||
otherwise
|
||||
error('side must be ''left'' or ''right''');
|
||||
end
|
||||
|
||||
% --- Evaluate scattering length over Bscan ---
|
||||
aScan = a_of_B(Bscan);
|
||||
|
||||
% --- Filter out any NaNs or extreme values ---
|
||||
valid = isfinite(aScan) & abs(aScan) < 1e4;
|
||||
Bscan = Bscan(valid);
|
||||
aScan = aScan(valid);
|
||||
|
||||
% --- Sort aScan and Bscan by ascending aScan ---
|
||||
[aSorted, idx] = sort(aScan);
|
||||
BSorted = Bscan(idx);
|
||||
|
||||
% --- Remove duplicates in aSorted for interp1 ---
|
||||
[uniqueA, iau] = unique(aSorted);
|
||||
uniqueB = BSorted(iau);
|
||||
|
||||
% --- Clip requested a_range to LUT limits ---
|
||||
requestedMin = min(a_range);
|
||||
requestedMax = max(a_range);
|
||||
|
||||
lutMin = uniqueA(1);
|
||||
lutMax = uniqueA(end);
|
||||
|
||||
if requestedMin < lutMin
|
||||
warning('Requested a_range min clipped to LUT minimum.');
|
||||
requestedMin = lutMin;
|
||||
end
|
||||
if requestedMax > lutMax
|
||||
warning('Requested a_range max clipped to LUT maximum.');
|
||||
requestedMax = lutMax;
|
||||
end
|
||||
|
||||
% --- Construct linear ramp in a ---
|
||||
% Use original order of a_range endpoints to respect increasing or decreasing ramp
|
||||
if a_range(1) < a_range(2)
|
||||
% ascending ramp
|
||||
a_lin = linspace(requestedMin, requestedMax, N_points);
|
||||
else
|
||||
% descending ramp
|
||||
a_lin = linspace(requestedMax, requestedMin, N_points);
|
||||
end
|
||||
|
||||
% --- Interpolate B from LUT ---
|
||||
B_lin = interp1(uniqueA, uniqueB, a_lin, 'pchip');
|
||||
|
||||
% --- Output time and ramp ---
|
||||
timeGrid = linspace(0, duration, N_points)';
|
||||
B_ramp = B_lin(:);
|
||||
a_ramp = a_lin(:);
|
||||
|
||||
% --- Optional plot ---
|
||||
if nargin >= 8 && doPlot
|
||||
Bvis = linspace(B0 - 20*wB, B0 + 20*wB, 2000);
|
||||
a_of_B = @(B) arrayfun(@(b) ...
|
||||
a_bkg * prod(1 - resonancewB ./ (b - resonanceB)), B);
|
||||
avis = a_of_B(Bvis);
|
||||
figure;
|
||||
plot(Bvis, avis, 'k-', 'DisplayName', 'Full a(B)');
|
||||
hold on;
|
||||
plot(B_ramp, a_ramp, 'r-', 'LineWidth', 2, 'DisplayName', 'Ramp a(t)');
|
||||
xlabel('B (G)');
|
||||
ylabel('a ($a_0$)', 'Interpreter', 'latex');
|
||||
title(sprintf('Feshbach Resonance #%d at B0=%.3f G (%s side)', resIndex, B0, side), 'Interpreter', 'tex');
|
||||
legend('Location','best');
|
||||
grid on;
|
||||
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
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user