Script to compute a non-linear ramp of the B field that results in a linear scan of scattering lengths
This commit is contained in:
Karthik
parent
ef2eab89ec
commit
99a9163a53
143
Estimations/LinearizeScatteringLengthScan.m
Normal file
143
Estimations/LinearizeScatteringLengthScan.m
Normal file
@ -0,0 +1,143 @@
|
|||||||
|
% --- User setup ---
|
||||||
|
resIndex = 3; % 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
|
||||||
|
|
||||||
|
% --- Generate ramp ---
|
||||||
|
[t, B_t, a_t] = makeLinearScatteringRamp(resIndex, duration, N_points, ...
|
||||||
|
FR_choice, ABKG_choice, side, a_range, doPlot);
|
||||||
|
|
||||||
|
% --- Plot results ---
|
||||||
|
figure;
|
||||||
|
subplot(2,1,1);
|
||||||
|
plot(t, B_t, 'b-', 'LineWidth', 1.5);
|
||||||
|
ylabel('B (G)');
|
||||||
|
title('Magnetic field B(t)');
|
||||||
|
|
||||||
|
subplot(2,1,2);
|
||||||
|
plot(t, a_t, 'r-', 'LineWidth', 1.5);
|
||||||
|
ylabel('Scattering length a (a_0)');
|
||||||
|
xlabel('Time (s)');
|
||||||
|
title('Linear ramp of a(t)');
|
||||||
|
|
||||||
|
|
||||||
|
%% 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.
|
||||||
|
|
||||||
|
% --- 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');
|
||||||
|
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');
|
||||||
|
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];
|
||||||
|
end
|
||||||
|
|
||||||
|
% --- Scattering length formula ---
|
||||||
|
a_of_B = @(B) arrayfun(@(b) ...
|
||||||
|
a_bkg * prod(1 - resonancewB ./ (b - resonanceB)), ...
|
||||||
|
B);
|
||||||
|
|
||||||
|
% --- 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);
|
||||||
|
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)');
|
||||||
|
title(sprintf('Feshbach Resonance #%d at B0=%.3f G (%s side)', resIndex, B0, side));
|
||||||
|
legend('Location','best');
|
||||||
|
grid on;
|
||||||
|
end
|
||||||
|
|
||||||
|
end
|
||||||
Loading…
Reference in New Issue
Block a user