diff --git a/Data-Analyzer/+Analyzer/fitTwoGaussianCurves.m b/Data-Analyzer/+Analyzer/fitTwoGaussianCurves.m index 8277d41..460ea36 100644 --- a/Data-Analyzer/+Analyzer/fitTwoGaussianCurves.m +++ b/Data-Analyzer/+Analyzer/fitTwoGaussianCurves.m @@ -1,4 +1,4 @@ -function fitResults = fitTwoGaussianCurves(S_theta_all, theta_vals, varargin) +function [fitResults, rawCurves] = fitTwoGaussianCurves(S_theta_all, theta_vals, varargin) %% fitTwoGaussianCurves % 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. @@ -13,7 +13,7 @@ function fitResults = fitTwoGaussianCurves(S_theta_all, theta_vals, varargin) % % Optional name-value pairs: % 'MaxTheta' - Maximum theta for search (default: pi/2) -% 'DeviationLimit' - Max relative deviation allowed (default: 0.3) +% 'ResidualThreshold' - 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) @@ -25,7 +25,7 @@ function fitResults = fitTwoGaussianCurves(S_theta_all, theta_vals, varargin) % --- 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, '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)); addParameter(p, 'RecenterCurves', true, @islogical); @@ -34,7 +34,9 @@ function fitResults = fitTwoGaussianCurves(S_theta_all, theta_vals, varargin) Ncurves = numel(S_theta_all); fitResults = struct('pFit',[],'thetaFit',[],'xFit',[],'yFit',[],... - 'thetaFine',[],'isValid',[],'fitMaxTheta',[]); + 'thetaFine',[],'isValid',[],'fitMaxTheta',[]); + + rawCurves = struct('x', cell(1, Ncurves), 'theta', cell(1, Ncurves)); if opts.RecenterCurves % --- Preprocess curves: shift first peak to zero without wrapping --- @@ -69,128 +71,187 @@ function fitResults = fitTwoGaussianCurves(S_theta_all, theta_vals, varargin) end for k = 1:Ncurves + + 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); + + x = x / x(1); + + % --- Store raw data for plotting --- + rawCurves(k).x = x; + rawCurves(k).theta = theta; + try - x = S_theta_all_shifted{k}; - theta = theta_recentered; + % --- Filter 1: Ensure curve covers full 0–MaxTheta range --- + thetaRange = max(theta) - min(theta); + if thetaRange < 0.9 * opts.MaxTheta + warning('Curve %d excluded: does not span full 0–MaxTheta range.', k); + fitResults(k) = fillZeroStruct(theta, opts.MaxTheta); + continue; + end - 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(); + fitResults(k) = fillNaNStruct(); continue; end + + xSmooth = smooth(x, 5, 'moving'); + curveAmp = max(xSmooth) - min(xSmooth); + + %% 1) Attempt primary detection via findpeaks + minProm = opts.AmplitudeThreshold * curveAmp; + [pk, locIdx] = findpeaks(xSmooth, 'MinPeakProminence', minProm); + thetaPeaks = theta(locIdx); + + foundSecondary = false; + fitMaxTheta = []; + secondaryAmp = []; + + secondaryIdx = find(thetaPeaks > opts.PositionThreshold, 1, 'first'); + if ~isempty(secondaryIdx) + foundSecondary = true; + fitMaxTheta = thetaPeaks(secondaryIdx); + secondaryAmp = pk(secondaryIdx); + end + + %% 2) If secondary found, fit constrained two-Gaussian directly + if foundSecondary + 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); + + mu1_guess = 0; + sigma1_guess = 0.1 * fitMaxTheta; + A2_guess = secondaryAmp; + mu2_guess = fitMaxTheta; + sigma2_guess = fitMaxTheta / 3; + p0 = [mu1_guess, sigma1_guess, A2_guess, mu2_guess, sigma2_guess]; + + 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); + + try + pFit = lsqcurvefit(twoGauss, p0, theta, x, lb, ub, optsLSQ); - x = x / x(1); - xSmooth = smooth(x, 5, 'moving'); - curveAmp = max(xSmooth) - min(xSmooth); + % --- Handle missing or weak secondary peak, discard poor fits --- + if pFit(3) <= opts.AmplitudeThreshold + % amplitude too small → discard secondary + pFit(3) = 0; + pFit(4:5) = NaN; + else + % --- Compute RMS of fit vs raw curve in [0, pi/4] --- + maskRMS = theta >= 0 & theta <= pi/4; + if any(maskRMS) + x_raw_subset = x(maskRMS); + y_fit_subset = twoGauss(pFit, theta(maskRMS)); + rms_val = sqrt(mean((x_raw_subset - y_fit_subset).^2)); + if rms_val > opts.ResidualThreshold + % residual too large → discard secondary + pFit(3) = NaN; + pFit(4:5) = NaN; + end + end + end + + % --- Map to canonical format with primary amplitude = 1 --- + pCanonical = [1.0, pFit(1), pFit(2), pFit(3), pFit(4), pFit(5)]; - % --- Robust peak detection --- - [pk, locIdx] = findpeaks(xSmooth, 'MinPeakProminence', opts.AmplitudeThreshold * curveAmp); - thetaPeaks = theta(locIdx); + thetaFine = linspace(0, opts.MaxTheta, 500); + yFit = 1.0*exp(-0.5*((thetaFine - pCanonical(2))/max(pCanonical(3),1e-6)).^2) + ... + pCanonical(4)*exp(-0.5*((thetaFine - pCanonical(5))/max(pCanonical(6),1e-6)).^2); - if isempty(pk) - warning('Curve %d has no significant peaks.', k); - fitResults(k) = fillZeroStruct(theta, opts.MaxTheta); - continue; + fitResults(k).pFit = pCanonical; + fitResults(k).thetaFit = theta; + fitResults(k).xFit = x; + fitResults(k).thetaFine = thetaFine; + fitResults(k).yFit = yFit; + fitResults(k).isValid = true; + fitResults(k).fitMaxTheta = fitMaxTheta; + continue; % success → next curve + catch + warning('Curve %d constrained fit failed, switching to generic fit.', k); + end 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 --- + + %% 3) Generic two-Gaussian fallback (no or failed secondary detection) try - pFit = lsqcurvefit(twoGauss, p0, thetaFit, xFit, lb, ub, optsLSQ); + twoGaussFree = @(p,th) ... + p(1)*exp(-0.5*((th - p(2))/max(p(3),1e-6)).^2) + ... + p(4)*exp(-0.5*((th - p(5))/max(p(6),1e-6)).^2); + + p0 = [1, 0, 0.1, 0.3, opts.MaxTheta/3, 0.1]; + lb = [0, 0, 1e-6, 0, opts.PositionThreshold, 1e-6]; + ub = [2, opts.MaxTheta/2, opts.MaxTheta, 2, opts.MaxTheta, opts.MaxTheta]; + + optsLSQfree = optimoptions('lsqcurvefit', 'Display', 'off', ... + 'MaxFunctionEvaluations', 2e4, 'MaxIterations', 2e4); + + pFree = lsqcurvefit(twoGaussFree, p0, theta(:), xSmooth(:), lb, ub, optsLSQfree); + + % --- Handle missing secondary peak --- + if pFree(4) <= 1e-3 + pFree(4) = 0; % secondary amplitude = 0 + pFree(5:6) = NaN; % secondary position & width = NaN + end + + fitMaxTheta = max(pFree(5), pFree(2)); + thetaFine = linspace(0, opts.MaxTheta, 500); + yFit = twoGaussFree(pFree, thetaFine); + + fitResults(k).pFit = pFree; + fitResults(k).thetaFit = theta; + fitResults(k).xFit = x; + fitResults(k).thetaFine = thetaFine; + fitResults(k).yFit = yFit; + fitResults(k).isValid = true; + fitResults(k).fitMaxTheta = fitMaxTheta; catch - warning('Curve %d fit failed.', k); + warning('Curve %d: generic two-Gaussian 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 completely.', k); - fitResults(k) = fillNaNStruct(); + fitResults(k) = fillNaNStruct(); end 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); + 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); + 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/+Plotter/plotFitParameterPDF.m b/Data-Analyzer/+Plotter/plotFitParameterPDF.m index 2afd49d..4b6a123 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.1 +% Version: 1.2 % % Description: % Plots 2D PDF (heatmap) of any parameter from two-Gaussian fit results @@ -12,7 +12,7 @@ function plotFitParameterPDF(fitResults, scanValues, paramName, varargin) % fitResults - struct array from fitTwoGaussianCurves % scanValues - vector of scan parameter values % paramName - string specifying the parameter to plot: -% 'mu1', 'sigma1', 'A2', 'mu2', 'sigma2' +% 'A1', 'mu1', 'sigma1', 'A2', 'mu2', 'sigma2' % % Optional name-value pairs (same as before, plus OverlayMeanSEM): % 'OverlayMeanSEM' - logical, overlay mean ± SEM (default: true) @@ -31,6 +31,7 @@ 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); @@ -40,9 +41,9 @@ parse(p, varargin{:}); opts = p.Results; % --- Map paramName to index --- -paramMap = struct('mu1',1,'sigma1',2,'A2',3,'mu2',4,'sigma2',5); +paramMap = struct('A1',1,'mu1',2,'sigma1',3,'A2',4,'mu2',5,'sigma2',6); if ~isfield(paramMap,paramName) - error('Invalid paramName. Must be one of: mu1, sigma1, A2, mu2, sigma2'); + error('Invalid paramName. Must be one of: A1, mu1, sigma1, A2, mu2, sigma2'); end paramIdx = paramMap.(paramName); @@ -149,6 +150,10 @@ if ~isempty(opts.XLim) xlim(opts.XLim); end +if ~isempty(opts.YLim) + ylim(opts.YLim); +end + % --- Overlay mean ± SEM if requested --- if opts.OverlayMeanSEM meanParam = nanmean(paramValues,1); diff --git a/Data-Analyzer/+Scripts/BECToDropletsToStripes/runSpectralDistributionAnalysis.m b/Data-Analyzer/+Scripts/BECToDropletsToStripes/runSpectralDistributionAnalysis.m index 0de0e5a..87545ba 100644 --- a/Data-Analyzer/+Scripts/BECToDropletsToStripes/runSpectralDistributionAnalysis.m +++ b/Data-Analyzer/+Scripts/BECToDropletsToStripes/runSpectralDistributionAnalysis.m @@ -165,13 +165,79 @@ Plotter.plotSpectralDistributionCumulants(results, ... 'SaveFileName', 'SpectralCumulants.fig'); %% ------------------ 4. Fit shifted Angular Spectral Distribution Curves ------------------ -fitResults = Analyzer.fitTwoGaussianCurves(... +[fitResults, rawCurves] = Analyzer.fitTwoGaussianCurves(... spectral_analysis_results.S_theta_norm_all, ... spectral_analysis_results.theta_vals, ... 'MaxTheta', pi/2, ... - 'DeviationLimit', 0.30, ... + 'ResidualThreshold', 0.15, ... 'PositionThreshold', pi/15, ... - 'AmplitudeThreshold', 0.02); + 'AmplitudeThreshold', 0.15); + +%{ +% --- Function call --- +plotTwoGaussianFitsOnRaw(fitResults, rawCurves, 8, 12); % 8×12 subplots per page + +% --- Function definition --- +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) +% rawCurves - struct array with fields: +% .x - raw curve +% .theta - corresponding theta values +% nRows - number of subplot rows per page (default: 4) +% nCols - number of subplot columns per page (default: 6) + + if nargin < 3, nRows = 4; end + if nargin < 4, nCols = 6; end + + Ncurves = numel(rawCurves); + plotsPerPage = nRows * nCols; + pageNum = 1; + + for k = 1:Ncurves + % --- New figure/page if needed --- + if mod(k-1, plotsPerPage) == 0 + figure('Name', sprintf('Curves Page %d', pageNum), ... + 'NumberTitle', 'off', 'Color', 'w'); + pageNum = pageNum + 1; + end + + % --- Select subplot --- + subplot(nRows, nCols, mod(k-1, plotsPerPage)+1); + hold on; grid on; box on; + + % --- Plot raw curve --- + xRaw = rawCurves(k).x; + thetaRaw = rawCurves(k).theta; + if ~isempty(xRaw) && ~isempty(thetaRaw) && all(isfinite(xRaw)) + plot(thetaRaw, xRaw, 'k.-', 'LineWidth', 1, 'DisplayName', 'Raw data'); + end + + % --- Overlay fit if valid --- + if k <= numel(fitResults) + fit = fitResults(k); + if isfield(fit, 'isValid') && fit.isValid ... + && ~isempty(fit.yFit) && ~isempty(fit.thetaFine) + plot(fit.thetaFine, fit.yFit, 'r-', 'LineWidth', 1.2, 'DisplayName', 'Two-Gaussian fit'); + end + end + + % --- Formatting --- + xlabel('\theta (rad)'); + ylabel('Normalized amplitude'); + title(sprintf('Curve %d', k), 'FontSize', 10); + + % --- Legend once per page --- + if mod(k-1, plotsPerPage) == 0 + legend('Location', 'best', 'FontSize', 8); + end + + hold off; + end +end +%} %% ------------------ 5. Plot fit parameters - amplitude ------------------ Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'A2', ... @@ -187,7 +253,8 @@ Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'A2', ... 'NumberOfBins', 20, ... 'NormalizeHistogram', true, ... 'Colormap', @Colormaps.coolwarm, ... - 'XLim', [min(scan_reference_values), max(scan_reference_values)]); + 'XLim', [min(scan_reference_values), max(scan_reference_values)], ... + 'YLim', [0, 1.6]); %% ------------------ 6. Plot fit parameters - position ------------------ Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'mu2', ... @@ -203,7 +270,8 @@ Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'mu2', ... 'NumberOfBins', 20, ... 'NormalizeHistogram', true, ... 'Colormap', @Colormaps.coolwarm, ... - 'XLim', [min(scan_reference_values), max(scan_reference_values)]); + 'XLim', [min(scan_reference_values), max(scan_reference_values)], ... + 'YLim', [0, 1.6]); %% ------------------ 7. Plot fit parameters - width ------------------ Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'sigma2', ... @@ -219,7 +287,8 @@ Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'sigma2', ... 'NumberOfBins', 20, ... 'NormalizeHistogram', true, ... 'Colormap', @Colormaps.coolwarm, ... - 'XLim', [min(scan_reference_values), max(scan_reference_values)]); + 'XLim', [min(scan_reference_values), max(scan_reference_values)], ... + 'YLim', [0, 1.6]); %% ------------------ 8. Plot fit parameters of mean shifted Angular Spectral Distribution Curves ------------------ @@ -242,13 +311,13 @@ for i = 1:numel(results.curves) end % --- Fit two-Gaussian model to mean curves --- -fitResultsMean = Analyzer.fitTwoGaussianCurves(... +[fitResultsMean, ~] = Analyzer.fitTwoGaussianCurves(... results.curves_mean, ... results.x_values*pi, ... 'MaxTheta', pi/2, ... - 'DeviationLimit', 0.30, ... + 'ResidualThreshold', 0.15, ... 'PositionThreshold', pi/15, ... - 'AmplitudeThreshold', 0.02, ... + 'AmplitudeThreshold', 0.15, ... 'RecenterCurves', false); % --- Prepare parameter values --- @@ -259,18 +328,11 @@ 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)) + if all(~isnan(pFit)) % Successful fit → use fitted values - amp2_vals(i) = pFit(3); - mu2_vals(i) = pFit(4); - sigma2_vals(i) = pFit(5); + amp2_vals(i) = pFit(4); + mu2_vals(i) = pFit(5); + sigma2_vals(i) = pFit(6); else % Fit failed → leave as NaN (skipped automatically) @@ -326,17 +388,17 @@ for i = 1:numel(results.curves) 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 +paramIdx = 1; % <-- 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( ... +[fitResults, rawCurves] = Analyzer.fitTwoGaussianCurves(... curves_cell, ... results.x_values*pi, ... - 'MaxTheta', pi/2, ... - 'DeviationLimit', 0.30, ... + 'MaxTheta', pi/2, ... + 'ResidualThreshold', 0.15, ... 'PositionThreshold', pi/15, ... - 'AmplitudeThreshold', 0.02, ... + 'AmplitudeThreshold', 0.15, ... 'RecenterCurves', false); diff --git a/Data-Analyzer/+Scripts/BECToStripesToDroplets/runSpectralDistributionAnalysis.m b/Data-Analyzer/+Scripts/BECToStripesToDroplets/runSpectralDistributionAnalysis.m new file mode 100644 index 0000000..75aba95 --- /dev/null +++ b/Data-Analyzer/+Scripts/BECToStripesToDroplets/runSpectralDistributionAnalysis.m @@ -0,0 +1,397 @@ +%% ===== BEC-Stripes-Droplets Settings ===== + +% Specify data location to run analysis on +dataSources = { + struct('sequence', 'TwoDGas', ... + 'date', '2025/06/24', ... + 'runs', [1]) % 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 = 'StripesToDroplets'; +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 = 24.6; +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 = true; +options.showProgressBar = true; + +% Extras +options.font = 'Bahnschrift'; +switch options.measurementName + case 'BECToDroplets' + options.scan_parameter = 'rot_mag_field'; + options.flipSortOrder = false; + options.scanParameterUnits = 'gauss'; + options.titleString = 'BEC to Droplets'; + case 'BECToStripes' + options.scan_parameter = 'rot_mag_field'; + options.flipSortOrder = false; + 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, rawCurves] = Analyzer.fitTwoGaussianCurves(... + spectral_analysis_results.S_theta_norm_all, ... + spectral_analysis_results.theta_vals, ... + 'MaxTheta', pi/2, ... + 'ResidualThreshold', 0.15, ... + 'PositionThreshold', pi/15, ... + 'AmplitudeThreshold', 0.15); + +%{ +% --- Function call --- +plotTwoGaussianFitsOnRaw(fitResults, rawCurves, 8, 12); % 8×12 subplots per page + +% --- Function definition --- +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) +% rawCurves - struct array with fields: +% .x - raw curve +% .theta - corresponding theta values +% nRows - number of subplot rows per page (default: 4) +% nCols - number of subplot columns per page (default: 6) + + if nargin < 3, nRows = 4; end + if nargin < 4, nCols = 6; end + + Ncurves = numel(rawCurves); + plotsPerPage = nRows * nCols; + pageNum = 1; + + for k = 1:Ncurves + % --- New figure/page if needed --- + if mod(k-1, plotsPerPage) == 0 + figure('Name', sprintf('Curves Page %d', pageNum), ... + 'NumberTitle', 'off', 'Color', 'w'); + pageNum = pageNum + 1; + end + + % --- Select subplot --- + subplot(nRows, nCols, mod(k-1, plotsPerPage)+1); + hold on; grid on; box on; + + % --- Plot raw curve --- + xRaw = rawCurves(k).x; + thetaRaw = rawCurves(k).theta; + if ~isempty(xRaw) && ~isempty(thetaRaw) && all(isfinite(xRaw)) + plot(thetaRaw, xRaw, 'k.-', 'LineWidth', 1, 'DisplayName', 'Raw data'); + end + + % --- Overlay fit if valid --- + if k <= numel(fitResults) + fit = fitResults(k); + if isfield(fit, 'isValid') && fit.isValid ... + && ~isempty(fit.yFit) && ~isempty(fit.thetaFine) + plot(fit.thetaFine, fit.yFit, 'r-', 'LineWidth', 1.2, 'DisplayName', 'Two-Gaussian fit'); + end + end + + % --- Formatting --- + xlabel('\theta (rad)'); + ylabel('Normalized amplitude'); + title(sprintf('Curve %d', k), 'FontSize', 10); + + % --- Legend once per page --- + if mod(k-1, plotsPerPage) == 0 + legend('Location', 'best', 'FontSize', 8); + end + + hold off; + end +end +%} + +%% ------------------ 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)], ... + 'YLim', [0, 1.6]); + +%% ------------------ 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)], ... + 'YLim', [0, 1.6]); + +%% ------------------ 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)], ... + 'YLim', [0, 1.6]); + +%% ------------------ 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, ... + 'ResidualThreshold', 0.15, ... + 'PositionThreshold', pi/15, ... + 'AmplitudeThreshold', 0.15, ... + '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(~isnan(pFit)) + % Successful fit → use fitted values + amp2_vals(i) = pFit(4); + mu2_vals(i) = pFit(5); + sigma2_vals(i) = pFit(6); + + 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 = 1; % <-- 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 --- +[fitResults, rawCurves] = Analyzer.fitTwoGaussianCurves(... + curves_cell, ... + results.x_values*pi, ... + 'MaxTheta', pi/2, ... + 'ResidualThreshold', 0.15, ... + 'PositionThreshold', pi/15, ... + 'AmplitudeThreshold', 0.15, ... + 'RecenterCurves', false);