function plotPhaseDiagramWithBoundaries(PhaseDiagramMatrix, NumberOfAtoms, ScatteringLengths, PhaseBoundary_Untilted, TitleString)
% plotPhaseDiagramWithBoundaries - Plots a phase diagram with overlaid interpolated boundaries
%
% Inputs:
%   PhaseDiagramMatrix     - 2D matrix of phase values
%   NumberOfAtoms          - Vector of atom numbers (x-axis)
%   ScatteringLengths      - Vector of scattering lengths (y-axis)
%   PhaseBoundary_Untilted - Object containing selected boundary points
%   TitleString            - Title string for the plot
%
% This function displays the interpolated phase diagram and overlays
% boundary curves extracted from user-selected points.

    % Extract boundary point sets
    rawPoints = PhaseBoundary_Untilted.SelectedPoints.values;
    if iscell(rawPoints)
        pointSets = rawPoints;
    else
        pointSets = {rawPoints};
    end

    % Interpolate phase diagram
    [X, Y]   = meshgrid(NumberOfAtoms, ScatteringLengths);
    [Xq, Yq] = meshgrid(linspace(min(NumberOfAtoms), max(NumberOfAtoms), 500), ...
                        linspace(min(ScatteringLengths), max(ScatteringLengths), 500));
    Mq       = interp2(X, Y, PhaseDiagramMatrix, Xq, Yq, 'nearest');

    % Plot
    figure('Color','w');
    clf
    set(gcf, 'Position', [100, 100, 950, 800]);
    set(gcf, 'Renderer', 'opengl');
    imagesc([min(NumberOfAtoms), max(NumberOfAtoms)], ...
            [min(ScatteringLengths), max(ScatteringLengths)], Mq);
    set(gca, 'YDir', 'normal');
    colormap([
        0.8 0.8 0.8;  % Unmodulated
        0.2 0.6 1.0;  % SSD
        0.2 0.8 0.2;  % Stripe
        1.0 0.6 0.2;  % Labyrinth
        0.8 0.2 0.8   % Honeycomb
    ]);
    cb = colorbar('Ticks', 0:4, ...
             'TickLabels', {'Unmodulated','SSD','Stripe','Labyrinth','Honeycomb'}, ...
             'FontSize', 12);
    cb.Color = 'k';
    clim([0 4]);
    xlabel('Number of Atoms');
    ylabel('Scattering Length a_s (a_0)');
    t = title(TitleString);
    t.Color = 'k';

    % Overlay curves
    hold on;
    for i = 1:numel(pointSets)
        pts = pointSets{i};
        x = pts(:,1) * 1e5;  % convert back to # of atoms scale
        y = pts(:,2);

        % Interpolation
        [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, 1000);
        x_smooth = interp1(t, x, t_fine, 'makima');
        y_smooth = interp1(t, y, t_fine, 'makima');

        % Plot smooth curve
        plot(x_smooth, y_smooth, 'k-', 'LineWidth', 2);
    end

    set(gca, 'FontSize', 16, ...
            'XColor', 'k', 'YColor', 'k', ...
            'Color', 'none');
    axis tight;
    grid on;

end