From 99a9163a532c1a182644d98667c52b34300345ca Mon Sep 17 00:00:00 2001
From: Karthik <arecibokck@gmail.com>
Date: Mon, 9 Jun 2025 11:45:36 +0200
Subject: [PATCH] Script to compute a non-linear ramp of the B field that
 results in a linear scan of scattering lengths

---
 Estimations/LinearizeScatteringLengthScan.m | 143 ++++++++++++++++++++
 1 file changed, 143 insertions(+)
 create mode 100644 Estimations/LinearizeScatteringLengthScan.m

diff --git a/Estimations/LinearizeScatteringLengthScan.m b/Estimations/LinearizeScatteringLengthScan.m
new file mode 100644
index 0000000..7ef6933
--- /dev/null
+++ b/Estimations/LinearizeScatteringLengthScan.m
@@ -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