/* @(#)dsplint.c	17.1.1.1 (ESO-IPG) 01/25/02 17:51:58 */
/*

  Christian Levin, Thomas Szeifert
  
  splint.c

  spline interpolation based on Hermite polynomials

 */

/* system includes */

#include <stdio.h>
#include <math.h>

#ifdef __STDC__
   double dsplint( double xp, double *x, float *y, int n, int *istart )
#else
   double dsplint( xp, x, y, n, istart )
   double xp,*x;
   float *y;
   int n,*istart;
#endif

/* double xp;	 x-value to interpolate */
/* double *x;	 x-array [1..n] */
/* float  *y;	 y-array [1..n] */
/* int *istart; initial index   */

{
double yp1, yp2, yp;
double xpi, xpi1, l1, l2, lp1, lp2;
int i;

if (xp < x[1] || xp > x[n]) return(0.0);

for (i=*istart; i<=n && xp>=x[i]; i++)
   ;

*istart = i;
i--;
    
lp1 = 1.0 / (x[i] - x[i+1]);
lp2 = -lp1;

if ( i == 1 )
   yp1 = (y[2] - y[1]) / (x[2] - x[1]);
else
   yp1 = (y[i+1] - y[i-1]) / (x[i+1] - x[i-1]);

if ( i >= n - 1 )
   yp2 = (y[n] - y[n-1]) / (x[n] - x[n-1]);
else
   yp2 = (y[i+2] - y[i]) / (x[i+2] - x[i]);

xpi1 = xp - x[i+1];
xpi  = xp - x[i];
l1   = xpi1*lp1;
l2   = xpi*lp2;

yp = y[i]*(1 - 2.0*lp1*xpi)*l1*l1 + 
     y[i+1]*(1 - 2.0*lp2*xpi1)*l2*l2 + 
     yp1*xpi*l1*l1 + yp2*xpi1*l2*l2;

return(yp);
}