function [distance, angle] = calculateGeodist(pPosition1, pPosition2, pUnit)
%CALCULATEGEODIST Calculate the distance and bearing between the two positions.
% Calculates the distance and bearing between pos1 and pos2. The optional
% parameter 'pUnit' specifies if the positions are in degrees ('degree') or
% in radiant ('rad'). If unit is not given, 'degree' would be the default.
%
% Goedist returns distance from pPosition1 to pPosition2 (dist) and the bearing
% from pPosition1 to pPosition2 (ang) in degrees rel. to north. 
%
%   Parameters:
%       pPosition1 (double) 1D-Vector with two values, latitude and
%                           longitude.
%       pPosition2 (double) 1D-Vector with two values, latitude and
%                           longitude.
%       pUnit (char)        Specifies if positions are in degrees or
%                           radiant.
%
%   Returns:
%       distance (double)   The calculated distance. [km]
%       angle (double)      The calculated bearing in degrees rel. to
%                           north.

%% Parameter check
if nargin < 3
  pUnit='degree';
end

%% Convert input
rlat = [pPosition1(1), pPosition2(1)];
rlon = [pPosition1(2), pPosition2(2)];

if strcmp(pUnit,'degree')
                    % rad = degree/180*pi !
  rlat = rlat / 180 * pi;
  rlon = rlon / 180 * pi;
end
  
rlat_size = size(rlat);
if rlat_size(1) > 1
                 % make sure that rlat and rlon have more rows than columns
   if rlat_size(1) < rlat_size(2)    
      rlat = rlat';
      rlon = rlon';
   end
end

st = [rlat(:,1) rlon(:,1)];
sc = [rlat(:,2) rlon(:,2)];

[distance, angle] = artoa.data.calculateEllipk2(st, sc);

end