///////////////////////////////////////////////////////////////////////////
//                                                                       //
// Program file name: Spline.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 creating cubic-spline approximation of
// a discrete function fi=f(xi). 

import java.lang.*;
import java.io.*;
import java.util.*;
public class Spline {
  final static int n = 20;
  final static int m = 100;
  public static void main(String argv[]) throws
    IOException {
    double xi[] = new double[n+1];
    double fi[] = new double[n+1];
    double p2[] = new double[n+1];

 // Read in data points xi and and data fi
   BufferedReader r = new BufferedReader
      (new InputStreamReader
        (new FileInputStream("xy.data")));
    for (int i=0; i<=n; ++i) {
      StringTokenizer
        s = new StringTokenizer(r.readLine());
      xi[i] = Double.parseDouble(s.nextToken());
      fi[i] = Double.parseDouble(s.nextToken());
    }
    p2 = cubicSpline(xi, fi);

 // Find the approximation of the function
    double h = (xi[n]-xi[0])/m;
    double x = xi[0];
    for (int i=1; i<m; ++i) {
      x += h;

   // Find the interval that x resides
      int k = 0;
      double dx = x-xi[0];
      while (dx >= 0) { 
        ++k;
        dx = x-xi[k];
      }
      --k;
    
   // Find the value of function f(x)
      dx = xi[k+1]-xi[k];
      double alpha = p2[k+1]/(6*dx);
      double beta = -p2[k]/(6*dx);
      double gamma = fi[k+1]/dx-dx*p2[k+1]/6;
      double eta = dx*p2[k]/6 - fi[k]/dx;
      double f = alpha*(x-xi[k])*(x-xi[k])*(x-xi[k])
                +beta*(x-xi[k+1])*(x-xi[k+1])*(x-xi[k+1])
                +gamma*(x-xi[k])+eta*(x-xi[k+1]);
      System.out.println(x + " " + f);
    }
  }

// Method to carry out the cubic-spline approximation
// with the second-order derivatives returned.

  public static double[] cubicSpline(double x[],
    double f[]) {
    int n = x.length-1;
    double p[] = new double [n+1];
    double d[] = new double [n-1];
    double b[] = new double [n-1];
    double c[] = new double [n-1];
    double g[] = new double [n];
    double h[] = new double [n];

 // Assign the intervals and function differences
    for (int i=0; i<n; ++i) {
      h[i] = x[i+1]-x[i];
      g[i] = f[i+1]-f[i];
    }

 // Evaluate the coefficient matrix elements
    for (int i=0; i<n-1; ++i) {
      d[i] = 2*(h[i+1]+h[i]);
      b[i] = 6*(g[i+1]/h[i+1]-g[i]/h[i]);
      c[i] = h[i+1];
    }

 // Obtain the second-order derivatives
    g = tridiagonalLinearEq(d, c, c, b);
    for (int i=1; i<n; ++i) p[i] = g[i-1];
    return p;
  }

// Method to solve the tridiagonal linear equation set.

  public static double[] tridiagonalLinearEq(double d[],
    double e[], double c[], double b[]) {
    int m = b.length;
    double w[] = new double[m];
    double y[] = new double[m];
    double z[] = new double[m];
    double v[] = new double[m-1];
    double t[] = new double[m-1];

 // Evaluate the elements in the LU decomposition
    w[0] = d[0];
    v[0]  = c[0];
    t[0]  = e[0]/w[0];
    for (int i=1; i<m-1; ++i) {
      w[i]  = d[i]-v[i-1]*t[i-1];
      v[i]  = c[i];
      t[i]  = e[i]/w[i];
    }
    w[m-1]  = d[m-1]-v[m-2]*t[m-2];

 // Forward substitution to obtain y
    y[0] = b[0]/w[0];
    for (int i=1; i<m; ++i)
      y[i] = (b[i]-v[i-1]*y[i-1])/w[i];

 // Backward substitution to obtain z
    z[m-1] = y[m-1];
    for (int i=m-2; i>=0; --i)
      z[i] = y[i]-t[i]*z[i+1];
     
    return z;
  }
}