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

Added script to identify and count droplets in the results of the numerics.

Browse Source
This commit is contained in:
Karthik 2025-05-16 18:44:04 +02:00
parent 883a5bc453
commit 8ae97eb267
2 changed files with 193 additions and 34 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

84
Dipolar-Gas-Simulator/+Scripts/identifyAndCountDroplets.m Normal file
View File

@ -0,0 +1,84 @@
function [DropletNumber] = identifyAndCountDroplets(folder_path, run_index, radius, minPeakHeight, SuppressPlotFlag)
set(0,'defaulttextInterpreter','latex')
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
format long
% Load data
Data = load(sprintf(horzcat(folder_path, '/Run_%03i/psi_gs.mat'),run_index),'psi','Params','Transf','Observ');
Params = Data.Params;
Transf = Data.Transf;
Observ = Data.Observ;
if isgpuarray(Data.psi)
psi = gather(Data.psi);
else
psi = Data.psi;
end
if isgpuarray(Data.Observ.residual)
Observ.residual = gather(Data.Observ.residual);
else
Observ.residual = Data.Observ.residual;
end
% Axes scaling and coordinates in micrometers
x = Transf.x * Params.l0 * 1e6;
y = Transf.y * Params.l0 * 1e6;
z = Transf.z * Params.l0 * 1e6;
dz = z(2)-z(1);
% Compute probability density |psi|^2
n = abs(psi).^2;
nxy = squeeze(trapz(n*dz,3));
state = nxy';
peaks = extractDroplets(state, radius, minPeakHeight);
DropletNumber = numel(peaks(peaks>0));
[row, col] = find(peaks);
if ~SuppressPlotFlag
figure(1);
clf
set(gcf,'Position', [100, 100, 900, 900])
plotxy = pcolor(x,y,state);
set(plotxy, 'EdgeColor', 'none');
hold on;
plot(x(col), y(row), 'w+', 'MarkerSize', 8, 'LineWidth', 1.5);
cbar1 = colorbar;
cbar1.Label.Interpreter = 'latex';
% cbar1.Ticks = []; % Disable the ticks
colormap(gca, Helper.Colormaps.plasma())
xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16)
ylabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16)
title(['N_d = ' num2str(DropletNumber)], ...
'Interpreter', 'tex', 'FontSize', 16)
set(gca, 'FontSize', 16); % For tick labels only
end
end
function [peaks] = extractDroplets(image, radius, minPeakHeight)
peaks = imregionalmax(image);
peaks(image < minPeakHeight) = 0;
[row, col] = find(peaks); % row and col are the coordinates of the detected peaks
unique_peaks = true(size(row)); % Assume all peaks are unique initially
% Check distances between each peak and remove duplicates
for i = 1:length(row)
if unique_peaks(i) % If this peak hasn't been excluded
% Find the distance from this peak to all other peaks
dist = sqrt((row(i) - row).^2 + (col(i) - col).^2); % Euclidean distance
% Exclude any peak within the radius
unique_peaks(dist < radius & dist > 0) = false; % Exclude other peaks in radius
end
end
% Save the peaks
peaks_mask = false(size(image));
peaks_mask(sub2ind(size(image), row(unique_peaks), col(unique_peaks))) = true;
peaks = peaks_mask; % Store the unique peaks mask
end

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

@ -33,7 +33,7 @@ OptionsStruct.PlotLive = true;
OptionsStruct.JobNumber = 0;
OptionsStruct.RunOnGPU = false;
OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = './Results/Data_3D/GradientDescent'; % './Results/Data_3D/AnisotropicTrap/Tilted0'
OptionsStruct.SaveDirectory = './Results/Data_Full3D/GradientDescent'; % './Results/Data_Full3D/AnisotropicTrap/Tilted0'
options = Helper.convertstruct2cell(OptionsStruct);
sim = Simulator.DipolarGas(options{:});
@ -85,7 +85,7 @@ OptionsStruct.PlotLive = true;
OptionsStruct.JobNumber = 0;
OptionsStruct.RunOnGPU = false;
OptionsStruct.SaveData = false;
OptionsStruct.SaveDirectory = './Results/Data_3D/RealTimeDynamics';
OptionsStruct.SaveDirectory = './Results/Data_Full3D/RealTimeDynamics';
options = Helper.convertstruct2cell(OptionsStruct);
sim = Simulator.DipolarGas(options{:});
@ -484,78 +484,153 @@ JobNumber = 2; % 82
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction_in_plane_trap(SaveDirectory, JobNumber);
%% - Plot GS wavefunction
% SaveDirectory = './Results/Data_3D/CompleteLHY/AspectRatio2_8';
% SaveDirectory = './Results/Data_3D/CompleteLHY/AspectRatio3_7';
SaveDirectory = './Results/Data_3D/CompleteLHY/AspectRatio3_8';
% SaveDirectory = './Results/Data_Full3D/CompleteLHY/AspectRatio2_8';
% SaveDirectory = './Results/Data_Full3D/CompleteLHY/AspectRatio3_7';
SaveDirectory = './Results/Data_Full3D/CompleteLHY/AspectRatio3_8';
JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
% SaveDirectory = './Results/Data_3D/CompleteLHY/BeyondSSD_SSD';
% SaveDirectory = './Results/Data_3D/CompleteLHY/BeyondSSD_Labyrinth';
SaveDirectory = './Results/Data_3D/CompleteLHY/BeyondSSD_Honeycomb';
% SaveDirectory = './Results/Data_Full3D/CompleteLHY/BeyondSSD_SSD';
% SaveDirectory = './Results/Data_Full3D/CompleteLHY/BeyondSSD_Labyrinth';
SaveDirectory = './Results/Data_Full3D/CompleteLHY/BeyondSSD_Honeycomb';
JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%% - Plot GS wavefunction
SaveDirectory = './Results/Data_3D/ApproximateLHY/AspectRatio2_8';
% SaveDirectory = './Results/Data_3D/ApproximateLHY/AspectRatio3_7';
% SaveDirectory = './Results/Data_3D/ApproximateLHY/AspectRatio3_8';
% SaveDirectory = './Results/Data_3D/ApproximateLHY/AspectRatio3_9';
SaveDirectory = './Results/Data_Full3D/ApproximateLHY/AspectRatio2_8';
% SaveDirectory = './Results/Data_Full3D/ApproximateLHY/AspectRatio3_7';
% SaveDirectory = './Results/Data_Full3D/ApproximateLHY/AspectRatio3_8';
% SaveDirectory = './Results/Data_Full3D/ApproximateLHY/AspectRatio3_9';
JobNumber = 3;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
% SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_SSD';
% SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_Labyrinth';
SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_Honeycomb';
% SaveDirectory = './Results/Data_Full3D/ApproximateLHY/BeyondSSD_SSD';
% SaveDirectory = './Results/Data_Full3D/ApproximateLHY/BeyondSSD_Labyrinth';
SaveDirectory = './Results/Data_Full3D/ApproximateLHY/BeyondSSD_Honeycomb';
JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
SaveDirectory = './Results/Data_3D/TiltedDipoles0';
SaveDirectory = './Results/Data_Full3D/TiltedDipoles0';
JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
SaveDirectory = './Results/Data_3D/TiltedDipoles15';
SaveDirectory = './Results/Data_Full3D/TiltedDipoles15';
JobNumber = 2;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
SaveDirectory = './Results/Data_3D/TiltedDipoles30';
SaveDirectory = './Results/Data_Full3D/TiltedDipoles30';
JobNumber = 2;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
SaveDirectory = './Results/Data_3D/TiltedDipoles45';
SaveDirectory = './Results/Data_Full3D/TiltedDipoles45';
JobNumber = 2;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
SaveDirectory = './Results/Data_3D/RealTimeDynamics';
SaveDirectory = './Results/Data_Full3D/RealTimeDynamics';
JobNumber = 0;
Plotter.makeMovie(SaveDirectory, JobNumber)
%%
SaveDirectory = './Results/Data_3D/AnisotropicTrap/TiltedDipoles15';
SaveDirectory = './Results/Data_Full3D/AnisotropicTrap/TiltedDipoles15';
JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
SaveDirectory = './Results/Data_3D/AnisotropicTrap/TiltedDipoles30';
SaveDirectory = './Results/Data_Full3D/AnisotropicTrap/TiltedDipoles30';
JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
SaveDirectory = './Results/Data_3D/AnisotropicTrap/TiltedDipoles35';
SaveDirectory = './Results/Data_Full3D/AnisotropicTrap/TiltedDipoles35';
JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
SaveDirectory = './Results/Data_3D/AnisotropicTrap/TiltedDipoles40';
SaveDirectory = './Results/Data_Full3D/AnisotropicTrap/TiltedDipoles40';
JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
SaveDirectory = './Results/Data_3D/AnisotropicTrap/TiltedDipoles45';
SaveDirectory = './Results/Data_Full3D/AnisotropicTrap/TiltedDipoles45';
JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
SaveDirectory = './Results/Data_3D/PhaseDiagram/GradientDescent/';
SaveDirectory = './Results/Data_Full3D/PhaseDiagram/GradientDescent/';
JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
SaveDirectory = './Results/Data_3D/PhaseDiagram/ImagTimePropagation/aS_8.833000e+01_theta_000_phi_000_N_81818';
SaveDirectory = './Results/Data_Full3D/PhaseDiagram/ImagTimePropagation/aS_8.312000e+01_theta_020_phi_000_N_100000';
JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
SaveDirectory = 'D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_8.729000e+01_theta_000_phi_000_N_100000';
JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%% Identify and count droplets
Radius = 2; % The radius within which peaks will be considered duplicates
PeakThreshold = 6E3;
SaveDirectory = 'D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_8.729000e+01_theta_000_phi_000_N_2550000';
JobNumber = 0;
SuppressPlotFlag = false;
DropletNumber = Scripts.identifyAndCountDroplets(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag);
%% Plot number of droplets
% Parameters
Radius = 2;
PeakThreshold = 6E3;
JobNumber = 0;
TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 0";
SuppressPlotFlag = true;
SCATTERING_LENGTH_RANGE = [85.21 86.25 87.29 88.33 89.38];
NUM_ATOMS_LIST = [100000 304167 508333 712500 916667 1120833 1325000 ...
1529167 1733333 1937500 2141667 2345833 2550000 ...
2754167 2958333 3162500 3366667 3570833];
% Prepare figure
figure(1);
clf;
set(gcf,'Position', [100, 100, 1000, 700]);
hold on;
% Color order for better visibility
colors = lines(length(SCATTERING_LENGTH_RANGE));
% Store legend labels
legendEntries = cell(1, length(SCATTERING_LENGTH_RANGE));
% Loop over scattering lengths
for j = 1:length(SCATTERING_LENGTH_RANGE)
aS = SCATTERING_LENGTH_RANGE(j);
% Format scattering length in scientific notation with 6 decimal places
aS_string = sprintf('%.6e', aS);
% Construct base directory for this aS
baseDir = ['D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_' ...
aS_string '_theta_000_phi_000_N_'];
% Preallocate results for this curve
DropletNumbers = zeros(size(NUM_ATOMS_LIST));
% Loop over all atom numbers
for i = 1:length(NUM_ATOMS_LIST)
N = NUM_ATOMS_LIST(i);
SaveDirectory = [baseDir num2str(N)];
% Call your droplet counting function
DropletNumbers(i) = Scripts.identifyAndCountDroplets(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag);
end
% Plot curve
plot(NUM_ATOMS_LIST, DropletNumbers, 'o-', ...
'Color', colors(j, :), 'LineWidth', 1.5);
% Store legend entry
legendEntries{j} = ['a_s = ' num2str(aS) 'a_o'];
end
% Finalize plot
xlabel('Number of Atoms', 'Interpreter', 'tex', 'FontSize', 16);
ylabel('Number of Droplets', 'Interpreter', 'tex', 'FontSize', 16);
title(TitleString, 'Interpreter', 'tex', 'FontSize', 18);
legend(legendEntries, 'Interpreter', 'tex', 'FontSize', 12, 'Location', 'bestoutside');
grid on;
hold off;
%% Generate lists
@ -613,7 +688,7 @@ for i = 1:length(SCATTERING_LENGTH_RANGE)
end
%% Visualize phase diagram
load('./Results/PhaseDiagramUntilted.mat')
load('./Results/Data_Full3D/PhaseDiagramUntilted.mat')
PhaseDiagramMatrix = M;
ScatteringLengths = SCATTERING_LENGTH_RANGE;
NumberOfAtoms = round(NUM_ATOMS_LIST * 1E-5,2);
@ -640,23 +715,23 @@ ylabel('Scattering Length (\times a_o)', 'Interpreter', 'tex', 'FontSize', 16);
grid on;
%% Edit phase diagram
load('./Results/PhaseDiagramUntilted.mat')
load('./Results/Data_Full3D/PhaseDiagramUntilted.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/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";
Scripts.plotSmoothedPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, titleString);
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";
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';
% SaveDirectory = './Results/Data_Full3D/GradientDescent/Phi020/aS_090_theta_020_phi_000_N_100000';
% JobNumber = 0;
SaveDirectory = './Results/Data_3D/PhaseDiagram';
SaveDirectory = './Results/Data_Full3D/PhaseDiagram';
JobNumber = 1;
% Load data