diff --git a/Estimations/LinearizeScatteringLengthScan.m b/Estimations/LinearizeScatteringLengthScan.m
index 10e455b..271d2e5 100644
--- a/Estimations/LinearizeScatteringLengthScan.m
+++ b/Estimations/LinearizeScatteringLengthScan.m
@@ -42,10 +42,13 @@ 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);
+[B_all, a_all] = getFullFeschbachSpectrum(FR_choice, ABKG_choice);
 
 figure(2);
 clf
-plot(B_curve, a_curve, 'k-', 'LineWidth', 1); hold on;
+plot(B_all, a_all , 'g-', 'LineWidth', 1); hold on;
+plot(B_curve, a_curve, 'k-', 'LineWidth', 1); 
+% Plot all unfiltered resonances in black
 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
@@ -56,7 +59,8 @@ yline(max(a_check), '--r', 'a_{max}');
 xlabel('B field (G)');
 ylabel('Scattering length a_s (a_0)');
 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');
+legend('Full a(B)', 'Modified a(B)', 'Zoomed Region', 'Ramp a_s(t)', 'Ramp Endpoints', 'Location', 'NorthWest');
+xlim([1.15 2.7])
 ylim([0 150])
 set(gca, 'FontSize', 14)
 grid on;
@@ -178,6 +182,7 @@ function [a_bkg, resonanceB, resonancewB] = getResonanceParams(FR_choice, ABKG_c
         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);
@@ -228,3 +233,53 @@ function [t, B_ramp, a_check] = generateLinearBRampUsingLUT(FR_choice, ABKG_choi
     a_check = a_of_B(B_ramp);
 
 end
+
+function [B, a_s] = getFullFeschbachSpectrum(FR_choice, ABKG_choice)
+    % plotScatteringLength(FR_choice, ABKG_choice)
+    % Plots the Dy-164 scattering length vs B field based on PhysRevX.9.021012 Suppl. Fig. S5
+    % 
+    % Inputs:
+    %   FR_choice: 1 = new resonance parameters, 2 = old
+    %   ABKG_choice: 1 = lower a_bkg, 2 = mid, 3 = upper
+    %
+    % Example:
+    %   plotScatteringLength(1, 1);
+
+    if nargin < 1, FR_choice = 1; end
+    if nargin < 2, ABKG_choice = 1; end
+
+    % Choose background scattering length
+    if FR_choice == 1  % New values
+        switch ABKG_choice
+            case 1, a_bkg = 85.5;
+            case 2, a_bkg = 93.5;
+            case 3, a_bkg = 77.5;
+        end
+
+        % Resonance positions (G) and widths (G)
+        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.0010, 0.0005, 0.021, ...
+                       0.015, 0.043, 0.0005, 0.130, 0.024, 0.0005, 0.036, 3.1];
+
+    else  % Old values
+        switch ABKG_choice
+            case 1, a_bkg = 87.2;
+            case 2, a_bkg = 95.2;
+            case 3, a_bkg = 79.2;
+        end
+
+        % Resonance positions (G) and widths (G)
+        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
+
+    % Magnetic field range for plotting (G)
+    B = linspace(0.5, 8, 2000);
+
+    % Compute scattering length using product formula
+    a_s = a_bkg * ones(size(B));
+    for j = 1:length(resonanceB)
+        a_s = a_s .* (1 - resonancewB(j) ./ (B - resonanceB(j)));
+    end
+end