From 40915f39138c2f838ec2dcede660094493976d11 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Thu, 15 May 2025 17:57:09 +0200 Subject: [PATCH] Added scripts to edit phase diagram matrix created manually using createPhaseDiagram function and another function to plot with colours, contours to define the boundaries between phases. --- .../+Scripts/createPhaseDiagram.m | 11 +- .../+Scripts/editPhaseDiagram.m | 74 +++++++++++ .../+Scripts/plotSmoothedPhaseDiagram.m | 39 ++++++ Dipolar-Gas-Simulator/+Scripts/run_locally.m | 115 +++++------------- 4 files changed, 150 insertions(+), 89 deletions(-) create mode 100644 Dipolar-Gas-Simulator/+Scripts/editPhaseDiagram.m create mode 100644 Dipolar-Gas-Simulator/+Scripts/plotSmoothedPhaseDiagram.m diff --git a/Dipolar-Gas-Simulator/+Scripts/createPhaseDiagram.m b/Dipolar-Gas-Simulator/+Scripts/createPhaseDiagram.m index 0f5e796..d6af174 100644 --- a/Dipolar-Gas-Simulator/+Scripts/createPhaseDiagram.m +++ b/Dipolar-Gas-Simulator/+Scripts/createPhaseDiagram.m @@ -1,17 +1,20 @@ function createPhaseDiagram() % Define axis values - SCATTERING_LENGTH_RANGE = [80.00 81.04 82.08 83.12 84.17 85.21 86.25 87.29 88.33 89.38 90.42 91.46 92.50 93.54 94.58 95.62 96.67 97.71 98.75 99.79 100.83 101.88 102.92 103.96 105.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]; + SCATTERING_LENGTH_RANGE = [75.00 76.09 77.18 78.27 79.36 80.00 81.04 82.08 83.12 84.17 85.21 86.25 87.29 88.33 89.38 90.42 91.46 92.50 93.54 94.58 95.62 96.67 97.71 98.75 99.79 100.83 101.88 102.92 103.96 105.00]; + NUM_ATOMS_LIST = [50000 54545 59091 63636 68182 72727 77273 81818 86364 90909 95455 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]; xlen = length(NUM_ATOMS_LIST); ylen = length(SCATTERING_LENGTH_RANGE); - M = zeros(xlen,ylen); % Phase diagram matrix + M = zeros(ylen,xlen); % Phase diagram matrix fig = figure('Name', 'Manual Phase Diagram Input', ... 'NumberTitle', 'off', ... 'KeyPressFcn', @keyCallback); - + fig = figure(1); + fig.WindowState = 'maximized'; + clf + set(gcf,'Position', [100, 100, 1600, 900]) hImg = imagesc(M); set(gca, 'YDir', 'normal'); colormap(parula); diff --git a/Dipolar-Gas-Simulator/+Scripts/editPhaseDiagram.m b/Dipolar-Gas-Simulator/+Scripts/editPhaseDiagram.m new file mode 100644 index 0000000..8f26492 --- /dev/null +++ b/Dipolar-Gas-Simulator/+Scripts/editPhaseDiagram.m @@ -0,0 +1,74 @@ +function editPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST) + % Validate dimensions + [ylen, xlen] = size(M); + if length(NUM_ATOMS_LIST) ~= xlen || length(SCATTERING_LENGTH_RANGE) ~= ylen + error('Matrix dimensions do not match axis lengths.'); + end + + % Create figure + fig = figure('Name', 'Manual Phase Diagram Editor', ... + 'NumberTitle', 'off', ... + 'KeyPressFcn', @keyCallback); + fig.WindowState = 'maximized'; + clf + set(gcf,'Position', [100, 100, 1600, 900]) + + hImg = imagesc(M); + 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)); + xlabel('Number of Atoms'); + ylabel('Scattering Length (a₀)'); + grid on; + title('Click a cell to modify value. Press "q" to quit.'); + + set(fig, 'WindowButtonDownFcn', @clickCallback); + uiwait(fig); % Wait until user quits + + % Prompt for saving + [filename, pathname] = uiputfile('edited_phase_diagram.mat', 'Save Edited Matrix As'); + if ischar(filename) + save(fullfile(pathname, filename), 'M', 'SCATTERING_LENGTH_RANGE', 'NUM_ATOMS_LIST'); + disp(['Matrix saved to ' fullfile(pathname, filename)]); + end + + % Callback: Mouse click to edit value + function clickCallback(~, ~) + C = get(gca, 'CurrentPoint'); + col = round(C(1,1)); + row = round(C(1,2)); + if row >= 1 && row <= ylen && col >= 1 && col <= xlen + prompt = sprintf('Enter new value for a = %.2f, N = %d:', ... + SCATTERING_LENGTH_RANGE(row), NUM_ATOMS_LIST(col)); + answer = inputdlg(prompt, 'Edit Cell Value', 1, {num2str(M(row,col))}); + if ~isempty(answer) + val = str2double(answer{1}); + if ~isnan(val) + M(row, col) = val; + refreshDisplay(); + else + warndlg('Invalid input. Please enter a numeric value.', 'Invalid Input'); + end + end + end + end + + % Callback: Quit on 'q' + function keyCallback(~, event) + if strcmp(event.Key, 'q') + uiresume(fig); + close(fig); + end + end + + % Update image + function refreshDisplay() + set(hImg, 'CData', M); + drawnow; + end +end diff --git a/Dipolar-Gas-Simulator/+Scripts/plotSmoothedPhaseDiagram.m b/Dipolar-Gas-Simulator/+Scripts/plotSmoothedPhaseDiagram.m new file mode 100644 index 0000000..72d8761 --- /dev/null +++ b/Dipolar-Gas-Simulator/+Scripts/plotSmoothedPhaseDiagram.m @@ -0,0 +1,39 @@ +function plotSmoothedPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, titleString) + % Region labels and corresponding colors + regionNames = ["Unmodulated", "SSD", "Stripe", "Labyrinth", "Honeycomb"]; + regionColors = [ + 0.8 0.8 0.8; % unmodulated - light gray + 0.2 0.6 1.0; % SSD - blue + 0.2 0.8 0.2; % Stripe - green + 1.0 0.6 0.2; % Labyrinth - orange + 0.8 0.2 0.8 % Honeycomb - purple + ]; + + % Interpolate the discrete-valued matrix using 'nearest' for smoothing + [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'); + + % Plot + figure('Color','w'); % Correctly assign figure handle + clf + set(gcf, 'Position', [100, 100, 950, 800]); + imagesc([min(NUM_ATOMS_LIST), max(NUM_ATOMS_LIST)], ... + [min(SCATTERING_LENGTH_RANGE), max(SCATTERING_LENGTH_RANGE)], ... + Mq); + % Overlay smooth contour lines to show boundaries + hold on; + 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); + 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); + 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 09c0ff4..91ed875 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -553,9 +553,10 @@ SaveDirectory = './Results/Data_3D/PhaseDiagram/GradientDescent/'; JobNumber = 0; Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) %% -SaveDirectory = './Results/Data_3D/PhaseDiagram/ImagTimePropagation/'; +SaveDirectory = './Results/Data_3D/PhaseDiagram/ImagTimePropagation/aS_8.833000e+01_theta_000_phi_000_N_81818'; JobNumber = 0; Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) + %% Generate lists % Set display format @@ -581,9 +582,12 @@ disp(row2) SCATTERING_LENGTH_RANGE = "[80.00 81.04 82.08 83.12 84.17 85.21 86.25 87.29 88.33 89.38 90.42 91.46 92.50 93.54 94.58 95.62 96.67 97.71 98.75 99.79 100.83 101.88 102.92 103.96 105.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]"; -%% Explore the phase diagram -SCATTERING_LENGTH_RANGE = [91.46 92.50 93.54 94.58 95.62 96.67 97.71]; -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]; +SCATTERING_LENGTH_RANGE ="[75.00 76.09 77.18 78.27 79.36 80.00 81.04 82.08 83.12 84.17 85.21 86.25 87.29]"; +NUM_ATOMS_LIST ="[50000 54545 59091 63636 68182 72727 77273 81818 86364 90909 95455]"; + +%% Explore the states +SCATTERING_LENGTH_RANGE = []; +NUM_ATOMS_LIST = []; JobNumber = 0; @@ -592,7 +596,7 @@ for i = 1:length(SCATTERING_LENGTH_RANGE) for j = 1:length(NUM_ATOMS_LIST) N = NUM_ATOMS_LIST(j); - SaveDirectory = sprintf('D:/Results-Numerics/PhaseDiagram/ImagTimePropagation/aS_%.6e_theta_000_phi_000_N_%d', aS, N); + SaveDirectory = sprintf('./aS_%.6e_theta_000_phi_000_N_%d', aS, N); fprintf('Processing JobNumber %d: %s\n', JobNumber, SaveDirectory); % Call the plotting function @@ -608,88 +612,14 @@ for i = 1:length(SCATTERING_LENGTH_RANGE) end end -%% Parameters that need to be run again -% Cell array of the input strings -strs = { - 'aS_9.875000e+01_theta_000_phi_000_N_2958333' - 'aS_9.562000e+01_theta_000_phi_000_N_4183333' - 'aS_9.562000e+01_theta_000_phi_000_N_1325000' - 'aS_9.458000e+01_theta_000_phi_000_N_508333' - 'aS_9.458000e+01_theta_000_phi_000_N_5000000' - 'aS_9.458000e+01_theta_000_phi_000_N_1937500' - 'aS_9.354000e+01_theta_000_phi_000_N_4795833' - 'aS_9.250000e+01_theta_000_phi_000_N_5000000' - 'aS_9.250000e+01_theta_000_phi_000_N_4183333' - 'aS_9.146000e+01_theta_000_phi_000_N_1937500' - 'aS_8.833000e+01_theta_000_phi_000_N_2754167' - 'aS_8.833000e+01_theta_000_phi_000_N_1325000' - 'aS_8.729000e+01_theta_000_phi_000_N_3775000' - 'aS_8.625000e+01_theta_000_phi_000_N_712500' - 'aS_8.625000e+01_theta_000_phi_000_N_4795833' - 'aS_8.625000e+01_theta_000_phi_000_N_3979167' - 'aS_8.521000e+01_theta_000_phi_000_N_1937500' - 'aS_8.208000e+01_theta_000_phi_000_N_3979167' - 'aS_105_theta_000_phi_000_N_916667' - 'aS_1.039600e+02_theta_000_phi_000_N_3366667' - 'aS_1.008300e+02_theta_000_phi_000_N_3979167' - 'aS_080_theta_000_phi_000_N_5000000' -}; - -% Extract values using regex -tokens = regexp(strs, 'aS_([0-9.e+-]+)_theta_000_phi_000_N_([0-9]+)', 'tokens'); - -% Convert to numeric array -pairs = cellfun(@(x) [str2double(x{1}), str2double(x{2})], [tokens{:}], 'UniformOutput', false); -pairs = vertcat(pairs{:}); - -% Sort by a_s (first column) -pairs = sortrows(pairs, 1); - -% Display as a table -disp(array2table(pairs, 'VariableNames', {'a_s', 'N'})); - -%% Phase diagram for untilted case -N = [1E5, 3.04E5, 5.08E5, 7.125E5, 9.16E5, 1.12E6, 1.325E6, 1.529E6, ... - 1.733E6, 1.9375E6, 2.141E6, 2.345E6, 2.55E6, 2.75E6, 2.95E6, ... - 3.1625E6, 3.367E6, 3.57E6, 3.775E6, 3.979E6, 4.183E6, 4.3875E6, ... - 4.591E6, 4.795E6, 5E6]; - -as_UB = [88.33, 94.58, 95.62, 96.67, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, ... - 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71]; - -as_LB = [87.29, 93.54, 94.58, 95.62, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, ... - 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67]; - -% Filter only rows with non-NaN data -valid_idx = ~isnan(as_LB) & ~isnan(as_UB); -N_valid = N(valid_idx); -LB_valid = as_LB(valid_idx); -UB_valid = as_UB(valid_idx); - -% Create shaded area between LB and UB -x_fill = [N_valid, fliplr(N_valid)]; -y_fill = [LB_valid, fliplr(UB_valid)]; - -% Plot settings -figure(1); -set(gcf,'Position',[100 100 950 750]) -fill(x_fill, y_fill, [0.8 0.8 1], 'EdgeColor', 'none'); % Light blue shade -% Axes settings -set(gca, 'XScale', 'log'); -xlim([1E4, 1E7]); -ylim([79, 106]); -xlabel('Atom number', 'Interpreter', 'latex', 'FontSize', 16); -ylabel('Scattering length $a_s$ ($a_0$)', 'Interpreter', 'latex', 'FontSize', 16); -grid on; -set(gca,'FontSize',16,'Box','On','Linewidth',2); - %% Visualize phase diagram -load('./Results/Data_3D/GradientDescent/phase_diagram_matrix_theta_30.mat') +% load('./Results/phase_diagram_matrix.mat') +load('./Results/edited_phase_diagram.mat') PhaseDiagramMatrix = M; ScatteringLengths = SCATTERING_LENGTH_RANGE; -NumberOfAtoms = NUM_ATOMS_LIST * 1E-5; -xlen = length(ScatteringLengths); -ylen = length(NumberOfAtoms); +NumberOfAtoms = round(NUM_ATOMS_LIST * 1E-5,2); +xlen = length(NumberOfAtoms); +ylen = length(ScatteringLengths); figure(1); set(gcf,'Position',[100 100 950 750]) @@ -704,11 +634,24 @@ xticks(1:xlen); yticks(1:ylen); xticklabels(string(NumberOfAtoms)); yticklabels(string(ScatteringLengths)); -xlabel('Number of Atoms (x 10^5)', 'Interpreter', 'tex', FontSize=16); +xlabel('Number of Atoms (x 1E5)', FontSize=16); ylabel('Scattering Length (a0)', FontSize=16); % title('Zero-temperature Phase Diagram for \theta = 0', 'Interpreter', 'tex', FontSize=16); grid on; +%% Edit phase diagram +% load('./Results/phase_diagram_matrix.mat') +load('./Results/edited_phase_diagram.mat') +PhaseDiagramMatrix = M; +ScatteringLengths = SCATTERING_LENGTH_RANGE; +NumberOfAtoms = NUM_ATOMS_LIST; +Scripts.editPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST) + +%% Smoothen phase diagram +load('./Results/edited_phase_diagram.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"; +Scripts.plotSmoothedPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, titleString); + %% Density modulation determination % SaveDirectory = './Results/Data_3D/GradientDescent/Phi020/aS_090_theta_020_phi_000_N_100000'; @@ -742,3 +685,5 @@ if ModulationFlag else disp('The state is not modulated'); end + +