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

Compare commits

base: karthik:ea14ccddf4c466758c0e011ce304929fec3703fd
Branches Tags
karthik:master
...
compare: karthik:3115554b27d36a5fa5586732b39ec0cb8ddc92d3
Branches Tags
karthik:master

6 Commits
ea14ccddf4 ... 3115554b27

Author SHA1 Message Date
Karthik Chandrashekara
3115554b27 MAJOR restructuring of Data-Analyzer: code re-written to be highly modular. 2025-08-18 10:24:55 +02:00
Karthik Chandrashekara
be84ba8bd5 Revert "MAJOR restructuring: Code rewritten to be highly modular." 
This reverts commit f61ab4deff.
 2025-08-18 00:27:27 +02:00
Karthik Chandrashekara
f61ab4deff MAJOR restructuring: Code rewritten to be highly modular. 2025-08-18 00:14:18 +02:00
Karthik Chandrashekara
ea9b40793f Revert "MAJOR restructuring: code rewritten to be highly modular." 
This reverts commit 9011919780.
 2025-08-18 00:10:17 +02:00
Karthik Chandrashekara
9011919780 MAJOR restructuring: code rewritten to be highly modular. 2025-08-17 23:40:06 +02:00
Karthik Chandrashekara
3519f0de82 MAJOR restructuring: code now rewritten to be highly modular. 2025-08-17 23:37:38 +02:00
52 changed files with 2599 additions and 5890 deletions
Show Stats Expand all files Collapse all files

111
Data-Analyzer/+Analyzer/conductPCA.m Normal file
View File

@ -0,0 +1,111 @@
function conductPCA(od_imgs, scan_reference_values, scan_parameter_values, doPlot, doSave, saveDir)
%% Performs PCA on optical density images, visualizes and optionally saves results.
%
% Inputs:
% od_imgs - cell array of OD images
% scan_reference_values - array of unique control parameter values
% scan_parameter_values - array mapping each image to a control parameter
% doPlot - logical, true to plot figures
% doSave - logical, true to save figures
% saveDir - directory to save figures if doSave is true
%
% Requires:
% +Calculator/computeCumulants.m
if nargin < 4, doPlot = true; end
if nargin < 5, doSave = false; end
if nargin < 6, saveDir = pwd; end
%% PCA computation
allImgs3D = cat(3, od_imgs{:});
[Nx, Ny] = size(allImgs3D(:,:,1));
Xall = reshape(allImgs3D, [], numel(od_imgs))';
[coeff, score, ~, ~, explained] = pca(Xall);
figCount = 1;
%% --- Figure 1: PC1 Image ---
if doPlot
pc1_vector = coeff(:,1);
pc1_image = reshape(pc1_vector, Nx, Ny);
figure(figCount); clf; set(gcf, 'Color', 'w', 'Position', [100 100 950 750]);
imagesc(pc1_image); axis image off; colormap(Colormaps.coolwarm()); colorbar;
title(sprintf('First Principal Component (PC1) Image - Explains %.2f%% Variance', explained(1)));
if doSave, saveas(gcf, fullfile(saveDir, 'PC1_Image.png')); end
figCount = figCount + 1;
end
%% --- Figure 2: PC1 Scores Scatter ---
if doPlot
numGroups = numel(scan_reference_values);
colors = lines(numGroups);
figure(figCount); clf; set(gcf, 'Color', 'w', 'Position', [100 100 950 750]); hold on;
for g = 1:numGroups
idx = scan_parameter_values == scan_reference_values(g);
scatter(repmat(scan_reference_values(g), sum(idx),1), score(idx,1), 36, colors(g,:), 'filled');
end
xlabel('Control Parameter'); ylabel('PC1 Score');
title('Evolution of PC1 Scores'); grid on;
if doSave, saveas(gcf, fullfile(saveDir, 'PC1_Scatter.png')); end
figCount = figCount + 1;
end
%% --- Figure 3: PC1 Distributions ---
if doPlot
figure(figCount); clf; set(gcf, 'Color', 'w', 'Position', [100 100 950 750]);
hold on;
for g = 1:numGroups
idx = scan_parameter_values == scan_reference_values(g);
data = score(idx,1);
histogram(data, 'Normalization', 'pdf', 'FaceColor', colors(g,:), 'FaceAlpha', 0.3);
[f, xi] = ksdensity(data);
plot(xi, f, 'Color', colors(g,:), 'LineWidth', 2);
end
xlabel('PC1 Score'); ylabel('Probability Density');
title('PC1 Score Distributions by Group');
legend(arrayfun(@num2str, scan_reference_values, 'UniformOutput', false), 'Location', 'Best');
grid on;
if doSave, saveas(gcf, fullfile(saveDir, 'PC1_Distributions.png')); end
figCount = figCount + 1;
end
%% --- Figure 4: Boxplot of PC1 Scores ---
if doPlot
figure(figCount); clf; set(gcf, 'Color', 'w', 'Position', [100 100 950 750]);
boxplot(score(:,1), scan_parameter_values);
xlabel('Control Parameter'); ylabel('PC1 Score');
title('PC1 Score Boxplots by Group'); grid on;
if doSave, saveas(gcf, fullfile(saveDir, 'PC1_Boxplot.png')); end
figCount = figCount + 1;
end
%% --- Figure 5: Mean ± SEM of PC1 Scores ---
if doPlot
meanScores = arrayfun(@(g) mean(score(scan_parameter_values == g,1)), scan_reference_values);
semScores = arrayfun(@(g) std(score(scan_parameter_values == g,1))/sqrt(sum(scan_parameter_values == g)), scan_reference_values);
figure(figCount); clf; set(gcf, 'Color', 'w', 'Position', [100 100 950 750]);
errorbar(scan_reference_values, meanScores, semScores, '-o', 'LineWidth', 2);
xlabel('Control Parameter'); ylabel('Mean PC1 Score ± SEM');
title('Mean ± SEM of PC1 Scores'); grid on;
if doSave, saveas(gcf, fullfile(saveDir, 'PC1_MeanSEM.png')); end
figCount = figCount + 1;
end
%% --- Figure 6: Binder Cumulant ---
if doPlot
binderVals = arrayfun(@(g) ...
Calculator.computeCumulants(score(scan_parameter_values == g,1)), ...
scan_reference_values);
figure(figCount); clf; set(gcf, 'Color', 'w', 'Position', [100 100 950 750]);
plot(scan_reference_values, binderVals, '-o', 'LineWidth', 2);
xlabel('Control Parameter'); ylabel('Binder Cumulant (PC1)');
title('Binder Cumulant of PC1 Scores'); grid on;
if doSave, saveas(gcf, fullfile(saveDir, 'PC1_BinderCumulant.png')); end
end
%% --- ANOVA Test ---
p = anova1(score(:,1), arrayfun(@num2str, scan_parameter_values, 'UniformOutput', false), 'off');
fprintf('ANOVA p-value for PC1 score differences between groups: %.4e\n', p);
end

296
Data-Analyzer/+Analyzer/conductSpectralAnalysis.m Normal file
View File

@ -0,0 +1,296 @@
function results = conductSpectralAnalysis(od_imgs, scan_parameter_values, options)
%% Performs Fourier analysis on a set of optical density (OD) images.
% Computes radial and angular spectral distributions, optionally plots
% results, saves figures, and can render a video of the analysis.
%
% Inputs:
% od_imgs - cell array of OD images
% scan_parameter_values - array of scan parameter values corresponding to each image
% OPTIONS -
% saveDirectory - directory to save files
% savefileName - base filename for saved figures/video
% skipMovieRender - skip creating the video of analysis
% skipSaveFigures - skip saving plots
% skipSaveOD - skip saving OD images as .mat
% skipPreprocessing - skip preprocessing of images before FFT
% skipMasking - skip masking of OD images
% skipIntensityThresholding- skip thresholding of intensity
% skipBinarization - skip binarization of OD images
% skipNormalization - skip normalization when plotting angular spectrum
% skipLivePlot = skip live plotting of figures
% pixel_size - physical pixel size of camera sensor (m)
% magnification - imaging magnification
% zoom_size - number of pixels to crop around FFT center
% k_min, k_max - min/max wavenumber for spectral contrast
% N_angular_bins - number of angular bins for S(θ)
% Angular_Threshold - threshold parameter for angular spectrum
% Angular_Sigma - Gaussian smoothing width for angular spectrum
% theta_min, theta_max - angular range for radial spectrum integration
% N_radial_bins - number of radial bins for S(k)
% Radial_WindowSize - window size for smoothing radial spectrum
% scan_parameter - string, type of scan parameter (used in plot text)
% font - font name for plots
%% Unpack struct arguments
pixel_size = options.pixel_size;
magnification = options.magnification;
zoom_size = options.zoom_size;
k_min = options.k_min;
k_max = options.k_max;
N_angular_bins = options.N_angular_bins;
Angular_Threshold = options.Angular_Threshold;
Angular_Sigma = options.Angular_Sigma;
theta_min = options.theta_min;
theta_max = options.theta_max;
N_radial_bins = options.N_radial_bins;
Radial_WindowSize = options.Radial_WindowSize;
skipNormalization = options.skipNormalization;
skipPreprocessing = options.skipPreprocessing;
skipMasking = options.skipMasking;
skipIntensityThresholding = options.skipIntensityThresholding;
skipBinarization = options.skipBinarization;
skipLivePlot = options.skipLivePlot;
skipMovieRender = options.skipMovieRender;
skipSaveFigures = options.skipSaveFigures;
skipSaveOD = options.skipSaveOD;
savefileName = options.savefileName;
saveDirectory = options.saveDirectory;
scan_parameter = options.scan_parameter;
font = options.font;
%% ===== Initialization =====
N_shots = length(od_imgs); % total number of images
fft_imgs = cell(1, N_shots); % FFT of each image
angular_spectral_distribution = cell(1, N_shots); % S(θ) angular spectrum
radial_spectral_contrast = zeros(1, N_shots); % radial contrast metric
angular_spectral_weight = zeros(1, N_shots); % integrated angular weight
S_theta_all = cell(1, N_shots);
S_k_all = cell(1, N_shots);
S_k_smoothed_all = cell(1, N_shots);
S_theta_norm_all = cell(1, N_shots);
% Optional save directory override
if ~isempty(saveDirectory)
savefileName = fullfile(saveDirectory, savefileName);
end
% Prepare video if enabled
if ~skipMovieRender
videoFile = VideoWriter([savefileName '.mp4'], 'MPEG-4');
videoFile.Quality = 100;
videoFile.FrameRate = 2;
open(videoFile);
end
% Prepare folder to save figures
if ~skipSaveFigures
saveFolder = [savefileName '_SavedFigures'];
if ~exist(saveFolder, 'dir')
mkdir(saveFolder);
end
end
% Initialize lists for power spectra and radial spectra
PS_all = cell(1, N_shots); % 2D FFT power spectrum |F(kx,ky)|^2
%% ===== Main loop over images =====
for k = 1:N_shots
IMG = od_imgs{k};
% Skip FFT if image is empty or has low intensity
if ~(max(IMG(:)) > 1)
IMGFFT = NaN(size(IMG));
else
% Compute FFT with optional preprocessing
[IMGFFT, ~] = Calculator.computeFourierTransform(IMG, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization);
end
% Image size
[Ny, Nx] = size(IMG);
% Real-space pixel size (meters)
dx = pixel_size / magnification;
dy = dx; % assume square pixels
% Real-space axes in µm
x = ((1:Nx) - ceil(Nx/2)) * dx * 1E6;
y = ((1:Ny) - ceil(Ny/2)) * dy * 1E6;
% Reciprocal space increments
dvx = 1 / (Nx * dx);
dvy = 1 / (Ny * dy);
% Frequency axes
vx = (-floor(Nx/2):ceil(Nx/2)-1) * dvx;
vy = (-floor(Ny/2):ceil(Ny/2)-1) * dvy;
% Wavenumber axes (µm⁻¹)
kx_full = 2 * pi * vx * 1E-6;
ky_full = 2 * pi * vy * 1E-6;
% Crop FFT image around center
mid_x = floor(Nx/2);
mid_y = floor(Ny/2);
fft_imgs{k} = IMGFFT(mid_y-zoom_size:mid_y+zoom_size, mid_x-zoom_size:mid_x+zoom_size);
% Crop wavenumber axes to match cropped FFT
kx = kx_full(mid_x - zoom_size : mid_x + zoom_size);
ky = ky_full(mid_y - zoom_size : mid_y + zoom_size);
%% ===== Spectral analysis =====
% Angular spectrum
[theta_vals, S_theta] = Calculator.computeAngularSpectralDistribution(fft_imgs{k}, kx, ky, k_min, k_max, N_angular_bins, Angular_Threshold, Angular_Sigma, []);
% Radial spectrum
[k_rho_vals, S_k] = Calculator.computeRadialSpectralDistribution(fft_imgs{k}, kx, ky, theta_min, theta_max, N_radial_bins);
% Smooth radial spectrum
S_k_smoothed = movmean(S_k, Radial_WindowSize);
% Store results
angular_spectral_distribution{k} = S_theta;
radial_spectral_contrast(k) = Calculator.computeRadialSpectralContrast(k_rho_vals, S_k_smoothed, k_min, k_max);
% Normalize angular spectrum and compute weight
S_theta_norm = S_theta / max(S_theta);
angular_spectral_weight(k) = trapz(theta_vals, S_theta_norm);
% Store results
S_theta_all{k} = S_theta;
S_k_all{k} = S_k;
S_k_smoothed_all{k} = S_k_smoothed;
S_theta_norm_all{k} = S_theta_norm;
PS_all{k} = abs(fft_imgs{k}).^2;
%% ===== Plotting =====
if ~skipLivePlot
figure(1); clf
set(gcf,'Position',[500 100 1000 800])
tiledlayout(2, 2, 'TileSpacing', 'compact', 'Padding', 'compact');
% OD image
ax1 = nexttile;
imagesc(x, y, IMG)
hold on;
Helper.drawODOverlays(x(1), y(1), x(end), y(end));
Helper.drawODOverlays(x(end), y(1), x(1), y(end));
hold off;
axis equal tight;
set(gca, 'FontSize', 14, 'YDir', 'normal')
colormap(ax1, Colormaps.inferno());
hcb = colorbar;
ylabel(hcb, 'Optical Density', 'Rotation', -90, 'FontSize', 14, 'FontName', font);
xlabel('x (\mum)', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font);
ylabel('y (\mum)', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font);
title('OD Image', 'FontSize', 16, 'FontWeight', 'bold', 'Interpreter', 'tex', 'FontName', font);
% Annotate scan parameter
if strcmp(scan_parameter, 'ps_rot_mag_fin_pol_angle')
text(0.975, 0.975, sprintf('%.1f^\\circ', scan_parameter_values(k)), ...
'Color', 'white', 'FontWeight', 'bold', 'FontSize', 14, ...
'Interpreter', 'tex', 'Units', 'normalized', ...
'HorizontalAlignment', 'right', 'VerticalAlignment', 'top');
else
text(0.975, 0.975, sprintf('%.2f G', scan_parameter_values(k)), ...
'Color', 'white', 'FontWeight', 'bold', 'FontSize', 14, ...
'Interpreter', 'tex', 'Units', 'normalized', ...
'HorizontalAlignment', 'right', 'VerticalAlignment', 'top');
end
% FFT power spectrum
ax2 = nexttile;
imagesc(kx, ky, log(1 + PS_all{k}));
hold on;
Helper.drawPSOverlays(kx, ky, k_min, k_max)
% Restrict axes strictly to image limits
xlim([min(kx), max(kx)]);
ylim([min(ky), max(ky)]);
axis image; % preserves aspect ratio
set(gca, 'FontSize', 14, 'YDir', 'normal')
xlabel('k_x [\mum^{-1}]', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font);
ylabel('k_y [\mum^{-1}]', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font);
title('Power Spectrum - S(k_x,k_y)', 'Interpreter', 'tex', ...
'FontSize', 16, 'FontWeight', 'bold', 'FontName', font);
colorbar;
colormap(ax2, Colormaps.coolwarm());
% Radial distribution
nexttile;
plot(k_rho_vals, S_k_smoothed, 'LineWidth', 2);
set(gca, 'FontSize', 14, 'YScale', 'log', 'XLim', [min(k_rho_vals), max(k_rho_vals)]);
xlabel('k_\rho [\mum^{-1}]', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font);
ylabel('Magnitude (a.u.)', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font);
title('Radial Spectral Distribution - S(k_\rho)', 'Interpreter', 'tex', ...
'FontSize', 16, 'FontWeight', 'bold', 'FontName', font);
grid on;
% Angular distribution
nexttile;
if ~skipNormalization
plot(theta_vals/pi, S_theta_norm, 'LineWidth', 2);
set(gca, 'FontSize', 14, 'YLim', [0, 1]);
else
plot(theta_vals/pi, S_theta, 'LineWidth', 2);
set(gca, 'FontSize', 14, 'YScale', 'log', 'YLim', [1E4, 1E7]);
end
xlabel('\theta/\pi [rad]', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font);
ylabel('Magnitude (a.u.)', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font);
title('Angular Spectral Distribution - S(\theta)', 'Interpreter', 'tex', ...
'FontSize', 16, 'FontWeight', 'bold', 'FontName', font);
grid on;
ax = gca;
ax.MinorGridLineStyle = ':';
ax.MinorGridColor = [0.7 0.7 0.7];
ax.MinorGridAlpha = 0.5;
ax.XMinorGrid = 'on';
ax.YMinorGrid = 'on';
end
%% ===== Save outputs =====
if ~skipMovieRender
frame = getframe(gcf);
writeVideo(videoFile, frame);
end
if ~skipSaveFigures
fileNamePNG = fullfile(saveFolder, sprintf('fft_analysis_img_%03d.png', k));
print(gcf, fileNamePNG, '-dpng', '-r100');
end
if ~skipSaveOD
odDataStruct = struct();
odDataStruct.IMG = IMG;
odDataStruct.x = x;
odDataStruct.y = y;
odDataStruct.scan_parameter_value = scan_parameter_values(k);
save(fullfile(saveFolder, sprintf('od_image_%03d.mat', k)), '-struct', 'odDataStruct');
end
if skipMovieRender && skipSaveFigures
pause(0.5);
end
end
% Package results into struct
results = struct();
results.kx = kx;
results.ky = ky;
results.PS_all = PS_all;
results.theta_vals = theta_vals;
results.S_theta_all = S_theta_all;
results.k_rho_vals = k_rho_vals;
results.S_k_all = S_k_all;
results.angular_spectral_distribution = angular_spectral_distribution;
results.S_k_smoothed_all = S_k_smoothed_all;
results.radial_spectral_contrast = radial_spectral_contrast;
results.S_theta_norm_all = S_theta_norm_all;
results.angular_spectral_weight = angular_spectral_weight;
if ~skipMovieRender
close(videoFile);
end
end

49
Data-Analyzer/+Analyzer/extractAutocorrelation.m Normal file
View File

@ -0,0 +1,49 @@
function results = extractAutocorrelation(theta_values, angular_spectral_distribution, scan_parameter_values, N_shots, N_angular_bins)
%% Extract g² (autocorrelation) from experimental images
% Computes angular autocorrelation g² for a set of experimental images.
% Uses conductSpectralAnalysis to compute S(θ) and θ-values, then groups
% images by scan parameter and computes normalized autocorrelations.
% ===== Convert spectral distributions to matrix =====
delta_nkr_all = zeros(N_shots, N_angular_bins);
for k = 1:N_shots
delta_nkr_all(k, :) = angular_spectral_distribution{k};
end
% ===== Group images by scan parameter values =====
[unique_scan_parameter_values, ~, idx] = unique(scan_parameter_values);
N_params = length(unique_scan_parameter_values);
% ===== Preallocate output arrays =====
g2_all = zeros(N_params, N_angular_bins);
g2_error_all = zeros(N_params, N_angular_bins);
% ===== Compute g²(θ) for each scan parameter group =====
for i = 1:N_params
group_idx = find(idx == i);
group_data = delta_nkr_all(group_idx, :);
for dtheta = 0:N_angular_bins-1
temp = zeros(length(group_idx), 1);
for j = 1:length(group_idx)
profile = group_data(j, :);
profile_shifted = circshift(profile, -dtheta, 2);
num = mean(profile .* profile_shifted);
denom = mean(profile.^2);
temp(j) = num / denom;
end
g2_all(i, dtheta+1) = mean(temp, 'omitnan');
g2_error_all(i, dtheta+1) = std(temp, 'omitnan') / sqrt(length(group_idx));
end
end
% ===== Package results =====
results = struct();
results.g2_all = g2_all;
results.g2_error_all = g2_error_all;
results.theta_values = theta_values;
end

107
Data-Analyzer/+Analyzer/extractCustomCorrelation.m Normal file
View File

@ -0,0 +1,107 @@
function results = extractCustomCorrelation(angular_spectral_distribution, scan_parameter_values, N_shots, N_angular_bins)
%% Extracts correlation of a single (highest) peak with possible secondary peak (50-70°)
%
% Inputs:
% od_imgs - Cell array of images
% scan_parameter_values - Vector of scan parameters corresponding to images
% pixel_size - Camera pixel size in meters
% magnification - Imaging magnification
% zoom_size - Half-size of FFT crop around center
% r_min, r_max - Radial bounds for angular spectral distribution
% N_angular_bins - Number of angular bins
% Angular_Threshold, Angular_Sigma - Parameters for angular weighting
% skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization - flags for FFT preprocessing
%
% Output:
% results - Struct containing g2 correlation and cumulant statistics per scan parameter
% ===== Convert spectral distributions to matrix (N_shots x N_angular_bins) =====
delta_nkr_all = zeros(N_shots, N_angular_bins);
for k = 1:N_shots
delta_nkr_all(k, :) = angular_spectral_distribution{k};
end
% ===== Group images by scan parameter values =====
[unique_scan_parameter_values, ~, idx] = unique(scan_parameter_values);
N_params = length(unique_scan_parameter_values);
% ===== Angular settings =====
angle_range = 180;
angle_per_bin = angle_range / N_angular_bins;
max_peak_bin = round(180 / angle_per_bin);
window_size = 10;
angle_threshold = 100;
% ===== Preallocate result arrays =====
mean_max_g2_values = zeros(1, N_params);
skew_max_g2_angle_values = zeros(1, N_params);
var_max_g2_values = zeros(1, N_params);
fourth_order_cumulant_max_g2_angle_values = zeros(1, N_params);
max_g2_all_per_scan_parameter_value = cell(1, N_params);
std_error_g2_values = zeros(1, N_params);
% ===== Compute correlations and cumulants per group =====
for i = 1:N_params
group_idx = find(idx == i);
group_data = delta_nkr_all(group_idx, :);
N_reps = size(group_data, 1);
g2_values = zeros(1, N_reps);
for j = 1:N_reps
profile = group_data(j, :);
% Find highest peak in 0–180° range
restricted_profile = profile(1:max_peak_bin);
[~, peak_idx_rel] = max(restricted_profile);
peak_idx = peak_idx_rel;
peak_angle = (peak_idx - 1) * angle_per_bin;
% Determine offsets for secondary peak correlation
if peak_angle < angle_threshold
offsets = round(50 / angle_per_bin) : round(70 / angle_per_bin);
else
offsets = -round(70 / angle_per_bin) : -round(50 / angle_per_bin);
end
ref_window = mod((peak_idx - window_size):(peak_idx + window_size) - 1, N_angular_bins) + 1;
ref = profile(ref_window);
correlations = zeros(size(offsets));
for k_off = 1:length(offsets)
shifted_idx = mod(peak_idx + offsets(k_off) - 1, N_angular_bins) + 1;
sec_window = mod((shifted_idx - window_size):(shifted_idx + window_size) - 1, N_angular_bins) + 1;
sec = profile(sec_window);
correlations(k_off) = mean(ref .* sec) / mean(ref.^2);
end
[max_corr, ~] = max(correlations);
g2_values(j) = max_corr;
end
% Store raw values
max_g2_all_per_scan_parameter_value{i} = g2_values;
% Compute cumulants
kappa = Calculator.computeCumulants(g2_values(:), 4);
mean_max_g2_values(i) = kappa(1);
var_max_g2_values(i) = kappa(2);
skew_max_g2_angle_values(i) = kappa(3);
fourth_order_cumulant_max_g2_angle_values(i) = kappa(4);
N_eff = sum(~isnan(g2_values));
std_error_g2_values(i) = sqrt(kappa(2)) / sqrt(N_eff);
end
% ===== Package results into struct =====
results = struct();
results.mean_max_g2 = mean_max_g2_values;
results.var_max_g2 = var_max_g2_values;
results.skew_max_g2_angle = skew_max_g2_angle_values;
results.fourth_order_cumulant_max_g2 = fourth_order_cumulant_max_g2_angle_values;
results.std_error_g2 = std_error_g2_values;
results.max_g2_all_per_scan_parameter_value = max_g2_all_per_scan_parameter_value;
end

78
Data-Analyzer/+Analyzer/performAnalysis.m Normal file
View File

@ -0,0 +1,78 @@
function results = performAnalysis(options)
arguments
options.scan_parameter (1,:) char
options.scan_reference_values (1,:) double
options.cam (1,1) double
options.angle (1,1) double
options.center (1,2) double
options.span (1,2) double
options.fraction (1,2) double
options.ImagingMode (1,:) char
options.PulseDuration (1,1) double
options.removeFringes (1,1) logical
options.skipUnshuffling (1,1) logical
options.pixel_size (1,1) double
options.magnification (1,1) double
options.zoom_size (1,1) double
options.N_angular_bins (1,1) double
options.Angular_Threshold (1,1) double
options.Angular_Sigma (1,1) double
options.Angular_WindowSize (1,1) double
options.theta_min (1,1) double
options.theta_max (1,1) double
options.N_radial_bins (1,1) double
options.Radial_Sigma (1,1) double
options.Radial_WindowSize (1,1) double
options.k_min (1,1) double
options.k_max (1,1) double
options.skipPreprocessing (1,1) logical
options.skipMasking (1,1) logical
options.skipIntensityThresholding (1,1) logical
options.skipBinarization (1,1) logical
options.skipNormalization (1,1) logical
options.skipLivePlot (1,1) logical
options.skipMovieRender (1,1) logical
options.skipSaveFigures (1,1) logical
options.skipSaveOD (1,1) logical
options.showProgressBar (1,1) logical
options.savefileName (1,:) char
options.folderPath (1,:) char
options.baseDataFolder (1,:) char
options.saveDirectory (1,:) char
options.titleString (1,:) char
options.font (1,:) char
end
% Collect OD images
[od_imgs, scan_parameter_values, ~] = Helper.collectODImages(options);
% Conduct spectral analysis
fprintf('\nInitiating spectral analysis...\n');
spectral_analysis_results = Analyzer.conductSpectralAnalysis(od_imgs, scan_parameter_values, options);
N_shots = length(od_imgs);
% Extract angular correlations
full_g2_results = Analyzer.extractAutocorrelation(...
spectral_analysis_results.theta_vals, ...
spectral_analysis_results.angular_spectral_distribution, ...
scan_parameter_values, N_shots, options.N_angular_bins);
custom_g_results = Analyzer.extractCustomCorrelation(...
spectral_analysis_results.angular_spectral_distribution, ...
scan_parameter_values, N_shots, options.N_angular_bins);
fprintf('\nSpectral analysis complete!\n');
% PCA
% Lattice Reconstruction
% Package results into struct
results = struct();
results.spectral_analysis_results = spectral_analysis_results;
results.full_g2_results = full_g2_results;
results.custom_g_results = custom_g_results;
end

96
Data-Analyzer/+Analyzer/runInteractiveODImageViewer.m Normal file
View File

@ -0,0 +1,96 @@
function runInteractiveODImageViewer(od_imgs, scan_parameter_values, file_list, options)
%% Interactive OD Image Viewer
% od_imgs : cell array of 2D OD images
% scan_parameter_values : array of corresponding scan parameter values
% file_list : cell array of corresponding filenames
% options : struct with fields
% .pixel_size, .magnification, .center, .font, .zoom_size, .scan_parameter
% Figure
hFig = figure('Name', 'OD Image Viewer', 'NumberTitle', 'off', 'Position', [50 50 1000 800]);
% Get image size
[Ny, Nx] = size(od_imgs{1});
% Pixel size and axes in μm
dx = options.pixel_size / options.magnification;
dy = dx; % square pixels
x = ((1:Nx) - (Nx+1)/2) * dx * 1e6;
y = ((1:Ny) - (Ny+1)/2) * dy * 1e6;
% Display first image
hAx = axes('Parent', hFig);
hImg = imagesc(hAx, x, y, od_imgs{1});
axis(hAx, 'equal', 'tight')
colormap(hAx, Colormaps.inferno());
set(hAx, 'FontSize', 14, 'YDir', 'normal');
xlabel(hAx, 'x (\mum)', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', options.font);
ylabel(hAx, 'y (\mum)', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', options.font);
title(hAx, ['Measurement: ', options.titleString], 'FontSize', 16, ...
'FontWeight', 'bold', 'Interpreter', 'tex', 'FontName', options.font);
colorbarHandle = colorbar(hAx);
ylabel(colorbarHandle, 'Optical Density', 'Rotation', -90, 'FontSize', 14, 'FontName', options.font);
hold(hAx, 'on')
% Draw diagonal overlays once
Helper.drawODOverlays(x(1), y(1), x(end), y(end));
Helper.drawODOverlays(x(end), y(1), x(1), y(end));
hold(hAx, 'off')
txtHandle = text(hAx, 0.975, 0.975, '', ...
'Color', 'white', 'FontWeight', 'bold', ...
'FontSize', 24, 'Interpreter', 'tex', ...
'Units', 'normalized', ...
'HorizontalAlignment', 'right', ...
'VerticalAlignment', 'top');
% Slider
sliderHandle = uicontrol('Style', 'slider', ...
'Min', 1, 'Max', length(od_imgs), 'Value', 1, ...
'SliderStep', [1/(length(od_imgs)-1), 10/(length(od_imgs)-1)], ...
'Position', [150 5 700 20], ...
'Callback', @(src, ~) updateImage(round(src.Value)));
% Initialize
currentIdx = 1;
updateImage(currentIdx);
% Arrow key callback
set(hFig, 'KeyPressFcn', @(src, event) keyPressCallback(event));
%% --- Nested Functions ---
function updateImage(idx)
currentIdx = idx;
hImg.CData = od_imgs{idx};
% Extract only filename (without path)
[~, fname, ext] = fileparts(file_list{idx});
shortName = [fname, ext];
% Update figure title with shot + filename
if strcmp(options.scan_parameter, 'rot_mag_fin_pol_angle')
hFig.Name = sprintf('Shot %d | %s', idx, shortName);
txtHandle.String = sprintf('%.1f^\\circ', scan_parameter_values(idx));
else
hFig.Name = sprintf('Shot %d | %s', idx, shortName);
txtHandle.String = sprintf('%.2f G', scan_parameter_values(idx));
end
sliderHandle.Value = idx;
drawnow;
end
function keyPressCallback(event)
switch event.Key
case 'rightarrow'
if currentIdx < length(od_imgs)
updateImage(currentIdx + 1);
end
case 'leftarrow'
if currentIdx > 1
updateImage(currentIdx - 1);
end
end
end
end

46
Data-Analyzer/+Calculator/computeAngularSpectralDistribution.m Normal file
View File

@ -0,0 +1,46 @@
function [theta_vals, S_theta] = computeAngularSpectralDistribution(IMGFFT, kx, ky, k_min, k_max, num_bins, threshold, sigma, windowSize)
% Apply threshold to isolate strong peaks
IMGFFT(IMGFFT < threshold) = 0;
% Create wavenumber meshgrid
[KX, KY] = meshgrid(kx, ky);
Kmag = sqrt(KX.^2 + KY.^2); % radial wavenumber magnitude
Theta = atan2(KY, KX); % range [-pi, pi]
% Restrict to radial band in wavenumber space
radial_mask = (Kmag >= k_min) & (Kmag <= k_max);
% Initialize angular structure factor
S_theta = zeros(1, num_bins);
theta_vals = linspace(0, pi, num_bins); % only 0 to pi due to symmetry
% Loop over angular bins
for i = 1:num_bins
angle_start = (i - 1) * pi / num_bins;
angle_end = i * pi / num_bins;
angle_mask = (Theta >= angle_start) & (Theta < angle_end);
bin_mask = radial_mask & angle_mask;
fft_angle = IMGFFT .* bin_mask;
S_theta(i) = sum(sum(abs(fft_angle).^2));
end
% Optional smoothing
if exist('sigma', 'var') && ~isempty(sigma)
% Gaussian smoothing
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);
% Circular convolution
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);
elseif exist('windowSize', 'var') && ~isempty(windowSize)
% Moving average smoothing
pad = floor(windowSize / 2);
kernel = ones(1, windowSize) / windowSize;
S_theta = conv([S_theta(end - pad + 1:end), S_theta, S_theta(1:pad)], kernel, 'same');
S_theta = S_theta(pad + 1:end - pad);
end
end

0
Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/computeCumulants.m → Data-Analyzer/+Calculator/computeCumulants.m
View File

70
Data-Analyzer/+Calculator/computeFourierTransform.m Normal file
View File

@ -0,0 +1,70 @@
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

28
Data-Analyzer/+Calculator/computeRadialSpectralContrast.m Normal file
View File

@ -0,0 +1,28 @@
function contrast = computeRadialSpectralContrast(k_rho_vals, S_k_smoothed, k_min, k_max)
% Computes the ratio of the peak in S_k_smoothed within [k_min, k_max]
% to the value at (or near) k = 0.
% Ensure inputs are column vectors
k_rho_vals = k_rho_vals(:);
S_k_smoothed = S_k_smoothed(:);
% Step 1: Find index of k ≈ 0
[~, idx_k0] = min(abs(k_rho_vals)); % Closest to zero
S_k0 = S_k_smoothed(idx_k0);
% Step 2: Find indices in specified k-range
in_range = (k_rho_vals >= k_min) & (k_rho_vals <= k_max);
if ~any(in_range)
warning('No values found in the specified k-range. Returning NaN.');
contrast = NaN;
return;
end
% Step 3: Find peak value in the specified k-range
S_k_peak = max(S_k_smoothed(in_range));
% Step 4: Compute contrast
contrast = S_k_peak / S_k0;
end

33
Data-Analyzer/+Calculator/computeRadialSpectralDistribution.m Normal file
View File

@ -0,0 +1,33 @@
function [k_rho_vals, S_radial] = computeRadialSpectralDistribution(IMGFFT, kx, ky, thetamin, thetamax, num_bins)
% IMGFFT : 2D FFT image (fftshifted and cropped)
% kx, ky : 1D physical wavenumber axes [μm⁻¹] matching FFT size
% thetamin : Minimum angle (in radians)
% thetamax : Maximum angle (in radians)
% num_bins : Number of radial bins
[KX, KY] = meshgrid(kx, ky);
K_rho = sqrt(KX.^2 + KY.^2);
Theta = atan2(KY, KX);
if thetamin < thetamax
angle_mask = (Theta >= thetamin) & (Theta <= thetamax);
else
angle_mask = (Theta >= thetamin) | (Theta <= thetamax);
end
power_spectrum = abs(IMGFFT).^2;
r_min = min(K_rho(angle_mask));
r_max = max(K_rho(angle_mask));
r_edges = linspace(r_min, r_max, num_bins + 1);
k_rho_vals = 0.5 * (r_edges(1:end-1) + r_edges(2:end));
S_radial = zeros(1, num_bins);
for i = 1:num_bins
r_low = r_edges(i);
r_high = r_edges(i + 1);
radial_mask = (K_rho >= r_low) & (K_rho < r_high);
full_mask = radial_mask & angle_mask;
S_radial(i) = sum(power_spectrum(full_mask));
end
end

0
Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/+Colormaps/coolwarm.m → Data-Analyzer/+Colormaps/coolwarm.m
View File

0
Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/+Colormaps/fake_parula.m → Data-Analyzer/+Colormaps/fake_parula.m
View File

0
Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/+Colormaps/inferno.m → Data-Analyzer/+Colormaps/inferno.m
View File

0
Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/+Colormaps/magma.m → Data-Analyzer/+Colormaps/magma.m
View File

0
Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/+Colormaps/plasma.m → Data-Analyzer/+Colormaps/plasma.m
View File

0
Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/+Colormaps/viridis.m → Data-Analyzer/+Colormaps/viridis.m
View File

36
Data-Analyzer/+Helper/PhysicsConstants.m Normal file
View File

@ -0,0 +1,36 @@
classdef PhysicsConstants < handle
properties (Constant)
% CODATA
PlanckConstant=6.62607015E-34;
PlanckConstantReduced=6.62607015E-34/(2*pi);
FineStructureConstant=7.2973525698E-3;
ElectronMass=9.10938291E-31;
GravitationalConstant=6.67384E-11;
ProtonMass=1.672621777E-27;
AtomicMassUnit=1.660539066E-27;
BohrRadius=5.2917721067E-11;
BohrMagneton=9.274009994E-24;
BoltzmannConstant=1.38064852E-23;
StandardGravityAcceleration=9.80665;
SpeedOfLight=299792458;
StefanBoltzmannConstant=5.670373E-8;
ElectronCharge=1.602176634E-19;
VacuumPermeability=1.25663706212E-6;
DielectricConstant=8.8541878128E-12;
ElectronGyromagneticFactor=-2.00231930436153;
AvogadroConstant=6.02214076E23;
ZeroKelvin = 273.15;
GravitationalAcceleration = 9.80553;
% Dy specific constants
Dy164Mass = 163.929174751*1.660539066E-27;
Dy164IsotopicAbundance = 0.2826;
DyMagneticMoment = 9.93*9.274009994E-24;
end
methods
function pc = PhysicsConstants()
end
end
end

68
Data-Analyzer/+Helper/ProgressBar.m Normal file
View File

@ -0,0 +1,68 @@
classdef ProgressBar < handle
% class for command-line progress-bar notification.
properties
strPercentageLength;
strDotsMaximum;
end
methods
%--- constructor
function this = ProgressBar()
%% Initialization
% Vizualization parameters
this.strPercentageLength = 10; % Length of percentage string (must be >5)
this.strDotsMaximum = 10; % The total number of dots in a progress bar
end
%--- print method
function run(this, msg)
% This function creates a text progress bar. It should be called with a
% STRING argument to initialize and terminate. Otherwise the number corresponding
% to progress in % should be supplied.
% INPUTS: C Either: Text string to initialize or terminate
% Percentage number to show progress
% OUTPUTS: N/A
% Example: Please refer to demo_textprogressbar.m
% Author: Paul Proteus (e-mail: proteus.paul (at) yahoo (dot) com)
% Version: 1.0
% Changes tracker: 29.06.2010 - First version
% Inspired by: http://blogs.mathworks.com/loren/2007/08/01/monitoring-progress-of-a-calculation/
%% Main
persistent strCR; % Carriage return pesistent variable
if isempty(strCR) && ~ischar(msg)
% Progress bar must be initialized with a string
error('The text progress must be initialized with a string!');
elseif isempty(strCR) && ischar(msg)
% Progress bar - initialization
fprintf('\n%s', msg);
strCR = -1;
elseif ~isempty(strCR) && ischar(msg)
% Progress bar - termination
strCR = [];
fprintf([msg '\n']);
elseif isnumeric(msg)
% Progress bar - normal progress
msg = floor(msg);
percentageOut = [num2str(msg) '%%'];
percentageOut = [percentageOut repmat(' ',1,this.strPercentageLength-length(percentageOut)-1)];
nDots = floor(msg/100*this.strDotsMaximum);
dotOut = ['[' repmat('.',1,nDots) repmat(' ',1,this.strDotsMaximum-nDots) ']'];
strOut = [percentageOut dotOut];
% Print it on the screen
if strCR == -1
% Don't do carriage return during first run
fprintf(strOut);
else
% Do it during all the other runs
fprintf([strCR strOut]);
end
% Update carriage return
strCR = repmat('\b',1,length(strOut)-1);
else
% Any other unexpected input
error('Unsupported argument type');
end
end
end
end

74
Data-Analyzer/+Helper/batchAnalyze.m Normal file
View File

@ -0,0 +1,74 @@
function results_all = batchAnalyze(dataSources, options)
arguments
dataSources (1,:) cell
options struct
end
% Default base folder if not specified
if ~isfield(options, 'baseDataFolder')
options.baseDataFolder = '//DyLabNAS/Data';
end
results_all = struct([]); % one element per folder
for i = 1:numel(dataSources)
ds = dataSources{i};
% Use per-sequence baseFolder if present, otherwise default from options
if isfield(ds, 'baseFolder') && ~isempty(ds.baseFolder)
baseFolder = fullfile(ds.baseFolder, ds.sequence, ds.date);
else
baseFolder = fullfile(options.baseDataFolder, ds.sequence, ds.date);
end
for j = 1:numel(ds.runs)
runItem = ds.runs(j);
% Convert numeric or char arrays to a string with leading zeros if needed
if isnumeric(runItem)
runID = sprintf('%04d', runItem); % adjust padding as needed
elseif isstring(runItem)
runID = runItem;
elseif ischar(runItem)
runID = string(runItem);
elseif iscell(runItem)
runID = string(runItem{1}); % handles cell of char
else
error('Unsupported type for run entry: %s', class(runItem));
end
% Build folder path
folderPath = fullfile(baseFolder, runID);
if ~endsWith(folderPath, filesep)
folderPath = [char(folderPath) filesep];
else
folderPath = char(folderPath);
end
options.folderPath = folderPath;
try
% Convert struct -> name-value args
args = [fieldnames(options), struct2cell(options)]';
args = args(:)';
% Perform analysis
analysisResults = Analyzer.performAnalysis(args{:});
% Store flat struct with metadata + results
result = struct();
result.sequence = ds.sequence;
result.date = ds.date;
result.run = runID;
result.path = folderPath;
result.results = analysisResults;
% Append to output
results_all(end+1,1) = result;
catch ME
warning("Error processing %s/%s/%s: %s", ...
ds.sequence, ds.date, runID, ME.message);
end
end
end
end

46
Data-Analyzer/+Helper/calculateODImage.m Normal file
View File

@ -0,0 +1,46 @@
function imageOD = calculateODImage(imageAtom, imageBackground, imageDark, mode, exposureTime)
%CALCULATEODIMAGE Calculates the optical density (OD) image for absorption imaging.
%
% imageOD = calculateODImage(imageAtom, imageBackground, imageDark, mode, exposureTime)
%
% Inputs:
% imageAtom - Image with atoms
% imageBackground - Image without atoms
% imageDark - Image without light
% mode - 'LowIntensity' (default) or 'HighIntensity'
% exposureTime - Required only for 'HighIntensity' [in seconds]
%
% Output:
% imageOD - Computed OD image
%
arguments
imageAtom (:,:) {mustBeNumeric}
imageBackground (:,:) {mustBeNumeric}
imageDark (:,:) {mustBeNumeric}
mode char {mustBeMember(mode, {'LowIntensity', 'HighIntensity'})} = 'LowIntensity'
exposureTime double = NaN
end
% Compute numerator and denominator
numerator = imageBackground - imageDark;
denominator = imageAtom - imageDark;
% Avoid division by zero
numerator(numerator == 0) = 1;
denominator(denominator == 0) = 1;
% Calculate OD based on mode
switch mode
case 'LowIntensity'
imageOD = -log(abs(denominator ./ numerator));
case 'HighIntensity'
if isnan(exposureTime)
error('Exposure time must be provided for HighIntensity mode.');
end
imageOD = abs(denominator ./ numerator);
imageOD = -log(imageOD) + (numerator - denominator) ./ (7000 * (exposureTime / 5e-6));
end
end

137
Data-Analyzer/+Helper/collectODImages.m Normal file
View File

@ -0,0 +1,137 @@
function [ordered_od_imgs, ordered_scan_parameter_values, ordered_file_list] = collectODImages(options)
%% Applies cropping, background subtraction, and optional fringe removal, optional unshuffling on OD image dataset
% Automatically reuses in-memory full dataset if available;
% otherwise, reads and processes raw HDF5 data.
%
% Inputs:
% options - structure containing processing options:
% .folderPath : path to raw HDF5 files
% .saveDirectory : path to save cache (if needed)
% .cam, .angle : camera selection and rotation angle
% .ImagingMode, .PulseDuration : imaging parameters
% .scan_parameter : name of scan parameter
% .center, .span : cropping settings
% .fraction : background subtraction fraction
% .removeFringes : logical flag for fringe removal
% .skipUnshuffling : logical flag to skip unshuffling
% .scan_reference_values: reference values for unshuffling
%
% Outputs:
% ordered_od_imgs : cell array of processed OD images (ordered)
% ordered_scan_parameter_values: vector of scan parameter values (ordered)
% ordered_file_list : cell array of file names (ordered)
% --- Check if the full OD dataset and scan parameters exist in workspace ---
fullDataExists = evalin('base', 'exist(''full_od_imgs'', ''var'')') && ...
evalin('base', 'exist(''full_bkg_imgs'', ''var'')') && ...
evalin('base', 'exist(''raw_scan_parameter_values'', ''var'')') && ...
evalin('base', 'exist(''raw_file_list'', ''var'')');
if fullDataExists
% Both required datasets exist, use them directly
fprintf('\nReusing full OD image dataset and scan parameters from memory.\n');
full_od_imgs = evalin('base', 'full_od_imgs');
full_bkg_imgs = evalin('base', 'full_bkg_imgs');
raw_scan_parameter_values = evalin('base', 'raw_scan_parameter_values');
raw_file_list = evalin('base', 'raw_file_list');
else
% Either dataset is missing, process raw HDF5 files completely
fprintf('\nFull OD image dataset or scan parameters not found in memory.\n');
[full_od_imgs, full_bkg_imgs, raw_scan_parameter_values, raw_file_list] = Helper.processRawData(options);
% Optionally save the full dataset into workspace for future reuse
assignin('base', 'full_od_imgs', full_od_imgs);
assignin('base', 'full_bkg_imgs', full_bkg_imgs);
assignin('base', 'raw_scan_parameter_values', raw_scan_parameter_values);
assignin('base', 'raw_file_list', raw_file_list);
fprintf('\nCompleted computing OD images. Stored in workspace for reuse.\n');
end
nFiles = size(full_od_imgs, 3);
% --- Preallocate arrays for processed images ---
absimages = zeros(options.span(1)+1, options.span(2)+1, nFiles, 'single');
refimages = zeros(options.span(1)+1, options.span(2)+1, nFiles, 'single');
% --- Process each image: crop and subtract background ---
for k = 1:nFiles
od_img = full_od_imgs(:,:,k); % original full OD image, never modified
bkg_img = full_bkg_imgs(:,:,k); % original full background image, never modified
if any(isnan(od_img(:)))
absimages(:,:,k) = nan(options.span(1)+1, options.span(2)+1, 'single');
continue
end
if any(isnan(bkg_img(:)))
refimages(:,:,k) = nan(options.span(1)+1, options.span(2)+1, 'single');
continue
end
% Crop image around the region of interest
cropped_absimage = Helper.cropODImage(od_img, options.center, options.span);
cropped_refimage = Helper.cropODImage(bkg_img, options.center, options.span);
% Subtract background offset based on fraction
processed_absimage = Helper.subtractBackgroundOffset(cropped_absimage, options.fraction);
processed_refimage = Helper.subtractBackgroundOffset(cropped_refimage, options.fraction);
% Store processed image (transpose to match orientation)
absimages(:,:,k) = processed_absimage';
refimages(:,:,k) = processed_refimage';
end
% --- Optional fringe removal ---
if isfield(options, 'removeFringes') && options.removeFringes
fprintf('\nApplying fringe removal to processed images...\n');
optrefimages = Helper.removeFringesInImage(absimages, refimages);
absimages_fringe_removed = absimages - optrefimages;
processed_od_imgs = arrayfun(@(i) absimages_fringe_removed(:,:,i), 1:nFiles, 'UniformOutput', false);
fprintf('\nFringe removal completed.\n');
else
processed_od_imgs = arrayfun(@(i) absimages(:,:,i), 1:nFiles, 'UniformOutput', false);
end
% --- Optional unshuffling based on scan reference values ---
if isfield(options, 'skipUnshuffling') && ~options.skipUnshuffling
fprintf('\nReordering images according to scan parameter reference values...\n');
n_values = length(options.scan_reference_values);
n_total = length(raw_scan_parameter_values);
n_reps = n_total / n_values;
ordered_scan_parameter_values = zeros(1, n_total);
ordered_od_imgs = cell(1, n_total);
ordered_file_list = cell(1, n_total);
counter = 1;
temp_scan_values = raw_scan_parameter_values; % copy for indexing
temp_od_imgs = processed_od_imgs;
temp_file_list = raw_file_list; % copy original file paths for reordering
for rep = 1:n_reps
for val = options.scan_reference_values
idx = find(temp_scan_values == val, 1, 'first');
if isempty(idx), continue; end
ordered_scan_parameter_values(counter) = temp_scan_values(idx);
ordered_od_imgs{counter} = temp_od_imgs{idx};
ordered_file_list{counter} = temp_file_list{idx}; % reorder file list
temp_scan_values(idx) = NaN; % mark as used
temp_od_imgs{idx} = [];
temp_file_list{idx} = [];
counter = counter + 1;
end
end
fprintf('\nImage reordering completed.\n');
else
% No unshuffling: keep original order
ordered_od_imgs = processed_od_imgs;
ordered_scan_parameter_values = raw_scan_parameter_values;
ordered_file_list = raw_file_list;
end
% Optionally save the full dataset into workspace for future reuse
assignin('base', 'od_imgs', ordered_od_imgs);
assignin('base', 'scan_parameter_values', ordered_scan_parameter_values);
fprintf('\nOD image dataset ready for further analysis.\n');
end

18
Data-Analyzer/+Helper/cropODImage.m Normal file
View File

@ -0,0 +1,18 @@
function ret = cropODImage(img, center, span)
% Crop the image according to the region of interest (ROI).
% :param dataSet: The images
% :type dataSet: xarray DataArray or DataSet
% :param center: The center of region of interest (ROI)
% :type center: tuple
% :param span: The span of region of interest (ROI)
% :type span: tuple
% :return: The cropped images
% :rtype: xarray DataArray or DataSet
x_start = floor(center(1) - span(1) / 2);
x_end = floor(center(1) + span(1) / 2);
y_start = floor(center(2) - span(2) / 2);
y_end = floor(center(2) + span(2) / 2);
ret = img(y_start:y_end, x_start:x_end);
end

58
Data-Analyzer/+Helper/drawODOverlays.m Normal file
View File

@ -0,0 +1,58 @@
function drawODOverlays(x1, y1, x2, y2)
% Parameters
tick_spacing = 10; % µm between ticks
tick_length = 2; % µm tick mark length
line_color = [0.5 0.5 0.5];
tick_color = [0.5 0.5 0.5];
font_size = 10;
% Vector from start to end
dx = x2 - x1;
dy = y2 - y1;
L = sqrt(dx^2 + dy^2);
% Unit direction vector along diagonal
ux = dx / L;
uy = dy / L;
% Perpendicular unit vector for ticks
perp_ux = -uy;
perp_uy = ux;
% Midpoint (center)
xc = (x1 + x2) / 2;
yc = (y1 + y2) / 2;
% Number of positive and negative ticks
n_ticks = floor(L / (2 * tick_spacing));
% Draw main diagonal line
plot([x1 x2], [y1 y2], '--', 'Color', line_color, 'LineWidth', 1.2);
for i = -n_ticks:n_ticks
d = i * tick_spacing;
xt = xc + d * ux;
yt = yc + d * uy;
% Tick line endpoints
xt1 = xt - 0.5 * tick_length * perp_ux;
yt1 = yt - 0.5 * tick_length * perp_uy;
xt2 = xt + 0.5 * tick_length * perp_ux;
yt2 = yt + 0.5 * tick_length * perp_uy;
% Draw tick
plot([xt1 xt2], [yt1 yt2], '--', 'Color', tick_color, 'LineWidth', 1);
% Label: centered at tick, offset slightly along diagonal
if d ~= 0
text(xt, yt, sprintf('%+d', d), ...
'Color', tick_color, ...
'FontSize', font_size, ...
'HorizontalAlignment', 'center', ...
'VerticalAlignment', 'bottom', ...
'Rotation', atan2d(dy, dx));
end
end
end

102
Data-Analyzer/+Helper/drawPSOverlays.m Normal file
View File

@ -0,0 +1,102 @@
function drawPSOverlays(kx, ky, k_min, k_max)
% drawPSOverlays - Draw overlays on existing FFT plot:
% - Radial lines every 30°
% - Annular highlight with white (upper half) and gray (lower half) circles at k_min and k_max
% - Horizontal white bands at ky=0 between k_min and k_max
% - Scale ticks and labels every 1 μm⁻¹ along each radial line
%
% Inputs:
% kx, ky - reciprocal space vectors (μm⁻¹)
% k_min - inner annulus radius (μm⁻¹)
% k_max - outer annulus radius (μm⁻¹)
hold on
% === Overlay Radial Lines + Scales ===
max_kx = max(abs(kx));
max_ky = max(abs(ky));
for angle = 0 : pi/6 : pi
x_line = [0, max_kx] * cos(angle);
y_line = [0, max_ky] * sin(angle);
% Plot radial lines
plot(x_line, y_line, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.2);
plot(x_line, -y_line, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.2);
% Draw scale ticks along both lines
drawTicksAlongLine(0,0, x_line(2), y_line(2));
drawTicksAlongLine(0,0, x_line(2), -y_line(2));
end
% === Overlay Annular Highlight ===
theta_full = linspace(0, 2*pi, 500);
% Upper half: white dashed circles
plot(k_min * cos(theta_full(theta_full <= pi)), ...
k_min * sin(theta_full(theta_full <= pi)), 'k--', 'LineWidth', 1.2);
plot(k_max * cos(theta_full(theta_full <= pi)), ...
k_max * sin(theta_full(theta_full <= pi)), 'k--', 'LineWidth', 1.2);
% Lower half: gray dashed circles
plot(k_min * cos(theta_full(theta_full > pi)), ...
k_min * sin(theta_full(theta_full > pi)), '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.0);
plot(k_max * cos(theta_full(theta_full > pi)), ...
k_max * sin(theta_full(theta_full > pi)), '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.0);
% === Highlight horizontal band across k_y = 0 ===
x_vals = kx;
xW1 = x_vals((x_vals >= -k_max) & (x_vals < -k_min));
xW2 = x_vals((x_vals > k_min) & (x_vals <= k_max));
plot(xW1, zeros(size(xW1)), 'k--', 'LineWidth', 1.2);
plot(xW2, zeros(size(xW2)), 'k--', 'LineWidth', 1.2);
hold off
% --- Nested helper function to draw ticks along a radial line ---
function drawTicksAlongLine(x_start, y_start, x_end, y_end)
% Tick parameters
tick_spacing = 1; % spacing between ticks in μm⁻¹
tick_length = 0.05 * sqrt((x_end - x_start)^2 + (y_end - y_start)^2);
tick_color = [0.5 0.5 0.5];
font_size = 8;
% Vector along the line
dx = x_end - x_start;
dy = y_end - y_start;
L = sqrt(dx^2 + dy^2);
ux = dx / L;
uy = dy / L;
% Perpendicular vector for ticks
perp_ux = -uy;
perp_uy = ux;
% Number of ticks
n_ticks = floor(L / tick_spacing);
for i = 1:n_ticks
xt = x_start + i * tick_spacing * ux;
yt = y_start + i * tick_spacing * uy;
% Tick endpoints
xt1 = xt - 0.5 * tick_length * perp_ux;
yt1 = yt - 0.5 * tick_length * perp_uy;
xt2 = xt + 0.5 * tick_length * perp_ux;
yt2 = yt + 0.5 * tick_length * perp_uy;
% Draw tick
plot([xt1 xt2], [yt1 yt2], '-', 'Color', tick_color, 'LineWidth', 1);
% Label
text(xt, yt, sprintf('%d', i), ...
'Color', tick_color, ...
'FontSize', font_size, ...
'HorizontalAlignment', 'center', ...
'VerticalAlignment', 'bottom', ...
'Rotation', atan2d(dy, dx));
end
end
end

11
Data-Analyzer/+Helper/getBkgOffsetFromCorners.m Normal file
View File

@ -0,0 +1,11 @@
function ret = getBkgOffsetFromCorners(img, x_fraction, y_fraction)
% image must be a 2D numerical array
[dim1, dim2] = size(img);
s1 = img(1:round(dim1 * y_fraction), 1:round(dim2 * x_fraction));
s2 = img(1:round(dim1 * y_fraction), round(dim2 - dim2 * x_fraction):dim2);
s3 = img(round(dim1 - dim1 * y_fraction):dim1, 1:round(dim2 * x_fraction));
s4 = img(round(dim1 - dim1 * y_fraction):dim1, round(dim2 - dim2 * x_fraction):dim2);
ret = mean([mean(s1(:)), mean(s2(:)), mean(s3(:)), mean(s4(:))]);
end

90
Data-Analyzer/+Helper/processRawData.m Normal file
View File

@ -0,0 +1,90 @@
function [full_od_imgs, full_bkg_imgs, raw_scan_parameter_values, raw_file_list] = processRawData(options)
%% Reads HDF5 files, computes OD images
%
% Inputs: options.folderPath, options.cam, options.angle, ImagingMode, PulseDuration, scan_parameter, etc.
%
% Returns the OD images and scan parameters immediately in memory.
% This function does NOT do cropping or fringe removal.
fprintf('\nProcessing raw data files at %s ...\n', options.folderPath);
% ===== Group paths in HDF5 files =====
groupList = ["/images/MOT_3D_Camera/in_situ_absorption", ...
"/images/ODT_1_Axis_Camera/in_situ_absorption", ...
"/images/ODT_2_Axis_Camera/in_situ_absorption", ...
"/images/Horizontal_Axis_Camera/in_situ_absorption", ...
"/images/Vertical_Axis_Camera/in_situ_absorption"];
% ===== Find files =====
files = dir(fullfile(options.folderPath, '*.h5'));
nFiles = length(files);
if nFiles == 0
error('\nNo HDF5 files found in %s', options.folderPath);
end
% Determine image size from first file
testFile = fullfile(files(1).folder, files(1).name);
atm_test = double(imrotate(h5read(testFile, append(groupList(options.cam), "/atoms")), options.angle));
[ny, nx] = size(atm_test);
% --- Preallocate in-memory arrays ---
full_od_imgs = nan(ny, nx, nFiles, 'single');
full_bkg_imgs = nan(ny, nx, nFiles, 'single');
raw_scan_parameter_values = zeros(1, nFiles);
% --- Progress bar ---
if isfield(options, 'showProgressBar') && options.showProgressBar
pb = Helper.ProgressBar();
pb.run('Computing OD images | Progress: ');
end
raw_file_list = strings(1, nFiles); % store full file paths
% ===== Loop over files =====
for k = 1:nFiles
fullFileName = fullfile(files(k).folder, files(k).name);
raw_file_list(k) = fullFileName; % track original file
if ~isfield(options, 'showProgressBar') || ~options.showProgressBar
fprintf('Now reading %s\n', fullFileName);
end
try
atm_img = double(imrotate(h5read(fullFileName, append(groupList(options.cam), "/atoms")), options.angle));
bkg_img = double(imrotate(h5read(fullFileName, append(groupList(options.cam), "/background")), options.angle));
dark_img = double(imrotate(h5read(fullFileName, append(groupList(options.cam), "/dark")), options.angle));
od_img = Helper.calculateODImage(atm_img, bkg_img, dark_img, options.ImagingMode, options.PulseDuration);
full_od_imgs(:, :, k) = single(od_img);
full_bkg_imgs(:, :, k) = single(bkg_img);
catch
warning('Missing data in %s, storing NaNs.', fullFileName);
full_od_imgs(:, :, k) = nan(ny, nx, 1, 'single');
full_bkg_imgs(:, :, k) = nan(ny, nx, 1, 'single');
continue;
end
% Extract scan parameter
info = h5info(fullFileName, '/globals');
for i = 1:length(info.Attributes)
if strcmp(info.Attributes(i).Name, options.scan_parameter)
if strcmp(options.scan_parameter, 'ps_rot_mag_fin_pol_angle')
raw_scan_parameter_values(k) = 180 - info.Attributes(i).Value;
else
raw_scan_parameter_values(k) = info.Attributes(i).Value;
end
end
end
% Update progress bar
if isfield(options, 'showProgressBar') && options.showProgressBar
progressPercent = round(k / nFiles * 100);
pb.run(progressPercent);
end
end
% Finish progress bar
if isfield(options, 'showProgressBar') && options.showProgressBar
pb.run(' Done!');
end
end

70
Data-Analyzer/+Helper/removeFringesInImage.m Normal file
View File

@ -0,0 +1,70 @@
function [optrefimages] = removeFringesInImage(absimages, refimages, bgmask)
% removefringesInImage - Fringe removal and noise reduction from absorption images.
% Creates an optimal reference image for each absorption image in a set as
% a linear combination of reference images, with coefficients chosen to
% minimize the least-squares residuals between each absorption image and
% the optimal reference image. The coefficients are obtained by solving a
% linear set of equations using matrix inverse by LU decomposition.
%
% Application of the algorithm is described in C. F. Ockeloen et al, Improved
% detection of small atom numbers through image processing, arXiv:1007.2136 (2010).
%
% Syntax:
% [optrefimages] = removefringesInImage(absimages,refimages,bgmask);
%
% Required inputs:
% absimages - Absorption image data,
% typically 16 bit grayscale images
% refimages - Raw reference image data
% absimages and refimages are both cell arrays containing
% 2D array data. The number of refimages can differ from the
% number of absimages.
%
% Optional inputs:
% bgmask - Array specifying background region used,
% 1=background, 0=data. Defaults to all ones.
% Outputs:
% optrefimages - Cell array of optimal reference images,
% equal in size to absimages.
%
% Dependencies: none
%
% Authors: Shannon Whitlock, Caspar Ockeloen
% Reference: C. F. Ockeloen, A. F. Tauschinsky, R. J. C. Spreeuw, and
% S. Whitlock, Improved detection of small atom numbers through
% image processing, arXiv:1007.2136
% Email:
% May 2009; Last revision: 11 August 2010
% Process inputs
% Set variables, and flatten absorption and reference images
nimgs = size(absimages,3);
nimgsR = size(refimages,3);
xdim = size(absimages(:,:,1),2);
ydim = size(absimages(:,:,1),1);
R = single(reshape(refimages,xdim*ydim,nimgsR));
A = single(reshape(absimages,xdim*ydim,nimgs));
optrefimages=zeros(size(absimages)); % preallocate
if not(exist('bgmask','var')); bgmask=ones(ydim,xdim); end
k = find(bgmask(:)==1); % Index k specifying background region
% Ensure there are no duplicate reference images
% R=unique(R','rows')'; % comment this line if you run out of memory
% Decompose B = R*R' using singular value or LU decomposition
[L,U,p] = lu(R(k,:)'*R(k,:),'vector'); % LU decomposition
for j=1:nimgs
b=R(k,:)'*A(k,j);
% Obtain coefficients c which minimise least-square residuals
lower.LT = true; upper.UT = true;
c = linsolve(U,linsolve(L,b(p,:),lower),upper);
% Compute optimised reference image
optrefimages(:,:,j)=reshape(R*c,[ydim xdim]);
end
end

16
Data-Analyzer/+Helper/subtractBackgroundOffset.m Normal file
View File

@ -0,0 +1,16 @@
function ret = subtractBackgroundOffset(img, fraction)
% Remove the background from the image.
% :param dataArray: The image
% :type dataArray: xarray DataArray
% :param x_fraction: The fraction of the pixels used in x axis
% :type x_fraction: float
% :param y_fraction: The fraction of the pixels used in y axis
% :type y_fraction: float
% :return: The image after removing background
% :rtype: xarray DataArray
x_fraction = fraction(1);
y_fraction = fraction(2);
offset = Helper.getBkgOffsetFromCorners(img, x_fraction, y_fraction);
ret = img - offset;
end

84
Data-Analyzer/+Plotter/compareMultipleDatasets.m Normal file
View File

@ -0,0 +1,84 @@
function compareMultipleDatasets(scanValsCell, meanValsCell, stderrValsCell, varargin)
% compareMultipleDatasets compares multiple datasets with error bars.
%
% Inputs:
% scanValsCell - cell array of x-values for each dataset
% meanValsCell - cell array of mean y-values for each dataset
% stderrValsCell - cell array of std/error values for each dataset
%
% Name-Value Pair Arguments:
% 'FigNum', 'FontName', 'MarkerSize', 'LineWidth', 'CapSize',
% 'YLim', 'Labels', 'Title', 'XLabel', 'YLabel',
% 'SkipSaveFigures', 'SaveFileName', 'SaveDirectory'
% --- Parse inputs ---
p = inputParser;
addParameter(p, 'FigNum', 1, @isnumeric);
addParameter(p, 'FontName', 'Arial', @ischar);
addParameter(p, 'MarkerSize', 6, @isnumeric);
addParameter(p, 'LineWidth', 1.5, @isnumeric);
addParameter(p, 'CapSize', 5, @isnumeric);
addParameter(p, 'YLim', [], @isnumeric);
addParameter(p, 'Labels', {}, @iscell);
addParameter(p, 'Title', '', @ischar);
addParameter(p, 'XLabel', '', @ischar);
addParameter(p, 'YLabel', '', @ischar);
addParameter(p, 'SkipSaveFigures', true, @islogical);
addParameter(p, 'SaveDirectory', pwd, @ischar);
addParameter(p, 'SaveFileName', 'figure.fig', @ischar);
parse(p, varargin{:});
opts = p.Results;
% --- Default labels ---
nDatasets = numel(scanValsCell);
if isempty(opts.Labels)
opts.Labels = arrayfun(@(i) sprintf('Dataset %d',i), 1:nDatasets, 'UniformOutput', false);
end
% --- Marker/line style cycle ---
markerList = {'o', 's', 'd', '^', 'v', '>', '<', 'p', 'h', '*', '+'};
lineList = {'-', '--', ':', '-.'};
% --- Plot ---
fig = figure(opts.FigNum); clf;
set(fig, 'Color', 'w', 'Position', [100 100 950 750]);
hold on;
for i = 1:nDatasets
marker = markerList{mod(i-1, numel(markerList)) + 1};
lineStyle = lineList{mod(i-1, numel(lineList)) + 1};
styleStr = [marker lineStyle];
if isempty(stderrValsCell{i})
plot(scanValsCell{i}, meanValsCell{i}, styleStr, ...
'MarkerSize', opts.MarkerSize, 'LineWidth', opts.LineWidth, ...
'DisplayName', opts.Labels{i});
else
errorbar(scanValsCell{i}, meanValsCell{i}, stderrValsCell{i}, styleStr, ...
'MarkerSize', opts.MarkerSize, 'LineWidth', opts.LineWidth, 'CapSize', opts.CapSize, ...
'DisplayName', opts.Labels{i});
end
end
hold off;
ax = gca;
axisFontSize = 14;
titleFontSize = 16;
set(ax, 'FontName', opts.FontName, 'FontSize', axisFontSize);
if ~isempty(opts.YLim)
ylim(opts.YLim);
end
xlabel(opts.XLabel, 'Interpreter', 'latex', 'FontSize', axisFontSize);
ylabel(opts.YLabel, 'Interpreter', 'latex', 'FontSize', axisFontSize);
title(opts.Title, 'Interpreter', 'latex', 'FontSize', titleFontSize);
legend('Location', 'best');
grid on;
% --- Save figure ---
Plotter.saveFigure(fig, ...
'SaveFileName', opts.SaveFileName, ...
'SaveDirectory', opts.SaveDirectory, ...
'SkipSaveFigures', opts.SkipSaveFigures);
end

126
Data-Analyzer/+Plotter/plotAverageSpectra.m Normal file
View File

@ -0,0 +1,126 @@
function plotAverageSpectra(scan_parameter_values, spectral_analysis_results, varargin)
%% plotAverageSpectra: Plot averaged power, radial, and angular spectra for a scan
%
% Inputs:
% scan_parameter_values - array of scan parameter values
% spectral_analysis_results - struct with fields:
% kx, ky, PS_all, k_rho_vals, S_k_all, theta_vals, S_theta_all
%
% Name-Value Pair Arguments:
% 'ScanParameterName', 'FigNum', 'ColormapPS', 'Font',
% 'SaveFileName', 'SaveDirectory', 'SkipSaveFigures'
% --- Extract data from struct ---
kx = spectral_analysis_results.kx;
ky = spectral_analysis_results.ky;
ps_list = spectral_analysis_results.PS_all;
k_rho_vals = spectral_analysis_results.k_rho_vals;
s_k_list = spectral_analysis_results.S_k_all;
theta_vals = spectral_analysis_results.theta_vals;
s_theta_list = spectral_analysis_results.S_theta_all;
% --- Parse optional parameters ---
p = inputParser;
addParameter(p, 'ScanParameterName', 'ScanParameter', @ischar);
addParameter(p, 'FigNum', 1, @(x) isnumeric(x) && isscalar(x));
addParameter(p, 'ColormapPS', Colormaps.coolwarm(), @(x) isnumeric(x) || ismatrix(x));
addParameter(p, 'Font', 'Arial', @ischar);
addParameter(p, 'SaveFileName', 'figure.fig', @ischar);
addParameter(p, 'SaveDirectory', pwd, @ischar);
addParameter(p, 'SkipSaveFigures', false, @islogical);
parse(p, varargin{:});
opts = p.Results;
scanParam = opts.ScanParameterName;
figNum = opts.FigNum;
colormapPS = opts.ColormapPS;
fontName = opts.Font;
saveFileName = opts.SaveFileName;
saveDirectory = opts.SaveDirectory;
skipSaveFigures = opts.SkipSaveFigures;
% --- Unique scan parameters ---
[uniqueParams, ~, idx] = unique(scan_parameter_values);
nParams = numel(uniqueParams);
% --- Loop over each unique parameter ---
for pIdx = 1:nParams
currentParam = uniqueParams(pIdx);
shotIndices = find(idx == pIdx);
nShots = numel(shotIndices);
% --- Initialize accumulators ---
avgPS = 0;
avgS_k = 0;
avgS_theta = 0;
% --- Sum over shots ---
for j = 1:nShots
avgPS = avgPS + ps_list{shotIndices(j)};
avgS_k = avgS_k + s_k_list{shotIndices(j)};
avgS_theta = avgS_theta + s_theta_list{shotIndices(j)};
end
% --- Average ---
avgPS = avgPS / nShots;
avgS_k = avgS_k / nShots;
avgS_theta = avgS_theta / nShots;
% ==== Plot ====
fig = figure(figNum); clf;
set(fig, 'Color', 'w', 'Position', [400 200 1200 400]);
tLayout = tiledlayout(1,3,'TileSpacing','compact','Padding','compact');
axisFontSize = 14;
titleFontSize = 16;
% --- 1. Power Spectrum ---
nexttile;
imagesc(kx, ky, log(1 + avgPS));
axis image;
set(gca, 'FontSize', axisFontSize, 'YDir', 'normal');
xlabel('k_x [\mum^{-1}]','Interpreter','tex','FontSize',axisFontSize,'FontName',fontName);
ylabel('k_y [\mum^{-1}]','Interpreter','tex','FontSize',axisFontSize,'FontName',fontName);
title('Average Power Spectrum','FontSize',titleFontSize,'FontWeight','bold');
colormap(colormapPS);
colorbar;
% --- Annotate scan parameter ---
if strcmp(scanParam,'ps_rot_mag_fin_pol_angle')
txt = sprintf('%.1f^\\circ', currentParam);
else
txt = sprintf('%.2f G', currentParam);
end
text(0.975,0.975,txt,'Color','white','FontWeight','bold','FontSize',axisFontSize, ...
'Interpreter','tex','Units','normalized','HorizontalAlignment','right','VerticalAlignment','top');
% --- 2. Radial Spectrum ---
nexttile;
plot(k_rho_vals, avgS_k, 'LineWidth', 2);
xlabel('k_\rho [\mum^{-1}]','Interpreter','tex','FontSize',axisFontSize);
ylabel('Magnitude (a.u.)','Interpreter','tex','FontSize',axisFontSize);
title('Average S(k_\rho)','FontSize',titleFontSize,'FontWeight','bold');
set(gca,'FontSize',axisFontSize,'YScale','log','XLim',[min(k_rho_vals), max(k_rho_vals)]);
grid on;
% --- 3. Angular Spectrum ---
nexttile;
plot(theta_vals/pi, avgS_theta, 'LineWidth', 2);
xlabel('\theta/\pi [rad]','Interpreter','tex','FontSize',axisFontSize);
ylabel('Magnitude (a.u.)','Interpreter','tex','FontSize',axisFontSize);
title('Average S(\theta)','FontSize',titleFontSize,'FontWeight','bold');
set(gca,'FontSize',axisFontSize,'YScale','log','YLim',[1e4, 1e7]);
ax = gca;
ax.XMinorGrid = 'on';
ax.YMinorGrid = 'on';
grid on;
drawnow;
% --- Save figure ---
saveFigure(fig, ...
'SaveFileName', saveFileName, ...
'SaveDirectory', saveDirectory, ...
'SkipSaveFigures', skipSaveFigures);
end
end

93
Data-Analyzer/+Plotter/plotCumulants.m Normal file
View File

@ -0,0 +1,93 @@
function plotCumulants(scan_vals, cumulant_data, varargin)
%% plotCumulants: Plots the first four cumulants vs. a scan parameter
%
% Usage:
% plotCumulants(scan_vals, {mean_vals, var_vals, skew_vals, fourth_order_vals}, ...
% 'Title', 'My Title', ...
% 'FigNum', 1, ...
% 'FontName', 'Arial', ...
% 'MarkerSize', 6, ...
% 'LineWidth', 1.5, ...
% 'SkipSaveFigures', false, ...
% 'SaveFileName', 'cumulants.fig', ...
% 'SaveDirectory', pwd);
% --- Parse optional name-value pairs ---
p = inputParser;
addParameter(p, 'Title', '', @ischar);
addParameter(p, 'XLabel', 'Scan Parameter', @ischar);
addParameter(p, 'FigNum', 1, @(x) isnumeric(x) && isscalar(x));
addParameter(p, 'FontName', 'Arial', @ischar);
addParameter(p, 'MarkerSize', 6, @isnumeric);
addParameter(p, 'LineWidth', 1.5, @isnumeric);
addParameter(p, 'SkipSaveFigures', false, @islogical);
addParameter(p, 'SaveFileName', 'cumulants.fig', @ischar);
addParameter(p, 'SaveDirectory', pwd, @ischar);
parse(p, varargin{:});
opts = p.Results;
% --- Extract cumulant data ---
mean_vals = cumulant_data{1};
var_vals = cumulant_data{2};
skew_vals = cumulant_data{3};
fourth_order_vals = cumulant_data{4};
% --- Figure setup ---
fig = figure(opts.FigNum); clf;
set(fig, 'Color', 'w', 'Position', [100 100 950 750]);
axisFontSize = 14;
labelFontSize = 16;
titleFontSize = 16;
tLayout = tiledlayout(2,2,'TileSpacing','compact','Padding','compact');
% --- Mean ---
nexttile;
errorbar(scan_vals, mean_vals, sqrt(var_vals), 'o-', ...
'LineWidth', opts.LineWidth, 'MarkerSize', opts.MarkerSize);
title('Mean', 'FontSize', titleFontSize, 'FontWeight', 'bold');
xlabel(opts.XLabel, 'FontSize', labelFontSize);
ylabel('\kappa_1', 'FontSize', labelFontSize);
set(gca, 'FontSize', axisFontSize, 'FontName', opts.FontName);
grid on;
% --- Variance ---
nexttile;
plot(scan_vals, var_vals, 's-', 'LineWidth', opts.LineWidth, 'MarkerSize', opts.MarkerSize);
title('Variance', 'FontSize', titleFontSize, 'FontWeight', 'bold');
xlabel(opts.XLabel, 'FontSize', labelFontSize);
ylabel('\kappa_2', 'FontSize', labelFontSize);
set(gca, 'FontSize', axisFontSize, 'FontName', opts.FontName);
grid on;
% --- Skewness ---
nexttile;
plot(scan_vals, skew_vals, 'd-', 'LineWidth', opts.LineWidth, 'MarkerSize', opts.MarkerSize);
title('Skewness', 'FontSize', titleFontSize, 'FontWeight', 'bold');
xlabel(opts.XLabel, 'FontSize', labelFontSize);
ylabel('\kappa_3', 'FontSize', labelFontSize);
set(gca, 'FontSize', axisFontSize, 'FontName', opts.FontName);
grid on;
% --- Binder Cumulant ---
nexttile;
plot(scan_vals, fourth_order_vals, '^-', 'LineWidth', opts.LineWidth, 'MarkerSize', opts.MarkerSize);
title('Binder Cumulant', 'FontSize', titleFontSize, 'FontWeight', 'bold');
xlabel(opts.XLabel, 'FontSize', labelFontSize);
ylabel('\kappa_4', 'FontSize', labelFontSize);
set(gca, 'FontSize', axisFontSize, 'FontName', opts.FontName);
grid on;
% --- Super title ---
if ~isempty(opts.Title)
sgtitle(opts.Title, 'FontWeight', 'bold', 'FontSize', titleFontSize, 'Interpreter', 'latex');
end
% --- Save figure ---
Plotter.saveFigure(fig, ...
'SaveFileName', opts.SaveFileName, ...
'SaveDirectory', opts.SaveDirectory, ...
'SkipSaveFigures', opts.SkipSaveFigures);
end

72
Data-Analyzer/+Plotter/plotG2.m Normal file
View File

@ -0,0 +1,72 @@
function plotG2(g2_all, g2_error_all, theta_values, scan_parameter_values, scan_parameter, varargin)
%% plotG2: Plots g2 angular correlations with optional parameters
%
% Usage:
% plotG2(g2_all, g2_error_all, theta_values, unique_scan_parameter_values, scan_parameter, ...
% 'Title', 'My Title', 'XLabel', 'B (G)', 'YLabel', '$g^{(2)}$', ...
% 'FigNum', 1, 'FontName', 'Arial', 'Colormap', @Colormaps.coolwarm, ...
% 'SaveFileName', 'myplot.fig', 'SaveDirectory', 'results')
% --- Parse name-value pairs ---
p = inputParser;
addParameter(p, 'Title', 'g^{(2)}(\delta\theta) vs \delta\theta', @(x) ischar(x) || isstring(x));
addParameter(p, 'XLabel', '$\delta\theta / \pi$', @(x) ischar(x) || isstring(x));
addParameter(p, 'YLabel', '$g^{(2)}(\delta\theta)$', @(x) ischar(x) || isstring(x));
addParameter(p, 'FontName', 'Arial', @ischar);
addParameter(p, 'FontSize', 14, @isnumeric);
addParameter(p, 'Colormap', @parula);
addParameter(p, 'FigNum', [], @(x) isempty(x) || (isnumeric(x) && isscalar(x)));
addParameter(p, 'SkipSaveFigures', false, @islogical);
addParameter(p, 'SaveFileName', 'figure.fig', @ischar);
addParameter(p, 'SaveDirectory', pwd, @ischar);
addParameter(p, 'YLim', [0 1], @isnumeric);
parse(p, varargin{:});
opts = p.Results;
nParams = size(g2_all, 1);
% --- Create figure ---
if isempty(opts.FigNum)
fig = figure;
else
fig = figure(opts.FigNum);
end
clf(fig);
set(fig, 'Color', 'w', 'Position', [100 100 950 750]);
hold on;
% --- Colormap ---
cmap = opts.Colormap(nParams);
% --- Plot data with errorbars ---
legend_entries = cell(nParams, 1);
for i = 1:nParams
errorbar(theta_values/pi, g2_all(i,:), g2_error_all(i,:), ...
'o', 'Color', cmap(i,:), 'MarkerSize', 4, 'MarkerFaceColor', cmap(i,:), 'CapSize', 4);
switch scan_parameter
case 'ps_rot_mag_fin_pol_angle'
legend_entries{i} = sprintf('$\\alpha = %g^\\circ$', scan_parameter_values(i));
case 'rot_mag_field'
legend_entries{i} = sprintf('B = %.2f G', scan_parameter_values(i));
otherwise
legend_entries{i} = sprintf('%g', scan_parameter_values(i));
end
end
% --- Formatting ---
xlabel(opts.XLabel, 'Interpreter', 'latex', 'FontName', opts.FontName, 'FontSize', opts.FontSize);
ylabel(opts.YLabel, 'Interpreter', 'latex', 'FontName', opts.FontName, 'FontSize', opts.FontSize);
title(opts.Title, 'Interpreter', 'latex', 'FontName', opts.FontName, 'FontSize', opts.FontSize + 2);
legend(legend_entries, 'Interpreter', 'latex', 'Location', 'bestoutside');
set(gca, 'FontName', opts.FontName, 'FontSize', opts.FontSize);
ylim(opts.YLim);
grid on;
% --- Save figure ---
Plotter.saveFigure(fig, ...
'SaveFileName', opts.SaveFileName, ...
'SaveDirectory', opts.SaveDirectory, ...
'SkipSaveFigures', opts.SkipSaveFigures);
end

69
Data-Analyzer/+Plotter/plotHeatmap.m Normal file
View File

@ -0,0 +1,69 @@
function plotHeatmap(results_all, x_values, y_values, fieldName, varargin)
%% plotHeatmap: Plots a heatmap for a field in a struct array.
%
% Usage:
% plotHeatmap(results_all, x_values, y_values, fieldName, ...
% 'FigNum', 1, 'Colormap', parula, 'CLim', [0 1], ...
% 'XLabel', '\alpha (degrees)', 'YLabel', 'BField (G)', ...
% 'Title', 'My Title', 'SaveFileName', 'heatmap.fig', ...
% 'SaveDirectory', 'results', 'SkipSaveFigures', false);
% --- Parse optional inputs ---
p = inputParser;
addParameter(p, 'FigNum', []);
addParameter(p, 'Colormap', parula);
addParameter(p, 'CLim', []);
addParameter(p, 'XLabel', '\alpha (degrees)');
addParameter(p, 'YLabel', 'BField (G)');
addParameter(p, 'Title', fieldName);
addParameter(p, 'FontName', 'Arial');
addParameter(p, 'FontSize', 14);
addParameter(p, 'SkipSaveFigures', false, @islogical);
addParameter(p, 'SaveFileName', 'heatmap.fig', @ischar);
addParameter(p, 'SaveDirectory', pwd, @ischar);
parse(p, varargin{:});
opts = p.Results;
N_y = length(results_all);
N_x = length(x_values);
% --- Preallocate data matrix ---
data_matrix = NaN(N_y, N_x);
for i = 1:N_y
if isfield(results_all(i), fieldName)
data_matrix(i, :) = results_all(i).(fieldName);
else
warning('Field "%s" does not exist in results_all(%d). Filling with NaN.', fieldName, i);
end
end
% --- Create figure ---
if isempty(opts.FigNum)
fig = figure;
else
fig = figure(opts.FigNum);
end
clf(fig);
set(fig, 'Color', 'w', 'Position', [50 50 950 750]);
% --- Plot heatmap ---
imagesc(x_values, y_values, data_matrix);
colormap(opts.Colormap);
if ~isempty(opts.CLim)
caxis(opts.CLim);
end
set(gca, 'FontName', opts.FontName, 'FontSize', opts.FontSize, 'YDir', 'normal');
% --- Labels and title ---
xlabel(opts.XLabel, 'Interpreter', 'tex', 'FontName', opts.FontName, 'FontSize', opts.FontSize);
ylabel(opts.YLabel, 'Interpreter', 'tex', 'FontName', opts.FontName, 'FontSize', opts.FontSize);
title(opts.Title, 'Interpreter', 'latex', 'FontSize', opts.FontSize + 2, 'FontWeight', 'bold');
colorbar;
% --- Save figure ---
Plotter.saveFigure(fig, ...
'SaveFileName', opts.SaveFileName, ...
'SaveDirectory', opts.SaveDirectory, ...
'SkipSaveFigures', opts.SkipSaveFigures);
end

70
Data-Analyzer/+Plotter/plotMeanWithSE.m Normal file
View File

@ -0,0 +1,70 @@
function plotMeanWithSE(scan_values, data_values, varargin)
%% plotMeanWithSE: Plots mean ± standard error vs a scan parameter.
%
% Usage:
% plotMeanWithSE(scan_values, data_values, ...
% 'Title', 'My Title', 'XLabel', 'Parameter', 'YLabel', 'Mean Value', ...
% 'FigNum', 1, 'FontName', 'Arial', 'YLim', [0 1], ...
% 'SaveFileName', 'mean_with_se.fig', 'SaveDirectory', 'results', ...
% 'SkipSaveFigures', false);
% --- Parse optional name-value pairs ---
p = inputParser;
addParameter(p, 'Title', '', @(x) ischar(x) || isstring(x));
addParameter(p, 'XLabel', '', @(x) ischar(x) || isstring(x));
addParameter(p, 'YLabel', '', @(x) ischar(x) || isstring(x));
addParameter(p, 'FigNum', [], @(x) isempty(x) || (isnumeric(x) && isscalar(x)));
addParameter(p, 'FontName', 'Arial', @ischar);
addParameter(p, 'FontSize', 14, @isnumeric);
addParameter(p, 'YLim', [], @(x) isempty(x) || isnumeric(x));
addParameter(p, 'SkipSaveFigures', false, @islogical);
addParameter(p, 'SaveFileName', 'mean_with_se.fig', @ischar);
addParameter(p, 'SaveDirectory', pwd, @ischar);
parse(p, varargin{:});
opts = p.Results;
% --- Compute mean and standard error ---
[unique_vals, ~, idx] = unique(scan_values);
mean_vals = zeros(size(unique_vals));
stderr_vals = zeros(size(unique_vals));
for i = 1:length(unique_vals)
group = data_values(idx == i);
if iscell(group)
groupVals = [group{:}];
else
groupVals = group;
end
mean_vals(i) = mean(groupVals);
stderr_vals(i) = std(groupVals) / sqrt(length(groupVals));
end
% --- Create figure ---
if isempty(opts.FigNum)
fig = figure;
else
fig = figure(opts.FigNum);
end
clf(fig);
set(fig, 'Color', 'w', 'Position', [100 100 950 750]);
% --- Plot error bars ---
errorbar(unique_vals, mean_vals, stderr_vals, 'o--', ...
'LineWidth', 1.8, 'MarkerSize', 6, 'CapSize', 5);
% --- Axis formatting ---
set(gca, 'FontName', opts.FontName, 'FontSize', opts.FontSize);
if ~isempty(opts.YLim)
ylim(opts.YLim);
end
xlabel(opts.XLabel, 'Interpreter', 'latex', 'FontName', opts.FontName, 'FontSize', opts.FontSize);
ylabel(opts.YLabel, 'Interpreter', 'latex', 'FontName', opts.FontName, 'FontSize', opts.FontSize);
title(opts.Title, 'Interpreter', 'latex', 'FontSize', opts.FontSize + 2, 'FontWeight', 'bold');
grid on;
% --- Save figure ---
Plotter.saveFigure(fig, ...
'SaveFileName', opts.SaveFileName, ...
'SaveDirectory', opts.SaveDirectory, ...
'SkipSaveFigures', opts.SkipSaveFigures);
end

81
Data-Analyzer/+Plotter/plotPDF.m Normal file
View File

@ -0,0 +1,81 @@
function plotPDF(dataCell, referenceValues, varargin)
%% plotPDF: Plots 2D heatmap of PDFs for grouped data
%
% Usage:
% Plotter.plotPDF(dataCell, referenceValues, ...
% 'Title', 'My Title', ...
% 'XLabel', 'Scan Parameter', ...
% 'YLabel', 'Data Values', ...
% 'FigNum', 1, ...
% 'FontName', 'Arial', ...
% 'SkipSaveFigures', true, ...
% 'SaveFileName', 'SavedPDFs', ...
% 'SaveDirectory', 'results', ...
% 'NumPoints', 200, ...
% 'DataRange', [min max], ...
% 'XLim', [xmin xmax], ...
% 'Colormap', @jet);
% --- Parse optional inputs ---
p = inputParser;
addParameter(p, 'Title', '', @(x) ischar(x) || isstring(x));
addParameter(p, 'XLabel', '', @(x) ischar(x) || isstring(x));
addParameter(p, 'YLabel', '', @(x) ischar(x) || isstring(x));
addParameter(p, 'FigNum', 1, @(x) isscalar(x));
addParameter(p, 'FontName', 'Arial', @ischar);
addParameter(p, 'FontSize', 14, @isnumeric);
addParameter(p, 'SkipSaveFigures', false, @islogical);
addParameter(p, 'SaveFileName', 'pdf.fig', @ischar);
addParameter(p, 'SaveDirectory', pwd, @ischar);
addParameter(p, 'NumPoints', 200, @(x) isscalar(x));
addParameter(p, 'DataRange', [], @(x) isempty(x) || numel(x)==2);
addParameter(p, 'XLim', [], @(x) isempty(x) || numel(x)==2);
addParameter(p, 'Colormap', @jet);
parse(p, varargin{:});
opts = p.Results;
N_params = numel(referenceValues);
% --- Determine y-grid for PDF ---
if isempty(opts.DataRange)
allData = cell2mat(dataCell(:));
y_grid = linspace(min(allData), max(allData), opts.NumPoints);
else
y_grid = linspace(opts.DataRange(1), opts.DataRange(2), opts.NumPoints);
end
pdf_matrix = zeros(numel(y_grid), N_params);
% --- Compute PDFs ---
for i = 1:N_params
data = dataCell{i};
data = data(~isnan(data));
if isempty(data), continue; end
f = ksdensity(data, y_grid);
pdf_matrix(:, i) = f;
end
% --- Plot heatmap ---
fig = figure(opts.FigNum); clf(fig);
set(fig, 'Color', 'w', 'Position',[100 100 950 750]);
imagesc(referenceValues, y_grid, pdf_matrix);
set(gca, 'YDir', 'normal', 'FontName', opts.FontName, 'FontSize', opts.FontSize);
xlabel(opts.XLabel, 'Interpreter', 'latex', 'FontSize', opts.FontSize, 'FontName', opts.FontName);
ylabel(opts.YLabel, 'Interpreter', 'latex', 'FontSize', opts.FontSize, 'FontName', opts.FontName);
title(opts.Title, 'Interpreter', 'latex', 'FontSize', opts.FontSize + 2, 'FontWeight', 'bold');
colormap(feval(opts.Colormap));
c = colorbar;
ylabel(c, 'PDF', 'Interpreter', 'latex', 'FontSize', opts.FontSize);
if ~isempty(opts.XLim)
xlim(opts.XLim);
end
% --- Save figure ---
Plotter.saveFigure(fig, ...
'SaveFileName', opts.SaveFileName, ...
'SaveDirectory', opts.SaveDirectory, ...
'SkipSaveFigures', opts.SkipSaveFigures);
end

50
Data-Analyzer/+Plotter/saveFigure.m Normal file
View File

@ -0,0 +1,50 @@
function saveFigure(fig, varargin)
%% saveFigure saves a MATLAB figure as a .fig file in a specified directory.
%
% Usage:
% saveFigure(fig)
% saveFigure(fig, 'SaveFileName', 'myplot.fig', 'SaveDirectory', 'results', 'SkipSaveFigures', false)
%
% Inputs:
% fig - Figure handle to save
%
% Optional Parameters:
% 'SaveFileName' - Name of the file (default: 'figure.fig')
% 'SaveDirectory' - Directory to save into (default: current working directory)
% 'SkipSaveFigures' - If true, skips saving (default: false)
%
% Example:
% fig = figure;
% plot(1:10, rand(1,10));
% saveFigure(fig, 'SaveFileName', 'test.fig', 'SaveDirectory', 'plots');
% --- Defaults ---
p = inputParser;
addParameter(p, 'SaveFileName', 'figure.fig');
addParameter(p, 'SaveDirectory', pwd);
addParameter(p, 'SkipSaveFigures', false);
parse(p, varargin{:});
opts = p.Results;
if opts.SkipSaveFigures
return; % Do nothing
end
% --- Ensure directory exists ---
if ~exist(opts.SaveDirectory, 'dir')
mkdir(opts.SaveDirectory);
end
% --- Ensure .fig extension ---
[~, name, ext] = fileparts(opts.SaveFileName);
if isempty(ext)
ext = '.fig';
elseif ~strcmpi(ext, '.fig')
warning('Overriding extension to .fig (was %s).', ext);
ext = '.fig';
end
saveFullPath = fullfile(opts.SaveDirectory, [name ext]);
savefig(fig, saveFullPath);
fprintf('Figure saved as MATLAB .fig: %s\n', saveFullPath);
end

159
Data-Analyzer/+Scripts/BECToDroplets/plotAnalysisResults.m Normal file
View File

@ -0,0 +1,159 @@
%% ------------------ 1. Mean ± Std Plots ------------------
% Plot Radial Spectral Contrast
Plotter.plotMeanWithSE(scan_parameter_values, results_all.spectral_analysis_results.radial_spectral_contrast, ...
'Title', options.titleString, ...
'XLabel', 'B (G)', ...
'YLabel', 'Radial Spectral Contrast', ...
'FigNum', 1, ...
'FontName', options.font, ...
'SaveFileName', 'RadialSpectralContrast.fig', ...
'SaveDirectory', [options.saveDirectory '/Results'], ...
'SkipSaveFigures', options.skipSaveFigures);
% Plot Angular Spectral Weight
Plotter.plotMeanWithSE(scan_parameter_values, results_all.spectral_analysis_results.angular_spectral_weight, ...
'Title', options.titleString, ...
'XLabel', 'B (G)', ...
'YLabel', 'Angular Spectral Weight', ...
'FigNum', 2, ...
'FontName', options.font, ...
'SaveFileName', 'AngularSpectralWeight.fig', ...
'SaveDirectory', [options.saveDirectory '/Results'], ...
'SkipSaveFigures', options.skipSaveFigures);
% Plot Peak Offset Angular Correlation
Plotter.plotMeanWithSE(options.scan_reference_values, results_all.custom_g_results.max_g2_all_per_scan_parameter_value, ...
'Title', options.titleString, ...
'XLabel', 'B (G)', ...
'YLabel', '$\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', ...
'FigNum', 3, ...
'YLim', [0 1], ...
'FontName', options.font, ...
'SaveFileName', 'PeakOffsetAngularCorrelation.fig', ...
'SaveDirectory', [options.saveDirectory '/Results'], ...
'SkipSaveFigures', options.skipSaveFigures);
%% ------------------ 2. g²(θ) across transition ------------------
Plotter.plotG2(results_all.full_g2_results.g2_all, ...
results_all.full_g2_results.g2_error_all, ...
results_all.full_g2_results.theta_values, ...
options.scan_reference_values, ...
'rot_mag_field', ...
'Title', options.titleString, ...
'XLabel', '$\delta\theta / \pi$', ...
'YLabel', '$g^{(2)}(\delta\theta)$', ...
'FigNum', 4, ...
'FontName', options.font, ...
'SkipSaveFigures', options.skipSaveFigures, ...
'SaveFileName', 'G2ThetaAcrossTransition.fig', ...
'SaveDirectory', [options.saveDirectory '/Results'], ...
'Colormap', @Colormaps.coolwarm);
%% ------------------ 3. PDF of max g² across transition ------------------
Plotter.plotPDF(results_all.custom_g_results.max_g2_all_per_scan_parameter_value, options.scan_reference_values, ...
'Title', options.titleString, ...
'XLabel', 'B (G)', ...
'YLabel', '$\mathrm{max}[g^{(2)}]$', ...
'FigNum', 5, ...
'FontName', options.font, ...
'SkipSaveFigures', options.skipSaveFigures, ...
'SaveFileName', 'PDF_MaxG2AcrossTransition.fig', ...
'SaveDirectory', [options.saveDirectory '/Results'], ...
'NumPoints', 200, ...
'DataRange', [0 1.5], ...
'Colormap', @Colormaps.coolwarm, ...
'XLim', [min(options.scan_reference_values) max(options.scan_reference_values)]);
%% ------------------ 4. Cumulants across transition ------------------
Plotter.plotCumulants(options.scan_reference_values, ...
{results_all.custom_g_results.mean_max_g2, results_all.custom_g_results.var_max_g2, results_all.custom_g_results.skew_max_g2_angle, results_all.custom_g_results.fourth_order_cumulant_max_g2}, ...
'Title', 'Cumulants of Peak Offset Angular Correlation', ...
'XLabel', 'B (G)', ...
'FigNum', 6, ...
'FontName', options.font, ...
'MarkerSize', 6, ...
'LineWidth', 1.5, ...
'SkipSaveFigures', options.skipSaveFigures, ...
'SaveFileName', 'CumulantOfPeakOffsetAngularCorrelation.fig', ...
'SaveDirectory', [options.saveDirectory '/Results']);
%{
%% ------------------ 6. Average of Spectra Plots ------------------
Plotter.plotAverageSpectra(scan_parameter_values, ...
spectral_analysis_results, ...
'ScanParameterName', scan_parameter, ...
'FigNum', 7, ...
'ColormapPS', Colormaps.coolwarm(), ...
'Font', 'Bahnschrift', ...
'SaveFileName', 'avgSpectra.fig', ...
'SaveDirectory', [options.saveDirectory '/Results'], ...
'SkipSaveFigures', options.skipSaveFigures);
%% ------------------ 7. Compare quantities ------------------
% Load Droplets → Stripes data
Data = load(dtsFile, ...
'unique_scan_parameter_values', ...
'mean_max_g2_values', ...
'std_error_g2_values');
dts_scan_parameter_values = Data.unique_scan_parameter_values;
dts_mean_mg2 = Data.mean_max_g2_values;
dts_stderr_mg2 = Data.std_error_g2_values;
% Load Stripes → Droplets data
Data = load(stdFile, ...
'unique_scan_parameter_values', ...
'mean_max_g2_values', ...
'std_error_g2_values');
std_scan_parameter_values = Data.unique_scan_parameter_values;
std_mean_mg2 = Data.mean_max_g2_values;
std_stderr_mg2 = Data.std_error_g2_values;
% Prepare cell arrays for multiple datasets
scanValsCell = {dts_scan_parameter_values, std_scan_parameter_values};
meanValsCell = {dts_mean_mg2, std_mean_mg2};
stderrValsCell = {dts_stderr_mg2, std_stderr_mg2};
% Compare datasets
compareMultipleDatasets(scanValsCell, meanValsCell, stderrValsCell, ...
'FigNum', 8, ...
'FontName', 'Bahnschrift', ...
'MarkerSize', 6, ...
'LineWidth', 1.5, ...
'CapSize', 5, ...
'YLim', [0 1], ...
'Labels', {'Droplets → Stripes', 'Stripes → Droplets'}, ...
'Title', 'AngularCorrelation_Comparison', ...
'XLabel', 'B (G)', ...
'YLabel', '$\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', ...
'SkipSaveFigures', options.skipSaveFigures, ...
'SaveDirectory', [options.saveDirectory '/Results'], ...
'SaveFileName', 'AngularCorrelation_Comparison.fig');
%% ------------------ 8. Heatmaps ------------------
BFields = [2.35, 2.15, 2.0, 1.85, 1.7, 1.55, 1.4, 1.35];
% Heatmap of mean_max_g2_values
Plotter.plotHeatmap(results_all, options.scan_groups, BFields, 'mean_max_g2_values', ...
'Colormap', @sky, ...
'CLim', [0 1], ...
'XLabel', '\alpha (degrees)', ...
'YLabel', 'BField (G)', ...
'Title', '$\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', ...
'FigNum', 9, ...
'SaveFileName', 'Heatmap_MaxG2.fig', ...
'SaveDirectory', options.resultsDir);
% Heatmap of radial_spectral_contrast
Plotter.plotHeatmap(results_all, options.scan_groups, BFields, 'radial_spectral_contrast', ...
'Colormap', @sky, ...
'CLim', [0 0.008], ...
'XLabel', '\alpha (degrees)', ...
'YLabel', 'BField (G)', ...
'Title', 'Radial Spectral Contrast', ...
'FigNum', 10, ...
'SaveFileName', 'Heatmap_RadialSpectralContrast.fig', ...
'SaveDirectory', options.resultsDir);
%}

75
Data-Analyzer/+Scripts/BECToDroplets/plotImages.m Normal file
View File

@ -0,0 +1,75 @@
%% ===== BEC-Droplets Settings =====
options = struct();
% File / paths
options.folderPath = "//DyLabNAS/Data/StructuralPhaseTransition/2025/08/13/0062";
options.savefileName = 'BECToDroplets';
options.saveDirectory = "Z:/Users/Karthik/Data-Analyzer/+Scripts";
% Camera / imaging
options.cam = 5;
options.angle = 0;
options.center = [1420, 2050];
options.span = [200, 200];
options.fraction = [0.1, 0.1];
options.pixel_size = 5.86e-6; % in meters
options.magnification = 23.94;
options.removeFringes = false;
options.ImagingMode = 'HighIntensity';
options.PulseDuration = 5e-6; % in s
% Fourier analysis settings
% Radial Spectral Distribution
options.theta_min = deg2rad(0);
options.theta_max = deg2rad(180);
options.N_radial_bins = 500;
options.Radial_Sigma = 2;
options.Radial_WindowSize = 5; % odd number for centered moving avg
% Angular Spectral Distribution
options.k_min = 1.2771; % in μm⁻¹
options.k_max = 2.5541; % in μm⁻¹
options.N_angular_bins = 180;
options.Angular_Threshold = 75;
options.Angular_Sigma = 2;
options.Angular_WindowSize = 5;
options.zoom_size = 50; % zoomed-in region around center
% Scan parameter
options.scan_parameter = 'rot_mag_field';
if strcmp(options.savefileName, 'BECToDroplets')
options.scan_reference_values = [2.40, 2.39, 2.38, 2.37, 2.35, 2.34, 2.32, 2.30, 2.28, 2.26, 2.24, 2.22, 2.2, 2.15, 2.10, 2.05, 2, 1.95, 1.90, 1.85, 1.8];
options.titleString = 'BEC to Droplets';
elseif strcmp(options.savefileName, 'BECToStripes')
options.scan_reference_values = [2.45, 2.44, 2.43, 2.42, 2.4, 2.39, 2.38, 2.37, 2.36, 2.35, 2.34, 2.32, 2.3, 2.28, 2.25, 2.2, 2.15, 2.10, 2.0, 1.90, 1.8];
options.titleString = 'BEC to Stripes';
elseif strcmp(options.savefileName, 'DropletsToStripes')
options.scan_reference_values = [0, 5, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 35, 40];
options.titleString = 'Droplets to Stripes';
elseif strcmp(options.savefileName, 'StripesToDroplets')
options.scan_reference_values = fliplr([0, 5, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 35, 40]);
options.titleString = 'Stripes to Droplets';
end
% Flags
options.skipNormalization = true;
options.skipUnshuffling = false;
options.skipPreprocessing = true;
options.skipMasking = true;
options.skipIntensityThresholding = true;
options.skipBinarization = true;
options.skipMovieRender = true;
options.skipSaveFigures = true;
options.skipSaveOD = true;
options.skipLivePlot = false;
options.showProgressBar = true;
% Optional extras
options.font = 'Bahnschrift';
%%
[od_imgs, scan_parameter_values, file_list] = Helper.collectODImages(options);
%%
Analyzer.runInteractiveODImageViewer(od_imgs, scan_parameter_values, file_list, options);

80
Data-Analyzer/+Scripts/BECToDroplets/runFullAnalysis.m Normal file
View File

@ -0,0 +1,80 @@
%% ===== BEC-Droplets Settings =====
% Specify data location to run analysis on
dataSources = {
struct('sequence', 'StructuralPhaseTransition', ...
'date', '2025/08/13', ...
'runs', [62]) % specify run numbers as a string in "" or just as a numeric value
};
options = struct();
% File / paths
options.baseDataFolder = '//DyLabNAS/Data';
options.savefileName = 'BECToDroplets';
scriptFullPath = mfilename('fullpath');
options.saveDirectory = fileparts(scriptFullPath);
% Camera / imaging
options.cam = 5;
options.angle = 0;
options.center = [1420, 2050];
options.span = [200, 200];
options.fraction = [0.1, 0.1];
options.pixel_size = 5.86e-6; % in meters
options.magnification = 23.94;
options.removeFringes = false;
options.ImagingMode = 'HighIntensity';
options.PulseDuration = 5e-6; % in s
% Fourier analysis settings
options.theta_min = deg2rad(0);
options.theta_max = deg2rad(180);
options.N_radial_bins = 500;
options.Radial_Sigma = 2;
options.Radial_WindowSize = 5; % odd number
options.k_min = 1.2771; % μm⁻¹
options.k_max = 2.5541; % μm⁻¹
options.N_angular_bins = 180;
options.Angular_Threshold = 75;
options.Angular_Sigma = 2;
options.Angular_WindowSize = 5;
options.zoom_size = 50;
% Scan parameter
options.scan_parameter = 'rot_mag_field';
switch options.savefileName
case 'BECToDroplets'
options.scan_reference_values = [2.40, 2.39, 2.38, 2.37, 2.35, 2.34, 2.32, 2.30, 2.28, 2.26, 2.24, 2.22, 2.2, 2.15, 2.10, 2.05, 2, 1.95, 1.90, 1.85, 1.8];
options.titleString = 'BEC to Droplets';
case 'BECToStripes'
options.scan_reference_values = [2.45, 2.44, 2.43, 2.42, 2.4, 2.39, 2.38, 2.37, 2.36, 2.35, 2.34, 2.32, 2.3, 2.28, 2.25, 2.2, 2.15, 2.10, 2.0, 1.90, 1.8];
options.titleString = 'BEC to Stripes';
case 'DropletsToStripes'
options.scan_reference_values = [0, 5, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 35, 40];
options.titleString = 'Droplets to Stripes';
case 'StripesToDroplets'
options.scan_reference_values = fliplr([0, 5, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 35, 40]);
options.titleString = 'Stripes to Droplets';
end
% Flags
options.skipNormalization = false;
options.skipUnshuffling = false;
options.skipPreprocessing = true;
options.skipMasking = true;
options.skipIntensityThresholding = true;
options.skipBinarization = true;
options.skipMovieRender = true;
options.skipSaveFigures = true;
options.skipSaveOD = true;
options.skipLivePlot = false;
options.showProgressBar = true;
% Extras
options.font = 'Bahnschrift';
%% ===== Run Batch Analysis =====
results_all = Helper.batchAnalyze(dataSources, options);

711
Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/analyzeFolder.m
View File

@ -1,711 +0,0 @@
function results = analyzeFolder(options)
% Ensure required fields are defined in options
arguments
options.scan_parameter (1,:) char
options.scan_groups (1,:) double
options.cam (1,1) double
options.angle (1,1) double
options.center (1,2) double
options.span (1,2) double
options.fraction (1,2) double
options.ImagingMode (1,:) char
options.PulseDuration (1,1) double
options.removeFringes (1,1) logical
options.skipUnshuffling (1,1) logical
options.pixel_size (1,1) double
options.magnification (1,1) double
options.zoom_size (1,1) double
options.r_min (1,1) double
options.r_max (1,1) double
options.N_angular_bins (1,1) double
options.Angular_Threshold (1,1) double
options.Angular_Sigma (1,1) double
options.Angular_WindowSize (1,1) double
options.theta_min (1,1) double
options.theta_max (1,1) double
options.N_radial_bins (1,1) double
options.Radial_Sigma (1,1) double
options.Radial_WindowSize (1,1) double
options.k_min (1,1) double
options.k_max (1,1) double
options.skipPreprocessing (1,1) logical
options.skipMasking (1,1) logical
options.skipIntensityThresholding (1,1) logical
options.skipBinarization (1,1) logical
options.folderPath (1,:) char
end
% Assign variables from options
scan_parameter = options.scan_parameter;
scan_groups = options.scan_groups;
folderPath = options.folderPath;
center = options.center;
span = options.span;
fraction = options.fraction;
ImagingMode = options.ImagingMode;
PulseDuration = options.PulseDuration;
removeFringes = options.removeFringes;
skipUnshuffling = options.skipUnshuffling;
pixel_size = options.pixel_size;
magnification = options.magnification;
zoom_size = options.zoom_size;
r_min = options.r_min;
r_max = options.r_max;
N_angular_bins = options.N_angular_bins;
Angular_Threshold = options.Angular_Threshold;
Angular_Sigma = options.Angular_Sigma;
Angular_WindowSize = options.Angular_WindowSize;
theta_min = options.theta_min;
theta_max = options.theta_max;
N_radial_bins = options.N_radial_bins;
Radial_Sigma = options.Radial_Sigma;
Radial_WindowSize = options.Radial_WindowSize;
k_min = options.k_min;
k_max = options.k_max;
skipPreprocessing = options.skipPreprocessing;
skipMasking = options.skipMasking;
skipIntensityThresholding = options.skipIntensityThresholding;
skipBinarization = options.skipBinarization;
cam = options.cam;
angle = options.angle;
% Load images and analyze them (keep using the cleaned body of your original function)
% Fix the incorrect usage of 'cam' and 'angle' not defined locally
groupList = ["/images/MOT_3D_Camera/in_situ_absorption",
"/images/ODT_1_Axis_Camera/in_situ_absorption",
"/images/ODT_2_Axis_Camera/in_situ_absorption",
"/images/Horizontal_Axis_Camera/in_situ_absorption",
"/images/Vertical_Axis_Camera/in_situ_absorption"];
filePattern = fullfile(folderPath, '*.h5');
files = dir(filePattern);
refimages = zeros(span(1) + 1, span(2) + 1, length(files));
absimages = zeros(span(1) + 1, span(2) + 1, length(files));
for k = 1 : length(files)
baseFileName = files(k).name;
fullFileName = fullfile(files(k).folder, baseFileName);
fprintf(1, 'Now reading %s\n', fullFileName);
atm_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/atoms")), angle));
bkg_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/background")), angle));
dark_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/dark")), angle));
if (isempty(atm_img) && isa(atm_img, 'double')) || ...
(isempty(bkg_img) && isa(bkg_img, 'double')) || ...
(isempty(dark_img) && isa(dark_img, 'double'))
refimages(:,:,k) = nan(size(refimages(:,:,k))); % fill with NaNs
absimages(:,:,k) = nan(size(absimages(:,:,k)));
else
refimages(:,:,k) = subtractBackgroundOffset(cropODImage(bkg_img, center, span), fraction)';
absimages(:,:,k) = subtractBackgroundOffset(cropODImage(calculateODImage(atm_img, bkg_img, dark_img, ImagingMode, PulseDuration), center, span), fraction)';
end
end
% ===== Fringe removal =====
if removeFringes
optrefimages = removefringesInImage(absimages, refimages);
absimages_fringe_removed = absimages(:, :, :) - optrefimages(:, :, :);
nimgs = size(absimages_fringe_removed,3);
od_imgs = cell(1, nimgs);
for i = 1:nimgs
od_imgs{i} = absimages_fringe_removed(:, :, i);
end
else
nimgs = size(absimages(:, :, :),3);
od_imgs = cell(1, nimgs);
for i = 1:nimgs
od_imgs{i} = absimages(:, :, i);
end
end
% ===== Get rotation angles =====
scan_parameter_values = zeros(1, length(files));
% Get information about the '/globals' group
for k = 1 : length(files)
baseFileName = files(k).name;
fullFileName = fullfile(files(k).folder, baseFileName);
info = h5info(fullFileName, '/globals');
for i = 1:length(info.Attributes)
if strcmp(info.Attributes(i).Name, scan_parameter)
if strcmp(scan_parameter, 'rot_mag_fin_pol_angle')
scan_parameter_values(k) = 180 - info.Attributes(i).Value;
else
scan_parameter_values(k) = info.Attributes(i).Value;
end
end
end
end
% ===== Unshuffle if necessary to do so =====
if ~skipUnshuffling
n_values = length(scan_groups);
n_total = length(scan_parameter_values);
% Infer number of repetitions
n_reps = n_total / n_values;
% Preallocate ordered arrays
ordered_scan_values = zeros(1, n_total);
ordered_od_imgs = cell(1, n_total);
counter = 1;
for rep = 1:n_reps
for val = scan_groups
% Find the next unused match for this val
idx = find(scan_parameter_values == val, 1, 'first');
% Assign and remove from list to avoid duplicates
ordered_scan_values(counter) = scan_parameter_values(idx);
ordered_od_imgs{counter} = od_imgs{idx};
% Mark as used by removing
scan_parameter_values(idx) = NaN; % NaN is safe since original values are 0:5:45
od_imgs{idx} = []; % empty cell so it won't be matched again
counter = counter + 1;
end
end
% Now assign back
scan_parameter_values = ordered_scan_values;
od_imgs = ordered_od_imgs;
end
% Extract quantities
fft_imgs = cell(1, nimgs);
spectral_distribution = cell(1, nimgs);
theta_values = cell(1, nimgs);
radial_spectral_contrast = zeros(1, nimgs);
angular_spectral_weight = zeros(1, nimgs);
N_shots = length(od_imgs);
for k = 1:N_shots
IMG = od_imgs{k};
if ~(max(IMG(:)) > 1)
IMGFFT = NaN(size(IMG));
else
[IMGFFT, IMGPR] = computeFourierTransform(IMG, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization);
end
% Size of original image (in pixels)
[Ny, Nx] = size(IMG);
% Real-space pixel size in micrometers after magnification
dx = pixel_size / magnification;
dy = dx; % assuming square pixels
% Real-space axes
x = ((1:Nx) - ceil(Nx/2)) * dx * 1E6;
y = ((1:Ny) - ceil(Ny/2)) * dy * 1E6;
% Reciprocal space increments (frequency domain, μm⁻¹)
dvx = 1 / (Nx * dx);
dvy = 1 / (Ny * dy);
% Frequency axes
vx = (-floor(Nx/2):ceil(Nx/2)-1) * dvx;
vy = (-floor(Ny/2):ceil(Ny/2)-1) * dvy;
% Wavenumber axes
kx_full = 2 * pi * vx * 1E-6; % μm⁻¹
ky_full = 2 * pi * vy * 1E-6;
% Crop FFT image around center
mid_x = floor(Nx/2);
mid_y = floor(Ny/2);
fft_imgs{k} = IMGFFT(mid_y-zoom_size:mid_y+zoom_size, mid_x-zoom_size:mid_x+zoom_size);
% Crop wavenumber axes to match fft_imgs{k}
kx = kx_full(mid_x - zoom_size : mid_x + zoom_size);
ky = ky_full(mid_y - zoom_size : mid_y + zoom_size);
[theta_vals, S_theta] = computeAngularSpectralDistribution(fft_imgs{k}, kx, ky, k_min, k_max, N_angular_bins, Angular_Threshold, Angular_Sigma, []);
[k_rho_vals, S_k] = computeRadialSpectralDistribution(fft_imgs{k}, kx, ky, theta_min, theta_max, N_radial_bins);
S_k_smoothed = movmean(S_k, Radial_WindowSize); % Compute moving average (use convolution) or use conv for more control
spectral_distribution{k} = S_theta;
theta_values{k} = theta_vals;
radial_spectral_contrast(k) = computeRadialSpectralContrast(k_rho_vals, S_k_smoothed, k_min, k_max);
S_theta_norm = S_theta / max(S_theta); % Normalize to 1
angular_spectral_weight(k) = trapz(theta_vals, S_theta_norm);
end
% 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_rsc = zeros(size(unique_scan_parameter_values));
stderr_rsc = 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 = radial_spectral_contrast(idx == i);
mean_rsc(i) = mean(group_vals, 'omitnan');
stderr_rsc(i) = std(group_vals, 'omitnan') / sqrt(length(group_vals)); % standard error = std / sqrt(N)
end
% Preallocate arrays
mean_asw = zeros(size(unique_scan_parameter_values));
stderr_asw = 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 = angular_spectral_weight(idx == i);
mean_asw(i) = mean(group_vals, 'omitnan');
stderr_asw(i) = std(group_vals, 'omitnan') / sqrt(length(group_vals)); % standard error = std / sqrt(N)
end
% Convert spectral distribution to matrix (N_shots x N_angular_bins)
delta_nkr_all = zeros(N_shots, N_angular_bins);
for k = 1:N_shots
delta_nkr_all(k, :) = spectral_distribution{k};
end
% Group by scan parameter values (e.g., alpha, angle, etc.)
[unique_scan_parameter_values, ~, idx] = unique(scan_parameter_values);
N_params = length(unique_scan_parameter_values);
% Define angular range and conversion
angle_range = 180;
angle_per_bin = angle_range / N_angular_bins;
max_peak_angle = 180;
max_peak_bin = round(max_peak_angle / angle_per_bin);
% Parameters for search
window_size = 10;
angle_threshold = 100;
% Initialize containers for final results
mean_max_g2_values = zeros(1, N_params);
mean_max_g2_angle_values = zeros(1, N_params);
var_max_g2_values = zeros(1, N_params);
var_max_g2_angle_values = zeros(1, N_params);
std_error_g2_values = zeros(1, N_params);
% Also store raw data per group
g2_all_per_group = cell(1, N_params);
angle_all_per_group = cell(1, N_params);
for i = 1:N_params
group_idx = find(idx == i);
group_data = delta_nkr_all(group_idx, :);
N_reps = size(group_data, 1);
g2_values = zeros(1, N_reps);
angle_at_max_g2 = zeros(1, N_reps);
for j = 1:N_reps
profile = group_data(j, :);
% Restrict search to 0–60° for highest peak
restricted_profile = profile(1:max_peak_bin);
[~, peak_idx_rel] = max(restricted_profile);
peak_idx = peak_idx_rel;
peak_angle = (peak_idx - 1) * angle_per_bin;
if peak_angle < angle_threshold
offsets = round(50 / angle_per_bin) : round(70 / angle_per_bin);
else
offsets = -round(70 / angle_per_bin) : -round(50 / angle_per_bin);
end
ref_window = mod((peak_idx - window_size):(peak_idx + window_size) - 1, N_angular_bins) + 1;
ref = profile(ref_window);
correlations = zeros(size(offsets));
angles = zeros(size(offsets));
for k = 1:length(offsets)
shifted_idx = mod(peak_idx + offsets(k) - 1, N_angular_bins) + 1;
sec_window = mod((shifted_idx - window_size):(shifted_idx + window_size) - 1, N_angular_bins) + 1;
sec = profile(sec_window);
num = mean(ref .* sec, 'omitnan');
denom = mean(ref.^2, 'omitnan');
g2 = num / denom;
correlations(k) = g2;
angles(k) = mod((peak_idx - 1 + offsets(k)) * angle_per_bin, angle_range);
end
[max_corr, max_idx] = max(correlations);
g2_values(j) = max_corr;
angle_at_max_g2(j) = angles(max_idx);
end
% Store raw values
g2_all_per_group{i} = g2_values;
angle_all_per_group{i} = angle_at_max_g2;
% Final stats
mean_max_g2_values(i) = mean(g2_values, 'omitnan');
var_max_g2_values(i) = var(g2_values, 0, 'omitnan');
mean_max_g2_angle_values(i)= mean(angle_at_max_g2, 'omitnan');
var_max_g2_angle_values(i) = var(angle_at_max_g2, 0, 'omitnan');
n_i = numel(g2_all_per_group{i}); % Number of repetitions for this param
std_error_g2_values(i) = sqrt(var_max_g2_values(i) / n_i);
end
results.folderPath = folderPath;
results.scan_parameter = scan_parameter;
results.scan_groups = scan_groups;
results.mean_max_g2_values = mean_max_g2_values;
results.std_error_g2_values = std_error_g2_values;
results.mean_max_g2_angle = mean_max_g2_angle_values;
results.radial_spectral_contrast= mean_rsc;
results.angular_spectral_weight = mean_asw;
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 [k_rho_vals, S_radial] = computeRadialSpectralDistribution(IMGFFT, kx, ky, thetamin, thetamax, num_bins)
% IMGFFT : 2D FFT image (fftshifted and cropped)
% kx, ky : 1D physical wavenumber axes [μm⁻¹] matching FFT size
% thetamin : Minimum angle (in radians)
% thetamax : Maximum angle (in radians)
% num_bins : Number of radial bins
[KX, KY] = meshgrid(kx, ky);
K_rho = sqrt(KX.^2 + KY.^2);
Theta = atan2(KY, KX);
if thetamin < thetamax
angle_mask = (Theta >= thetamin) & (Theta <= thetamax);
else
angle_mask = (Theta >= thetamin) | (Theta <= thetamax);
end
power_spectrum = abs(IMGFFT).^2;
r_min = min(K_rho(angle_mask));
r_max = max(K_rho(angle_mask));
r_edges = linspace(r_min, r_max, num_bins + 1);
k_rho_vals = 0.5 * (r_edges(1:end-1) + r_edges(2:end));
S_radial = zeros(1, num_bins);
for i = 1:num_bins
r_low = r_edges(i);
r_high = r_edges(i + 1);
radial_mask = (K_rho >= r_low) & (K_rho < r_high);
full_mask = radial_mask & angle_mask;
S_radial(i) = sum(power_spectrum(full_mask));
end
end
function [theta_vals, S_theta] = computeAngularSpectralDistribution(IMGFFT, kx, ky, k_min, k_max, num_bins, threshold, sigma, windowSize)
% Apply threshold to isolate strong peaks
IMGFFT(IMGFFT < threshold) = 0;
% Create wavenumber meshgrid
[KX, KY] = meshgrid(kx, ky);
Kmag = sqrt(KX.^2 + KY.^2); % radial wavenumber magnitude
Theta = atan2(KY, KX); % range [-pi, pi]
% Restrict to radial band in wavenumber space
radial_mask = (Kmag >= k_min) & (Kmag <= k_max);
% Initialize angular structure factor
S_theta = zeros(1, num_bins);
theta_vals = linspace(0, pi, num_bins); % only 0 to pi due to symmetry
% Loop over angular bins
for i = 1:num_bins
angle_start = (i - 1) * pi / num_bins;
angle_end = i * pi / num_bins;
angle_mask = (Theta >= angle_start) & (Theta < angle_end);
bin_mask = radial_mask & angle_mask;
fft_angle = IMGFFT .* bin_mask;
S_theta(i) = sum(sum(abs(fft_angle).^2));
end
% Optional smoothing
if exist('sigma', 'var') && ~isempty(sigma)
% Gaussian smoothing
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);
% Circular convolution
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);
elseif exist('windowSize', 'var') && ~isempty(windowSize)
% Moving average smoothing
pad = floor(windowSize / 2);
kernel = ones(1, windowSize) / windowSize;
S_theta = conv([S_theta(end - pad + 1:end), S_theta, S_theta(1:pad)], kernel, 'same');
S_theta = S_theta(pad + 1:end - pad);
end
end
function contrast = computeRadialSpectralContrast(k_rho_vals, S_k_smoothed, k_min, k_max)
% Computes the ratio of the peak in S_k_smoothed within [k_min, k_max]
% to the value at (or near) k = 0.
% Ensure inputs are column vectors
k_rho_vals = k_rho_vals(:);
S_k_smoothed = S_k_smoothed(:);
% Step 1: Find index of k ≈ 0
[~, idx_k0] = min(abs(k_rho_vals)); % Closest to zero
S_k0 = S_k_smoothed(idx_k0);
% Step 2: Find indices in specified k-range
in_range = (k_rho_vals >= k_min) & (k_rho_vals <= k_max);
if ~any(in_range)
warning('No values found in the specified k-range. Returning NaN.');
contrast = NaN;
return;
end
% Step 3: Find peak value in the specified k-range
S_k_peak = max(S_k_smoothed(in_range));
% Step 4: Compute contrast
contrast = S_k_peak / S_k0;
end
function ret = getBkgOffsetFromCorners(img, x_fraction, y_fraction)
% image must be a 2D numerical array
[dim1, dim2] = size(img);
s1 = img(1:round(dim1 * y_fraction), 1:round(dim2 * x_fraction));
s2 = img(1:round(dim1 * y_fraction), round(dim2 - dim2 * x_fraction):dim2);
s3 = img(round(dim1 - dim1 * y_fraction):dim1, 1:round(dim2 * x_fraction));
s4 = img(round(dim1 - dim1 * y_fraction):dim1, round(dim2 - dim2 * x_fraction):dim2);
ret = mean([mean(s1(:)), mean(s2(:)), mean(s3(:)), mean(s4(:))]);
end
function ret = subtractBackgroundOffset(img, fraction)
% Remove the background from the image.
% :param dataArray: The image
% :type dataArray: xarray DataArray
% :param x_fraction: The fraction of the pixels used in x axis
% :type x_fraction: float
% :param y_fraction: The fraction of the pixels used in y axis
% :type y_fraction: float
% :return: The image after removing background
% :rtype: xarray DataArray
x_fraction = fraction(1);
y_fraction = fraction(2);
offset = getBkgOffsetFromCorners(img, x_fraction, y_fraction);
ret = img - offset;
end
function ret = cropODImage(img, center, span)
% Crop the image according to the region of interest (ROI).
% :param dataSet: The images
% :type dataSet: xarray DataArray or DataSet
% :param center: The center of region of interest (ROI)
% :type center: tuple
% :param span: The span of region of interest (ROI)
% :type span: tuple
% :return: The cropped images
% :rtype: xarray DataArray or DataSet
x_start = floor(center(1) - span(1) / 2);
x_end = floor(center(1) + span(1) / 2);
y_start = floor(center(2) - span(2) / 2);
y_end = floor(center(2) + span(2) / 2);
ret = img(y_start:y_end, x_start:x_end);
end
function imageOD = calculateODImage(imageAtom, imageBackground, imageDark, mode, exposureTime)
%CALCULATEODIMAGE Calculates the optical density (OD) image for absorption imaging.
%
% imageOD = calculateODImage(imageAtom, imageBackground, imageDark, mode, exposureTime)
%
% Inputs:
% imageAtom - Image with atoms
% imageBackground - Image without atoms
% imageDark - Image without light
% mode - 'LowIntensity' (default) or 'HighIntensity'
% exposureTime - Required only for 'HighIntensity' [in seconds]
%
% Output:
% imageOD - Computed OD image
%
arguments
imageAtom (:,:) {mustBeNumeric}
imageBackground (:,:) {mustBeNumeric}
imageDark (:,:) {mustBeNumeric}
mode char {mustBeMember(mode, {'LowIntensity', 'HighIntensity'})} = 'LowIntensity'
exposureTime double = NaN
end
% Compute numerator and denominator
numerator = imageBackground - imageDark;
denominator = imageAtom - imageDark;
% Avoid division by zero
numerator(numerator == 0) = 1;
denominator(denominator == 0) = 1;
% Calculate OD based on mode
switch mode
case 'LowIntensity'
imageOD = -log(abs(denominator ./ numerator));
case 'HighIntensity'
if isnan(exposureTime)
error('Exposure time must be provided for HighIntensity mode.');
end
imageOD = abs(denominator ./ numerator);
imageOD = -log(imageOD) + (numerator - denominator) ./ (7000 * (exposureTime / 5e-6));
end
end
function [optrefimages] = removefringesInImage(absimages, refimages, bgmask)
% removefringesInImage - Fringe removal and noise reduction from absorption images.
% Creates an optimal reference image for each absorption image in a set as
% a linear combination of reference images, with coefficients chosen to
% minimize the least-squares residuals between each absorption image and
% the optimal reference image. The coefficients are obtained by solving a
% linear set of equations using matrix inverse by LU decomposition.
%
% Application of the algorithm is described in C. F. Ockeloen et al, Improved
% detection of small atom numbers through image processing, arXiv:1007.2136 (2010).
%
% Syntax:
% [optrefimages] = removefringesInImage(absimages,refimages,bgmask);
%
% Required inputs:
% absimages - Absorption image data,
% typically 16 bit grayscale images
% refimages - Raw reference image data
% absimages and refimages are both cell arrays containing
% 2D array data. The number of refimages can differ from the
% number of absimages.
%
% Optional inputs:
% bgmask - Array specifying background region used,
% 1=background, 0=data. Defaults to all ones.
% Outputs:
% optrefimages - Cell array of optimal reference images,
% equal in size to absimages.
%
% Dependencies: none
%
% Authors: Shannon Whitlock, Caspar Ockeloen
% Reference: C. F. Ockeloen, A. F. Tauschinsky, R. J. C. Spreeuw, and
% S. Whitlock, Improved detection of small atom numbers through
% image processing, arXiv:1007.2136
% Email:
% May 2009; Last revision: 11 August 2010
% Process inputs
% Set variables, and flatten absorption and reference images
nimgs = size(absimages,3);
nimgsR = size(refimages,3);
xdim = size(absimages(:,:,1),2);
ydim = size(absimages(:,:,1),1);
R = single(reshape(refimages,xdim*ydim,nimgsR));
A = single(reshape(absimages,xdim*ydim,nimgs));
optrefimages=zeros(size(absimages)); % preallocate
if not(exist('bgmask','var')); bgmask=ones(ydim,xdim); end
k = find(bgmask(:)==1); % Index k specifying background region
% Ensure there are no duplicate reference images
% R=unique(R','rows')'; % comment this line if you run out of memory
% Decompose B = R*R' using singular value or LU decomposition
[L,U,p] = lu(R(k,:)'*R(k,:),'vector'); % LU decomposition
for j=1:nimgs
b=R(k,:)'*A(k,j);
% Obtain coefficients c which minimise least-square residuals
lower.LT = true; upper.UT = true;
c = linsolve(U,linsolve(L,b(p,:),lower),upper);
% Compute optimised reference image
optrefimages(:,:,j)=reshape(R*c,[ydim xdim]);
end
end

576
Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/analyzewithPCA.m
View File

@ -1,576 +0,0 @@
%% Extract Images
clear; close all; clc;
%% ===== D-S Settings =====
groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1_Axis_Camera/in_situ_absorption", ...
"/images/ODT_2_Axis_Camera/in_situ_absorption", "/images/Horizontal_Axis_Camera/in_situ_absorption", ...
"/images/Vertical_Axis_Camera/in_situ_absorption"];
folderPath = "//DyLabNAS/Data/TwoDGas/2025/06/23/";
run = '0300';
folderPath = strcat(folderPath, run);
cam = 5;
angle = 0;
center = [1410, 2030];
span = [200, 200];
fraction = [0.1, 0.1];
pixel_size = 5.86e-6; % in meters
magnification = 23.94;
removeFringes = false;
ImagingMode = 'HighIntensity';
PulseDuration = 5e-6; % in s
% Fourier analysis settings
% Radial Spectral Distribution
theta_min = deg2rad(0);
theta_max = deg2rad(180);
N_radial_bins = 500;
Radial_Sigma = 2;
Radial_WindowSize = 5; % Choose an odd number for a centered moving average
% Angular Spectral Distribution
r_min = 10;
r_max = 20;
N_angular_bins = 180;
Angular_Threshold = 75;
Angular_Sigma = 2;
Angular_WindowSize = 5;
zoom_size = 50; % Zoomed-in region around center
% Plotting and saving
scan_parameter = 'ps_rot_mag_fin_pol_angle';
% scan_parameter = 'rot_mag_field';
savefileName = 'DropletsToStripes';
font = 'Bahnschrift';
if strcmp(savefileName, 'DropletsToStripes')
scan_groups = 0:5:45;
titleString = 'Droplets to Stripes';
elseif strcmp(savefileName, 'StripesToDroplets')
scan_groups = 45:-5:0;
titleString = 'Stripes to Droplets';
end
% Flags
skipNormalization = true;
skipUnshuffling = true;
skipPreprocessing = true;
skipMasking = true;
skipIntensityThresholding = true;
skipBinarization = true;
skipMovieRender = true;
skipSaveFigures = true;
skipSaveOD = true;
%% ===== S-D Settings =====
groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1_Axis_Camera/in_situ_absorption", ...
"/images/ODT_2_Axis_Camera/in_situ_absorption", "/images/Horizontal_Axis_Camera/in_situ_absorption", ...
"/images/Vertical_Axis_Camera/in_situ_absorption"];
folderPath = "//DyLabNAS/Data/TwoDGas/2025/06/24/";
run = '0001';
folderPath = strcat(folderPath, run);
cam = 5;
angle = 0;
center = [1410, 2030];
span = [200, 200];
fraction = [0.1, 0.1];
pixel_size = 5.86e-6; % in meters
magnification = 23.94;
removeFringes = false;
ImagingMode = 'HighIntensity';
PulseDuration = 5e-6; % in s
% Fourier analysis settings
% Radial Spectral Distribution
theta_min = deg2rad(0);
theta_max = deg2rad(180);
N_radial_bins = 500;
Radial_Sigma = 2;
Radial_WindowSize = 5; % Choose an odd number for a centered moving average
% Angular Spectral Distribution
r_min = 10;
r_max = 20;
N_angular_bins = 180;
Angular_Threshold = 75;
Angular_Sigma = 2;
Angular_WindowSize = 5;
zoom_size = 50; % Zoomed-in region around center
% Plotting and saving
scan_parameter = 'ps_rot_mag_fin_pol_angle';
% scan_parameter = 'rot_mag_field';
savefileName = 'StripesToDroplets';
font = 'Bahnschrift';
if strcmp(savefileName, 'DropletsToStripes')
scan_groups = 0:5:45
titleString = 'Droplets to Stripes';
elseif strcmp(savefileName, 'StripesToDroplets')
scan_groups = 45:-5:0;
titleString = 'Stripes to Droplets';
end
% Flags
skipNormalization = true;
skipUnshuffling = false;
skipPreprocessing = true;
skipMasking = true;
skipIntensityThresholding = true;
skipBinarization = true;
skipMovieRender = true;
skipSaveFigures = true;
skipSaveOD = true;
%% ===== Load and compute OD image, rotate and extract ROI for analysis =====
% Get a list of all files in the folder with the desired file name pattern.
filePattern = fullfile(folderPath, '*.h5');
files = dir(filePattern);
refimages = zeros(span(1) + 1, span(2) + 1, length(files));
absimages = zeros(span(1) + 1, span(2) + 1, length(files));
for k = 1 : length(files)
baseFileName = files(k).name;
fullFileName = fullfile(files(k).folder, baseFileName);
fprintf(1, 'Now reading %s\n', fullFileName);
atm_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/atoms")), angle));
bkg_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/background")), angle));
dark_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/dark")), angle));
if (isempty(atm_img) && isa(atm_img, 'double')) || ...
(isempty(bkg_img) && isa(bkg_img, 'double')) || ...
(isempty(dark_img) && isa(dark_img, 'double'))
refimages(:,:,k) = nan(size(refimages(:,:,k))); % fill with NaNs
absimages(:,:,k) = nan(size(absimages(:,:,k)));
else
refimages(:,:,k) = subtractBackgroundOffset(cropODImage(bkg_img, center, span), fraction)';
absimages(:,:,k) = subtractBackgroundOffset(cropODImage(calculateODImage(atm_img, bkg_img, dark_img, ImagingMode, PulseDuration), center, span), fraction)';
end
end
% ===== Fringe removal =====
if removeFringes
optrefimages = removefringesInImage(absimages, refimages);
absimages_fringe_removed = absimages(:, :, :) - optrefimages(:, :, :);
nimgs = size(absimages_fringe_removed,3);
od_imgs = cell(1, nimgs);
for i = 1:nimgs
od_imgs{i} = absimages_fringe_removed(:, :, i);
end
else
nimgs = size(absimages(:, :, :),3);
od_imgs = cell(1, nimgs);
for i = 1:nimgs
od_imgs{i} = absimages(:, :, i);
end
end
% ===== Get rotation angles =====
scan_parameter_values = zeros(1, length(files));
% Get information about the '/globals' group
for k = 1 : length(files)
baseFileName = files(k).name;
fullFileName = fullfile(files(k).folder, baseFileName);
info = h5info(fullFileName, '/globals');
for i = 1:length(info.Attributes)
if strcmp(info.Attributes(i).Name, scan_parameter)
if strcmp(scan_parameter, 'ps_rot_mag_fin_pol_angle')
scan_parameter_values(k) = 180 - info.Attributes(i).Value;
else
scan_parameter_values(k) = info.Attributes(i).Value;
end
end
end
end
% ===== Unshuffle if necessary to do so =====
if ~skipUnshuffling
n_values = length(scan_groups);
n_total = length(scan_parameter_values);
% Infer number of repetitions
n_reps = n_total / n_values;
% Preallocate ordered arrays
ordered_scan_values = zeros(1, n_total);
ordered_od_imgs = cell(1, n_total);
counter = 1;
for rep = 1:n_reps
for val = scan_groups
% Find the next unused match for this val
idx = find(scan_parameter_values == val, 1, 'first');
% Assign and remove from list to avoid duplicates
ordered_scan_values(counter) = scan_parameter_values(idx);
ordered_od_imgs{counter} = od_imgs{idx};
% Mark as used by removing
scan_parameter_values(idx) = NaN; % NaN is safe since original values are 0:5:45
od_imgs{idx} = []; % empty cell so it won't be matched again
counter = counter + 1;
end
end
% Now assign back
scan_parameter_values = ordered_scan_values;
od_imgs = ordered_od_imgs;
end
%% Carry out PCA
numPCs = 5;
% Stack all 600 images into one data matrix [nImages x nPixels]
allImgs3D = cat(3, od_imgs{:});
[Nx, Ny] = size(allImgs3D(:,:,1));
Xall = reshape(allImgs3D, [], numel(od_imgs))'; % [600 x (Nx*Ny)]
% Global PCA
[coeff, score, ~, ~, explained] = pca(Xall);
%% Visualize PC1
% Extract the first principal component vector (eigenimage)
pc1_vector = coeff(:,1);
% Reshape back to original image dimensions
pc1_image = reshape(pc1_vector, Nx, Ny);
% Plot the PC1 image
figure(1); clf; set(gcf, 'Color', 'w', 'Position', [100 100 950 750]);
imagesc(pc1_image);
axis image off;
colormap(Colormaps.coolwarm()); % or use 'jet', 'parula', etc.
colorbar;
title(sprintf('First Principal Component (PC1) Image - Explains %.2f%% Variance', explained(1)));
%% Distribution scatter plot
numGroups = numel(scan_groups);
colors = lines(numGroups);
figure(2); clf; set(gcf, 'Color', 'w', 'Position', [100 100 950 750]); hold on;
for g = 1:numGroups
idx = scan_parameter_values == scan_groups(g);
scatter(repmat(scan_groups(g), sum(idx),1), score(idx,1), 36, colors(g,:), 'filled');
end
xlabel('Control Parameter');
ylabel('PC1 Score');
title('Evolution of PC1 Scores');
grid on;
%% Distribution Histogram plot
numGroups = length(scan_groups);
colors = lines(numGroups);
% Define number of bins globally
numBins = 20;
% Define common bin edges based on global PC1 score range
minScore = min(score(:,1));
maxScore = max(score(:,1));
binEdges = linspace(minScore, maxScore, numBins+1); % +1 because edges are one more than bins
binWidth = binEdges(2) - binEdges(1); % for scaling KDE
figure(3);
clf; set(gcf, 'Color', 'w', 'Position', [100 100 950 750]);
tiledlayout(ceil(numGroups/2), 2, 'TileSpacing', 'compact', 'Padding', 'compact');
for g = 1:numGroups
groupVal = scan_groups(g);
idx = scan_parameter_values == groupVal;
groupPC1 = score(idx,1);
nexttile;
% Plot histogram
histogram(groupPC1, 'Normalization', 'probability', ...
'FaceColor', colors(g,:), 'EdgeColor', 'none', ...
'BinEdges', binEdges);
hold on;
% Compute KDE
[f, xi] = ksdensity(groupPC1, 'NumPoints', 1000);
% Scale KDE to histogram probability scale
f_scaled = f * binWidth;
% Overlay KDE curve
plot(xi, f_scaled, 'k', 'LineWidth', 1.5);
% Vertical line at median
med = median(groupPC1);
yl = ylim;
plot([med med], yl, 'k--', 'LineWidth', 1);
xlabel('PC1 Score');
ylabel('Probability');
title(sprintf('Control Parameter = %d', groupVal));
grid on;
hold off;
end
sgtitle('PC1 Score Distributions');
%% Box plot for PC1 scores by group
groupLabels = cell(size(score,1),1);
for g = 1:numGroups
idx = scan_parameter_values == scan_groups(g);
groupLabels(idx) = {sprintf('%d', scan_groups(g))};
end
figure(4);
clf; set(gcf, 'Color', 'w', 'Position', [100 100 950 750]);
boxplot(score(:,1), groupLabels);
xlabel('Control Parameter');
ylabel('PC1 Score');
title('Evolution of PC1 Scores');
grid on;
%% Mean and SEM plot for PC1 scores
numGroups = length(scan_groups);
meanPC1Scores = zeros(numGroups,1);
semPC1Scores = zeros(numGroups,1);
for g = 1:numGroups
groupVal = scan_groups(g);
idx = scan_parameter_values == groupVal;
groupPC1 = score(idx,1); % PC1 scores for this group
meanPC1Scores(g) = mean(groupPC1);
semPC1Scores(g) = std(groupPC1)/sqrt(sum(idx)); % Standard error of mean
end
% Plot mean ± SEM with error bars
figure(5);
clf; set(gcf, 'Color', 'w', 'Position', [100 100 950 750]);
errorbar(scan_groups, meanPC1Scores, semPC1Scores, 'o-', ...
'LineWidth', 1.5, 'MarkerSize', 8, 'MarkerFaceColor', 'b');
xlabel('Control Parameter');
ylabel('Mean PC1 Score ± SEM');
title('Evolution of PC1 Scores');
grid on;
%% Plot Binder Cumulant
maxOrder = 4; % We only need up to order 4 here
numGroups = length(scan_groups);
kappa4 = NaN(1, numGroups);
for g = 1:numGroups
groupVal = scan_groups(g);
idx = scan_parameter_values == groupVal;
groupPC1 = score(idx, 1);
cumulants = computeCumulants(groupPC1, maxOrder);
kappa4(g) = cumulants(4); % 4th-order cumulant
end
% Plot
figure(6);
clf; set(gcf, 'Color', 'w', 'Position', [100 100 950 750]);
plot(scan_groups, kappa4 * 1E-5, '-o', 'LineWidth', 1.5, 'MarkerFaceColor', 'b');
ylim([-12 12])
xlabel('Control Parameter');
ylabel('\kappa_4 (\times 10^{5})');
grid on;
title('Evolution of Binder Cumulant of PC1 Score');
%% --- ANOVA test ---
p = anova1(score(:,1), groupLabels, 'off');
fprintf('ANOVA p-value for PC1 score differences between groups: %.4e\n', p);
%% Helper Functions
function ret = getBkgOffsetFromCorners(img, x_fraction, y_fraction)
% image must be a 2D numerical array
[dim1, dim2] = size(img);
s1 = img(1:round(dim1 * y_fraction), 1:round(dim2 * x_fraction));
s2 = img(1:round(dim1 * y_fraction), round(dim2 - dim2 * x_fraction):dim2);
s3 = img(round(dim1 - dim1 * y_fraction):dim1, 1:round(dim2 * x_fraction));
s4 = img(round(dim1 - dim1 * y_fraction):dim1, round(dim2 - dim2 * x_fraction):dim2);
ret = mean([mean(s1(:)), mean(s2(:)), mean(s3(:)), mean(s4(:))]);
end
function ret = subtractBackgroundOffset(img, fraction)
% Remove the background from the image.
% :param dataArray: The image
% :type dataArray: xarray DataArray
% :param x_fraction: The fraction of the pixels used in x axis
% :type x_fraction: float
% :param y_fraction: The fraction of the pixels used in y axis
% :type y_fraction: float
% :return: The image after removing background
% :rtype: xarray DataArray
x_fraction = fraction(1);
y_fraction = fraction(2);
offset = getBkgOffsetFromCorners(img, x_fraction, y_fraction);
ret = img - offset;
end
function ret = cropODImage(img, center, span)
% Crop the image according to the region of interest (ROI).
% :param dataSet: The images
% :type dataSet: xarray DataArray or DataSet
% :param center: The center of region of interest (ROI)
% :type center: tuple
% :param span: The span of region of interest (ROI)
% :type span: tuple
% :return: The cropped images
% :rtype: xarray DataArray or DataSet
x_start = floor(center(1) - span(1) / 2);
x_end = floor(center(1) + span(1) / 2);
y_start = floor(center(2) - span(2) / 2);
y_end = floor(center(2) + span(2) / 2);
ret = img(y_start:y_end, x_start:x_end);
end
function imageOD = calculateODImage(imageAtom, imageBackground, imageDark, mode, exposureTime)
%CALCULATEODIMAGE Calculates the optical density (OD) image for absorption imaging.
%
% imageOD = calculateODImage(imageAtom, imageBackground, imageDark, mode, exposureTime)
%
% Inputs:
% imageAtom - Image with atoms
% imageBackground - Image without atoms
% imageDark - Image without light
% mode - 'LowIntensity' (default) or 'HighIntensity'
% exposureTime - Required only for 'HighIntensity' [in seconds]
%
% Output:
% imageOD - Computed OD image
%
arguments
imageAtom (:,:) {mustBeNumeric}
imageBackground (:,:) {mustBeNumeric}
imageDark (:,:) {mustBeNumeric}
mode char {mustBeMember(mode, {'LowIntensity', 'HighIntensity'})} = 'LowIntensity'
exposureTime double = NaN
end
% Compute numerator and denominator
numerator = imageBackground - imageDark;
denominator = imageAtom - imageDark;
% Avoid division by zero
numerator(numerator == 0) = 1;
denominator(denominator == 0) = 1;
% Calculate OD based on mode
switch mode
case 'LowIntensity'
imageOD = -log(abs(denominator ./ numerator));
case 'HighIntensity'
if isnan(exposureTime)
error('Exposure time must be provided for HighIntensity mode.');
end
imageOD = abs(denominator ./ numerator);
imageOD = -log(imageOD) + (numerator - denominator) ./ (7000 * (exposureTime / 5e-6));
end
end
function [optrefimages] = removefringesInImage(absimages, refimages, bgmask)
% removefringesInImage - Fringe removal and noise reduction from absorption images.
% Creates an optimal reference image for each absorption image in a set as
% a linear combination of reference images, with coefficients chosen to
% minimize the least-squares residuals between each absorption image and
% the optimal reference image. The coefficients are obtained by solving a
% linear set of equations using matrix inverse by LU decomposition.
%
% Application of the algorithm is described in C. F. Ockeloen et al, Improved
% detection of small atom numbers through image processing, arXiv:1007.2136 (2010).
%
% Syntax:
% [optrefimages] = removefringesInImage(absimages,refimages,bgmask);
%
% Required inputs:
% absimages - Absorption image data,
% typically 16 bit grayscale images
% refimages - Raw reference image data
% absimages and refimages are both cell arrays containing
% 2D array data. The number of refimages can differ from the
% number of absimages.
%
% Optional inputs:
% bgmask - Array specifying background region used,
% 1=background, 0=data. Defaults to all ones.
% Outputs:
% optrefimages - Cell array of optimal reference images,
% equal in size to absimages.
%
% Dependencies: none
%
% Authors: Shannon Whitlock, Caspar Ockeloen
% Reference: C. F. Ockeloen, A. F. Tauschinsky, R. J. C. Spreeuw, and
% S. Whitlock, Improved detection of small atom numbers through
% image processing, arXiv:1007.2136
% Email:
% May 2009; Last revision: 11 August 2010
% Process inputs
% Set variables, and flatten absorption and reference images
nimgs = size(absimages,3);
nimgsR = size(refimages,3);
xdim = size(absimages(:,:,1),2);
ydim = size(absimages(:,:,1),1);
R = single(reshape(refimages,xdim*ydim,nimgsR));
A = single(reshape(absimages,xdim*ydim,nimgs));
optrefimages=zeros(size(absimages)); % preallocate
if not(exist('bgmask','var')); bgmask=ones(ydim,xdim); end
k = find(bgmask(:)==1); % Index k specifying background region
% Ensure there are no duplicate reference images
% R=unique(R','rows')'; % comment this line if you run out of memory
% Decompose B = R*R' using singular value or LU decomposition
[L,U,p] = lu(R(k,:)'*R(k,:),'vector'); % LU decomposition
for j=1:nimgs
b=R(k,:)'*A(k,j);
% Obtain coefficients c which minimise least-square residuals
lower.LT = true; upper.UT = true;
c = linsolve(U,linsolve(L,b(p,:),lower),upper);
% Compute optimised reference image
optrefimages(:,:,j)=reshape(R*c,[ydim xdim]);
end
end

41
Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/bootstrapCumulants.m
View File

@ -1,41 +0,0 @@
function [cumulants_mean, cumulants_ci, bootstrap_samples] = bootstrapCumulants(x, maxOrder, nBoot)
% bootstrapCumulants - compute bootstrap estimates of cumulants and confidence intervals
%
% Syntax:
% [meanC, ciC, allC] = bootstrapCumulants(x, maxOrder, nBoot)
%
% Inputs:
% x - Data vector (may contain NaNs)
% maxOrder - Max cumulant order (default: 6)
% nBoot - Number of bootstrap samples (default: 1000)
%
% Outputs:
% cumulants_mean - Mean of bootstrap cumulants
% cumulants_ci - 95% confidence intervals [2.5th; 97.5th] percentile
% bootstrap_samples - All bootstrap cumulants (nBoot x maxOrder)
if nargin < 2, maxOrder = 6; end
if nargin < 3, nBoot = 1000; end
x = x(:);
x = x(~isnan(x)); % Remove NaNs
if isempty(x)
cumulants_mean = NaN(1, maxOrder);
cumulants_ci = NaN(2, maxOrder);
bootstrap_samples = NaN(nBoot, maxOrder);
return;
end
N = numel(x);
bootstrap_samples = zeros(nBoot, maxOrder);
for b = 1:nBoot
xb = x(randi(N, [N, 1])); % Resample with replacement
bootstrap_samples(b, :) = computeCumulants(xb, maxOrder);
end
cumulants_mean = mean(bootstrap_samples, 1);
cumulants_ci = prctile(bootstrap_samples, [2.5, 97.5]);
end

39
Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/compareAngularCorrelation.m
View File

@ -1,39 +0,0 @@
%% Track spectral weight across the transition
set(0,'defaulttextInterpreter','latex')
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
format long
font = 'Bahnschrift';
% Load data
Data = load('C:/Users/Karthik/Documents/GitRepositories/Calculations/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/Max_g2_DropletsToStripes.mat', 'unique_scan_parameter_values', 'mean_max_g2_values', 'std_error_g2_values');
dts_scan_parameter_values = Data.unique_scan_parameter_values;
dts_mean_mg2 = Data.mean_max_g2_values;
dts_stderr_mg2 = Data.std_error_g2_values;
Data = load('C:/Users/Karthik/Documents/GitRepositories/Calculations/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/Max_g2_StripesToDroplets.mat', 'unique_scan_parameter_values', 'mean_max_g2_values', 'std_error_g2_values');
std_scan_parameter_values = Data.unique_scan_parameter_values;
std_mean_mg2 = Data.mean_max_g2_values;
std_stderr_mg2 = Data.std_error_g2_values;
figure(1);
set(gcf,'Position',[100 100 950 750])
errorbar(dts_scan_parameter_values, dts_mean_mg2, dts_stderr_mg2, 'o--', ...
'LineWidth', 1.5, 'MarkerSize', 6, 'CapSize', 5, 'DisplayName' , 'Droplets to Stripes');
hold on
errorbar(std_scan_parameter_values, std_mean_mg2, std_stderr_mg2, 'o--', ...
'LineWidth', 1.5, 'MarkerSize', 6, 'CapSize', 5, 'DisplayName', 'Stripes to Droplets');
set(gca, 'FontSize', 14, 'YLim', [0, 1]);
hXLabel = xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex');
hYLabel = ylabel('$\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', 'Interpreter', 'latex');
% hTitle = title('B = 2.42 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
%%

1372
Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/conductSpectralAnalysis.m
View File

File diff suppressed because it is too large Load Diff

545
Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/extractAutocorrelation.m
View File

@ -1,545 +0,0 @@
%% ===== Settings =====
groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1_Axis_Camera/in_situ_absorption", ...
"/images/ODT_2_Axis_Camera/in_situ_absorption", "/images/Horizontal_Axis_Camera/in_situ_absorption", ...
"/images/Vertical_Axis_Camera/in_situ_absorption"];
folderPath = "//DyLabNAS/Data/TwoDGas/2025/06/23/";
run = '0300';
folderPath = strcat(folderPath, run);
cam = 5;
angle = 0;
center = [1410, 2030];
span = [200, 200];
fraction = [0.1, 0.1];
pixel_size = 5.86e-6; % in meters
magnification = 23.94;
removeFringes = false;
ImagingMode = 'HighIntensity';
PulseDuration = 5e-6; % in s
% Fourier analysis settings
% Radial Spectral Distribution
theta_min = deg2rad(0);
theta_max = deg2rad(180);
N_radial_bins = 500;
Radial_Sigma = 2;
Radial_WindowSize = 5; % Choose an odd number for a centered moving average
% Angular Spectral Distribution
r_min = 10;
r_max = 20;
N_angular_bins = 180;
Angular_Threshold = 75;
Angular_Sigma = 2;
Angular_WindowSize = 5;
zoom_size = 50; % Zoomed-in region around center
% Plotting and saving
scan_parameter = 'ps_rot_mag_fin_pol_angle';
% scan_parameter = 'rot_mag_field';
scan_parameter_text = 'Angle = ';
% scan_parameter_text = 'BField = ';
savefileName = 'DropletsToStripes';
font = 'Bahnschrift';
if strcmp(savefileName, 'DropletsToStripes')
scan_groups = 0:5:45;
titleString = 'Droplets to Stripes';
elseif strcmp(savefileName, 'StripesToDroplets')
scan_groups = 45:-5:0;
titleString = 'Stripes to Droplets';
end
% Flags
skipUnshuffling = true;
skipPreprocessing = true;
skipMasking = true;
skipIntensityThresholding = true;
skipBinarization = true;
skipMovieRender = true;
skipSaveFigures = true;
%% ===== Load and compute OD image, rotate and extract ROI for analysis =====
% Get a list of all files in the folder with the desired file name pattern.
filePattern = fullfile(folderPath, '*.h5');
files = dir(filePattern);
refimages = zeros(span(1) + 1, span(2) + 1, length(files));
absimages = zeros(span(1) + 1, span(2) + 1, length(files));
for k = 1 : length(files)
baseFileName = files(k).name;
fullFileName = fullfile(files(k).folder, baseFileName);
fprintf(1, 'Now reading %s\n', fullFileName);
atm_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/atoms")), angle));
bkg_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/background")), angle));
dark_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/dark")), angle));
if (isempty(atm_img) && isa(atm_img, 'double')) || ...
(isempty(bkg_img) && isa(bkg_img, 'double')) || ...
(isempty(dark_img) && isa(dark_img, 'double'))
refimages(:,:,k) = nan(size(refimages(:,:,k))); % fill with NaNs
absimages(:,:,k) = nan(size(absimages(:,:,k)));
else
refimages(:,:,k) = subtractBackgroundOffset(cropODImage(bkg_img, center, span), fraction)';
absimages(:,:,k) = subtractBackgroundOffset(cropODImage(calculateODImage(atm_img, bkg_img, dark_img, ImagingMode, PulseDuration), center, span), fraction)';
end
end
%% ===== Fringe removal =====
if removeFringes
optrefimages = removefringesInImage(absimages, refimages);
absimages_fringe_removed = absimages(:, :, :) - optrefimages(:, :, :);
nimgs = size(absimages_fringe_removed,3);
od_imgs = cell(1, nimgs);
for i = 1:nimgs
od_imgs{i} = absimages_fringe_removed(:, :, i);
end
else
nimgs = size(absimages(:, :, :),3);
od_imgs = cell(1, nimgs);
for i = 1:nimgs
od_imgs{i} = absimages(:, :, i);
end
end
%% ===== Get rotation angles =====
scan_parameter_values = zeros(1, length(files));
% Get information about the '/globals' group
for k = 1 : length(files)
baseFileName = files(k).name;
fullFileName = fullfile(files(k).folder, baseFileName);
info = h5info(fullFileName, '/globals');
for i = 1:length(info.Attributes)
if strcmp(info.Attributes(i).Name, scan_parameter)
if strcmp(scan_parameter, 'ps_rot_mag_fin_pol_angle')
scan_parameter_values(k) = 180 - info.Attributes(i).Value;
else
scan_parameter_values(k) = info.Attributes(i).Value;
end
end
end
end
%% ===== Extract g2 from experiment data =====
fft_imgs = cell(1, nimgs);
spectral_distribution = cell(1, nimgs);
theta_values = cell(1, nimgs);
N_shots = length(od_imgs);
% Compute FFT
for k = 1:N_shots
IMG = od_imgs{k};
[IMGFFT, IMGPR] = computeFourierTransform(IMG, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization);
% Size of original image (in pixels)
[Ny, Nx] = size(IMG);
% Real-space pixel size in micrometers after magnification
dx = pixel_size / magnification;
dy = dx; % assuming square pixels
% Real-space axes
x = ((1:Nx) - ceil(Nx/2)) * dx * 1E6;
y = ((1:Ny) - ceil(Ny/2)) * dy * 1E6;
% Reciprocal space increments (frequency domain, μm⁻¹)
dvx = 1 / (Nx * dx);
dvy = 1 / (Ny * dy);
% Frequency axes
vx = (-floor(Nx/2):ceil(Nx/2)-1) * dvx;
vy = (-floor(Ny/2):ceil(Ny/2)-1) * dvy;
% Wavenumber axes
kx_full = 2 * pi * vx * 1E-6; % μm⁻¹
ky_full = 2 * pi * vy * 1E-6;
% Crop FFT image around center
mid_x = floor(Nx/2);
mid_y = floor(Ny/2);
fft_imgs{k} = IMGFFT(mid_y-zoom_size:mid_y+zoom_size, mid_x-zoom_size:mid_x+zoom_size);
% Crop wavenumber axes to match fft_imgs{k}
kx = kx_full(mid_x - zoom_size : mid_x + zoom_size);
ky = ky_full(mid_y - zoom_size : mid_y + zoom_size);
[theta_values, S_theta] = computeAngularSpectralDistribution(fft_imgs{k}, r_min, r_max, N_angular_bins, Angular_Threshold, Angular_Sigma, []);
spectral_distribution{k} = S_theta;
end
% Create matrix of shape (N_shots x N_angular_bins)
delta_nkr_all = zeros(N_shots, N_angular_bins);
for k = 1:N_shots
delta_nkr_all(k, :) = spectral_distribution{k};
end
% Grouping by scan parameter value (e.g., alpha)
[unique_scan_parameter_values, ~, idx] = unique(scan_parameter_values);
% Number of unique parameter values
N_params = length(unique_scan_parameter_values);
% Preallocate result arrays
g2_all = zeros(N_params, N_angular_bins);
g2_error_all = zeros(N_params, N_angular_bins);
% Compute g2
for i = 1:N_params
group_idx = find(idx == i);
group_data = delta_nkr_all(group_idx, :);
for dtheta = 0:N_angular_bins-1
temp = zeros(length(group_idx), 1);
for j = 1:length(group_idx)
profile = group_data(j, :);
profile_shifted = circshift(profile, -dtheta, 2);
num = mean(profile .* profile_shifted);
denom = mean(profile.^2);
temp(j) = num / denom;
end
g2_all(i, dtheta+1) = mean(temp, 'omitnan');
g2_error_all(i, dtheta+1) = std(temp, 'omitnan') / sqrt(length(group_idx)); % Standard error
end
end
% Number of unique parameter values
nParams = size(g2_all, 1);
% Generate a colormap with enough unique colors
cmap = sky(nParams); % You can also try 'jet', 'turbo', 'hot', etc.
figure(1);
clf;
set(gcf, 'Color', 'w', 'Position',[100 100 950 750])
hold on;
legend_entries = cell(nParams, 1);
for i = 1:nParams
errorbar(theta_values/pi, g2_all(i, :), g2_error_all(i, :), ...
'o', 'Color', cmap(i,:), ...
'MarkerSize', 3, 'MarkerFaceColor', cmap(i,:), ...
'CapSize', 4);
if strcmp(scan_parameter, 'ps_rot_mag_fin_pol_angle')
legend_entries{i} = sprintf('$\\alpha = %g^\\circ$', unique_scan_parameter_values(i));
elseif strcmp(scan_parameter, 'rot_mag_field')
legend_entries{i} = sprintf('B = %.2f G', unique_scan_parameter_values(i));
end
end
ylim([0.0 1.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;
%% 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] = computeAngularSpectralDistribution(IMGFFT, r_min, r_max, num_bins, threshold, sigma, windowSize)
% 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 angular structure factor
S_theta = zeros(1, num_bins);
theta_vals = linspace(0, pi, num_bins);
% Loop through angle bins
for i = 1:num_bins
angle_start = (i-1) * pi / num_bins;
angle_end = i * pi / num_bins;
angle_mask = (Theta >= angle_start & Theta < angle_end);
bin_mask = radial_mask & angle_mask;
fft_angle = IMGFFT .* bin_mask;
S_theta(i) = sum(sum(abs(fft_angle).^2));
end
% Smooth using either Gaussian or moving average
if exist('sigma', 'var') && ~isempty(sigma)
% Gaussian convolution
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);
% Circular convolution
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);
elseif exist('windowSize', 'var') && ~isempty(windowSize)
% Moving average via convolution (circular)
pad = floor(windowSize / 2);
kernel = ones(1, windowSize) / windowSize;
S_theta = conv([S_theta(end-pad+1:end), S_theta, S_theta(1:pad)], kernel, 'same');
S_theta = S_theta(pad+1:end-pad);
end
end
function ret = getBkgOffsetFromCorners(img, x_fraction, y_fraction)
% image must be a 2D numerical array
[dim1, dim2] = size(img);
s1 = img(1:round(dim1 * y_fraction), 1:round(dim2 * x_fraction));
s2 = img(1:round(dim1 * y_fraction), round(dim2 - dim2 * x_fraction):dim2);
s3 = img(round(dim1 - dim1 * y_fraction):dim1, 1:round(dim2 * x_fraction));
s4 = img(round(dim1 - dim1 * y_fraction):dim1, round(dim2 - dim2 * x_fraction):dim2);
ret = mean([mean(s1(:)), mean(s2(:)), mean(s3(:)), mean(s4(:))]);
end
function ret = subtractBackgroundOffset(img, fraction)
% Remove the background from the image.
% :param dataArray: The image
% :type dataArray: xarray DataArray
% :param x_fraction: The fraction of the pixels used in x axis
% :type x_fraction: float
% :param y_fraction: The fraction of the pixels used in y axis
% :type y_fraction: float
% :return: The image after removing background
% :rtype: xarray DataArray
x_fraction = fraction(1);
y_fraction = fraction(2);
offset = getBkgOffsetFromCorners(img, x_fraction, y_fraction);
ret = img - offset;
end
function ret = cropODImage(img, center, span)
% Crop the image according to the region of interest (ROI).
% :param dataSet: The images
% :type dataSet: xarray DataArray or DataSet
% :param center: The center of region of interest (ROI)
% :type center: tuple
% :param span: The span of region of interest (ROI)
% :type span: tuple
% :return: The cropped images
% :rtype: xarray DataArray or DataSet
x_start = floor(center(1) - span(1) / 2);
x_end = floor(center(1) + span(1) / 2);
y_start = floor(center(2) - span(2) / 2);
y_end = floor(center(2) + span(2) / 2);
ret = img(y_start:y_end, x_start:x_end);
end
function imageOD = calculateODImage(imageAtom, imageBackground, imageDark, mode, exposureTime)
%CALCULATEODIMAGE Calculates the optical density (OD) image for absorption imaging.
%
% imageOD = calculateODImage(imageAtom, imageBackground, imageDark, mode, exposureTime)
%
% Inputs:
% imageAtom - Image with atoms
% imageBackground - Image without atoms
% imageDark - Image without light
% mode - 'LowIntensity' (default) or 'HighIntensity'
% exposureTime - Required only for 'HighIntensity' [in seconds]
%
% Output:
% imageOD - Computed OD image
%
arguments
imageAtom (:,:) {mustBeNumeric}
imageBackground (:,:) {mustBeNumeric}
imageDark (:,:) {mustBeNumeric}
mode char {mustBeMember(mode, {'LowIntensity', 'HighIntensity'})} = 'LowIntensity'
exposureTime double = NaN
end
% Compute numerator and denominator
numerator = imageBackground - imageDark;
denominator = imageAtom - imageDark;
% Avoid division by zero
numerator(numerator == 0) = 1;
denominator(denominator == 0) = 1;
% Calculate OD based on mode
switch mode
case 'LowIntensity'
imageOD = -log(abs(denominator ./ numerator));
case 'HighIntensity'
if isnan(exposureTime)
error('Exposure time must be provided for HighIntensity mode.');
end
imageOD = abs(denominator ./ numerator);
imageOD = -log(imageOD) + (numerator - denominator) ./ (7000 * (exposureTime / 5e-6));
end
end
function [optrefimages] = removefringesInImage(absimages, refimages, bgmask)
% removefringesInImage - Fringe removal and noise reduction from absorption images.
% Creates an optimal reference image for each absorption image in a set as
% a linear combination of reference images, with coefficients chosen to
% minimize the least-squares residuals between each absorption image and
% the optimal reference image. The coefficients are obtained by solving a
% linear set of equations using matrix inverse by LU decomposition.
%
% Application of the algorithm is described in C. F. Ockeloen et al, Improved
% detection of small atom numbers through image processing, arXiv:1007.2136 (2010).
%
% Syntax:
% [optrefimages] = removefringesInImage(absimages,refimages,bgmask);
%
% Required inputs:
% absimages - Absorption image data,
% typically 16 bit grayscale images
% refimages - Raw reference image data
% absimages and refimages are both cell arrays containing
% 2D array data. The number of refimages can differ from the
% number of absimages.
%
% Optional inputs:
% bgmask - Array specifying background region used,
% 1=background, 0=data. Defaults to all ones.
% Outputs:
% optrefimages - Cell array of optimal reference images,
% equal in size to absimages.
%
% Dependencies: none
%
% Authors: Shannon Whitlock, Caspar Ockeloen
% Reference: C. F. Ockeloen, A. F. Tauschinsky, R. J. C. Spreeuw, and
% S. Whitlock, Improved detection of small atom numbers through
% image processing, arXiv:1007.2136
% Email:
% May 2009; Last revision: 11 August 2010
% Process inputs
% Set variables, and flatten absorption and reference images
nimgs = size(absimages,3);
nimgsR = size(refimages,3);
xdim = size(absimages(:,:,1),2);
ydim = size(absimages(:,:,1),1);
R = single(reshape(refimages,xdim*ydim,nimgsR));
A = single(reshape(absimages,xdim*ydim,nimgs));
optrefimages=zeros(size(absimages)); % preallocate
if not(exist('bgmask','var')); bgmask=ones(ydim,xdim); end
k = find(bgmask(:)==1); % Index k specifying background region
% Ensure there are no duplicate reference images
% R=unique(R','rows')'; % comment this line if you run out of memory
% Decompose B = R*R' using singular value or LU decomposition
[L,U,p] = lu(R(k,:)'*R(k,:),'vector'); % LU decomposition
for j=1:nimgs
b=R(k,:)'*A(k,j);
% Obtain coefficients c which minimise least-square residuals
lower.LT = true; upper.UT = true;
c = linsolve(U,linsolve(L,b(p,:),lower),upper);
% Compute optimised reference image
optrefimages(:,:,j)=reshape(R*c,[ydim xdim]);
end
end

738
Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/extractCustomCorrelation.m
View File

@ -1,738 +0,0 @@
%% ===== D-S Settings =====
groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1_Axis_Camera/in_situ_absorption", ...
"/images/ODT_2_Axis_Camera/in_situ_absorption", "/images/Horizontal_Axis_Camera/in_situ_absorption", ...
"/images/Vertical_Axis_Camera/in_situ_absorption"];
folderPath = "//DyLabNAS/Data/TwoDGas/2025/07/22/";
run = '0021';
folderPath = strcat(folderPath, run);
cam = 5;
angle = 0;
center = [1410, 2030];
span = [200, 200];
fraction = [0.1, 0.1];
pixel_size = 5.86e-6; % in meters
magnification = 23.94;
removeFringes = false;
ImagingMode = 'HighIntensity';
PulseDuration = 5e-6; % in s
% Fourier analysis settings
% Radial Spectral Distribution
theta_min = deg2rad(0);
theta_max = deg2rad(180);
N_radial_bins = 500;
Radial_Sigma = 2;
Radial_WindowSize = 5; % Choose an odd number for a centered moving average
% Angular Spectral Distribution
r_min = 10;
r_max = 20;
N_angular_bins = 180;
Angular_Threshold = 75;
Angular_Sigma = 2;
Angular_WindowSize = 5;
zoom_size = 50; % Zoomed-in region around center
% Plotting and saving
scan_parameter = 'ps_rot_mag_fin_pol_angle';
% scan_parameter = 'rot_mag_field';
savefileName = 'DropletsToStripes';
font = 'Bahnschrift';
if strcmp(savefileName, 'DropletsToStripes')
scan_groups = 0:1:40;
titleString = 'Droplets to Stripes';
elseif strcmp(savefileName, 'StripesToDroplets')
scan_groups = 40:-1:0;
titleString = 'Stripes to Droplets';
end
% Flags
skipNormalization = true;
skipUnshuffling = false;
skipPreprocessing = true;
skipMasking = true;
skipIntensityThresholding = true;
skipBinarization = true;
skipMovieRender = true;
skipSaveFigures = false;
skipSaveOD = true;
%% ===== Load and compute OD image, rotate and extract ROI for analysis =====
% Get a list of all files in the folder with the desired file name pattern.
filePattern = fullfile(folderPath, '*.h5');
files = dir(filePattern);
refimages = zeros(span(1) + 1, span(2) + 1, length(files));
absimages = zeros(span(1) + 1, span(2) + 1, length(files));
for k = 1 : length(files)
baseFileName = files(k).name;
fullFileName = fullfile(files(k).folder, baseFileName);
fprintf(1, 'Now reading %s\n', fullFileName);
atm_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/atoms")), angle));
bkg_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/background")), angle));
dark_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/dark")), angle));
if (isempty(atm_img) && isa(atm_img, 'double')) || ...
(isempty(bkg_img) && isa(bkg_img, 'double')) || ...
(isempty(dark_img) && isa(dark_img, 'double'))
refimages(:,:,k) = nan(size(refimages(:,:,k))); % fill with NaNs
absimages(:,:,k) = nan(size(absimages(:,:,k)));
else
refimages(:,:,k) = subtractBackgroundOffset(cropODImage(bkg_img, center, span), fraction)';
absimages(:,:,k) = subtractBackgroundOffset(cropODImage(calculateODImage(atm_img, bkg_img, dark_img, ImagingMode, PulseDuration), center, span), fraction)';
end
end
% ===== Fringe removal =====
if removeFringes
optrefimages = removefringesInImage(absimages, refimages);
absimages_fringe_removed = absimages(:, :, :) - optrefimages(:, :, :);
nimgs = size(absimages_fringe_removed,3);
od_imgs = cell(1, nimgs);
for i = 1:nimgs
od_imgs{i} = absimages_fringe_removed(:, :, i);
end
else
nimgs = size(absimages(:, :, :),3);
od_imgs = cell(1, nimgs);
for i = 1:nimgs
od_imgs{i} = absimages(:, :, i);
end
end
%% ===== Get rotation angles =====
scan_parameter_values = zeros(1, length(files));
% Get information about the '/globals' group
for k = 1 : length(files)
baseFileName = files(k).name;
fullFileName = fullfile(files(k).folder, baseFileName);
info = h5info(fullFileName, '/globals');
for i = 1:length(info.Attributes)
if strcmp(info.Attributes(i).Name, scan_parameter)
if strcmp(scan_parameter, 'ps_rot_mag_fin_pol_angle')
scan_parameter_values(k) = 180 - info.Attributes(i).Value;
else
scan_parameter_values(k) = info.Attributes(i).Value;
end
end
end
end
%% ===== Correlation of a single (highest) peak with a possible peak between 50-70 degrees from experiment data =====
fft_imgs = cell(1, nimgs);
spectral_distribution = cell(1, nimgs);
theta_values = cell(1, nimgs);
N_shots = length(od_imgs);
% Compute FFT for all images
for k = 1:N_shots
IMG = od_imgs{k};
[IMGFFT, IMGPR] = computeFourierTransform(IMG, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization);
[Ny, Nx] = size(IMG);
dx = pixel_size / magnification;
dy = dx; % assuming square pixels
x = ((1:Nx) - ceil(Nx/2)) * dx * 1E6;
y = ((1:Ny) - ceil(Ny/2)) * dy * 1E6;
dvx = 1 / (Nx * dx);
dvy = 1 / (Ny * dy);
vx = (-floor(Nx/2):ceil(Nx/2)-1) * dvx;
vy = (-floor(Ny/2):ceil(Ny/2)-1) * dvy;
kx_full = 2 * pi * vx * 1E-6;
ky_full = 2 * pi * vy * 1E-6;
mid_x = floor(Nx/2);
mid_y = floor(Ny/2);
fft_imgs{k} = IMGFFT(mid_y-zoom_size:mid_y+zoom_size, mid_x-zoom_size:mid_x+zoom_size);
kx = kx_full(mid_x - zoom_size : mid_x + zoom_size);
ky = ky_full(mid_y - zoom_size : mid_y + zoom_size);
[theta_vals, S_theta] = computeAngularSpectralDistribution(fft_imgs{k}, r_min, r_max, N_angular_bins, Angular_Threshold, Angular_Sigma, []);
spectral_distribution{k} = S_theta;
theta_values{k} = theta_vals;
end
% Convert spectral distribution to matrix (N_shots x N_angular_bins)
delta_nkr_all = zeros(N_shots, N_angular_bins);
for k = 1:N_shots
delta_nkr_all(k, :) = spectral_distribution{k};
end
% Group by scan parameter values (e.g., alpha, angle, etc.)
[unique_scan_parameter_values, ~, idx] = unique(scan_parameter_values);
N_params = length(unique_scan_parameter_values);
% Define angular range and conversion
angle_range = 180;
angle_per_bin = angle_range / N_angular_bins;
max_peak_angle = 180;
max_peak_bin = round(max_peak_angle / angle_per_bin);
% Parameters for search
window_size = 10;
angle_threshold = 100;
% Initialize containers for final results
mean_max_g2_values = zeros(1, N_params);
skew_max_g2_angle_values = zeros(1, N_params);
var_max_g2_values = zeros(1, N_params);
fourth_order_cumulant_max_g2_angle_values= zeros(1, N_params);
fifth_order_cumulant_max_g2_angle_values = zeros(1, N_params);
sixth_order_cumulant_max_g2_angle_values = zeros(1, N_params);
% Also store raw data per group
max_g2_all_per_group = cell(1, N_params);
std_error_g2_values = zeros(1, N_params);
for i = 1:N_params
group_idx = find(idx == i);
group_data = delta_nkr_all(group_idx, :);
N_reps = size(group_data, 1);
g2_values = zeros(1, N_reps);
for j = 1:N_reps
profile = group_data(j, :);
% Restrict search to 0–60° for highest peak
restricted_profile = profile(1:max_peak_bin);
[~, peak_idx_rel] = max(restricted_profile);
peak_idx = peak_idx_rel;
peak_angle = (peak_idx - 1) * angle_per_bin;
if peak_angle < angle_threshold
offsets = round(50 / angle_per_bin) : round(70 / angle_per_bin);
else
offsets = -round(70 / angle_per_bin) : -round(50 / angle_per_bin);
end
ref_window = mod((peak_idx - window_size):(peak_idx + window_size) - 1, N_angular_bins) + 1;
ref = profile(ref_window);
correlations = zeros(size(offsets));
for k = 1:length(offsets)
shifted_idx = mod(peak_idx + offsets(k) - 1, N_angular_bins) + 1;
sec_window = mod((shifted_idx - window_size):(shifted_idx + window_size) - 1, N_angular_bins) + 1;
sec = profile(sec_window);
num = mean(ref .* sec);
denom = mean(ref.^2);
g2 = num / denom;
correlations(k) = g2;
end
[max_corr, max_idx] = max(correlations);
g2_values(j) = max_corr;
end
% Store raw values
max_g2_all_per_group{i} = g2_values;
% Compute cumulants
kappa = computeCumulants(g2_values(:), 6);
% Final stats
mean_max_g2_values(i) = kappa(1);
var_max_g2_values(i) = kappa(2);
N_eff = sum(~isnan(g2_values));
std_error_g2_values(i) = sqrt(kappa(2)) / sqrt(N_eff);
skew_max_g2_angle_values(i) = kappa(3);
fourth_order_cumulant_max_g2_angle_values(i)= kappa(4);
fifth_order_cumulant_max_g2_angle_values(i) = kappa(5);
sixth_order_cumulant_max_g2_angle_values(i) = kappa(6);
end
%% Plot PDF of order parameter
if ~skipSaveFigures
% Define folder for saving images
saveFolder = [savefileName '_SavedFigures'];
if ~exist(saveFolder, 'dir')
mkdir(saveFolder);
end
end
figure(2); % one persistent figure
set(gcf, 'Color', 'w', 'Position', [100 100 950 750])
for val = scan_groups
% Find the index i that matches this scan parameter value
i = find(unique_scan_parameter_values == val, 1);
% Skip if not found (sanity check)
if isempty(i)
continue;
end
g2_vals = max_g2_all_per_group{i};
g2_vals = g2_vals(~isnan(g2_vals));
if isempty(g2_vals)
continue;
end
% KDE estimation
[f, xi] = ksdensity(g2_vals, 'NumPoints', 200);
clf;
histogram(g2_vals, 'Normalization', 'pdf', ...
'NumBins', 10, ...
'FaceAlpha', 0.3, ...
'EdgeColor', 'none', ...
'FaceColor', [0.3 0.5 0.8]);
hold on;
plot(xi, f, 'LineWidth', 2, 'Color', [0 0.2 0.6]);
set(gca, 'FontSize', 16);
title(sprintf('%s: \\boldmath$\\alpha = %.1f^{\\circ}$', titleString, val), ...
'FontSize', 16, 'Interpreter', 'latex');
xlabel('$\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', 'Interpreter', 'latex', 'FontSize', 14);
ylabel('PDF', 'FontSize', 14);
xlim([0.0, 1.5]);
grid on;
drawnow;
% ==== Save Figure ====
if ~skipSaveFigures
% Create a filename for each averaged plot
fileNamePNG = fullfile(saveFolder, sprintf('max_g2_analysis_param_%03d.png', val));
% Save current figure as PNG with high resolution
print(gcf, fileNamePNG, '-dpng', '-r300'); % 300 dpi for high quality
else
pause(0.5)
end
end
%% Plot all cumulants
figure(3)
set(gcf, 'Color', 'w', 'Position', [100 100 950 750])
scan_vals = unique_scan_parameter_values;
% Define font style for consistency
axis_fontsize = 14;
label_fontsize = 16;
title_fontsize = 16;
% 1. Mean with error bars
subplot(3,2,1);
errorbar(scan_vals, mean_max_g2_values, std_error_g2_values, 'o-', ...
'LineWidth', 1.5, 'MarkerSize', 6);
title('Mean of $\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', ...
'Interpreter', 'latex', 'FontSize', title_fontsize);
xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', ...
'FontSize', label_fontsize);
ylabel('$\kappa_1$', 'Interpreter', 'latex', ...
'FontSize', label_fontsize);
set(gca, 'FontSize', axis_fontsize);
grid on;
% 2. Variance
subplot(3,2,2);
plot(scan_vals, var_max_g2_values, 's-', 'LineWidth', 1.5, 'MarkerSize', 6);
title('Variance of $\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', ...
'Interpreter', 'latex', 'FontSize', title_fontsize);
xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', ...
'FontSize', label_fontsize);
ylabel('$\kappa_2$', 'Interpreter', 'latex', ...
'FontSize', label_fontsize);
set(gca, 'FontSize', axis_fontsize);
grid on;
% 3. Skewness
subplot(3,2,3);
plot(scan_vals, skew_max_g2_angle_values, 'd-', 'LineWidth', 1.5, 'MarkerSize', 6);
title('Skewness of $\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', ...
'Interpreter', 'latex', 'FontSize', title_fontsize);
xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', ...
'FontSize', label_fontsize);
ylabel('$\kappa_3$', 'Interpreter', 'latex', ...
'FontSize', label_fontsize);
set(gca, 'FontSize', axis_fontsize);
grid on;
% 4. Binder Cumulant
subplot(3,2,4);
plot(scan_vals, fourth_order_cumulant_max_g2_angle_values, '^-', 'LineWidth', 1.5, 'MarkerSize', 6);
title('Binder Cumulant of $\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', ...
'Interpreter', 'latex', 'FontSize', title_fontsize);
xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', ...
'FontSize', label_fontsize);
ylabel('$\kappa_4$', 'Interpreter', 'latex', ...
'FontSize', label_fontsize);
set(gca, 'FontSize', axis_fontsize);
grid on;
% 5. 5th-order cumulant
subplot(3,2,5);
plot(scan_vals, fifth_order_cumulant_max_g2_angle_values, 'v-', 'LineWidth', 1.5, 'MarkerSize', 6);
title('Fifth-order cumulant of $\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', ...
'Interpreter', 'latex', 'FontSize', title_fontsize);
xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', ...
'FontSize', label_fontsize);
ylabel('$\kappa_5$', 'Interpreter', 'latex', ...
'FontSize', label_fontsize);
set(gca, 'FontSize', axis_fontsize);
grid on;
% 6. 6th-order cumulant
subplot(3,2,6);
plot(scan_vals, sixth_order_cumulant_max_g2_angle_values, '>-', 'LineWidth', 1.5, 'MarkerSize', 6);
title('Sixth-order cumulant of $\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', ...
'Interpreter', 'latex', 'FontSize', title_fontsize);
xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', ...
'FontSize', label_fontsize);
ylabel('$\kappa_6$', 'Interpreter', 'latex', ...
'FontSize', label_fontsize);
set(gca, 'FontSize', axis_fontsize);
grid on;
% Super title
sgtitle(sprintf('Cumulants of Peak Offset Angular Correlation - %s', titleString), ...
'FontWeight', 'bold', 'FontSize', 16, 'Interpreter', 'latex');
%% ── Mean ± Std vs. scan parameter ──────────────────────────────────────
% Plot mean ± SEM
figure(1);
set(gcf, 'Color', 'w', 'Position',[100 100 950 750])
set(gca, 'FontSize', 14); % For tick labels only
errorbar(unique_scan_parameter_values, ... % x-axis
mean_max_g2_values, ... % y-axis (mean)
std_error_g2_values, ... % ± SEM
'--o', 'LineWidth', 1.8, 'MarkerSize', 6 );
set(gca, 'FontSize', 14, 'YLim', [0, 1]);
hXLabel = xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex');
hYLabel = ylabel('$\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', 'Interpreter', 'latex');
hTitle = title(titleString, 'Interpreter', 'tex');
% set([hXLabel, hYLabel], 'FontName', font);
set([hXLabel, hYLabel], 'FontSize', 14);
set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold');
grid on;
% Define folder for saving images
saveFolder = [savefileName '_SavedFigures'];
if ~exist(saveFolder, 'dir')
mkdir(saveFolder);
end
save([saveFolder savefileName '.mat'], 'unique_scan_parameter_values', 'mean_max_g2_values', 'std_error_g2_values');
%% 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] = computeAngularSpectralDistribution(IMGFFT, r_min, r_max, num_bins, threshold, sigma, windowSize)
% 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 angular structure factor
S_theta = zeros(1, num_bins);
theta_vals = linspace(0, pi, num_bins);
% Loop through angle bins
for i = 1:num_bins
angle_start = (i-1) * pi / num_bins;
angle_end = i * pi / num_bins;
angle_mask = (Theta >= angle_start & Theta < angle_end);
bin_mask = radial_mask & angle_mask;
fft_angle = IMGFFT .* bin_mask;
S_theta(i) = sum(sum(abs(fft_angle).^2));
end
% Smooth using either Gaussian or moving average
if exist('sigma', 'var') && ~isempty(sigma)
% Gaussian convolution
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);
% Circular convolution
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);
elseif exist('windowSize', 'var') && ~isempty(windowSize)
% Moving average via convolution (circular)
pad = floor(windowSize / 2);
kernel = ones(1, windowSize) / windowSize;
S_theta = conv([S_theta(end-pad+1:end), S_theta, S_theta(1:pad)], kernel, 'same');
S_theta = S_theta(pad+1:end-pad);
end
end
function ret = getBkgOffsetFromCorners(img, x_fraction, y_fraction)
% image must be a 2D numerical array
[dim1, dim2] = size(img);
s1 = img(1:round(dim1 * y_fraction), 1:round(dim2 * x_fraction));
s2 = img(1:round(dim1 * y_fraction), round(dim2 - dim2 * x_fraction):dim2);
s3 = img(round(dim1 - dim1 * y_fraction):dim1, 1:round(dim2 * x_fraction));
s4 = img(round(dim1 - dim1 * y_fraction):dim1, round(dim2 - dim2 * x_fraction):dim2);
ret = mean([mean(s1(:)), mean(s2(:)), mean(s3(:)), mean(s4(:))]);
end
function ret = subtractBackgroundOffset(img, fraction)
% Remove the background from the image.
% :param dataArray: The image
% :type dataArray: xarray DataArray
% :param x_fraction: The fraction of the pixels used in x axis
% :type x_fraction: float
% :param y_fraction: The fraction of the pixels used in y axis
% :type y_fraction: float
% :return: The image after removing background
% :rtype: xarray DataArray
x_fraction = fraction(1);
y_fraction = fraction(2);
offset = getBkgOffsetFromCorners(img, x_fraction, y_fraction);
ret = img - offset;
end
function ret = cropODImage(img, center, span)
% Crop the image according to the region of interest (ROI).
% :param dataSet: The images
% :type dataSet: xarray DataArray or DataSet
% :param center: The center of region of interest (ROI)
% :type center: tuple
% :param span: The span of region of interest (ROI)
% :type span: tuple
% :return: The cropped images
% :rtype: xarray DataArray or DataSet
x_start = floor(center(1) - span(1) / 2);
x_end = floor(center(1) + span(1) / 2);
y_start = floor(center(2) - span(2) / 2);
y_end = floor(center(2) + span(2) / 2);
ret = img(y_start:y_end, x_start:x_end);
end
function imageOD = calculateODImage(imageAtom, imageBackground, imageDark, mode, exposureTime)
%CALCULATEODIMAGE Calculates the optical density (OD) image for absorption imaging.
%
% imageOD = calculateODImage(imageAtom, imageBackground, imageDark, mode, exposureTime)
%
% Inputs:
% imageAtom - Image with atoms
% imageBackground - Image without atoms
% imageDark - Image without light
% mode - 'LowIntensity' (default) or 'HighIntensity'
% exposureTime - Required only for 'HighIntensity' [in seconds]
%
% Output:
% imageOD - Computed OD image
%
arguments
imageAtom (:,:) {mustBeNumeric}
imageBackground (:,:) {mustBeNumeric}
imageDark (:,:) {mustBeNumeric}
mode char {mustBeMember(mode, {'LowIntensity', 'HighIntensity'})} = 'LowIntensity'
exposureTime double = NaN
end
% Compute numerator and denominator
numerator = imageBackground - imageDark;
denominator = imageAtom - imageDark;
% Avoid division by zero
numerator(numerator == 0) = 1;
denominator(denominator == 0) = 1;
% Calculate OD based on mode
switch mode
case 'LowIntensity'
imageOD = -log(abs(denominator ./ numerator));
case 'HighIntensity'
if isnan(exposureTime)
error('Exposure time must be provided for HighIntensity mode.');
end
imageOD = abs(denominator ./ numerator);
imageOD = -log(imageOD) + (numerator - denominator) ./ (7000 * (exposureTime / 5e-6));
end
end
function [optrefimages] = removefringesInImage(absimages, refimages, bgmask)
% removefringesInImage - Fringe removal and noise reduction from absorption images.
% Creates an optimal reference image for each absorption image in a set as
% a linear combination of reference images, with coefficients chosen to
% minimize the least-squares residuals between each absorption image and
% the optimal reference image. The coefficients are obtained by solving a
% linear set of equations using matrix inverse by LU decomposition.
%
% Application of the algorithm is described in C. F. Ockeloen et al, Improved
% detection of small atom numbers through image processing, arXiv:1007.2136 (2010).
%
% Syntax:
% [optrefimages] = removefringesInImage(absimages,refimages,bgmask);
%
% Required inputs:
% absimages - Absorption image data,
% typically 16 bit grayscale images
% refimages - Raw reference image data
% absimages and refimages are both cell arrays containing
% 2D array data. The number of refimages can differ from the
% number of absimages.
%
% Optional inputs:
% bgmask - Array specifying background region used,
% 1=background, 0=data. Defaults to all ones.
% Outputs:
% optrefimages - Cell array of optimal reference images,
% equal in size to absimages.
%
% Dependencies: none
%
% Authors: Shannon Whitlock, Caspar Ockeloen
% Reference: C. F. Ockeloen, A. F. Tauschinsky, R. J. C. Spreeuw, and
% S. Whitlock, Improved detection of small atom numbers through
% image processing, arXiv:1007.2136
% Email:
% May 2009; Last revision: 11 August 2010
% Process inputs
% Set variables, and flatten absorption and reference images
nimgs = size(absimages,3);
nimgsR = size(refimages,3);
xdim = size(absimages(:,:,1),2);
ydim = size(absimages(:,:,1),1);
R = single(reshape(refimages,xdim*ydim,nimgsR));
A = single(reshape(absimages,xdim*ydim,nimgs));
optrefimages=zeros(size(absimages)); % preallocate
if not(exist('bgmask','var')); bgmask=ones(ydim,xdim); end
k = find(bgmask(:)==1); % Index k specifying background region
% Ensure there are no duplicate reference images
% R=unique(R','rows')'; % comment this line if you run out of memory
% Decompose B = R*R' using singular value or LU decomposition
[L,U,p] = lu(R(k,:)'*R(k,:),'vector'); % LU decomposition
for j=1:nimgs
b=R(k,:)'*A(k,j);
% Obtain coefficients c which minimise least-square residuals
lower.LT = true; upper.UT = true;
c = linsolve(U,linsolve(L,b(p,:),lower),upper);
% Compute optimised reference image
optrefimages(:,:,j)=reshape(R*c,[ydim xdim]);
end
end

158
Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/extractQuantities.m
View File

@ -1,158 +0,0 @@
%% Parameters
% === Define folders and settings ===
baseFolder = '//DyLabNAS/Data/TwoDGas/2025/04/';
dates = ["01", "02"]; % Example: three folders
runs = {
["0059", "0060", "0061"],
["0007", "0008", "0009", "0010", "0011"]
};
options.scan_parameter = 'rot_mag_fin_pol_angle';
options.scan_groups = 0:10:50;
options.cam = 5;
% Image cropping and alignment
options.angle = 0;
options.center = [1285, 2100];
options.span = [200, 200];
options.fraction = [0.1, 0.1];
% Imaging and calibration parameters
options.pixel_size = 5.86e-6; % in meters
options.magnification = 23.94;
options.removeFringes = false;
options.ImagingMode = 'HighIntensity';
options.PulseDuration = 5e-6;
% Fourier analysis: Radial
options.theta_min = deg2rad(0);
options.theta_max = deg2rad(180);
options.N_radial_bins = 500;
options.Radial_Sigma = 2;
options.Radial_WindowSize = 5; % Must be odd
% Fourier analysis: Angular
options.r_min = 10;
options.r_max = 20;
options.k_min = 1.2771; % in μm⁻¹
options.k_max = 2.5541; % in μm⁻¹
options.N_angular_bins = 180;
options.Angular_Threshold = 75;
options.Angular_Sigma = 2;
options.Angular_WindowSize = 5;
% Optional visualization / zooming
options.zoom_size = 50;
% Optional flags or settings struct
options.skipUnshuffling = false;
options.skipPreprocessing = true;
options.skipMasking = true;
options.skipIntensityThresholding = true;
options.skipBinarization = true;
% === Loop through folders and collect results ===
results_all = [];
assert(length(dates) == length(runs), ...
'Each entry in `dates` must correspond to a cell in `runs`.');
for i = 1:length(dates)
currentDate = dates(i);
currentRuns = runs{i};
for j = 1:length(currentRuns)
runID = currentRuns(j);
folderPath = fullfile(baseFolder, currentDate, runID);
if ~endsWith(folderPath, filesep)
options.folderPath = [char(folderPath) filesep];
else
options.folderPath = char(folderPath);
end
try
% Unpack options struct into name-value pairs
args = [fieldnames(options), struct2cell(options)]';
args = args(:)';
results = analyzeFolder(args{:});
results_all = [results_all; results];
catch ME
warning("Error processing %s/%s: %s", currentDate, runID, ME.message);
end
end
end
%% Plotting heatmap of mean_max_g2_values
N_x = length(options.scan_groups);
N_y = length(results_all);
BFields = [2.35, 2.15, 2.0, 1.85, 1.7, 1.55, 1.4, 1.35];
% Preallocate
g2_matrix = zeros(N_y, N_x);
for i = 1:N_y
for j = 1:N_x
g2_matrix(i, j) = results_all(i).mean_max_g2_values(j);
end
end
% Plot heatmap
font = 'Bahnschrift';
figure(1)
clf
set(gcf,'Position',[50 50 950 750])
imagesc(options.scan_groups, BFields, g2_matrix);
colormap(sky);
clim([0, 1])
set(gca, 'FontSize', 14, 'YDir', 'normal');
hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex');
hYLabel = ylabel('BField (G)', 'Interpreter', 'tex');
hTitle = title('$\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', 'Interpreter', 'latex');
set([hXLabel, hYLabel], 'FontSize', 14)
set(hTitle, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title
colorbar;
%% Heat map of radial spectral contrast
N_x = length(options.scan_groups);
N_y = length(results_all);
BFields = [2.35, 2.15, 2.0, 1.85, 1.7, 1.55, 1.4, 1.35];
% Preallocate
radial_spectral_contrast_matrix = zeros(N_y, N_x);
for i = 1:N_y
for j = 1:N_x
radial_spectral_contrast_matrix(i, j) = results_all(i).radial_spectral_contrast(j);
end
end
% Plot heatmap
font = 'Bahnschrift';
figure(3)
clf
set(gcf,'Position',[50 50 950 750])
imagesc(options.scan_groups, BFields, radial_spectral_contrast_matrix);
colormap(sky);
clim([0 0.008])
set(gca, 'FontSize', 14, 'YDir', 'normal');
hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex');
hYLabel = ylabel('BField (G)', 'Interpreter', 'tex');
hTitle = title('Radial Spectral Contrast');
set([hXLabel, hYLabel], 'FontSize', 14)
set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title
colorbar;

416
Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/plotImages.m
View File

@ -1,416 +0,0 @@
%% Parameters
groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1_Axis_Camera/in_situ_absorption", ...
"/images/ODT_2_Axis_Camera/in_situ_absorption", "/images/Horizontal_Axis_Camera/in_situ_absorption", ...
"/images/Vertical_Axis_Camera/in_situ_absorption"];
folderPath = "//DyLabNAS/Data/TwoDGas/2025/07/16/";
run = '0002';
folderPath = strcat(folderPath, run);
cam = 5;
angle = 0;
center = [1430, 2025];
span = [200, 200];
fraction = [0.1, 0.1];
pixel_size = 5.86e-6; % in meters
magnification = 23.94;
removeFringes = false;
ImagingMode = 'HighIntensity';
PulseDuration = 5e-6;
% Plotting and saving
scan_parameter = 'evap_rot_mag_field';
scan_groups = 0:10:50;
savefileName = 'Droplets';
font = 'Bahnschrift';
% Flags
skipUnshuffling = true;
%% ===== Load and compute OD image, rotate and extract ROI for analysis =====
% Get a list of all files in the folder with the desired file name pattern.
filePattern = fullfile(folderPath, '*.h5');
files = dir(filePattern);
refimages = zeros(span(1) + 1, span(2) + 1, length(files));
absimages = zeros(span(1) + 1, span(2) + 1, length(files));
for k = 1 : length(files)
baseFileName = files(k).name;
fullFileName = fullfile(files(k).folder, baseFileName);
fprintf(1, 'Now reading %s\n', fullFileName);
atm_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/atoms")), angle));
bkg_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/background")), angle));
dark_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/dark")), angle));
refimages(:,:,k) = subtractBackgroundOffset(cropODImage(bkg_img, center, span), fraction)';
absimages(:,:,k) = subtractBackgroundOffset(cropODImage(calculateODImage(atm_img, bkg_img, dark_img, ImagingMode, PulseDuration), center, span), fraction)';
end
% ===== Fringe removal =====
if removeFringes
optrefimages = removefringesInImage(absimages, refimages);
absimages_fringe_removed = absimages(:, :, :) - optrefimages(:, :, :);
nimgs = size(absimages_fringe_removed,3);
od_imgs = cell(1, nimgs);
for i = 1:nimgs
od_imgs{i} = absimages_fringe_removed(:, :, i);
end
else
nimgs = size(absimages(:, :, :),3);
od_imgs = cell(1, nimgs);
for i = 1:nimgs
od_imgs{i} = absimages(:, :, i);
end
end
% ===== Get rotation angles =====
scan_parameter_values = zeros(1, length(files));
% Get information about the '/globals' group
for k = 1 : length(files)
baseFileName = files(k).name;
fullFileName = fullfile(files(k).folder, baseFileName);
info = h5info(fullFileName, '/globals');
for i = 1:length(info.Attributes)
if strcmp(info.Attributes(i).Name, scan_parameter)
if strcmp(scan_parameter, 'rot_mag_fin_pol_angle')
scan_parameter_values(k) = 180 - info.Attributes(i).Value;
else
scan_parameter_values(k) = info.Attributes(i).Value;
end
end
end
end
% ===== Unshuffle if necessary to do so =====
if ~skipUnshuffling
n_values = length(scan_groups);
n_total = length(scan_parameter_values);
% Infer number of repetitions
n_reps = n_total / n_values;
% Preallocate ordered arrays
ordered_scan_values = zeros(1, n_total);
ordered_od_imgs = cell(1, n_total);
counter = 1;
for rep = 1:n_reps
for val = scan_groups
% Find the next unused match for this val
idx = find(scan_parameter_values == val, 1, 'first');
% Assign and remove from list to avoid duplicates
ordered_scan_values(counter) = scan_parameter_values(idx);
ordered_od_imgs{counter} = od_imgs{idx};
% Mark as used by removing
scan_parameter_values(idx) = NaN; % NaN is safe since original values are 0:5:45
od_imgs{idx} = []; % empty cell so it won't be matched again
counter = counter + 1;
end
end
% Now assign back
scan_parameter_values = ordered_scan_values;
od_imgs = ordered_od_imgs;
end
%% Display Images
figure(1)
clf
set(gcf,'Position',[50 50 950 750])
% Get image size in pixels
[Ny, Nx] = size(od_imgs{1});
% Define pixel size and magnification (if not already defined earlier)
dx = pixel_size / magnification; % e.g. in meters
dy = dx; % assuming square pixels
% Define x and y axes in μm (centered at image center)
x = ((1:Nx) - ceil(Nx/2)) * dx * 1E6; % micrometers
y = ((1:Ny) - ceil(Ny/2)) * dy * 1E6;
% Display the cropped image
for k = 1 : length(od_imgs)
imagesc(x, y, od_imgs{k});
hold on;
% Convert pixel grid to µm (already done: x and y axes)
% Draw ↘ diagonal (top-left to bottom-right)
drawODOverlays(x(1), y(1), x(end), y(end));
% Draw ↙ diagonal (top-right to bottom-left)
drawODOverlays(x(end), y(1), x(1), y(end));
hold off;
axis equal tight;
colormap(Colormaps.inferno());
set(gca, 'FontSize', 14, 'YDir', 'normal');
if strcmp(scan_parameter, 'rot_mag_fin_pol_angle')
text(0.975, 0.975, [num2str(scan_parameter_values(k), '%.1f^\\circ')], ...
'Color', 'white', 'FontWeight', 'bold', 'FontSize', 24, ...
'Interpreter', 'tex', 'Units', 'normalized', ...
'HorizontalAlignment', 'right', 'VerticalAlignment', 'top');
else
text(0.975, 0.975, [num2str(scan_parameter_values(k), '%.2f'), ' G'], ...
'Color', 'white', 'FontWeight', 'bold', 'FontSize', 24, ...
'Interpreter', 'tex', 'Units', 'normalized', ...
'HorizontalAlignment', 'right', 'VerticalAlignment', 'top');
end
colorbarHandle = colorbar;
ylabel(colorbarHandle, 'Optical Density', 'Rotation', -90, 'FontSize', 14, 'FontName', font);
xlabel('x (\mum)', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font);
ylabel('y (\mum)', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font);
title('OD Image', 'FontSize', 16, 'FontWeight', 'bold', 'Interpreter', 'tex', 'FontName', font);
drawnow;
pause(0.5);
end
%% Helper Functions
function ret = getBkgOffsetFromCorners(img, x_fraction, y_fraction)
% image must be a 2D numerical array
[dim1, dim2] = size(img);
s1 = img(1:round(dim1 * y_fraction), 1:round(dim2 * x_fraction));
s2 = img(1:round(dim1 * y_fraction), round(dim2 - dim2 * x_fraction):dim2);
s3 = img(round(dim1 - dim1 * y_fraction):dim1, 1:round(dim2 * x_fraction));
s4 = img(round(dim1 - dim1 * y_fraction):dim1, round(dim2 - dim2 * x_fraction):dim2);
ret = mean([mean(s1(:)), mean(s2(:)), mean(s3(:)), mean(s4(:))]);
end
function ret = subtractBackgroundOffset(img, fraction)
% Remove the background from the image.
% :param dataArray: The image
% :type dataArray: xarray DataArray
% :param x_fraction: The fraction of the pixels used in x axis
% :type x_fraction: float
% :param y_fraction: The fraction of the pixels used in y axis
% :type y_fraction: float
% :return: The image after removing background
% :rtype: xarray DataArray
x_fraction = fraction(1);
y_fraction = fraction(2);
offset = getBkgOffsetFromCorners(img, x_fraction, y_fraction);
ret = img - offset;
end
function ret = cropODImage(img, center, span)
% Crop the image according to the region of interest (ROI).
% :param dataSet: The images
% :type dataSet: xarray DataArray or DataSet
% :param center: The center of region of interest (ROI)
% :type center: tuple
% :param span: The span of region of interest (ROI)
% :type span: tuple
% :return: The cropped images
% :rtype: xarray DataArray or DataSet
x_start = floor(center(1) - span(1) / 2);
x_end = floor(center(1) + span(1) / 2);
y_start = floor(center(2) - span(2) / 2);
y_end = floor(center(2) + span(2) / 2);
ret = img(y_start:y_end, x_start:x_end);
end
function imageOD = calculateODImage(imageAtom, imageBackground, imageDark, mode, exposureTime)
%CALCULATEODIMAGE Calculates the optical density (OD) image for absorption imaging.
%
% imageOD = calculateODImage(imageAtom, imageBackground, imageDark, mode, exposureTime)
%
% Inputs:
% imageAtom - Image with atoms
% imageBackground - Image without atoms
% imageDark - Image without light
% mode - 'LowIntensity' (default) or 'HighIntensity'
% exposureTime - Required only for 'HighIntensity' [in seconds]
%
% Output:
% imageOD - Computed OD image
%
arguments
imageAtom (:,:) {mustBeNumeric}
imageBackground (:,:) {mustBeNumeric}
imageDark (:,:) {mustBeNumeric}
mode char {mustBeMember(mode, {'LowIntensity', 'HighIntensity'})} = 'LowIntensity'
exposureTime double = NaN
end
% Compute numerator and denominator
numerator = imageBackground - imageDark;
denominator = imageAtom - imageDark;
% Avoid division by zero
numerator(numerator == 0) = 1;
denominator(denominator == 0) = 1;
% Calculate OD based on mode
switch mode
case 'LowIntensity'
imageOD = -log(abs(denominator ./ numerator));
case 'HighIntensity'
if isnan(exposureTime)
error('Exposure time must be provided for HighIntensity mode.');
end
imageOD = abs(denominator ./ numerator);
imageOD = -log(imageOD) + (numerator - denominator) ./ (7000 * (exposureTime / 5e-6));
end
end
function drawODOverlays(x1, y1, x2, y2)
% Parameters
tick_spacing = 10; % µm between ticks
tick_length = 2; % µm tick mark length
line_color = [0.5 0.5 0.5];
tick_color = [0.5 0.5 0.5];
font_size = 10;
% Vector from start to end
dx = x2 - x1;
dy = y2 - y1;
L = sqrt(dx^2 + dy^2);
% Unit direction vector along diagonal
ux = dx / L;
uy = dy / L;
% Perpendicular unit vector for ticks
perp_ux = -uy;
perp_uy = ux;
% Midpoint (center)
xc = (x1 + x2) / 2;
yc = (y1 + y2) / 2;
% Number of positive and negative ticks
n_ticks = floor(L / (2 * tick_spacing));
% Draw main diagonal line
plot([x1 x2], [y1 y2], '--', 'Color', line_color, 'LineWidth', 1.2);
for i = -n_ticks:n_ticks
d = i * tick_spacing;
xt = xc + d * ux;
yt = yc + d * uy;
% Tick line endpoints
xt1 = xt - 0.5 * tick_length * perp_ux;
yt1 = yt - 0.5 * tick_length * perp_uy;
xt2 = xt + 0.5 * tick_length * perp_ux;
yt2 = yt + 0.5 * tick_length * perp_uy;
% Draw tick
plot([xt1 xt2], [yt1 yt2], '--', 'Color', tick_color, 'LineWidth', 1);
% Label: centered at tick, offset slightly along diagonal
if d ~= 0
text(xt, yt, sprintf('%+d', d), ...
'Color', tick_color, ...
'FontSize', font_size, ...
'HorizontalAlignment', 'center', ...
'VerticalAlignment', 'bottom', ...
'Rotation', atan2d(dy, dx));
end
end
end
function [optrefimages] = removefringesInImage(absimages, refimages, bgmask)
% removefringesInImage - Fringe removal and noise reduction from absorption images.
% Creates an optimal reference image for each absorption image in a set as
% a linear combination of reference images, with coefficients chosen to
% minimize the least-squares residuals between each absorption image and
% the optimal reference image. The coefficients are obtained by solving a
% linear set of equations using matrix inverse by LU decomposition.
%
% Application of the algorithm is described in C. F. Ockeloen et al, Improved
% detection of small atom numbers through image processing, arXiv:1007.2136 (2010).
%
% Syntax:
% [optrefimages] = removefringesInImage(absimages,refimages,bgmask);
%
% Required inputs:
% absimages - Absorption image data,
% typically 16 bit grayscale images
% refimages - Raw reference image data
% absimages and refimages are both cell arrays containing
% 2D array data. The number of refimages can differ from the
% number of absimages.
%
% Optional inputs:
% bgmask - Array specifying background region used,
% 1=background, 0=data. Defaults to all ones.
% Outputs:
% optrefimages - Cell array of optimal reference images,
% equal in size to absimages.
%
% Dependencies: none
%
% Authors: Shannon Whitlock, Caspar Ockeloen
% Reference: C. F. Ockeloen, A. F. Tauschinsky, R. J. C. Spreeuw, and
% S. Whitlock, Improved detection of small atom numbers through
% image processing, arXiv:1007.2136
% Email:
% May 2009; Last revision: 11 August 2010
% Process inputs
% Set variables, and flatten absorption and reference images
nimgs = size(absimages,3);
nimgsR = size(refimages,3);
xdim = size(absimages(:,:,1),2);
ydim = size(absimages(:,:,1),1);
R = single(reshape(refimages,xdim*ydim,nimgsR));
A = single(reshape(absimages,xdim*ydim,nimgs));
optrefimages=zeros(size(absimages)); % preallocate
if not(exist('bgmask','var')); bgmask=ones(ydim,xdim); end
k = find(bgmask(:)==1); % Index k specifying background region
% Ensure there are no duplicate reference images
% R=unique(R','rows')'; % comment this line if you run out of memory
% Decompose B = R*R' using singular value or LU decomposition
[L,U,p] = lu(R(k,:)'*R(k,:),'vector'); % LU decomposition
for j=1:nimgs
b=R(k,:)'*A(k,j);
% Obtain coefficients c which minimise least-square residuals
lower.LT = true; upper.UT = true;
c = linsolve(U,linsolve(L,b(p,:),lower),upper);
% Compute optimised reference image
optrefimages(:,:,j)=reshape(R*c,[ydim xdim]);
end
end

767
Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/plotPhaseDiagram.m
View File

@ -1,767 +0,0 @@
% === Parameters ===
baseFolder = '//DyLabNAS/Data/TwoDGas/2025/04/';
dates = ["01", "02"];
runs = {
["0059", "0060", "0061"],
["0007", "0008", "0009", "0010", "0011"]
};
scan_groups = 0:10:50;
scan_parameter = 'rot_mag_fin_pol_angle';
cam = 5;
% Image cropping and alignment
angle = 0;
center = [1285, 2100];
span = [200, 200];
fraction = [0.1, 0.1];
% Imaging and calibration parameters
pixel_size = 5.86e-6; % in meters
magnification = 23.94;
removeFringes = false;
ImagingMode = 'LowIntensity';
PulseDuration = 5e-6;
% Optional visualization / zooming
options.zoom_size = 50;
% Optional flags or settings struct
skipUnshuffling = false;
skipPreprocessing = true;
skipMasking = true;
skipIntensityThresholding = true;
skipBinarization = true;
groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1_Axis_Camera/in_situ_absorption", ...
"/images/ODT_2_Axis_Camera/in_situ_absorption", "/images/Horizontal_Axis_Camera/in_situ_absorption", ...
"/images/Vertical_Axis_Camera/in_situ_absorption"];
%%
allData = {}; % now a growing list of structs per B field
dataCounter = 1;
for i = 1:length(dates)
dateStr = dates(i);
runList = runs{i};
for j = 1:length(runList)
folderPath = fullfile(baseFolder, dateStr, runList{j});
filePattern = fullfile(folderPath, '*.h5');
files = dir(filePattern);
refimages = zeros(span(1) + 1, span(2) + 1, length(files));
absimages = zeros(span(1) + 1, span(2) + 1, length(files));
for k = 1 : length(files)
baseFileName = files(k).name;
fullFileName = fullfile(files(k).folder, baseFileName);
fprintf(1, 'Now reading %s\n', fullFileName);
atm_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/atoms")), angle));
bkg_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/background")), angle));
dark_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/dark")), angle));
refimages(:,:,k) = subtractBackgroundOffset(cropODImage(bkg_img, center, span), fraction)';
absimages(:,:,k) = subtractBackgroundOffset(cropODImage(calculateODImage(atm_img, bkg_img, dark_img, ImagingMode, PulseDuration), center, span), fraction)';
end
% ===== Fringe removal =====
if removeFringes
optrefimages = removefringesInImage(absimages, refimages);
absimages_fringe_removed = absimages(:, :, :) - optrefimages(:, :, :);
nimgs = size(absimages_fringe_removed,3);
od_imgs = cell(1, nimgs);
for k = 1:nimgs
od_imgs{k} = absimages_fringe_removed(:, :, k);
end
else
nimgs = size(absimages(:, :, :),3);
od_imgs = cell(1, nimgs);
for k = 1:nimgs
od_imgs{k} = absimages(:, :, k);
end
end
%% ===== Get rotation angles =====
scan_parameter_values = zeros(1, length(files));
% Get information about the '/globals' group
for k = 1 : length(files)
baseFileName = files(k).name;
fullFileName = fullfile(files(k).folder, baseFileName);
info = h5info(fullFileName, '/globals');
for i = 1:length(info.Attributes)
if strcmp(info.Attributes(i).Name, scan_parameter)
if strcmp(scan_parameter, 'rot_mag_fin_pol_angle')
scan_parameter_values(k) = 180 - info.Attributes(i).Value;
else
scan_parameter_values(k) = info.Attributes(i).Value;
end
end
if strcmp(info.Attributes(i).Name, "rot_mag_field")
B = info.Attributes(i).Value;
end
end
end
% ===== Unshuffle if necessary to do so =====
if ~skipUnshuffling
n_values = length(scan_groups);
n_total = length(scan_parameter_values);
% Infer number of repetitions
n_reps = n_total / n_values;
% Preallocate ordered arrays
ordered_scan_values = zeros(1, n_total);
ordered_od_imgs = cell(1, n_total);
counter = 1;
for rep = 1:n_reps
for val = scan_groups
% Find the next unused match for this val
idx = find(scan_parameter_values == val, 1, 'first');
% Assign and remove from list to avoid duplicates
ordered_scan_values(counter) = scan_parameter_values(idx);
ordered_od_imgs{counter} = od_imgs{idx};
% Mark as used by removing
scan_parameter_values(idx) = NaN; % NaN is safe since original values are 0:5:45
od_imgs{idx} = []; % empty cell so it won't be matched again
counter = counter + 1;
end
end
% Now assign back
scan_parameter_values = ordered_scan_values;
od_imgs = ordered_od_imgs;
end
% === Reshape ===
od_imgs_reshaped = reshape(od_imgs, [length(scan_groups), n_reps]);
% === Store ===
allData{dataCounter} = struct(...
'B', B, ...
'theta_vals', scan_groups, ...
'od_imgs', od_imgs_reshaped ...
);
dataCounter = dataCounter + 1;
end
end
%% === % Plot PD - 1st rep of each θ per B-field ===
[theta_vals, ~, idx] = unique(scan_parameter_values);
nB = numel(allData);
nTheta = numel(theta_vals);
% Select every 2nd B-field index
idxToPlot = 1:2:nB; % indices 1, 3, 5, ...
% Update number of B-fields to plot
nB_new = numel(idxToPlot);
figure(101); clf;
% Make the figure wider to fit the colorbar comfortably
set(gcf, 'Position', [100, 100, 1300, 800]);
% Create tiled layout with some right padding to reserve space for colorbar
t = tiledlayout(nB_new, nTheta, 'TileSpacing', 'compact', 'Padding', 'compact');
font = 'Bahnschrift';
allAxes = gobjects(nB_new, nTheta);
for new_i = 1:nB_new
i = idxToPlot(new_i); % original index in allData
data = allData{i};
for j = 1:nTheta
ax = nexttile((new_i-1)*nTheta + j);
allAxes(new_i,j) = ax;
od = data(j).od_imgs;
imagesc(od, 'Parent', ax);
set(ax, 'YDir', 'normal');
axis(ax, 'image');
ax.XTick = [];
ax.YTick = [];
colormap(ax, Colormaps.inferno());
end
end
% Use colorbar associated with the last image tile
cb = colorbar('Location', 'eastoutside');
cb.Layout.Tile = 'east'; % Attach it to the layout edge
cb.FontName = font;
cb.FontSize = 18;
cb.Label.FontSize = 20;
cb.Label.Rotation = 90;
cb.Label.VerticalAlignment = 'bottom';
cb.Label.HorizontalAlignment = 'center';
cb.Direction = 'normal'; % Ensure ticks go bottom-to-top
% Add x and y tick labels along bottom and left
% Use bottom row for θ ticks
for j = 1:nTheta
ax = allAxes(end, j);
ax.XTick = size(od,2)/2;
ax.XTickLabel = sprintf('%d°', theta_vals(j));
ax.XTickLabelRotation = 0;
ax.FontName = font;
ax.FontSize = 20;
end
% Use first column for B ticks (only the plotted subset)
for new_i = 1:nB_new
i = idxToPlot(new_i);
ax = allAxes(new_i, 1);
ax.YTick = size(od,1)/2;
ax.YTickLabel = sprintf('%.2f G', allData{i}(1).B);
ax.FontName = font;
ax.FontSize = 20;
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 [k_rho_vals, S_radial] = computeRadialSpectralDistribution(IMGFFT, kx, ky, thetamin, thetamax, num_bins)
% IMGFFT : 2D FFT image (fftshifted and cropped)
% kx, ky : 1D physical wavenumber axes [μm⁻¹] matching FFT size
% thetamin : Minimum angle (in radians)
% thetamax : Maximum angle (in radians)
% num_bins : Number of radial bins
[KX, KY] = meshgrid(kx, ky);
K_rho = sqrt(KX.^2 + KY.^2);
Theta = atan2(KY, KX);
if thetamin < thetamax
angle_mask = (Theta >= thetamin) & (Theta <= thetamax);
else
angle_mask = (Theta >= thetamin) | (Theta <= thetamax);
end
power_spectrum = abs(IMGFFT).^2;
r_min = min(K_rho(angle_mask));
r_max = max(K_rho(angle_mask));
r_edges = linspace(r_min, r_max, num_bins + 1);
k_rho_vals = 0.5 * (r_edges(1:end-1) + r_edges(2:end));
S_radial = zeros(1, num_bins);
for i = 1:num_bins
r_low = r_edges(i);
r_high = r_edges(i + 1);
radial_mask = (K_rho >= r_low) & (K_rho < r_high);
full_mask = radial_mask & angle_mask;
S_radial(i) = sum(power_spectrum(full_mask));
end
end
function [theta_vals, S_theta] = computeAngularSpectralDistribution(IMGFFT, r_min, r_max, num_bins, threshold, sigma, windowSize)
% 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 angular structure factor
S_theta = zeros(1, num_bins);
theta_vals = linspace(0, pi, num_bins);
% Loop through angle bins
for i = 1:num_bins
angle_start = (i-1) * pi / num_bins;
angle_end = i * pi / num_bins;
angle_mask = (Theta >= angle_start & Theta < angle_end);
bin_mask = radial_mask & angle_mask;
fft_angle = IMGFFT .* bin_mask;
S_theta(i) = sum(sum(abs(fft_angle).^2));
end
% Smooth using either Gaussian or moving average
if exist('sigma', 'var') && ~isempty(sigma)
% Gaussian convolution
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);
% Circular convolution
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);
elseif exist('windowSize', 'var') && ~isempty(windowSize)
% Moving average via convolution (circular)
pad = floor(windowSize / 2);
kernel = ones(1, windowSize) / windowSize;
S_theta = conv([S_theta(end-pad+1:end), S_theta, S_theta(1:pad)], kernel, 'same');
S_theta = S_theta(pad+1:end-pad);
end
end
function contrast = computeRadialSpectralContrast(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);
s1 = img(1:round(dim1 * y_fraction), 1:round(dim2 * x_fraction));
s2 = img(1:round(dim1 * y_fraction), round(dim2 - dim2 * x_fraction):dim2);
s3 = img(round(dim1 - dim1 * y_fraction):dim1, 1:round(dim2 * x_fraction));
s4 = img(round(dim1 - dim1 * y_fraction):dim1, round(dim2 - dim2 * x_fraction):dim2);
ret = mean([mean(s1(:)), mean(s2(:)), mean(s3(:)), mean(s4(:))]);
end
function ret = subtractBackgroundOffset(img, fraction)
% Remove the background from the image.
% :param dataArray: The image
% :type dataArray: xarray DataArray
% :param x_fraction: The fraction of the pixels used in x axis
% :type x_fraction: float
% :param y_fraction: The fraction of the pixels used in y axis
% :type y_fraction: float
% :return: The image after removing background
% :rtype: xarray DataArray
x_fraction = fraction(1);
y_fraction = fraction(2);
offset = getBkgOffsetFromCorners(img, x_fraction, y_fraction);
ret = img - offset;
end
function ret = cropODImage(img, center, span)
% Crop the image according to the region of interest (ROI).
% :param dataSet: The images
% :type dataSet: xarray DataArray or DataSet
% :param center: The center of region of interest (ROI)
% :type center: tuple
% :param span: The span of region of interest (ROI)
% :type span: tuple
% :return: The cropped images
% :rtype: xarray DataArray or DataSet
x_start = floor(center(1) - span(1) / 2);
x_end = floor(center(1) + span(1) / 2);
y_start = floor(center(2) - span(2) / 2);
y_end = floor(center(2) + span(2) / 2);
ret = img(y_start:y_end, x_start:x_end);
end
function imageOD = calculateODImage(imageAtom, imageBackground, imageDark, mode, exposureTime)
%CALCULATEODIMAGE Calculates the optical density (OD) image for absorption imaging.
%
% imageOD = calculateODImage(imageAtom, imageBackground, imageDark, mode, exposureTime)
%
% Inputs:
% imageAtom - Image with atoms
% imageBackground - Image without atoms
% imageDark - Image without light
% mode - 'LowIntensity' (default) or 'HighIntensity'
% exposureTime - Required only for 'HighIntensity' [in seconds]
%
% Output:
% imageOD - Computed OD image
%
arguments
imageAtom (:,:) {mustBeNumeric}
imageBackground (:,:) {mustBeNumeric}
imageDark (:,:) {mustBeNumeric}
mode char {mustBeMember(mode, {'LowIntensity', 'HighIntensity'})} = 'LowIntensity'
exposureTime double = NaN
end
% Compute numerator and denominator
numerator = imageBackground - imageDark;
denominator = imageAtom - imageDark;
% Avoid division by zero
numerator(numerator == 0) = 1;
denominator(denominator == 0) = 1;
% Calculate OD based on mode
switch mode
case 'LowIntensity'
imageOD = -log(abs(denominator ./ numerator));
case 'HighIntensity'
if isnan(exposureTime)
error('Exposure time must be provided for HighIntensity mode.');
end
imageOD = abs(denominator ./ numerator);
imageOD = -log(imageOD) + (numerator - denominator) ./ (7000 * (exposureTime / 5e-6));
end
end
function drawODOverlays(x1, y1, x2, y2)
% Parameters
tick_spacing = 10; % µm between ticks
tick_length = 2; % µm tick mark length
line_color = [0.5 0.5 0.5];
tick_color = [0.5 0.5 0.5];
font_size = 10;
% Vector from start to end
dx = x2 - x1;
dy = y2 - y1;
L = sqrt(dx^2 + dy^2);
% Unit direction vector along diagonal
ux = dx / L;
uy = dy / L;
% Perpendicular unit vector for ticks
perp_ux = -uy;
perp_uy = ux;
% Midpoint (center)
xc = (x1 + x2) / 2;
yc = (y1 + y2) / 2;
% Number of positive and negative ticks
n_ticks = floor(L / (2 * tick_spacing));
% Draw main diagonal line
plot([x1 x2], [y1 y2], '--', 'Color', line_color, 'LineWidth', 1.2);
for i = -n_ticks:n_ticks
d = i * tick_spacing;
xt = xc + d * ux;
yt = yc + d * uy;
% Tick line endpoints
xt1 = xt - 0.5 * tick_length * perp_ux;
yt1 = yt - 0.5 * tick_length * perp_uy;
xt2 = xt + 0.5 * tick_length * perp_ux;
yt2 = yt + 0.5 * tick_length * perp_uy;
% Draw tick
plot([xt1 xt2], [yt1 yt2], '--', 'Color', tick_color, 'LineWidth', 1);
% Label: centered at tick, offset slightly along diagonal
if d ~= 0
text(xt, yt, sprintf('%+d', d), ...
'Color', tick_color, ...
'FontSize', font_size, ...
'HorizontalAlignment', 'center', ...
'VerticalAlignment', 'bottom', ...
'Rotation', atan2d(dy, dx));
end
end
end
function drawPSOverlays(kx, ky, r_min, r_max)
% drawFFTOverlays - Draw overlays on existing FFT plot:
% - Radial lines every 30°
% - Annular highlight with white (upper half) and gray (lower half) circles between r_min and r_max
% - Horizontal white bands at ky=0 in annulus region
% - Scale ticks and labels every 1 μm⁻¹ along each radial line
%
% Inputs:
% kx, ky - reciprocal space vectors (μm⁻¹)
% r_min - inner annulus radius offset index (integer)
% r_max - outer annulus radius offset index (integer)
%
% Example:
% hold on;
% drawFFTOverlays(kx, ky, 10, 30);
hold on
% === Overlay Radial Lines + Scales ===
[kx_grid, ky_grid] = meshgrid(kx, ky);
[~, kr_grid] = cart2pol(kx_grid, ky_grid); % kr_grid in μm⁻¹
max_kx = max(kx);
max_ky = max(ky);
for angle = 0 : pi/6 : pi
x_line = [0, max_kx] * cos(angle);
y_line = [0, max_ky] * sin(angle);
% Plot radial lines
plot(x_line, y_line, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.2);
plot(x_line, -y_line, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.2);
% Draw scale ticks along positive radial line
drawTicksAlongLine(0, 0, x_line(2), y_line(2));
% Draw scale ticks along negative radial line (reflect y)
drawTicksAlongLine(0, 0, x_line(2), -y_line(2));
end
% === Overlay Annular Highlight: White (r_min to r_max), Gray elsewhere ===
theta_full = linspace(0, 2*pi, 500);
center_x = ceil(size(kr_grid, 2) / 2);
center_y = ceil(size(kr_grid, 1) / 2);
k_min = kr_grid(center_y, center_x + r_min);
k_max = kr_grid(center_y, center_x + r_max);
% Upper half: white dashed circles
x1_upper = k_min * cos(theta_full(theta_full <= pi));
y1_upper = k_min * sin(theta_full(theta_full <= pi));
x2_upper = k_max * cos(theta_full(theta_full <= pi));
y2_upper = k_max * sin(theta_full(theta_full <= pi));
plot(x1_upper, y1_upper, 'k--', 'LineWidth', 1.2);
plot(x2_upper, y2_upper, 'k--', 'LineWidth', 1.2);
% Lower half: gray dashed circles
x1_lower = k_min * cos(theta_full(theta_full > pi));
y1_lower = k_min * sin(theta_full(theta_full > pi));
x2_lower = k_max * cos(theta_full(theta_full > pi));
y2_lower = k_max * sin(theta_full(theta_full > pi));
plot(x1_lower, y1_lower, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.0);
plot(x2_lower, y2_lower, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.0);
% === Highlight horizontal band across k_y = 0 ===
x_vals = kx;
xW1 = x_vals((x_vals >= -k_max) & (x_vals < -k_min));
xW2 = x_vals((x_vals > k_min) & (x_vals <= k_max));
plot(xW1, zeros(size(xW1)), 'k--', 'LineWidth', 1.2);
plot(xW2, zeros(size(xW2)), 'k--', 'LineWidth', 1.2);
hold off
% --- Nested helper function to draw ticks along a radial line ---
function drawTicksAlongLine(x_start, y_start, x_end, y_end)
% Tick parameters
tick_spacing = 1; % spacing between ticks in μm⁻¹
tick_length = 0.05 * sqrt((x_end - x_start)^2 + (y_end - y_start)^2); % relative tick length
line_color = [0.5 0.5 0.5];
tick_color = [0.5 0.5 0.5];
font_size = 8;
% Vector along the line
dx = x_end - x_start;
dy = y_end - y_start;
L = sqrt(dx^2 + dy^2);
ux = dx / L;
uy = dy / L;
% Perpendicular vector for ticks
perp_ux = -uy;
perp_uy = ux;
% Number of ticks (from 0 up to max length)
n_ticks = floor(L / tick_spacing);
for i = 1:n_ticks
% Position of tick along the line
xt = x_start + i * tick_spacing * ux;
yt = y_start + i * tick_spacing * uy;
% Tick endpoints
xt1 = xt - 0.5 * tick_length * perp_ux;
yt1 = yt - 0.5 * tick_length * perp_uy;
xt2 = xt + 0.5 * tick_length * perp_ux;
yt2 = yt + 0.5 * tick_length * perp_uy;
% Draw tick
plot([xt1 xt2], [yt1 yt2], '-', 'Color', tick_color, 'LineWidth', 1);
% Label with distance (integer)
text(xt, yt, sprintf('%d', i), ...
'Color', tick_color, ...
'FontSize', font_size, ...
'HorizontalAlignment', 'center', ...
'VerticalAlignment', 'bottom', ...
'Rotation', atan2d(dy, dx));
end
end
end
function [optrefimages] = removefringesInImage(absimages, refimages, bgmask)
% removefringesInImage - Fringe removal and noise reduction from absorption images.
% Creates an optimal reference image for each absorption image in a set as
% a linear combination of reference images, with coefficients chosen to
% minimize the least-squares residuals between each absorption image and
% the optimal reference image. The coefficients are obtained by solving a
% linear set of equations using matrix inverse by LU decomposition.
%
% Application of the algorithm is described in C. F. Ockeloen et al, Improved
% detection of small atom numbers through image processing, arXiv:1007.2136 (2010).
%
% Syntax:
% [optrefimages] = removefringesInImage(absimages,refimages,bgmask);
%
% Required inputs:
% absimages - Absorption image data,
% typically 16 bit grayscale images
% refimages - Raw reference image data
% absimages and refimages are both cell arrays containing
% 2D array data. The number of refimages can differ from the
% number of absimages.
%
% Optional inputs:
% bgmask - Array specifying background region used,
% 1=background, 0=data. Defaults to all ones.
% Outputs:
% optrefimages - Cell array of optimal reference images,
% equal in size to absimages.
%
% Dependencies: none
%
% Authors: Shannon Whitlock, Caspar Ockeloen
% Reference: C. F. Ockeloen, A. F. Tauschinsky, R. J. C. Spreeuw, and
% S. Whitlock, Improved detection of small atom numbers through
% image processing, arXiv:1007.2136
% Email:
% May 2009; Last revision: 11 August 2010
% Process inputs
% Set variables, and flatten absorption and reference images
nimgs = size(absimages,3);
nimgsR = size(refimages,3);
xdim = size(absimages(:,:,1),2);
ydim = size(absimages(:,:,1),1);
R = single(reshape(refimages,xdim*ydim,nimgsR));
A = single(reshape(absimages,xdim*ydim,nimgs));
optrefimages=zeros(size(absimages)); % preallocate
if not(exist('bgmask','var')); bgmask=ones(ydim,xdim); end
k = find(bgmask(:)==1); % Index k specifying background region
% Ensure there are no duplicate reference images
% R=unique(R','rows')'; % comment this line if you run out of memory
% Decompose B = R*R' using singular value or LU decomposition
[L,U,p] = lu(R(k,:)'*R(k,:),'vector'); % LU decomposition
for j=1:nimgs
b=R(k,:)'*A(k,j);
% Obtain coefficients c which minimise least-square residuals
lower.LT = true; upper.UT = true;
c = linsolve(U,linsolve(L,b(p,:),lower),upper);
% Compute optimised reference image
optrefimages(:,:,j)=reshape(R*c,[ydim xdim]);
end
end

304
Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/simulateDistribution.m
View File

@ -1,304 +0,0 @@
%% Evolve across a second-order-like transition
clear; clc;
N_params = 50;
N_reps = 50;
alpha_values = linspace(0, 45, N_params);
all_data = cell(1, N_params);
% Transition control
alpha_start = 5; % where sigma starts changing
alpha_widen_end = 15; % when sigma finishes first change
alpha_shift_start = 15; % when mean starts shifting
alpha_end = 40; % when mean finishes shifting and sigma narrows
mu_start = 1.2; % high initial mean
mu_end = 0.2; % low final mean
sigma_start = 0.25; % wide std at start
sigma_mid = 0.15; % mid-range std in middle
sigma_end = 0.07; % narrow std at end
max_skew = 5; % peak skew strength
for i = 1:N_params
alpha = alpha_values(i);
% === Sigma evolution (variance large -> small) ===
if alpha < alpha_start
sigma = sigma_start; % wide at start
elseif alpha < alpha_widen_end
% Smooth transition from wide to mid
t_sigma = (alpha - alpha_start) / (alpha_widen_end - alpha_start);
sigma = sigma_start * (1 - t_sigma) + sigma_mid * t_sigma;
elseif alpha < alpha_end
% Smooth transition from mid to narrow
t_sigma = (alpha - alpha_widen_end) / (alpha_end - alpha_widen_end);
sigma = sigma_mid * (1 - t_sigma) + sigma_end * t_sigma;
else
sigma = sigma_end; % narrow at end
end
% === Mean evolution ===
if alpha < alpha_shift_start
mu = mu_start; % fixed at high initially
elseif alpha <= alpha_end
% Smooth cosine shift
t_mu = (alpha - alpha_shift_start) / (alpha_end - alpha_shift_start);
smooth_t_mu = (1 - cos(pi * t_mu)) / 2;
mu = mu_start * (1 - smooth_t_mu) + mu_end * smooth_t_mu;
else
mu = mu_end;
end
% === Skew evolution ===
if alpha < alpha_end
t_skew = (alpha - alpha_start) / (alpha_end - alpha_start);
skew_strength = max_skew * (1 - t_skew); % fade out
else
skew_strength = 0;
end
% Generate data
if abs(skew_strength) < 1e-2
data = normrnd(mu, sigma, [N_reps, 1]);
else
data = skewnormrnd(mu, sigma, skew_strength, N_reps);
end
all_data{i} = data;
% Cumulants
kappa = computeCumulants(data, 6);
mean_vals(i) = kappa(1);
var_vals(i) = kappa(2);
skew_vals(i) = kappa(3);
kappa4_vals(i) = kappa(4);
kappa5_vals(i) = kappa(5);
kappa6_vals(i) = kappa(6);
end
%% Evolve across a first-order-like transition
% First-order-like distribution evolution with significant bimodality
clear; clc;
N_params = 50;
N_reps = 50;
alpha_values = linspace(0, 45, N_params);
all_data = cell(1, N_params);
% Define transition regions
skewed_start = 10;
bimodal_start = 20;
bimodal_end = 35;
final_narrow_start = 40;
% Peak positions and widths
mu_high = 1.2; % Initial metastable peak
mu_low = 0.2; % Final stable peak
mu_new_peak = 0.8; % New peak appears slightly lower
sigma_initial = 0.08;
for i = 1:N_params
alpha = alpha_values(i);
if alpha < skewed_start
% Region I: Narrow unimodal at high mean
data = normrnd(mu_high, sigma_initial, [N_reps, 1]);
elseif alpha < bimodal_start
% Region II: Slightly skewed
t_skew = (alpha - skewed_start) / (bimodal_start - skewed_start);
mu = mu_high - 0.15 * t_skew;
sigma = sigma_initial + 0.02 * t_skew;
skew_strength = 3 * t_skew;
data = skewnormrnd(mu, sigma, skew_strength, N_reps);
elseif alpha < bimodal_end
% Region III: Bimodal with fixed or slowly drifting peak positions
t = (alpha - bimodal_start) / (bimodal_end - bimodal_start); % t in [0, 1]
% Increased separation between peaks
drift_amount = 0.3; % larger = more drift toward final mean
sep_offset = 0.25; % larger = more initial separation between peaks
% Peaks start separated and move toward mu_low
mu1 = mu_high * (1 - t)^drift_amount + mu_low * (1 - (1 - t)^drift_amount); % Right peak drifts to left
mu2 = (mu_new_peak - sep_offset) * (1 - t)^drift_amount + mu_low * (1 - (1 - t)^drift_amount); % Left peak moves slightly
sigma1 = sigma_initial + 0.02 * (1 - abs(0.5 - t) * 2);
sigma2 = sigma1;
% Weight shift: right peak dies out, left peak grows
w2 = 0.5 + 0.5 * t; % left peak grows: 0.5 → 1
w1 = 1 - w2; % right peak fades: 0.5 → 0
N1 = round(N_reps * w1);
N2 = N_reps - N1;
mode1 = normrnd(mu1, sigma1, [N1, 1]);
mode2 = normrnd(mu2, sigma2, [N2, 1]);
data = [mode1; mode2];
data = data(randperm(length(data)));
else
% Region IV: Final stable low-mean Gaussian
data = normrnd(mu_low, sigma_initial, [N_reps, 1]);
end
% Store data and compute cumulants
all_data{i} = data;
kappa = computeCumulants(data, 6);
mean_vals(i) = kappa(1);
var_vals(i) = kappa(2);
skew_vals(i) = kappa(3);
kappa4_vals(i) = kappa(4);
kappa5_vals(i) = kappa(5);
kappa6_vals(i) = kappa(6);
end
%% === Compute 2D PDF heatmap: f(x, alpha) ===
x_grid = linspace(0.0, 1.8, 200); % max[g²] values on y-axis
pdf_matrix = zeros(numel(x_grid), N_params); % Now: rows = y, columns = alpha
for i = 1:N_params
data = all_data{i};
f = ksdensity(data, x_grid, 'Bandwidth', 0.025);
pdf_matrix(:, i) = f; % Transpose for y-axis to be vertical
end
% === Plot PDF vs. alpha heatmap ===
figure(2); clf;
set(gcf, 'Color', 'w', 'Position',[100 100 950 750])
imagesc(alpha_values, x_grid, pdf_matrix);
set(gca, 'YDir', 'normal'); % Flip y-axis to normal orientation
xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', 'FontSize', 14);
ylabel('$\mathrm{max}[g^{(2)}]$', 'Interpreter', 'latex', 'FontSize', 14);
title('Evolving PDF of $\mathrm{max}[g^{(2)}]$', ...
'Interpreter', 'latex', 'FontSize', 16);
colormap(Colormaps.coolwarm()); % More aesthetic than default
colorbar;
c = colorbar;
ylabel(c, 'PDF', 'FontSize', 14, 'Interpreter', 'latex');
set(gca, 'FontSize', 14);
%% Animate evolving distribution and cumulant value
figure(1); clf;
set(gcf, 'Color', 'w', 'Position',[100 100 1300 750])
for i = 1:N_params
clf;
% PDF
subplot(1,2,1); cla; hold on;
data = all_data{i};
% Plot histogram with normalized PDF
histogram(data, 'Normalization', 'pdf', 'BinWidth', 0.03, ...
'FaceColor', [0.3 0.5 0.8], 'EdgeColor', 'k', 'FaceAlpha', 0.6);
title(sprintf('Histogram at $\\alpha = %.1f^\\circ$', alpha_values(i)), ...
'Interpreter', 'latex', 'FontSize', 16);
xlabel('$\mathrm{max}[g^{(2)}]$', 'Interpreter', 'latex', 'FontSize', 14);
ylabel('PDF', 'FontSize', 14);
set(gca, 'FontSize', 12); grid on;
xlim([0.0, 2.0]);
% Cumulant evolution
subplot(1,2,2); hold on;
plot(alpha_values(1:i), kappa4_vals(1:i), 'bo-', 'LineWidth', 2);
title('Binder Cumulant Tracking', 'Interpreter', 'latex', 'FontSize', 16);
xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', 'FontSize', 14);
ylabel('$\kappa_4$', 'Interpreter', 'latex', 'FontSize', 14);
xlim([0, 45]); grid on;
set(gca, 'FontSize', 12);
pause(0.3);
end
%% === Plotting ===
figure(1)
set(gcf, 'Color', 'w', 'Position', [100 100 950 750])
t = tiledlayout(2, 2, 'TileSpacing', 'compact', 'Padding', 'compact');
scan_vals = alpha_values; % your parameter sweep values
% Define font style for consistency
axis_fontsize = 14;
label_fontsize = 16;
title_fontsize = 16;
% 1. Mean with error bars (if you have error data, else just plot)
% If no error, replace errorbar with plot or omit error data
% For now, no error bars assumed
nexttile;
plot(scan_vals, mean_vals, 'o-', 'LineWidth', 1.5, 'MarkerSize', 6);
title('Mean', 'FontSize', title_fontsize, 'Interpreter', 'latex');
xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', 'FontSize', label_fontsize);
ylabel('$\kappa_1$', 'Interpreter', 'latex', 'FontSize', label_fontsize);
set(gca, 'FontSize', axis_fontsize);
grid on;
% 2. Variance
nexttile;
plot(scan_vals, var_vals, 's-', 'LineWidth', 1.5, 'MarkerSize', 6);
title('Variance', 'FontSize', title_fontsize, 'Interpreter', 'latex');
xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', 'FontSize', label_fontsize);
ylabel('$\kappa_2$', 'Interpreter', 'latex', 'FontSize', label_fontsize);
set(gca, 'FontSize', axis_fontsize);
grid on;
% 3. Skewness
nexttile;
plot(scan_vals, skew_vals, 'd-', 'LineWidth', 1.5, 'MarkerSize', 6);
title('Skewness', 'FontSize', title_fontsize, 'Interpreter', 'latex');
xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', 'FontSize', label_fontsize);
ylabel('$\kappa_3$', 'Interpreter', 'latex', 'FontSize', label_fontsize);
set(gca, 'FontSize', axis_fontsize);
grid on;
% 4. Binder Cumulant
nexttile;
plot(scan_vals, kappa4_vals, '^-', 'LineWidth', 1.5, 'MarkerSize', 6);
title('Binder Cumulant', 'FontSize', title_fontsize, 'Interpreter', 'latex');
xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', 'FontSize', label_fontsize);
ylabel('$\kappa_4$', 'Interpreter', 'latex', 'FontSize', label_fontsize);
set(gca, 'FontSize', axis_fontsize);
grid on;
% Super title (you can customize the string)
sgtitle('Cumulants of a simulated evolving distribution', ...
'FontWeight', 'bold', 'FontSize', 18, 'Interpreter', 'latex');
%% === Helper: Cumulant Calculation ===
function kappa = computeCumulants(data, max_order)
data = data(:);
mu = mean(data);
c = zeros(1, max_order);
centered = data - mu;
for n = 1:max_order
c(n) = mean(centered.^n);
end
kappa = zeros(1, max_order);
kappa(1) = mu;
kappa(2) = c(2);
kappa(3) = c(3);
kappa(4) = c(4) - 3*c(2)^2;
kappa(5) = c(5) - 10*c(3)*c(2);
kappa(6) = c(6) - 15*c(4)*c(2) - 10*c(3)^2 + 30*c(2)^3;
end
%% === Helper: Skewed Normal Distribution ===
function x = skewnormrnd(mu, sigma, alpha, n)
% Skew-normal using Azzalini's method
delta = alpha / sqrt(1 + alpha^2);
u0 = randn(n,1);
v = randn(n,1);
u1 = delta * u0 + sqrt(1 - delta^2) * v;
x = mu + sigma * u1 .* sign(u0);
end

223
Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/understandingCumulants.m
View File

@ -1,223 +0,0 @@
%% Main script: Sweep different parameter pairs
% Default parameters
defaults.mu1 = 0.5;
defaults.mu2 = 1.0;
defaults.sigma1 = 0.1;
defaults.sigma2 = 0.1;
defaults.weight1 = 0.5;
% Parameter pair definitions
%{
param_pairs = {
'mu1', linspace(0.7, 1.0, 40), ...
'mu2', linspace(1.0, 1.3, 40);
'mu1', linspace(0.7, 1.0, 40), ...
'weight1', linspace(0.2, 0.8, 40);
'sigma1', linspace(0.05, 0.2, 40), ...
'sigma2', linspace(0.05, 0.2, 40);
'mu1', linspace(0.7, 1.0, 40), ...
'sigma1', linspace(0.05, 0.2, 40);
'mu2', linspace(1.0, 1.3, 40), ...
'weight1', linspace(0.2, 0.8, 40);
};
%}
param_pairs = {
'mu1', linspace(0.1, 1.5, 40), ...
'mu2', linspace(0.1, 1.5, 40);
};
% Cumulant index to visualize (2=variance, 3=skewness, 4=kurtosis)
cumulant_to_plot = 4;
% Run sweep for each pair
for i = 1:size(param_pairs,1)
param1_name = param_pairs{i,1};
param1_vals = param_pairs{i,2};
param2_name = param_pairs{i,3};
param2_vals = param_pairs{i,4};
fprintf('Sweeping %s and %s...\n', param1_name, param2_name);
Z = sweepBimodalCumulants(param1_name, param1_vals, ...
param2_name, param2_vals, ...
defaults, cumulant_to_plot);
end
%%
% Parameters
N_total = 10000;
mu1 = 0.5;
sigma1 = 0.1;
mu2 = 1.0;
sigma2 = 0.1;
weight1 = 0.7;
% Generate data
data = generateBimodalDistribution(N_total, mu1, mu2, sigma1, sigma2, weight1);
% Plot histogram
figure(3);
clf
set(gcf, 'Color', 'w', 'Position', [100 100 950 750]);
histogram(data, 'Normalization', 'pdf', 'EdgeColor', 'none', 'FaceAlpha', 0.5);
hold on;
% Overlay smooth density estimate
[xi, f] = ksdensity(data);
plot(f, xi, 'r-', 'LineWidth', 2);
% Labels and title
xlabel('Value');
ylabel('Probability Density');
legend('Histogram', 'Smoothed Density');
grid on;
%%
N = size(Z, 1);
main_diag_values = diag(Z);
anti_diag_values = diag(flipud(Z));
param1_diag_main = param1_vals;
param2_diag_main = param2_vals;
param1_diag_anti = param1_vals;
param2_diag_anti = flip(param2_vals);
% For example, plot the main diagonal cumulants:
figure(4);
set(gcf, 'Color', 'w', 'Position', [100 100 950 750]);
plot(1:N, main_diag_values, '-o');
xlabel('Index along diagonal');
ylabel('$\kappa_4$', 'Interpreter', 'latex');
title('$\kappa_4$ along anti-diagonal', 'Interpreter', 'latex');
% Plot anti-diagonal cumulants:
figure(5);
set(gcf, 'Color', 'w', 'Position', [100 100 950 750]);
plot(1:N, anti_diag_values, '-o');
xlabel('Index along anti-diagonal');
ylabel('$\kappa_4$', 'Interpreter', 'latex');
title('$\kappa_4$ along anti-diagonal', 'Interpreter', 'latex');
%% === Helper: Bimodal Distribution ===
function data = generateBimodalDistribution(N_total, mu1, mu2, sigma1, sigma2, weight1)
%GENERATEBIMODALDISTRIBUTION Generates a single bimodal distribution.
%
% data = generateBimodalDistribution(N_total, mu1, mu2, sigma1, sigma2, weight1)
%
% Inputs:
% N_total - total number of samples
% mu1, mu2 - means of the two modes
% sigma1, sigma2 - standard deviations of the two modes
% weight1 - fraction of samples from mode 1 (between 0 and 1)
%
% Output:
% data - shuffled samples from the bimodal distribution
% Validate weight
weight1 = min(max(weight1, 0), 1);
weight2 = 1 - weight1;
% Determine number of samples for each mode
N1 = round(N_total * weight1);
N2 = N_total - N1;
% Generate samples
mode1_samples = normrnd(mu1, sigma1, [N1, 1]);
mode2_samples = normrnd(mu2, sigma2, [N2, 1]);
% Combine and shuffle
data = [mode1_samples; mode2_samples];
data = data(randperm(length(data)));
end
%% === Helper: Cumulant Calculation ===
function kappa = computeCumulants(data, max_order)
data = data(:);
mu = mean(data);
centered = data - mu;
% Preallocate
c = zeros(1, max_order);
kappa = zeros(1, max_order);
% Compute central moments up to max_order
for n = 1:max_order
c(n) = mean(centered.^n);
end
% Assign cumulants based on available order
if max_order >= 1, kappa(1) = mu; end
if max_order >= 2, kappa(2) = c(2); end
if max_order >= 3, kappa(3) = c(3); end
if max_order >= 4, kappa(4) = c(4) - 3*c(2)^2; end
if max_order >= 5, kappa(5) = c(5) - 10*c(3)*c(2); end
if max_order >= 6
kappa(6) = c(6) - 15*c(4)*c(2) - 10*c(3)^2 + 30*c(2)^3;
end
end
%% === Helper: Cumulant Calculation ===
function Z = sweepBimodalCumulants(param1_name, param1_vals, ...
param2_name, param2_vals, ...
fixed_params, ...
cumulant_index)
%SWEEPBIMODALCUMULANTS Sweep 2 parameters and return a chosen cumulant.
%
% Z = sweepBimodalCumulants(...)
% Returns a matrix Z of cumulant values at each grid point.
% Setup grid
[P1, P2] = meshgrid(param1_vals, param2_vals);
Z = zeros(size(P1));
N_samples = 1000;
maxOrder = max(4, cumulant_index);
for i = 1:numel(P1)
% Copy fixed parameters
params = fixed_params;
% Override swept parameters
params.(param1_name) = P1(i);
params.(param2_name) = P2(i);
% Generate and compute cumulants
data = generateBimodalDistribution(N_samples, ...
params.mu1, params.mu2, ...
params.sigma1, params.sigma2, ...
params.weight1);
kappa = computeCumulants(data, maxOrder);
Z(i) = kappa(cumulant_index);
end
% Plot full heatmap
figure;
set(gcf, 'Color', 'w', 'Position', [100 100 950 750]);
imagesc(param1_vals, param2_vals, Z);
set(gca, 'YDir', 'normal');
xlabel(param1_name);
ylabel(param2_name);
title(['Cumulant \kappa_', num2str(cumulant_index)]);
colorbar;
axis tight;
% Optional binary colormap (red = ≥0, blue = <0)
figure;
set(gcf, 'Color', 'w', 'Position', [100 100 950 750]);
imagesc(param1_vals, param2_vals, Z);
set(gca, 'YDir', 'normal');
xlabel(param1_name);
ylabel(param2_name);
title(['Binary color split of \kappa_', num2str(cumulant_index)]);
clim([-1 1]);
colormap([0 0 1; 1 0 0]); % Blue (neg), Red (pos & zero)
colorbar;
axis tight;
end