---

polyinterp.h

Go to the documentation of this file.

#ifndef _POLYINTERP_H_
#define _POLYINTERP_H_

#include <stdio.h>



class polyinterp {

   private:
   double fourPointTable[4][4];
   double twoPointTable[4][4];

  void lagrange( double x, int N, double c[] )
  {
    double num, denom;

    for (int i = 0; i < N; i++) {
      num = 1;
      denom = 1;
      for (int k = 0; k < N; k++) {
        if (i != k) {
          num = num * (x - k);
          denom = denom * (i - k);
        }
      } // for k
      c[i] = num / denom;
    } // for i
  } // lagrange


  void fillTable( const int N, double table[4][4] )
  {
    double x;
    double n = N;
    int i = 0;

    for (x = 0.5; x < n; x = x + 1.0) {
      lagrange( x, N, table[i] );
      i++;
    }
  } // fillTable


  void printTable( double table[4][4], int N)
  {
    printf("%d-point interpolation table:\n", N);
    double x = 0.5;
    for (int i = 0; i < N; i++) {
      printf("%4.2f: ", x);
      for (int j = 0; j < N; j++) {
        printf("%6.4f", table[i][j] );
        if (j < N-1)
          printf(", ");
      }
      printf("\n");
      x = x + 1.0;
    }    
  }


  void printTables()
  {
    printTable( fourPointTable, 4 );
    printTable( twoPointTable, 2 );
    printf("\n");
  } // printTables


  void getCoef( double x, int n, double c[])
  {
    int j = (int)x;

    if (j < 0 || j >= n) {
      printf("poly::getCoef: n = %d, bad x value = %f\n", n, x);
    }

    if (n == 2) {
      c[2] = 0.0;
      c[3] = 0.0;
    }
    else if (n != 4) {
      printf("poly::getCoef: bad value for N\n");
      return;
    }

    for (int i = 0; i < n; i++) {
      if (n == 4) {
        c[i] = fourPointTable[j][i];
      }
      else { // n == 2
        c[i] = twoPointTable[j][i];
      }
    }

  } // getCoef


   public:

   polyinterp() 
   {
    // Fill in the 4-point polynomial interplation table
    // for the points 0.5, 1.5, 2.5, 3.5
    fillTable( 4, fourPointTable );
    
    // Fill in the 2-point polynomial interpolation table
    // for 0.5 and 1.5
    fillTable( 2, twoPointTable );
   } // polyinterp class constructor


  double interpPoint( double x, int N, double d[])
  {
    double c[4];
    double point = 0;
    
    int n = 4;
    if (N < 4)
      n = N;

    getCoef( x, n, c );

    if (n == 4) {
      point = c[0]*d[0] + c[1]*d[1] + c[2]*d[2] + c[3]*d[3]; 
    }
    else if (n == 2) {
      point = c[0]*d[0] + c[1]*d[1];
    }

    return point;
  } // interpPoint

};  // class interpolation


#endif

---