function plotPCAResults(pcaResults, scan_parameter_values, scan_reference_values,  varargin)
%% plotPCAResults: Plots PCA results in a style consistent with plotG2
%
% Inputs:
%   pcaResults - struct returned by computePCAfromImages
%   scan_parameter_values, scan_reference_values
%   varargin   - name-value pairs (same as plotG2 plus 'FigNumRange')

    % --- Parse name-value pairs ---
    p = inputParser;
    addParameter(p, 'FontName', 'Arial', @ischar);
    addParameter(p, 'FontSize', 14, @isnumeric);
    addParameter(p, 'Colormap', @Colormaps.coolwarm);
    addParameter(p, 'SkipSaveFigures', false, @islogical);
    addParameter(p, 'SaveDirectory', pwd, @ischar);
    addParameter(p, 'FigNumRange', [], @(x) isnumeric(x) && all(x>0));
    parse(p, varargin{:});
    opts = p.Results;

    Nx         = pcaResults.Nx;
    Ny         = pcaResults.Ny;
    coeff      = pcaResults.coeff;
    score      = pcaResults.score;
    explained  = pcaResults.explained;

    raw_scan_param_vals    = scan_parameter_values;
    unique_scan_param_vals = scan_reference_values;
    numGroups              = numel(unique_scan_param_vals);
    colors                 = lines(numGroups);

    % --- Figure numbering setup ---
    if isempty(opts.FigNumRange)
        figCount = 1;
        figNums = [];
    else
        figNums = opts.FigNumRange;
        figCount = 1;
    end

    figPos = [100 100 950 750];

    %% --- Figure 1: PC1 Image ---
    pc1_image = reshape(coeff(:,1), Nx, Ny);
    if ~isempty(figNums)
        fig = figure(figNums(figCount)); clf;
    else
        fig = figure; clf;
    end
    set(fig, 'Color', 'w', 'Position', figPos);
    imagesc(pc1_image); axis image off; colormap(opts.Colormap()); colorbar;
    title(sprintf('First Principal Component (PC1) Image - Explains %.2f%% Variance', explained(1)), ...
        'FontName', opts.FontName, 'FontSize', opts.FontSize + 2);
    if ~opts.SkipSaveFigures
        Plotter.saveFigure(fig, 'SaveFileName', 'PC1_Image.fig', 'SaveDirectory', opts.SaveDirectory);
    end
    figCount = figCount + 1;

    %% --- Figure 2: PC1 Scores Distribution Scatterplot ---
    if ~isempty(figNums)
        fig = figure(figNums(figCount)); clf;
    else
        fig = figure; clf;
    end
    set(fig, 'Color', 'w', 'Position', figPos); hold on;
    for g = 1:numGroups
        idx = raw_scan_param_vals == unique_scan_param_vals(g);
        scatter(repmat(unique_scan_param_vals(g), sum(idx),1), score(idx,1), 36, colors(g,:), 'filled');
    end
    xlabel('Control Parameter', 'FontName', opts.FontName, 'FontSize', opts.FontSize);
    ylabel('PC1 Score', 'FontName', opts.FontName, 'FontSize', opts.FontSize);
    title('Evolution of PC1 Scores', 'FontName', opts.FontName, 'FontSize', opts.FontSize + 2);
    grid on;
    set(gca,'XDir','reverse');  % Ensure decreasing scan_reference_values go left-to-right
    if ~opts.SkipSaveFigures
        Plotter.saveFigure(fig, 'SaveFileName', 'PC1_Scatter.fig', 'SaveDirectory', opts.SaveDirectory);
    end
    figCount = figCount + 1;

    %% --- Figure 3: PC1 Scores Distribution Histograms ---
    numTiles = min(6, numGroups);  % show up to 6 groups
    tileIndices = round(linspace(1, numGroups, numTiles));

    if ~isempty(figNums)
        fig = figure(figNums(figCount)); clf;
    else
        fig = figure; clf;
    end
    set(fig, 'Color', 'w', 'Position', figPos);
    tLayout = tiledlayout(3,2,'TileSpacing','compact','Padding','compact');

    for t = 1:numTiles
        g = tileIndices(t);
        idx = raw_scan_param_vals == unique_scan_param_vals(g);
        data = score(idx,1);

        nexttile;
        histogram(data, 'Normalization', 'pdf', 'FaceColor', colors(g,:), 'FaceAlpha', 0.3);
        hold on;
        [f, xi] = ksdensity(data);
        plot(xi, f, 'Color', colors(g,:), 'LineWidth', 2);
        yl = ylim;
        plot([median(data) median(data)], yl, 'k--', 'LineWidth', 1);
        xlabel('PC1 Score');
        ylabel('Probability');
        title(sprintf('Control = %g', unique_scan_param_vals(g)));
        grid on;
    end
    sgtitle('PC1 Score Distributions');
    if ~opts.SkipSaveFigures
        Plotter.saveFigure(fig, 'SaveFileName', 'PC1_Distributions.fig', 'SaveDirectory', opts.SaveDirectory);
    end
    figCount = figCount + 1;

    %% --- Figure 4: PC1 Scores Distribution Boxplot ---
    % Construct group labels explicitly
    groupLabels = arrayfun(@num2str, raw_scan_param_vals, 'UniformOutput', false);
    
    % Create categorical variable with specified order
    groupCats = categorical(groupLabels, ...
                            arrayfun(@num2str, unique_scan_param_vals, 'UniformOutput', false), ...
                            'Ordinal', true);
    
    if ~isempty(figNums)
        fig = figure(figNums(figCount)); clf;
    else
        fig = figure; clf;
    end
    set(fig, 'Color', 'w', 'Position', figPos);
    
    % Plot boxplot with categorical groups
    boxplot(score(:,1), groupCats);
    
    xlabel('Control Parameter', 'FontName', opts.FontName, 'FontSize', opts.FontSize);
    ylabel('PC1 Score', 'FontName', opts.FontName, 'FontSize', opts.FontSize);
    title('PC1 Score Boxplots by Group', 'FontName', opts.FontName, 'FontSize', opts.FontSize + 2);
    grid on;
    
    if ~opts.SkipSaveFigures
        Plotter.saveFigure(fig, 'SaveFileName', 'PC1_Boxplot.fig', 'SaveDirectory', opts.SaveDirectory);
    end
    figCount = figCount + 1;


    %% --- Figure 5: PC1 Scores Distribution Mean ± SEM ---
    if ~isempty(figNums)
        fig = figure(figNums(figCount)); clf;
    else
        fig = figure; clf;
    end
    set(fig, 'Color', 'w', 'Position', figPos);
    meanScores = arrayfun(@(g) mean(score(raw_scan_param_vals == g,1)), unique_scan_param_vals);
    semScores  = arrayfun(@(g) std(score(raw_scan_param_vals == g,1))/sqrt(sum(raw_scan_param_vals == g)), unique_scan_param_vals);
    errorbar(unique_scan_param_vals, meanScores, semScores, '-o', 'LineWidth', 2);
    xlabel('Control Parameter', 'FontName', opts.FontName, 'FontSize', opts.FontSize);
    ylabel('Mean PC1 Score ± SEM', 'FontName', opts.FontName, 'FontSize', opts.FontSize);
    title('Mean ± SEM of PC1 Scores', 'FontName', opts.FontName, 'FontSize', opts.FontSize + 2);
    grid on;
    set(gca,'XDir','reverse');  % consistent ordering
    if ~opts.SkipSaveFigures
        Plotter.saveFigure(fig, 'SaveFileName', 'PC1_MeanSEM.fig', 'SaveDirectory', opts.SaveDirectory);
    end
    figCount = figCount + 1;

    %% --- Figure 6: PC1 Scores Distribution Binder Cumulant ---
    cumulantsAll = cell2mat(arrayfun(@(g) {Calculator.computeCumulants(score(raw_scan_param_vals == g,1), 4)}, unique_scan_param_vals));
    cumulantsAll = reshape(cumulantsAll, 4, numGroups);
    binderVals   = cumulantsAll(4,:);

    if ~isempty(figNums)
        fig = figure(figNums(figCount)); clf;
    else
        fig = figure; clf;
    end
    set(fig, 'Color', 'w', 'Position', figPos);
    plot(unique_scan_param_vals, binderVals*1E-5, '-o', 'LineWidth', 2);  % scale like older code
    xlabel('Control Parameter', 'FontName', opts.FontName, 'FontSize', opts.FontSize);
    ylabel('\kappa (\times 10^5)', 'FontName', opts.FontName, 'FontSize', opts.FontSize);
    title('Binder Cumulant of PC1 Scores', 'FontName', opts.FontName, 'FontSize', opts.FontSize + 2);
    grid on;
    set(gca,'XDir','reverse');  % consistent ordering
    if ~opts.SkipSaveFigures
        Plotter.saveFigure(fig, 'SaveFileName', 'PC1_BinderCumulant.fig', 'SaveDirectory', opts.SaveDirectory);
    end

    %% --- ANOVA Test ---
    p = anova1(score(:,1), groupLabels, 'off');
    fprintf('[INFO] ANOVA p-value for PC1 score differences between groups: %.4e\n', p);

end