PolynomialInterpolation.h

Go to the documentation of this file.

#ifndef _utl_PolynomialInterpolation_h_
#define _utl_PolynomialInterpolation_h_

#include <cmath>
#include <utl/ErrorLogger.h>


namespace utl {

  inline
  double
  PolynomialInterpolation(const unsigned int n, const double px[], const double py[],
                          const double x, double& dy)
  {
    if (!n)
      return 0;

    double minDist = std::fabs(x - px[0]);
    unsigned int k = 0;
    double c[n];
    double d[n];
    c[0] = py[0];
    d[0] = py[0];
    double oldx = px[0];
    for (unsigned int i = 1; i < n; ++i) {
      const double newx = px[i];
      if (oldx == newx) {
        ERROR("x values should be distinct.");
        return 0;
      }
      oldx = newx;
      const double dist = std::fabs(x - newx);
      if (dist < minDist) {
        minDist = dist;
        k = i;
      }
      c[i] = d[i] = py[i];
    }

    double y = py[k];

    for (unsigned int m = 1; m < n; ++m) {
      const unsigned int nm = n - m;
      for (unsigned int i = 0; i < nm; ++i) {
        const double xa = px[i];
        const double xb = px[i+m];
        const double dx = xb - xa;
        const double cd = (c[i+1] - d[i]) / dx;
        c[i] = (x-xa)*cd;
        d[i] = (x-xb)*cd;
      }
      dy = (2*k < nm) ? c[k] : d[--k];
      y += dy;
    }

    return y;
  }


  inline
  double
  PolynomialInterpolation(const unsigned int n, const double px[], const double py[],
                          const double x)
  {
    double dummy;
    return PolynomialInterpolation(n, px, py, x, dummy);
  }

}

#endif