diff --git a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/analyzeFolder.m b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/analyzeFolder.m index a6a980e..ed483b3 100644 --- a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/analyzeFolder.m +++ b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/analyzeFolder.m @@ -92,8 +92,16 @@ function results = analyzeFolder(options) 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)'; + if (isempty(atm_img) && isa(atm_img, 'double')) || ... + (isempty(bkg_img) && isa(bkg_img, 'double')) || ... + (isempty(dark_img) && isa(dark_img, 'double')) + + refimages = nan(size(refimages)); % fill with NaNs + absimages = nan(size(absimages)); + 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 ===== diff --git a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/bootstrapCumulants.m b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/bootstrapCumulants.m new file mode 100644 index 0000000..e59ec25 --- /dev/null +++ b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/bootstrapCumulants.m @@ -0,0 +1,41 @@ +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 diff --git a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/compareAngularCorrelation.m b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/compareAngularCorrelation.m index a2574a0..cb0030d 100644 --- a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/compareAngularCorrelation.m +++ b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/compareAngularCorrelation.m @@ -8,13 +8,13 @@ format long font = 'Bahnschrift'; % Load data -Data = load('C:/Users/Karthik/Documents/GitRepositories/Calculations/Data-Analyzer/Comparison/Max_g2_DropletsToStripes.mat', 'unique_scan_parameter_values', 'mean_max_g2_values', 'std_error_g2_values'); +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/Comparison/Max_g2_StripesToDroplets.mat', 'unique_scan_parameter_values', 'mean_max_g2_values', '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; @@ -30,10 +30,10 @@ errorbar(std_scan_parameter_values, std_mean_mg2, std_stderr_mg2, 'o--', ... 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'); +% 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 +% set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title grid on %% \ No newline at end of file diff --git a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/computeCumulants.m b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/computeCumulants.m new file mode 100644 index 0000000..52ed552 --- /dev/null +++ b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/computeCumulants.m @@ -0,0 +1,52 @@ +function cumulants = computeCumulants(x, maxOrder) +% computeCumulants - compute cumulants up to specified order from data vector x +% +% Syntax: cumulants = computeCumulants(x, maxOrder) +% +% Inputs: +% x - 1D numeric vector (may contain NaNs) +% maxOrder - maximum order of cumulants to compute (default: 6) +% +% Output: +% cumulants - vector [kappa_1, ..., kappa_maxOrder] + + if nargin < 2 + maxOrder = 6; + end + + x = x(:); + x = x(~isnan(x)); % Remove NaNs + + if isempty(x) + cumulants = NaN(1, maxOrder); + return; + end + + mu1 = mean(x, 'omitnan'); + x_centered = x - mu1; + + cumulants = zeros(1, maxOrder); + cumulants(1) = mu1; + + mu = zeros(1, maxOrder); + for k = 2:maxOrder + mu(k) = mean(x_centered.^k, 'omitnan'); + end + + if maxOrder >= 2 + cumulants(2) = mu(2); + end + if maxOrder >= 3 + cumulants(3) = mu(3); + end + if maxOrder >= 4 + cumulants(4) = mu(4) - 3 * mu(2)^2; + end + if maxOrder >= 5 + cumulants(5) = mu(5) - 10 * mu(3) * mu(2); + end + if maxOrder >= 6 + cumulants(6) = mu(6) - 15 * mu(4) * mu(2) - 10 * mu(3)^2 + 30 * mu(2)^3; + end + +end diff --git a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/conductSpectralAnalysis.m b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/conductSpectralAnalysis.m index 839cf77..e0fb551 100644 --- a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/conductSpectralAnalysis.m +++ b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/conductSpectralAnalysis.m @@ -3,9 +3,9 @@ groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1 "/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/"; +folderPath = "//DyLabNAS/Data/TwoDGas/2025/07/22/"; -run = '0300'; +run = '0021'; folderPath = strcat(folderPath, run); @@ -50,32 +50,32 @@ savefileName = 'DropletsToStripes'; font = 'Bahnschrift'; if strcmp(savefileName, 'DropletsToStripes') - scan_groups = 0:5:45; + scan_groups = 0:1:40; titleString = 'Droplets to Stripes'; elseif strcmp(savefileName, 'StripesToDroplets') - scan_groups = 45:-5:0; + scan_groups = 40:-1:0; titleString = 'Stripes to Droplets'; end % Flags -skipNormalization = false; -skipUnshuffling = true; +skipNormalization = true; +skipUnshuffling = false; skipPreprocessing = true; skipMasking = true; skipIntensityThresholding = true; skipBinarization = true; skipMovieRender = true; -skipSaveFigures = false; -skipSaveOD = false; +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/"; +folderPath = "//DyLabNAS/Data/TwoDGas/2025/07/23/"; -run = '0001'; +run = '0055'; folderPath = strcat(folderPath, run); @@ -120,10 +120,10 @@ savefileName = 'StripesToDroplets'; font = 'Bahnschrift'; if strcmp(savefileName, 'DropletsToStripes') - scan_groups = 0:5:45; + scan_groups = 0:1:40 titleString = 'Droplets to Stripes'; elseif strcmp(savefileName, 'StripesToDroplets') - scan_groups = 45:-5:0; + scan_groups = 40:-1:0; titleString = 'Stripes to Droplets'; end @@ -135,8 +135,8 @@ skipMasking = true; skipIntensityThresholding = true; skipBinarization = true; skipMovieRender = true; -skipSaveFigures = false; -skipSaveOD = false; +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. @@ -156,8 +156,16 @@ for k = 1 : length(files) 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)'; + if (isempty(atm_img) && isa(atm_img, 'double')) || ... + (isempty(bkg_img) && isa(bkg_img, 'double')) || ... + (isempty(dark_img) && isa(dark_img, 'double')) + + refimages = nan(size(refimages)); % fill with NaNs + absimages = nan(size(absimages)); + 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 ===== diff --git a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/extractAutocorrelation.m b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/extractAutocorrelation.m index 477f812..2085fab 100644 --- a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/extractAutocorrelation.m +++ b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/extractAutocorrelation.m @@ -86,8 +86,16 @@ for k = 1 : length(files) 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)'; + if (isempty(atm_img) && isa(atm_img, 'double')) || ... + (isempty(bkg_img) && isa(bkg_img, 'double')) || ... + (isempty(dark_img) && isa(dark_img, 'double')) + + refimages = nan(size(refimages)); % fill with NaNs + absimages = nan(size(absimages)); + 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 ===== @@ -209,8 +217,8 @@ for i = 1:N_params temp(j) = num / denom; end - g2_all(i, dtheta+1) = mean(temp); - g2_error_all(i, dtheta+1) = std(temp) / sqrt(length(group_idx)); % Standard error + 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 @@ -222,7 +230,7 @@ cmap = sky(nParams); % You can also try 'jet', 'turbo', 'hot', etc. figure(1); clf; -set(gcf,'Position',[100 100 950 750]) +set(gcf, 'Color', 'w', 'Position',[100 100 950 750]) hold on; legend_entries = cell(nParams, 1); diff --git a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/extractCustomCorrelation.m b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/extractCustomCorrelation.m index 66ca5a4..f17b79b 100644 --- a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/extractCustomCorrelation.m +++ b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/extractCustomCorrelation.m @@ -3,9 +3,9 @@ groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1 "/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/"; +folderPath = "//DyLabNAS/Data/TwoDGas/2025/07/22/"; -run = '0300'; +run = '0021'; folderPath = strcat(folderPath, run); @@ -86,8 +86,16 @@ for k = 1 : length(files) 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)'; + if (isempty(atm_img) && isa(atm_img, 'double')) || ... + (isempty(bkg_img) && isa(bkg_img, 'double')) || ... + (isempty(dark_img) && isa(dark_img, 'double')) + + refimages = nan(size(refimages)); % fill with NaNs + absimages = nan(size(absimages)); + 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 ===== @@ -176,7 +184,7 @@ 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); +N_params = length(unique_scan_parameter_values); % Define angular range and conversion angle_range = 180; @@ -189,23 +197,24 @@ 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); +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); +kurt_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 -g2_all_per_group = cell(1, N_params); -angle_all_per_group = cell(1, N_params); +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); - angle_at_max_g2 = zeros(1, N_reps); + 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, :); @@ -225,8 +234,7 @@ for i = 1:N_params 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; @@ -237,36 +245,191 @@ for i = 1:N_params 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; + max_g2_all_per_group{i} = g2_values; + + % Compute cumulants + kappa = computeCumulants(g2_values(:), 6); % 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'); + 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); + kurt_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 -%% ── Mean ± Std vs. scan parameter ────────────────────────────────────── -% Compute standard error instead of standard deviation -std_error_g2_values = zeros(1, N_params); -for i = 1:N_params - 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); +%% 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 = 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. Kurtosis +subplot(3,2,4); +plot(scan_vals, kurt_max_g2_angle_values, '^-', 'LineWidth', 1.5, 'MarkerSize', 6); +title('Kurtosis 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,'Position',[100 100 950 750]) +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) @@ -289,75 +452,6 @@ if ~exist(saveFolder, 'dir') end save([saveFolder savefileName '.mat'], 'unique_scan_parameter_values', 'mean_max_g2_values', 'std_error_g2_values'); -%{ -%% Plot histograms within 0-180 degrees only -figure(1); -hold on; - -bin_edges = 0:10:180; - -h1 = histogram(ref_peak_angles, 'BinEdges', bin_edges, ... - 'FaceColor', [0.3 0.7 0.9], 'EdgeColor', 'none', 'FaceAlpha', 0.6); - -h2 = histogram(angle_at_max_g2, 'BinEdges', bin_edges, ... - 'FaceColor', [0.9 0.4 0.4], 'EdgeColor', 'none', 'FaceAlpha', 0.6); - -h1.Normalization = 'probability'; -h2.Normalization = 'probability'; - -xlabel('Angle (degrees)', 'FontSize', 12); -ylabel('Probability', 'FontSize', 12); -legend({'Reference Peak Angle', 'Angle at Max g₂'}, 'FontSize', 12); -title('Comparison of Reference Peak and Max g₂ Angles', 'FontSize', 14); -grid on; -xlim([0 180]); -hold off; - -% Assume ref_peak_angles and angle_at_max_g2 are row or column vectors of angles in [0,180] - -% Define fine angle grid for KDE evaluation -angle_grid = linspace(0, 180, 1000); - -% KDE for reference peak angles -[f_ref, xi_ref] = ksdensity(ref_peak_angles, angle_grid, 'Bandwidth', 5); - -% KDE for max g2 angles -[f_g2, xi_g2] = ksdensity(angle_at_max_g2, angle_grid, 'Bandwidth', 5); - -% Plot KDEs -figure(2); -plot(xi_ref, f_ref, 'LineWidth', 2, 'DisplayName', 'Reference Peak Angles'); -hold on; -plot(xi_g2, f_g2, 'LineWidth', 2, 'DisplayName', 'Max g_2 Angles'); -xlabel('Angle (degrees)'); -ylabel('Probability Density'); -title('KDE of Angle Distributions'); -legend; -grid on; - -% Find modes (angle at max KDE value) -[~, mode_idx_ref] = max(f_ref); -mode_ref = xi_ref(mode_idx_ref); - -[~, mode_idx_g2] = max(f_g2); -mode_g2 = xi_g2(mode_idx_g2); - -% Calculate difference in mode -mode_diff = abs(mode_ref - mode_g2); -fprintf('Mode difference between distributions: %.2f degrees\n', mode_diff); - -% Add vertical dashed lines at mode positions -yl = ylim; % get y-axis limits for text positioning - -% Reference peak mode line and label -xline(mode_ref, 'k--', 'LineWidth', 1.5, 'DisplayName', sprintf('Ref Mode: %.1f°', mode_ref)); -text(mode_ref, yl(2)*0.9, sprintf('%.1f°', mode_ref), 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 12, 'Color', 'k'); - -% Max g2 mode line and label -xline(mode_g2, 'r--', 'LineWidth', 1.5, 'DisplayName', sprintf('g_2 Mode: %.1f°', mode_g2)); -text(mode_g2, yl(2)*0.75, sprintf('%.1f°', mode_g2), 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 12, 'Color', 'r'); -%} - %% Helper Functions function [IMGFFT, IMGPR] = computeFourierTransform(I, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization) % computeFourierSpectrum - Computes the 2D Fourier power spectrum