/*************************************************************************/
/* Program file name: Ex2_2.java                                         */
/*                                                                       */
/* © Tao Pang 2006                                                       */
/*                                                                       */
/* Last modified: June 9, 2006                                           */
/*                                                                       */
/* (1) This Java 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.     */
/*                                                                       */
/*************************************************************************/

// Program for Exercise 2.2.  An direct implementation of the Lagrange
// interpolation scheme on the data given.

import java.lang.*;
public class Ex2_2 {
  public static void main(String argv[]) {
    double xi[] = {0, 0.4, 0.8, 1.2, 1.6};
    double fi[] = {0, 0.428392, 0.742101, 0.910314, 0.970348};

    double x = 0.3;
    double f = lagrange(x, xi, fi);
    double df = Math.abs(f-0.328627);
    System.out.println("f(0.3) = " + f + "+-" + df);

    x = 0.5;
    f = lagrange(x, xi, fi);
    df = Math.abs(f-0.520500);
    System.out.println("f(0.5) = " + f + "+-" + df);
  }

// Method to carry out the Lagrange interpolation directly.

  public static double lagrange(double x, double xi[],
    double fi[]) {
    int n = xi.length-1;
    double p[] = new double[n+1];
    double f = 0;

    for (int j=0; j<=n; ++j) {
      p[j] = 1;
      for (int k=0; k<j; ++k)
        p[j] *=(x-xi[k])/(xi[j]-xi[k]);
      for (int k=j+1; k<=n; ++k)
        p[j] *=(x-xi[k])/(xi[j]-xi[k]);
      f += fi[j]*p[j];
    }
    return f;
  }
}