seval3.f

C******************** START FILE SEVAL3.FOR ; GROUP TRKRLIB ************
C
      REAL FUNCTION SEVAL3(N, U, X, Y, B, C, D, deriv, deriv2)
      INTEGER N
      REAL  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    SEVAL3 = 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  and the 2nd derivative
C
C    deriv2 = 2*C(i) + 6*D(i)*(U-X(i))
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 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
        seval3 = y(i) + dx*(b(i) + dx*(c(i) + dx*d(i)))
        deriv = b(i) + dx*(2.0*c(i) + dx*3.0*d(i))
        deriv2 = 2.0*c(i) + dx*6.0*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
        seval3=y(i)+dx*zslop
        deriv=zslop
        deriv2=0.0
      ENDIF
C
      RETURN
      END
C******************** END FILE SEVAL3.FOR ; GROUP TRKRLIB **************