function results = extractAutocorrelation(theta_values, angular_spectral_distribution, scan_parameter_values, N_shots, N_angular_bins)
    %% Extract g² (autocorrelation) from experimental images
    % Computes angular autocorrelation g² for a set of experimental images.
    % Returns both all repetition curves and their mean + error for each scan parameter.

    % ===== Convert spectral distributions to matrix =====
    delta_nkr_all = zeros(N_shots, N_angular_bins);
    for k = 1:N_shots
        delta_nkr_all(k, :) = angular_spectral_distribution{k};
    end

    % ===== Group images by scan parameter values =====
    [unique_scan_parameter_values, ~, idx] = unique(scan_parameter_values);
    N_params = length(unique_scan_parameter_values);

    % ===== Preallocate outputs =====
    g2_curves   = cell(1, N_params);   % each cell: [N_reps × N_angular_bins]
    g2_mean     = zeros(N_params, N_angular_bins);
    g2_error    = zeros(N_params, N_angular_bins);

    % ===== Compute g²(θ) for each scan parameter group =====
    for i = 1:N_params
        group_idx  = find(idx == i);                  
        group_data = delta_nkr_all(group_idx, :);      
        N_reps     = length(group_idx);

        g2_matrix = zeros(N_reps, N_angular_bins);

        for j = 1:N_reps
            profile = group_data(j, :);

            for dtheta = 0:N_angular_bins-1
                profile_shifted = circshift(profile, -dtheta, 2);

                num   = mean(profile .* profile_shifted);
                denom = mean(profile.^2);

                g2_matrix(j, dtheta+1) = num / denom;
            end
        end

        % Store per-repetition curves
        g2_curves{i} = g2_matrix;

        % Compute mean and standard error
        g2_mean(i,:)  = mean(g2_matrix, 1, 'omitnan');
        g2_error(i,:) = std(g2_matrix, 0, 1, 'omitnan') ./ sqrt(N_reps);
    end

    % ===== Package results =====
    results = struct();
    results.g2_curves             = g2_curves;                  % raw [N_reps × N_angular_bins] per group
    results.g2_mean               = g2_mean;                    % mean per group
    results.g2_error              = g2_error;                   % standard error per group
    results.theta_values          = theta_values;
    results.scan_parameter_values = unique_scan_parameter_values;
end