ccccccccccccccccccccccccc     Program 2.A     ccccccccccccccccccccccccccc
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                       c
c Please Note:                                                          c
c                                                                       c
c (1) This computer program is written by Tao Pang in conjunction with  c
c     his book, "An Introduction to Computational Physics," published   c
c     by Cambridge University Press in 1997.                            c
c                                                                       c
c (2) No warranties, express or implied, are made for this program.     c
c                                                                       c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      PROGRAM INTERPOLATION2
C
C Main program for the Lagrange interpolation with the
C upward and downward correction method.
C Copyright (c) Tao Pang 1997.
C
      PARAMETER (N=5)
      DIMENSION XI(N),FI(N)
      DATA XI/0.0,0.5,1.0,1.5,2.0/,
     *     FI/1.0,0.938470,0.765198,0.511828,0.223891/
C
      X = 0.9
      CALL UPDOWN (N,XI,FI,X,F,DF)
      WRITE (6, 999) X,F,DF
      STOP
  999 FORMAT (3F16.8)
      END
C
      SUBROUTINE UPDOWN (N,XI,FI,X,F,DF)
C
C Subroutine performing the Lagrange interpolation with the
C upward and downward correction method.  F: interpolated
C value.  DF: error estimated.
C Copyright (c) Tao Pang 1997.
C
      PARAMETER (NMAX=21) 
      DIMENSION XI(N),FI(N),DP(NMAX,NMAX),DM(NMAX,NMAX)
C
      IF (N.GT.NMAX) STOP 'Dimension too large'
      DX = ABS(XI(N)-XI(1))
      DO     100 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
        ELSE
        END IF
  100 CONTINUE
      J0 = I0
C
C Evaluate correction matrices
C
      DO     200 I = 1, N-1
        DO   150 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)
  150   CONTINUE
  200 CONTINUE
C
C Update the approximation
C
      F = FI(I0)
      IT = 0
      IF(X.LT.XI(I0)) IT = 1
      DO     300 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
  300 CONTINUE
      DF = ABS(DF)
      RETURN
      END