function [plotHandles] = plotOffsetVariations(pAxesHandle, pVariationResults, 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)';

%% 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)
        continue;
    end
    
    if mod(i, 10) == 0
        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), ...
                scatter(currentResults.observations(:, 2), currentResults.observations(:, 1), 10, 'k', '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