From e256caa0dcb1e3e69b7180052016fcb8229128fa Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Mon, 2 Jun 2025 17:12:37 +0200 Subject: [PATCH 1/4] Adjusted parameters, new script to run Fourier analysis on simulated data --- .../+Scripts/conductFourierAnalysis.m | 296 ++++++++++++++++++ .../+Scripts/run_hybrid_worker.m | 4 +- Dipolar-Gas-Simulator/+Scripts/run_locally.m | 66 ++++ .../+Scripts/run_on_cluster.m | 4 +- 4 files changed, 366 insertions(+), 4 deletions(-) create mode 100644 Dipolar-Gas-Simulator/+Scripts/conductFourierAnalysis.m diff --git a/Dipolar-Gas-Simulator/+Scripts/conductFourierAnalysis.m b/Dipolar-Gas-Simulator/+Scripts/conductFourierAnalysis.m new file mode 100644 index 0000000..0a49e39 --- /dev/null +++ b/Dipolar-Gas-Simulator/+Scripts/conductFourierAnalysis.m @@ -0,0 +1,296 @@ +function [spectral_weight, g2, theta_vals] = conductFourierAnalysis(folder_path, run_index, N_bins, Threshold, Sigma, SuppressPlotFlag) + + arguments + folder_path (1,:) char + run_index (1,:) {mustBeNumeric,mustBeReal} + N_bins (1,:) {mustBeNumeric,mustBeReal} + Threshold (1,:) {mustBeNumeric,mustBeReal} + Sigma (1,:) {mustBeNumeric,mustBeReal} + SuppressPlotFlag (1,:) logical = true + end + + set(0,'defaulttextInterpreter','latex') + set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex'); + + % Load data + Data = load(fullfile(fullfile(folder_path, sprintf('Run_%03i', run_index)), 'psi_gs.mat'), 'psi', 'Transf', 'Observ', 'Params'); + + 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 + + alpha = Params.theta; + + % 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); + + % Calculate frequency increment (frequency axes) + Nx = length(x); % grid size along X + Ny = length(y); % grid size along Y + dx = mean(diff(x)); % real space increment in the X direction (in micrometers) + dy = mean(diff(y)); % real space increment in the Y direction (in micrometers) + dvx = 1 / (Nx * dx); % reciprocal space increment in the X direction (in micrometers^-1) + dvy = 1 / (Ny * dy); % reciprocal space increment in the Y direction (in micrometers^-1) + + % Create the frequency axes + vx = (-Nx/2:Nx/2-1) * dvx; % Frequency axis in X (micrometers^-1) + vy = (-Ny/2:Ny/2-1) * dvy; % Frequency axis in Y (micrometers^-1) + + % Calculate maximum frequencies + % kx_max = pi / dx; + % ky_max = pi / dy; + + % Generate reciprocal axes + % kx = linspace(-kx_max, kx_max * (Nx-2)/Nx, Nx); + % ky = linspace(-ky_max, ky_max * (Ny-2)/Ny, Ny); + + % Create the Wavenumber axes + kx = 2*pi*vx; % Wavenumber axis in X + ky = 2*pi*vy; % Wavenumber axis in Y + + % Compute probability density |psi|^2 + n = abs(psi).^2; + + nxy = squeeze(trapz(n*dz,3)); + + skipPreprocessing = true; + skipMasking = true; + skipIntensityThresholding = true; + skipBinarization = true; + + %% Extract Spectral Weight and g2 + + IMG = nxy; + + [IMGFFT, ~] = computeFourierTransform(IMG, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization); + + [theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(IMGFFT, 10, 35, N_bins, Threshold, Sigma); + + spectral_weight = trapz(theta_vals, S_theta); + + g2 = zeros(1, N_bins); % Preallocate + + for dtheta = 0:N_bins-1 + profile = S_theta; + profile_shifted = circshift(profile, -dtheta, 2); + + num = mean(profile .* profile_shifted); + denom = mean(profile)^2; + + g2(dtheta+1) = num / denom - 1; + end + + if ~SuppressPlotFlag + figure(1); + clf + set(gcf,'Position',[500 100 1000 800]) + t = tiledlayout(2, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 1x4 grid + font = 'Bahnschrift'; + % Display the cropped OD image + ax1 = nexttile; + plotxy = pcolor(x,y,IMG'); + set(plotxy, 'EdgeColor', 'none'); + % Define normalized positions (relative to axis limits) + x_offset = 0.025; % 5% offset from the edges + y_offset = 0.025; % 5% offset from the edges + % Top-right corner (normalized axis coordinates) + hText = text(1 - x_offset, 1 - y_offset, ... + ['\alpha = ', num2str(rad2deg(alpha), '%.1f'), '^\circ'], ... + 'Color', 'white', 'FontWeight', 'bold', ... + 'Interpreter', 'tex', 'FontSize', 20, ... + 'Units', 'normalized', ... + 'HorizontalAlignment', 'right', ... + 'VerticalAlignment', 'top'); + axis square; + hcb = colorbar; + hcb.Label.Interpreter = 'latex'; + colormap(gca, Helper.Colormaps.plasma()) + set(gca, 'FontSize', 14); % For tick labels only + hL = ylabel(hcb, 'Column Density'); + hXLabel = xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14); + hYLabel = ylabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14); + hTitle = title('$|\Psi(x,y)|^2$', 'Interpreter', 'latex', 'FontSize', 14) ; + set(hText, 'FontName', font) + set([hXLabel, hYLabel, hL], 'FontSize', 14) + set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title + + + % Plot the power spectrum + nexttile; + imagesc(kx, ky, log(1 + abs(IMGFFT).^2)); + axis square; + hcb = colorbar; + colormap(gca, Helper.Colormaps.plasma()) + set(gca, 'FontSize', 14); % For tick labels only + set(gca,'YDir','normal') + hXLabel = xlabel('k_x', 'Interpreter', 'tex'); + hYLabel = ylabel('k_y', 'Interpreter', 'tex'); + hTitle = title('Power Spectrum - S(k_x,k_y)', 'Interpreter', 'tex'); + set([hXLabel, hYLabel], 'FontName', font) + set([hXLabel, hYLabel], 'FontSize', 14) + set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title + + % Plot the angular distribution + nexttile + plot(theta_vals/pi, S_theta,'Linewidth',2); + set(gca, 'FontSize', 14); % For tick labels only + hXLabel = xlabel('\theta/\pi [rad]', 'Interpreter', 'tex'); + hYLabel = ylabel('Normalized magnitude (a.u.)', 'Interpreter', 'tex'); + hTitle = title('Angular Spectral Distribution - S(\theta)', 'Interpreter', 'tex'); + set([hXLabel, hYLabel], 'FontName', font) + set([hXLabel, hYLabel], 'FontSize', 14) + set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title + grid on + + nexttile + plot(theta_vals/pi, g2, 'o-', 'LineWidth', 1.2, 'MarkerSize', 5); + set(gca, 'FontSize', 14); + ylim([-1.5 3.0]); % Set y-axis limits here + hXLabel = xlabel('$\delta\theta / \pi$', 'Interpreter', 'latex'); + hYLabel = ylabel('$g^{(2)}(\delta\theta)$', 'Interpreter', 'latex'); + hTitle = title('Autocorrelation', 'Interpreter', 'tex'); + set([hXLabel, hYLabel], 'FontName', font) + set([hXLabel, hYLabel], 'FontSize', 14) + set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title + grid on; + end +end + +%% Helper Functions +function [IMGFFT, IMGPR] = computeFourierTransform(I, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization) + % computeFourierSpectrum - Computes the 2D Fourier power spectrum + % of binarized and enhanced lattice image features, with optional central mask. + % + % Inputs: + % I - Grayscale or RGB image matrix + % + % Output: + % F_mag - 2D Fourier power spectrum (shifted) + + if ~skipPreprocessing + % Preprocessing: Denoise + filtered = imgaussfilt(I, 10); + IMGPR = I - filtered; % adjust sigma as needed + else + IMGPR = I; + end + + if ~skipMasking + [rows, cols] = size(IMGPR); + [X, Y] = meshgrid(1:cols, 1:rows); + % Elliptical mask parameters + cx = cols / 2; + cy = rows / 2; + + % Shifted coordinates + x = X - cx; + y = Y - cy; + + % Ellipse semi-axes + rx = 0.4 * cols; + ry = 0.2 * rows; + + % Rotation angle in degrees -> radians + theta_deg = 30; % Adjust as needed + theta = deg2rad(theta_deg); + + % Rotated ellipse equation + cos_t = cos(theta); + sin_t = sin(theta); + + x_rot = (x * cos_t + y * sin_t); + y_rot = (-x * sin_t + y * cos_t); + + ellipseMask = (x_rot.^2) / rx^2 + (y_rot.^2) / ry^2 <= 1; + + % Apply cutout mask + IMGPR = IMGPR .* ellipseMask; + end + + if ~skipIntensityThresholding + % Apply global intensity threshold mask + intensity_thresh = 0.20; + intensity_mask = IMGPR > intensity_thresh; + IMGPR = IMGPR .* intensity_mask; + end + + if ~skipBinarization + % Adaptive binarization and cleanup + IMGPR = imbinarize(IMGPR, 'adaptive', 'Sensitivity', 0.0); + IMGPR = imdilate(IMGPR, strel('disk', 2)); + IMGPR = imerode(IMGPR, strel('disk', 1)); + IMGPR = imfill(IMGPR, 'holes'); + F = fft2(double(IMGPR)); % Compute 2D Fourier Transform + IMGFFT = abs(fftshift(F))'; % Shift zero frequency to center + else + F = fft2(double(IMGPR)); % Compute 2D Fourier Transform + IMGFFT = abs(fftshift(F))'; % Shift zero frequency to center + end +end + +function [theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(IMGFFT, r_min, r_max, num_bins, threshold, sigma) + % Apply threshold to isolate strong peaks + IMGFFT(IMGFFT < threshold) = 0; + + % Prepare polar coordinates + [ny, nx] = size(IMGFFT); + [X, Y] = meshgrid(1:nx, 1:ny); + cx = ceil(nx/2); + cy = ceil(ny/2); + R = sqrt((X - cx).^2 + (Y - cy).^2); + Theta = atan2(Y - cy, X - cx); % range [-pi, pi] + + % Choose radial band + radial_mask = (R >= r_min) & (R <= r_max); + + % Initialize the angular structure factor array + S_theta = zeros(1, num_bins); % Pre-allocate for 180 angle bins + % Define the angle values for the x-axis + theta_vals = linspace(0, pi, num_bins); + + % Loop through each angle bin + for i = 1:num_bins + angle_start = (i-1) * pi / num_bins; + angle_end = i * pi / num_bins; + + % Define a mask for the given angle range + angle_mask = (Theta >= angle_start & Theta < angle_end); + + bin_mask = radial_mask & angle_mask; + + % Extract the Fourier components for the given angle + fft_angle = IMGFFT .* bin_mask; + + % Integrate the Fourier components over the radius at the angle + S_theta(i) = sum(sum(abs(fft_angle).^2)); % sum of squared magnitudes + end + + % Create a 1D Gaussian kernel + half_width = ceil(3 * sigma); + x = -half_width:half_width; + gauss_kernel = exp(-x.^2 / (2 * sigma^2)); + gauss_kernel = gauss_kernel / sum(gauss_kernel); % normalize + + % Apply convolution (circular padding to preserve periodicity) + S_theta = conv([S_theta(end-half_width+1:end), S_theta, S_theta(1:half_width)], gauss_kernel, 'same'); + S_theta = S_theta(half_width+1:end-half_width); % crop back to original size + + % Normalize to 1 + S_theta = S_theta / max(S_theta); +end + diff --git a/Dipolar-Gas-Simulator/+Scripts/run_hybrid_worker.m b/Dipolar-Gas-Simulator/+Scripts/run_hybrid_worker.m index 8cf74e6..11fa382 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_hybrid_worker.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_hybrid_worker.m @@ -48,7 +48,7 @@ function run_hybrid_worker(batchParams, batchIdx) OptionsStruct.MaxIterationsForGD = 15000; OptionsStruct.TimeStepSize = 1E-3; % in s OptionsStruct.MinimumTimeStepSize = 1E-6; % in s - OptionsStruct.TimeCutOff = 2E5; % in s + OptionsStruct.TimeCutOff = 1E6; % in s OptionsStruct.EnergyTolerance = 5E-10; OptionsStruct.ResidualTolerance = 1E-05; OptionsStruct.NoiseScaleFactor = 0.010; @@ -77,4 +77,4 @@ function run_hybrid_worker(batchParams, batchIdx) end delete(pool); -end +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 bcb7984..8362af4 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -1066,3 +1066,69 @@ else disp('The state is not modulated'); end +%% Plot Spectral weight and g2 across transition for simulated data +N_atoms = 5E5; +a_s = 95.62; +phi_deg = 0.0; +alpha_vals = [0.0 5.0 10.0 15.0 20.0 25.0 30.0 35.0 40.0 45.0]; + +N_alpha = length(alpha_vals); +N_bins = 90; +Threshold = 75; +Sigma = 2; + +spectral_weights = zeros(size(alpha_vals)); +g2_all = zeros(N_alpha, N_bins); +theta_vals_all = zeros(N_alpha, N_bins); + +JobNumber = 1; +SuppressPlotFlag = true; + +for i = 1:N_alpha + folderName = sprintf('aS_%03d_theta_%03d_phi_%03d_N_%d', a_s, alpha_vals(i), phi_deg, N_atoms); + SaveDirectory = fullfile('D:/Results - Numerics/Data_Full3D/PhaseTransition/DTS/', folderName); + [spectral_weight, g2, theta_vals] = Scripts.conductFourierAnalysis(SaveDirectory, JobNumber, N_bins, Threshold, Sigma, SuppressPlotFlag); + spectral_weights(i) = spectral_weight; % Store the spectral weight for the current alpha value + g2_all(i, :) = g2; % Store the g2 results for the current alpha value + theta_vals_all(i, :) = theta_vals; +end + +figure(2); +set(gcf,'Position',[100 100 950 750]) +font = 'Bahnschrift'; +plot(alpha_vals, spectral_weights, 'o--', ... + 'LineWidth', 1.5, 'MarkerSize', 6); +set(gca, 'FontSize', 14); % For tick labels only +hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex'); +hYLabel = ylabel('Spectral Weight', 'Interpreter', 'tex'); +hTitle = title('Change across transition', 'Interpreter', 'tex'); +set([hXLabel, hYLabel], 'FontName', font) +set([hXLabel, hYLabel], 'FontSize', 14) +set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title +grid on + +figure(3); +set(gcf,'Position',[100 100 950 750]) +font = 'Bahnschrift'; +legend_entries = cell(N_alpha, 1); +theta_vals = theta_vals_all(1, :); +cmap = sky(N_alpha); % Generate a colormap with enough unique colors +hold on + +for i = 1:N_alpha + plot(theta_vals/pi, g2_all(i, :), ... + 'o-', 'Color', cmap(i,:), 'LineWidth', 1.2, ... + 'MarkerSize', 5); + legend_entries{i} = sprintf('$\\alpha = %g^\\circ$', alpha_vals(i)); +end +set(gca, 'FontSize', 14); +ylim([-1.5 10.0]); % Set y-axis limits here +hXLabel = xlabel('$\delta\theta / \pi$', 'Interpreter', 'latex'); +hYLabel = ylabel('$g^{(2)}(\delta\theta)$', 'Interpreter', 'latex'); +hTitle = title('Transition from Droplets to Stripes', 'Interpreter', 'tex'); +legend(legend_entries, 'Interpreter', 'latex', 'Location', 'bestoutside'); +set([hXLabel, hYLabel], 'FontName', font) +set([hXLabel, hYLabel], 'FontSize', 14) +set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title +grid on; + diff --git a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m index 5a16eab..62efda6 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m @@ -46,13 +46,13 @@ function run_on_cluster(batchParams, batchIdx) OptionsStruct.MaxIterationsForGD = 15000; OptionsStruct.TimeStepSize = 1E-3; OptionsStruct.MinimumTimeStepSize = 1E-6; - OptionsStruct.TimeCutOff = 5E5; + OptionsStruct.TimeCutOff = 2E5; OptionsStruct.EnergyTolerance = 5E-08; OptionsStruct.ResidualTolerance = 1E-05; OptionsStruct.NoiseScaleFactor = 0.010; OptionsStruct.PlotLive = false; - OptionsStruct.JobNumber = k; + OptionsStruct.JobNumber = 0; OptionsStruct.RunOnGPU = true; OptionsStruct.SaveData = true; OptionsStruct.SaveDirectory = saveDir; From 4416bc58b649de113c917fcdfaffd7ef2ac61463 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Tue, 3 Jun 2025 12:38:11 +0200 Subject: [PATCH 2/4] New analysis script - characterization of lattice geometries in real space --- Data-Analyzer/computeBondOrderParameters.m | 179 +++++++++++++++++++++ Data-Analyzer/execution_scripts.m | 63 ++++++++ Data-Analyzer/extractHexaticOrder.m | 98 +++++++++++ 3 files changed, 340 insertions(+) create mode 100644 Data-Analyzer/computeBondOrderParameters.m create mode 100644 Data-Analyzer/execution_scripts.m create mode 100644 Data-Analyzer/extractHexaticOrder.m diff --git a/Data-Analyzer/computeBondOrderParameters.m b/Data-Analyzer/computeBondOrderParameters.m new file mode 100644 index 0000000..fe3ecee --- /dev/null +++ b/Data-Analyzer/computeBondOrderParameters.m @@ -0,0 +1,179 @@ +function [psi_global, order_ratio, Npoints] = computeBondOrderParameters(I) +% Analyze bond-orientational order in a 2D image, restricted to a circular region +% +% Input: +% I - 2D grayscale image (double or uint8) +% +% Outputs: +% psi_global - struct with fields psi2, psi4, psi6 (⟨|ψₙ|⟩ values) +% order_ratio - scalar = ⟨|ψ₆|⟩ / ⟨|ψ₂|⟩ +% Npoints - number of detected points used + + %% Step 1: Create mask + [Ny, Nx] = size(I); + % Parameters for elliptical mask + a = Nx / 2.25; % Semi-major axis (x-direction) + b = Ny / 2.25; % Semi-minor axis (y-direction) + theta_deg = 0; % Rotation angle in degrees + theta = deg2rad(theta_deg); % Convert to radians + + % Meshgrid coordinates + [Ny, Nx] = size(I); + [X, Y] = meshgrid(1:Nx, 1:Ny); + cx = Nx / 2; + cy = Ny / 2; + + % Shifted coordinates + Xc = X - cx; + Yc = Y - cy; + + % Rotate coordinates + Xr = cos(theta) * Xc + sin(theta) * Yc; + Yr = -sin(theta) * Xc + cos(theta) * Yc; + + % Elliptical mask equation + mask = (Xr / a).^2 + (Yr / b).^2 <= 1; + + % Apply mask to image + I = double(I) .* mask; + + %% Step 2: Detect local maxima within mask + I_smooth = imgaussfilt(I, 2); + peaks = imregionalmax(I_smooth); + [y, x] = find(peaks); + Npoints = length(x); + + if Npoints < 6 + warning('Too few points (%d) to compute bond order meaningfully.', Npoints); + psi_global = struct('psi2', NaN, 'psi4', NaN, 'psi6', NaN); + order_ratio = NaN; + return; + end + + %% Step 3: Try Delaunay triangulation + use_kNN = false; + try + tri = delaunay(x, y); + catch + warning('Delaunay triangulation failed - falling back to kNN neighbor search'); + use_kNN = true; + end + + % 3) Find neighbors + neighbors = cell(Npoints, 1); + if ~use_kNN + % Use Delaunay neighbors + for i = 1:size(tri,1) + for j = 1:3 + idx = tri(i,j); + nbrs = tri(i,[1 2 3] ~= j); + neighbors{idx} = unique([neighbors{idx}, nbrs]); + end + end + else + % Use kNN neighbors (e.g. k=6) + k = min(6, Npoints-1); + pts = [x, y]; + [idx, ~] = knnsearch(pts, pts, 'K', k+1); % +1 for self + for j = 1:Npoints + neighbors{j} = idx(j, 2:end); + end + end + + % 4) Compute bond order parameters + psi2 = zeros(Npoints,1); + psi4 = zeros(Npoints,1); + psi6 = zeros(Npoints,1); + + for j = 1:Npoints + nbrs = neighbors{j}; + nbrs(nbrs == j) = []; + if isempty(nbrs) + continue; + end + angles = atan2(y(nbrs) - y(j), x(nbrs) - x(j)); + psi2(j) = abs(mean(exp(1i * 2 * angles))); + psi4(j) = abs(mean(exp(1i * 4 * angles))); + psi6(j) = abs(mean(exp(1i * 6 * angles))); + end + + psi_global.psi2 = mean(psi2); + psi_global.psi4 = mean(psi4); + psi_global.psi6 = mean(psi6); + + order_ratio = psi_global.psi6 / psi_global.psi2; + + % --- After neighbors are found and (x,y) extracted --- + + % 1) Angular distribution of all bonds pooled from all points + all_angles = []; + for j = 1:Npoints + nbrs = neighbors{j}; + nbrs(nbrs == j) = []; + if isempty(nbrs), continue; end + angles = atan2(y(nbrs)-y(j), x(nbrs)-x(j)); + all_angles = [all_angles; angles(:)]; + end + + % Histogram of bond angles over [-pi, pi] + nbins = 36; % 10 degree bins + edges = linspace(-pi, pi, nbins+1); + counts = histcounts(all_angles, edges, 'Normalization','probability'); + + % Circular variance of bond angles + R = abs(mean(exp(1i*all_angles))); + circ_variance = 1 - R; % 0 means angles clustered; 1 means uniform + + % 2) Positional correlation function g(r) estimate + + % Compute pairwise distances + coords = [x y]; + D = pdist(coords); + max_r = min([Nx, Ny])/2; % max radius for g(r) + + % Compute histogram of distances + dr = 1; % bin width in pixels (adjust) + r_edges = 0:dr:max_r; + [counts_r, r_bins] = histcounts(D, r_edges); + + % Normalize g(r) by ideal gas distribution (circle annulus area) + r_centers = r_edges(1:end-1) + dr/2; + area_annulus = pi*((r_edges(2:end)).^2 - (r_edges(1:end-1)).^2); + density = Npoints / (Nx*Ny); + + g_r = counts_r ./ (density * area_annulus * Npoints); + + figure(1); + clf + set(gcf,'Position',[500 100 1250 500]) + t = tiledlayout(1, 3, 'TileSpacing', 'compact', 'Padding', 'compact'); % 1x4 grid + + nexttile + imshow(I, []); + hold on; + scatter(x, y, 30, 'w', 'filled'); + hold off; + colormap(Helper.Colormaps.plasma()); % Example colormap for original image plot + cbar1 = colorbar; + cbar1.Label.Interpreter = 'latex'; + xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14); + ylabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14); + title('$|\Psi(x,y)|^2$', 'Interpreter', 'latex', 'FontSize', 14); + axis square; + + % --- Plot angular distribution --- + nexttile + polarhistogram(all_angles, nbins, 'Normalization', 'probability'); + title(sprintf('Bond angle distribution (circ. variance = %.3f)', circ_variance)); + + % --- Plot g(r) --- + nexttile + plot(r_centers, g_r, 'LineWidth', 2); + xlabel('r (pixels)'); + ylabel('g(r)'); + title('Radial Pair Correlation Function g(r)'); + grid on; + axis square + + +end \ No newline at end of file diff --git a/Data-Analyzer/execution_scripts.m b/Data-Analyzer/execution_scripts.m new file mode 100644 index 0000000..653cfda --- /dev/null +++ b/Data-Analyzer/execution_scripts.m @@ -0,0 +1,63 @@ +%% + +Data = load('D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_712500/Run_000/psi_gs.mat','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)); +IMG = nxy; + +extractHexaticOrder(IMG) + +%% + +Data = load('D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_712500/Run_000/psi_gs.mat','psi','Params','Transf','Observ'); +% Data = load('D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta40/HighN/aS_9.562000e+01_theta_040_phi_000_N_508333/Run_000/psi_gs.mat','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)); +IMG = nxy; + +% + +[psi, ratio, N] = computeBondOrderParameters(IMG); + +fprintf('Points: %d\n⟨|ψ₂|⟩ = %.3f, ⟨|ψ₄|⟩ = %.3f, ⟨|ψ₆|⟩ = %.3f\n', N, psi.psi2, psi.psi4, psi.psi6); +fprintf('(⟨|ψ₆|⟩ / ⟨|ψ₂|⟩) = %.3f\n', ratio); diff --git a/Data-Analyzer/extractHexaticOrder.m b/Data-Analyzer/extractHexaticOrder.m new file mode 100644 index 0000000..6c81b0c --- /dev/null +++ b/Data-Analyzer/extractHexaticOrder.m @@ -0,0 +1,98 @@ +function extractHexaticOrder(I) + % Input: + % I = 2D image matrix (grayscale, double) + % Output: + % Plots original image + hexatic order + histogram + + %% 1) Detect dots (local maxima) + I_smooth = imgaussfilt(I, 2); + bw = imregionalmax(I_smooth); + [y, x] = find(bw); + + if length(x) < 6 + error('Too few detected dots (%d) for meaningful hexatic order.', length(x)); + end + + %% 2) Compute hexatic order parameter psi6 + + % Delaunay triangulation to find neighbors + tri = delaunay(x,y); + N = length(x); + neighbors = cell(N,1); + for i = 1:size(tri,1) + neighbors{tri(i,1)} = unique([neighbors{tri(i,1)} tri(i,2) tri(i,3)]); + neighbors{tri(i,2)} = unique([neighbors{tri(i,2)} tri(i,1) tri(i,3)]); + neighbors{tri(i,3)} = unique([neighbors{tri(i,3)} tri(i,1) tri(i,2)]); + end + + psi6_vals = zeros(N,1); + for j = 1:N + nbrs = neighbors{j}; + nbrs(nbrs == j) = []; + if isempty(nbrs) + psi6_vals(j) = 0; + continue; + end + angles = atan2(y(nbrs)-y(j), x(nbrs)-x(j)); + psi6_vals(j) = abs(mean(exp(1i*6*angles))); + end + psi6_global = mean(psi6_vals); + + %% 3) Plotting + + figure('Name','Hexatic Order Analysis','NumberTitle','off', 'Units','normalized', 'Position',[0.2 0.2 0.6 0.7]); + t = tiledlayout(2,2, 'Padding','none', 'TileSpacing','compact'); + + % Top left: Original image with detected dots (square) + ax1 = nexttile(1); + imshow(I, []); + hold on; + scatter(x, y, 30, 'w', 'filled'); + hold off; + colormap(ax1, Helper.Colormaps.plasma()); % Example colormap for original image plot + cbar1 = colorbar(ax1); + cbar1.Label.Interpreter = 'latex'; + xlabel(ax1, '$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14); + ylabel(ax1, '$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14); + title(ax1, '$|\Psi(x,y)|^2$', 'Interpreter', 'latex', 'FontSize', 14); + axis(ax1, 'square'); + + % Top right: Histogram of local |psi6| (square) + ax2 = nexttile(2); + histogram(psi6_vals, 20, 'FaceColor', 'b'); + xlabel(ax2, '$|\psi_6|$', 'Interpreter', 'latex', 'FontSize', 14); + ylabel(ax2, 'Count'); + title(ax2, 'Histogram of Local Hexatic Order'); + xlim(ax2, [0 1]); + grid(ax2, 'on'); + axis(ax2, 'square'); + + % Bottom: Hexatic order magnitude on dots with bond orientation arrows (span 2 tiles) + ax3 = nexttile([1 2]); + scatter_handle = scatter(x, y, 50, psi6_vals, 'filled'); + colormap(ax3, 'jet'); % Different colormap for hexatic order magnitude + cbar3 = colorbar(ax3); + cbar3.Label.String = '|$\psi_6$|'; + cbar3.Label.Interpreter = 'latex'; + caxis(ax3, [0 1]); + axis(ax3, 'equal'); + axis(ax3, 'tight'); + title(ax3, sprintf('Local Hexatic Order |\\psi_6|, Global = %.3f', psi6_global)); + hold(ax3, 'on'); + + for j = 1:N + nbrs = neighbors{j}; + nbrs(nbrs == j) = []; + if isempty(nbrs) + continue; + end + angles = atan2(y(nbrs)-y(j), x(nbrs)-x(j)); + psi6_complex = mean(exp(1i*6*angles)); + local_angle = angle(psi6_complex)/6; + arrow_length = 5; + quiver(ax3, x(j), y(j), arrow_length*cos(local_angle), arrow_length*sin(local_angle), ... + 'k', 'MaxHeadSize', 2, 'AutoScale', 'off'); + end + hold(ax3, 'off'); + +end From 63904f2a95d2cf0fbf61b1c14ae9960ac173b76c Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Tue, 3 Jun 2025 16:48:11 +0200 Subject: [PATCH 3/4] 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 From 047a1b2d854993cda3eb58fd055cec2a01f65e30 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Sat, 7 Jun 2025 16:39:00 +0200 Subject: [PATCH 4/4] Added scripts to analyse simulation data, now possible to extract spectral contrast, weight and autocorrelation. --- .../analyzePhaseTransitionSimulation.m | 554 ++++++++++++++++++ Data-Analyzer/compareSpectralWeights.m | 55 +- Data-Analyzer/execution_scripts.m | 2 +- Data-Analyzer/extractAutocorrelation.m | 144 ----- Data-Analyzer/extractSpectralWeight.m | 81 ++- Data-Analyzer/matchTFRadii.m | 164 ++++++ Dipolar-Gas-Simulator/+Scripts/run_locally.m | 47 +- 7 files changed, 868 insertions(+), 179 deletions(-) create mode 100644 Data-Analyzer/analyzePhaseTransitionSimulation.m create mode 100644 Data-Analyzer/matchTFRadii.m diff --git a/Data-Analyzer/analyzePhaseTransitionSimulation.m b/Data-Analyzer/analyzePhaseTransitionSimulation.m new file mode 100644 index 0000000..dc94bb6 --- /dev/null +++ b/Data-Analyzer/analyzePhaseTransitionSimulation.m @@ -0,0 +1,554 @@ +%% Extract Images + +baseDir = 'D:/Results - Numerics/Data_Full3D/PhaseTransition/DTS/'; +JobNumber = 0; +runFolder = sprintf('Run_%03d', JobNumber); +movieFileName = 'DropletsToStripes.mp4'; % Output file name +datafileName = './DropletsToStripes.mat'; +reverseOrder = false; % Set this to true to reverse the theta ordering + +TitleString = 'Change across transition: Droplets To Stripes'; + +%% +baseDir = 'D:/Results - Numerics/Data_Full3D/PhaseTransition/STD/'; +JobNumber = 0; +runFolder = sprintf('Run_%03d', JobNumber); +movieFileName = 'StripesToDroplets.mp4'; % Output file name +datafileName = './StripesToDroplets.mat'; +reverseOrder = true; % Set this to true to reverse the theta ordering + + + TitleString = 'Change across transition: Stripes To Droplets'; + +%% +folderList = dir(baseDir); +isValid = [folderList.isdir] & ~ismember({folderList.name}, {'.', '..'}); +folderNames = {folderList(isValid).name}; +nimgs = numel(folderNames); + +% Extract theta values from folder names +PolarAngleVals = zeros(1, nimgs); +for k = 1:nimgs + tokens = regexp(folderNames{k}, 'theta_(\d{3})', 'tokens'); + if isempty(tokens) + warning('No theta found in folder name: %s', folderNames{k}); + PolarAngleVals(k) = NaN; + else + PolarAngleVals(k) = str2double(tokens{1}{1}); + end +end + +% Choose sort direction +sortDirection = 'ascend'; +if reverseOrder + sortDirection = 'descend'; +end + +% Sort folderNames based on polar angle +[~, sortIdx] = sort(PolarAngleVals, sortDirection); +folderNames = folderNames(sortIdx); +PolarAngleVals = PolarAngleVals(sortIdx); % Optional: if you still want sorted list + +imgs = cell(1, nimgs); +alphas = zeros(1, nimgs); + +for k = 1:nimgs + folderName = folderNames{k}; + SaveDirectory = fullfile(baseDir, folderName, runFolder); + + % Extract alpha (theta) again from folder name + tokens = regexp(folderName, 'theta_(\d{3})', 'tokens'); + alpha_val = str2double(tokens{1}{1}); + alphas(k) = alpha_val; + + matPath = fullfile(SaveDirectory, 'psi_gs.mat'); + if ~isfile(matPath) + warning('Missing psi_gs.mat in %s', SaveDirectory); + continue; + end + + try + Data = load(matPath, 'psi', 'Params', 'Transf', 'Observ'); + catch ME + warning('Failed to load %s: %s', matPath, ME.message); + continue; + end + + Params = Data.Params; + Transf = Data.Transf; + Observ = Data.Observ; + + psi = Data.psi; + if isgpuarray(psi) + psi = gather(psi); + end + if isgpuarray(Observ.residual) + Observ.residual = gather(Observ.residual); + end + + % Axes and projection + x = Transf.x * Params.l0 * 1e6; + y = Transf.y * Params.l0 * 1e6; + z = Transf.z * Params.l0 * 1e6; + dx = x(2)-x(1); dy = y(2)-y(1); dz = z(2)-z(1); + + % Calculate frequency increment (frequency axes) + Nx = length(x); % grid size along X + Ny = length(y); % grid size along Y + dx = mean(diff(x)); % real space increment in the X direction (in micrometers) + dy = mean(diff(y)); % real space increment in the Y direction (in micrometers) + dvx = 1 / (Nx * dx); % reciprocal space increment in the X direction (in micrometers^-1) + dvy = 1 / (Ny * dy); % reciprocal space increment in the Y direction (in micrometers^-1) + + % Create the frequency axes + vx = (-Nx/2:Nx/2-1) * dvx; % Frequency axis in X (micrometers^-1) + vy = (-Ny/2:Ny/2-1) * dvy; % Frequency axis in Y (micrometers^-1) + + % Calculate maximum frequencies + % kx_max = pi / dx; + % ky_max = pi / dy; + + % Generate reciprocal axes + % kx = linspace(-kx_max, kx_max * (Nx-2)/Nx, Nx); + % ky = linspace(-ky_max, ky_max * (Ny-2)/Ny, Ny); + + % Create the Wavenumber axes + kx = 2*pi*vx; % Wavenumber axis in X + ky = 2*pi*vy; % Wavenumber axis in Y + + n = abs(psi).^2; + nxy = squeeze(trapz(n * dz, 3)); + + imgs{k} = nxy; +end + +%% Analyze Images + + +makeMovie = true; % Set to false to disable movie creation +font = 'Bahnschrift'; + +skipPreprocessing = true; +skipMasking = true; +skipIntensityThresholding = true; +skipBinarization = true; + +% Run Fourier analysis over images + +fft_imgs = cell(1, nimgs); +spectral_contrast = zeros(1, nimgs); +spectral_weight = zeros(1, nimgs); +g2_all = cell(1, nimgs); +theta_values_all = cell(1, nimgs); + +N_bins = 180; +Threshold = 25; +Sigma = 2; + +if makeMovie + % Create VideoWriter object for movie + videoFile = VideoWriter(movieFileName, 'MPEG-4'); + videoFile.Quality = 100; % Set quality to maximum (0–100) + videoFile.FrameRate = 2; % Set the frame rate (frames per second) + open(videoFile); % Open the video file to write +end + +% Display the cropped image +for k = 1:nimgs + IMG = imgs{k}; + [IMGFFT, IMGPR] = computeFourierTransform(IMG, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization); + [theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(IMGFFT, 10, 35, N_bins, Threshold, Sigma); + + g2 = zeros(1, N_bins); % Preallocate + + for dtheta = 0:N_bins-1 + profile = S_theta; + profile_shifted = circshift(profile, -dtheta, 2); + + num = mean(profile .* profile_shifted); + denom = mean(profile)^2; + + g2(dtheta+1) = num / denom - 1; + end + + g2_all{k} = g2; + theta_values_all{k} = theta_vals; + + figure(1); + clf + set(gcf,'Position',[500 100 1000 800]) + t = tiledlayout(2, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 1x4 grid + + y_min = min(y); + y_max = max(y); + x_min = min(x); + x_max = max(x); + + % Display the cropped OD image + ax1 = nexttile; + imagesc(x, y, IMG') + % Define normalized positions (relative to axis limits) + x_offset = 0.025; % 5% offset from the edges + y_offset = 0.025; % 5% offset from the edges + % Top-right corner (normalized axis coordinates) + hText = text(1 - x_offset, 1 - y_offset, ['Angle = ', num2str(alphas(k), '%.1f')], ... + 'Color', 'white', 'FontWeight', 'bold', 'Interpreter', 'tex', 'FontSize', 20, 'Units', 'normalized', 'HorizontalAlignment', 'right', 'VerticalAlignment', 'top'); + axis square; + hcb = colorbar; + colormap(ax1, 'jet'); + set(gca, 'FontSize', 14); % For tick labels only + hL = ylabel(hcb, 'Optical Density'); + set(hL,'Rotation',-90); + set(gca,'YDir','normal') + % set(gca, 'YTick', linspace(y_min, y_max, 5)); % Define y ticks + % set(gca, 'YTickLabel', flip(linspace(y_min, y_max, 5))); % Flip only the labels + hXLabel = xlabel('x (pixels)', 'Interpreter', 'tex'); + hYLabel = ylabel('y (pixels)', 'Interpreter', 'tex'); + hTitle = title('Density', 'Interpreter', 'tex'); + set([hXLabel, hYLabel, hL, hText], 'FontName', font) + set([hXLabel, hYLabel, hL], 'FontSize', 14) + set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title + + % Plot the power spectrum + ax2 = nexttile; + [rows, cols] = size(IMGFFT); + zoom_size = 50; % Zoomed-in region around center + mid_x = floor(cols/2); + mid_y = floor(rows/2); + fft_imgs{k} = IMGFFT(mid_y-zoom_size:mid_y+zoom_size, mid_x-zoom_size:mid_x+zoom_size); + imagesc(log(1 + abs(fft_imgs{k}).^2)); + % Define normalized positions (relative to axis limits) + x_offset = 0.025; % 5% offset from the edges + y_offset = 0.025; % 5% offset from the edges + axis square; + hcb = colorbar; + colormap(ax2, 'jet'); + set(gca, 'FontSize', 14); % For tick labels only + set(gca,'YDir','normal') + hXLabel = xlabel('k_x', 'Interpreter', 'tex'); + hYLabel = ylabel('k_y', 'Interpreter', 'tex'); + hTitle = title('Power Spectrum - S(k_x,k_y)', 'Interpreter', 'tex'); + set([hXLabel, hYLabel, hText], 'FontName', font) + set([hXLabel, hYLabel], 'FontSize', 14) + set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title + + % Plot the angular distribution + nexttile + spectral_contrast(k) = computeSpectralContrast(fft_imgs{k}, 10, 25, Threshold); + [theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(fft_imgs{k}, 10, 25, N_bins, Threshold, Sigma); + spectral_weight(k) = trapz(theta_vals, S_theta); + plot(theta_vals/pi, S_theta,'Linewidth',2); + axis square; + set(gca, 'FontSize', 14); % For tick labels only + hXLabel = xlabel('\theta/\pi [rad]', 'Interpreter', 'tex'); + hYLabel = ylabel('Normalized magnitude (a.u.)', 'Interpreter', 'tex'); + hTitle = title('Angular Spectral Distribution - S(\theta)', 'Interpreter', 'tex'); + set([hXLabel, hYLabel, hText], 'FontName', font) + set([hXLabel, hYLabel], 'FontSize', 14) + set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title + grid on + + nexttile + plot(theta_vals/pi, g2, 'o-', 'LineWidth', 1.2, 'MarkerSize', 5); + set(gca, 'FontSize', 14); + ylim([-1.5 3.0]); % Set y-axis limits here + hXLabel = xlabel('$\delta\theta / \pi$', 'Interpreter', 'latex'); + hYLabel = ylabel('$g^{(2)}(\delta\theta)$', 'Interpreter', 'latex'); + hTitle = title('Autocorrelation', 'Interpreter', 'tex'); + set([hXLabel, hYLabel], 'FontName', font) + set([hXLabel, hYLabel], 'FontSize', 14) + set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title + grid on; + + if makeMovie + frame = getframe(gcf); + writeVideo(videoFile, frame); + else + pause(0.5); % Only pause when not recording + end + +end + +if makeMovie + close(videoFile); + disp(['Movie saved to ', movieFileName]); +end + +%% Track across the transition + +figure(2); +set(gcf,'Position',[100 100 950 750]) +plot(alphas, spectral_contrast, 'o--', ... + 'LineWidth', 1.5, 'MarkerSize', 6); +set(gca, 'FontSize', 14); % For tick labels only +hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex'); +% hXLabel = xlabel('B_z (G)', 'Interpreter', 'tex'); +hYLabel = ylabel('Spectral Contrast', 'Interpreter', 'tex'); +hTitle = title(TitleString, 'Interpreter', 'tex'); +set([hXLabel, hYLabel], 'FontName', font) +set([hXLabel, hYLabel], 'FontSize', 14) +set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title +grid on + +figure(3); +set(gcf,'Position',[100 100 950 750]) +plot(alphas, spectral_weight, 'o--', ... + 'LineWidth', 1.5, 'MarkerSize', 6); +set(gca, 'FontSize', 14); % For tick labels only +hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex'); +% hXLabel = xlabel('B_z (G)', 'Interpreter', 'tex'); +hYLabel = ylabel('Spectral Weight', 'Interpreter', 'tex'); +hTitle = title(TitleString, 'Interpreter', 'tex'); +set([hXLabel, hYLabel], 'FontName', font) +set([hXLabel, hYLabel], 'FontSize', 14) +set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title +grid on + +save(datafileName, 'alphas', 'spectral_contrast', 'spectral_weight'); + +figure(4); +clf; +set(gcf,'Position',[100 100 950 750]) +hold on; + +% Reconstruct theta axis from any one of the stored values +theta_vals = theta_values_all{1}; % assuming it's in radians + +legend_entries = cell(nimgs, 1); + +% Generate a colormap with enough unique colors +cmap = sky(nimgs); % You can also try 'jet', 'turbo', 'hot', etc. + +for i = 1:nimgs + plot(theta_vals/pi, g2_all{i}, ... + 'o-', 'Color', cmap(i,:), 'LineWidth', 1.2, ... + 'MarkerSize', 5); + legend_entries{i} = sprintf('$\\alpha = %g^\\circ$', alphas(i)); +end + +ylim([-1.5 3.0]); % Set y-axis limits here +set(gca, 'FontSize', 14); +hXLabel = xlabel('$\delta\theta / \pi$', 'Interpreter', 'latex'); +hYLabel = ylabel('$g^{(2)}(\delta\theta)$', 'Interpreter', 'latex'); +hTitle = title(TitleString, 'Interpreter', 'tex'); +legend(legend_entries, 'Interpreter', 'latex', 'Location', 'bestoutside'); +set([hXLabel, hYLabel], 'FontName', font) +set([hXLabel, hYLabel], 'FontSize', 14) +set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title +grid on; + +%% Track across the transition + +set(0,'defaulttextInterpreter','latex') +set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex'); + +format long + +font = 'Bahnschrift'; + +% Load data +Data = load('./DropletsToStripes.mat', 'alphas', 'spectral_contrast', 'spectral_weight'); +dts_alphas = Data.alphas; +dts_sc = Data.spectral_contrast; +dts_sw = Data.spectral_weight; + +Data = load('./StripesToDroplets.mat', 'alphas', 'spectral_contrast', 'spectral_weight'); +std_alphas = Data.alphas; +std_sc = Data.spectral_contrast; +std_sw = Data.spectral_weight; + +% Normalize dts data +dts_min = min(dts_sw); +dts_max = max(dts_sw); +dts_range = dts_max - dts_min; +dts_sf_norm = (dts_sw - dts_min) / dts_range; + +% Normalize std data +std_min = min(std_sw); +std_max = max(std_sw); +std_range = std_max - std_min; +std_sf_norm = (std_sw - std_min) / std_range; + +figure(5); +set(gcf,'Position',[100 100 950 750]) +plot(dts_alphas, dts_sc, 'o--', 'LineWidth', 1.5, 'MarkerSize', 6, 'DisplayName' , 'Droplets to Stripes'); +hold on +plot(std_alphas, std_sc, 'o--', 'LineWidth', 1.5, 'MarkerSize', 6, 'DisplayName' , 'Stripes to Droplets'); +set(gca, 'FontSize', 14); % For tick labels only +hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex'); +hYLabel = ylabel('Spectral Contrast', 'Interpreter', 'tex'); +hTitle = title('Change across transition', 'Interpreter', 'tex'); +legend +set([hXLabel, hYLabel], 'FontName', font) +set([hXLabel, hYLabel], 'FontSize', 14) +set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title +grid on + +figure(6); +set(gcf,'Position',[100 100 950 750]) +plot(dts_alphas, dts_sw, 'o--', 'LineWidth', 1.5, 'MarkerSize', 6, 'DisplayName' , 'Droplets to Stripes'); +hold on +plot(std_alphas, std_sw, 'o--', 'LineWidth', 1.5, 'MarkerSize', 6, 'DisplayName' , 'Stripes to Droplets'); +set(gca, 'FontSize', 14); % For tick labels only +hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex'); +hYLabel = ylabel('Spectral Weight', 'Interpreter', 'tex'); +hTitle = title('Change across transition', 'Interpreter', 'tex'); +legend +set([hXLabel, hYLabel], 'FontName', font) +set([hXLabel, hYLabel], 'FontSize', 14) +set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title +grid on + +%% +function [IMGFFT, IMGPR] = computeFourierTransform(I, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization) + % computeFourierSpectrum - Computes the 2D Fourier power spectrum + % of binarized and enhanced lattice image features, with optional central mask. + % + % Inputs: + % I - Grayscale or RGB image matrix + % + % Output: + % F_mag - 2D Fourier power spectrum (shifted) + + if ~skipPreprocessing + % Preprocessing: Denoise + filtered = imgaussfilt(I, 10); + IMGPR = I - filtered; % adjust sigma as needed + else + IMGPR = I; + end + + if ~skipMasking + [rows, cols] = size(IMGPR); + [X, Y] = meshgrid(1:cols, 1:rows); + % Elliptical mask parameters + cx = cols / 2; + cy = rows / 2; + + % Shifted coordinates + x = X - cx; + y = Y - cy; + + % Ellipse semi-axes + rx = 0.4 * cols; + ry = 0.2 * rows; + + % Rotation angle in degrees -> radians + theta_deg = 30; % Adjust as needed + theta = deg2rad(theta_deg); + + % Rotated ellipse equation + cos_t = cos(theta); + sin_t = sin(theta); + + x_rot = (x * cos_t + y * sin_t); + y_rot = (-x * sin_t + y * cos_t); + + ellipseMask = (x_rot.^2) / rx^2 + (y_rot.^2) / ry^2 <= 1; + + % Apply cutout mask + IMGPR = IMGPR .* ellipseMask; + end + + if ~skipIntensityThresholding + % Apply global intensity threshold mask + intensity_thresh = 0.20; + intensity_mask = IMGPR > intensity_thresh; + IMGPR = IMGPR .* intensity_mask; + end + + if ~skipBinarization + % Adaptive binarization and cleanup + IMGPR = imbinarize(IMGPR, 'adaptive', 'Sensitivity', 0.0); + IMGPR = imdilate(IMGPR, strel('disk', 2)); + IMGPR = imerode(IMGPR, strel('disk', 1)); + IMGPR = imfill(IMGPR, 'holes'); + F = fft2(double(IMGPR)); % Compute 2D Fourier Transform + IMGFFT = abs(fftshift(F))'; % Shift zero frequency to center + else + F = fft2(double(IMGPR)); % Compute 2D Fourier Transform + IMGFFT = abs(fftshift(F))'; % Shift zero frequency to center + end +end + +function [theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(IMGFFT, r_min, r_max, num_bins, threshold, sigma) + % Apply threshold to isolate strong peaks + IMGFFT(IMGFFT < threshold) = 0; + + % Prepare polar coordinates + [ny, nx] = size(IMGFFT); + [X, Y] = meshgrid(1:nx, 1:ny); + cx = ceil(nx/2); + cy = ceil(ny/2); + R = sqrt((X - cx).^2 + (Y - cy).^2); + Theta = atan2(Y - cy, X - cx); % range [-pi, pi] + + % Choose radial band + radial_mask = (R >= r_min) & (R <= r_max); + + % Initialize the angular structure factor array + S_theta = zeros(1, num_bins); % Pre-allocate for 180 angle bins + % Define the angle values for the x-axis + theta_vals = linspace(0, pi, num_bins); + + % Loop through each angle bin + for i = 1:num_bins + angle_start = (i-1) * pi / num_bins; + angle_end = i * pi / num_bins; + + % Define a mask for the given angle range + angle_mask = (Theta >= angle_start & Theta < angle_end); + + bin_mask = radial_mask & angle_mask; + + % Extract the Fourier components for the given angle + fft_angle = IMGFFT .* bin_mask; + + % Integrate the Fourier components over the radius at the angle + S_theta(i) = sum(sum(abs(fft_angle).^2)); % sum of squared magnitudes + end + + % Create a 1D Gaussian kernel + half_width = ceil(3 * sigma); + x = -half_width:half_width; + gauss_kernel = exp(-x.^2 / (2 * sigma^2)); + gauss_kernel = gauss_kernel / sum(gauss_kernel); % normalize + + % Apply convolution (circular padding to preserve periodicity) + S_theta = conv([S_theta(end-half_width+1:end), S_theta, S_theta(1:half_width)], gauss_kernel, 'same'); + S_theta = S_theta(half_width+1:end-half_width); % crop back to original size + + % Normalize to 1 + S_theta = S_theta / max(S_theta); +end + +function contrast = computeSpectralContrast(IMGFFT, r_min, r_max, threshold) + % Apply threshold to isolate strong peaks + IMGFFT(IMGFFT < threshold) = 0; + + % Prepare polar coordinates + [ny, nx] = size(IMGFFT); + [X, Y] = meshgrid(1:nx, 1:ny); + cx = ceil(nx/2); + cy = ceil(ny/2); + R = sqrt((X - cx).^2 + (Y - cy).^2); + + % Ring region (annulus) mask + ring_mask = (R >= r_min) & (R <= r_max); + + % Squared magnitude in the ring + ring_power = abs(IMGFFT).^2 .* ring_mask; + + % Maximum power in the ring + ring_max = max(ring_power(:)); + + % Power at the DC component + dc_power = abs(IMGFFT(cy, cx))^2; + + % Avoid division by zero + if dc_power == 0 + contrast = Inf; % or NaN or 0, depending on how you want to handle this + else + contrast = ring_max / dc_power; + end +end diff --git a/Data-Analyzer/compareSpectralWeights.m b/Data-Analyzer/compareSpectralWeights.m index 94f075e..fd1e852 100644 --- a/Data-Analyzer/compareSpectralWeights.m +++ b/Data-Analyzer/compareSpectralWeights.m @@ -8,42 +8,63 @@ format long font = 'Bahnschrift'; % Load data -Data = load('C:/Users/Karthik/Documents/GitRepositories/Calculations/Data-Analyzer/B2.45G/DropletsToStripes.mat', 'unique_scan_parameter_values', 'mean_sf', 'stderr_sf'); +Data = load('D:/Results - Experiment/B2.45G/DropletsToStripes.mat', 'unique_scan_parameter_values', 'mean_sc', 'stderr_sc', 'mean_sw', 'stderr_sw'); dts_scan_parameter_values = Data.unique_scan_parameter_values; -dts_mean_sf = Data.mean_sf; -dts_stderr_sf = Data.stderr_sf; +dts_mean_sc = Data.mean_sc; +dts_stderr_sc = Data.stderr_sc; +dts_mean_sw = Data.mean_sw; +dts_stderr_sw = Data.stderr_sw; -Data = load('C:/Users/Karthik/Documents/GitRepositories/Calculations/Data-Analyzer/B2.45G/StripesToDroplets.mat', 'unique_scan_parameter_values', 'mean_sf', 'stderr_sf'); +Data = load('D:/Results - Experiment/B2.45G/StripesToDroplets.mat', 'unique_scan_parameter_values', 'mean_sc', 'stderr_sc', 'mean_sw', 'stderr_sw'); std_scan_parameter_values = Data.unique_scan_parameter_values; -std_mean_sf = Data.mean_sf; -std_stderr_sf = Data.stderr_sf; +std_mean_sw = Data.mean_sw; +std_stderr_sw = Data.stderr_sw; +std_mean_sc = Data.mean_sc; +std_stderr_sc = Data.stderr_sc; % Normalize dts data -dts_min = min(dts_mean_sf); -dts_max = max(dts_mean_sf); +dts_min = min(dts_mean_sw); +dts_max = max(dts_mean_sw); dts_range = dts_max - dts_min; -dts_mean_sf_norm = (dts_mean_sf - dts_min) / dts_range; -dts_stderr_sf_norm = dts_stderr_sf / dts_range; +dts_mean_sw_norm = (dts_mean_sw - dts_min) / dts_range; +dts_stderr_sw_norm = dts_stderr_sw / dts_range; % Normalize std data -std_min = min(std_mean_sf); -std_max = max(std_mean_sf); +std_min = min(std_mean_sw); +std_max = max(std_mean_sw); std_range = std_max - std_min; -std_mean_sf_norm = (std_mean_sf - std_min) / std_range; -std_stderr_sf_norm = std_stderr_sf / std_range; +std_mean_sw_norm = (std_mean_sw - std_min) / std_range; +std_stderr_sw_norm = std_stderr_sw / std_range; figure(1); set(gcf,'Position',[100 100 950 750]) -errorbar(dts_scan_parameter_values, dts_mean_sf_norm, dts_stderr_sf_norm, 'o--', ... +errorbar(dts_scan_parameter_values, dts_mean_sc, dts_stderr_sc, 'o--', ... 'LineWidth', 1.5, 'MarkerSize', 6, 'CapSize', 5, 'DisplayName' , 'Droplets to Stripes'); hold on -errorbar(std_scan_parameter_values, std_mean_sf_norm, std_stderr_sf_norm, 'o--', ... +errorbar(std_scan_parameter_values, std_mean_sc, std_stderr_sc, 'o--', ... 'LineWidth', 1.5, 'MarkerSize', 6, 'CapSize', 5, 'DisplayName', 'Stripes to Droplets'); set(gca, 'FontSize', 14); % For tick labels only hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex'); -hYLabel = ylabel('Normalized Spectral Weight', 'Interpreter', 'tex'); +hYLabel = ylabel('Spectral Contrast', 'Interpreter', 'tex'); +hTitle = title('B = 2.45 G', 'Interpreter', 'tex'); +legend +set([hXLabel, hYLabel], 'FontName', font) +set([hXLabel, hYLabel], 'FontSize', 14) +set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title +grid on + +figure(2); +set(gcf,'Position',[100 100 950 750]) +errorbar(dts_scan_parameter_values, dts_mean_sw, dts_stderr_sw, 'o--', ... + 'LineWidth', 1.5, 'MarkerSize', 6, 'CapSize', 5, 'DisplayName' , 'Droplets to Stripes'); +hold on +errorbar(std_scan_parameter_values, std_mean_sw, std_stderr_sw, 'o--', ... + 'LineWidth', 1.5, 'MarkerSize', 6, 'CapSize', 5, 'DisplayName', 'Stripes to Droplets'); +set(gca, 'FontSize', 14); % For tick labels only +hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex'); +hYLabel = ylabel('Spectral Weight', 'Interpreter', 'tex'); hTitle = title('B = 2.45 G', 'Interpreter', 'tex'); legend set([hXLabel, hYLabel], 'FontName', font) diff --git a/Data-Analyzer/execution_scripts.m b/Data-Analyzer/execution_scripts.m index 653cfda..2fb647e 100644 --- a/Data-Analyzer/execution_scripts.m +++ b/Data-Analyzer/execution_scripts.m @@ -60,4 +60,4 @@ IMG = nxy; [psi, ratio, N] = computeBondOrderParameters(IMG); fprintf('Points: %d\n⟨|ψ₂|⟩ = %.3f, ⟨|ψ₄|⟩ = %.3f, ⟨|ψ₆|⟩ = %.3f\n', N, psi.psi2, psi.psi4, psi.psi6); -fprintf('(⟨|ψ₆|⟩ / ⟨|ψ₂|⟩) = %.3f\n', ratio); +fprintf('(⟨|ψ₆|⟩ / ⟨|ψ₂|⟩) = %.3f\n', ratio); \ No newline at end of file diff --git a/Data-Analyzer/extractAutocorrelation.m b/Data-Analyzer/extractAutocorrelation.m index 7e66f76..073ab11 100644 --- a/Data-Analyzer/extractAutocorrelation.m +++ b/Data-Analyzer/extractAutocorrelation.m @@ -198,150 +198,6 @@ set([hXLabel, hYLabel], 'FontSize', 14) set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title grid on; -%% Extract g2 from simulation data - -Data = load('E:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_712500/Run_000/psi_gs.mat','psi','Params','Transf','Observ'); - -% Data = load('E:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta40/HighN/aS_9.562000e+01_theta_040_phi_000_N_508333/Run_000/psi_gs.mat','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; - -dx = x(2)-x(1); dy = y(2)-y(1); dz = z(2)-z(1); - -% Calculate frequency increment (frequency axes) -Nx = length(x); % grid size along X -Ny = length(y); % grid size along Y -dx = mean(diff(x)); % real space increment in the X direction (in micrometers) -dy = mean(diff(y)); % real space increment in the Y direction (in micrometers) -dvx = 1 / (Nx * dx); % reciprocal space increment in the X direction (in micrometers^-1) -dvy = 1 / (Ny * dy); % reciprocal space increment in the Y direction (in micrometers^-1) - -% Create the frequency axes -vx = (-Nx/2:Nx/2-1) * dvx; % Frequency axis in X (micrometers^-1) -vy = (-Ny/2:Ny/2-1) * dvy; % Frequency axis in Y (micrometers^-1) - -% Calculate maximum frequencies -% kx_max = pi / dx; -% ky_max = pi / dy; - -% Generate reciprocal axes -% kx = linspace(-kx_max, kx_max * (Nx-2)/Nx, Nx); -% ky = linspace(-ky_max, ky_max * (Ny-2)/Ny, Ny); - -% Create the Wavenumber axes -kx = 2*pi*vx; % Wavenumber axis in X -ky = 2*pi*vy; % Wavenumber axis in Y - -% Compute probability density |psi|^2 -n = abs(psi).^2; - -nxz = squeeze(trapz(n*dy,2)); -nyz = squeeze(trapz(n*dx,1)); -nxy = squeeze(trapz(n*dz,3)); - -skipPreprocessing = true; -skipMasking = true; -skipIntensityThresholding = true; -skipBinarization = true; - -font = 'Bahnschrift'; -% Extract g2 - -N_bins = 90; -Threshold = 75; -Sigma = 2; - -IMG = nxy; - -[IMGFFT, IMGPR] = computeFourierTransform(IMG, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization); - -[theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(IMGFFT, 10, 35, N_bins, Threshold, Sigma); - -g2_all = zeros(1, N_bins); % Preallocate - -for dtheta = 0:N_bins-1 - profile = S_theta; - profile_shifted = circshift(profile, -dtheta, 2); - - num = mean(profile .* profile_shifted); - denom = mean(profile)^2; - - g2_all(dtheta+1) = num / denom - 1; -end - -figure(2); -clf -set(gcf,'Position',[500 100 1000 800]) -t = tiledlayout(2, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 1x4 grid - -% Display the cropped OD image -nexttile -plotxy = pcolor(x,y,IMG'); -set(plotxy, 'EdgeColor', 'none'); -cbar1 = colorbar; -cbar1.Label.Interpreter = 'latex'; -colormap(gca, Helper.Colormaps.plasma()) -xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14) -ylabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14) -title('$|\Psi(x,y)|^2$', 'Interpreter', 'latex', 'FontSize', 14) - -% Plot the power spectrum -nexttile; -imagesc(kx, ky, log(1 + abs(IMGFFT).^2)); -axis square; -hcb = colorbar; -colormap(gca, Helper.Colormaps.plasma()) -set(gca, 'FontSize', 14); % For tick labels only -set(gca,'YDir','normal') -hXLabel = xlabel('k_x', 'Interpreter', 'tex'); -hYLabel = ylabel('k_y', 'Interpreter', 'tex'); -hTitle = title('Power Spectrum - S(k_x,k_y)', 'Interpreter', 'tex'); -set([hXLabel, hYLabel], 'FontName', font) -set([hXLabel, hYLabel], 'FontSize', 14) -set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title - -% Plot the angular distribution -nexttile -plot(theta_vals/pi, S_theta,'Linewidth',2); -set(gca, 'FontSize', 14); % For tick labels only -hXLabel = xlabel('\theta/\pi [rad]', 'Interpreter', 'tex'); -hYLabel = ylabel('Normalized magnitude (a.u.)', 'Interpreter', 'tex'); -hTitle = title('Angular Spectral Distribution - S(\theta)', 'Interpreter', 'tex'); -set([hXLabel, hYLabel], 'FontName', font) -set([hXLabel, hYLabel], 'FontSize', 14) -set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title -grid on - -nexttile -plot(theta_vals/pi, g2_all, 'o-', 'LineWidth', 1.2, 'MarkerSize', 5); -set(gca, 'FontSize', 14); -ylim([-1.5 3.0]); % Set y-axis limits here -hXLabel = xlabel('$\delta\theta / \pi$', 'Interpreter', 'latex'); -hYLabel = ylabel('$g^{(2)}(\delta\theta)$', 'Interpreter', 'latex'); -hTitle = title('Autocorrelation', 'Interpreter', 'tex'); -set([hXLabel, hYLabel], 'FontName', font) -set([hXLabel, hYLabel], 'FontSize', 14) -set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title -grid on; - %% Helper Functions function [IMGFFT, IMGPR] = computeFourierTransform(I, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization) % computeFourierSpectrum - Computes the 2D Fourier power spectrum diff --git a/Data-Analyzer/extractSpectralWeight.m b/Data-Analyzer/extractSpectralWeight.m index 526838a..accaf48 100644 --- a/Data-Analyzer/extractSpectralWeight.m +++ b/Data-Analyzer/extractSpectralWeight.m @@ -6,7 +6,7 @@ groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1_Axi folderPath = "D:/Data - Experiment/2025/05/22/"; -run = '0078'; +run = '0079'; folderPath = strcat(folderPath, run); @@ -25,7 +25,9 @@ scan_parameter = 'rot_mag_fin_pol_angle'; scan_parameter_text = 'Angle = '; % scan_parameter_text = 'BField = '; -font = 'Bahnschrift'; +savefodlerPath = 'D:/Results - Experiment/B2.45G/'; +savefileName = 'StripesToDroplets.mat'; +font = 'Bahnschrift'; skipPreprocessing = true; skipMasking = true; @@ -96,6 +98,7 @@ end %% Run Fourier analysis over images fft_imgs = cell(1, nimgs); +spectral_contrast = zeros(1, nimgs); spectral_weight = zeros(1, nimgs); N_bins = 180; @@ -197,6 +200,7 @@ for k = 1:N_shots % Plot the angular distribution nexttile + spectral_contrast(k) = computeSpectralContrast(fft_imgs{k}, 10, 25, Threshold); [theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(fft_imgs{k}, 10, 20, N_bins, Threshold, Sigma); spectral_weight(k) = trapz(theta_vals, S_theta); plot(theta_vals/pi, S_theta,'Linewidth',2); @@ -220,25 +224,51 @@ end % Close the video file close(videoFile); -%% Track spectral weight across the transition +%% Track across the transition % Assuming scan_parameter_values and spectral_weight are column vectors (or row vectors of same length) [unique_scan_parameter_values, ~, idx] = unique(scan_parameter_values); % Preallocate arrays -mean_sf = zeros(size(unique_scan_parameter_values)); -stderr_sf = zeros(size(unique_scan_parameter_values)); +mean_sc = zeros(size(unique_scan_parameter_values)); +stderr_sc = zeros(size(unique_scan_parameter_values)); % Loop through each unique theta and compute mean and standard error for i = 1:length(unique_scan_parameter_values) - group_vals = spectral_weight(idx == i); - mean_sf(i) = mean(group_vals); - stderr_sf(i) = std(group_vals) / sqrt(length(group_vals)); % standard error = std / sqrt(N) + group_vals = spectral_contrast(idx == i); + mean_sc(i) = mean(group_vals); + stderr_sc(i) = std(group_vals) / sqrt(length(group_vals)); % standard error = std / sqrt(N) end figure(2); set(gcf,'Position',[100 100 950 750]) -errorbar(unique_scan_parameter_values, mean_sf, stderr_sf, 'o--', ... +errorbar(unique_scan_parameter_values, mean_sc, stderr_sc, 'o--', ... + 'LineWidth', 1.5, 'MarkerSize', 6, 'CapSize', 5); +set(gca, 'FontSize', 14); % For tick labels only +hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex'); +% hXLabel = xlabel('B_z (G)', 'Interpreter', 'tex'); +hYLabel = ylabel('Spectral Contrast', 'Interpreter', 'tex'); +hTitle = title('Change across transition', 'Interpreter', 'tex'); +set([hXLabel, hYLabel], 'FontName', font) +set([hXLabel, hYLabel], 'FontSize', 14) +set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title +grid on + + +% Preallocate arrays +mean_sw = zeros(size(unique_scan_parameter_values)); +stderr_sw = zeros(size(unique_scan_parameter_values)); + +% Loop through each unique theta and compute mean and standard error +for i = 1:length(unique_scan_parameter_values) + group_vals = spectral_weight(idx == i); + mean_sw(i) = mean(group_vals); + stderr_sw(i) = std(group_vals) / sqrt(length(group_vals)); % standard error = std / sqrt(N) +end + +figure(3); +set(gcf,'Position',[100 100 950 750]) +errorbar(unique_scan_parameter_values, mean_sw, stderr_sw, 'o--', ... 'LineWidth', 1.5, 'MarkerSize', 6, 'CapSize', 5); set(gca, 'FontSize', 14); % For tick labels only hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex'); @@ -250,6 +280,8 @@ set([hXLabel, hYLabel], 'FontSize', 14) set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title grid on +save([savefolderPath savefileName], 'unique_scan_parameter_values', 'mean_sc', 'stderr_sc', 'mean_sw', 'stderr_sw'); + %% k-means Clustering % Reshape to column vector @@ -433,6 +465,37 @@ function [theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(IM S_theta = S_theta / max(S_theta); end +function contrast = computeSpectralContrast(IMGFFT, r_min, r_max, threshold) + % Apply threshold to isolate strong peaks + IMGFFT(IMGFFT < threshold) = 0; + + % Prepare polar coordinates + [ny, nx] = size(IMGFFT); + [X, Y] = meshgrid(1:nx, 1:ny); + cx = ceil(nx/2); + cy = ceil(ny/2); + R = sqrt((X - cx).^2 + (Y - cy).^2); + + % Ring region (annulus) mask + ring_mask = (R >= r_min) & (R <= r_max); + + % Squared magnitude in the ring + ring_power = abs(IMGFFT).^2 .* ring_mask; + + % Maximum power in the ring + ring_max = max(ring_power(:)); + + % Power at the DC component + dc_power = abs(IMGFFT(cy, cx))^2; + + % Avoid division by zero + if dc_power == 0 + contrast = Inf; % or NaN or 0, depending on how you want to handle this + else + contrast = ring_max / dc_power; + end +end + function ret = getBkgOffsetFromCorners(img, x_fraction, y_fraction) % image must be a 2D numerical array [dim1, dim2] = size(img); diff --git a/Data-Analyzer/matchTFRadii.m b/Data-Analyzer/matchTFRadii.m new file mode 100644 index 0000000..8b5d4d2 --- /dev/null +++ b/Data-Analyzer/matchTFRadii.m @@ -0,0 +1,164 @@ +% ---------------------------- +% Experimental data +% ---------------------------- +PixelSize = 5.86; % microns + +AtomNumbers = [9.14950197, 8.6845907 , 8.82521245, 8.5899089 , 8.21675841, ... + 8.96234044, 8.8636914 , 8.70332154, 8.82930706, 8.91869919, ... + 8.58553165, 8.73391981, 8.71943552, 8.11717678, 8.59490351, ... + 8.57514491, 8.81628891, 8.37211343, 8.76077699, 8.71297796, ... + 9.17634469, 8.81424285, 8.61176745, 8.40555897, 8.97137861, ... + 8.88393124, 8.66625724, 8.30688943, 9.02338201, 8.57729816, ... + 8.50333458, 8.67617084, 8.8936879 , 9.02031475, 8.98459233, ... + 8.76525048, 8.76801503, 8.58302559, 8.4617431 , 8.74479855, ... + 8.83882896, 8.69091377, 8.79282459, 8.51785483, 8.75629649, ... + 8.58994308, 8.36816564, 9.2429294 , 8.6583425 , 8.55827961]; + +EstimatedAtomNumber = mean(AtomNumbers) * 1E4; + +TF_Radii_X_pixels = [48.44308968, 46.01326593, 46.45950681, 46.41644117, 45.56176919, ... + 46.60816438, 46.85307478, 47.61086543, 46.66687703, 46.17986721, ... + 46.67877165, 46.07789481, 46.42285497, 46.22167708, 45.95144492, ... + 47.05400117, 49.03005788, 45.84588659, 46.85742777, 45.9824117 , ... + 47.14731188, 47.4984484 , 45.9055646 , 47.31804553, 47.52321888, ... + 47.76823968, 46.459749 , 45.4498851 , 45.38339308, 46.68736642, ... + 45.76607233, 48.1796053 , 46.94291541, 47.54092708, 48.26130406, ... + 47.44092616, 48.73463214, 46.39356452, 48.74120217, 45.57014182, ... + 47.56467835, 46.62867035, 46.62322802, 46.03032919, 44.78559832, ... + 46.31282562, 46.83537518, 47.68015029, 47.71093571, 47.34079816]; + +TF_Radii_Y_pixels = [113.52610841, 113.68862761, 112.84031747, 114.22062324, ... + 112.45378285, 114.53863928, 111.39181472, 112.67024271, ... + 113.65387448, 113.57576769, 110.22576589, 110.45091803, ... + 109.97966067, 112.84785553, 109.3836049 , 111.22290862, ... + 111.17028727, 110.71088554, 111.72973603, 112.39623635, ... + 113.18160954, 112.00016346, 109.66542778, 111.98705097, ... + 112.35983901, 110.21703075, 112.14565939, 111.2029942 , ... + 110.74296 , 112.56607054, 112.58015318, 111.93031032, ... + 111.59774288, 112.30723266, 112.79543793, 111.08288891, ... + 113.85269603, 111.77349741, 113.58639434, 111.28694353, ... + 112.1993445 , 111.72215918, 111.93271101, 112.17593036, ... + 110.82246602, 113.61806907, 114.13693144, 114.27245731, ... + 114.24223538, 112.61704196]; + +% ---------------------------- +% Load simulation data for fitting +% ---------------------------- +% baseDir = 'D:\Results - Numerics\Data_Full3D\PhaseDiagram\TFRadii\LowN'; +baseDir = 'C:\Users\Karthik\Documents\GitRepositories\Calculations\Estimations\ThomasFermiRadius'; +refData = load(fullfile(baseDir, 'TFFermi_Theta0.mat')); + +% ---------------------------- +% Find simulation values at the estimated atom number +% ---------------------------- +[~, idxClosest] = min(abs(refData.NUM_ATOMS_LIST - EstimatedAtomNumber)); + +if ndims(refData.TF_Radii) > 2 + TF_Radii = squeeze(refData.TF_Radii); + TF_X_target = TF_Radii(idxClosest, 1); + TF_Y_target = TF_Radii(idxClosest, 2); +else + TF_X_target = refData.TF_Radii(idxClosest, 1); + TF_Y_target = refData.TF_Radii(idxClosest, 2); +end + +fprintf('Target radii from simulation at N = %.0f:\n', refData.NUM_ATOMS_LIST(idxClosest)); +fprintf('TF_X_target = %.2f µm\n', TF_X_target); +fprintf('TF_Y_target = %.2f µm\n', TF_Y_target); + +% ---------------------------- +% Find optimal magnification +% ---------------------------- +errorFunc = @(mag) ... + (mean(TF_Radii_X_pixels)*PixelSize/mag - TF_X_target)^2 + ... + (mean(TF_Radii_Y_pixels)*PixelSize/mag - TF_Y_target)^2; + +Magnification = fminbnd(errorFunc, 10, 50); +fprintf('Best-fit magnification: %.4f\n', Magnification); + +% ---------------------------- +% Convert to real space and get stats +% ---------------------------- +TF_Radii_X_Real = TF_Radii_X_pixels * PixelSize / Magnification; +TF_Radii_Y_Real = TF_Radii_Y_pixels * PixelSize / Magnification; + +Avg_X = mean(TF_Radii_X_Real); +Avg_Y = mean(TF_Radii_Y_Real); +Std_X = std(TF_Radii_X_Real); +Std_Y = std(TF_Radii_Y_Real); + +fprintf('TF Radius X = %.2f ± %.2f µm\n', Avg_X, Std_X); +fprintf('TF Radius Y = %.2f ± %.2f µm\n', Avg_Y, Std_Y); + +% ---------------------------- +% Plotting +% ---------------------------- +fileList = {'TFFermi_Theta0.mat'}; +thetaLabels = {'\theta = 0^\circ', '\theta = 20^\circ', '\theta = 40^\circ'}; + +fig = figure(1); clf; +set(gcf,'Position', [100, 100, 1200, 500]) +t = tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); +colors = lines(length(fileList)); +legendEntries = cell(1, length(fileList)); + +for j = 1:length(fileList) + data = load(fullfile(baseDir, fileList{j})); + + aS = data.SCATTERING_LENGTH_RANGE; + NUM_ATOMS_LIST = data.NUM_ATOMS_LIST; + if ndims(data.TF_Radii) > 2 + TF_Radii = squeeze(data.TF_Radii); + else + TF_Radii = data.TF_Radii; + end + + legendEntries{j} = sprintf('%s, a_s = %.2f a_0', thetaLabels{j}, aS); + + % Rx + nexttile(1); + plot(NUM_ATOMS_LIST, TF_Radii(:,1), '-', ... + 'Color', colors(j,:), 'LineWidth', 1.5, ... + 'DisplayName', legendEntries{j}); hold on; + + % Ry + nexttile(2); + plot(NUM_ATOMS_LIST, TF_Radii(:,2), '-', ... + 'Color', colors(j,:), 'LineWidth', 1.5, ... + 'DisplayName', legendEntries{j}); hold on; +end + +% ---------------------------- +% Add experimental point w/ error bars and annotation +% ---------------------------- +% TF Radius X +nexttile(1); +errorbar(EstimatedAtomNumber, Avg_X, Std_X, ... + 'd', 'MarkerSize', 8, 'MarkerFaceColor', [0.2 0.2 0.8], ... + 'MarkerEdgeColor', 'k', 'Color', 'k', 'LineWidth', 1.2, ... + 'DisplayName', '\theta = 0^\circ, Experimental Value'); + +% TF Radius Y +nexttile(2); +errorbar(EstimatedAtomNumber, Avg_Y, Std_Y, ... + 'd', 'MarkerSize', 8, 'MarkerFaceColor', [0.2 0.2 0.8], ... + 'MarkerEdgeColor', 'k', 'Color', 'k', 'LineWidth', 1.2, ... + 'DisplayName', '\theta = 0^\circ, Experimental Value'); + +% ---------------------------- +% Finalize +% ---------------------------- +nexttile(1); +xlabel('Number of Atoms', 'FontSize', 16); +ylabel('TF Radius - X ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16); +legend('FontSize', 12, 'Interpreter', 'tex', 'Location', 'bestoutside'); +axis square; grid on; + +nexttile(2); +xlabel('Number of Atoms', 'FontSize', 16); +ylabel('TF Radius - Y ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16); +legend('FontSize', 12, 'Interpreter', 'tex', 'Location', 'bestoutside'); +axis square; grid on; + +sgtitle('[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz', ... + 'Interpreter', 'tex', 'FontSize', 18); diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m index 75ff163..c41a1aa 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -573,9 +573,14 @@ JobNumber = 0; Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) %% -SaveDirectory = 'D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_712500'; +SaveDirectory = 'D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_1733333'; JobNumber = 0; Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) +%% +SaveDirectory = 'D:/Results - Numerics/Data_Full3D/PhaseTransition/STD/aS_9.562000e+01_theta_025_phi_000_N_500000'; +JobNumber = 0; +Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) + %% Identify and count droplets Radius = 2; % The radius within which peaks will be considered duplicates @@ -591,7 +596,7 @@ 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); +AveragePCD = Scripts.extractAveragePeakColumnDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag); %% Extract average unit cell density - Droplets Radius = 2; % The radius within which peaks will be considered duplicates @@ -604,7 +609,7 @@ UCD = Scripts.extractAverageUnitCellDensity(SaveDirectory, J %% 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'; +SaveDirectory = 'D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_1325000'; JobNumber = 0; SuppressPlotFlag = false; UCD = Scripts.extractAverageUnitCellDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag); @@ -766,10 +771,10 @@ 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 +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"; @@ -777,7 +782,7 @@ 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)); +UCD_values = zeros(length(SCATTERING_LENGTH_RANGE), length(NUM_ATOMS_LIST)); % Prepare figure figure(1); @@ -824,6 +829,32 @@ legend(arrayfun(@(aS) sprintf('a_s = %.2f a_0', aS), SCATTERING_LENGTH_RANGE, .. set(gca, 'FontSize', 14); grid on; +% Physical constants +PlanckConstant = 6.62607015E-34; +PlanckConstantReduced = 6.62607015E-34/(2*pi); +AtomicMassUnit = 1.660539066E-27; +BohrMagneton = 9.274009994E-24; + +% Dy specific constants +Dy164Mass = 163.929174751*AtomicMassUnit; +Dy164IsotopicAbundance = 0.2826; +DyMagneticMoment = 9.93*BohrMagneton; + +add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length +nadd2s = 0.01:0.01:0.25; +ppmum = nadd2s.*(1E12*add^2)^-1; +% +figure(2); +clf; +set(gcf,'Position', [100, 100, 850, 700]); +hold on +plot(nadd2s, ppmum, 'o-', 'LineWidth', 1.5) +xlabel('na_{dd}^2', 'Interpreter', 'tex', 'FontSize', 16); +ylabel('Unit Cell Density (UCD) [$\mu m^{-2}$]', 'Interpreter', 'latex', 'FontSize', 16); +title("[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 0, 0, 72.4 ] Hz; \theta = 0^\circ", 'Interpreter', 'tex', 'FontSize', 18); +set(gca, 'FontSize', 14); +grid on; + %% Plot TF radii of unmodulated states % Parameters JobNumber = 0;