Minor modifications to scripts.

Data-Analyzer/+Scripts/BECToDroplets/runAngularSpectralDistributionAnalysis.m

@@ -254,7 +254,6 @@ Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'sigma2', ...
     '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, ...
@@ -339,6 +338,98 @@ Plotter.plotMeanWithSE(scan_reference_values, sigma2_vals, ...
     'SaveFileName', 'SecondaryPeakWidthMeanWithSEM.fig', ...
     'SaveDirectory', options.saveDirectory);
 
+%% ------------------ 9. 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 (\theta, rad)', '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;
+    ylim([0, 1.1])
+    xlabel('\alpha (degrees)', '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
+
+%% ------------------ 10. Save Fit Results ------------------
+outputDir = fullfile(options.saveDirectory, 'BECToD_ASDFitData');
+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');
+
 %% Inspect individual realizations of a single parameter
 
 % --- Recenter curves first ---

Data-Analyzer/+Scripts/BECToDroplets/runRadialSpectralDistributionAnalysis.m

@@ -122,10 +122,31 @@ Plotter.plotSpectralCurves( ...
     'PositionThreshold', 0.5, ...                 % minimum separation between peaks
     'AmplitudeThreshold', 0.015);                 % minimum secondary peak fraction
 
+% --- Post-fit diagnostics ---
+pFit_all    = vertcat(fitResults.pFit);   % Ncurves x Nparams
+isValid_all = [fitResults.isValid];       % logical vector
+
+% Extract last three parameters (usually secondary peak amplitude, position, width)
+pTail       = pFit_all(:, end-2:end);
+
+% Count curves where *any* of these last three values are NaN or zero
+numWithNaN  = sum(any(isnan(pTail), 2));
+numWithZero = sum(any(pTail == 0, 2));
+
+% Count of valid fits
+numValid    = sum(isValid_all);
+
+% Display summary
+fprintf('\n--- Fit diagnostics (last 3 params) ---\n');
+fprintf('Curves with ≥1 NaN in last three pFit elements: %d\n', numWithNaN);
+fprintf('Curves with ≥1 zero in last three pFit elements: %d\n', numWithZero);
+fprintf('Curves marked as isValid = true: %d\n', numValid);
+fprintf('----------------------------------------\n\n');
+
 %% ------------------ Plot Fits on Raw ------------------
 options.skipLivePlot    = true;
 options.skipSaveFigures = false; 
-Plotter.plotTwoModeGaussianFitsOnRaw(fitResults, rawCurves, 4, 6, ...
+Plotter.plotTwoModeGaussianFitsOnRawRSD(fitResults, rawCurves, 4, 6, ...
     'SkipLivePlot', options.skipLivePlot, ...
     'SkipSaveFigures', options.skipSaveFigures, ...
     'SaveDirectory', options.saveDirectory);
@@ -140,7 +161,7 @@ Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'A2', ...
     'YLabel', 'Secondary peak amplitude', ...
     'FontName', options.font, ...
     'FontSize', 16, ...
-    'FigNum', 4, ...
+    'FigNum', 2, ...
     'SkipSaveFigures', options.skipSaveFigures, ...
     'SaveFileName', 'SecondaryPeakAmplitudePDF.fig', ...
     'SaveDirectory', options.saveDirectory, ...
@@ -158,7 +179,7 @@ Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'mu2', ...
     'YLabel', 'Secondary peak position (\theta, rad)', ...
     'FontName', options.font, ...
     'FontSize', 16, ...
-    'FigNum', 5, ...
+    'FigNum', 3, ...
     'SkipSaveFigures', options.skipSaveFigures, ...
     'SaveFileName', 'SecondaryPeakPositionPDF.fig', ...
     'SaveDirectory', options.saveDirectory, ...
@@ -175,7 +196,7 @@ Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'sigma2', ...
     'YLabel', 'Secondary peak width (\sigma, rad)', ...
     'FontName', options.font, ...
     'FontSize', 16, ...
-    'FigNum', 6, ...
+    'FigNum', 4, ...
     'SkipSaveFigures', options.skipSaveFigures, ...
     'SaveFileName', 'SecondaryPeakWidthPDF.fig', ...
     'SaveDirectory', options.saveDirectory, ...
@@ -193,7 +214,7 @@ Plotter.plotSecondaryGaussianRange(fitResults, scan_reference_values, ...
     'YLabel', '\mu_2 \pm \sigma_2 (rad)', ...
     'FontName', options.font, ...
     'FontSize', 16, ...
-    'FigNum', 8, ...
+    'FigNum', 5, ...
     'NumberOfBins', 20, ...
     'NormalizeHistogram', true, ...
     'Colormap', @Colormaps.coolwarm, ...
@@ -202,4 +223,159 @@ Plotter.plotSecondaryGaussianRange(fitResults, scan_reference_values, ...
     '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.15, ...                % 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 (\theta, rad)', ...
+    '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 (\theta, rad)', '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('\alpha (degrees)', '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, 'BECToD_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');

Data-Analyzer/+Scripts/BECToDropletsToStripes/runAngularSpectralDistributionAnalysis.m

@@ -340,7 +340,6 @@ Plotter.plotMeanWithSE(scan_reference_values, sigma2_vals, ...
     'SaveDirectory', options.saveDirectory);
 
 %% ------------------ 9. 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 (\theta, rad)', 'Secondary peak width (\sigma, rad)'};
@@ -398,8 +397,8 @@ for p = 1:numel(paramNames)
     end
 end
 
-%% ------------------ 5. Save Fit Results ------------------
+%% ------------------ 10. Save Fit Results ------------------
-outputDir = fullfile(options.saveDirectory, 'DToS_FitData');
+outputDir = fullfile(options.saveDirectory, 'DToS_ASDFitData');
 if ~exist(outputDir, 'dir')
     mkdir(outputDir);
 end

Data-Analyzer/+Scripts/BECToDropletsToStripes/runRadialSpectralDistributionAnalysis.m

@@ -122,10 +122,31 @@ Plotter.plotSpectralCurves( ...
     'PositionThreshold', 0.5, ...                 % minimum separation between peaks
     'AmplitudeThreshold', 0.015);                 % minimum secondary peak fraction
 
+% --- Post-fit diagnostics ---
+pFit_all    = vertcat(fitResults.pFit);   % Ncurves x Nparams
+isValid_all = [fitResults.isValid];       % logical vector
+
+% Extract last three parameters (usually secondary peak amplitude, position, width)
+pTail       = pFit_all(:, end-2:end);
+
+% Count curves where *any* of these last three values are NaN or zero
+numWithNaN  = sum(any(isnan(pTail), 2));
+numWithZero = sum(any(pTail == 0, 2));
+
+% Count of valid fits
+numValid    = sum(isValid_all);
+
+% Display summary
+fprintf('\n--- Fit diagnostics (last 3 params) ---\n');
+fprintf('Curves with ≥1 NaN in last three pFit elements: %d\n', numWithNaN);
+fprintf('Curves with ≥1 zero in last three pFit elements: %d\n', numWithZero);
+fprintf('Curves marked as isValid = true: %d\n', numValid);
+fprintf('----------------------------------------\n\n');
+
 %% ------------------ Plot Fits on Raw ------------------
 options.skipLivePlot    = true;
 options.skipSaveFigures = false; 
-Plotter.plotTwoModeGaussianFitsOnRaw(fitResults, rawCurves, 4, 6, ...
+Plotter.plotTwoModeGaussianFitsOnRawRSD(fitResults, rawCurves, 4, 6, ...
     'SkipLivePlot', options.skipLivePlot, ...
     'SkipSaveFigures', options.skipSaveFigures, ...
     'SaveDirectory', options.saveDirectory);
@@ -140,7 +161,7 @@ Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'A2', ...
     'YLabel', 'Secondary peak amplitude', ...
     'FontName', options.font, ...
     'FontSize', 16, ...
-    'FigNum', 4, ...
+    'FigNum', 2, ...
     'SkipSaveFigures', options.skipSaveFigures, ...
     'SaveFileName', 'SecondaryPeakAmplitudePDF.fig', ...
     'SaveDirectory', options.saveDirectory, ...
@@ -158,7 +179,7 @@ Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'mu2', ...
     'YLabel', 'Secondary peak position (\theta, rad)', ...
     'FontName', options.font, ...
     'FontSize', 16, ...
-    'FigNum', 5, ...
+    'FigNum', 3, ...
     'SkipSaveFigures', options.skipSaveFigures, ...
     'SaveFileName', 'SecondaryPeakPositionPDF.fig', ...
     'SaveDirectory', options.saveDirectory, ...
@@ -175,7 +196,7 @@ Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'sigma2', ...
     'YLabel', 'Secondary peak width (\sigma, rad)', ...
     'FontName', options.font, ...
     'FontSize', 16, ...
-    'FigNum', 6, ...
+    'FigNum', 4, ...
     'SkipSaveFigures', options.skipSaveFigures, ...
     'SaveFileName', 'SecondaryPeakWidthPDF.fig', ...
     'SaveDirectory', options.saveDirectory, ...
@@ -193,7 +214,7 @@ Plotter.plotSecondaryGaussianRange(fitResults, scan_reference_values, ...
     'YLabel', '\mu_2 \pm \sigma_2 (rad)', ...
     'FontName', options.font, ...
     'FontSize', 16, ...
-    'FigNum', 8, ...
+    'FigNum', 5, ...
     'NumberOfBins', 20, ...
     'NormalizeHistogram', true, ...
     'Colormap', @Colormaps.coolwarm, ...
@@ -202,4 +223,159 @@ Plotter.plotSecondaryGaussianRange(fitResults, scan_reference_values, ...
     '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.15, ...                % 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', '\alpha (degrees)', ...
+    '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', '\alpha (degrees)', ...
+    'YLabel', 'Secondary peak position (\theta, rad)', ...
+    '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', '\alpha (degrees)', ...
+    '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 (\theta, rad)', '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('\alpha (degrees)', '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, 'DToS_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');

Data-Analyzer/+Scripts/BECToStripes/runAngularSpectralDistributionAnalysis.m

@@ -254,7 +254,6 @@ Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'sigma2', ...
     '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, ...
@@ -339,6 +338,98 @@ Plotter.plotMeanWithSE(scan_reference_values, sigma2_vals, ...
     'SaveFileName', 'SecondaryPeakWidthMeanWithSEM.fig', ...
     'SaveDirectory', options.saveDirectory);
 
+%% ------------------ 9. 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 (\theta, rad)', '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;
+    ylim([0, 1.1])
+    xlabel('\alpha (degrees)', '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
+
+%% ------------------ 10. Save Fit Results ------------------
+outputDir = fullfile(options.saveDirectory, 'BECToS_ASDFitData');
+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');
+
 %% Inspect individual realizations of a single parameter
 
 % --- Recenter curves first ---

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

@@ -122,10 +122,31 @@ Plotter.plotSpectralCurves( ...
     'PositionThreshold', 0.5, ...                 % minimum separation between peaks
     'AmplitudeThreshold', 0.015);                 % minimum secondary peak fraction
 
+% --- Post-fit diagnostics ---
+pFit_all    = vertcat(fitResults.pFit);   % Ncurves x Nparams
+isValid_all = [fitResults.isValid];       % logical vector
+
+% Extract last three parameters (usually secondary peak amplitude, position, width)
+pTail       = pFit_all(:, end-2:end);
+
+% Count curves where *any* of these last three values are NaN or zero
+numWithNaN  = sum(any(isnan(pTail), 2));
+numWithZero = sum(any(pTail == 0, 2));
+
+% Count of valid fits
+numValid    = sum(isValid_all);
+
+% Display summary
+fprintf('\n--- Fit diagnostics (last 3 params) ---\n');
+fprintf('Curves with ≥1 NaN in last three pFit elements: %d\n', numWithNaN);
+fprintf('Curves with ≥1 zero in last three pFit elements: %d\n', numWithZero);
+fprintf('Curves marked as isValid = true: %d\n', numValid);
+fprintf('----------------------------------------\n\n');
+
 %% ------------------ Plot Fits on Raw ------------------
 options.skipLivePlot    = true;
 options.skipSaveFigures = false; 
-Plotter.plotTwoModeGaussianFitsOnRaw(fitResults, rawCurves, 4, 6, ...
+Plotter.plotTwoModeGaussianFitsOnRawRSD(fitResults, rawCurves, 4, 6, ...
     'SkipLivePlot', options.skipLivePlot, ...
     'SkipSaveFigures', options.skipSaveFigures, ...
     'SaveDirectory', options.saveDirectory);
@@ -140,7 +161,7 @@ Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'A2', ...
     'YLabel', 'Secondary peak amplitude', ...
     'FontName', options.font, ...
     'FontSize', 16, ...
-    'FigNum', 4, ...
+    'FigNum', 2, ...
     'SkipSaveFigures', options.skipSaveFigures, ...
     'SaveFileName', 'SecondaryPeakAmplitudePDF.fig', ...
     'SaveDirectory', options.saveDirectory, ...
@@ -158,7 +179,7 @@ Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'mu2', ...
     'YLabel', 'Secondary peak position (\theta, rad)', ...
     'FontName', options.font, ...
     'FontSize', 16, ...
-    'FigNum', 5, ...
+    'FigNum', 3, ...
     'SkipSaveFigures', options.skipSaveFigures, ...
     'SaveFileName', 'SecondaryPeakPositionPDF.fig', ...
     'SaveDirectory', options.saveDirectory, ...
@@ -175,7 +196,7 @@ Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'sigma2', ...
     'YLabel', 'Secondary peak width (\sigma, rad)', ...
     'FontName', options.font, ...
     'FontSize', 16, ...
-    'FigNum', 6, ...
+    'FigNum', 4, ...
     'SkipSaveFigures', options.skipSaveFigures, ...
     'SaveFileName', 'SecondaryPeakWidthPDF.fig', ...
     'SaveDirectory', options.saveDirectory, ...
@@ -193,7 +214,7 @@ Plotter.plotSecondaryGaussianRange(fitResults, scan_reference_values, ...
     'YLabel', '\mu_2 \pm \sigma_2 (rad)', ...
     'FontName', options.font, ...
     'FontSize', 16, ...
-    'FigNum', 8, ...
+    'FigNum', 5, ...
     'NumberOfBins', 20, ...
     'NormalizeHistogram', true, ...
     'Colormap', @Colormaps.coolwarm, ...
@@ -202,4 +223,159 @@ Plotter.plotSecondaryGaussianRange(fitResults, scan_reference_values, ...
     '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.15, ...                % 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 (\theta, rad)', ...
+    '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 (\theta, rad)', '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('\alpha (degrees)', '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');

Data-Analyzer/+Scripts/BECToStripesToDroplets/runAngularSpectralDistributionAnalysis.m

@@ -254,7 +254,6 @@ Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'sigma2', ...
     '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, ...
@@ -340,7 +339,6 @@ Plotter.plotMeanWithSE(scan_reference_values, sigma2_vals, ...
     'SaveDirectory', options.saveDirectory);
 
 %% ------------------ 9. 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 (\theta, rad)', 'Secondary peak width (\sigma, rad)'};
@@ -398,8 +396,8 @@ for p = 1:numel(paramNames)
     end
 end
 
-%% ------------------ 5. Save Fit Results ------------------
+%% ------------------ 10. Save Fit Results ------------------
-outputDir = fullfile(options.saveDirectory, 'SToD_FitData');
+outputDir = fullfile(options.saveDirectory, 'SToD_ASDFitData');
 if ~exist(outputDir, 'dir')
     mkdir(outputDir);
 end

Data-Analyzer/+Scripts/BECToStripesToDroplets/runRadialSpectralDistributionAnalysis.m

@@ -122,10 +122,31 @@ Plotter.plotSpectralCurves( ...
     'PositionThreshold', 0.5, ...                 % minimum separation between peaks
     'AmplitudeThreshold', 0.015);                 % minimum secondary peak fraction
 
+% --- Post-fit diagnostics ---
+pFit_all    = vertcat(fitResults.pFit);   % Ncurves x Nparams
+isValid_all = [fitResults.isValid];       % logical vector
+
+% Extract last three parameters (usually secondary peak amplitude, position, width)
+pTail       = pFit_all(:, end-2:end);
+
+% Count curves where *any* of these last three values are NaN or zero
+numWithNaN  = sum(any(isnan(pTail), 2));
+numWithZero = sum(any(pTail == 0, 2));
+
+% Count of valid fits
+numValid    = sum(isValid_all);
+
+% Display summary
+fprintf('\n--- Fit diagnostics (last 3 params) ---\n');
+fprintf('Curves with ≥1 NaN in last three pFit elements: %d\n', numWithNaN);
+fprintf('Curves with ≥1 zero in last three pFit elements: %d\n', numWithZero);
+fprintf('Curves marked as isValid = true: %d\n', numValid);
+fprintf('----------------------------------------\n\n');
+
 %% ------------------ Plot Fits on Raw ------------------
 options.skipLivePlot    = true;
 options.skipSaveFigures = false; 
-Plotter.plotTwoModeGaussianFitsOnRaw(fitResults, rawCurves, 4, 6, ...
+Plotter.plotTwoModeGaussianFitsOnRawRSD(fitResults, rawCurves, 4, 6, ...
     'SkipLivePlot', options.skipLivePlot, ...
     'SkipSaveFigures', options.skipSaveFigures, ...
     'SaveDirectory', options.saveDirectory);
@@ -140,7 +161,7 @@ Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'A2', ...
     'YLabel', 'Secondary peak amplitude', ...
     'FontName', options.font, ...
     'FontSize', 16, ...
-    'FigNum', 4, ...
+    'FigNum', 2, ...
     'SkipSaveFigures', options.skipSaveFigures, ...
     'SaveFileName', 'SecondaryPeakAmplitudePDF.fig', ...
     'SaveDirectory', options.saveDirectory, ...
@@ -158,7 +179,7 @@ Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'mu2', ...
     'YLabel', 'Secondary peak position (\theta, rad)', ...
     'FontName', options.font, ...
     'FontSize', 16, ...
-    'FigNum', 5, ...
+    'FigNum', 3, ...
     'SkipSaveFigures', options.skipSaveFigures, ...
     'SaveFileName', 'SecondaryPeakPositionPDF.fig', ...
     'SaveDirectory', options.saveDirectory, ...
@@ -175,7 +196,7 @@ Plotter.plotFitParameterPDF(fitResults, scan_reference_values, 'sigma2', ...
     'YLabel', 'Secondary peak width (\sigma, rad)', ...
     'FontName', options.font, ...
     'FontSize', 16, ...
-    'FigNum', 6, ...
+    'FigNum', 4, ...
     'SkipSaveFigures', options.skipSaveFigures, ...
     'SaveFileName', 'SecondaryPeakWidthPDF.fig', ...
     'SaveDirectory', options.saveDirectory, ...
@@ -193,7 +214,7 @@ Plotter.plotSecondaryGaussianRange(fitResults, scan_reference_values, ...
     'YLabel', '\mu_2 \pm \sigma_2 (rad)', ...
     'FontName', options.font, ...
     'FontSize', 16, ...
-    'FigNum', 8, ...
+    'FigNum', 5, ...
     'NumberOfBins', 20, ...
     'NormalizeHistogram', true, ...
     'Colormap', @Colormaps.coolwarm, ...
@@ -202,4 +223,159 @@ Plotter.plotSecondaryGaussianRange(fitResults, scan_reference_values, ...
     '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.15, ...                % 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', '\alpha (degrees)', ...
+    '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', '\alpha (degrees)', ...
+    'YLabel', 'Secondary peak position (\theta, rad)', ...
+    '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', '\alpha (degrees)', ...
+    '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 (\theta, rad)', '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('\alpha (degrees)', '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, 'SToD_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');