function [D,bearing]=latlon2dist(lat1,lon1,lat2,lon2)
% latlon2dist = finds distance & bearing between pairs of points from lat/lon
%
% [D,bearing]=latlon2dist(lat1,lon1,lat2,lon2)
%
% INPUT
%  lat1    : vector of latitudes for start points (decimal degrees)
%  lon1    : vector of longitudes for start points (decimal degrees)
%  lat2    : vector of latitudes for end points (decimal degrees)
%  lon2    : vector of longitudes for end points (decimal degrees)
% 
% OUTPUT
%  D       : distance from start to end points (km)
%  bearing : bearing of end points from start points (degrees)
%
% IMB : 23/7/2009

R = 6371; % radius of earth (km)

lat1r = lat1*pi/180;
lon1r = lon1*pi/180;
lat2r = lat2*pi/180;
lon2r = lon2*pi/180;
dlonr = (lon2-lon1)*pi/180;

% distance between point 1 and point 2 (km) determined using the spherical 
% law of cosines. With 64-bit floating point numbers should be good for
% separations down to a few metres
D = R*acos( sin(lat1r).*sin(lat2r) + cos(lat1r).*cos(lat2r).*cos(lon2r-lon1r));

% bearing of point 2 from point 1 (degrees)
bearing = atan2( sin(dlonr).*cos(lat2r),...
                 cos(lat1r).*sin(lat2r) - sin(lat1r).*cos(lat2r).*cos(dlonr));
bearing = (bearing*180/pi);
ii = bearing < 0;
bearing(ii)=bearing(ii)+360;