karthik/Calculations
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity
972d356a88
Calculations/Estimations/LinearizeScatteringLengthScan.m

231 lines
8.5 KiB
Matlab
Raw Blame History

FR_choice = 1;
ABKG_choice = 1;
a_start = 86; % initial scattering length (a_0)
a_end = 85; % final scattering length (a_0)
T = 0.030; % ramp duration (s)
Nt = 1000; % number of time points
ResonanceRange = [1.2, 2.6];
% Options for ramp generation
options = struct( ...
'smoothingMethod', 'sgolay', ... % 'sgolay', 'lowpass', or 'none'
'sgolayOrder', 3, ...
'sgolayFrameLength', 51, ...
'maxRampRate', 5, ...
'Bmin', 0.5, ...
'Bmax', 3.5, ...
'rampShape', 'linear' ... % Choose: 'linear', 'exponential', 'sigmoid'
);
[t, B_ramp, a_check] = generateSmoothBRamp(FR_choice, ABKG_choice, a_start, a_end, ResonanceRange, T, Nt, options);
% ----- Plot results -----
figure(1); clf
subplot(2,1,1);
plot(t, B_ramp, 'b-', 'LineWidth', 2);
xlabel('Time (s)');
ylabel('Magnetic field B(t) (G)');
title('Generated magnetic field ramp B(t)');
grid on;
set(gca, 'FontSize', 14);
subplot(2,1,2);
plot(t, a_check, 'r-', 'LineWidth', 2);
xlabel('Time (s)');
ylabel('Scattering length a_s(t) (a_0)');
legend('Resulting a_s(t)');
title('Scattering length ramp');
grid on;
set(gca, 'FontSize', 14);
% 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);
clf
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('Feshbach spectrum 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])
set(gca, 'FontSize', 14)
grid on;
%% 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
Reference in New Issue View Git Blame Copy Permalink