stellinstall/r8spleen_8f_source.html

C-------------------------------------------

C  derived from Forsyth/Malcolm/Moler SPLINE.FOR -- different BC:

C  y'=0 at RHS of data interval (high index side)

C

CCCCCCCCCCCCCCC

      SUBROUTINE r8spleen (N, X, Y, B, C, D)

      IMPLICIT NONE

      INTEGER, PARAMETER :: R8=selected_real_kind(12,100)

      INTEGER N

      real*8 x(n), y(n), b(n), c(n), d(n)

C

C  THE COEFFICIENTS B(I), C(I), AND D(I), I=1,2,...,N ARE COMPUTED

C  FOR A CUBIC INTERPOLATING SPLINE FOR WHICH S'(XN)=0.

C

C    S(X) = Y(I) + B(I)*(X-X(I)) + C(I)*(X-X(I))**2 + D(I)*(X-X(I))**3

C

C    FOR  X(I) .LE. X .LE. X(I+1)

C

C  INPUT..

C

C    N = THE NUMBER OF DATA POINTS OR KNOTS (N.GE.2)

C    X = THE ABSCISSAS OF THE KNOTS IN STRICTLY INCREASING ORDER

C    Y = THE ORDINATES OF THE KNOTS

C

C  OUTPUT..

C

C    B, C, D  = ARRAYS OF SPLINE COEFFICIENTS AS DEFINED ABOVE.

C

C  USING  P  TO DENOTE DIFFERENTIATION,

C

C    Y(I) = S(X(I))

C    B(I) = SP(X(I))

C    C(I) = SPP(X(I))/2

C    D(I) = SPPP(X(I))/6  (DERIVATIVE FROM THE RIGHT)

C

CCCCCCCCCCCCCCC

C  THE ACCOMPANYING FUNCTION SUBPROGRAM  SEVAL  CAN BE USED

C  TO EVALUATE THE SPLINE.

C

C

      INTEGER NM1, IB, I

      real*8 t

C

      nm1 = n-1

      IF ( n .LT. 2 ) RETURN

      IF ( n .LT. 3 ) GO TO 50

C

C  SET UP TRIDIAGONAL SYSTEM

C

C  B = DIAGONAL, D = OFFDIAGONAL, C = RIGHT HAND SIDE.

C

      d(1) = x(2) - x(1)

      c(2) = (y(2) - y(1))/d(1)

      DO 10 i = 2, nm1

         d(i) = x(i+1) - x(i)

         b(i) = 2._r8*(d(i-1) + d(i))

         c(i+1) = (y(i+1) - y(i))/d(i)

         c(i) = c(i+1) - c(i)

   10 CONTINUE

C

C  END CONDITIONS.  THIRD DERIVATIVES AT  X(1)  AND  X(N)

C  OBTAINED FROM DIVIDED DIFFERENCES

C

      b(1) = -d(1)

      b(n) = +2._r8*d(n-1)

      c(1) = 0._r8

      c(n) = 0._r8

      IF ( n .EQ. 3 ) GO TO 15

      c(1) = c(3)/(x(4)-x(2)) - c(2)/(x(3)-x(1))

      c(1) = c(1)*d(1)**2/(x(4)-x(1))

        c(n)= -(y(n)-y(n-1))/d(n-1)

C

C  FORWARD ELIMINATION

C

   15 DO 20 i = 2, n

         t = d(i-1)/b(i-1)

         b(i) = b(i) - t*d(i-1)

         c(i) = c(i) - t*c(i-1)

   20 CONTINUE

C

C  BACK SUBSTITUTION

C

      c(n) = c(n)/b(n)

      DO 30 ib = 1, nm1

         i = n-ib

         c(i) = (c(i) - d(i)*c(i+1))/b(i)

   30 CONTINUE

C

C  C(I) IS NOW THE SIGMA(I) OF THE TEXT

C

C  COMPUTE POLYNOMIAL COEFFICIENTS

C

      b(n) = (y(n) - y(nm1))/d(nm1) + d(nm1)*(c(nm1) + 2._r8*c(n))

      DO 40 i = 1, nm1

         b(i) = (y(i+1) - y(i))/d(i) - d(i)*(c(i+1) + 2._r8*c(i))

         d(i) = (c(i+1) - c(i))/d(i)

         c(i) = 3._r8*c(i)

   40 CONTINUE

      c(n) = 3._r8*c(n)

      d(n) = d(n-1)

      RETURN

C

   50 b(1) = (y(2)-y(1))/(x(2)-x(1))

      c(1) = 0._r8

      d(1) = 0._r8

      b(2) = b(1)

      c(2) = 0._r8

      d(2) = 0._r8

      RETURN

      END