function beautifyBERplot(options)
% BEAUTIFYBERPLOT  Enhances BER-style plots for publication-quality figures.
% Supports automatic smoothing and trend-line overlay.
%
% Usage examples:
%   beautifyBERplot;                                 % default
%   beautifyBERplot("polyfit",1);                    % add polynomial fit
%   beautifyBERplot("polyfit",1,"fitmethod","pchip") % piecewise cubic fit
%
% Supported fitmethod options: 'polyfit', 'smoothingspline', 'loess', 'pchip'

arguments
    options.logscale (1,1) logical = 1
    options.polyfit (1,1) logical = 0
    options.polyorder (1,1) double = 2
    options.fitmethod (1,1) string = "polyfit"  % choose fit type
end

% --- find all line objects in current axes
lines = findall(gca, 'Type', 'Line');
markers = {'o', 's', 'd', '^', 'v', '>', '<', 'p', 'h'}; 
num_markers = length(markers);

% --- style all lines consistently
for i = 1:length(lines)
    lines(i).LineWidth = 1;
    % lines(i).LineStyle = '-';
    if string(lines(i).Marker) == "none"
        lines(i).Marker = markers{mod(i-1, num_markers) + 1};
    end
    lines(i).MarkerSize = 4;
    lines(i).MarkerFaceColor = lines(i).Color;
end

% --- optional smoothing/fitting overlay
if options.polyfit
    hold on
    for i = 1:length(lines)
        x = lines(i).XData;
        y = lines(i).YData;
        valid = isfinite(x) & isfinite(y);
        if sum(valid) < options.polyorder + 1
            continue;
        end

        xf = linspace(min(x(valid)), max(x(valid)), 200);

        % ----- choose fitting method -----
        switch lower(options.fitmethod)
            case "polyfit"
                p = polyfit(x(valid), y(valid), options.polyorder);
                yf = polyval(p, xf);

            case "smoothingspline"
                try
                    f = fit(x(valid)', y(valid)', 'smoothingspline');
                    yf = feval(f, xf);
                catch
                    yf = interp1(x(valid), y(valid), xf, 'pchip');
                end

            case "loess"
                yf = smooth(x(valid), y(valid), 0.2, 'loess');
                yf = interp1(x(valid), yf, xf, 'linear', 'extrap');

            case "pchip"
                yf = interp1(x(valid), y(valid), xf, 'pchip');

            otherwise
                warning('Unknown fitmethod "%s". Using polyfit.', options.fitmethod);
                p = polyfit(x(valid), y(valid), options.polyorder);
                yf = polyval(p, xf);
        end

        % --- lightened color for fit overlay
        lightcol = lines(i).Color + 0.4 * (1 - lines(i).Color);
        lightcol(lightcol > 1) = 1;

        plot(xf, yf, '-', 'Color', lightcol, ...
            'LineWidth', 0.7, 'Marker', 'none', ...
            'HandleVisibility','off');
    end
    hold off
end

% --- axis scaling and cosmetics
if options.logscale
    set(gca, 'YScale', 'log');
end

% --- Figure size in centimeters ---
% width_pt  = 500;
% height_pt = 300;
% 
% pt2cm = 0.03514598;   % TeX point → cm
% width_cm  = width_pt  * pt2cm;   % = 8.85 cm
% height_cm = height_pt * pt2cm;   % = 2.81 cm
% 
% set(gcf, 'Units', 'centimeters', 'Position', [2 2 width_cm height_cm]);
% set(gcf, 'PaperUnits', 'centimeters', 'PaperPosition', [0 0 width_cm height_cm]);

% --- Formatting ---
set(findall(gca, '-property', 'Interpreter'), 'Interpreter', 'latex');
set(gcf, 'Color', 'w');
set(gca, 'Box', 'on', 'LineWidth', 0.8);
grid on;
set(gca, 'FontSize', 10, 'FontName', 'Latin Modern Roman');


end