
function y0 = lagrange(tab,x0,N);
%LAGRANGE  Lagrange interpolation of arbitrary order. Y = LAGRANGE(TAB,X0,N)
%          returns an N-th order interpolated value from table TAB, looking
%          up X0 in the first column of TAB.

%          NOTE: TAB's 1st column is checked for monotonicity. It is an
%          error to request a value outside the range of the first column
%          of TAB for X0.

%          Michael F. Saucier 10-16-87

if (nargin ~= 3), error('Wrong number of input arguments.'), end

dx = diff(tab(:,1));
sig = sign(dx(1));
if any(sign(dx)-sig),
  error('First column of the table must be monotonic.')
end

i = find(tab(:,1) == x0);
%if i ~= 0, y0 = tab(i,2); return, end
if (~ isempty(i)), y0 = tab(i,2); return, end %DT

[m,n] = size(tab);
jmin = min(max(min(find(tab(:,1) > x0)) - fix((N+1)/2),1),m-N);
tab2 = tab(jmin:jmin+N,:);
jj = 1:N+1;

seq = x0*ones(1,N+1) - tab2(jj,1)';
lnum = prod(seq) ./ seq;

lden = ones(1,N+1);
for i=jj,
  for j=jj,
    if j~=i, lden(i) = lden(i) * (tab2(i,1)-tab2(j,1)); end
  end
end

y0 = sum(lnum' ./ lden' .* tab2(jj,2));

