From e671d129695e7233a27296645b6ab1077edc1b85 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Wed, 15 Oct 2025 17:55:05 +0200 Subject: [PATCH] Added routines to analyze radial spectral distributions and corrected other bugs. --- .../+Analyzer/computeSpectralAverages.m | 50 ++++ .../+Analyzer/conductSpectralAnalysis.m | 4 +- .../extractFullRadialSpectralDistribution.m | 95 +++++++ ...sianCurvesToAngularSpectralDistribution.m} | 4 +- ...ussianCurvesToRadialSpectralDistribution.m | 146 ++++++++++ .../computeAngularSpectralDistribution.m | 24 +- .../computeRadialSpectralDistribution.m | 18 +- .../+Plotter/plotSecondaryGaussianRange.m | 166 +++++++++++ ...verageSpectra.m => plotSpectralAverages.m} | 79 +++--- Data-Analyzer/+Plotter/plotSpectralCurves.m | 39 ++- .../plotAnalysisResults.m | 17 +- ... runAngularSpectralDistributionAnalysis.m} | 15 +- .../runCorrelationAnalysis.m | 5 +- .../runRadialSpectralDistributionAnalysis.m | 257 ++++++++++++++++++ .../plotAnalysisResults.m | 16 +- ... runAngularSpectralDistributionAnalysis.m} | 8 +- .../runCorrelationAnalysis.m | 3 +- .../runRadialSpectralDistributionAnalysis.m | 257 ++++++++++++++++++ 18 files changed, 1074 insertions(+), 129 deletions(-) create mode 100644 Data-Analyzer/+Analyzer/computeSpectralAverages.m create mode 100644 Data-Analyzer/+Analyzer/extractFullRadialSpectralDistribution.m rename Data-Analyzer/+Analyzer/{fitTwoGaussianCurves.m => fitTwoGaussianCurvesToAngularSpectralDistribution.m} (98%) create mode 100644 Data-Analyzer/+Analyzer/fitTwoGaussianCurvesToRadialSpectralDistribution.m create mode 100644 Data-Analyzer/+Plotter/plotSecondaryGaussianRange.m rename Data-Analyzer/+Plotter/{plotAverageSpectra.m => plotSpectralAverages.m} (65%) rename Data-Analyzer/+Scripts/BECToDropletsToStripes/{runSpectralDistributionAnalysis.m => runAngularSpectralDistributionAnalysis.m} (97%) create mode 100644 Data-Analyzer/+Scripts/BECToDropletsToStripes/runRadialSpectralDistributionAnalysis.m rename Data-Analyzer/+Scripts/BECToStripesToDroplets/{runSpectralDistributionAnalysis.m => runAngularSpectralDistributionAnalysis.m} (97%) create mode 100644 Data-Analyzer/+Scripts/BECToStripesToDroplets/runRadialSpectralDistributionAnalysis.m diff --git a/Data-Analyzer/+Analyzer/computeSpectralAverages.m b/Data-Analyzer/+Analyzer/computeSpectralAverages.m new file mode 100644 index 0000000..85daf6b --- /dev/null +++ b/Data-Analyzer/+Analyzer/computeSpectralAverages.m @@ -0,0 +1,50 @@ +function results = computeSpectralAverages(S_all, scan_reference_values) +%% computeSpectralAverages +% Compute mean and SEM of spectra grouped by scan parameter. +% Handles both 1D curves and 2D matrices. + + [uniqueParams, ~, idx] = unique(scan_reference_values); + nParams = numel(uniqueParams); + nTotal = numel(S_all); + nReps = nTotal / nParams; + + % Determine data dimensionality + firstData = S_all{1}; + dataSize = size(firstData); + isMatrix = (numel(size(firstData)) == 2) && all(size(firstData) > 1); + + curves_all = cell(1, nParams); + curves_mean = cell(1, nParams); + curves_err = cell(1, nParams); + + for i = 1:nParams + if isMatrix + % --- Case: 2D matrix (e.g., power spectra) --- + data_sum = 0; + for r = 1:nReps + idx_r = (r-1)*nParams + i; + data_sum = data_sum + S_all{idx_r}; + end + avg = data_sum / nReps; + curves_all{i} = []; % not used for 2D + curves_mean{i} = avg; + curves_err{i} = []; + else + % --- Case: 1D curve --- + nPoints = numel(S_all{1}); + curves = zeros(nReps, nPoints); + for r = 1:nReps + idx_r = (r-1)*nParams + i; + curves(r,:) = S_all{idx_r}; + end + curves_all{i} = curves; + curves_mean{i} = mean(curves,1); + curves_err{i} = std(curves,0,1)/sqrt(nReps); + end + end + + results.curves_all = curves_all; + results.curves_mean = curves_mean; + results.curves_error = curves_err; + results.scan_parameters = uniqueParams; +end diff --git a/Data-Analyzer/+Analyzer/conductSpectralAnalysis.m b/Data-Analyzer/+Analyzer/conductSpectralAnalysis.m index fcb1a75..2a3aadf 100644 --- a/Data-Analyzer/+Analyzer/conductSpectralAnalysis.m +++ b/Data-Analyzer/+Analyzer/conductSpectralAnalysis.m @@ -147,11 +147,11 @@ function results = conductSpectralAnalysis(od_imgs, scan_parameter_values, optio % Store results 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; diff --git a/Data-Analyzer/+Analyzer/extractFullRadialSpectralDistribution.m b/Data-Analyzer/+Analyzer/extractFullRadialSpectralDistribution.m new file mode 100644 index 0000000..50f3095 --- /dev/null +++ b/Data-Analyzer/+Analyzer/extractFullRadialSpectralDistribution.m @@ -0,0 +1,95 @@ +function results = extractFullRadialSpectralDistribution(od_imgs, options) +%% extractFullRadialSpectralDistribution +% Author: Karthik +% Date: 2025-10-15 +% Version: 1.0 +% +% Description: +% Computes the radial spectral distribution S(k) from OD images +% +% Inputs: +% od_imgs - cell array of OD images +% options - struct containing relevant parameters: +% pixel_size, magnification, zoom_size, +% theta_min, theta_max, N_radial_bins, Radial_WindowSize +% +% Outputs: +% results - struct containing k_rho values and radial spectra +% +% Notes: +% This function is a minimal variant of conductSpectralAnalysis, +% stopping right after radial spectral dsitribution extraction. + + %% ===== Unpack options ===== + pixel_size = options.pixel_size; + magnification = options.magnification; + zoom_size = options.zoom_size; + theta_min = options.theta_min; + theta_max = options.theta_max; + N_radial_bins = options.N_radial_bins; + radial_window_size = options.Radial_WindowSize; + skipPreprocessing = options.skipPreprocessing; + skipMasking = options.skipMasking; + skipIntensityThresholding = options.skipIntensityThresholding; + skipBinarization = options.skipBinarization; + + %% ===== Initialization ===== + N_shots = length(od_imgs); + fft_imgs = cell(1, N_shots); + S_k_all = cell(1, N_shots); + S_k_smoothed_all = cell(1, N_shots); + k_rho_vals = []; + + %% ===== Main loop ===== + for k = 1:N_shots + IMG = od_imgs{k}; + + % Skip empty or low-intensity images + if ~(max(IMG(:)) > 1) + IMGFFT = NaN(size(IMG)); + else + % Compute FFT (with same preprocessing pipeline) + [IMGFFT, ~] = Calculator.computeFourierTransform(IMG, ... + skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization); + end + + % Image dimensions + [Ny, Nx] = size(IMG); + dx = pixel_size / magnification; + dy = dx; + + % Reciprocal-space sampling + 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; + + % Convert to wavenumber axes [µm⁻¹] + kx_full = 2 * pi * vx * 1E-6; + ky_full = 2 * pi * vy * 1E-6; + + % Crop FFT and wavenumber axes 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); + kx = kx_full(mid_x - zoom_size : mid_x + zoom_size); + ky = ky_full(mid_y - zoom_size : mid_y + zoom_size); + + % ===== Compute 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_window_size); + + % Store results + S_k_all{k} = S_k; + S_k_smoothed_all{k} = S_k_smoothed; + end + + %% ===== Package results ===== + results = struct(); + results.k_rho_vals = k_rho_vals; + results.S_k_all = S_k_all; + results.S_k_smoothed_all = S_k_smoothed_all; +end diff --git a/Data-Analyzer/+Analyzer/fitTwoGaussianCurves.m b/Data-Analyzer/+Analyzer/fitTwoGaussianCurvesToAngularSpectralDistribution.m similarity index 98% rename from Data-Analyzer/+Analyzer/fitTwoGaussianCurves.m rename to Data-Analyzer/+Analyzer/fitTwoGaussianCurvesToAngularSpectralDistribution.m index 460ea36..5428b1b 100644 --- a/Data-Analyzer/+Analyzer/fitTwoGaussianCurves.m +++ b/Data-Analyzer/+Analyzer/fitTwoGaussianCurvesToAngularSpectralDistribution.m @@ -1,5 +1,5 @@ -function [fitResults, rawCurves] = fitTwoGaussianCurves(S_theta_all, theta_vals, varargin) -%% fitTwoGaussianCurves +function [fitResults, rawCurves] = fitTwoGaussianCurvesToAngularSpectralDistribution(S_theta_all, theta_vals, varargin) +%% fitTwoGaussianCurvesToAngularSpectralDistribution % Fits a two-Gaussian model to multiple spectral curves from θ=0 up to % the first secondary peak ≥50% of the primary peak, within 0–π/2 radians. % diff --git a/Data-Analyzer/+Analyzer/fitTwoGaussianCurvesToRadialSpectralDistribution.m b/Data-Analyzer/+Analyzer/fitTwoGaussianCurvesToRadialSpectralDistribution.m new file mode 100644 index 0000000..1ee4417 --- /dev/null +++ b/Data-Analyzer/+Analyzer/fitTwoGaussianCurvesToRadialSpectralDistribution.m @@ -0,0 +1,146 @@ +function [fitResults, rawCurves] = fitTwoGaussianCurvesToRadialSpectralDistribution(S_k_all, k_rho_vals, varargin) +%% fitTwoGaussianCurvesToRadialSpectralDistribution +% Fits a two-Gaussian model to multiple radial spectral curves (S(k_rho)), +% truncating each curve to user-specified k_rho and amplitude ranges. +% +% Author: Karthik +% Date: 2025-10-15 +% Version: 1.0 +% +% Inputs: +% S_k_all - 1xN cell array of radial spectral curves +% k_rho_vals - vector of corresponding k_rho values +% +% Optional name-value pairs: +% 'KRhoRange' - [k_min, k_max] (default: [min(k_rho_vals), max(k_rho_vals)]) +% 'AmplitudeRange' - [y_min, y_max] (default: [0, Inf]) +% 'ResidualThreshold', 'PositionThreshold', 'AmplitudeThreshold' as before +% +% Outputs: +% fitResults - struct array with fields: +% .pFit, .kFit, .xFit, .yFit, .kFine, .isValid, .fitMaxKRho +% rawCurves - struct array of truncated raw curves for plotting + + % --- Parse optional inputs --- + p = inputParser; + addParameter(p, 'KRhoRange', [min(k_rho_vals), max(k_rho_vals)], @(x) isnumeric(x) && numel(x)==2); + addParameter(p, 'AmplitudeRange', [0, Inf], @(x) isnumeric(x) && numel(x)==2); + addParameter(p, 'ResidualThreshold', 0.3, @(x) isnumeric(x) && isscalar(x)); + addParameter(p, 'PositionThreshold', 0.1, @(x) isnumeric(x) && isscalar(x)); + addParameter(p, 'AmplitudeThreshold', 0.5, @(x) isnumeric(x) && isscalar(x)); + parse(p, varargin{:}); + opts = p.Results; + + Ncurves = numel(S_k_all); + fitResults = struct('pFit',[],'kFit',[],'xFit',[],'yFit',[],... + 'kFine',[],'isValid',[],'fitMaxKRho',[]); + rawCurves = struct('x', cell(1, Ncurves), 'k', cell(1, Ncurves)); + + for k = 1:Ncurves + % --- Truncate curve to k_rho and amplitude ranges --- + xRaw = S_k_all{k}; + % --- Normalize and smooth --- + xNorm = xRaw / max(xRaw); + xSmooth = smooth(xNorm, 5, 'moving'); + + mask = (k_rho_vals(:) >= opts.KRhoRange(1)) & (k_rho_vals(:) <= opts.KRhoRange(2)) & ... + (xSmooth(:) >= opts.AmplitudeRange(1)) & (xSmooth(:) <= opts.AmplitudeRange(2)); + + x = xSmooth(mask); + kVals = k_rho_vals(mask); + + rawCurves(k).x = x; + rawCurves(k).k = kVals; + + if isempty(x) + warning('Curve %d is empty after truncation.', k); + fitResults(k) = fillZeroStructRadial(kVals, opts.KRhoRange(2)); + continue; + elseif numel(x) < 5 + warning('Curve %d has too few points after truncation, skipping.', k); + fitResults(k) = fillNaNStructRadial(); + continue; + end + + try + % --- Detect peaks --- + curveAmp = max(x) - min(x); + minProm = opts.AmplitudeThreshold * curveAmp; + [pk, locIdx] = findpeaks(x, 'MinPeakProminence', minProm); + kPeaks = kVals(locIdx); + + %% Generic two-Gaussian + twoGauss = @(p,k) ... + p(1)*exp(-0.5*((k - p(2))/max(p(3),1e-6)).^2) + ... + p(4)*exp(-0.5*((k - p(5))/max(p(6),1e-6)).^2); + + secondaryIdx = find(kPeaks > opts.PositionThreshold, 1, 'first'); + if ~isempty(secondaryIdx) + kSecondary = kPeaks(secondaryIdx); + secondaryAmp = pk(secondaryIdx); + A1_guess = opts.AmplitudeRange(2); + mu1_guess = 0.0; + sigma1_guess = 0.1 * kSecondary; + A2_guess = secondaryAmp; + mu2_guess = kSecondary; + sigma2_guess = 0.1 * kSecondary; + p0 = [A1_guess, mu1_guess, sigma1_guess, A2_guess, mu2_guess, sigma2_guess]; + else + p0 = [opts.AmplitudeRange(2), 0, 0.1, 0.01, opts.KRhoRange(2)/3, 0.1]; + end + + lb = [0, 0, 1e-6, 0, opts.PositionThreshold, 1e-6]; + ub = [2, opts.KRhoRange(2)/2, opts.KRhoRange(2), 2, opts.KRhoRange(2), opts.KRhoRange(2)]; + optsLSQ = optimoptions('lsqcurvefit','Display','off', ... + 'MaxFunctionEvaluations',1e4,'MaxIterations',1e4); + pFit = lsqcurvefit(twoGauss, p0, kVals(:), x(:), lb, ub, optsLSQ); + + if pFit(4) <= 1e-3 + pFit(4) = 0; + pFit(5:6) = NaN; + end + + fitMaxKRho = max(pFit(5), pFit(2)); + kFine = linspace(min(kVals), max(kVals), 500); + yFit = twoGauss(pFit, kFine); + + fitResults(k).pFit = pFit; + fitResults(k).kFit = kVals; + fitResults(k).xFit = x; + fitResults(k).kFine = kFine; + fitResults(k).yFit = yFit; + fitResults(k).isValid = true; + fitResults(k).fitMaxKRho = fitMaxKRho; + + catch + warning('Curve %d fit failed completely.', k); + fitResults(k) = fillNaNStructRadial(); + end + end +end + +%% --- Helper: NaN struct for failed fits --- +function s = fillNaNStructRadial() + s = struct( ... + 'pFit', nan(1,6), ... + 'kFit', nan, ... + 'xFit', nan, ... + 'yFit', nan, ... + 'kFine', nan, ... + 'isValid', false, ... + 'fitMaxKRho', nan ... + ); +end + +%% --- Helper: Zero struct for empty curves --- +function s = fillZeroStructRadial(kVals, MaxKRho) + s = struct( ... + 'pFit', zeros(1,6), ... + 'kFit', kVals, ... + 'xFit', zeros(size(kVals)), ... + 'yFit', zeros(size(kVals)), ... + 'kFine', zeros(1,500), ... + 'isValid', false, ... + 'fitMaxKRho', MaxKRho ... + ); +end diff --git a/Data-Analyzer/+Calculator/computeAngularSpectralDistribution.m b/Data-Analyzer/+Calculator/computeAngularSpectralDistribution.m index 10feea7..cc63aba 100644 --- a/Data-Analyzer/+Calculator/computeAngularSpectralDistribution.m +++ b/Data-Analyzer/+Calculator/computeAngularSpectralDistribution.m @@ -16,28 +16,28 @@ function [theta_vals, S_theta] = computeAngularSpectralDistribution(IMGFFT, kx, 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] + [KX, KY] = meshgrid(kx, ky); + Kmag = sqrt(KX.^2 + KY.^2); % radial wavenumber magnitude + Theta = atan2(KY, KX); % range [-pi, pi] % Map Theta into [0, 2*pi) Theta(Theta < 0) = Theta(Theta < 0) + 2*pi; % Restrict to radial band in wavenumber space - radial_mask = (Kmag >= k_min) & (Kmag <= k_max); + radial_mask = (Kmag >= k_min) & (Kmag <= k_max); % Initialize angular structure factor - S_theta = zeros(1, num_bins); - theta_vals = linspace(0, theta_max, num_bins); + S_theta = zeros(1, num_bins); + theta_vals = linspace(0, theta_max, num_bins); % Loop over angular bins for i = 1:num_bins - angle_start = (i - 1) * theta_max / num_bins; - angle_end = i * theta_max / 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)); + angle_start = (i - 1) * theta_max / num_bins; + angle_end = i * theta_max / 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 diff --git a/Data-Analyzer/+Calculator/computeRadialSpectralDistribution.m b/Data-Analyzer/+Calculator/computeRadialSpectralDistribution.m index 58c49e1..e357471 100644 --- a/Data-Analyzer/+Calculator/computeRadialSpectralDistribution.m +++ b/Data-Analyzer/+Calculator/computeRadialSpectralDistribution.m @@ -22,8 +22,8 @@ function [k_rho_vals, S_radial] = computeRadialSpectralDistribution(IMGFFT, kx, % Optional notes, references. [KX, KY] = meshgrid(kx, ky); - K_rho = sqrt(KX.^2 + KY.^2); - Theta = atan2(KY, KX); + K_rho = sqrt(KX.^2 + KY.^2); + Theta = atan2(KY, KX); if thetamin < thetamax angle_mask = (Theta >= thetamin) & (Theta <= thetamax); @@ -33,17 +33,17 @@ function [k_rho_vals, S_radial] = computeRadialSpectralDistribution(IMGFFT, kx, 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); + 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); + S_radial = zeros(1, num_bins); for i = 1:num_bins - r_low = r_edges(i); - r_high = r_edges(i + 1); + 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; + full_mask = radial_mask & angle_mask; S_radial(i) = sum(power_spectrum(full_mask)); end end \ No newline at end of file diff --git a/Data-Analyzer/+Plotter/plotSecondaryGaussianRange.m b/Data-Analyzer/+Plotter/plotSecondaryGaussianRange.m new file mode 100644 index 0000000..5931dc4 --- /dev/null +++ b/Data-Analyzer/+Plotter/plotSecondaryGaussianRange.m @@ -0,0 +1,166 @@ +function plotSecondaryGaussianRange(fitResults, scanValues, varargin) +%% plotSecondaryGaussianRange +% Plots secondary Gaussian ranges (mu2 ± sigma2) in a heatmap style +% exactly like plotFitParameterPDF, with dashed lines for distribution +% edges and solid line for mean mu2. + +% --- 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', 'SecondaryGaussianRange.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, 'YLim', [], @(x) isempty(x) || numel(x)==2); +addParameter(p, 'Colormap', @jet); +addParameter(p, 'PlotType', 'histogram', @(x) any(validatestring(x,{'kde','histogram'}))); +addParameter(p, 'NumberOfBins', 50, @isscalar); +addParameter(p, 'NormalizeHistogram', true, @islogical); +addParameter(p, 'OverlayMeanSEM', true, @islogical); +parse(p, varargin{:}); +opts = p.Results; + +% --- Determine repetitions and scan parameters --- +N_params = numel(scanValues); +N_total = numel(fitResults); +N_reps = N_total / N_params; + +% --- Extract mu2 and sigma2 --- +mu2_values = nan(N_reps, N_params); +sigma2_values = nan(N_reps, N_params); +for k = 1:N_total + paramIdxScan = mod(k-1, N_params) + 1; + repIdx = floor((k-1)/N_params) + 1; + mu2_values(repIdx, paramIdxScan) = fitResults(k).pFit(5); + sigma2_values(repIdx, paramIdxScan) = fitResults(k).pFit(6); +end + +% --- Mask invalid fits --- +for i = 1:N_total + paramIdxScan = mod(i-1, N_params) + 1; + repIdx = floor((i-1)/N_params) + 1; + if ~fitResults(i).isValid && all(fitResults(i).pFit == 0) + mu2_values(repIdx, paramIdxScan) = NaN; + sigma2_values(repIdx, paramIdxScan) = NaN; + end +end + +% --- Compute lower and upper bounds --- +lowerBound = mu2_values - sigma2_values; +upperBound = mu2_values + sigma2_values; + +% Define how much to trim (e.g., 2% on each side) +trimFraction = 0.01; % can tune this (0.01–0.05 typical) + +% Compute robust global limits +globalLower = quantile(lowerBound(:), trimFraction); +globalUpper = quantile(upperBound(:), 1 - trimFraction); + +fprintf('Global lower bound (%.0f%%): %.4f\n', trimFraction*100, globalLower); +fprintf('Global upper bound (%.0f%%): %.4f\n', (1-trimFraction)*100, globalUpper); + +% --- Prepare data per scan parameter --- +dataCell = cell(N_params,1); +for i = 1:N_params + combinedData = []; + for j = 1:N_reps + if ~isnan(mu2_values(j,i)) && ~isnan(sigma2_values(j,i)) + combinedData = [combinedData; linspace(lowerBound(j,i), upperBound(j,i), 3)']; %#ok + end + end + dataCell{i} = combinedData; +end + +% --- Determine y-range --- +if isempty(opts.DataRange) + allData = cell2mat(dataCell(:)); + y_min = min(allData); + y_max = max(allData); +else + y_min = opts.DataRange(1); + y_max = opts.DataRange(2); +end + +% --- Prepare PDF matrix --- +if strcmpi(opts.PlotType,'kde') + y_grid = linspace(y_min, y_max, opts.NumPoints); + pdf_matrix = zeros(numel(y_grid), N_params); + for i = 1:N_params + data = dataCell{i}; + if isempty(data), continue; end + pdf_matrix(:,i) = ksdensity(data, y_grid); + end +else + edges = linspace(y_min, y_max, opts.NumberOfBins+1); + binCenters = (edges(1:end-1)+edges(2:end))/2; + pdf_matrix = zeros(numel(binCenters), N_params); + for i = 1:N_params + data = dataCell{i}; + if isempty(data), continue; end + counts = histcounts(data, edges); + if opts.NormalizeHistogram + counts = counts / sum(counts); + end + pdf_matrix(:,i) = counts(:); + end +end + +% --- Mask out non-data regions --- +maskPerScan = cellfun(@(c) ~isempty(c), dataCell); +pdf_matrix(:, ~maskPerScan) = NaN; + +% --- Plot heatmap --- +fig = figure(opts.FigNum); clf(fig); +set(fig, 'Color', 'w', 'Position', [100 100 950 750]); + +if strcmpi(opts.PlotType,'kde') + h = imagesc(scanValues, y_grid, pdf_matrix); + set(h, 'AlphaData', ~isnan(pdf_matrix)); +else + h = imagesc(scanValues, binCenters, pdf_matrix); + set(h, 'AlphaData', ~isnan(pdf_matrix)); +end + +colormap([1 1 1; feval(opts.Colormap)]); % white for NaNs +c = colorbar; +ylabel(c, 'Probability Density', 'Rotation', -90, 'FontName', opts.FontName, 'FontSize', opts.FontSize); + +set(gca, 'YDir','normal', 'FontName', opts.FontName, 'FontSize', opts.FontSize); +xlabel(opts.XLabel, 'FontName', opts.FontName, 'FontSize', opts.FontSize); +ylabel(opts.YLabel, 'FontName', opts.FontName, 'FontSize', opts.FontSize); +title(opts.Title, 'FontName', opts.FontName, 'FontSize', opts.FontSize+2, 'FontWeight','bold'); +grid on; hold on; + +% --- Overlay mu2 and distribution edges --- +meanMu2 = nanmean(mu2_values,1); +minLower = min(lowerBound,[],1); +maxUpper = max(upperBound,[],1); + +plot(scanValues, meanMu2, 'k-', 'LineWidth', 2); % solid mu2 +plot(scanValues, minLower, 'k--', 'LineWidth', 1.5); % dashed lower edge +plot(scanValues, maxUpper, 'k--', 'LineWidth', 1.5); % dashed upper edge + +% --- Overlay mean ± SEM if requested --- +if opts.OverlayMeanSEM + semMu2 = nanstd(mu2_values,0,1) ./ sqrt(sum(~isnan(mu2_values),1)); + xVec = reshape(scanValues,1,[]); + fill([xVec, fliplr(xVec)], [meanMu2 - semMu2, fliplr(meanMu2 + semMu2)], ... + [0.2 0.4 0.8], 'FaceAlpha',0.2, 'EdgeColor','none'); +end + +% --- Apply axis limits if specified --- +if ~isempty(opts.XLim) + xlim(opts.XLim); +end +if ~isempty(opts.YLim) + ylim(opts.YLim); +end + +end diff --git a/Data-Analyzer/+Plotter/plotAverageSpectra.m b/Data-Analyzer/+Plotter/plotSpectralAverages.m similarity index 65% rename from Data-Analyzer/+Plotter/plotAverageSpectra.m rename to Data-Analyzer/+Plotter/plotSpectralAverages.m index 2515794..eb3effa 100644 --- a/Data-Analyzer/+Plotter/plotAverageSpectra.m +++ b/Data-Analyzer/+Plotter/plotSpectralAverages.m @@ -1,5 +1,5 @@ -function plotAverageSpectra(scan_parameter_values, spectral_analysis_results, varargin) -%% plotAverageSpectra +function plotSpectralAverages(scan_parameter_values, results_full, varargin) +%% plotSpectralAverages % Author: Karthik % Date: 2025-09-12 % Version: 1.0 @@ -13,7 +13,7 @@ function plotAverageSpectra(scan_parameter_values, spectral_analysis_results, va % kx, ky, PS_all, k_rho_vals, S_k_all, theta_vals, S_theta_all % % Name-Value Pair Arguments: -% 'ScanParameterName', 'FigNum', 'ColormapPS', 'Font', +% 'FigNum', 'ColormapPS', 'Font', % 'SaveFileName', 'SaveDirectory', 'SkipSaveFigures' % % Notes: @@ -23,78 +23,69 @@ function plotAverageSpectra(scan_parameter_values, spectral_analysis_results, va % Dataset colors are distinct % --- 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; + kx = results_full.kx; + ky = results_full.ky; + ps_list = results_full.PS_all; + k_rho_vals = results_full.k_rho_vals; + s_k_list = results_full.S_k_all; + theta_vals = results_full.theta_vals; + s_theta_list = results_full.S_theta_all; % --- Parse optional parameters --- p = inputParser; - addParameter(p, 'ScanParameterName', 'ScanParameter', @ischar); + addParameter(p, 'Title', 'Average Curves', @(x) ischar(x) || isstring(x)); + addParameter(p, 'AnnotationSuffix', '%.1f^\\circ', @(x) ischar(x) || isstring(x)); + addParameter(p, 'ColormapPS', [], @(x) isnumeric(x) || ismatrix(x)); + addParameter(p, 'FontName', 'Arial', @ischar); + addParameter(p, 'FontSize', 14, @isnumeric); 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', 'avgspectra.fig', @ischar); addParameter(p, 'SaveDirectory', pwd, @ischar); addParameter(p, 'SkipSaveFigures', false, @islogical); parse(p, varargin{:}); opts = p.Results; + + axisFontSize = opts.FontSize; + titleFontSize = opts.FontSize+2; + + % --- Compute averaged data using shared helper --- + PS_data = Analyzer.computeSpectralAverages(ps_list, scan_parameter_values); + S_k_data = Analyzer.computeSpectralAverages(s_k_list, scan_parameter_values); + S_theta_data = Analyzer.computeSpectralAverages(s_theta_list, scan_parameter_values); % --- Unique scan parameters --- - [uniqueParams, ~, idx] = unique(scan_parameter_values); + uniqueParams = S_k_data.scan_parameters; 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; + avgPS = PS_data.curves_mean{pIdx}; + avgS_k = S_k_data.curves_mean{pIdx}; + avgS_theta = S_theta_data.curves_mean{pIdx}; % ==== Plot ==== fig = figure(opts.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',opts.Font); - ylabel('k_y [\mum^{-1}]','Interpreter','tex','FontSize',axisFontSize,'FontName',opts.Font); + xlabel('k_x [\mum^{-1}]','Interpreter','tex','FontSize',axisFontSize,'FontName',opts.FontName); + ylabel('k_y [\mum^{-1}]','Interpreter','tex','FontSize',axisFontSize,'FontName',opts.FontName); title('Average Power Spectrum','FontSize',titleFontSize,'FontWeight','bold'); + if isempty(opts.ColormapPS) + opts.ColormapPS = Colormaps.coolwarm(); + end colormap(opts.ColormapPS); colorbar; % --- Annotate scan parameter --- - if strcmp(opts.ScanParameterName,'ps_rot_mag_fin_pol_angle') - txt = sprintf('%.1f^\\circ', currentParam); - else - txt = sprintf('%.2f G', currentParam); - end + txt = sprintf(opts.AnnotationSuffix, currentParam); text(0.975,0.975,txt,'Color','white','FontWeight','bold','FontSize',axisFontSize, ... 'Interpreter','tex','Units','normalized','HorizontalAlignment','right','VerticalAlignment','top'); @@ -118,11 +109,13 @@ function plotAverageSpectra(scan_parameter_values, spectral_analysis_results, va ax.XMinorGrid = 'on'; ax.YMinorGrid = 'on'; grid on; + + sgtitle(opts.Title, 'FontName', opts.FontName, 'FontSize', titleFontSize+2); drawnow; % --- Save figure --- - saveFigure(fig, ... + Plotter.saveFigure(fig, ... 'SaveFileName', opts.SaveFileName, ... 'SaveDirectory', opts.SaveDirectory, ... 'SkipSaveFigures', opts.SkipSaveFigures); diff --git a/Data-Analyzer/+Plotter/plotSpectralCurves.m b/Data-Analyzer/+Plotter/plotSpectralCurves.m index 8f7c57a..fdf6b32 100644 --- a/Data-Analyzer/+Plotter/plotSpectralCurves.m +++ b/Data-Analyzer/+Plotter/plotSpectralCurves.m @@ -1,4 +1,4 @@ -function results = plotSpectralCurves(S_theta_all, theta_vals, scan_reference_values, varargin) +function results = plotSpectralCurves(S_all, x_vals, scan_reference_values, varargin) %% plotSpectralCurves % Author: Karthik % Date: 2025-10-01 @@ -33,6 +33,9 @@ function results = plotSpectralCurves(S_theta_all, theta_vals, scan_reference_va addParameter(p, 'Title', 'Spectral Curves', @(x) ischar(x) || isstring(x)); addParameter(p, 'XLabel', '\theta / \pi', @(x) ischar(x) || isstring(x)); addParameter(p, 'YLabel', 'S(\theta)', @(x) ischar(x) || isstring(x)); + addParameter(p, 'YScale', 'linear', @(x) ismember(x, {'linear','log'})); + addParameter(p, 'XLim', [], @(x) isempty(x) || (isnumeric(x) && numel(x)==2)); + addParameter(p, 'YLim', [], @(x) isempty(x) || (isnumeric(x) && numel(x)==2)); addParameter(p, 'FontName', 'Arial', @ischar); addParameter(p, 'FontSize', 14, @isnumeric); addParameter(p, 'FigNum', [], @(x) isempty(x) || (isnumeric(x) && isscalar(x))); @@ -45,31 +48,19 @@ function results = plotSpectralCurves(S_theta_all, theta_vals, scan_reference_va parse(p, varargin{:}); opts = p.Results; - % --- Determine number of scan parameters and repetitions --- - N_params = numel(scan_reference_values); - Ntotal = numel(S_theta_all); - Nreps = Ntotal / N_params; + % --- Compute averages using computeSpectralAverages --- + avgData = Analyzer.computeSpectralAverages(S_all, scan_reference_values); - % --- Prepare curves, mean, SEM --- - curves_all = cell(1, N_params); - curves_mean = cell(1, N_params); - curves_err = cell(1, N_params); - Npoints = numel(S_theta_all{1}); - - for i = 1:N_params - curves = zeros(Nreps, Npoints); - for r = 1:Nreps - idx = (r-1)*N_params + i; % correct interleaved indexing - curves(r,:) = S_theta_all{idx}; - end - curves_all{i} = curves; - curves_mean{i} = mean(curves,1); - curves_err{i} = std(curves,0,1)/sqrt(Nreps); - end + % --- Extract averaged curves --- + curves_all = avgData.curves_all; + curves_mean = avgData.curves_mean; + curves_err = avgData.curves_error; + N_params = numel(avgData.scan_parameters); + Npoints = numel(S_all{1}); % --- Create results struct compatible with plotting --- results.curves = curves_all; - results.x_values = theta_vals; % generic x-axis + results.x_values = x_vals; % generic x-axis results.curves_mean = curves_mean; results.curves_error = curves_err; results.scan_parameter_values = scan_reference_values; @@ -111,7 +102,9 @@ function results = plotSpectralCurves(S_theta_all, theta_vals, scan_reference_va ylabel(ax, opts.YLabel, 'FontName', opts.FontName, 'FontSize', opts.FontSize); title(ax, sprintf('%s=%.3g%s', opts.TileTitlePrefix, results.scan_parameter_values(i), opts.TileTitleSuffix), ... 'FontName', opts.FontName, 'FontSize', opts.FontSize); - set(ax, 'FontName', opts.FontName, 'FontSize', opts.FontSize); + set(ax, 'FontName', opts.FontName, 'FontSize', opts.FontSize, 'YScale', opts.YScale); + if ~isempty(opts.XLim), xlim(ax, opts.XLim); end + if ~isempty(opts.YLim), ylim(ax, opts.YLim); end end % --- Save figure --- diff --git a/Data-Analyzer/+Scripts/BECToDropletsToStripes/plotAnalysisResults.m b/Data-Analyzer/+Scripts/BECToDropletsToStripes/plotAnalysisResults.m index 967114f..6d432fa 100644 --- a/Data-Analyzer/+Scripts/BECToDropletsToStripes/plotAnalysisResults.m +++ b/Data-Analyzer/+Scripts/BECToDropletsToStripes/plotAnalysisResults.m @@ -187,16 +187,15 @@ Plotter.plotMultiplePCAResults(compiled_results.pca_results, scan_parameter_valu 'SkipSaveFigures', options.skipSaveFigures, ... 'SaveDirectory', figSaveDir); -%{ %% ------------------ 14. Average of Spectra Plots ------------------ - -Plotter.plotAverageSpectra(scan_parameter_values, ... - spectral_analysis_results, ... - 'ScanParameterName', scan_parameter, ... +Plotter.plotSpectralAverages(scan_parameter_values, ... + compiled_results.spectral_analysis_results, ... + 'Title', options.titleString, ... + 'AnnotationSuffix', '%.1f^\\circ', ... 'FigNum', 20, ... 'ColormapPS', Colormaps.coolwarm(), ... - 'Font', 'Bahnschrift', ... + 'FontName', 'Bahnschrift', ... 'SaveFileName', 'avgSpectra.fig', ... - 'SaveDirectory', figSaveDir, ... - 'SkipSaveFigures', options.skipSaveFigures); -%} \ No newline at end of file + 'SaveDirectory', pwd, ... + 'SkipSaveFigures', true); + diff --git a/Data-Analyzer/+Scripts/BECToDropletsToStripes/runSpectralDistributionAnalysis.m b/Data-Analyzer/+Scripts/BECToDropletsToStripes/runAngularSpectralDistributionAnalysis.m similarity index 97% rename from Data-Analyzer/+Scripts/BECToDropletsToStripes/runSpectralDistributionAnalysis.m rename to Data-Analyzer/+Scripts/BECToDropletsToStripes/runAngularSpectralDistributionAnalysis.m index 87545ba..d4d1d83 100644 --- a/Data-Analyzer/+Scripts/BECToDropletsToStripes/runSpectralDistributionAnalysis.m +++ b/Data-Analyzer/+Scripts/BECToDropletsToStripes/runAngularSpectralDistributionAnalysis.m @@ -42,13 +42,6 @@ options.Angular_Sigma = 2; options.Angular_WindowSize = 5; options.zoom_size = 50; -% -options.maximumShift = 8; -options.Radial_Theta = deg2rad(45); -options.Radial_Minimum = 2; -options.Radial_Maximum = 6; -options.skipLivePlot = false; - % Flags options.skipUnshuffling = false; options.skipNormalization = false; @@ -165,7 +158,7 @@ Plotter.plotSpectralDistributionCumulants(results, ... 'SaveFileName', 'SpectralCumulants.fig'); %% ------------------ 4. Fit shifted Angular Spectral Distribution Curves ------------------ -[fitResults, rawCurves] = Analyzer.fitTwoGaussianCurves(... +[fitResults, rawCurves] = Analyzer.fitTwoGaussianCurvesToAngularSpectralDistribution(... spectral_analysis_results.S_theta_norm_all, ... spectral_analysis_results.theta_vals, ... 'MaxTheta', pi/2, ... @@ -182,7 +175,7 @@ function plotTwoGaussianFitsOnRaw(fitResults, rawCurves, nRows, nCols) % plotTwoGaussianFitsOnRaw - Plots raw curves and overlays valid two-Gaussian fits. % % Inputs: -% fitResults - struct array from fitTwoGaussianCurves (may contain fits) +% fitResults - struct array from fitTwoGaussianCurvesToAngularSpectralDistribution (may contain fits) % rawCurves - struct array with fields: % .x - raw curve % .theta - corresponding theta values @@ -311,7 +304,7 @@ for i = 1:numel(results.curves) end % --- Fit two-Gaussian model to mean curves --- -[fitResultsMean, ~] = Analyzer.fitTwoGaussianCurves(... +[fitResultsMean, ~] = Analyzer.fitTwoGaussianCurvesToAngularSpectralDistribution(... results.curves_mean, ... results.x_values*pi, ... 'MaxTheta', pi/2, ... @@ -394,7 +387,7 @@ curves_cell = num2cell(curves_matrix, 2); % 1x60 cell array curves_cell = cellfun(@(x) x(:).', curves_cell, 'UniformOutput', false); % ensure 1xN row vectors % --- Fit two-Gaussian model to these curves --- -[fitResults, rawCurves] = Analyzer.fitTwoGaussianCurves(... +[fitResults, rawCurves] = Analyzer.fitTwoGaussianCurvesToAngularSpectralDistribution(... curves_cell, ... results.x_values*pi, ... 'MaxTheta', pi/2, ... diff --git a/Data-Analyzer/+Scripts/BECToDropletsToStripes/runCorrelationAnalysis.m b/Data-Analyzer/+Scripts/BECToDropletsToStripes/runCorrelationAnalysis.m index 092907f..c27e3d2 100644 --- a/Data-Analyzer/+Scripts/BECToDropletsToStripes/runCorrelationAnalysis.m +++ b/Data-Analyzer/+Scripts/BECToDropletsToStripes/runCorrelationAnalysis.m @@ -42,12 +42,11 @@ options.Angular_Sigma = 2; options.Angular_WindowSize = 5; options.zoom_size = 50; -% +% Real-Space g2 extraction and analysis options.maximumShift = 8; options.Radial_Theta = deg2rad(45); options.Radial_Minimum = 2; options.Radial_Maximum = 6; -options.skipLivePlot = false; % Flags options.skipUnshuffling = false; @@ -117,7 +116,7 @@ analysis_results = Analyzer.analyzeG2Structures(g2_analysis_results, o %% Plot raw OD images and the corresponding real space g2 matrix options.skipLivePlot = true; -options.skipSaveFigures = false; +options.skipSaveFigures = true; saveDirectory = 'C:\Users\Karthik-OfficePC\Documents\GitRepositories\Calculations\Data-Analyzer\+Scripts'; Plotter.plotODG2withAnalysis(od_imgs, scan_parameter_values, g2_analysis_results, analysis_results, options, ... diff --git a/Data-Analyzer/+Scripts/BECToDropletsToStripes/runRadialSpectralDistributionAnalysis.m b/Data-Analyzer/+Scripts/BECToDropletsToStripes/runRadialSpectralDistributionAnalysis.m new file mode 100644 index 0000000..f637d95 --- /dev/null +++ b/Data-Analyzer/+Scripts/BECToDropletsToStripes/runRadialSpectralDistributionAnalysis.m @@ -0,0 +1,257 @@ +%% ===== BEC-Droplets-Stripes Settings ===== + +% Specify data location to run analysis on +dataSources = { + struct('sequence', 'TwoDGas', ... + 'date', '2025/06/23', ... + 'runs', [300]) % specify run numbers as a string in "" or just as a numeric value +}; + +options = struct(); + +% File paths +options.baseDataFolder = '//DyLabNAS/Data'; +options.FullODImagesFolder = 'E:/Data - Experiment/FullODImages/202506'; +options.measurementName = 'DropletsToStripes'; +scriptFullPath = mfilename('fullpath'); +options.saveDirectory = fileparts(scriptFullPath); + +% Camera / imaging settings +options.cam = 4; % 1 - ODT_1_Axis_Camera; 2 - ODT_2_Axis_Camera; 3 - Horizontal_Axis_Camera;, 4 - Vertical_Axis_Camera; +options.angle = 0; % angle by which image will be rotated +options.center = [1410, 2030]; +options.span = [200, 200]; +options.fraction = [0.1, 0.1]; +options.pixel_size = 5.86e-6; % in meters +options.magnification = 23.94; +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 = 150; +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 = 360; +options.Angular_Threshold = 75; +options.Angular_Sigma = 2; +options.Angular_WindowSize = 5; +options.zoom_size = 50; + +% Flags +options.skipUnshuffling = false; +options.skipNormalization = false; + +options.skipFringeRemoval = true; +options.skipPreprocessing = true; +options.skipMasking = true; +options.skipIntensityThresholding = true; +options.skipBinarization = true; + +options.skipFullODImagesFolderUse = false; +options.skipSaveData = false; +options.skipSaveFigures = true; +options.skipSaveProcessedOD = true; +options.skipLivePlot = false; +options.showProgressBar = true; + +% Extras +options.font = 'Bahnschrift'; +switch options.measurementName + case 'BECToDroplets' + options.scan_parameter = 'rot_mag_field'; + options.flipSortOrder = true; + options.scanParameterUnits = 'gauss'; + options.titleString = 'BEC to Droplets'; + case 'BECToStripes' + options.scan_parameter = 'rot_mag_field'; + options.flipSortOrder = true; + options.scanParameterUnits = 'gauss'; + options.titleString = 'BEC to Stripes'; + case 'DropletsToStripes' + options.scan_parameter = 'ps_rot_mag_fin_pol_angle'; + options.flipSortOrder = false; + options.scanParameterUnits = 'degrees'; + options.titleString = 'Droplets to Stripes'; + case 'StripesToDroplets' + options.scan_parameter = 'ps_rot_mag_fin_pol_angle'; + options.flipSortOrder = false; + options.scanParameterUnits = 'degrees'; + options.titleString = 'Stripes to Droplets'; +end + +%% ===== Collect Images and Launch Viewer ===== + +[options.selectedPath, options.folderPath] = Helper.selectDataSourcePath(dataSources, options); + +[od_imgs, scan_parameter_values, scan_reference_values, file_list] = Helper.collectODImages(options); + +%% Conduct spectral analysis +spectral_analysis_results = Analyzer.extractFullRadialSpectralDistribution(od_imgs, options); + +%% ------------------ 1. Plot of all Radial Spectral Distribution Curves ------------------ +Plotter.plotSpectralCurves( ... + spectral_analysis_results.S_k_all, ... + spectral_analysis_results.k_rho_vals, ... + scan_reference_values, ... % correct scan params + 'Title', options.titleString, ... + 'XLabel', 'k_\rho', ... + 'YLabel', 'S(k_\rho)', ... + 'YScale', 'log', ... + 'XLim', [min(spectral_analysis_results.k_rho_vals), max(spectral_analysis_results.k_rho_vals)], ... + 'HighlightEvery', 10, ... % highlight every 10th repetition + 'FontName', options.font, ... + 'FigNum', 1, ... + 'TileTitlePrefix', '\alpha', ... % user-defined tile prefix + 'TileTitleSuffix', '^\circ', ... % user-defined suffix (e.g., degrees symbol) + 'SkipSaveFigures', options.skipSaveFigures, ... + 'SaveFileName', 'SpectralCurves.fig', ... + 'SaveDirectory', options.saveDirectory); + +%% ------------------ 2. Fit truncated Radial Spectral Distribution Curves ------------------ +[fitResults, rawCurves] = Analyzer.fitTwoGaussianCurvesToRadialSpectralDistribution(... + spectral_analysis_results.S_k_all, ... % radial spectral curves + spectral_analysis_results.k_rho_vals, ... % corresponding k_rho values + 'KRhoRange', [0, 3], ... % truncate curves to 0 <= k_rho <= 3 + 'AmplitudeRange', [0, 0.5], ... % truncate curves to y >= 0 (all amplitudes) + 'ResidualThreshold', 0.15, ... % maximum allowed residual + 'PositionThreshold', 0.5, ... % minimum separation between peaks + 'AmplitudeThreshold', 0.015); % minimum secondary peak fraction +%% +%{ +% --- Function call --- +plotTwoGaussianFitsOnRaw(fitResults, rawCurves, 8, 12); % 8×12 subplots per page + +% --- Function definition --- +function plotTwoGaussianFitsOnRaw(fitResults, rawCurves, nRows, nCols) +%% plotTwoGaussianFitsOnRaw +% Plots raw radial spectral curves with their two-Gaussian fits. +% +% Inputs: +% fitResults - struct array from fitTwoGaussianCurvesToRadialSpectralDistribution +% rawCurves - struct array with fields: +% .x - raw normalized amplitudes +% .k - corresponding k_rho values +% nRows - number of subplot rows per page (default: 8) +% nCols - number of subplot columns per page (default: 12) + + if nargin < 3, nRows = 8; end + if nargin < 4, nCols = 12; end + + Ncurves = numel(rawCurves); + plotsPerPage = nRows * nCols; + pageNum = 1; + + for k = 1:Ncurves + % --- Create new figure page if needed --- + if mod(k-1, plotsPerPage) == 0 + figure('Name', sprintf('Radial Spectra Page %d', pageNum), ... + 'NumberTitle', 'off', 'Color', 'w', ... + 'Position', [50 50 1600 900]); + pageNum = pageNum + 1; + end + + % --- Subplot selection --- + subplot(nRows, nCols, mod(k-1, plotsPerPage) + 1); + hold on; grid on; box on; + + % --- Plot raw curve --- + if isfield(rawCurves, 'x') && isfield(rawCurves, 'k') && ... + ~isempty(rawCurves(k).x) && all(isfinite(rawCurves(k).x)) + plot(rawCurves(k).k, rawCurves(k).x, 'k.-', ... + 'LineWidth', 1, 'DisplayName', 'Raw data'); + end + + % --- Overlay fit if valid --- + if k <= numel(fitResults) + fit = fitResults(k); + if isfield(fit, 'isValid') && fit.isValid && ... + isfield(fit, 'kFine') && isfield(fit, 'yFit') && ... + ~isempty(fit.kFine) && all(isfinite(fit.yFit)) + plot(fit.kFine, fit.yFit, 'r-', ... + 'LineWidth', 1.2, 'DisplayName', 'Two-Gaussian fit'); + end + end + + % --- Labels and title --- + xlabel('k_\rho', 'FontSize', 10); + ylabel('Normalized amplitude', 'FontSize', 10); + title(sprintf('Curve %d', k), 'FontSize', 9, 'Interpreter', 'none'); + + % --- Legend on first subplot of each page --- + if mod(k-1, plotsPerPage) == 0 + legend('Location', 'best', 'FontSize', 8); + end + + hold off; + end +end +%} +%% ------------------ 3. Plot fit parameters - amplitude ------------------ +Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'A2', ... + 'Title', options.titleString, ... + 'XLabel', '\alpha (degrees)', ... + 'YLabel', 'Secondary peak amplitude', ... + 'FontName', options.font, ... + 'FontSize', 16, ... + 'FigNum', 4, ... + 'SkipSaveFigures', options.skipSaveFigures, ... + 'SaveFileName', 'SecondaryPeakAmplitudePDF.fig', ... + 'PlotType', 'histogram', ... + 'NumberOfBins', 20, ... + 'NormalizeHistogram', true, ... + 'Colormap', @Colormaps.coolwarm, ... + 'XLim', [min(scan_reference_values), max(scan_reference_values)], ... + 'YLim', [0, 0.06]); + +%% ------------------ 4. Plot fit parameters - position ------------------ +Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'mu2', ... + 'Title', options.titleString, ... + 'XLabel', '\alpha (degrees)', ... + 'YLabel', 'Secondary peak position (\theta, rad)', ... + 'FontName', options.font, ... + 'FontSize', 16, ... + 'FigNum', 5, ... + 'SkipSaveFigures', options.skipSaveFigures, ... + 'SaveFileName', 'SecondaryPeakPositionPDF.fig', ... + 'PlotType', 'histogram', ... + 'NumberOfBins', 20, ... + 'NormalizeHistogram', true, ... + 'Colormap', @Colormaps.coolwarm, ... + 'XLim', [min(scan_reference_values), max(scan_reference_values)], ... + 'YLim', [1, 2]); +%% ------------------ 5. Plot fit parameters - width ------------------ +Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'sigma2', ... + 'Title', options.titleString, ... + 'XLabel', '\alpha (degrees)', ... + 'YLabel', 'Secondary peak width (\sigma, rad)', ... + 'FontName', options.font, ... + 'FontSize', 16, ... + 'FigNum', 6, ... + 'SkipSaveFigures', options.skipSaveFigures, ... + 'SaveFileName', 'SecondaryPeakWidthPDF.fig', ... + 'PlotType', 'histogram', ... + 'NumberOfBins', 20, ... + 'NormalizeHistogram', true, ... + 'Colormap', @Colormaps.coolwarm, ... + 'XLim', [min(scan_reference_values), max(scan_reference_values)], ... + 'YLim', [0, 1.5]); + +%% ------------------ 6. Plot secondary Gaussian range ------------------ +Plotter.plotSecondaryGaussianRange(fitResults, scan_reference_values, ... + 'Title', options.titleString, ... + 'XLabel', '\alpha (degrees)', ... + 'YLabel', '\mu_2 \pm \sigma_2 (rad)', ... + 'FontName', options.font, ... + 'FontSize', 16, ... + 'FigNum', 8, ... + 'NumberOfBins', 20, ... + 'NormalizeHistogram', true, ... + 'Colormap', @Colormaps.coolwarm, ... + 'OverlayMeanSEM', true); + + diff --git a/Data-Analyzer/+Scripts/BECToStripesToDroplets/plotAnalysisResults.m b/Data-Analyzer/+Scripts/BECToStripesToDroplets/plotAnalysisResults.m index bda2f3a..bb7fbef 100644 --- a/Data-Analyzer/+Scripts/BECToStripesToDroplets/plotAnalysisResults.m +++ b/Data-Analyzer/+Scripts/BECToStripesToDroplets/plotAnalysisResults.m @@ -187,16 +187,14 @@ Plotter.plotMultiplePCAResults(compiled_results.pca_results, scan_parameter_valu 'SkipSaveFigures', options.skipSaveFigures, ... 'SaveDirectory', figSaveDir); -%{ %% ------------------ 14. Average of Spectra Plots ------------------ - -Plotter.plotAverageSpectra(scan_parameter_values, ... - spectral_analysis_results, ... - 'ScanParameterName', scan_parameter, ... +Plotter.plotSpectralAverages(scan_parameter_values, ... + compiled_results.spectral_analysis_results, ... + 'Title', options.titleString, ... + 'AnnotationSuffix', '%.1f^\\circ', ... 'FigNum', 20, ... 'ColormapPS', Colormaps.coolwarm(), ... - 'Font', 'Bahnschrift', ... + 'FontName', 'Bahnschrift', ... 'SaveFileName', 'avgSpectra.fig', ... - 'SaveDirectory', figSaveDir, ... - 'SkipSaveFigures', options.skipSaveFigures); -%} \ No newline at end of file + 'SaveDirectory', pwd, ... + 'SkipSaveFigures', true); \ No newline at end of file diff --git a/Data-Analyzer/+Scripts/BECToStripesToDroplets/runSpectralDistributionAnalysis.m b/Data-Analyzer/+Scripts/BECToStripesToDroplets/runAngularSpectralDistributionAnalysis.m similarity index 97% rename from Data-Analyzer/+Scripts/BECToStripesToDroplets/runSpectralDistributionAnalysis.m rename to Data-Analyzer/+Scripts/BECToStripesToDroplets/runAngularSpectralDistributionAnalysis.m index 75aba95..09c9e68 100644 --- a/Data-Analyzer/+Scripts/BECToStripesToDroplets/runSpectralDistributionAnalysis.m +++ b/Data-Analyzer/+Scripts/BECToStripesToDroplets/runAngularSpectralDistributionAnalysis.m @@ -158,7 +158,7 @@ Plotter.plotSpectralDistributionCumulants(results, ... 'SaveFileName', 'SpectralCumulants.fig'); %% ------------------ 4. Fit shifted Angular Spectral Distribution Curves ------------------ -[fitResults, rawCurves] = Analyzer.fitTwoGaussianCurves(... +[fitResults, rawCurves] = Analyzer.fitTwoGaussianCurvesToAngularSpectralDistribution(... spectral_analysis_results.S_theta_norm_all, ... spectral_analysis_results.theta_vals, ... 'MaxTheta', pi/2, ... @@ -175,7 +175,7 @@ function plotTwoGaussianFitsOnRaw(fitResults, rawCurves, nRows, nCols) % plotTwoGaussianFitsOnRaw - Plots raw curves and overlays valid two-Gaussian fits. % % Inputs: -% fitResults - struct array from fitTwoGaussianCurves (may contain fits) +% fitResults - struct array from fitTwoGaussianCurvesToAngularSpectralDistribution (may contain fits) % rawCurves - struct array with fields: % .x - raw curve % .theta - corresponding theta values @@ -304,7 +304,7 @@ for i = 1:numel(results.curves) end % --- Fit two-Gaussian model to mean curves --- -[fitResultsMean, ~] = Analyzer.fitTwoGaussianCurves(... +[fitResultsMean, ~] = Analyzer.fitTwoGaussianCurvesToAngularSpectralDistribution(... results.curves_mean, ... results.x_values*pi, ... 'MaxTheta', pi/2, ... @@ -387,7 +387,7 @@ curves_cell = num2cell(curves_matrix, 2); % 1x60 cell array curves_cell = cellfun(@(x) x(:).', curves_cell, 'UniformOutput', false); % ensure 1xN row vectors % --- Fit two-Gaussian model to these curves --- -[fitResults, rawCurves] = Analyzer.fitTwoGaussianCurves(... +[fitResults, rawCurves] = Analyzer.fitTwoGaussianCurvesToAngularSpectralDistribution(... curves_cell, ... results.x_values*pi, ... 'MaxTheta', pi/2, ... diff --git a/Data-Analyzer/+Scripts/BECToStripesToDroplets/runCorrelationAnalysis.m b/Data-Analyzer/+Scripts/BECToStripesToDroplets/runCorrelationAnalysis.m index 9600754..d447b5b 100644 --- a/Data-Analyzer/+Scripts/BECToStripesToDroplets/runCorrelationAnalysis.m +++ b/Data-Analyzer/+Scripts/BECToStripesToDroplets/runCorrelationAnalysis.m @@ -42,12 +42,11 @@ options.Angular_Sigma = 2; options.Angular_WindowSize = 5; options.zoom_size = 50; -% +% Real-Space g2 extraction and analysis options.maximumShift = 8; options.Radial_Theta = deg2rad(45); options.Radial_Minimum = 2; options.Radial_Maximum = 6; -options.skipLivePlot = false; % Flags options.skipUnshuffling = false; diff --git a/Data-Analyzer/+Scripts/BECToStripesToDroplets/runRadialSpectralDistributionAnalysis.m b/Data-Analyzer/+Scripts/BECToStripesToDroplets/runRadialSpectralDistributionAnalysis.m new file mode 100644 index 0000000..4caf1fe --- /dev/null +++ b/Data-Analyzer/+Scripts/BECToStripesToDroplets/runRadialSpectralDistributionAnalysis.m @@ -0,0 +1,257 @@ +%% ===== BEC-Droplets-Stripes Settings ===== + +% Specify data location to run analysis on +dataSources = { + struct('sequence', 'TwoDGas', ... + 'date', '2025/06/23', ... + 'runs', [300]) % specify run numbers as a string in "" or just as a numeric value +}; + +options = struct(); + +% File paths +options.baseDataFolder = '//DyLabNAS/Data'; +options.FullODImagesFolder = 'E:/Data - Experiment/FullODImages/202506'; +options.measurementName = 'DropletsToStripes'; +scriptFullPath = mfilename('fullpath'); +options.saveDirectory = fileparts(scriptFullPath); + +% Camera / imaging settings +options.cam = 4; % 1 - ODT_1_Axis_Camera; 2 - ODT_2_Axis_Camera; 3 - Horizontal_Axis_Camera;, 4 - Vertical_Axis_Camera; +options.angle = 0; % angle by which image will be rotated +options.center = [1410, 2030]; +options.span = [200, 200]; +options.fraction = [0.1, 0.1]; +options.pixel_size = 5.86e-6; % in meters +options.magnification = 23.94; +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 = 360; +options.Angular_Threshold = 75; +options.Angular_Sigma = 2; +options.Angular_WindowSize = 5; +options.zoom_size = 50; + +% Flags +options.skipUnshuffling = false; +options.skipNormalization = false; + +options.skipFringeRemoval = true; +options.skipPreprocessing = true; +options.skipMasking = true; +options.skipIntensityThresholding = true; +options.skipBinarization = true; + +options.skipFullODImagesFolderUse = false; +options.skipSaveData = false; +options.skipSaveFigures = true; +options.skipSaveProcessedOD = true; +options.skipLivePlot = false; +options.showProgressBar = true; + +% Extras +options.font = 'Bahnschrift'; +switch options.measurementName + case 'BECToDroplets' + options.scan_parameter = 'rot_mag_field'; + options.flipSortOrder = true; + options.scanParameterUnits = 'gauss'; + options.titleString = 'BEC to Droplets'; + case 'BECToStripes' + options.scan_parameter = 'rot_mag_field'; + options.flipSortOrder = true; + options.scanParameterUnits = 'gauss'; + options.titleString = 'BEC to Stripes'; + case 'DropletsToStripes' + options.scan_parameter = 'ps_rot_mag_fin_pol_angle'; + options.flipSortOrder = false; + options.scanParameterUnits = 'degrees'; + options.titleString = 'Droplets to Stripes'; + case 'StripesToDroplets' + options.scan_parameter = 'ps_rot_mag_fin_pol_angle'; + options.flipSortOrder = false; + options.scanParameterUnits = 'degrees'; + options.titleString = 'Stripes to Droplets'; +end + +%% ===== Collect Images and Launch Viewer ===== + +[options.selectedPath, options.folderPath] = Helper.selectDataSourcePath(dataSources, options); + +[od_imgs, scan_parameter_values, scan_reference_values, file_list] = Helper.collectODImages(options); + +%% Conduct spectral analysis +spectral_analysis_results = Analyzer.extractFullRadialSpectralDistribution(od_imgs, options); + +%% ------------------ 1. Plot of all Radial Spectral Distribution Curves ------------------ +Plotter.plotSpectralCurves( ... + spectral_analysis_results.S_k_all, ... + spectral_analysis_results.k_rho_vals, ... + scan_reference_values, ... % correct scan params + 'Title', options.titleString, ... + 'XLabel', 'k_\rho', ... + 'YLabel', 'S(k_\rho)', ... + 'YScale', 'log', ... + 'XLim', [min(spectral_analysis_results.k_rho_vals), max(spectral_analysis_results.k_rho_vals)], ... + 'HighlightEvery', 10, ... % highlight every 10th repetition + 'FontName', options.font, ... + 'FigNum', 1, ... + 'TileTitlePrefix', '\alpha', ... % user-defined tile prefix + 'TileTitleSuffix', '^\circ', ... % user-defined suffix (e.g., degrees symbol) + 'SkipSaveFigures', options.skipSaveFigures, ... + 'SaveFileName', 'SpectralCurves.fig', ... + 'SaveDirectory', options.saveDirectory); + +%% ------------------ 2. Fit truncated Radial Spectral Distribution Curves ------------------ +[fitResults, rawCurves] = Analyzer.fitTwoGaussianCurvesToRadialSpectralDistribution(... + spectral_analysis_results.S_k_all, ... % radial spectral curves + spectral_analysis_results.k_rho_vals, ... % corresponding k_rho values + 'KRhoRange', [0, 3], ... % truncate curves to 0 <= k_rho <= 3 + 'AmplitudeRange', [0, 0.5], ... % truncate curves to y >= 0 (all amplitudes) + 'ResidualThreshold', 0.15, ... % maximum allowed residual + 'PositionThreshold', 0.5, ... % minimum separation between peaks + 'AmplitudeThreshold', 0.015); % minimum secondary peak fraction +%% +%{ +% --- Function call --- +plotTwoGaussianFitsOnRaw(fitResults, rawCurves, 8, 12); % 8×12 subplots per page + +% --- Function definition --- +function plotTwoGaussianFitsOnRaw(fitResults, rawCurves, nRows, nCols) +%% plotTwoGaussianFitsOnRaw +% Plots raw radial spectral curves with their two-Gaussian fits. +% +% Inputs: +% fitResults - struct array from fitTwoGaussianCurvesToRadialSpectralDistribution +% rawCurves - struct array with fields: +% .x - raw normalized amplitudes +% .k - corresponding k_rho values +% nRows - number of subplot rows per page (default: 8) +% nCols - number of subplot columns per page (default: 12) + + if nargin < 3, nRows = 8; end + if nargin < 4, nCols = 12; end + + Ncurves = numel(rawCurves); + plotsPerPage = nRows * nCols; + pageNum = 1; + + for k = 1:Ncurves + % --- Create new figure page if needed --- + if mod(k-1, plotsPerPage) == 0 + figure('Name', sprintf('Radial Spectra Page %d', pageNum), ... + 'NumberTitle', 'off', 'Color', 'w', ... + 'Position', [50 50 1600 900]); + pageNum = pageNum + 1; + end + + % --- Subplot selection --- + subplot(nRows, nCols, mod(k-1, plotsPerPage) + 1); + hold on; grid on; box on; + + % --- Plot raw curve --- + if isfield(rawCurves, 'x') && isfield(rawCurves, 'k') && ... + ~isempty(rawCurves(k).x) && all(isfinite(rawCurves(k).x)) + plot(rawCurves(k).k, rawCurves(k).x, 'k.-', ... + 'LineWidth', 1, 'DisplayName', 'Raw data'); + end + + % --- Overlay fit if valid --- + if k <= numel(fitResults) + fit = fitResults(k); + if isfield(fit, 'isValid') && fit.isValid && ... + isfield(fit, 'kFine') && isfield(fit, 'yFit') && ... + ~isempty(fit.kFine) && all(isfinite(fit.yFit)) + plot(fit.kFine, fit.yFit, 'r-', ... + 'LineWidth', 1.2, 'DisplayName', 'Two-Gaussian fit'); + end + end + + % --- Labels and title --- + xlabel('k_\rho', 'FontSize', 10); + ylabel('Normalized amplitude', 'FontSize', 10); + title(sprintf('Curve %d', k), 'FontSize', 9, 'Interpreter', 'none'); + + % --- Legend on first subplot of each page --- + if mod(k-1, plotsPerPage) == 0 + legend('Location', 'best', 'FontSize', 8); + end + + hold off; + end +end +%} +%% ------------------ 3. Plot fit parameters - amplitude ------------------ +Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'A2', ... + 'Title', options.titleString, ... + 'XLabel', '\alpha (degrees)', ... + 'YLabel', 'Secondary peak amplitude', ... + 'FontName', options.font, ... + 'FontSize', 16, ... + 'FigNum', 4, ... + 'SkipSaveFigures', options.skipSaveFigures, ... + 'SaveFileName', 'SecondaryPeakAmplitudePDF.fig', ... + 'PlotType', 'histogram', ... + 'NumberOfBins', 20, ... + 'NormalizeHistogram', true, ... + 'Colormap', @Colormaps.coolwarm, ... + 'XLim', [min(scan_reference_values), max(scan_reference_values)], ... + 'YLim', [0, 0.06]); + +%% ------------------ 4. Plot fit parameters - position ------------------ +Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'mu2', ... + 'Title', options.titleString, ... + 'XLabel', '\alpha (degrees)', ... + 'YLabel', 'Secondary peak position (\theta, rad)', ... + 'FontName', options.font, ... + 'FontSize', 16, ... + 'FigNum', 5, ... + 'SkipSaveFigures', options.skipSaveFigures, ... + 'SaveFileName', 'SecondaryPeakPositionPDF.fig', ... + 'PlotType', 'histogram', ... + 'NumberOfBins', 20, ... + 'NormalizeHistogram', true, ... + 'Colormap', @Colormaps.coolwarm, ... + 'XLim', [min(scan_reference_values), max(scan_reference_values)], ... + 'YLim', [1, 2]); +%% ------------------ 5. Plot fit parameters - width ------------------ +Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'sigma2', ... + 'Title', options.titleString, ... + 'XLabel', '\alpha (degrees)', ... + 'YLabel', 'Secondary peak width (\sigma, rad)', ... + 'FontName', options.font, ... + 'FontSize', 16, ... + 'FigNum', 6, ... + 'SkipSaveFigures', options.skipSaveFigures, ... + 'SaveFileName', 'SecondaryPeakWidthPDF.fig', ... + 'PlotType', 'histogram', ... + 'NumberOfBins', 20, ... + 'NormalizeHistogram', true, ... + 'Colormap', @Colormaps.coolwarm, ... + 'XLim', [min(scan_reference_values), max(scan_reference_values)], ... + 'YLim', [0, 1.5]); + +%% ------------------ 6. Plot secondary Gaussian range ------------------ +Plotter.plotSecondaryGaussianRange(fitResults, scan_reference_values, ... + 'Title', options.titleString, ... + 'XLabel', '\alpha (degrees)', ... + 'YLabel', '\mu_2 \pm \sigma_2 (rad)', ... + 'FontName', options.font, ... + 'FontSize', 16, ... + 'FigNum', 8, ... + 'NumberOfBins', 20, ... + 'NormalizeHistogram', true, ... + 'Colormap', @Colormaps.coolwarm, ... + 'OverlayMeanSEM', true); + +