diff --git a/Estimations/LinearizeScatteringLengthScan.m b/Estimations/LinearizeScatteringLengthScan.m index 7a48f71..10e455b 100644 --- a/Estimations/LinearizeScatteringLengthScan.m +++ b/Estimations/LinearizeScatteringLengthScan.m @@ -1,45 +1,50 @@ -% --- Main script --- FR_choice = 1; ABKG_choice = 1; - -% Ramp parameters -a_start = 91; % initial scattering length (a_0) -a_end = 89; % final scattering length (a_0) -T = 0.030; % ramp duration (s) +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 = [2.3, 3.0]; -% options.rampShape = @(t) a_start + (a_end - a_start) * sin(pi*t/T/2).^2; +ResonanceRange = [1.2, 2.6]; % Options for ramp generation -options = struct(... - 'smoothingMethod', 'sgolay', ... +options = struct( ... + 'smoothingMethod', 'sgolay', ... % 'sgolay', 'lowpass', or 'none' 'sgolayOrder', 3, ... 'sgolayFrameLength', 51, ... 'maxRampRate', 5, ... 'Bmin', 0.5, ... 'Bmax', 3.5, ... - 'rampShape', 'linear' ... % 'linear', 'exponential', 'sigmoid', or function handle + '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 -figure(1); +% ----- Plot results ----- +figure(1); clf + subplot(2,1,1); plot(t, B_ramp, 'b-', 'LineWidth', 2); -xlabel('Time (s)'); ylabel('B(t) (G)'); -title('Generated magnetic field ramp'); grid on; +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('a_s(B(t)) (a_0)'); -title('Resulting ramp in scattering length a_s(t)'); grid on; +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 @@ -50,41 +55,25 @@ 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'); +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; - -%% --- generateSmoothBRamp --- +%% 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); - % --- 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.'); + 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 - a_target = a_start + (a_end - a_start) * base; + % 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); @@ -129,7 +118,30 @@ function [t, B_ramp, a_check] = generateSmoothBRamp(FR_choice, ABKG_choice, a_st a_check = a_of_B(B_ramp); end -%% --- extractBetweenResonances --- +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'); @@ -186,3 +198,33 @@ function [a_bkg, resonanceB, resonancewB] = getResonanceParams(FR_choice, ABKG_c 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