/********************************************************************************
 * COD potential, bound state and dynamics code                     T.K. 2008   *
 *                                                                              *
 * ops-spline.cc                                                                *
 *                                                                              *
 * spline initialization (calculation of moments, natural cubic  spline):       *
 *                                                                              *
 *             spline_init(....)                                                *
 *                                                                              *
 * actual spline interpolation:                                                 *
 *                                                                              *
 *             spline_interpolation(....)                                       *
 *                                                                              *
 * mostly adopted from the numerical recipes.                                   *
 ********************************************************************************/
#include <iostream>
#include <math.h>

using namespace std;

//Intern functions
void spline_init(double* values_X, double* values_FUNC_X,  
		 int nr_points,  double* spl_coef);
void spline_interpolation(double* values_X, double* values_FUNC_X,  
          		  double* spl_coef, int nr_points, double X, double *vout, double *dvdout);


/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* 
 *     spline_init                                                               * 
 *                                                                               * 
 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/


void spline_init(double* values_X, double* values_FUNC_X,  
		 int nr_points,  double* spl_coef){      
  double p,qn,sig,un;
  static double *dum;
  static int    alt_nr_points = -1;
  
  //Allocation of temporary space
  if(nr_points != alt_nr_points){
    if(dum != NULL)  delete(dum);
    dum = new double[nr_points];
  }
  
  //natural cubic spline
  spl_coef[0]=dum[0]=0.0;
  
  for (int i=1;i<=nr_points-2;i++) {
    sig=(values_X[i]-values_X[i-1])/(values_X[i+1]-values_X[i-1]);
    p=sig*spl_coef[i-1]+2.0;
    spl_coef[i]=(sig-1.0)/p;
    dum[i]=(values_FUNC_X[i+1]-values_FUNC_X[i])/(values_X[i+1]-values_X[i]) - 
      (values_FUNC_X[i]-values_FUNC_X[i-1])/(values_X[i]-values_X[i-1]);
    dum[i]=(6.0*dum[i]/(values_X[i+1]-values_X[i-1])-sig*dum[i-1])/p;
  }
  
  //natural cubic spline 
  qn=un=0.0;
  
  spl_coef[nr_points-1]=(un-qn*dum[nr_points-2])/(qn*spl_coef[nr_points-2]+1.0);
  for (int k=nr_points-2;k>=0;k--)
    spl_coef[k]=spl_coef[k]*spl_coef[k+1]+dum[k];
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* 
 *     spline_interpolation                                                      * 
 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

void spline_interpolation(double* values_X, double* values_FUNC_X,  
			    double* spl_coef, int nr_points, double X, double* vout, double* dvdout){
  double y, dy;
  int klo,khi,k;
  double h,b,a;

  //check boundaries
  if(X < values_X[0]){ cerr << "X to low for spine " << X  << "\n";} //return(values_FUNC_X[0]); }
  if(X > values_X[nr_points-1] ){ cerr << "X to high for spine " << X << " " << values_X[nr_points-1] << "\n"; }//return(values_FUNC_X[nr_points-1]); }

  	
  klo=0;
  khi=nr_points-1;
  while (khi-klo > 1) {
    k=(khi+klo) >> 1;
    if (values_X[k] > X) khi=k;
    else klo=k;
  }
  h=values_X[khi]-values_X[klo];
//  if(X < values_X[0]) return(values_FUNC_X[0]);
  if (h == 0.0) cerr << "Value out of splines region!";
  a=(values_X[khi]-X)/h;
  b=(X-values_X[klo])/h;
  y=a*values_FUNC_X[klo]+b*values_FUNC_X[khi]+((a*a*a-a)*spl_coef[klo]+
                                   (b*b*b-b)*spl_coef[khi])*(h*h)/6.0;

  dy= (values_FUNC_X[khi]-values_FUNC_X[klo])/h + ( (3.*b*b - 1.)*spl_coef[khi]
                                                   +(1. - 3.*a*a)*spl_coef[klo])*h/6.  ; 
  *dvdout = dy;
  *vout = y;	
//  return(y);
}