Calculations/Data-Analyzer/+Scripts/BECToStripes/runRadialSpectralDistributionAnalysis.m

448 lines
18 KiB
Matlab
Raw Permalink Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

%% ===== BEC-Stripes Settings =====
% Specify data location to run analysis on
%%
dataSources = {
struct('sequence', 'StructuralPhaseTransition', ...
'date', '2025/10/11', ...
'runs', [1]) % specify run numbers as a string in "" or just as a numeric value
};
%%
%{
dataSources = {
struct('sequence', 'StructuralPhaseTransition', ...
'date', '2025/10/10', ...
'runs', [2]) % 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/202510';
options.measurementName = 'BECToStripes';
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 = [1420, 2040];
options.span = [200, 200];
options.fraction = [0.1, 0.1];
options.pixel_size = 5.86e-6; % in meters
options.magnification = 23.2;
options.ImagingMode = 'HighIntensity';
options.PulseDuration = 5e-6; % in s
% Fourier analysis settings
options.theta_min = deg2rad(0);
options.theta_max = deg2rad(180);
options.N_radial_bins = 500;
options.Radial_Sigma = 2;
options.Radial_WindowSize = 5; % odd number
options.k_min = 1.2771; % μm⁻¹
options.k_max = 2.5541; % μm⁻¹
options.N_angular_bins = 180;
options.Angular_Threshold = 75;
options.Angular_Sigma = 2;
options.Angular_WindowSize = 5;
options.zoom_size = 50;
% 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.extractFullRadialSpectralDistribution(od_imgs, options);
%% ------------------ 1. Plot of all Radial Spectral Distribution Curves ------------------
Plotter.plotSpectralCurves( ...
spectral_analysis_results.S_k_all, ...
spectral_analysis_results.k_rho_vals, ...
scan_reference_values, ... % correct scan params
'Title', options.titleString, ...
'XLabel', 'k_\rho', ...
'YLabel', 'S(k_\rho)', ...
'YScale', 'log', ...
'XLim', [min(spectral_analysis_results.k_rho_vals), max(spectral_analysis_results.k_rho_vals)], ...
'HighlightEvery', 10, ... % highlight every 10th repetition
'FontName', options.font, ...
'FigNum', 1, ...
'TileTitlePrefix', 'B', ... % user-defined tile prefix
'TileTitleSuffix', 'G', ... % user-defined suffix (e.g., degrees symbol)
'SkipSaveFigures', options.skipSaveFigures, ...
'SaveFileName', 'SpectralCurves.fig', ...
'SaveDirectory', options.saveDirectory);
%% ------------------ 2. Fit truncated Radial Spectral Distribution Curves ------------------
[fitResults, rawCurves] = Analyzer.fitTwoGaussianCurvesToRadialSpectralDistribution(...
spectral_analysis_results.S_k_all, ... % radial spectral curves
spectral_analysis_results.k_rho_vals, ... % corresponding k_rho values
'KRhoRange', [0, 3], ... % truncate curves to 0 <= k_rho <= 3
'AmplitudeRange', [0, 0.3], ... % truncate curves to y >= 0 (all amplitudes)
'ResidualThreshold', 0.175, ... % maximum allowed residual
'PositionThreshold', 0.5, ... % minimum separation between peaks
'AmplitudeThreshold', 0.015); ... % minimum secondary peak fraction
%% Post-fit diagnostics
% --- Determine max number of parameters across all fits ---
maxParams = max(cellfun(@numel, {fitResults.pFit}));
% --- Preallocate padded array ---
Ncurves = numel(fitResults);
pFit_all = nan(Ncurves, maxParams);
% --- Fill padded array ---
for k = 1:Ncurves
pk = fitResults(k).pFit;
pFit_all(k,1:numel(pk)) = pk;
end
isValid_all = [fitResults.isValid]; % logical vector
% --- Extract parameters of the SECOND Gaussian only ---
% For both 2G and 3G, second Gaussian = indices 4:6
pSecond = pFit_all(:, 4:6);
% --- Diagnostics ---
numWithNaN = sum(any(isnan(pSecond), 2));
numWithZero = sum(any(pSecond == 0, 2));
numValid = sum(isValid_all);
% Display summary
fprintf('\n--- Fit diagnostics (second Gaussian only) ---\n');
fprintf('Curves with ≥1 NaN in second Gaussian params: %d\n', numWithNaN);
fprintf('Curves with ≥1 zero in second Gaussian params: %d\n', numWithZero);
fprintf('Curves marked as isValid = true: %d\n', numValid);
fprintf('----------------------------------------\n\n');
%% ------------------ Plot Fits on Raw ------------------
options.skipLivePlot = true;
Plotter.plotTwoModeGaussianFitsOnRawRSD(fitResults, rawCurves, 4, 6, ...
'SkipLivePlot', options.skipLivePlot, ...
'SkipSaveFigures', options.skipSaveFigures, ...
'SaveDirectory', options.saveDirectory);
%% ------------------ Interactive Fit Viewer ------------------
Analyzer.runInteractiveRSDTwoGaussianFitGUI(spectral_analysis_results);
%% ------------------ 3. Plot fit parameters - amplitude ------------------
Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'A2', ...
'Title', options.titleString, ...
'XLabel', 'B (G)', ...
'YLabel', 'Secondary peak amplitude', ...
'FontName', options.font, ...
'FontSize', 16, ...
'FigNum', 2, ...
'SkipSaveFigures', options.skipSaveFigures, ...
'SaveFileName', 'SecondaryPeakAmplitudePDF.fig', ...
'SaveDirectory', options.saveDirectory, ...
'PlotType', 'histogram', ...
'NumberOfBins', 20, ...
'NormalizeHistogram', true, ...
'Colormap', @Colormaps.coolwarm, ...
'XLim', [min(scan_reference_values), max(scan_reference_values)], ...
'YLim', [0, 0.01]);
%% ------------------ 4. Plot fit parameters - position ------------------
Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'mu2', ...
'Title', options.titleString, ...
'XLabel', 'B (G)', ...
'YLabel', 'Secondary peak position (k_\rho, \mum^{-1})', ...
'FontName', options.font, ...
'FontSize', 16, ...
'FigNum', 3, ...
'SkipSaveFigures', options.skipSaveFigures, ...
'SaveFileName', 'SecondaryPeakPositionPDF.fig', ...
'SaveDirectory', options.saveDirectory, ...
'PlotType', 'histogram', ...
'NumberOfBins', 20, ...
'NormalizeHistogram', true, ...
'Colormap', @Colormaps.coolwarm, ...
'XLim', [min(scan_reference_values), max(scan_reference_values)], ...
'YLim', [0, 3]);
%% ------------------ 5. Plot fit parameters - width ------------------
Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'sigma2', ...
'Title', options.titleString, ...
'XLabel', 'B (G)', ...
'YLabel', 'Secondary peak width (\sigma, rad)', ...
'FontName', options.font, ...
'FontSize', 16, ...
'FigNum', 4, ...
'SkipSaveFigures', options.skipSaveFigures, ...
'SaveFileName', 'SecondaryPeakWidthPDF.fig', ...
'SaveDirectory', options.saveDirectory, ...
'PlotType', 'histogram', ...
'NumberOfBins', 20, ...
'NormalizeHistogram', true, ...
'Colormap', @Colormaps.coolwarm, ...
'XLim', [min(scan_reference_values), max(scan_reference_values)], ...
'YLim', [0, 1.5]);
%% ------------------ 6. Plot secondary Gaussian range ------------------
Plotter.plotSecondaryGaussianRange(fitResults, scan_reference_values, ...
'Title', options.titleString, ...
'XLabel', 'B (G)', ...
'YLabel', '\mu_2 \pm \sigma_2 (rad)', ...
'YLim', [0, 5], ...
'FontName', options.font, ...
'FontSize', 16, ...
'FigNum', 5, ...
'NumberOfBins', 20, ...
'NormalizeHistogram', true, ...
'Colormap', @Colormaps.coolwarm, ...
'OverlayMeanSEM', true, ...
'SkipSaveFigures', options.skipSaveFigures, ...
'SaveFileName', 'SecondaryGaussianRange.fig', ...
'SaveDirectory', options.saveDirectory);
%% ------------------ 7. Plot fit parameters of mean shifted Angular Spectral Distribution Curves ------------------
avgData = Analyzer.computeSpectralAverages(spectral_analysis_results.S_k_all, scan_reference_values);
% Extract averaged curves (mean across repetitions)
curves_mean = avgData.curves_mean; % [N_scanValues × N_k_rho]
krho_values = spectral_analysis_results.k_rho_vals;
% ---------- Fit two-Gaussian model to mean curves ----------
[fitResultsMean, ~] = Analyzer.fitTwoGaussianCurvesToRadialSpectralDistribution(...
curves_mean, ...
krho_values, ...
'KRhoRange', [0, 3], ... % truncate curves to 0 <= k_rho <= 3
'AmplitudeRange', [0, 0.3], ... % truncate curves to y >= 0 (all amplitudes)
'ResidualThreshold', 0.175, ... % maximum allowed residual
'PositionThreshold', 0.5, ... % minimum separation between peaks
'AmplitudeThreshold', 0.015); % minimum secondary peak fraction
% ---------- Extract fit parameters ----------
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))
amp2_vals(i) = pFit(4);
mu2_vals(i) = pFit(5);
sigma2_vals(i) = pFit(6);
end
end
% ---------- Plot with SEM ----------
Plotter.plotMeanWithSE(scan_reference_values, amp2_vals, ...
'Title', options.titleString, ...
'XLabel', 'B (G)', ...
'YLabel', 'Secondary peak amplitude', ...
'FigNum', 6, ...
'FontName', options.font, ...
'FontSize', 16, ...
'SkipSaveFigures', options.skipSaveFigures, ...
'SaveFileName', 'SecondaryPeakAmplitudeMeanWithSEM.fig', ...
'SaveDirectory', options.saveDirectory);
Plotter.plotMeanWithSE(scan_reference_values, mu2_vals, ...
'Title', options.titleString, ...
'XLabel', 'B (G)', ...
'YLabel', 'Secondary peak position (k_\rho, \mum^{-1})', ...
'FigNum', 7, ...
'FontName', options.font, ...
'FontSize', 16, ...
'SkipSaveFigures', options.skipSaveFigures, ...
'SaveFileName', 'SecondaryPeakPositionMeanWithSEM.fig', ...
'SaveDirectory', options.saveDirectory);
Plotter.plotMeanWithSE(scan_reference_values, sigma2_vals, ...
'Title', options.titleString, ...
'XLabel', 'B (G)', ...
'YLabel', 'Secondary peak width (\sigma, rad)', ...
'FigNum', 8, ...
'FontName', options.font, ...
'FontSize', 16, ...
'SkipSaveFigures', options.skipSaveFigures, ...
'SaveFileName', 'SecondaryPeakWidthMeanWithSEM.fig', ...
'SaveDirectory', options.saveDirectory);
%% ------------------ 8. Compare mean of individual fits vs fit of mean ------------------
paramMap = struct('A1',1,'mu1',2,'sigma1',3,'A2',4,'mu2',5,'sigma2',6);
paramNames = {'A2', 'mu2', 'sigma2'};
paramLabels = {'Secondary peak amplitude', 'Secondary peak position (k_\rho, \mum^{-1})', 'Secondary peak width (\sigma, rad)'};
meanOfIndividualFits = struct();
fitOfMeanCurve = struct();
for p = 1:numel(paramNames)
paramName = paramNames{p};
paramLabel = paramLabels{p};
paramIdx = paramMap.(paramName);
% --- Mean of individual fits per scan_reference_values ---
meanVals = nan(size(scan_reference_values)); % already unique
for k = 1:numel(scan_reference_values)
thisAlpha = scan_reference_values(k);
mask = scan_parameter_values == thisAlpha; % pick all repetitions
pvals = arrayfun(@(f) f.pFit(paramIdx), fitResults(mask), 'UniformOutput', true);
meanVals(k) = mean(pvals(~isnan(pvals) & isfinite(pvals)), 'omitnan');
end
meanOfIndividualFits.(paramName) = meanVals;
% --- Fit of mean curve ---
switch paramName
case 'A2'
fitOfMeanCurve.(paramName) = amp2_vals;
case 'mu2'
fitOfMeanCurve.(paramName) = mu2_vals;
case 'sigma2'
fitOfMeanCurve.(paramName) = sigma2_vals;
end
% --- Plot comparison ---
fig = figure('Name', 'Comparison', 'NumberTitle', 'off');
set(fig, 'Color', 'w', 'Position', [100 100 950 750]);
hold on;
plot(scan_reference_values, meanOfIndividualFits.(paramName), 'o-', 'LineWidth', 1.8, ...
'DisplayName', 'Mean of Individual Fits');
plot(scan_reference_values, fitOfMeanCurve.(paramName), 's--', 'LineWidth', 1.8, ...
'DisplayName', 'Fit of Mean Curve');
hold off;
xlabel('B (G)', 'FontName', options.font, 'FontSize', 16);
ylabel(paramLabel, 'FontName', options.font, 'FontSize', 16);
title(options.titleString, 'FontName', options.font, 'FontSize', 16);
legend('Location', 'northeast');
grid on;
set(gca, 'FontSize', 14);
if ~options.skipSaveFigures
saveas(fig, fullfile(options.saveDirectory, ['Compare_' paramName '.fig']));
exportgraphics(fig, fullfile(options.saveDirectory, ['Compare_' paramName '.png']), 'Resolution', 300);
end
end
%% ------------------ 9. Save Fit Results ------------------
outputDir = fullfile(options.saveDirectory, 'BECToS_RSDFitData');
if ~exist(outputDir, 'dir')
mkdir(outputDir);
end
% --- Define filenames for each variable ---
file_scanParamVals = fullfile(outputDir, 'scan_parameter_values');
file_scanRefVals = fullfile(outputDir, 'scan_reference_values');
file_fitResults = fullfile(outputDir, 'fitResults.mat');
file_rawCurves = fullfile(outputDir, 'rawCurves.mat');
file_spectralResults = fullfile(outputDir, 'spectral_analysis_results.mat');
file_meanFits = fullfile(outputDir, 'fitOfMeanCurve.mat');
file_indivFits = fullfile(outputDir, 'meanOfIndividualFits.mat');
% --- Save each variable separately ---
save(file_scanParamVals, 'scan_parameter_values', '-v7.3');
save(file_scanRefVals, 'scan_reference_values', '-v7.3');
save(file_fitResults, 'fitResults', '-v7.3');
save(file_rawCurves, 'rawCurves', '-v7.3');
save(file_spectralResults, 'spectral_analysis_results', '-v7.3');
save(file_meanFits, 'fitOfMeanCurve', '-v7.3');
save(file_indivFits, 'meanOfIndividualFits', '-v7.3');
fprintf('\n--- Fit data saved successfully ---\n');
fprintf('scanParamVals: %s\n', file_scanParamVals);
fprintf('scanRefVals: %s\n', file_scanRefVals);
fprintf('fitResults: %s\n', file_fitResults);
fprintf('rawCurves: %s\n', file_rawCurves);
fprintf('spectral_analysis_results:%s\n', file_spectralResults);
fprintf('fitOfMeanCurve: %s\n', file_meanFits);
fprintf('meanOfIndividualFits: %s\n', file_indivFits);
fprintf('----------------------------------\n\n');
%{
%% ------------------ Single Gaussian Fit to truncated Radial Spectral Distribution Curves ------------------
[fitResults, rawCurves] = Analyzer.fitSingleGaussianCurvesToRadialSpectralDistribution(...
spectral_analysis_results.S_k_all, ... % radial spectral curves
spectral_analysis_results.k_rho_vals, ... % corresponding k_rho values
'KRhoRange', [0.5, 3], ... % truncate curves to 0 <= k_rho <= 3
'ResidualThreshold', 0.15); % maximum allowed residual
Plotter.plotSingleGaussianFitParameterPDF(fitResults, scan_reference_values, 'A', ...
'Title', options.titleString, ...
'XLabel', 'B (G)', ...
'YLabel', 'Peak amplitude', ...
'FontName', options.font, ...
'FontSize', 16, ...
'FigNum', 2, ...
'SkipSaveFigures', options.skipSaveFigures, ...
'SaveFileName', 'PeakAmplitudePDF.fig', ...
'SaveDirectory', options.saveDirectory, ...
'PlotType', 'histogram', ...
'NumberOfBins', 20, ...
'NormalizeHistogram', true, ...
'Colormap', @Colormaps.coolwarm, ...
'XLim', [min(scan_reference_values), max(scan_reference_values)], ...
'YLim', [0, 0.01]);
Plotter.plotSingleGaussianFitParameterPDF(fitResults, scan_reference_values, 'mu', ...
'Title', options.titleString, ...
'XLabel', 'B (G)', ...
'YLabel', 'Peak position', ...
'FontName', options.font, ...
'FontSize', 16, ...
'FigNum', 3, ... % changed figure number to avoid overlap
'SkipSaveFigures', options.skipSaveFigures, ...
'SaveFileName', 'PeakPositionPDF.fig', ...
'SaveDirectory', options.saveDirectory, ...
'PlotType', 'histogram', ...
'NumberOfBins', 20, ...
'NormalizeHistogram', true, ...
'Colormap', @Colormaps.coolwarm, ...
'XLim', [min(scan_reference_values), max(scan_reference_values)], ...
'YLim', [0, 3]);
%}