function a = rbf_matrix ( m, nd, xd, r0, phi )

%*****************************************************************************80
%
%% RBF_MATRIX computes the radial basis interpolation matrix.
%
%  Discussion:
%
%    Given a set of data FD = F(XD), the RBF weights can be found by
%      WD = A \ FD
%    after which the RBF interpolant can be evaluated at any point XI by
%      FI(XI) = WD' * PHI(XD-XI)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2015
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    William Press, Brian Flannery, Saul Teukolsky, William Vetterling,
%    Numerical Recipes in FORTRAN: The Art of Scientific Computing,
%    Third Edition,
%    Cambridge University Press, 2007,
%    ISBN13: 978-0-521-88068-8,
%    LC: QA297.N866.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer ND, the number of data points.
%
%    Input, real XD(M,ND), the data points.
%
%    Input, real R0, a scale factor.  R0 should be larger than the typical
%    separation between points, but smaller than the maximum separation.
%    The value of R0 has a significant effect on the resulting interpolant.
%
%    Input, function V = PHI ( R, R0 ), a function handle to evaluate the radial
%    basis functions.
%
%    Input, real FD(ND,1), the function values at the data points.
%
%    Output, real A(ND,ND), the matrix.
%
  a = zeros ( nd, nd );
 
  for i = 1 : nd
    if ( m == 1 )
      d = xd - xd(i);
      r = abs ( d );
    else
      d = xd - repmat ( xd(:,i), 1, nd );
      r = sqrt ( sum ( d.^2 ) );
    end
    a(i,1:nd) = phi ( r(1:nd), r0 );
  end

  return
end
