From 658ebb049665edcd8822328257db49b1b46e6e7b Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Mon, 16 Jun 2025 13:49:14 +0200 Subject: [PATCH] Latest script to choose resonance range, ignore resonances between two chosen ones and generate a B Field ramp that results in a linear ramp in scattering length. --- Estimations/LinearizeScatteringLengthScan.m | 297 +++++++++++--------- 1 file changed, 171 insertions(+), 126 deletions(-) diff --git a/Estimations/LinearizeScatteringLengthScan.m b/Estimations/LinearizeScatteringLengthScan.m index 32e6a85..7a48f71 100644 --- a/Estimations/LinearizeScatteringLengthScan.m +++ b/Estimations/LinearizeScatteringLengthScan.m @@ -1,143 +1,188 @@ -% --- 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 +% --- Main script --- +FR_choice = 1; +ABKG_choice = 1; -% --- Generate ramp --- -[t, B_t, a_t] = makeLinearScatteringRamp(resIndex, duration, N_points, ... - FR_choice, ABKG_choice, side, a_range, doPlot); +% Ramp parameters +a_start = 91; % initial scattering length (a_0) +a_end = 89; % final scattering length (a_0) +T = 0.030; % ramp duration (s) +Nt = 1000; % number of time points +ResonanceRange = [2.3, 3.0]; +% options.rampShape = @(t) a_start + (a_end - a_start) * sin(pi*t/T/2).^2; -% --- Plot results --- -figure; +% Options for ramp generation +options = struct(... + 'smoothingMethod', 'sgolay', ... + 'sgolayOrder', 3, ... + 'sgolayFrameLength', 51, ... + 'maxRampRate', 5, ... + 'Bmin', 0.5, ... + 'Bmax', 3.5, ... + 'rampShape', 'linear' ... % 'linear', 'exponential', 'sigmoid', or function handle +); + +[t, B_ramp, a_check] = generateSmoothBRamp(FR_choice, ABKG_choice, a_start, a_end, ResonanceRange, T, Nt, options); + +% Plot +figure(1); subplot(2,1,1); -plot(t, B_t, 'b-', 'LineWidth', 1.5); -ylabel('B (G)', 'Interpreter', 'tex'); -title('Magnetic field B(t)', 'Interpreter', 'tex'); +plot(t, B_ramp, 'b-', 'LineWidth', 2); +xlabel('Time (s)'); ylabel('B(t) (G)'); +title('Generated magnetic field ramp'); grid on; 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'); +plot(t, a_check, 'r--', 'LineWidth', 2); +xlabel('Time (s)'); ylabel('a_s(B(t)) (a_0)'); +title('Resulting ramp in scattering length a_s(t)'); grid on; + +% 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); +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('Resonance curve 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]) +grid on; -%% 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. +%% --- generateSmoothBRamp --- +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); -% --- 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'); + % --- Generate base ramp shape --- + ascending = (a_end > a_start); + switch lower(class(opts.rampShape)) + case 'char' + switch lower(opts.rampShape) + case 'linear' + base = t / T; + case 'exponential' + tau = T / 3; + base = (1 - exp(-t / tau)) / (1 - exp(-T / tau)); + 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))); + otherwise + error('Unknown ramp shape string: %s', opts.rampShape); + end + case 'function_handle' + base = opts.rampShape(t); + otherwise + error('Invalid type for rampShape. Use string or function handle.'); 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'); + + a_target = a_start + (a_end - a_start) * base; + + % --- 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 - 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]; + + % --- 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 -% --- Scattering length formula --- -a_of_B = @(B) arrayfun(@(b) ... - a_bkg * prod(1 - resonancewB ./ (b - resonanceB)), ... - B); +%% --- extractBetweenResonances --- +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); -% --- 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); + a_of_B = @(B) arrayfun(@(b) ... + a_bkg * prod(1 - resonancewB ./ (b - resonanceB)), B); 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; + 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 +