From 249d81ce4d5321e6cb4277b961a7098d0e0ba5ac Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Thu, 9 Oct 2025 20:02:52 +0200 Subject: [PATCH] Major update - correction to analysis of recentered angular spectral distribution curves, corrections to fitting and plotting of fit results. --- .../extractFullAngularSpectralDistribution.m | 99 +++++ .../+Analyzer/fitTwoGaussianCurves.m | 295 ++++++++------- .../+Analyzer/recenterSpectralCurves.m | 88 +++++ .../computeAngularSpectralDistribution.m | 25 +- Data-Analyzer/+Plotter/plotFitParameterPDF.m | 35 +- Data-Analyzer/+Plotter/plotG2Curves.m | 2 +- .../+Plotter/plotSpectralCurvesRecentered.m | 62 +--- .../plotSpectralDistributionCumulants.m | 2 +- .../plotAnalysisResults.m | 240 +----------- .../runSpectralDistributionAnalysis.m | 342 ++++++++++++++++++ .../plotAnalysisResults.m | 243 +------------ 11 files changed, 733 insertions(+), 700 deletions(-) create mode 100644 Data-Analyzer/+Analyzer/extractFullAngularSpectralDistribution.m create mode 100644 Data-Analyzer/+Analyzer/recenterSpectralCurves.m create mode 100644 Data-Analyzer/+Scripts/BECToDropletsToStripes/runSpectralDistributionAnalysis.m diff --git a/Data-Analyzer/+Analyzer/extractFullAngularSpectralDistribution.m b/Data-Analyzer/+Analyzer/extractFullAngularSpectralDistribution.m new file mode 100644 index 0000000..70792f3 --- /dev/null +++ b/Data-Analyzer/+Analyzer/extractFullAngularSpectralDistribution.m @@ -0,0 +1,99 @@ +function results = extractFullAngularSpectralDistribution(od_imgs, options) +%% extractFullAngularSpectralDistribution +% Author: Karthik +% Date: 2025-10-09 +% Version: 1.0 +% +% Description: +% Computes the angular spectral distribution S(θ) from OD images, +% extending θ from 0 → 2π using the computeAngularSpectralDistribution function. +% +% Inputs: +% od_imgs - cell array of OD images +% options - struct containing relevant parameters: +% pixel_size, magnification, zoom_size, +% k_min, k_max, N_angular_bins, +% Angular_Threshold, Angular_Sigma +% +% Outputs: +% results - struct containing θ values and angular spectra +% +% Notes: +% This function is a minimal variant of conductSpectralAnalysis, +% stopping right after angular spectral extraction. + + %% ===== Unpack options ===== + 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; + 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_theta_all = cell(1, N_shots); + S_theta_norm_all = cell(1, N_shots); + theta_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 angular spectrum (0 → 2π) ===== + [theta_vals, S_theta] = Calculator.computeAngularSpectralDistribution( ... + fft_imgs{k}, kx, ky, k_min, k_max, N_angular_bins, ... + Angular_Threshold, Angular_Sigma, [], 2*pi); + + % Normalize angular spectrum and compute weight + S_theta_norm = S_theta / max(S_theta); + + % Store results + S_theta_all{k} = S_theta; + S_theta_norm_all{k} = S_theta_norm; + end + + %% ===== Package results ===== + results = struct(); + results.theta_vals = theta_vals; + results.S_theta_all = S_theta_all; + results.S_theta_norm_all = S_theta_norm_all; +end diff --git a/Data-Analyzer/+Analyzer/fitTwoGaussianCurves.m b/Data-Analyzer/+Analyzer/fitTwoGaussianCurves.m index 1d04a65..8277d41 100644 --- a/Data-Analyzer/+Analyzer/fitTwoGaussianCurves.m +++ b/Data-Analyzer/+Analyzer/fitTwoGaussianCurves.m @@ -5,171 +5,192 @@ function fitResults = fitTwoGaussianCurves(S_theta_all, theta_vals, varargin) % % Author: Karthik % Date: 2025-10-06 -% Version: 1.0 +% Version: 1.1 % % Inputs: % S_theta_all - 1xN cell array of spectral curves % theta_vals - vector of corresponding theta values (radians) % % Optional name-value pairs: -% 'MaxTheta' - Maximum theta for search (default: pi/2) -% 'DeviationLimit' - Max relative deviation allowed (default: 0.3) +% 'MaxTheta' - Maximum theta for search (default: pi/2) +% 'DeviationLimit' - Max relative deviation allowed (default: 0.3) +% 'PositionThreshold' - Minimum angular separation between main and +% secondary peaks (default: 0.1 rad) +% 'AmplitudeThreshold' - Minimum fraction of primary peak amplitude for secondary peak (default: 0.5) % % Output: % fitResults - struct array with fields: -% .pFit - fitted parameters -% .thetaFit - theta values used for fit -% .xFit - curve values used for fit -% .yFit - fitted curve evaluated on thetaFine -% .thetaFine - fine theta for plotting -% .isValid - whether fit passed deviation threshold -% .fitMaxTheta - theta used as fit limit +% .pFit, .thetaFit, .xFit, .yFit, .thetaFine, .isValid, .fitMaxTheta % --- Parse optional inputs --- p = inputParser; addParameter(p, 'MaxTheta', pi/2, @(x) isnumeric(x) && isscalar(x)); addParameter(p, 'DeviationLimit', 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)); + addParameter(p, 'RecenterCurves', true, @islogical); parse(p, varargin{:}); opts = p.Results; Ncurves = numel(S_theta_all); fitResults = struct('pFit',[],'thetaFit',[],'xFit',[],'yFit',[],... 'thetaFine',[],'isValid',[],'fitMaxTheta',[]); - - % --- Preprocess curves: shift first peak to zero without wrapping --- - S_theta_all_shifted = cell(size(S_theta_all)); - for k = 1:Ncurves - curve = S_theta_all{k}; - - % --- Find first peak in 0–MaxTheta range --- - idx_range = find(theta_vals >= 0 & theta_vals <= opts.MaxTheta); - if isempty(idx_range) - [~, peak_idx] = max(curve); - else - [~, local_max_idx] = max(curve(idx_range)); - peak_idx = idx_range(local_max_idx); - end - - % Take only the part from the peak to the end - S_theta_all_shifted{k} = curve(peak_idx:end); - end - - % --- Pad curves to the same length --- - Npoints_shifted = max(cellfun(@numel, S_theta_all_shifted)); - for k = 1:Ncurves - len = numel(S_theta_all_shifted{k}); - if len < Npoints_shifted - S_theta_all_shifted{k} = [S_theta_all_shifted{k}, nan(1, Npoints_shifted-len)]; - end - end - - % --- Define recentered theta values --- - theta_recentered = theta_vals(1:Npoints_shifted) - theta_vals(1); - - for k = 1:Ncurves - x = S_theta_all_shifted{k}; % use recentered curve - theta = theta_recentered; % recentered x-values - - validIdx = ~isnan(x); - x = x(validIdx); - theta = theta(validIdx); - - % --- Restrict to 0–MaxTheta --- - mask = theta>=0 & theta<=opts.MaxTheta; - x = x(mask); - theta = theta(mask); - if numel(theta) < 8 - warning('Curve %d too short (<8 points), skipping.', k); - continue; - end - - % --- Normalize so S(0)=1 --- - x = x / x(1); - - % --- Smooth for stable peak detection --- - xSmooth = smooth(x, 5, 'moving'); - - % --- Find local maxima --- - [pk, locIdx] = findpeaks(xSmooth); - thetaPeaks = theta(locIdx); - - if isempty(pk) - warning('Curve %d has no significant peaks, restricting to within pi/4.', k); - fitMaxTheta = pi/4; - else - % --- Primary peak --- - [mainAmp, mainIdx] = max(pk); - mainTheta = thetaPeaks(mainIdx); - thetaPeaks = thetaPeaks(:); % column vector - pk = pk(:); + if opts.RecenterCurves + % --- Preprocess curves: shift first peak to zero without wrapping --- + S_theta_all_shifted = cell(size(S_theta_all)); + for k = 1:Ncurves + curve = S_theta_all{k}; - % --- Secondary peak ≥50% of primary after main --- - secondaryIdx = find((thetaPeaks>mainTheta) & (pk>=0.50*mainAmp), 1, 'first'); - if isempty(secondaryIdx) - fitMaxTheta = opts.MaxTheta; + idx_range = find(theta_vals >= 0 & theta_vals <= opts.MaxTheta); + if isempty(idx_range) + [~, peak_idx] = max(curve); else - fitMaxTheta = thetaPeaks(secondaryIdx); + [~, local_max_idx] = max(curve(idx_range)); + peak_idx = idx_range(local_max_idx); + end + + S_theta_all_shifted{k} = curve(peak_idx:end); + end + + % --- Pad curves to same length --- + Npoints_shifted = max(cellfun(@numel, S_theta_all_shifted)); + for k = 1:Ncurves + len = numel(S_theta_all_shifted{k}); + if len < Npoints_shifted + S_theta_all_shifted{k} = [S_theta_all_shifted{k}, nan(1, Npoints_shifted-len)]; end end + + theta_recentered = theta_vals(1:Npoints_shifted) - theta_vals(1); + else + S_theta_all_shifted = S_theta_all; + theta_recentered = theta_vals; + end - % --- Extract data up to secondary peak --- - fitMask = theta>=0 & theta<=fitMaxTheta; - thetaFit = theta(fitMask); - xFit = x(fitMask); - if numel(thetaFit) < 8 - warning('Curve %d has too few points for fitting, skipping.', k); - continue; - end - - % --- Two-Gaussian model --- - twoGauss = @(p,theta) ... - 1.0*exp(-0.5*((theta-p(1))/max(p(2),1e-6)).^2) + ... - p(3)*exp(-0.5*((theta-p(4))/max(p(5),1e-6)).^2); - % p = [mu1, sigma1, A2, mu2, sigma2] - - % --- Initial guesses --- - mu1_guess = 0; - sigma1_guess = 0.1*fitMaxTheta; - A2_guess = max(xFit)*0.8; - mu2_guess = fitMaxTheta/2; - sigma2_guess = fitMaxTheta/5; - p0 = [mu1_guess, sigma1_guess, A2_guess, mu2_guess, sigma2_guess]; - - % --- Bounds --- - lb = [0, 1e-6, 0, 0, 1e-6]; - ub = [fitMaxTheta/4, fitMaxTheta, 2, fitMaxTheta, fitMaxTheta]; - - optsLSQ = optimoptions('lsqcurvefit','Display','off', ... - 'MaxFunctionEvaluations',1e4,'MaxIterations',1e4); - - % --- Fit --- + for k = 1:Ncurves try - pFit = lsqcurvefit(twoGauss, p0, thetaFit, xFit, lb, ub, optsLSQ); + x = S_theta_all_shifted{k}; + theta = theta_recentered; + + validIdx = ~isnan(x); + x = x(validIdx); + theta = theta(validIdx); + + mask = theta>=0 & theta<=opts.MaxTheta; + x = x(mask); + theta = theta(mask); + if numel(theta) < 8 + warning('Curve %d too short (<8 points), skipping.', k); + fitResults(k) = fillNaNStruct(); + continue; + end + + x = x / x(1); + xSmooth = smooth(x, 5, 'moving'); + curveAmp = max(xSmooth) - min(xSmooth); + + % --- Robust peak detection --- + [pk, locIdx] = findpeaks(xSmooth, 'MinPeakProminence', opts.AmplitudeThreshold * curveAmp); + thetaPeaks = theta(locIdx); + + if isempty(pk) + warning('Curve %d has no significant peaks.', k); + fitResults(k) = fillZeroStruct(theta, opts.MaxTheta); + continue; + end + + thetaPeaks = thetaPeaks(:); + pk = pk(:); + + % --- Find secondary peak --- + secondaryIdx = find((thetaPeaks > opts.PositionThreshold) & ... + (pk >= opts.AmplitudeThreshold * curveAmp), 1, 'first'); + + if isempty(secondaryIdx) + % No secondary peak → bypass fit, fill with zeros + fitResults(k) = fillZeroStruct(theta, opts.MaxTheta); + continue; + else + fitMaxTheta = thetaPeaks(secondaryIdx); + secondaryAmp = pk(secondaryIdx); + end + + % --- Use the full curve for fitting --- + thetaFit = theta(:); + xFit = x(:); + + if numel(thetaFit) < 8 + warning('Curve %d has too few points for fitting.', k); + fitResults(k) = fillNaNStruct(); + continue; + end + + % --- Two-Gaussian model --- + twoGauss = @(p,theta) ... + 1.0 * exp(-0.5*((theta - p(1))/max(p(2),1e-6)).^2) + ... % main Gaussian + p(3) * exp(-0.5*((theta - p(4))/max(p(5),1e-6)).^2); % secondary Gaussian + + % --- Parameter guesses based on secondary peak --- + mu1_guess = 0; + sigma1_guess = 0.1 * fitMaxTheta; + A2_guess = secondaryAmp; % use detected secondary amplitude + mu2_guess = fitMaxTheta; % center at detected secondary peak + sigma2_guess = fitMaxTheta / 3; % broad guess + + p0 = [mu1_guess, sigma1_guess, A2_guess, mu2_guess, sigma2_guess]; + + % --- Bounds that enforce secondary peak constraints --- + lb = [0, 1e-6, 0, max(fitMaxTheta - 0.1, 0), 1e-6]; + ub = [fitMaxTheta/4, fitMaxTheta, 2, min(fitMaxTheta + 0.1, opts.MaxTheta), fitMaxTheta]; + + optsLSQ = optimoptions('lsqcurvefit','Display','off', ... + 'MaxFunctionEvaluations',1e4,'MaxIterations',1e4); + + % --- Fit --- + try + pFit = lsqcurvefit(twoGauss, p0, thetaFit, xFit, lb, ub, optsLSQ); + catch + warning('Curve %d fit failed.', k); + fitResults(k) = fillNaNStruct(); + continue; + end + + % --- Evaluate fitted curve --- + thetaFine = linspace(0, opts.MaxTheta, 500); + yFit = twoGauss(pFit, thetaFine); + + % --- Assess deviation --- + relDevAmplitude = abs(pFit(3) - secondaryAmp) / max(secondaryAmp,1e-6); + relDevPosition = abs(pFit(4) - fitMaxTheta) / max(fitMaxTheta, 1e-6); + combinedDeviation = sqrt(relDevAmplitude^2 + relDevPosition^2); + isValid = combinedDeviation <= opts.DeviationLimit; + + % --- Store results --- + fitResults(k).pFit = pFit; + fitResults(k).thetaFit = thetaFit; + fitResults(k).xFit = xFit; + fitResults(k).thetaFine = thetaFine; + fitResults(k).yFit = yFit; + fitResults(k).isValid = isValid; + fitResults(k).fitMaxTheta = fitMaxTheta; + catch - warning('Curve %d fit failed, using initial guess.', k); - pFit = p0; + warning('Curve %d fit failed completely.', k); + fitResults(k) = fillNaNStruct(); end - - % --- Evaluate fitted curve --- - thetaFine = linspace(0, opts.MaxTheta, 500); - yFit = twoGauss(pFit, thetaFine); - - % --- Compute relative deviation --- - yFitInterp = twoGauss(pFit, thetaFit); - relDeviation = abs(yFitInterp - xFit) ./ max(xFit, 1e-6); - maxRelDev = max(relDeviation); - - % --- Goodness-of-fit --- - isValid = maxRelDev <= opts.DeviationLimit; - - % --- Store results --- - fitResults(k).pFit = pFit; - fitResults(k).thetaFit = thetaFit; - fitResults(k).xFit = xFit; - fitResults(k).thetaFine = thetaFine; - fitResults(k).yFit = yFit; - fitResults(k).isValid = isValid; - fitResults(k).fitMaxTheta = fitMaxTheta; end end + +% --- Helper: return struct with NaNs (fit failed) --- +function s = fillNaNStruct() + s = struct('pFit',nan(1,5),'thetaFit',nan,'xFit',nan,'yFit',nan,'thetaFine',nan,... + 'isValid',false,'fitMaxTheta',nan); +end + +% --- Helper: return struct with zeros (no secondary peak found) --- +function s = fillZeroStruct(theta, fitMaxTheta) + s = struct('pFit',zeros(1,5),'thetaFit',theta,'xFit',zeros(size(theta)),... + 'yFit',zeros(size(theta)),'thetaFine',zeros(1,500),... + 'isValid',false,'fitMaxTheta',fitMaxTheta); +end diff --git a/Data-Analyzer/+Analyzer/recenterSpectralCurves.m b/Data-Analyzer/+Analyzer/recenterSpectralCurves.m new file mode 100644 index 0000000..d4cf65d --- /dev/null +++ b/Data-Analyzer/+Analyzer/recenterSpectralCurves.m @@ -0,0 +1,88 @@ +function results = recenterSpectralCurves(S_theta_all, theta_vals, scan_reference_values, varargin) +%% recenterSpectralCurves +% Author: Karthik +% Date: 2025-10-09 +% Version: 1.0 +% +% Description: +% Recenters each curve so its first peak is at zero (no wrap-around), +% pads shorter curves with NaN, and computes mean and SEM for each +% scan parameter. Returns a struct ready for plotting. +% +% Inputs: +% S_theta_all - 1x(N_reps*N_params) cell array of curves +% theta_vals - vector of theta values (radians) +% scan_reference_values - vector of unique scan parameters +% +% Optional name-value pairs: +% 'SearchRange' - [deg_min deg_max] range for first-peak search (default: [0 90]) +% +% Output struct fields: +% results.curves +% results.x_values +% results.curves_mean +% results.curves_error +% results.scan_parameter_values + + % --- Parse optional inputs --- + p = inputParser; + addParameter(p, 'SearchRange', [0 pi/2], @(x) isnumeric(x) && numel(x)==2); + parse(p, varargin{:}); + opts = p.Results; + + % --- Setup --- + N_params = numel(scan_reference_values); + Ntotal = numel(S_theta_all); + Nreps = Ntotal / N_params; + theta_min = opts.SearchRange(1); + theta_max = opts.SearchRange(2); + + % --- Shift each curve so its first peak is at zero --- + S_theta_all_shifted = cell(size(S_theta_all)); + for k = 1:Ntotal + curve = S_theta_all{k}; + + idx_range = find(theta_vals >= theta_min & theta_vals <= theta_max); + if isempty(idx_range) + [~, peak_idx] = max(curve); + else + [~, local_max_idx] = max(curve(idx_range)); + peak_idx = idx_range(local_max_idx); + end + + % From first peak onward + S_theta_all_shifted{k} = curve(peak_idx:end); + end + + % --- Pad shorter curves with NaN --- + Npoints_shifted = max(cellfun(@numel, S_theta_all_shifted)); + for k = 1:Ntotal + len = numel(S_theta_all_shifted{k}); + if len < Npoints_shifted + S_theta_all_shifted{k} = [S_theta_all_shifted{k}, nan(1, Npoints_shifted - len)]; + end + end + + % --- Compute mean and SEM --- + curves_all = cell(1, N_params); + curves_mean = cell(1, N_params); + curves_err = cell(1, N_params); + + for i = 1:N_params + curves = zeros(Nreps, Npoints_shifted); + for r = 1:Nreps + idx = (r-1)*N_params + i; % interleaved indexing + curves(r,:) = S_theta_all_shifted{idx}; + end + curves_all{i} = curves; + curves_mean{i} = nanmean(curves,1); + curves_err{i} = nanstd(curves,0,1)/sqrt(Nreps); + end + + % --- Construct results struct --- + results.curves = curves_all; + results.x_values = theta_vals(1:Npoints_shifted) - theta_vals(1); % shift start to zero + results.curves_mean = curves_mean; + results.curves_error = curves_err; + results.scan_reference_values = scan_reference_values; +end diff --git a/Data-Analyzer/+Calculator/computeAngularSpectralDistribution.m b/Data-Analyzer/+Calculator/computeAngularSpectralDistribution.m index 8e203f4..10feea7 100644 --- a/Data-Analyzer/+Calculator/computeAngularSpectralDistribution.m +++ b/Data-Analyzer/+Calculator/computeAngularSpectralDistribution.m @@ -1,14 +1,16 @@ -function [theta_vals, S_theta] = computeAngularSpectralDistribution(IMGFFT, kx, ky, k_min, k_max, num_bins, threshold, sigma, windowSize) +function [theta_vals, S_theta] = computeAngularSpectralDistribution(IMGFFT, kx, ky, k_min, k_max, num_bins, threshold, sigma, windowSize, theta_max) %% computeAngularSpectralDistribution % Author: Karthik % Date: 2025-09-12 -% Version: 1.0 +% Version: 1.1 % % Description: -% Brief description of the script functionality. -% -% Notes: -% Optional notes, references. +% Computes the angular spectral distribution within a specified radial band. +% Supports custom angular range up to 2π. + + if nargin < 10 || isempty(theta_max) + theta_max = pi; % default preserves current behavior + end % Apply threshold to isolate strong peaks IMGFFT(IMGFFT < threshold) = 0; @@ -18,17 +20,20 @@ function [theta_vals, S_theta] = computeAngularSpectralDistribution(IMGFFT, kx, Kmag = sqrt(KX.^2 + KY.^2); % radial wavenumber magnitude Theta = atan2(KY, KX); % range [-pi, pi] - % Restrict to radial band in wavenumber space + % 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); % Initialize angular structure factor S_theta = zeros(1, num_bins); - theta_vals = linspace(0, pi, num_bins); % only 0 to pi due to symmetry + theta_vals = linspace(0, theta_max, num_bins); % Loop over angular bins for i = 1:num_bins - angle_start = (i - 1) * pi / num_bins; - angle_end = i * pi / 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; diff --git a/Data-Analyzer/+Plotter/plotFitParameterPDF.m b/Data-Analyzer/+Plotter/plotFitParameterPDF.m index f956013..2afd49d 100644 --- a/Data-Analyzer/+Plotter/plotFitParameterPDF.m +++ b/Data-Analyzer/+Plotter/plotFitParameterPDF.m @@ -2,7 +2,7 @@ function plotFitParameterPDF(fitResults, scanValues, paramName, varargin) %% plotFitParameterPDF % Author: Karthik % Date: 2025-10-06 -% Version: 1.0 +% Version: 1.1 % % Description: % Plots 2D PDF (heatmap) of any parameter from two-Gaussian fit results @@ -56,8 +56,16 @@ paramValues = nan(N_reps, N_params); for k = 1:N_total paramIdxScan = mod(k-1, N_params) + 1; repIdx = floor((k-1)/N_params) + 1; - if fitResults(k).isValid - paramValues(repIdx, paramIdxScan) = fitResults(k).pFit(paramIdx); + paramValues(repIdx, paramIdxScan) = fitResults(k).pFit(paramIdx); +end + +% --- Create mask of true zeros (from fillZeroStruct) --- +trueZeroMask = false(size(paramValues)); +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) + trueZeroMask(repIdx, paramIdxScan) = true; end end @@ -105,13 +113,21 @@ for i = 1:N_params end end +% --- Mask out non-data regions (make them white) --- +dataMask = ~isnan(paramValues); % where valid fits exist +maskPerScan = any(dataMask,1); % scans that have any valid data +pdf_matrix(:, ~maskPerScan) = NaN; % make empty scans white + % --- Plot heatmap --- fig = figure(opts.FigNum); clf(fig); set(fig, 'Color', 'w', 'Position', [100 100 950 750]); + if strcmpi(opts.PlotType,'kde') - imagesc(scanValues, y_grid, pdf_matrix); + h = imagesc(scanValues, y_grid, pdf_matrix); + set(h, 'AlphaData', ~isnan(pdf_matrix)); % make NaNs transparent else - imagesc(scanValues, binCenters, pdf_matrix); + h = imagesc(scanValues, binCenters, pdf_matrix); + set(h, 'AlphaData', ~isnan(pdf_matrix)); % make NaNs transparent end set(gca, 'YDir', 'normal', 'FontName', opts.FontName, 'FontSize', opts.FontSize); @@ -119,7 +135,9 @@ xlabel(opts.XLabel, 'FontSize', opts.FontSize, 'FontName', opts.FontName); ylabel(opts.YLabel, 'FontSize', opts.FontSize, 'FontName', opts.FontName); title(opts.Title, 'FontSize', opts.FontSize+2, 'FontWeight', 'bold'); -colormap(feval(opts.Colormap)); +% --- Colormap and colorbar --- +cmap = feval(opts.Colormap); +colormap([1 1 1; cmap]); % prepend white for NaN regions c = colorbar; if strcmpi(opts.PlotType,'kde') || opts.NormalizeHistogram ylabel(c, 'Probability Density', 'Rotation', -90, 'FontName', opts.FontName, 'FontSize', opts.FontSize); @@ -142,9 +160,4 @@ if opts.OverlayMeanSEM plot(scanValues, meanParam, 'k-', 'LineWidth', 2); end -% --- Save figure --- -if ~opts.SkipSaveFigures - if ~exist(opts.SaveDirectory,'dir'), mkdir(opts.SaveDirectory); end - savefig(fig, fullfile(opts.SaveDirectory, opts.SaveFileName)); -end end diff --git a/Data-Analyzer/+Plotter/plotG2Curves.m b/Data-Analyzer/+Plotter/plotG2Curves.m index 5fd7678..b260d36 100644 --- a/Data-Analyzer/+Plotter/plotG2Curves.m +++ b/Data-Analyzer/+Plotter/plotG2Curves.m @@ -30,7 +30,7 @@ function plotG2Curves(results, varargin) addParameter(p, 'SaveFileName', 'G2Curves.fig', @ischar); addParameter(p, 'SaveDirectory', pwd, @ischar); addParameter(p, 'TileTitlePrefix', 'Control Parameter', @(x) ischar(x) || isstring(x)); - addParameter(p, 'TileTitleSuffix', '', @(x) ischar(x) || isstring(x)); % NEW + addParameter(p, 'TileTitleSuffix', '', @(x) ischar(x) || isstring(x)); addParameter(p, 'HighlightEvery', 10, @(x) isnumeric(x) && isscalar(x) && x>=1); parse(p, varargin{:}); opts = p.Results; diff --git a/Data-Analyzer/+Plotter/plotSpectralCurvesRecentered.m b/Data-Analyzer/+Plotter/plotSpectralCurvesRecentered.m index d9a1269..b3209ce 100644 --- a/Data-Analyzer/+Plotter/plotSpectralCurvesRecentered.m +++ b/Data-Analyzer/+Plotter/plotSpectralCurvesRecentered.m @@ -1,4 +1,4 @@ -function results = plotSpectralCurvesRecentered(S_theta_all, theta_vals, scan_reference_values, varargin) +function plotSpectralCurvesRecentered(results, scan_reference_values, varargin) %% plotSpectralCurves % Author: Karthik % Date: 2025-10-01 @@ -49,65 +49,7 @@ function results = plotSpectralCurvesRecentered(S_theta_all, theta_vals, scan_re % --- Determine number of scan parameters and repetitions --- N_params = numel(scan_reference_values); - Ntotal = numel(S_theta_all); - Nreps = Ntotal / N_params; - theta_min = deg2rad(0); - theta_max = deg2rad(90); - - % --- Preprocess curves: shift first peak to zero without wrapping --- - S_theta_all_shifted = cell(size(S_theta_all)); - for k = 1:Ntotal - curve = S_theta_all{k}; - - % --- Find first peak only in specified theta range --- - - idx_range = find(theta_vals >= theta_min & theta_vals <= theta_max); - if isempty(idx_range) - % fallback: search entire curve if range is empty - [~, peak_idx] = max(curve); - else - [~, local_max_idx] = max(curve(idx_range)); - peak_idx = idx_range(local_max_idx); - end - - % Take only the part from peak to end - S_theta_all_shifted{k} = curve(peak_idx:end); - end - - % --- Adjust Npoints for plotting --- - Npoints_shifted = max(cellfun(@numel, S_theta_all_shifted)); - for k = 1:Ntotal - % Pad shorter curves with NaN to keep sizes consistent - len = numel(S_theta_all_shifted{k}); - if len < Npoints_shifted - S_theta_all_shifted{k} = [S_theta_all_shifted{k}, nan(1, Npoints_shifted-len)]; - end - end - - % --- Prepare curves, mean, SEM --- - curves_all = cell(1, N_params); - curves_mean = cell(1, N_params); - curves_err = cell(1, N_params); - - for i = 1:N_params - curves = zeros(Nreps, Npoints_shifted); - for r = 1:Nreps - idx = (r-1)*N_params + i; % correct interleaved indexing - curves(r,:) = S_theta_all_shifted{idx}; - end - curves_all{i} = curves; - curves_mean{i} = nanmean(curves,1); - curves_err{i} = nanstd(curves,0,1)/sqrt(Nreps); - end - - % --- Create results struct compatible with plotting --- - results.curves = curves_all; - results.x_values = theta_vals(1:Npoints_shifted) - theta_vals(1); % center first peak at zero - results.curves_mean = curves_mean; - results.curves_error = curves_err; - results.scan_parameter_values = scan_reference_values; - % --- Create figure --- if isempty(opts.FigNum), fig = figure; else, fig = figure(opts.FigNum); end clf(fig); @@ -149,7 +91,7 @@ function results = plotSpectralCurvesRecentered(S_theta_all, theta_vals, scan_re grid(ax,'on'); xlabel(ax, opts.XLabel, 'FontName', opts.FontName, 'FontSize', opts.FontSize); 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), ... + title(ax, sprintf('%s=%.3g%s', opts.TileTitlePrefix, results.scan_reference_values(i), opts.TileTitleSuffix), ... 'FontName', opts.FontName, 'FontSize', opts.FontSize); set(ax, 'FontName', opts.FontName, 'FontSize', opts.FontSize); end diff --git a/Data-Analyzer/+Plotter/plotSpectralDistributionCumulants.m b/Data-Analyzer/+Plotter/plotSpectralDistributionCumulants.m index d7ba5c7..897dd3e 100644 --- a/Data-Analyzer/+Plotter/plotSpectralDistributionCumulants.m +++ b/Data-Analyzer/+Plotter/plotSpectralDistributionCumulants.m @@ -31,7 +31,7 @@ function plotSpectralDistributionCumulants(results, varargin) % --- Setup --- N_params = numel(results.curves); - xvals = results.scan_parameter_values; + xvals = results.scan_reference_values; thetaVals = results.x_values * pi; % --- Select specific theta indices (like in g2 case) --- diff --git a/Data-Analyzer/+Scripts/BECToDropletsToStripes/plotAnalysisResults.m b/Data-Analyzer/+Scripts/BECToDropletsToStripes/plotAnalysisResults.m index 3760845..967114f 100644 --- a/Data-Analyzer/+Scripts/BECToDropletsToStripes/plotAnalysisResults.m +++ b/Data-Analyzer/+Scripts/BECToDropletsToStripes/plotAnalysisResults.m @@ -187,251 +187,13 @@ Plotter.plotMultiplePCAResults(compiled_results.pca_results, scan_parameter_valu 'SkipSaveFigures', options.skipSaveFigures, ... 'SaveDirectory', figSaveDir); -%% ------------------ 7. Plot of all Angular Spectral Distribution Curves ------------------ -Plotter.plotSpectralCurves( ... - results_all{1}.results.spectral_analysis_results.S_theta_norm_all, ... - results_all{1}.results.spectral_analysis_results.theta_vals/pi, ... % correct θ values - results_all{1}.scan_reference_values, ... % correct scan params - 'Title', options.titleString, ... - 'XLabel', '\theta / \pi', ... - 'YLabel', 'S(\theta)', ... - 'HighlightEvery', 10, ... % highlight every 10th repetition - 'FontName', options.font, ... - 'FigNum', 20, ... - 'TileTitlePrefix', '\alpha', ... % user-defined tile prefix - 'TileTitleSuffix', '^\circ', ... % user-defined suffix (e.g., degrees symbol) - 'SkipSaveFigures', options.skipSaveFigures, ... - 'SaveFileName', 'SpectralCurves.fig', ... - 'SaveDirectory', figSaveDir); - -%% ------------------ 8. Plot of all Angular Spectral Distribution Curves shifted ------------------ -results = Plotter.plotSpectralCurvesRecentered( ... - results_all{1}.results.spectral_analysis_results.S_theta_norm_all, ... - results_all{1}.results.spectral_analysis_results.theta_vals/pi, ... % correct θ values - results_all{1}.scan_reference_values, ... % correct scan params - 'Title', options.titleString, ... - 'XLabel', '\theta / \pi', ... - 'YLabel', 'S(\theta)', ... - 'HighlightEvery', 10, ... % highlight every 10th repetition - 'FontName', options.font, ... - 'FigNum', 21, ... - 'TileTitlePrefix', '\alpha', ... % user-defined tile prefix - 'TileTitleSuffix', '^\circ', ... % user-defined suffix (e.g., degrees symbol) - 'SkipSaveFigures', options.skipSaveFigures, ... - 'SaveFileName', 'SpectralCurves.fig', ... - 'SaveDirectory', figSaveDir); -%% ------------------ 9. Plot cumulants from shifted Angular Spectral Distribution Curves ------------------ -Plotter.plotSpectralDistributionCumulants(results, ... - 'Title', options.titleString, ... - 'XLabel', '\alpha (degrees)', ... - 'FontName', options.font, ... - 'FontSize', 14, ... - 'FigNum', 22, ... - 'SkipSaveFigures', false, ... - 'SaveFileName', 'SpectralCumulants.fig'); - -%% ------------------ 10. Fit shifted Angular Spectral Distribution Curves ------------------ -fitResults = Analyzer.fitTwoGaussianCurves(... - results_all{1}.results.spectral_analysis_results.S_theta_norm_all, ... - results_all{1}.results.spectral_analysis_results.theta_vals, ... - 'MaxTheta', pi/2, ... - 'DeviationLimit', 1.00); - -%% ------------------ 11. Plot fit parameters - position ------------------ -Plotter.plotFitParameterPDF(fitResults, results_all{1}.scan_reference_values, 'mu2', ... - 'Title', options.titleString, ... - 'XLabel', '\alpha (degrees)', ... - 'YLabel', 'Secondary peak position (\theta, rad)', ... - 'FontName', options.font, ... - 'FontSize', 16, ... - 'FigNum', 23, ... - 'SkipSaveFigures', options.skipSaveFigures, ... - 'SaveFileName', 'SecondaryPeakPositionPDF.fig', ... - 'PlotType', 'histogram', ... - 'NumberOfBins', 20, ... - 'NormalizeHistogram', true, ... - 'Colormap', @Colormaps.coolwarm, ... - 'XLim', [min(results_all{1}.scan_reference_values), max(results_all{1}.scan_reference_values)]); - -%% ------------------ 12. Plot fit parameters - width ------------------ -Plotter.plotFitParameterPDF(fitResults, results_all{1}.scan_reference_values, 'sigma2', ... - 'Title', options.titleString, ... - 'XLabel', '\alpha (degrees)', ... - 'YLabel', 'Secondary peak width (\sigma, rad)', ... - 'FontName', options.font, ... - 'FontSize', 16, ... - 'FigNum', 24, ... - 'SkipSaveFigures', options.skipSaveFigures, ... - 'SaveFileName', 'SecondaryPeakWidthPDF.fig', ... - 'PlotType', 'histogram', ... - 'NumberOfBins', 20, ... - 'NormalizeHistogram', true, ... - 'Colormap', @Colormaps.coolwarm, ... - 'XLim', [min(results_all{1}.scan_reference_values), max(results_all{1}.scan_reference_values)]); - -%% ------------------ 13. Plot fit parameters of mean shifted Angular Spectral Distribution Curves ------------------ -S_theta_all = results_all{1}.results.spectral_analysis_results.S_theta_norm_all; -theta_vals = results_all{1}.results.spectral_analysis_results.theta_vals; -scanValues = results_all{1}.scan_reference_values; - -N_params = numel(scanValues); -N_total = numel(S_theta_all); -N_reps = N_total / N_params; - -theta_min = deg2rad(0); -theta_max = deg2rad(90); - -% --- Shift curves so first peak is at start --- -S_theta_all_shifted = cell(size(S_theta_all)); -for k = 1:N_total - curve = S_theta_all{k}; - idx_range = find(theta_vals >= theta_min & theta_vals <= theta_max); - if isempty(idx_range) - [~, peak_idx] = max(curve); - else - [~, local_max_idx] = max(curve(idx_range)); - peak_idx = idx_range(local_max_idx); - end - S_theta_all_shifted{k} = curve(peak_idx:end); -end - -% --- Pad shorter curves with NaN to match lengths --- -Npoints_shifted = max(cellfun(@numel, S_theta_all_shifted)); -for k = 1:N_total - len = numel(S_theta_all_shifted{k}); - if len < Npoints_shifted - S_theta_all_shifted{k} = [S_theta_all_shifted{k}, nan(1, Npoints_shifted-len)]; - end -end - -% --- Compute mean curves per scan parameter --- -meanCurves = cell(1, N_params); -for i = 1:N_params - curves = zeros(N_reps, Npoints_shifted); - for r = 1:N_reps - idx = (r-1)*N_params + i; % interleaved indexing - curves(r,:) = S_theta_all_shifted{idx}; - end - meanCurves{i} = nanmean(curves,1); % mean over repetitions -end - -% --- Fit two-Gaussian model to mean curves --- -fitResultsMean = Analyzer.fitTwoGaussianCurves(meanCurves, theta_vals(1:Npoints_shifted)-theta_vals(1), ... - 'MaxTheta', pi/2, ... - 'DeviationLimit', 1.0); - -% --- Scatter plot of secondary peak position (mu2) vs scan parameter --- -mu2_vals = nan(1, N_params); -sigma2_vals = nan(1, N_params); - -for i = 1:N_params - if fitResultsMean(i).isValid - mu2_vals(i) = fitResultsMean(i).pFit(4); % secondary peak position - sigma2_vals(i) = fitResultsMean(i).pFit(5); % secondary peak width - end -end - -% Secondary peak position -plotSecondaryPeakScatter(fitResultsMean, results_all{1}.scan_reference_values, 'mu2', ... - 'Title', options.titleString, ... - 'XLabel', '\alpha (degrees)', ... - 'YLabel', 'Secondary peak position (\theta, rad)', ... - 'FigNum', 23, ... - 'FontName', options.font, ... - 'FontSize', 16, ... - 'SkipSaveFigures', options.skipSaveFigures, ... - 'SaveFileName', 'SecondaryPeakPositionScatter.fig'); - -% Secondary peak width -plotSecondaryPeakScatter(fitResultsMean, results_all{1}.scan_reference_values, 'sigma2', ... - 'Title', options.titleString, ... - 'XLabel', '\alpha (degrees)', ... - 'YLabel', 'Secondary peak width (\sigma, rad)', ... - 'FigNum', 24, ... - 'FontName', options.font, ... - 'FontSize', 16, ... - 'SkipSaveFigures', options.skipSaveFigures, ... - 'SaveFileName', 'SecondaryPeakWidthScatter.fig'); - - -function plotSecondaryPeakScatter(fitResultsMean, scanValues, parameterName, varargin) -%% plotSecondaryPeakScatter -% Author: Karthik -% Date: 2025-10-06 -% Version: 1.0 -% -% Description: -% Scatter plot of secondary peak fit parameters (mu2 or sigma2) vs scan parameter -% in the same style as plotMeanWithSE. - - % --- 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', [], @(x) isempty(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', 'secondary_peak_scatter.fig', @ischar); - addParameter(p, 'SaveDirectory', pwd, @ischar); - parse(p, varargin{:}); - opts = p.Results; - - % --- Extract parameter values --- - N_params = numel(fitResultsMean); - paramVals = nan(1, N_params); - for i = 1:N_params - if fitResultsMean(i).isValid - switch parameterName - case 'mu2' - paramVals(i) = fitResultsMean(i).pFit(4); - case 'sigma2' - paramVals(i) = fitResultsMean(i).pFit(5); - otherwise - error('Unknown parameter name: %s', parameterName); - end - end - end - - % --- Prepare figure --- - if isempty(opts.FigNum) - fig = figure; - else - fig = figure(opts.FigNum); - clf(fig); - end - set(fig, 'Color', 'w', 'Position', [100 100 950 750]); - - % --- Plot as mean ± SE style (SE=0 here) --- - errorbar(scanValues, paramVals, zeros(size(paramVals)), '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', 'tex', 'FontName', opts.FontName, 'FontSize', opts.FontSize); - ylabel(opts.YLabel, 'Interpreter', 'tex', 'FontName', opts.FontName, 'FontSize', opts.FontSize); - title(opts.Title, 'FontName', opts.FontName, 'FontSize', opts.FontSize + 2, 'FontWeight', 'bold'); - grid on; - - % --- Save figure --- - if ~opts.SkipSaveFigures - if ~exist(opts.SaveDirectory, 'dir'), mkdir(opts.SaveDirectory); end - savefig(fig, fullfile(opts.SaveDirectory, opts.SaveFileName)); - end -end - %{ %% ------------------ 14. Average of Spectra Plots ------------------ Plotter.plotAverageSpectra(scan_parameter_values, ... spectral_analysis_results, ... 'ScanParameterName', scan_parameter, ... - 'FigNum', 25, ... + 'FigNum', 20, ... 'ColormapPS', Colormaps.coolwarm(), ... 'Font', 'Bahnschrift', ... 'SaveFileName', 'avgSpectra.fig', ... diff --git a/Data-Analyzer/+Scripts/BECToDropletsToStripes/runSpectralDistributionAnalysis.m b/Data-Analyzer/+Scripts/BECToDropletsToStripes/runSpectralDistributionAnalysis.m new file mode 100644 index 0000000..0de0e5a --- /dev/null +++ b/Data-Analyzer/+Scripts/BECToDropletsToStripes/runSpectralDistributionAnalysis.m @@ -0,0 +1,342 @@ +%% ===== 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; + +% +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; + +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.extractFullAngularSpectralDistribution(od_imgs, options); + +%% ------------------ 1. Plot of all Angular Spectral Distribution Curves ------------------ +Plotter.plotSpectralCurves( ... + spectral_analysis_results.S_theta_norm_all, ... + spectral_analysis_results.theta_vals/pi, ... % correct θ values + scan_reference_values, ... % correct scan params + 'Title', options.titleString, ... + 'XLabel', '\theta / \pi', ... + 'YLabel', 'S(\theta)', ... + '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. Plot of all Angular Spectral Distribution Curves shifted ------------------ + +% --- Recenter curves first --- +results = Analyzer.recenterSpectralCurves(spectral_analysis_results.S_theta_norm_all, ... + spectral_analysis_results.theta_vals/pi, ... + scan_reference_values, ... + 'SearchRange', [0 90]); % degrees + +% --- Restrict to desired theta range (e.g., 0 to 0.5*pi) --- +thetaMin = 0; % in units of pi (since you divided by pi) +thetaMax = 1; % corresponds to pi/2 + +mask = results.x_values >= thetaMin & results.x_values <= thetaMax; +results.x_values = results.x_values(mask); + +% Apply the same mask to each curve set +for i = 1:numel(results.curves) + results.curves{i} = results.curves{i}(:, mask); + results.curves_mean{i} = results.curves_mean{i}(mask); + results.curves_error{i}= results.curves_error{i}(mask); +end + +Plotter.plotSpectralCurvesRecentered( ... + results, ... + scan_reference_values, ... % correct scan params + 'Title', options.titleString, ... + 'XLabel', '\theta / \pi', ... + 'YLabel', 'S(\theta)', ... + 'HighlightEvery', 10, ... % highlight every 10th repetition + 'FontName', options.font, ... + 'FigNum', 2, ... + 'TileTitlePrefix', '\alpha', ... % user-defined tile prefix + 'TileTitleSuffix', '^\circ', ... % user-defined suffix (e.g., degrees symbol) + 'SkipSaveFigures', options.skipSaveFigures, ... + 'SaveFileName', 'SpectralCurves.fig', ... + 'SaveDirectory', options.saveDirectory); + +%% ------------------ 3. Plot cumulants from shifted Angular Spectral Distribution Curves ------------------ +Plotter.plotSpectralDistributionCumulants(results, ... + 'Title', options.titleString, ... + 'XLabel', '\alpha (degrees)', ... + 'FontName', options.font, ... + 'FontSize', 14, ... + 'FigNum', 3, ... + 'SkipSaveFigures', false, ... + 'SaveFileName', 'SpectralCumulants.fig'); + +%% ------------------ 4. Fit shifted Angular Spectral Distribution Curves ------------------ +fitResults = Analyzer.fitTwoGaussianCurves(... + spectral_analysis_results.S_theta_norm_all, ... + spectral_analysis_results.theta_vals, ... + 'MaxTheta', pi/2, ... + 'DeviationLimit', 0.30, ... + 'PositionThreshold', pi/15, ... + 'AmplitudeThreshold', 0.02); + +%% ------------------ 5. 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)]); + +%% ------------------ 6. 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)]); + +%% ------------------ 7. 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)]); + +%% ------------------ 8. Plot fit parameters of mean shifted Angular Spectral Distribution Curves ------------------ + +% --- Recenter curves first --- +results = Analyzer.recenterSpectralCurves(spectral_analysis_results.S_theta_norm_all, ... + spectral_analysis_results.theta_vals/pi, ... + scan_reference_values, ... + 'SearchRange', [0 90]); % degrees + +% --- Restrict to desired theta range (e.g., 0 to 0.5*pi) --- +thetaMin = 0; % in units of pi (since you divided by pi) +thetaMax = 1; % corresponds to pi/2 + +mask = results.x_values >= thetaMin & results.x_values <= thetaMax; +results.x_values = results.x_values(mask); + +% Apply the same mask to each curve set +for i = 1:numel(results.curves) + results.curves_mean{i} = results.curves_mean{i}(mask); +end + +% --- Fit two-Gaussian model to mean curves --- +fitResultsMean = Analyzer.fitTwoGaussianCurves(... + results.curves_mean, ... + results.x_values*pi, ... + 'MaxTheta', pi/2, ... + 'DeviationLimit', 0.30, ... + 'PositionThreshold', pi/15, ... + 'AmplitudeThreshold', 0.02, ... + 'RecenterCurves', false); + +% --- Prepare parameter values --- +N_params = numel(fitResultsMean); +amp2_vals = nan(1, N_params); +mu2_vals = nan(1, N_params); +sigma2_vals = nan(1, N_params); + +for i = 1:N_params + pFit = fitResultsMean(i).pFit; + + if all(pFit == 0) + % No secondary peak → plot as zero + amp2_vals(i) = 0; + mu2_vals(i) = 0; + sigma2_vals(i) = 0; + + elseif all(~isnan(pFit)) + % Successful fit → use fitted values + amp2_vals(i) = pFit(3); + mu2_vals(i) = pFit(4); + sigma2_vals(i) = pFit(5); + + else + % Fit failed → leave as NaN (skipped automatically) + continue; + end +end + +% --- Plot using plotMeanWithSE --- +Plotter.plotMeanWithSE(scan_reference_values, amp2_vals, ... + 'Title', options.titleString, ... + 'XLabel', '\alpha (degrees)', ... + 'YLabel', 'Secondary peak amplitude', ... + 'FigNum', 7, ... + 'FontName', options.font, ... + 'FontSize', 16); + +% --- Plot using plotMeanWithSE --- +Plotter.plotMeanWithSE(scan_reference_values, mu2_vals, ... + 'Title', options.titleString, ... + 'XLabel', '\alpha (degrees)', ... + 'YLabel', 'Secondary peak position (\theta, rad)', ... + 'FigNum', 8, ... + 'FontName', options.font, ... + 'FontSize', 16); + +Plotter.plotMeanWithSE(scan_reference_values, sigma2_vals, ... + 'Title', options.titleString, ... + 'XLabel', '\alpha (degrees)', ... + 'YLabel', 'Secondary peak width (\sigma, rad)', ... + 'FigNum', 9, ... + 'FontName', options.font, ... + 'FontSize', 16); + +%% Inspect individual realizations of a single parameter + +% --- Recenter curves first --- +results = Analyzer.recenterSpectralCurves( ... + spectral_analysis_results.S_theta_norm_all, ... + spectral_analysis_results.theta_vals/pi, ... + scan_reference_values, ... + 'SearchRange', [0 90]); % degrees + +% --- Restrict to desired theta range (e.g., 0 to 0.5*pi) --- +thetaMin = 0; % in units of pi (since you divided by pi) +thetaMax = 1; % corresponds to pi/2 + +mask = results.x_values >= thetaMin & results.x_values <= thetaMax; +results.x_values = results.x_values(mask); + +% --- Apply the same mask to each curve set (1x10 cell, each 60x180) --- +for i = 1:numel(results.curves) + results.curves{i} = results.curves{i}(:, mask); +end + +% --- Convert selected curve set (e.g., 5th) into 1x60 cell array of 1xN row vectors --- +paramIdx = 10; % <-- choose which scan point or curve set to analyze +curves_matrix = results.curves{paramIdx}; % 60xN numeric +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 --- +fitResultsMean = Analyzer.fitTwoGaussianCurves( ... + curves_cell, ... + results.x_values*pi, ... + 'MaxTheta', pi/2, ... + 'DeviationLimit', 0.30, ... + 'PositionThreshold', pi/15, ... + 'AmplitudeThreshold', 0.02, ... + 'RecenterCurves', false); diff --git a/Data-Analyzer/+Scripts/BECToStripesToDroplets/plotAnalysisResults.m b/Data-Analyzer/+Scripts/BECToStripesToDroplets/plotAnalysisResults.m index 74e8d41..bda2f3a 100644 --- a/Data-Analyzer/+Scripts/BECToStripesToDroplets/plotAnalysisResults.m +++ b/Data-Analyzer/+Scripts/BECToStripesToDroplets/plotAnalysisResults.m @@ -106,8 +106,7 @@ Plotter.plotG2Curves(compiled_results.full_g2_results, ... 'SaveFileName', 'G2Curves.fig', ... 'SaveDirectory', figSaveDir); - -Plotter.plotG2PDF(compiled_results.full_g2_results, pi/6, ... +Plotter.plotG2PDF(compiled_results.full_g2_results, pi/2, ... 'Title', options.titleString, ... 'XLabel', '\alpha (degrees)', ... 'YLabel', 'g^{(2)}(\delta\theta)', ... @@ -188,251 +187,13 @@ Plotter.plotMultiplePCAResults(compiled_results.pca_results, scan_parameter_valu 'SkipSaveFigures', options.skipSaveFigures, ... 'SaveDirectory', figSaveDir); -%% ------------------ 7. Plot of all Angular Spectral Distribution Curves ------------------ -Plotter.plotSpectralCurves( ... - results_all{1}.results.spectral_analysis_results.S_theta_norm_all, ... - results_all{1}.results.spectral_analysis_results.theta_vals/pi, ... % correct θ values - results_all{1}.scan_reference_values, ... % correct scan params - 'Title', options.titleString, ... - 'XLabel', '\theta / \pi', ... - 'YLabel', 'S(\theta)', ... - 'HighlightEvery', 10, ... % highlight every 10th repetition - 'FontName', options.font, ... - 'FigNum', 20, ... - 'TileTitlePrefix', '\alpha', ... % user-defined tile prefix - 'TileTitleSuffix', '^\circ', ... % user-defined suffix (e.g., degrees symbol) - 'SkipSaveFigures', options.skipSaveFigures, ... - 'SaveFileName', 'SpectralCurves.fig', ... - 'SaveDirectory', figSaveDir); - -%% ------------------ 8. Plot of all Angular Spectral Distribution Curves shifted ------------------ -results = Plotter.plotSpectralCurvesRecentered( ... - results_all{1}.results.spectral_analysis_results.S_theta_norm_all, ... - results_all{1}.results.spectral_analysis_results.theta_vals/pi, ... % correct θ values - results_all{1}.scan_reference_values, ... % correct scan params - 'Title', options.titleString, ... - 'XLabel', '\theta / \pi', ... - 'YLabel', 'S(\theta)', ... - 'HighlightEvery', 10, ... % highlight every 10th repetition - 'FontName', options.font, ... - 'FigNum', 21, ... - 'TileTitlePrefix', '\alpha', ... % user-defined tile prefix - 'TileTitleSuffix', '^\circ', ... % user-defined suffix (e.g., degrees symbol) - 'SkipSaveFigures', options.skipSaveFigures, ... - 'SaveFileName', 'SpectralCurves.fig', ... - 'SaveDirectory', figSaveDir); -%% ------------------ 9. Plot cumulants from shifted Angular Spectral Distribution Curves ------------------ -Plotter.plotSpectralDistributionCumulants(results, ... - 'Title', options.titleString, ... - 'XLabel', '\alpha (degrees)', ... - 'FontName', options.font, ... - 'FontSize', 14, ... - 'FigNum', 22, ... - 'SkipSaveFigures', false, ... - 'SaveFileName', 'SpectralCumulants.fig'); - -%% ------------------ 10. Fit shifted Angular Spectral Distribution Curves ------------------ -fitResults = Analyzer.fitTwoGaussianCurves(... - results_all{1}.results.spectral_analysis_results.S_theta_norm_all, ... - results_all{1}.results.spectral_analysis_results.theta_vals, ... - 'MaxTheta', pi/2, ... - 'DeviationLimit', 1.00); - -%% ------------------ 11. Plot fit parameters - position ------------------ -Plotter.plotFitParameterPDF(fitResults, results_all{1}.scan_reference_values, 'mu2', ... - 'Title', options.titleString, ... - 'XLabel', '\alpha (degrees)', ... - 'YLabel', 'Secondary peak position (\theta, rad)', ... - 'FontName', options.font, ... - 'FontSize', 16, ... - 'FigNum', 23, ... - 'SkipSaveFigures', options.skipSaveFigures, ... - 'SaveFileName', 'SecondaryPeakPositionPDF.fig', ... - 'PlotType', 'histogram', ... - 'NumberOfBins', 20, ... - 'NormalizeHistogram', true, ... - 'Colormap', @Colormaps.coolwarm, ... - 'XLim', [min(results_all{1}.scan_reference_values), max(results_all{1}.scan_reference_values)]); - -%% ------------------ 12. Plot fit parameters - width ------------------ -Plotter.plotFitParameterPDF(fitResults, results_all{1}.scan_reference_values, 'sigma2', ... - 'Title', options.titleString, ... - 'XLabel', '\alpha (degrees)', ... - 'YLabel', 'Secondary peak width (\sigma, rad)', ... - 'FontName', options.font, ... - 'FontSize', 16, ... - 'FigNum', 24, ... - 'SkipSaveFigures', options.skipSaveFigures, ... - 'SaveFileName', 'SecondaryPeakWidthPDF.fig', ... - 'PlotType', 'histogram', ... - 'NumberOfBins', 20, ... - 'NormalizeHistogram', true, ... - 'Colormap', @Colormaps.coolwarm, ... - 'XLim', [min(results_all{1}.scan_reference_values), max(results_all{1}.scan_reference_values)]); - -%% ------------------ 13. Plot fit parameters of mean shifted Angular Spectral Distribution Curves ------------------ -S_theta_all = results_all{1}.results.spectral_analysis_results.S_theta_norm_all; -theta_vals = results_all{1}.results.spectral_analysis_results.theta_vals; -scanValues = results_all{1}.scan_reference_values; - -N_params = numel(scanValues); -N_total = numel(S_theta_all); -N_reps = N_total / N_params; - -theta_min = deg2rad(0); -theta_max = deg2rad(90); - -% --- Shift curves so first peak is at start --- -S_theta_all_shifted = cell(size(S_theta_all)); -for k = 1:N_total - curve = S_theta_all{k}; - idx_range = find(theta_vals >= theta_min & theta_vals <= theta_max); - if isempty(idx_range) - [~, peak_idx] = max(curve); - else - [~, local_max_idx] = max(curve(idx_range)); - peak_idx = idx_range(local_max_idx); - end - S_theta_all_shifted{k} = curve(peak_idx:end); -end - -% --- Pad shorter curves with NaN to match lengths --- -Npoints_shifted = max(cellfun(@numel, S_theta_all_shifted)); -for k = 1:N_total - len = numel(S_theta_all_shifted{k}); - if len < Npoints_shifted - S_theta_all_shifted{k} = [S_theta_all_shifted{k}, nan(1, Npoints_shifted-len)]; - end -end - -% --- Compute mean curves per scan parameter --- -meanCurves = cell(1, N_params); -for i = 1:N_params - curves = zeros(N_reps, Npoints_shifted); - for r = 1:N_reps - idx = (r-1)*N_params + i; % interleaved indexing - curves(r,:) = S_theta_all_shifted{idx}; - end - meanCurves{i} = nanmean(curves,1); % mean over repetitions -end - -% --- Fit two-Gaussian model to mean curves --- -fitResultsMean = Analyzer.fitTwoGaussianCurves(meanCurves, theta_vals(1:Npoints_shifted)-theta_vals(1), ... - 'MaxTheta', pi/2, ... - 'DeviationLimit', 1.0); - -% --- Scatter plot of secondary peak position (mu2) vs scan parameter --- -mu2_vals = nan(1, N_params); -sigma2_vals = nan(1, N_params); - -for i = 1:N_params - if fitResultsMean(i).isValid - mu2_vals(i) = fitResultsMean(i).pFit(4); % secondary peak position - sigma2_vals(i) = fitResultsMean(i).pFit(5); % secondary peak width - end -end - -% Secondary peak position -plotSecondaryPeakScatter(fitResultsMean, results_all{1}.scan_reference_values, 'mu2', ... - 'Title', options.titleString, ... - 'XLabel', '\alpha (degrees)', ... - 'YLabel', 'Secondary peak position (\theta, rad)', ... - 'FigNum', 23, ... - 'FontName', options.font, ... - 'FontSize', 16, ... - 'SkipSaveFigures', options.skipSaveFigures, ... - 'SaveFileName', 'SecondaryPeakPositionScatter.fig'); - -% Secondary peak width -plotSecondaryPeakScatter(fitResultsMean, results_all{1}.scan_reference_values, 'sigma2', ... - 'Title', options.titleString, ... - 'XLabel', '\alpha (degrees)', ... - 'YLabel', 'Secondary peak width (\sigma, rad)', ... - 'FigNum', 24, ... - 'FontName', options.font, ... - 'FontSize', 16, ... - 'SkipSaveFigures', options.skipSaveFigures, ... - 'SaveFileName', 'SecondaryPeakWidthScatter.fig'); - - -function plotSecondaryPeakScatter(fitResultsMean, scanValues, parameterName, varargin) -%% plotSecondaryPeakScatter -% Author: Karthik -% Date: 2025-10-06 -% Version: 1.0 -% -% Description: -% Scatter plot of secondary peak fit parameters (mu2 or sigma2) vs scan parameter -% in the same style as plotMeanWithSE. - - % --- 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', [], @(x) isempty(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', 'secondary_peak_scatter.fig', @ischar); - addParameter(p, 'SaveDirectory', pwd, @ischar); - parse(p, varargin{:}); - opts = p.Results; - - % --- Extract parameter values --- - N_params = numel(fitResultsMean); - paramVals = nan(1, N_params); - for i = 1:N_params - if fitResultsMean(i).isValid - switch parameterName - case 'mu2' - paramVals(i) = fitResultsMean(i).pFit(4); - case 'sigma2' - paramVals(i) = fitResultsMean(i).pFit(5); - otherwise - error('Unknown parameter name: %s', parameterName); - end - end - end - - % --- Prepare figure --- - if isempty(opts.FigNum) - fig = figure; - else - fig = figure(opts.FigNum); - clf(fig); - end - set(fig, 'Color', 'w', 'Position', [100 100 950 750]); - - % --- Plot as mean ± SE style (SE=0 here) --- - errorbar(scanValues, paramVals, zeros(size(paramVals)), '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', 'tex', 'FontName', opts.FontName, 'FontSize', opts.FontSize); - ylabel(opts.YLabel, 'Interpreter', 'tex', 'FontName', opts.FontName, 'FontSize', opts.FontSize); - title(opts.Title, 'FontName', opts.FontName, 'FontSize', opts.FontSize + 2, 'FontWeight', 'bold'); - grid on; - - % --- Save figure --- - if ~opts.SkipSaveFigures - if ~exist(opts.SaveDirectory, 'dir'), mkdir(opts.SaveDirectory); end - savefig(fig, fullfile(opts.SaveDirectory, opts.SaveFileName)); - end -end - %{ %% ------------------ 14. Average of Spectra Plots ------------------ Plotter.plotAverageSpectra(scan_parameter_values, ... spectral_analysis_results, ... 'ScanParameterName', scan_parameter, ... - 'FigNum', 25, ... + 'FigNum', 20, ... 'ColormapPS', Colormaps.coolwarm(), ... 'Font', 'Bahnschrift', ... 'SaveFileName', 'avgSpectra.fig', ...