/*************************************************************************/
/* Program file name: lagrange2.c                                        */
/*                                                                       */
/* © Tao Pang 2006                                                       */
/*                                                                       */
/* Last modified: June 2, 2006                                           */
/*                                                                       */
/* (1) This C program is created for the book, "An Introduction to       */
/*     Computational Physics, 2nd Edition," written by Tao Pang and      */
/*     published by Cambridge University Press on January 19, 2006.      */
/*                                                                       */
/* (2) No warranties, express or implied, are made for this program.     */
/*                                                                       */
/*************************************************************************/

#include <stdio.h>
#include <math.h>
#define nmax 21

/* An example of completing the Lagrange interpolation
   with the upward and downward correction method. */

main() {
  int n = 5;
  double xi[] = {0,0.5,1,1.5,2};
  double fi[] = {1,0.938470,0.765198,0.511828,0.223891};
  double x = 0.9, f, df;
  void updown();

  updown(n, xi, fi, x, &f, &df);

  printf("%16.8lf %16.8lf %16.8lf\n", x,f,df);
}

/* Function to complete the interpolation via upward and
   downward corrections. */

void updown(n, xi, fi, x, f, df)
int n;
double x, *f, *df;
double xi[], fi[];
{
  int i, i0, j, j0, k, it;
  double dx0, dx1, dt;
  double dp[nmax][nmax], dm[nmax][nmax];

  if (n > nmax) {
    printf("Dimension too large\n");
    exit(1);
  }

/* Find the closest point to x */

  i0  = 0;
  dx0 = fabs(xi[n-1]-xi[0]);
  for (i=0; i<n; ++i) {
    dp[i][i] = fi[i];
    dm[i][i] = fi[i];
    dx1 = fabs(x-xi[i]);
    if (dx1 < dx0) {
      i0  = i;
      dx0 = dx1;
    }
  }
  j0 = i0;

/* Evaluate the rest of the corrections recursively */

  for (i=0; i<n-1; ++i) {
    for (j=0; j<n-i-1; ++j) {
      k = j+i+1;
      dt = (dp[j][k-1]-dm[j+1][k])/(xi[k]-xi[j]);
      dp[j][k] = dt*(xi[k]-x);
      dm[j][k] = dt*(xi[j]-x);
    }
  }

/* Update the interpolation with the corrections */

  *f = fi[i0];
  it = 0;
  if (x < xi[i0]) it = 1;
  for (i=0; i<n-1; ++i) {
    if ((it==1) || (j0==n-1)) {
      i0  = i0-1;
      *df = dp[i0][j0];
      *f  = *f+*df;
      it = 0;
      if (j0==n-1) it = 1;
    }
    else if ((it==0) || (i0==0)) {
      j0 = j0+1;
      *df = dm[i0][j0];
      *f  = *f+*df;
      it = 1;
      if (i0==0) it = 0;
    }
  }
  *df = fabs(*df);
}