function D=dist_gc(x1lat,x1lon,x2lat,x2lon)

% Works for x1 scalar and x2 vector, both scalars, or both vectors of the same length

% Great circle distances (Haversine formula)
A=sind((x1lat-x2lat)/2).^2+cosd(x2lat).*cosd(x1lat).*sind((x1lon-x2lon)/2).^2;
% First factor is volumetric mean earth radius in km
D=6371.0*2*atan2(sqrt(A),real(sqrt(1-A)));

