//-------------------------------------------------------------------//
//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.

//          If the last argument SKIPM is specified, the monotonicity
//          check will be skipped.

//          Original Author: Michael F. Saucier 10-16-87
//-------------------------------------------------------------------//

lagrange = function (tab, x0, N, skipm)
{
  local (dx, i, j, jj, jmin, lden, lnum, ...
         m, n, seq, sig, tab2, y0);

  if (!exist (N)) 
  { 
    error("Wrong number of input arguments."); 
  }

  if (!exist (skipm))
  {
    dx = diff (tab[;1]);
    sig = sign (dx[1]);
    if (any (sign (dx) - sig))
    {
      error("First column of the table must be monotonic.");
    }
  }

  // Check range of 1st column versus x0

  if (tab[1;1] > x0 || tab[tab.nr;1] < x0) 
  {
    error ("lagrange: x0 must be within range of tab");
  }

  i = find (tab[;1] == x0);
  if (i != []) 
  {
    return tab[i;2];
  }

  m = tab.nr; n = tab.nc;
  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 in jj)
  {
    for (j in jj)
    {
      if (j != i) 
      {
        lden[i] = lden[i] * (tab2[i;1] - tab2[j;1]);
      }
    }
  }
  
  y0 = sum (lnum' ./ lden' .* tab2[jj;2]);
  return y0;
};