#include "../iolib/message.h"

float *farray1();

float cubinterpe(yy,n,dx,x,ifnewy)
float *yy;	/* data samples, uniformly distributed on x */
int n;		/* number of data samples */
float dx;       /* sampling interval */
float x;        /* on which function value is returned */
int ifnewy;	/* spline will be calculated iff ifnewy!=0 */
{
static int dim=0;
static float *xx,*y2;
float yp1,ypn,y;
int i,nx;
	
	if(n<4) STOP(cubinterpe: n too small);
	if(x>(n-1)*dx) STOP(cubinterpe: x too large);
	
	if(ifnewy)
	{ if(dim) {free(&xx[1]); free(&y2[1]);}
	  xx=farray1(1,n);
	  y2=farray1(1,n);
	  dim=n;
	  for(i=1;i<=n;i++) xx[i]=dx*(i-1);
	  yp1=(-11.*yy[0]  +18.*yy[1]  -9.*yy[2]  +2.*yy[3]  )/(6.*dx);
	  ypn=( 11.*yy[n-1]-18.*yy[n-2]+9.*yy[n-3]-2.*yy[n-4])/(6.*dx);
	  spline(xx,yy-1,n,yp1,ypn,y2);	/* yy-1: C--FORTRAN */
	}
	
	nx=x/dx;
	if(x/dx-nx==0.0) return(yy[nx]);
	
	splint(xx,yy-1,y2,n,x,&y);
	return(y);
}