From 63904f2a95d2cf0fbf61b1c14ae9960ac173b76c Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Tue, 3 Jun 2025 16:48:11 +0200 Subject: [PATCH] MAJOR corrections to scripts to extract densities --- .../extractAveragePeakColumnDensity.m | 10 +- .../+Scripts/extractAverageUnitCellDensity.m | 241 ++++++++++++------ Dipolar-Gas-Simulator/+Scripts/run_locally.m | 153 ++++++++--- 3 files changed, 293 insertions(+), 111 deletions(-) diff --git a/Dipolar-Gas-Simulator/+Scripts/extractAveragePeakColumnDensity.m b/Dipolar-Gas-Simulator/+Scripts/extractAveragePeakColumnDensity.m index a840d90..a22c9a6 100644 --- a/Dipolar-Gas-Simulator/+Scripts/extractAveragePeakColumnDensity.m +++ b/Dipolar-Gas-Simulator/+Scripts/extractAveragePeakColumnDensity.m @@ -5,7 +5,15 @@ function [AveragePCD] = extractAveragePeakColumnDensity(folder_path, run_index, format long % Load data - Data = load(sprintf(horzcat(folder_path, '/Run_%03i/psi_gs.mat'),run_index),'psi','Params','Transf','Observ'); + filePath = sprintf(horzcat(folder_path, '/Run_%03i/psi_gs.mat'), run_index); + + try + Data = load(filePath, 'psi', 'Params', 'Transf', 'Observ'); + catch ME + warning('Failed to load file: %s\n%s', filePath, ME.message); + AveragePCD = NaN; + return; + end Params = Data.Params; Transf = Data.Transf; diff --git a/Dipolar-Gas-Simulator/+Scripts/extractAverageUnitCellDensity.m b/Dipolar-Gas-Simulator/+Scripts/extractAverageUnitCellDensity.m index 05b56ea..98a79a5 100644 --- a/Dipolar-Gas-Simulator/+Scripts/extractAverageUnitCellDensity.m +++ b/Dipolar-Gas-Simulator/+Scripts/extractAverageUnitCellDensity.m @@ -1,96 +1,177 @@ -function [UCD] = extractAverageUnitCellDensity(folder_path, run_index, radius, minPeakHeight, SuppressPlotFlag) +function [UCD] = extractAverageUnitCellDensity(folder_path, run_index, radius, minPeakHeight, SuppressPlotFlag) + set(0,'defaulttextInterpreter','latex') - set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex'); - - format long - + set(groot, 'defaultAxesTickLabelInterpreter','latex'); + set(groot, 'defaultLegendInterpreter','latex'); + % Load data - Data = load(sprintf(horzcat(folder_path, '/Run_%03i/psi_gs.mat'),run_index),'psi','Params','Transf','Observ'); - + filePath = sprintf(horzcat(folder_path, '/Run_%03i/psi_gs.mat'), run_index); + + try + Data = load(filePath, 'psi', 'Params', 'Transf', 'Observ'); + catch ME + warning('Failed to load file: %s\n%s', filePath, ME.message); + UCD = NaN; + return; + end + 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); + psi = gather(Data.psi); - % 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.'); + % Axes 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 integrated density + n = abs(psi).^2; + nxy = squeeze(trapz(n * dz, 3)); + state = nxy'; + + % Fourier spectrum for orientation detection + F = fftshift(abs(fft2(state - mean(state(:))))); + [nx, ny] = size(state); + kx = linspace(-pi / (x(2) - x(1)), pi / (x(2) - x(1)), nx); + ky = linspace(-pi / (y(2) - y(1)), pi / (y(2) - y(1)), ny); + [KX, KY] = meshgrid(kx, ky); + theta = atan2(KY, KX); + + % Remove center/DC + R = sqrt(KX.^2 + KY.^2); + F(R < 0.1 * max(kx(:))) = 0; + + % Angular histogram of power spectrum + nbins = 180; + angles = linspace(-pi, pi, nbins); + powerAngular = zeros(1, nbins); + for k = 1:nbins-1 + mask = theta >= angles(k) & theta < angles(k+1); + powerAngular(k) = sum(F(mask), 'all'); 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)); + powerAngular(end) = powerAngular(1); % wrap around - UCD = NumberOfParticlesInVoronoiCell / AreaOfVoronoiCell; + % Anisotropy measure + peakRatio = max(powerAngular) / mean(powerAngular); + % Threshold to distinguish stripe vs droplet + isStripe = peakRatio > 20; + + if ~isStripe + % === DROPLET MODE: Voronoi cell === + peaks = extractPeaks(state, radius, minPeakHeight); + [row, col] = find(peaks); + peakHeights = state(peaks); + [~, maxIdx] = max(peakHeights); + MaxPeakLocation = [row(maxIdx), col(maxIdx)]; + points = [col, row]; + + [V, C] = voronoin(points); + Vcell_indices = C{maxIdx}; + + 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'); + vx(end+1) = vx(1); + vy(end+1) = vy(1); + AreaOfVoronoiCell = polyarea(vx, vy); + + [X, Y] = meshgrid(x, y); + inCellMask = inpolygon(X, Y, vx, vy); + NumberOfParticles = sum(state(inCellMask)); + UCD = NumberOfParticles / AreaOfVoronoiCell; + else + warning('Voronoi cell for max peak is invalid.'); + UCD = NaN; + end + + else + % === STRIPE MODE: Use rectangular unit cell aligned with stripe direction === + [~, idxMax] = max(F(:)); + stripeAngle = theta(idxMax); % radians + + % Rotate image to align stripes horizontally + rotatedState = imrotate(state, -rad2deg(stripeAngle), 'bilinear', 'crop'); + + % Estimate vertical stripe spacing (stripe-normal direction) + stripeProfileY = mean(rotatedState, 2) - mean(rotatedState(:)); + fftY = abs(fft(stripeProfileY)); + fftY(1:3) = 0; + [~, fyIdx] = max(fftY(1:floor(end/2))); + spacingY = round(numel(stripeProfileY) / fyIdx); + + % Estimate horizontal spacing (along-stripe periodicity) + stripeProfileX = mean(rotatedState, 1) - mean(rotatedState(:)); + fftX = abs(fft(stripeProfileX)); + fftX(1:3) = 0; + [~, fxIdx] = max(fftX(1:floor(end/2))); + spacingX = round(numel(stripeProfileX) / fxIdx); + + % Find stripe center (vertical) + [~, centerY] = max(stripeProfileY); + rowRange = max(1, centerY - round(0.5 * spacingY)) : ... + min(size(rotatedState,1), centerY + round(0.5 * spacingY)); + + % Find repeat center along stripe direction (horizontal) + [~, centerX] = max(stripeProfileX); + colRange = max(1, centerX - round(0.5 * spacingX)) : ... + min(size(rotatedState,2), centerX + round(0.5 * spacingX)); + + % Extract unit cell region + unitCellRegion = rotatedState(rowRange, colRange); + NumberOfParticles = sum(unitCellRegion(:)); + AreaOfUnitCell = numel(unitCellRegion) * abs(x(2)-x(1)) * abs(y(2)-y(1)); + UCD = NumberOfParticles / AreaOfUnitCell; + end + + + % Optional plot 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 + if ~isStripe + plotxy = pcolor(x,y,state); + set(plotxy, 'EdgeColor', 'none'); + hold on; + plot(vx, vy, 'w-', 'LineWidth', 2); + 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 + else + imagesc(rotatedState); + axis image; + cbar1 = colorbar; + cbar1.Label.Interpreter = 'latex'; + colormap(gca, Helper.Colormaps.plasma()) + xlabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16) + ylabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16) + + % Title with UCD + title(['UCD = ' num2str(UCD, '%.3f') ' $\mu$m$^{-2}$'], ... + 'Interpreter', 'latex', 'FontSize', 16) + + % Highlight unit cell region + rectX = colRange(1) - 0.5; + rectY = rowRange(1) - 0.5; + rectW = length(colRange); + rectH = length(rowRange); + rectangle('Position', [rectX, rectY, rectW, rectH], ... + 'EdgeColor', 'w', 'LineWidth', 2); + title(sprintf('UCD = %.4f $\\mu$m$^{-2}$', UCD), ... + 'Interpreter', 'latex', 'FontSize', 16); + + + end end -end +end function [peaks_mask] = extractPeaks(image, radius, minPeakHeight) peaks = imregionalmax(image); diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m index 8362af4..75ff163 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -593,7 +593,7 @@ JobNumber = 0; SuppressPlotFlag = false; AveragePCD = Scripts.extractAveragePeakColumnDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag); -%% Extract average unit cell density +%% Extract average unit cell density - Droplets 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'; @@ -601,6 +601,19 @@ JobNumber = 0; SuppressPlotFlag = false; UCD = Scripts.extractAverageUnitCellDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag); +%% Extract average unit cell density - Stripes +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_1529167'; +JobNumber = 0; +SuppressPlotFlag = false; +UCD = Scripts.extractAverageUnitCellDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag); + +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]; + %% Plot number of droplets % Parameters Radius = 2; @@ -676,60 +689,140 @@ SuppressPlotFlag = true; 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_FULL = [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_INSET = [50000 54545 59091 63636 68182 72727 77273 81818 86364 90909 95455]; % Prepare figure figure(1); clf; set(gcf,'Position', [100, 100, 1000, 700]); -hold on; -% Color order for better visibility +% Color order colors = lines(length(SCATTERING_LENGTH_RANGE)); -% Store legend labels -legendEntries = cell(1, length(SCATTERING_LENGTH_RANGE)); +% Store data for both sets +AverageCDs_full = zeros(length(SCATTERING_LENGTH_RANGE), length(NUM_ATOMS_LIST_FULL)); +AverageCDs_inset = zeros(length(SCATTERING_LENGTH_RANGE), length(NUM_ATOMS_LIST_INSET)); + +% Main plot +main_ax = axes; +hold(main_ax, 'on'); -% 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/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); + % Full list + for i = 1:length(NUM_ATOMS_LIST_FULL) + N = NUM_ATOMS_LIST_FULL(i); SaveDirectory = [baseDir num2str(N)]; - - % Call your droplet counting function - AverageCDs(i) = Scripts.extractAveragePeakColumnDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag); + AverageCDs_full(j,i) = Scripts.extractAveragePeakColumnDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag); end - % Plot curve - plot(NUM_ATOMS_LIST, AverageCDs, 'o-', ... - 'Color', colors(j, :), 'LineWidth', 1.5); + baseDir = ['D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/LowN/aS_' ... + aS_string '_theta_000_phi_000_N_']; + + % Inset list + for i = 1:length(NUM_ATOMS_LIST_INSET) + N = NUM_ATOMS_LIST_INSET(i); + SaveDirectory = [baseDir num2str(N)]; + AverageCDs_inset(j,i) = Scripts.extractAveragePeakColumnDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag); + end - % Store legend entry - legendEntries{j} = ['a_s = ' num2str(aS) 'a_o']; + % Plot main curve + x = NUM_ATOMS_LIST_FULL; + y = AverageCDs_full(j,:); + valid = ~isnan(y); % logical index of valid points + plot(main_ax,x(valid), y(valid), 'o-', ... + 'Color', colors(j,:), 'LineWidth', 1.5); 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'); +legend(arrayfun(@(aS) sprintf('a_s = %.2f a_0', aS), SCATTERING_LENGTH_RANGE, ... + 'UniformOutput', false), ... + 'Interpreter', 'tex', 'FontSize', 12, 'Location', 'bestoutside'); +grid(main_ax, 'on'); + +% Inset plot +inset_ax = axes('Position', [0.45 0.18 0.28 0.28]); % Normalized position [x y w h] +box(inset_ax, 'on'); +hold(inset_ax, 'on'); + +for j = 1:length(SCATTERING_LENGTH_RANGE) + plot(inset_ax, NUM_ATOMS_LIST_INSET, AverageCDs_inset(j,:), 'o-', ... + 'Color', colors(j,:), 'LineWidth', 1.2); +end + +set(inset_ax, 'FontSize', 8); +title(inset_ax, 'Low-N', 'FontSize', 9); +grid(inset_ax, 'on'); +xlabel(inset_ax, 'N', 'FontSize', 9); +ylabel(inset_ax, 'CD', 'FontSize', 9); + +%% Plot average unit cell density +Radius = 2; +PeakThreshold = 3E3; +JobNumber = 0; +SuppressPlotFlag = true; % Suppress plots during batch processing +TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 0^\circ"; + + +SCATTERING_LENGTH_RANGE = [95.62]; + +NUM_ATOMS_LIST = [712500 916667 1120833 1325000 1529167 1733333 1937500 2141667 2345833 2550000 2754167 2958333 3162500 3366667 3570833]; + +UCD_values = zeros(length(SCATTERING_LENGTH_RANGE), length(NUM_ATOMS_LIST)); + +% Prepare figure +figure(1); +clf; +set(gcf,'Position', [100, 100, 1000, 700]); +hold on + +% Color order +colors = lines(length(SCATTERING_LENGTH_RANGE)); + +for j = 1:length(SCATTERING_LENGTH_RANGE) + aS = SCATTERING_LENGTH_RANGE(j); + aS_string = sprintf('%.6e', aS); + + baseDir = ['D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_' ... + aS_string '_theta_000_phi_000_N_']; + + for i = 1:length(NUM_ATOMS_LIST) + N = NUM_ATOMS_LIST(i); + % Construct folder path for this N + SaveDirectory = sprintf('%s%d', baseDir, N); + try + UCD_values(j,i) = Scripts.extractAverageUnitCellDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag); + catch ME + warning('Error processing N=%d: %s', N, ME.message); + UCD_values(j,i) = NaN; % mark as NaN on error + end + end + + x = NUM_ATOMS_LIST; + y = UCD_values(j,:); + + valid = ~isnan(y); % logical index of valid points + plot(x(valid), y(valid), 'o-', 'Color', colors(j,:), 'LineWidth', 1.5); + +end + +xlabel('Number of Atoms', 'Interpreter', 'latex', 'FontSize', 16); +ylabel('Unit Cell Density (UCD) [$\mu m^{-2}$]', 'Interpreter', 'latex', 'FontSize', 16); +title(TitleString, 'Interpreter', 'tex', 'FontSize', 18); +legend(arrayfun(@(aS) sprintf('a_s = %.2f a_0', aS), SCATTERING_LENGTH_RANGE, ... + 'UniformOutput', false), ... + 'Interpreter', 'tex', 'FontSize', 12, 'Location', 'bestoutside'); +set(gca, 'FontSize', 14); grid on; -hold off; %% Plot TF radii of unmodulated states % Parameters