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

MAJOR corrections to scripts to extract densities

Browse Source
This commit is contained in:
Karthik 2025-06-03 16:48:11 +02:00
parent 4416bc58b6
commit 63904f2a95
3 changed files with 293 additions and 111 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
Dipolar-Gas-Simulator/+Scripts/extractAveragePeakColumnDensity.m
View File

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

241
Dipolar-Gas-Simulator/+Scripts/extractAverageUnitCellDensity.m
View File

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

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

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