/*************************************************************************/
/* Program file name: spline.c                                           */
/*                                                                       */
/* © Tao Pang 2006                                                       */
/*                                                                       */
/* Last modified: June 22, 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 n 20
#define m 100

/* An example of creating cubic-spline approximation of
   a discrete function fi=f(xi). */

main() {
  FILE *fp;
  int i,k;
  double x, f, dx, h;
  double alpha, beta, gamma, eta;
  double xi[n+1], fi[n+1], p2[n+1];
  void cubicSpline();

// Read in data points xi and and data fi
  fp = fopen("xy.data", "r");
  for (i=0; i<=n; ++i) {
    fscanf(fp, "%lf %lf", &x, &f);
    xi[i] = x;
    fi[i] = f;
  }
  cubicSpline(xi, fi, p2);

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

// Find the interval that x resides
    k = 0;
    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];
    alpha = p2[k+1]/(6*dx);
    beta = -p2[k]/(6*dx);
    gamma = fi[k+1]/dx-dx*p2[k+1]/6;
    eta = dx*p2[k]/6 - fi[k]/dx;
    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]);
    printf("%lf %lf\n", x, f);
  }
}

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

void cubicSpline(x, f, p)
  double x[n+1], f[n+1], p[n+1];

{
  int i;
  double d[n-1], b[n-1], c[n-1], g[n], h[n];
  void tridiagonalLinearEq();

// Assign the intervals and function differences
  for (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 (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
  tridiagonalLinearEq(d, c, c, b, g);
  p[0] = 0;
  p[n] = 0;
  for (i=1; i<n; ++i) p[i] = g[i-1];
}

/* Functione to solve the tridiagonal linear equation set. */

void tridiagonalLinearEq(d, e, c, b, z)
  double d[n-1], e[n-1], c[n-1], b[n-1], z[n-1];
{
  int i;
  double y[n-1], w[n-1], v[n-2], t[n-2];

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

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

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