/* @(#)splint.c	17.1.1.1 (ES0-DMD) 01/25/02 17:27:00 */
/*===========================================================================
  Copyright (C) 1995 European Southern Observatory (ESO)
 
  This program is free software; you can redistribute it and/or 
  modify it under the terms of the GNU General Public License as 
  published by the Free Software Foundation; either version 2 of 
  the License, or (at your option) any later version.
 
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public 
  License along with this program; if not, write to the Free 
  Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
  MA 02139, USA.
 
  Corresponding concerning ESO-MIDAS should be addressed as follows:
	Internet e-mail: midas@eso.org
	Postal address: European Southern Observatory
			Data Management Division 
			Karl-Schwarzschild-Strasse 2
			D 85748 Garching bei Muenchen 
			GERMANY
===========================================================================*/

/* @(#)splint.c	17.1.1.1  (ESO)  01/25/02  17:27: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.
 *
 ***************************************************************/

float splint( xp, X, Y, n, istart )
float xp;	/* x-value to interpolate */
float *X;	/* x-array [0..n-1] */
float *Y;	/* y-array [0..n-1] */
int n;	
int *istart;	/* initial index */
{
    float yp1, yp2, yp;
    float xpi, xpi1, l1, l2, lp1, lp2;
    int i;

    if ( xp < X[0] || xp > X[n-1] )
	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 == 0 )
	yp1 = (Y[1] - Y[0]) / (X[1] - X[0]);
    else
	yp1 = (Y[i+1] - Y[i-1]) / (X[i+1] - X[i-1]);

    if ( i >= n - 2 )
	yp2 = (Y[n-1] - Y[n-2]) / (X[n-1] - X[n-2]);
    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((float)yp);
}

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

float splint2( xp, x, y, y2, n, kstart )
float xp;	/* x-value to interpolate */
float *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;
    float 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);
}