Estimations/LinearizeScatteringLengthScan.m

@@ -1,143 +0,0 @@
-% --- 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