!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Program file name: lagrange2.f90                                       !
!                                                                         !
!  © Tao Pang 2006                                                        !
!                                                                         !
!  Last modified: June 6, 2006                                            !
!                                                                         !
!  (1) This F90 program is created for the book, "An Introduction to      !
!      Computational Physics, 2nd Edition," written by Tao Pang and       !
!      published by Cambridge University Press on January 19, 2006.       !
!                                                                         !
!  (2) No warranties, express or implied, are made for this program.      !
!                                                                         !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
PROGRAM LAGRANGE2
!
! An example of completing the Lagrange interpolation
! with the upward and downward correction method.
!
  IMPLICIT NONE
  INTEGER, PARAMETER :: N=5
  REAL :: X, F, DF
  REAL, DIMENSION (N) :: XI, FI
  DATA XI/0.0, 0.5, 1.0, 1.5, 2.0/, &
       FI/1.0, 0.938470, 0.765198, 0.511828, 0.223891/
!
  X = 0.9
  CALL UPDOWN (N, XI, FI, X, F, DF)
  WRITE (6,"(3F16.8)") X, F, DF
END PROGRAM LAGRANGE2
!
SUBROUTINE UPDOWN (N, XI, FI, X, F, DF)
!
! Subroutine to complete the interpolation via upward and
! downward corrections.
!
  IMPLICIT NONE
  INTEGER, PARAMETER :: NMAX=21
  INTEGER, INTENT (IN) :: N
  INTEGER :: I, J, I0, J0, IT, K
  REAL, INTENT (IN) :: X
  REAL, INTENT (OUT) :: F, DF
  REAL :: DX, DXT, DT
  REAL, INTENT (IN), DIMENSION (N) :: XI, FI
  REAL, DIMENSION (NMAX,NMAX) :: DP, DM
!
  IF (N.GT.NMAX) STOP 'Dimension of the data set is too large.'
!
! Find the closest point to x
!
  DX = ABS(XI(N)-XI(1))
  DO  I = 1, N
    DP(I,I) = FI(I)
    DM(I,I) = FI(I)
    DXT = ABS(X-XI(I))
    IF (DXT.LT.DX) THEN
      I0 = I
      DX = DXT
    END IF
  END DO
  J0 = I0
!
! Evaluate the rest of the corrections recursively
!
  DO I = 1, N-1
    DO J = 1, N-I
      K = J+I
      DT =(DP(J,K-1)-DM(J+1,K))/(XI(K)-XI(J))
      DP(J,K) = DT*(XI(K)-X)
      DM(J,K) = DT*(XI(J)-X)
    END DO
  END DO
!
! Update the interpolation with the corrections
!
  F = FI(I0)
  IT = 0
  IF(X.LT.XI(I0)) IT = 1
  DO I = 1, N-1
    IF ((IT.EQ.1).OR.(J0.EQ.N)) THEN
      I0 = I0-1
      DF = DP(I0,J0)
      F  = F+DF
      IT = 0
      IF (J0.EQ.N) IT = 1
    ELSE IF ((IT.EQ.0).OR.(I0.EQ.1)) THEN
      J0 = J0+1
      DF = DM(I0,J0)
      F  = F+DF
      IT = 1
      IF (I0.EQ.1) IT = 0
    END IF
  END DO
  DF = ABS(DF)
END SUBROUTINE UPDOWN