r8seval.f

C******************** START FILE SEVAL.FOR ; GROUP TRKRLIB *************
C
      REAL*8 FUNCTION r8seval(N, U, X, Y, B, C, D)
      IMPLICIT NONE
      INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
      INTEGER N
      real*8  u, x(n), y(n), b(n), c(n), d(n)
C
CCCCCCCCCCCCCCC
CCCCCCCCCCCCCCC
C  THIS SUBROUTINE EVALUATES THE CUBIC SPLINE FUNCTION
C
C    SEVAL = Y(I) + B(I)*(U-X(I)) + C(I)*(U-X(I))**2 + D(I)*(U-X(I))**3
C
C    WHERE  X(I) .LT. U .LT. X(I+1), USING HORNER'S RULE
C
C  IF  U .LT. X(1) THEN  I = 1  IS USED.
C  IF  U .GE. X(N) THEN  I = N  IS USED.
C
C  INPUT..
C
C    N = THE NUMBER OF DATA POINTS
C    U = THE ABSCISSA AT WHICH THE SPLINE IS TO BE EVALUATED
C    X,Y = THE ARRAYS OF DATA ABSCISSAS AND ORDINATES
C    B,C,D = ARRAYS OF SPLINE COEFFICIENTS COMPUTED BY SPLINE
C
C  IF  U  IS NOT IN THE SAME INTERVAL AS THE PREVIOUS CALL, THEN A
C  BINARY SEARCH IS PERFORMED TO DETERMINE THE PROPER INTERVAL.
C
      INTEGER I, J, K
      real*8 dx
      DATA i/1/
      IF ( i .GE. n ) i = 1
      IF ( u .LT. x(i) ) GO TO 10
      IF ( u .LE. x(i+1) ) GO TO 30
C
C  BINARY SEARCH
C
   10 i = 1
      j = n+1
   20 k = (i+j)/2
      IF ( u .LT. x(k) ) j = k
      IF ( u .GE. x(k) ) i = k
      IF ( j .GT. i+1 ) GO TO 20
C
C  EVALUATE SPLINE
C
   30 dx = u - x(i)
      r8seval = y(i) + dx*(b(i) + dx*(c(i) + dx*d(i)))
C
      RETURN
      END
C******************** END FILE SEVAL.FOR ; GROUP TRKRLIB ***************