% --- 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

% --- 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)', 'Interpreter', 'tex');
title('Magnetic field B(t)', 'Interpreter', 'tex');

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');


%% 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$)', 'Interpreter', 'latex');
    title(sprintf('Feshbach Resonance #%d at B0=%.3f G (%s side)', resIndex, B0, side), 'Interpreter', 'tex');
    legend('Location','best');
    grid on;
end

end