From 8ae97eb26775ac6313176c85b1508762a7934f7b Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Fri, 16 May 2025 18:44:04 +0200 Subject: [PATCH] Added script to identify and count droplets in the results of the numerics. --- .../+Scripts/identifyAndCountDroplets.m | 84 ++++++++++ Dipolar-Gas-Simulator/+Scripts/run_locally.m | 143 +++++++++++++----- 2 files changed, 193 insertions(+), 34 deletions(-) create mode 100644 Dipolar-Gas-Simulator/+Scripts/identifyAndCountDroplets.m diff --git a/Dipolar-Gas-Simulator/+Scripts/identifyAndCountDroplets.m b/Dipolar-Gas-Simulator/+Scripts/identifyAndCountDroplets.m new file mode 100644 index 0000000..d07febe --- /dev/null +++ b/Dipolar-Gas-Simulator/+Scripts/identifyAndCountDroplets.m @@ -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 \ 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 29414a9..d8660c4 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -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