karthik/Calculations
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity

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.

Browse Source
This commit is contained in:
Karthik 2025-05-15 17:57:09 +02:00
parent 5a8b117e3e
commit 40915f3913
4 changed files with 150 additions and 89 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
Dipolar-Gas-Simulator/+Scripts/createPhaseDiagram.m
View File

@ -1,17 +1,20 @@
function createPhaseDiagram() function createPhaseDiagram()
% Define axis values % 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]; 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 = [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 = [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); xlen = length(NUM_ATOMS_LIST);
ylen = length(SCATTERING_LENGTH_RANGE); 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', ... fig = figure('Name', 'Manual Phase Diagram Input', ...
'NumberTitle', 'off', ... 'NumberTitle', 'off', ...
'KeyPressFcn', @keyCallback); 'KeyPressFcn', @keyCallback);
fig = figure(1);
fig.WindowState = 'maximized';
clf
set(gcf,'Position', [100, 100, 1600, 900])
hImg = imagesc(M); hImg = imagesc(M);
set(gca, 'YDir', 'normal'); set(gca, 'YDir', 'normal');
colormap(parula); colormap(parula);

74
Dipolar-Gas-Simulator/+Scripts/editPhaseDiagram.m Normal file
View File

@ -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

39
Dipolar-Gas-Simulator/+Scripts/plotSmoothedPhaseDiagram.m Normal file
View File

@ -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

115
Dipolar-Gas-Simulator/+Scripts/run_locally.m
View File

@ -553,9 +553,10 @@ SaveDirectory = './Results/Data_3D/PhaseDiagram/GradientDescent/';
JobNumber = 0; JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) 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; JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%% Generate lists %% Generate lists
% Set display format % 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]"; 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]"; 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 ="[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]";
SCATTERING_LENGTH_RANGE = [91.46 92.50 93.54 94.58 95.62 96.67 97.71]; NUM_ATOMS_LIST ="[50000 54545 59091 63636 68182 72727 77273 81818 86364 90909 95455]";
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 states
SCATTERING_LENGTH_RANGE = [];
NUM_ATOMS_LIST = [];
JobNumber = 0; JobNumber = 0;
@ -592,7 +596,7 @@ for i = 1:length(SCATTERING_LENGTH_RANGE)
for j = 1:length(NUM_ATOMS_LIST) for j = 1:length(NUM_ATOMS_LIST)
N = NUM_ATOMS_LIST(j); 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); fprintf('Processing JobNumber %d: %s\n', JobNumber, SaveDirectory);
% Call the plotting function % Call the plotting function
@ -608,88 +612,14 @@ for i = 1:length(SCATTERING_LENGTH_RANGE)
end end
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 %% 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; PhaseDiagramMatrix = M;
ScatteringLengths = SCATTERING_LENGTH_RANGE; ScatteringLengths = SCATTERING_LENGTH_RANGE;
NumberOfAtoms = NUM_ATOMS_LIST * 1E-5; NumberOfAtoms = round(NUM_ATOMS_LIST * 1E-5,2);
xlen = length(ScatteringLengths); xlen = length(NumberOfAtoms);
ylen = length(NumberOfAtoms); ylen = length(ScatteringLengths);
figure(1); figure(1);
set(gcf,'Position',[100 100 950 750]) set(gcf,'Position',[100 100 950 750])
@ -704,11 +634,24 @@ xticks(1:xlen);
yticks(1:ylen); yticks(1:ylen);
xticklabels(string(NumberOfAtoms)); xticklabels(string(NumberOfAtoms));
yticklabels(string(ScatteringLengths)); 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); ylabel('Scattering Length (a0)', FontSize=16);
% title('Zero-temperature Phase Diagram for \theta = 0', 'Interpreter', 'tex', FontSize=16); % title('Zero-temperature Phase Diagram for \theta = 0', 'Interpreter', 'tex', FontSize=16);
grid on; 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 %% Density modulation determination
% SaveDirectory = './Results/Data_3D/GradientDescent/Phi020/aS_090_theta_020_phi_000_N_100000'; % SaveDirectory = './Results/Data_3D/GradientDescent/Phi020/aS_090_theta_020_phi_000_N_100000';
@ -742,3 +685,5 @@ if ModulationFlag
else else
disp('The state is not modulated'); disp('The state is not modulated');
end end