From d34178246232f3d0254acd7820e8e49c95e73c6c Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Wed, 28 May 2025 18:11:16 +0200 Subject: [PATCH] Latest additions - scripts to calculate average peak column density and unit cell density. --- .../extractAveragePeakColumnDensity.m | 83 ++++++++++++ .../+Scripts/extractAverageUnitCellDensity.m | 116 ++++++++++++++++ Dipolar-Gas-Simulator/+Scripts/run_locally.m | 128 +++++++++++++++--- 3 files changed, 311 insertions(+), 16 deletions(-) create mode 100644 Dipolar-Gas-Simulator/+Scripts/extractAveragePeakColumnDensity.m create mode 100644 Dipolar-Gas-Simulator/+Scripts/extractAverageUnitCellDensity.m diff --git a/Dipolar-Gas-Simulator/+Scripts/extractAveragePeakColumnDensity.m b/Dipolar-Gas-Simulator/+Scripts/extractAveragePeakColumnDensity.m new file mode 100644 index 0000000..a840d90 --- /dev/null +++ b/Dipolar-Gas-Simulator/+Scripts/extractAveragePeakColumnDensity.m @@ -0,0 +1,83 @@ +function [AveragePCD] = extractAveragePeakColumnDensity(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 = extractPeaks(state, radius, minPeakHeight); + peakHeights = state(peaks); + AveragePCD = mean(peakHeights); + [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([' = ' num2str(AveragePCD)], ... + 'Interpreter', 'tex', 'FontSize', 16) + set(gca, 'FontSize', 16); % For tick labels only + end +end + +function [peaks_mask] = extractPeaks(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; +end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Scripts/extractAverageUnitCellDensity.m b/Dipolar-Gas-Simulator/+Scripts/extractAverageUnitCellDensity.m new file mode 100644 index 0000000..05b56ea --- /dev/null +++ b/Dipolar-Gas-Simulator/+Scripts/extractAverageUnitCellDensity.m @@ -0,0 +1,116 @@ +function [UCD] = extractAverageUnitCellDensity(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 = extractPeaks(state, radius, minPeakHeight); + peakHeights = state(peaks); + [row, col] = find(peaks); + [~, maxIdx] = max(peakHeights); + MaxPeakLocation = [row(maxIdx), col(maxIdx)]; + + % Voronoi diagram of peak positions + points = [col, row]; % Voronoi uses [x, y] + [V, C] = voronoin(points); + + % Voronoi cell of the max peak + Vcell_indices = C{maxIdx}; + + % Plot the Voronoi cell + if all(Vcell_indices > 0) && all(Vcell_indices <= size(V, 1)) && ~any(isinf(V(Vcell_indices, 1))) + vx = interp1(1:numel(x), x, V(Vcell_indices,1), 'linear', 'extrap'); + vy = interp1(1:numel(y), y, V(Vcell_indices,2), 'linear', 'extrap'); + % Close the polygon by repeating the first vertex + vx(end+1) = vx(1); + vy(end+1) = vy(1); + % Compute area of the Voronoi cell polygon using the shoelace formula + AreaOfVoronoiCell = polyarea(vx, vy); % Area of Voronoi cell around max peak in um^2 + else + warning('Voronoi cell for max peak is unbounded or invalid.'); + end + + % Create grid points mesh + [X, Y] = meshgrid(x, y); % Note: size(X) and size(Y) match size(state) + + % Determine points inside Voronoi polygon + inCellMask = inpolygon(X, Y, vx, vy); + + % Sum all state values inside the Voronoi cell polygon + NumberOfParticlesInVoronoiCell = sum(state(inCellMask)); + + UCD = NumberOfParticlesInVoronoiCell / AreaOfVoronoiCell; + + 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(MaxPeakLocation(2)), y(MaxPeakLocation(1)), 'w+', 'MarkerSize', 8, 'LineWidth', 1.5); + plot(vx, vy, 'w-', '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(['UCD = ' num2str(UCD) ' \mum^{-2}'], ... + 'Interpreter', 'tex', 'FontSize', 16) + set(gca, 'FontSize', 16); % For tick labels only + end +end + +function [peaks_mask] = extractPeaks(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; +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 7e85654..4d602a2 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -483,18 +483,21 @@ JobNumber = 2; % 82 % JobNumber = 3; % 83 Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber) [contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction_in_plane_trap(SaveDirectory, JobNumber); + %% - Plot GS wavefunction % 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_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_Full3D/ApproximateLHY/AspectRatio2_8'; % SaveDirectory = './Results/Data_Full3D/ApproximateLHY/AspectRatio3_7'; @@ -502,36 +505,44 @@ SaveDirectory = './Results/Data_Full3D/ApproximateLHY/AspectRatio2_8'; % SaveDirectory = './Results/Data_Full3D/ApproximateLHY/AspectRatio3_9'; JobNumber = 3; Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) + %% % 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_Full3D/TiltedDipoles0'; JobNumber = 0; Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) + %% SaveDirectory = './Results/Data_Full3D/TiltedDipoles15'; JobNumber = 2; Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) + %% SaveDirectory = './Results/Data_Full3D/TiltedDipoles30'; JobNumber = 2; Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) + %% SaveDirectory = './Results/Data_Full3D/TiltedDipoles45'; JobNumber = 2; Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) + %% SaveDirectory = './Results/Data_Full3D/RealTimeDynamics'; JobNumber = 0; Plotter.makeMovie(SaveDirectory, JobNumber) + %% SaveDirectory = './Results/Data_Full3D/AnisotropicTrap/TiltedDipoles15'; JobNumber = 0; Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) + %% SaveDirectory = './Results/Data_Full3D/AnisotropicTrap/TiltedDipoles30'; JobNumber = 0; @@ -540,33 +551,56 @@ Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) SaveDirectory = './Results/Data_Full3D/AnisotropicTrap/TiltedDipoles35'; JobNumber = 0; Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) + %% SaveDirectory = './Results/Data_Full3D/AnisotropicTrap/TiltedDipoles40'; JobNumber = 0; Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) + %% SaveDirectory = './Results/Data_Full3D/AnisotropicTrap/TiltedDipoles45'; JobNumber = 0; Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) + %% SaveDirectory = './Results/Data_Full3D/PhaseDiagram/GradientDescent/'; JobNumber = 0; Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) + %% 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/Theta40/HighN/aS_9.250000e+01_theta_040_phi_000_N_916667'; +SaveDirectory = 'D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_712500'; 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'; +PeakThreshold = 3E3; +SaveDirectory = 'D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_712500'; JobNumber = 0; SuppressPlotFlag = false; DropletNumber = Scripts.identifyAndCountDroplets(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag); + +%% Extract average peak column density +Radius = 2; % The radius within which peaks will be considered duplicates +PeakThreshold = 3E3; +SaveDirectory = 'D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_712500'; +JobNumber = 0; +SuppressPlotFlag = false; +AveragePCD = Scripts.extractAveragePeakColumnDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag); + +%% Extract average unit cell density +Radius = 2; % The radius within which peaks will be considered duplicates +PeakThreshold = 3E3; +SaveDirectory = 'D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_712500'; +JobNumber = 0; +SuppressPlotFlag = false; +UCD = Scripts.extractAverageUnitCellDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag); + %% Plot number of droplets % Parameters Radius = 2; @@ -632,15 +666,31 @@ legend(legendEntries, 'Interpreter', 'tex', 'FontSize', 12, 'Location', 'bestout grid on; hold off; -%% Plot TF radii of unmodulated states +%% Plot average peak column density % Parameters +Radius = 2; +PeakThreshold = 400; JobNumber = 0; -TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 40^\circ"; +TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 0^\circ"; SuppressPlotFlag = true; -SCATTERING_LENGTH_RANGE = [97.71]; +SCATTERING_LENGTH_RANGE = [95.62]; -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]; +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]; + +% 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) @@ -650,7 +700,57 @@ for j = 1:length(SCATTERING_LENGTH_RANGE) aS_string = sprintf('%.6e', aS); % Construct base directory for this aS - baseDir = ['D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta40/aS_' ... + baseDir = ['D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_' ... + aS_string '_theta_000_phi_000_N_']; + + % Preallocate results for this curve + AverageCDs = 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 + AverageCDs(i) = Scripts.extractAveragePeakColumnDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag); + end + + % Plot curve + plot(NUM_ATOMS_LIST, AverageCDs, '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('Average Peak Column Density', 'Interpreter', 'tex', 'FontSize', 16); +title(TitleString, 'Interpreter', 'tex', 'FontSize', 18); +legend(legendEntries, 'Interpreter', 'tex', 'FontSize', 12, 'Location', 'bestoutside'); +grid on; +hold off; + +%% Plot TF radii of unmodulated states +% Parameters +JobNumber = 0; +TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 40^\circ"; +SuppressPlotFlag = true; + +SCATTERING_LENGTH_RANGE = [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]; + +% 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/Theta40/LowN/aS_' ... aS_string '_theta_040_phi_000_N_']; % Preallocate results for this curve: rows = N, cols = [Rx, Ry] @@ -711,6 +811,7 @@ axis square grid on; sgtitle(TitleString, 'Interpreter', 'tex', 'FontSize', 18); + %% Overlay curves % File list and associated theta labels @@ -799,7 +900,7 @@ 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 = [97.71]; +SCATTERING_LENGTH_RANGE = [95.62]; 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]; @@ -810,7 +911,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/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta40/HighN/aS_%.6e_theta_040_phi_000_N_%d', aS, N); + SaveDirectory = sprintf('D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_%.6e_theta_000_phi_000_N_%d', aS, N); fprintf('Processing JobNumber %d: %s\n', JobNumber, SaveDirectory); % Call the plotting function @@ -828,12 +929,7 @@ end %% Missing files data = { - 'aS_8.208000e+01_theta_000_phi_000_N_90909' - 'aS_8.729000e+01_theta_000_phi_000_N_77273' - 'aS_8.729000e+01_theta_000_phi_000_N_95455' - 'aS_8.729000e+01_theta_000_phi_000_N_3775000' - 'aS_9.250000e+01_theta_000_phi_000_N_5000000' - 'aS_9.771000e+01_theta_020_phi_000_N_3775000' + '' }; % Preallocate arrays