/* interp.cpp
   Test interpolation routines. Enter number of points, and their
   values  xa's and ya's. Calculate the interpolation at a desired
   x. Remember: "avoid extrapolation. Former stock-market analysts
                 know it well ".
   C.A. Bertulani (04/19/2000)
*/

typedef double Number;
#include <math.h>
#include <iomanip.h>
#include <fstream.h>
#include "butil.h"

void main (){
	
	void polint(Number [], Number [], int, Number , Number *, Number *);
	void splint(Number [], Number [], Number [], int , Number , Number *);
	void spline(Number [], Number [], int, Number , Number , Number []);
	
	int n;
	Number *xa, *ya, a, b, x, *y2;
    
	/* constants for little random number generator  */
	int FM = 7875;
	int FA = 211;
	int FC = 1663;
	long int ifake=0L;

	n = 10;  // number of points
	xa=vector(n);
	ya=vector(n);
	y2=vector(n);
    ofstream out("interp.dat");              /* save data in interp.dat */

    cout << "\n Polynomial interpolation: ready to go?\n";

// generate points from an exponential funcion
	cout << "Generate table of xa vs. ya from ya=exp(-xa)" << endl;
	xa[1]=1.;
	ya[1]=10*exp(-xa[1]);
    for (int i = 2; i < n; i++) {
		/*  generate a random point xa  */
        ifake = (ifake*FA+FC) % FM;
		xa[i]= xa[i-1]+xa[i-1]*ifake/FM;  
        ya[i] = 10*exp(-xa[i]);
		out << xa[i] << "    " << ya[i] << endl;
    }

	cout  << endl << "Enter value of x" << endl;
	cin >> x;
		
		// call interpolating routine
	polint(xa, ya, n, x, &a, &b);

    out << " x = " << x << "     y = " << a << endl;
	out << " Compare with analytical result: y = " << 10*exp(-x) << endl << endl;

	// generate points from a  polynomial
	cout << "Now we try xa vs. ya from ya = 1+xa**2+xa**5.2" << endl;
	ya[1]=1.+pow(xa[1],2)+pow(xa[1],5.2);
    for (i = 2; i < n; i++) {
		ya[i]=1.+pow(xa[i],2)+pow(xa[i],5.2);
		out << xa[i] << "    " << ya[i] << endl;
    }

		
		// call interpolating routine
	polint(xa, ya, n, x, &a, &b);

    out << " x = " << x << "     y = " << a  << endl;
	out << " Compare with analytical result: y = " << 1+x*x+pow(x,5.2) << endl << endl;

    cout << "\nNow we repeat calculations with splint - spline interpolation\n";
// generate points from an exponential funcion
	ya[1]=10*exp(-xa[1]);
    for (i = 2; i < n; i++) {
        ya[i] = 10*exp(-xa[i]);
		out << xa[i] << "    " << ya[i] << endl;
    }


	// call interpolating routine
	spline(xa, ya, n, 1.e31, 1.e31, y2);
	splint(xa, ya, y2, n, x, &a);

    out << " x = " << x << "     y = " << a <<  endl;
	out << " Compare with analytical result: y = " << 10*exp(-x) << endl << endl;

		

	cout << endl << "data stored in interp.dat\n";

}

/***********************************************************************
Given arrays xa[1..n] and ya[1..n], and given a value of x, this
routine returns a value y, and an error estimate dy. If P(x) is the 
polynomial of degree N-1 such that P(xa_i) = ya_i, i=1,..n, then the
returned value y = P(x).
************************************************************************/

void polint(Number xa[], Number ya[], int n, Number x, Number *y, Number *dy)
{
	int i,m,ns=1;
	Number den,dif,dift,ho,hp,w;
	Number *c,*d;

	dif=fabs(x-xa[1]);
	c=vector(n);
	d=vector(n);
	for (i=1;i<=n;i++) {    // Here we find the index ns of the closest table entry,
		if ( (dift=fabs(x-xa[i])) < dif) {
			ns=i;
			dif=dift;
		}
		c[i]=ya[i];			// and initialize the tableau of c' and d's.
		d[i]=ya[i];
	}
	*y=ya[ns--];			// This is the initial approximation to y
	for (m=1;m<n;m++) {		// For each column of the tableau
		for (i=1;i<=n-m;i++) {	// we loop over the current c's and d's and update
			ho=xa[i]-x;			// them.
			hp=xa[i+m]-x;
			w=c[i+1]-d[i];
			try{if ( (den=ho-hp) == 0.0){
				throw "Error in routine polint";
			}
			}
			// This error can occur only if two input xa's are (to within roundoff) identical
			catch (char* message) {
            cout << message;
			}
			den=w/den;
			d[i]=hp*den;		// Here the c's and d'c are updated
			c[i]=ho*den;
		}
		*y += (*dy=(2*ns < (n-m) ? c[ns+1] : d[ns--]));
		// After each column in the tableau is completed, we decide which correction, c or d,
		// we want to add to our accumulating value of y, i.e., which path to take through
		// the tableau - forking up or down. We do this in such a way as to take the most
		// "straight line" route through the tableau to its apex, updating ns accordingly
		// to keep track of where we are. This route keeps the partial approximation
		// centered (insofar as possible) on the target x. The last dy added is thus the 
		// error indication.
	}
	delete []d;
	delete []c;
}

/******************************************************************************
Given arrays xa[1..n] and ya[1..n], which tabulate a function (with the xa's
in order), and given the array y2a[1..n], which is the output from spline 
above, and given a value of x this routine returns a cubic-spline interpolated
value y.
******************************************************************************/
void splint(Number xa[], Number ya[], Number y2a[], int n, Number x, Number *y)
{
	int klo,khi,k;
	Number h,b,a;

	klo=1;		// We will find the right place in the table by means of bisection.
	khi=n;		// This is optimal if sequential calls to this routine are at 
	while (khi-klo > 1) {	// random values of x. If sequential calls are in order, 
		k=(khi+klo) >> 1;	// and closely spaced, one would do better to store the
		if (xa[k] > x) khi=k;	// previous values of klo and khi and test if they
		else klo=k;			// remain appropriate on the next call.
	}						// klo and khi now bracket the input value of x.
	h=xa[khi]-xa[klo];
	try{if (h == 0.0) {
		throw "Bad xa input to routine splint";
	}
	}
	// This error can occur only if two input xa's are (to within roundoff) identical
	catch (char* message) {
		cout << message;
	}
	a=(xa[khi]-x)/h;
	b=(x-xa[klo])/h;
	*y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
}

/***************************************************************************
Given arrays xa[1..n] and ya[1..n], with x1 < x2 < ... < xn, and given  
values of yp1 and ypn for the first derivative of the interpolating function
at points 1 and n, respectively, this routine returns an array y2[1..n]
that contains the second derivatives of the interpolating function at the
tabulated points xi. If yp1 an/or ypn are equal to 1 x 10^30 or larger, the
routine is signaled to set the corresponding boundary condition for a
natural spline, with zero second derivative on that bounday.
****************************************************************************/
void spline(Number x[], Number y[], int n, Number yp1, Number ypn, Number y2[])
{
	int i,k;
	Number p,qn,sig,un,*u;

	u=vector(n-1);
	if (yp1 > 0.99e30)		// The lower boundary condition is set either to be
		y2[1]=u[1]=0.0;		// "natural" or else to have a specified first 
	else {					// derivative.
		y2[1] = -0.5;
		u[1]=(3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);
	}
	for (i=2;i<=n-1;i++) {	// This is the decomposition loop of the tridiagonal
		sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);	// algorithm. y2 and u are used for
		p=sig*y2[i-1]+2.0;		// temporary storage of the decomposed factors.
		y2[i]=(sig-1.0)/p;
		u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
		u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
	}
	if (ypn > 0.99e30)		// The lower boundary condition is set either to be
		qn=un=0.0;			// "natural" or else to have a specified first
	else {					// derivative.
		qn=0.5;
		un=(3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]));
	}
	y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
	for (k=n-1;k>=1;k--)	// This is the backsubstitution loop of the tridiagonal 		
		y2[k]=y2[k]*y2[k+1]+u[k];		//algorithm
	delete []u;
}