r8seval2.f

C******************** START FILE SEVAL2.FOR ; GROUP TRKRLIB ************
C
      REAL*8 FUNCTION r8seval2(N, U, X, Y, B, C, D, deriv)
!============
! idecl:  explicitize implicit REAL declarations:
      IMPLICIT NONE
      INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
      real*8 zslop
!============
      INTEGER N
      real*8  u, x(n), y(n), b(n), c(n), d(n), deriv
C
CCCCCCCCCCCCCCC
CCCCCCCCCCCCCCC
C  THIS SUBROUTINE EVALUATES THE CUBIC SPLINE FUNCTION
C
C    SEVAL2 = Y(I) + B(I)*(U-X(I)) + C(I)*(U-X(I))**2 + D(I)*(U-X(I))**3
C
C  and the derivative
C
C    deriv = B(i) + 2*C(i)*(U-X(i)) + 3*D(i)*(U-X(i))**2
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, ILIN
      real*8 dx
      DATA i/1/
      DATA ilin/0/
C
      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)
      IF(ilin.EQ.0) THEN
        r8seval2 = y(i) + dx*(b(i) + dx*(c(i) + dx*d(i)))
        deriv = b(i) + dx*(2.0_r8*c(i) + dx*3.0_r8*d(i))
      ELSE
        IF(i.EQ.n) THEN
          zslop=(y(n)-y(n-1))/(x(n)-x(n-1))
        ELSE
          zslop=(y(i+1)-y(i))/(x(i+1)-x(i))
        ENDIF
        r8seval2=y(i)+dx*zslop
        deriv=zslop
      ENDIF
C
      RETURN
      END
C******************** END FILE SEVAL2.FOR ; GROUP TRKRLIB **************