diff --git a/Dipolar-Gas-Simulator/+Scripts/plotPhaseDiagramWithBoundaries.m b/Dipolar-Gas-Simulator/+Scripts/plotPhaseDiagramWithBoundaries.m
new file mode 100644
index 0000000..2235d26
--- /dev/null
+++ b/Dipolar-Gas-Simulator/+Scripts/plotPhaseDiagramWithBoundaries.m
@@ -0,0 +1,79 @@
+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
diff --git a/Dipolar-Gas-Simulator/+Scripts/plotSmoothedContourOnPhaseDiagram.m b/Dipolar-Gas-Simulator/+Scripts/plotSmoothedContourOnPhaseDiagram.m
new file mode 100644
index 0000000..99cc8ed
--- /dev/null
+++ b/Dipolar-Gas-Simulator/+Scripts/plotSmoothedContourOnPhaseDiagram.m
@@ -0,0 +1,230 @@
+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);
+    title(ax, titleString, 'FontSize', 18, 'Interpreter', 'tex');
+    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
diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m
index 769f659..87c2bc8 100644
--- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m
+++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m
@@ -751,12 +751,47 @@ load('./Results/Data_Full3D/PhaseDiagramTilted_Theta_20.mat')
 PhaseDiagramMatrix   = M;
 ScatteringLengths    = SCATTERING_LENGTH_RANGE;
 NumberOfAtoms        = NUM_ATOMS_LIST;
+
 Scripts.editPhaseDiagram(PhaseDiagramMatrix, ScatteringLengths, NumberOfAtoms)
 
 %% Smoothen phase diagram
 load('./Results/Data_Full3D/PhaseDiagramTilted_Theta_20.mat');  % Load M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST
-TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 20";
+TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 20^\circ";
+
 Scripts.plotSmoothedPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, TitleString);
+%% Select boundary points from original phase diagram
+load('./Results/Data_Full3D/PhaseDiagramUntilted.mat')
+M = M;
+N_atoms = round(NUM_ATOMS_LIST * 1E-5, 2);
+scatt_lengths = SCATTERING_LENGTH_RANGE;
+savePath = './Results/Data_Full3D//BoundaryPoints/SelectedPoints_Untilted_1.mat';
+
+Scripts.selectPhaseDiagramPoints(M, N_atoms, scatt_lengths, savePath);
+%% Interactively modify selected boundary points
+load('./Results/Data_Full3D/PhaseDiagramUntilted.mat')
+TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 20^\circ";
+
+Scripts.plotSmoothedContourOnPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, ...
+    TitleString, ...
+    './Results/Data_Full3D/BoundaryPoints/SelectedPoints_Untilted_1.mat', './Results/Data_Full3D/BoundaryPoints/SelectedPoints_Untilted_2.mat', './Results/Data_Full3D/BoundaryPoints/SelectedPoints_Untilted_3.mat', ...
+    'interpMethod', 'pchip', 'showPoints', true);
+%% Edit modified points
+load('./Results/Data_Full3D/PhaseDiagramUntilted.mat')
+TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 20^\circ";
+
+Scripts.plotSmoothedContourOnPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, ...
+    TitleString, ...
+    './Results/Data_Full3D//BoundaryPoints/ModifiedPoints_Untilted_1.mat', './Results/Data_Full3D//BoundaryPoints/ModifiedPoints_Untilted_2.mat', './Results/Data_Full3D//BoundaryPoints/ModifiedPoints_Untilted_3.mat', ...
+    'interpMethod', 'pchip', 'showPoints', true);
+%% Plot with interpolated phase boundary
+load('./Results/Data_Full3D/PhaseDiagramUntilted.mat')
+PhaseDiagramMatrix   = M;
+ScatteringLengths    = SCATTERING_LENGTH_RANGE;
+NumberOfAtoms        = NUM_ATOMS_LIST;
+PhaseBoundary_Untilted = load("./Results/Data_Full3D/BoundaryPoints/PhaseBoundary_Untilted.mat");
+TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 0^\circ";
+
+Scripts.plotPhaseDiagramWithBoundaries(PhaseDiagramMatrix, NumberOfAtoms, ScatteringLengths, PhaseBoundary_Untilted, TitleString);
 
 %% Density modulation determination
 
diff --git a/Dipolar-Gas-Simulator/+Scripts/selectPhaseDiagramPoints.m b/Dipolar-Gas-Simulator/+Scripts/selectPhaseDiagramPoints.m
new file mode 100644
index 0000000..d8037b8
--- /dev/null
+++ b/Dipolar-Gas-Simulator/+Scripts/selectPhaseDiagramPoints.m
@@ -0,0 +1,99 @@
+function selectPhaseDiagramPoints(PhaseDiagramMatrix, NumberOfAtoms, ScatteringLengths, savePath)
+% Interactive point selector for phase diagram.
+% - Single click to select a point (plotted in red)
+% - Double click to remove nearest point
+% - Save points to a MAT file using savePath
+
+    xlen = length(NumberOfAtoms);
+    ylen = length(ScatteringLengths);
+
+    selectedIndices = [];
+    selectedValues = [];
+    plottedPoints = [];  % handles of plotted markers
+
+    % Create figure and axes
+    fig = figure('Name', 'Interactive Phase Diagram', 'Position', [100, 100, 1600, 900]);
+    ax = axes(fig);
+    imagesc(ax, PhaseDiagramMatrix);
+    set(ax, 'YDir', 'normal');
+    colormap(ax, parula);
+    colorbar(ax);
+    axis(ax, 'equal', 'tight');
+    hold(ax, 'on');  % This is crucial!
+
+    % Style
+    set(ax, 'FontSize', 16, 'Box', 'On', 'Linewidth', 2);
+    xticks(ax, 1:xlen);
+    yticks(ax, 1:ylen);
+    xticklabels(ax, string(NumberOfAtoms));
+    yticklabels(ax, string(ScatteringLengths));
+    xlabel(ax, 'Number of Atoms', 'FontSize', 16);
+    ylabel(ax, 'Scattering Length (\times a_o)', 'FontSize', 16);
+    grid(ax, 'on');
+
+    % Save Button
+    uicontrol('Style', 'pushbutton', ...
+              'String', 'Save Points', ...
+              'FontSize', 12, ...
+              'Position', [20 20 120 40], ...
+              'Callback', @savePoints);
+
+    % Set callback on axes, not on image
+    fig.WindowButtonDownFcn = @handleClick;
+    
+    selectedIndices = zeros(0, 2);  % ← ensure it's Nx2, not just []
+    selectedValues  = zeros(0, 2);
+    plottedPoints   = gobjects(0);  % empty graphics array
+
+
+    function handleClick(~, ~)
+        try
+            clickType = get(fig, 'SelectionType');
+            cp = get(ax, 'CurrentPoint');
+            xClick = round(cp(1, 1));
+            yClick = round(cp(1, 2));
+
+            if xClick < 1 || xClick > xlen || yClick < 1 || yClick > ylen
+                return;
+            end
+
+            if strcmp(clickType, 'normal')  % Single click = add point
+                alreadySelected = any(ismember(selectedIndices, [xClick, yClick], 'rows'));
+                if ~alreadySelected
+                    selectedIndices(end+1, :) = [xClick, yClick];
+                    selectedValues(end+1, :) = [NumberOfAtoms(xClick), ScatteringLengths(yClick)];
+                    h = plot(ax, xClick, yClick, 'r*', 'MarkerSize', 10, 'LineWidth', 1.5);
+                    plottedPoints(end+1) = h;
+                    fprintf('Selected: N = %.2f, a = %.2f\n', NumberOfAtoms(xClick), ScatteringLengths(yClick));
+                end
+
+            elseif strcmp(clickType, 'open')  % Double click = remove point
+                if isempty(selectedIndices), return; end
+                dists = vecnorm(selectedIndices - [xClick, yClick], 2, 2);
+                [minDist, idx] = min(dists);
+                if minDist <= 1
+                    fprintf('Removed: N = %.2f, a = %.2f\n', selectedValues(idx,1), selectedValues(idx,2));
+                    delete(plottedPoints(idx));
+                    selectedIndices(idx, :) = [];
+                    selectedValues(idx, :) = [];
+                    plottedPoints(idx) = [];
+                end
+            end
+
+        catch ME
+            fprintf('Error in handleClick:\n%s\n', getReport(ME, 'extended'));
+        end
+    end
+
+
+    function savePoints(~, ~)
+        if isempty(selectedValues)
+            warndlg('No points selected.', 'Save Warning');
+            return;
+        end
+        SelectedPoints.indices = selectedIndices;
+        SelectedPoints.values = selectedValues;
+        save(savePath, 'SelectedPoints');
+        msgbox(sprintf('Saved %d points to:\n%s', size(selectedValues,1), savePath), 'Saved');
+    end
+end