function d = dist(lat1,lon1,lat2,lon2)
% function d = dist(lat1,lon1,lat2,lon2)
%
% calculate distance between any two points on a sphere
% spherical coordinate in lat/lon decimal degrees

R = 6371; % Radius of Earth in km

ax = sin(lon1*pi/180)*cos(lat1*pi/180);
ay = cos(lon1*pi/180)*cos(lat1*pi/180);
az = sin(lat1*pi/180);

bx = sin(lon2*pi/180)*cos(lat2*pi/180);
by = cos(lon2*pi/180)*cos(lat2*pi/180);
bz = sin(lat2*pi/180);

a = [ax ay az];
b = [bx by bz];

theta = acos((a*b')/(norm(a)*norm(b)));
theta*180/pi;

d = R*theta;