function D = distances(lon,lat)
% computes the matrix of distances D between 
% the points in lon and lat vectors 
% the longitude and latitude of the points. 
%
 R = 6373; % radius of the Earth (km)
 n = length(lon);
 D = zeros(n,n);
 for k = 1:n
   for j = 1:k-1
     D(k,j) = R*acos(sind(lon(k))*sind(lon(j)) + ...
        cosd(lon(k))*cosd(lon(j))*cosd(lat(k)-lat(j)));
     D(j,k) = D(k,j);
   end
 end
end