function [results, initialResults, variationResults] = calculateByOffsetVariations(pPreparedData, pCurrentCombination, pFloatReferenceTime, pSoundVelocity, pTrackingParameter, pOffsetVariations)
%CALCULATEERRORELLIPSES Summary of this function goes here
%   Detailed explanation goes here

%% Create all combinations of offsets that should be calculated
offsetVariations = [ ...
    pOffsetVariations; ...
    0, 0 ...
    ];

variations = cell(4, 1);
variationResults = cell(4, 1);

soundsources = strsplit(pCurrentCombination.soundsources{1}, ' ');
for i = 1:length(offsetVariations)
    variations{i} = pPreparedData;
    % apply variations to the data
    for k = 1:length(soundsources)
        variations{i}.(soundsources{k}).toa = variations{i}.(soundsources{k}).toa + offsetVariations(i, k);
    end
    variationResults{i} = struct();
    [variationResults{i}.positions, variationResults{i}.dates] = helperCalculateCombinationSegment(variations{i});
    % create date vector that contains all common dates
    if (i == 1)
        intersectedDates = variationResults{i}.dates;
    else
        intersectedDates = intersect(intersectedDates, variationResults{i}.dates);
    end

end

%% Prepare return variables
results = cell(1, length(intersectedDates));
initialResults = variationResults(end);

%% Calculate covariance for every date of the initial results
for i = 1:length(variationResults{end}.dates)
    results{i} = struct();
    observations = NaN(length(variationResults), 2);
    for k = 1:length(variationResults)
        selection = variationResults{k}.dates == variationResults{end}.dates(i);
        if sum(selection) ~= 1
            continue;
        end
        observations(k, :) = variationResults{k}.positions(selection, :);
    end
    results{i}.date = variationResults{end}.dates(i);
    results{i}.observations = observations;
    % get covariance matrix
    results{i}.covarianceMatrix = cov(observations);
    % eigendecomposition
    [results{i}.V, results{i}.D] = eig(results{i}.covarianceMatrix);
end

variationResults = variationResults(1:end - 1);

%% Helper functions
    function [segmentPositions, segmentDates] = helperCalculateCombinationSegment(pCustomizedToaData)
        [ ...
            segmentPositions, ...
            segmentDates, ...
            ~, ...
            ~ ...
            ] = artoa.trajectory.calculateCombinationSegment( ...
            pCustomizedToaData, ...
            pCurrentCombination, ...
            pFloatReferenceTime, ...
            pSoundVelocity, ...
            pCurrentCombination.backwardTracking, ...
            false ...
        );

        % remove NaNs from segment
        validPositions = all(~isnan(segmentPositions), 2);
        segmentPositions = segmentPositions(validPositions, :);
        segmentDates = segmentDates(validPositions);
    end

end