ABACUS/src/FITTING/polint.cc

/**********************************************************

This software is part of J.-S. Caux's ABACUS library.

Copyright (c) J.-S. Caux.

-----------------------------------------------------------

File:  polint.cc

Purpose:  Polynomial interpolation

***********************************************************/

#include "ABACUS.h"
using namespace std;

void ABACUS::polint(Vect_DP& xa, Vect_DP& ya, const DP x, DP& y, DP& dy)
{
  // Polynomial interpolation/extrapolation, NR page 113.

  int i, m, ns=0;
  DP den, dif, dift, ho, hp, w;

  int n = xa.size();
  Vect_DP c(n), d(n);
  dif = fabs(x-xa[0]);
  for (i = 0; i < n; i++) {
    if ((dift = fabs(x-xa[i])) < dif) {
      ns = i;
      dif = dift;
    }
    c[i] = ya[i];
    d[i] = ya[i];
  }
  y = ya[ns--];
  for (m = 1; m < n; m++) {
    for (i = 0; i < n-m; i++) {
      ho = xa[i] - x;
      hp = xa[i+m] - x;
      w = c[i+1] - d[i];
      if ((den = ho-hp) == 0.0) ABACUSerror("Error in routine polint.");
      den = w/den;
      d[i] = hp * den;
      c[i] = ho * den;
    }
    y += (dy = (2 * (ns + 1) < (n-m) ? c[ns + 1] : d[ns--]));
  }
}