diff --git a/Dipolar-Gas-Simulator/+Scripts/plotSmoothedPhaseDiagram.m b/Dipolar-Gas-Simulator/+Scripts/plotSmoothedPhaseDiagram.m
index 72d8761..ee8bb98 100644
--- a/Dipolar-Gas-Simulator/+Scripts/plotSmoothedPhaseDiagram.m
+++ b/Dipolar-Gas-Simulator/+Scripts/plotSmoothedPhaseDiagram.m
@@ -19,6 +19,7 @@ function plotSmoothedPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, ti
     figure('Color','w');  % Correctly assign figure handle
     clf
     set(gcf, 'Position', [100, 100, 950, 800]);
+    set(gcf, 'Renderer', 'opengl');
     imagesc([min(NUM_ATOMS_LIST), max(NUM_ATOMS_LIST)], ...
             [min(SCATTERING_LENGTH_RANGE), max(SCATTERING_LENGTH_RANGE)], ...
             Mq);
@@ -27,13 +28,18 @@ function plotSmoothedPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, ti
     contour(Xq, Yq, Mq, 0.5:1:3.5, 'k', 'LineWidth', 1.2);  % 0.5 to 3.5 separates 0–1, 1–2, etc.
     set(gca, 'YDir', 'normal');  % Keep Y increasing upward
     colormap(regionColors);
-    colorbar('Ticks', 0:4, 'TickLabels', regionNames, 'FontSize', 12);
+    cb = colorbar('Ticks', 0:4, 'TickLabels', regionNames, 'FontSize', 12);
+    cb.Color = 'k';  % set colorbar text color
     clim([0 4]);  % Match value range to colormap rows
     % Axis formatting
     xlabel('Number of Atoms', 'Interpreter', 'tex', 'FontSize', 16);
     ylabel('Scattering Length (\times a_o)', 'Interpreter', 'tex', 'FontSize', 16);
-    title(titleString, 'Interpreter', 'tex', 'FontSize', 18);
-    set(gca, 'FontSize', 16);
+    t = title(titleString, 'Interpreter', 'tex', 'FontSize', 18);
+    t.Color = 'k';  % force title color to black
+
+    set(gca, 'FontSize', 16, ...
+        'XColor', 'k', 'YColor', 'k', ...
+        'Color', 'none');  % Ensure transparent axes background
     axis tight;
     grid on;
 end
\ No newline at end of file
diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m
index 428fcd0..cf0321f 100644
--- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m
+++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m
@@ -661,8 +661,9 @@ SCATTERING_LENGTH_RANGE ="[75.00 76.09 77.18 78.27 79.36 80.00 81.04 82.08 83.12
 NUM_ATOMS_LIST          ="[50000 54545 59091 63636 68182 72727 77273 81818 86364 90909 95455]";
 
 %% Explore the states
-SCATTERING_LENGTH_RANGE = [];
-NUM_ATOMS_LIST          = [];
+SCATTERING_LENGTH_RANGE = [80.00];
+NUM_ATOMS_LIST          = [100000 304167 508333 712500 916667 1120833 1325000 1529167 1733333 1937500 2141667 2345833 2550000 2754167 2958333 3162500 3366667 3570833 3775000 3979167 4183333 4387500 4591667 4795833 5000000];
+% NUM_ATOMS_LIST          = [4183333 4387500 4591667 4795833 5000000];
 
 JobNumber = 0;
 
@@ -671,7 +672,7 @@ for i = 1:length(SCATTERING_LENGTH_RANGE)
     for j = 1:length(NUM_ATOMS_LIST)
         N = NUM_ATOMS_LIST(j);
         
-        SaveDirectory = sprintf('./aS_%.6e_theta_000_phi_000_N_%d', aS, N);
+        SaveDirectory = sprintf('D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta20/HighN/aS_%.6e_theta_020_phi_000_N_%d', aS, N);
         fprintf('Processing JobNumber %d: %s\n', JobNumber, SaveDirectory);
 
         % Call the plotting function
@@ -718,7 +719,7 @@ T = table(aS_values, N_values, 'VariableNames', {'aS', 'N'});
 disp(T);
 
 %% Visualize phase diagram
-load('./Results/Data_Full3D/PhaseDiagramUntilted.mat')
+load('./Results/Data_Full3D/PhaseDiagramTilted_Theta_20.mat')
 PhaseDiagramMatrix   = M;
 ScatteringLengths    = SCATTERING_LENGTH_RANGE;
 NumberOfAtoms        = round(NUM_ATOMS_LIST * 1E-5,2);
@@ -731,29 +732,29 @@ fig.WindowState = 'maximized';
 clf
 set(gcf,'Position', [100, 100, 1600, 900])
 set(gca,'FontSize',16,'Box','On','Linewidth',2);
-hImg = imagesc(M);
+hImg = imagesc(PhaseDiagramMatrix);
 set(gca, 'YDir', 'normal');
 colormap(parula);
 colorbar;
 axis equal tight;
 xticks(1:xlen);
 yticks(1:ylen);
-xticklabels(string(NUM_ATOMS_LIST));
-yticklabels(string(SCATTERING_LENGTH_RANGE));
+xticklabels(string(NumberOfAtoms));
+yticklabels(string(ScatteringLengths));
 xlabel('Number of Atoms', 'Interpreter', 'tex', 'FontSize', 16);
 ylabel('Scattering Length (\times a_o)', 'Interpreter', 'tex', 'FontSize', 16);
 grid on;
 
 %% Edit phase diagram
-load('./Results/Data_Full3D/PhaseDiagramUntilted.mat')
+load('./Results/Data_Full3D/PhaseDiagramTilted_Theta_20.mat')
 PhaseDiagramMatrix   = M;
 ScatteringLengths    = SCATTERING_LENGTH_RANGE;
 NumberOfAtoms        = NUM_ATOMS_LIST;
-Scripts.editPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST)
+Scripts.editPhaseDiagram(PhaseDiagramMatrix, ScatteringLengths, NumberOfAtoms)
 
 %% Smoothen phase diagram
-load('./Results/Data_Full3D/PhaseDiagramUntilted.mat');  % Load M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST
-TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 0";
+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";
 Scripts.plotSmoothedPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, TitleString);
 
 %% Density modulation determination