!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Program file name: lagrange.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 LAGRANGE
!
! An example of extracting an approximate function
! value via the Lagrange interpolation scheme.
!
  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 AITKEN (N, XI, FI, X, F, DF)
  WRITE (6,"(3F16.8)") X, F, DF
END PROGRAM LAGRANGE
!
SUBROUTINE AITKEN (N, XI, FI, X, F, DF)
!
! Subroutine to carry out the Aitken recursions.
!
  INTEGER, PARAMETER :: NMAX=21
  INTEGER, INTENT (IN) :: N
  INTEGER :: I,J
  REAL, INTENT (IN) :: X
  REAL, INTENT (OUT) :: F, DF
  REAL :: X1, X2, F1, F2
  REAL, INTENT (IN), DIMENSION (N):: XI, FI
  REAL, DIMENSION (NMAX):: FT
!
  IF (N.GT.NMAX) STOP 'Dimension of the data is too large.'
  DO I = 1, N
    FT(I) = FI(I)
  END DO
!
  DO I = 1, N-1  
    DO J = 1, N-I
      X1 = XI(J)
      X2 = XI(J+I)
      F1 = FT(J)
      F2 = FT(J+1)
      FT(J) = (X-X1)/(X2-X1)*F2+(X-X2)/(X1-X2)*F1
    END DO
  END DO
  F = FT(1) 
  DF = (ABS(F-F1)+ABS(F-F2))/2.0
END SUBROUTINE AITKEN