function D = distances ( lon, lat )

%*****************************************************************************80
%
%% distances() computes a distance table from longitude and latitude data.
%
%  Discussion:
%
%    The distances are computed on the surface of a sphere representing
%    the earth.
%
%  Licensing:
%
%    Copyright 1984-2018 The MathWorks, Inc.
%
%  Modified:
%
%    13 March 2023
%
%  Author:
%
%    Cleve Moler
%    Modifications by John Burkardt
%
%  Reference:
%
%    Great circle distance matrix between pairs of cities.
%    https://en.wikipedia.org/wiki/Great-circle_distance.
%
%    Cleve Moler,
%    USA Traveling Salesman Tour,
%    Posted 17 September 2018,
%    logs.mathworks.com/cleve/2018/09/17/usa-traveling-salesman-tour/?s_tid=srchtitle_traveling%20salesman_1
%
%  Usage:
%
%    D = distances ( lon, lat )
%
%  Input:
%
%   real lon(*), lat(*): the longitude and latitude of several cities,
%   in degrees.
%
%  Output:
%
%    real D(*,*):  D(k,j) is the distance in miles between cities k and j.
%
%  Local:
%
%    real R: the radius of the earth in miles.
%
  R = 3959;

  n = length ( lon );
  D = zeros ( n, n );

  for k = 1 : n
    for j = 1 : k - 1
      D(k,j) = R * acos ( sind ( lat(k) ) * sind ( lat(j) ) + ...
        cosd ( lat(k) ) * cosd ( lat(j) ) * cosd ( lon(k) - lon(j) ) );
      D(j,k) = D(k,j);
    end
  end

  return
end
