** LAGRAN **************************************************************
*  Function to implement Lagrange interpolation. Local variables used: *
*     I, J  : subscripts                                               *
*     PROD  : a product in Lagrange's interpolation formula            *
*     SUM   : running sum of the products in Lagrange's formula        *
*                                                                      *
*  Accepts:  Arrays X and F, where F[I] = function value at X[I],      *
             the dimension LIMIT of these arrays, the degree N of the  *
*            interpolating polynomial, and an x value XSTAR at which   *
*            the function is to be approximated                        *
*  Returns:  Value of Lagrange polynomial at X                         *
************************************************************************

      FUNCTION LAGRAN(X, F, LIMIT, N, XSTAR)

      INTEGER LIMIT, N, I, J
      REAL LAGRAN, X(0:LIMIT), F(0:LIMIT), XSTAR, SUM, PROD

      SUM = 0
      DO 20 I = 0, N
         PROD = F(I)
         DO 10 J = 0, N
            IF (I .NE. J) THEN
               PROD = PROD * (XSTAR - X(J)) / (X(I) - X(J))
            END IF
10       CONTINUE
         SUM = SUM + PROD
20    CONTINUE

      LAGRAN = SUM
      END

