/* @(#)splint.c	8.3 (ESO-IPG) 01/03/95 14:09:20 */
/* @(#)splint.c	1.2 (ESO-IPG) 2/9/94 18:16:12 */
/* @(#)splint.c	5.1.1.1  (ESO)  4/5/93  17:52:13 */
/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* .COPYRIGHT   (C) 1993 European Southern Observatory     */
/* .IDENT       .prg                                       */
/* .AUTHORS     Pascal Ballester (ESO/Garching)            */
/*              Cristian Levin   (ESO/La Silla)            */
/* .KEYWORDS    Spectroscopy, Long-Slit                    */
/* .PURPOSE                                                */
/* .VERSION     1.0  Package Creation  17-MAR-1993         */
/* ------------------------------------------------------- */

/* @(#)splint.c    1.0.0.0 (ESO-La Silla) 10/08/91 12:00:00 */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.IDENT        splint.c
.MODULE       subroutines -- sprebin.exe
.LANGUAGE     C
.AUTHOR       Cristian Levin - ESO La Silla
.PURPOSE      This file contains two routines for making spline
              interpolation:
                  - using Hermite polynomials.
                  - using natural-cubic spline.
              The first one is currently used for rebinning.
.KEYWORDS     rebin, spline interpolation.
.VERSION 1.0  1-Apr-1991   Implementation
.ENVIRONMENT  UNIX
------------------------------------------------------------*/

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

/***************************************************************
 *
 * splint(): spline interpolation based on Hermite polynomials.
 *
 ***************************************************************/
#ifdef __STDC__
   double splint( double xp, double *x, float *y, int n, int *istart )
#else
   double splint( 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);
}

/***************************************************************
 *
 * splint2(): natural cubic-spline interpolation.
 *
 ***************************************************************/

#ifdef __STDC__
   double splint2( double xp, double *x, float *y, float *y2, int n, 
                  int *kstart )
#else
   double splint2( xp, x, y, y2,  n, kstart )
   double xp,*x;
   float  *y,*y2;
   int    n,*kstart;
#endif

/* double xp;	 x-value to interpolate */
/* double *x;	 x-array [1..n] */
/* float  *y;	 y-array [1..n] */
/* float  *y2;	 y2-array [1..n] (2-nd derivatives) */
/* int n;	 */
/* int *kstart;	 initial index */

{
    int klo, khi, k;
    double a, b, h, yp;

    klo = *kstart;
    khi = n;

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

    for ( k = klo; k < n && xp > x[k]; k++ )
	;

    klo = *kstart = k-1;
    khi = k;

    h = x[khi] - x[klo];
    if ( h == 0.0 )
	nrerror( "Error in spline interpolation" );

    a = (x[khi] - xp) / h;
    b = (xp - x[klo]) / h;

    yp = a*y[klo] + b*y[khi] + ((a*a*a - a)*y2[klo] + (b*b*b - b)*y2[khi])*
	 (h*h) / 6.0;

    return(yp);
}