!     Last change:  DOS  29 Jul 2000   12:22 pm
!                                       *** copyright 2000 ***
! *** filename splint.f95               *** John F. Monahan **
!                                       **********************
      program splint
!  splint -- demonstration of cubic spline interpolation
!
      implicit none
      integer, parameter :: n = 7
      integer, parameter :: npts = 100
      real, dimension(n) :: xk,z,mn,md,h
      real xi,s,fn,fd
      real splev
      integer i
!
  21  format(4f12.6)
  24  format(f9.6,2x,4f10.6)
!
      open( unit=6, file='splint.out')
      open( unit=9, file='splint.dat')
!                   get n points of curve
      do i = 1,n
!                   x's equally spaced on (-4,4)
      xi = 4.*real(2*i-1-n)/real(n-1)
      xk(i) = xi 
!                   function interpolated is simple-looking
      z(i) = 1./(1.+xi*xi)
      end do  ! loop on i
!                   find natural spline first
      call splstn(xk,z,mn,h,n)
      WRITE(6,21) mn
!                   now find spline with derivative conditions
      md(1) = -8./289. 
      md(n) =  8./289.  
!
      call splstd(xk,z,md,h,n)
      WRITE(6,21) md
!
!                  write out results for 2 interpolants
      do i = 1,npts
      xi = 4.*real(2*i-1-npts)/real(npts-1)
      s = 1./(1.+xi*xi)
      fn = splev(xi,xk,z,mn,h,n)
      fd = splev(xi,xk,z,md,h,n)
      write(9,24) xi,s,fn,fd
      end do  ! loop on i
      stop
      end program splint
      SUBROUTINE SPLSTN(XK,Y,M,H,N)
!  NATURAL CUBIC SPLINES WITH UNEQUALLY SPACED KNOTS
!  M(J) GIVES SECOND DERIV AT XKNOT(J)  J=1,...,N
!  J F MONAHAN    AUGUST 1981  DEPT OF STAT, NCSU, RALEIGH, NC 27650
!         RECODED OCTOBER 1999, APRIL 2000 FOR FORTRAN 95
!
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: N
      REAL, DIMENSION( N ), INTENT(IN) :: XK,Y
      REAL, DIMENSION( N ), INTENT(OUT) :: M,H
      REAL, DIMENSION(3*N) :: W
      REAL S
      INTEGER J,NM1
!
      NM1 = N-1
!                   COMPUTE DIFFERENCES H(J) AND FIRST DIFS OF Y'S
      DO J = 2,N
      H(J) = XK(J)-XK(J-1)
      M(J) = (Y(J)-Y(J-1))/H(J)         ! S(J)
      END DO  ! LOOP ON J
!                   FILL MATRIX TRIDIAGONAL MATRIX A IN VECTOR W
!      DIAGONAL W(J) = A(J,J) = 2
!   SUBDIAGONAL W(N+J) = A(J,J-1) = LAMBDA(J)
! SUPERDIAGONAL W(2*N+J) = A(J,J+1) = MU(J)
!
!  M(J) NOW HOLDS FIRST DIFFERENCES
!
!                   FIX TWO ENDS
      DO J = 2,NM1
      S = H(J)+H(J+1)
      W(2*N+J) = H(J+1)/S               ! SUPERDIAGONAL
      W(N+J) = H(J)/S                   ! SUBDIAGONAL
      W(J) = 2.                         ! DIAGONAL
      M(J) = 6.*(M(J+1)-M(J))/S         ! SECOND DIFFERENCES
      END DO  ! LOOP ON J
!                   FINISH ENDS
      W(1) = 2.
      W(N) = 2.
      W(2*N+1) = 0.
      W(2*N) = 0.
      M(1) = 0.
      M(N) = 0.
!                   SOLVE TRIDIAGONAL SYSTEM
      CALL STRIDS(W,M,N)                                                
      RETURN
      END SUBROUTINE SPLSTN
      SUBROUTINE SPLSTD(XK,Y,M,H,N)
!  CUBIC SPLINES WITH UNEQUALLY SPACED KNOTS WITH DERIVATIVE CONDITION
!  DERIVATIVES AT XK(1) & XK(N) ARE GIVEN IN M(1) & M(N) AS INPUT
!  M(J) GIVES SECOND DERIV AT XKNOT(J)  J=1,...,N
!  J F MONAHAN    AUGUST 1981  DEPT OF STAT, NCSU, RALEIGH, NC 27650
!         RECODED OCTOBER 1999, APRIL 2000 FOR FORTRAN 95
!
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: N
      REAL, DIMENSION( N ), INTENT(IN) :: XK,Y
      REAL, DIMENSION( N ), INTENT(OUT) :: M,H
      REAL, DIMENSION(3*N) :: W
      REAL S
      INTEGER J,NM1
!
      NM1 = N-1
!                   COMPUTE DIFFERENCES H(J) AND FIRST DIFS OF Y'S
      W(N) = M(N)                       ! HOLD THIS DERIV FOR LATER
      DO J = 2,N
      H(J) = XK(J)-XK(J-1)
      M(J) = (Y(J)-Y(J-1))/H(J)         ! S(J)
      END DO  ! LOOP ON J
!                   FILL MATRIX TRIDIAGONAL MATRIX A IN VECTOR W
!      DIAGONAL W(J) = A(J,J) = 2
!   SUBDIAGONAL W(N+J) = A(J,J-1) = LAMBDA(J)
! SUPERDIAGONAL W(2*N+J) = A(J,J+1) = MU(J)
!
!  M(J) NOW HOLDS FIRST DIFFERENCES
!
!                   FIX TWO ENDS
      M(1) = 6.*(M(2)-M(1))/H(2)        ! SECOND DIFFERENCE
      DO J = 2,NM1
      S = H(J)+H(J+1)
      W(2*N+J) = H(J+1)/S               ! SUPERDIAGONAL
      W(N+J) = H(J)/S                   ! SUBDIAGONAL
      W(J) = 2.                         ! DIAGONAL
      M(J) = 6.*(M(J+1)-M(J))/S         ! SECOND DIFFERENCES
      END DO  ! LOOP ON J
      M(N) = 6.*(W(N)-M(N))/H(N)        ! SECOND DIFFERENCE ** YES W(N)
!                   FINISH ENDS
      W(1) = 2.
      W(N) = 2.
      W(2*N+1) = 1.
      W(2*N) = 1.
!                   SOLVE TRIDIAGONAL SYSTEM
      CALL STRIDS(W,M,N)                                                
      RETURN                                                            
      END SUBROUTINE SPLSTD
      SUBROUTINE STRIDS(W,X,N)
!  SOLVES TRIDIAGONAL SYSTEM:  A*(X-OUT)=(X-IN)
!  WHERE DIAGONAL ELEMENTS IN W(I) = A(I,I)
!    SUBDIAGONAL ELEMENTS IN W(N+I) = A(I,I-1)
!   SUPERDIAGONAL ELEMENTS IN W(2*N+I) = A(I,I+1)
!
!  J F MONAHAN    AUGUST 1981  DEPT OF STAT, NCSU, RALEIGH, NC 27650
!       RECODED  OCTOBER 1999, APRIL 2000 FOR FORTRAN 95
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: N
      REAL, DIMENSION( 3*N ), INTENT(IN OUT) :: W
      REAL, DIMENSION( N ), INTENT(IN OUT) :: X
      INTEGER I,K
!                   COMPUTE LU FACTORIZATION (L IS UNIT LOWER)
      DO K = 2,N
      W(N+K) = W(N+K)/W(K-1)
      W(K) = W(K) - W(N+K)*W(2*N+K-1)
      END DO  ! LOOP ON K
!                   SOLVE LOWER TRIANGULAR (BIDIAGONAL) SYSTEM
      DO I = 2,N
      X(I) = X(I) - W(N+I)*X(I-1)
      END DO  ! LOOP ON I
!                   SOLVE UPPER TRIANGULAR SYSTEM (START AT BOTTOM)
      X(N) = X(N) / W(N)
      DO K = 2,N
      I = N+1-K
      X(I) = ( X(I) - W(2*N+I)*X(I+1) ) / W(I)
      END DO  ! LOOP ON K
      RETURN
      END SUBROUTINE STRIDS
      REAL FUNCTION SPLEV(X,XK,Y,M,H,N)
!  EVALUATES SPLINE AT X    XK(JC-1) LT X LE XK(JC)
!  CUBIC SPLINES WITH UNEQUALLY SPACED KNOTS
!  M(J) GIVES SECOND DERIV AT XKNOT(J)  J=1,...,N
!  J F MONAHAN    AUGUST 1981  DEPT OF STAT, NCSU, RALEIGH, NC 27650
!      RECODED   OCTOBER 1999, APRIL 2000 FOR FORTRAN 95
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: N
      REAL, INTENT(IN) :: X
      REAL, DIMENSION( N ), INTENT(IN) :: XK,Y,M,H
      REAL DA,DB
      INTEGER, SAVE :: J = 0
      INTEGER JFIND
!                   INITIAL VALUE -- RETURN THIS IF OUT OF RANGE
      SPLEV = 0.
!                   IF X OUT OF RANGE, SPLEV=0 AND QUIT
      IF( (X .LT. XK(1)) .OR. (X .GT. XK(N)) ) RETURN
!                   IF LAST CALL IN SAME INTERVAL, NO NEED TO SEARCH
      IF( ( J .LE. 1 ) .OR. ( J .GT. N ) ) THEN
             J = JFIND(X,XK,N)          ! FIND THE INTERVAL
                                           ELSE
             IF( (X .LE. XK(J-1)) .OR. (X .GT. XK(J)) ) THEN
             J = JFIND(X,XK,N)          ! FIND THE INTERVAL
                                           ENDIF
                                           ENDIF
!                   COMPUTE DIFFERENCES FROM NEARBY KNOTS
      DB = X-XK(J-1)
      DA = XK(J)-X
!                   EVALUATE CUBIC SPLINE
      SPLEV = ( (M(J-1)*DA**3 + M(J)*DB**3)/6.            &
     &      + (Y(J-1)-M(J-1)*H(J)*H(J)/6.)*DA           &
     &      + (Y(J)-M(J)*H(J)*H(J)/6.)*DB ) / H(J)
      RETURN
      END FUNCTION SPLEV
      INTEGER FUNCTION JFIND(X,XK,N)
!  FIND J SUCH THAT XK(J-1) LT X LE XK(J)
!  JFIND=2 IF X EQ XK(1)         JFIND=0 IF X IS OUT OF RANGE
!  J F MONAHAN    AUGUST 1981  DEPT OF STAT, NCSU, RALEIGH, NC 27650
!   RECODED      OCTOBER 1999, APRIL 2000 FOR FORTRAN 95
!
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: N
      REAL, INTENT(IN) :: X
      REAL, DIMENSION( N ), INTENT(IN) :: XK
      INTEGER JU,JL,I
      JFIND = 0
!                   IF X OUT OF RANGE, JFIND=0 AND QUIT
      IF( (X .LT. XK(1)) .OR. (X .GT. XK(N)) ) RETURN
!                   CHECK SPECIAL CASE
      IF( X .EQ. XK(1) ) THEN
                              JFIND = 2
                              RETURN
                         END IF
!                   INITIALIZE POINTERS
      JU = N
      JL = 1
      DO WHILE( JU-JL .GT. 1 )
!                   TEST THIS ONE
      I = ( JU + JL ) / 2
      IF( X .LE. XK(I) ) THEN
                              JU = I
                         ELSE
                              JL = I
                         END IF
      END DO  ! WHILE( JU-JL .GT. 1 )
      JFIND = JU
      RETURN
      END FUNCTION JFIND
! *** end of filename splint.f95        *********************