function L = lebesgue ( x )

%*****************************************************************************80
%
%% lebesgue() computes the Lebesgue coefficient for a set of nodes on an interval.
%
%  Author:
%
%    Greg von Winckel
%
%  References:
%
%    Jean-Paul Berrut, Lloyd N. Trefethen, 
%    Barycentric Lagrange Interpolation" 
%    http://web.comlab.ox.ac.uk/oucl/work/nick.trefethen/berrut.ps.gz
%
%    Walter Gaustschi, 
%    Numerical Analysis, An Introduction,
%    (1997) page 76+
%
%  Input:
%
%    real X, the node coordinates.
%
%  Output:
%
%    real L, the Lebesgue coefficient.
%

%
%  Get number of nodes
%
  n=length(x);
%
%  use this finer mesh for interpolating between nodes
%
  N=8*n+1;

  x=x(:);

  X=repmat(x,1,n);
%
%  Compute the weights
%
  w=1./prod(X-X'+eye(n),2);
%
%  Fine mesh for interpolating 
%
  xp=linspace(min(x),max(x),N)';

  xdiff=repmat(xp.',n,1)-repmat(x,1,N);
%
%  find all the points where the difference is zero
%
  zerodex=(xdiff==0); 
%
%  See eq. 3.1 in Ref (1)
%
  lfun=prod(xdiff,1);
%
%  kill zeros
%
  xdiff(zerodex)=eps;
%
%  Compute lebesgue function
%
  lebfun=sum(abs(diag(w)*repmat(lfun,n,1)./xdiff));
  plot(xp,lebfun);
  L=max(lebfun);

  return
end

