function [plotHandles] = plotOffsetVariations(pAxesHandle, pVariationResults, pOffsetVariations, pColor)
%PLOTOFFSETVARIATIONS Summary of this function goes here
%   Detailed explanation goes here

%% Parameter check
if nargin < 3
    pColor = [0 0 0];
end

%% Initialize return variables
plotHandles = {};

%% Prepare variables for plotting
% get plotting positions
lat = cellfun(@(x) x.observations(end, 1), pVariationResults)';
lon = cellfun(@(x) x.observations(end, 2), pVariationResults)';
dates = cellfun(@(x) x.date, pVariationResults)';

%% Find indices of 0 1, 0 -1 and 1 0, -1 0
indicesFirstDimension = [find(all(pOffsetVariations == [0, 1], 2)), find(all(pOffsetVariations == [0, -1], 2))];
indicesSecondDimension = [find(all(pOffsetVariations == [1, 0], 2)), find(all(pOffsetVariations == [-1, 0], 2))];

%% Plot the eigenvectors and ellipse for every position
hold(pAxesHandle, 'on');
for i = 1:length(dates)
    if isnan(dates(i))
        continue;
    end
    currentResults = pVariationResults{i};
    % find ellipse by fit_ellipse function
    [~, plot_values] = fit_ellipse(currentResults.observations(:, 2), currentResults.observations(:, 1));
    if isempty(plot_values) | all(structfun(@(x) isempty(x), plot_values))
        continue;
    end
    
    [ver_line_lat, ver_line_lon] = track2(currentResults.observations(indicesFirstDimension(1), 1), currentResults.observations(indicesFirstDimension(1), 2), currentResults.observations(indicesFirstDimension(2), 1), currentResults.observations(indicesFirstDimension(2), 2));
    [hor_line_lat, hor_line_lon] = track2(currentResults.observations(indicesSecondDimension(1), 1), currentResults.observations(indicesSecondDimension(1), 2), currentResults.observations(indicesSecondDimension(2), 1), currentResults.observations(indicesSecondDimension(2), 2));
    
    if mod(i, 10) == 0
        plotHandles = [ ...
            plotHandles, ...
            { ...
                line(ver_line_lon, ver_line_lat), ...line(plot_values.new_horz_line(1, :), plot_values.new_horz_line(2, :)), ...
                line(hor_line_lon, hor_line_lat), ...
                plot(plot_values.rotated_ellipse(1, :), plot_values.rotated_ellipse(2, :), 'Color', pColor), ...
                scatter(currentResults.observations(:, 2), currentResults.observations(:, 1), 10, [.6 .6 .6], 'filled') ...
            } ...
            ];
    else
        plotHandles = [ ...
            plotHandles, ...
            { ...
                line(plot_values.new_ver_line(1, :), plot_values.new_ver_line(2, :)), ...
                line(plot_values.new_horz_line(1, :), plot_values.new_horz_line(2, :)), ...
                plot(plot_values.rotated_ellipse(1, :), plot_values.rotated_ellipse(2, :), 'Color', pColor) ...
            } ...
            ];
    end
    
end
hold(pAxesHandle, 'off');

end