///////////////////////////////////////////////////////////////////////////
//                                                                       //
// Program file name: Deriv2.java                                        //
//                                                                       //
// © Tao Pang 2006                                                       //
//                                                                       //
// Last modified: January 18, 2006                                       //
//                                                                       //
// (1) This Java program is part of 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.     //
//                                                                       //
///////////////////////////////////////////////////////////////////////////

// An example of evaluating the 1st-order derivative with the
// nonuniform 3-point formula for f(x)=sin(x).

import java.lang.*;
public class Deriv2 {
  static final int n = 100, m = 5;
  public static void main(String argv[]) {
    double[] x = new double[n+1];
    double[] f = new double[n+1];
    double[] f1 = new double[n+1];
    double[] f2 = new double[n+1];

 // Assign constants, data points, and function
    int k = 2;
    double h = Math.PI/(2*n);
    for (int i=0; i<=n; ++i) {
      x[i]  = h*i;
      f[i]  = Math.sin(x[i]);
    }

 // Calculate 1st-order derivative of the data
    f1 = firstOrderDerivative2(x, f, k);

// Output the result in every m data points
    for (int i=0; i<=n; i+=m) {
      double df1 = f1[i] - Math.cos(x[i]);
      System.out.println("x = " + x[i]);
      System.out.println("f'(x) = " + f1[i]);
      System.out.println("Error in f'(x): " + df1);
    }
  }

// Method to calculate the 1st-order derivative with the
// nonuniform 3-point formula.  Extrapolations are made
// at the boundaries.

  public static double[] firstOrderDerivative2(double x[],
    double f[], int k) {
    int n = x.length-1;
    double[] y = new double[n+1];
    double[] xl = new double[k+1];
    double[] fl = new double[k+1];
    double[] xr = new double[k+1];
    double[] fr = new double[k+1];

 // Calculate the derivative at the field points
    double h0 = x[1]-x[0];
    double a0 = h0*h0;
    for (int i=1; i<n; ++i) {
      double h = x[i+1]-x[i];
      double a = h*h;
      double b = a-a0;
      double c = h*h0*(h+h0);
      y[i] = (a0*f[i+1]+b*f[i]-a*f[i-1])/c;
      h0 = h;
      a0 = a;
    }

 // Lagrange-extrapolate the boundary points
    for (int i=1; i<=(k+1); ++i) {
      xl[i-1] = x[i]-x[0];
      fl[i-1] = y[i];
      xr[i-1] = x[n]-x[n-i];
      fr[i-1] = y[n-i];
    }
    y[0] = aitken(0, xl, fl);
    y[n] = aitken(0, xr, fr);
    return y;
  }

// The Aitken method for the Lagrange interpolation.

  public static double aitken(double x, double xi[],
    double fi[]) {
    int n = xi.length-1;
    double ft[] = (double[]) fi.clone();
    for (int i=0; i<n; ++i) {
      for (int j=0; j<n-i; ++j) {
        ft[j] = (x-xi[j])/(xi[i+j+1]-xi[j])*ft[j+1]
               +(x-xi[i+j+1])/(xi[j]-xi[i+j+1])*ft[j];
      }
    }
    return ft[0];
  }
}