function plotSmoothedContourOnPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, TitleString, varargin)
    interpMethod = 'makima';
    showPoints = true;
    matFilePaths = {};
    i = 1;
    while i <= numel(varargin)
        arg = varargin{i};
        if ischar(arg) || isstring(arg)
            switch lower(arg)
                case {'interp', 'interpmethod'}
                    interpMethod = varargin{i+1};
                    i = i + 2;
                case 'showpoints'
                    showPoints = varargin{i+1};
                    i = i + 2;
                otherwise
                    matFilePaths{end+1} = char(arg);
                    i = i + 1;
            end
        else
            error('Unexpected argument type');
        end
    end

    regionNames = ["Unmodulated", "SSD", "Stripe", "Labyrinth", "Honeycomb"];
    regionColors = [
        0.8 0.8 0.8;
        0.2 0.6 1.0;
        0.2 0.8 0.2;
        1.0 0.6 0.2;
        0.8 0.2 0.8
    ];

    [X, Y] = meshgrid(NUM_ATOMS_LIST, SCATTERING_LENGTH_RANGE);
    [Xq, Yq] = meshgrid(linspace(min(NUM_ATOMS_LIST), max(NUM_ATOMS_LIST), 500), ...
                        linspace(min(SCATTERING_LENGTH_RANGE), max(SCATTERING_LENGTH_RANGE), 500));
    Mq = interp2(X, Y, M, Xq, Yq, 'nearest');

    fig = figure('Color', 'w'); clf;
    set(fig, 'Position', [100, 100, 1050, 800], 'Renderer', 'opengl');
    ax = axes('Parent', fig, 'Position', [0.1 0.15 0.75 0.8]);

    imagesc(ax, [min(NUM_ATOMS_LIST), max(NUM_ATOMS_LIST)], ...
                [min(SCATTERING_LENGTH_RANGE), max(SCATTERING_LENGTH_RANGE)], Mq);
    set(ax, 'YDir', 'normal');
    colormap(ax, regionColors);
    colorbar(ax, 'Ticks', 0:4, 'TickLabels', regionNames, 'FontSize', 12);
    clim(ax, [0 4]);
    hold(ax, 'on');
    contour(ax, Xq, Yq, Mq, 0.5:1:3.5, 'k', 'LineWidth', 1.2);

    xlabel(ax, 'Number of Atoms', 'FontSize', 16);
    ylabel(ax, 'Scattering Length (\times a_o)', 'FontSize', 16);
    t = title(ax, TitleString, 'FontSize', 18, 'Interpreter', 'tex');
    t.Color = 'k';  % force title color to black
    set(ax, 'FontSize', 16, 'Color', 'none');
    axis(ax, 'tight');
    grid(ax, 'on');

    % Storage
    pointData = {}; % {struct with .x, .y, .points, .curve, .orig}
    draggingPoint = [];
    currentSet = 1;

    % Load points
    for k = 1:numel(matFilePaths)
        data = load(matFilePaths{k});
        if ~isfield(data, 'SelectedPoints') || ~isfield(data.SelectedPoints, 'values')
            warning('Skipping %s: missing SelectedPoints.values', matFilePaths{k});
            continue;
        end
        rawPoints = data.SelectedPoints.values;
        if iscell(rawPoints)
            pointSets = rawPoints;
        else
            pointSets = {rawPoints};
        end

        for i = 1:numel(pointSets)
            pts = pointSets{i};
            x = pts(:,1) * 1e5;
            y = pts(:,2);
            pointData{end+1} = plotPointSet(ax, x, y, interpMethod, showPoints);
        end
    end

    % Add control panel
    uicontrol('Style', 'pushbutton', 'String', 'Save', 'FontSize', 12, ...
              'Position', [860, 720, 150, 40], 'Callback', @(~,~) savePoints());
    uicontrol('Style', 'pushbutton', 'String', 'Export As...', 'FontSize', 12, ...
              'Position', [860, 670, 150, 40], 'Callback', @(~,~) exportPoints());
    uicontrol('Style', 'pushbutton', 'String', 'Reset All', 'FontSize', 12, ...
              'Position', [860, 620, 150, 40], 'Callback', @(~,~) resetPoints());

    % Mouse interaction
    set(fig, 'WindowButtonDownFcn', @startInteraction);
    set(fig, 'WindowButtonUpFcn', @stopDrag);
    set(fig, 'WindowButtonMotionFcn', @doDrag);

    function s = plotPointSet(ax, x, y, method, showPts)
        [x, idx] = sort(x); y = y(idx);
        dist = sqrt(diff(x).^2 + diff(y).^2);
        t = [0; cumsum(dist)];
        t = t / t(end);
        t_fine = linspace(0,1,500);
        x_smooth = interp1(t, x, t_fine, method);
        y_smooth = interp1(t, y, t_fine, method);

        s.x = x; s.y = y;
        s.orig = [x(:), y(:)];
        s.points = gobjects(numel(x),1);
        if showPts
            for j = 1:numel(x)
                s.points(j) = plot(ax, x(j), y(j), 'ro', 'MarkerFaceColor', 'r', ...
                    'ButtonDownFcn', @(src,~) pointClick(src));
            end
        end
        s.curve = plot(ax, x_smooth, y_smooth, 'w-', 'LineWidth', 2);
    end

    function pointClick(hPoint)
        modifiers = get(fig, 'CurrentModifier'); % get current key modifiers
        if any(strcmp(modifiers, 'shift'))  % Shift+Click deletes point
            deletePoint(hPoint);
        else
            startDrag(hPoint);  % normal click starts drag
        end
    end

    function deletePoint(hPoint)
        for i = 1:numel(pointData)
            idx = find(pointData{i}.points == hPoint);
            if ~isempty(idx)
                % Remove point data
                s = pointData{i};
                delete(s.points(idx)); % delete point graphic
                s.points(idx) = [];
                s.x(idx) = [];
                s.y(idx) = [];
                updateCurve(s);
                pointData{i} = s;
                break;
            end
        end
    end

    function startInteraction(~, ~)
        click = get(ax, 'CurrentPoint');
        if strcmp(fig.SelectionType, 'alt') % right click to add point
            pt = click(1,1:2);
            s = pointData{currentSet};
            s.x(end+1) = pt(1); s.y(end+1) = pt(2);
            s.points(end+1) = plot(ax, pt(1), pt(2), 'ro', ...
                'MarkerFaceColor', 'r', 'ButtonDownFcn', @(src,~) pointClick(src));
            [s.x, I] = sort(s.x); s.y = s.y(I); s.points = s.points(I);
            updateCurve(s);
            pointData{currentSet} = s;
        end
    end

    function startDrag(h)
        draggingPoint = h;
    end

    function doDrag(~, ~)
        if isempty(draggingPoint), return; end
        pt = get(ax, 'CurrentPoint');
        draggingPoint.XData = pt(1,1);
        draggingPoint.YData = pt(1,2);
        updateAssociatedCurve(draggingPoint);
    end

    function stopDrag(~, ~)
        draggingPoint = [];
    end

    function updateAssociatedCurve(hPoint)
        for i = 1:numel(pointData)
            idx = find(pointData{i}.points == hPoint);
            if ~isempty(idx)
                pointData{i}.x(idx) = hPoint.XData;
                pointData{i}.y(idx) = hPoint.YData;
                updateCurve(pointData{i});
                break;
            end
        end
    end

    function updateCurve(s)
        [s.x, idx] = sort(s.x); s.y = s.y(idx); s.points = s.points(idx);
        dist = sqrt(diff(s.x).^2 + diff(s.y).^2);
        t = [0; cumsum(dist)];
        if t(end) == 0, return; end
        t = t / t(end);
        t_fine = linspace(0,1,500);
        xs = interp1(t, s.x, t_fine, interpMethod);
        ys = interp1(t, s.y, t_fine, interpMethod);
        set(s.curve, 'XData', xs, 'YData', ys);
    end

    function savePoints()
        for i = 1:numel(pointData)
            pts = [pointData{i}.x(:)/1e5, pointData{i}.y(:)];
            SelectedPoints.values = pts;
            save(sprintf('ModifiedSet_%d.mat', i), 'SelectedPoints');
        end
        disp('Saved all modified point sets.');
    end

    function exportPoints()
        [f,p] = uiputfile('*.mat', 'Save As...');
        if isequal(f,0), return; end
        pts = cell(1, numel(pointData));  % Initialize as cell array
        for i = 1:numel(pointData)
            pts{i} = [pointData{i}.x(:)/1e5, pointData{i}.y(:)];
        end
        SelectedPoints.values = pts;
        save(fullfile(p,f), 'SelectedPoints');
        disp(['Exported to ', fullfile(p,f)]);
    end


    function resetPoints()
        for i = 1:numel(pointData)
            s = pointData{i};
            delete(s.points); delete(s.curve);
            pointData{i} = plotPointSet(ax, s.orig(:,1), s.orig(:,2), interpMethod, true);
        end
    end
end