*-----------------------------------------------------------------------
* poly - find an interpolating polynomial by the Bjorck-Pereyra algorithm
*-----------------------------------------------------------------------
* refer to the Algol:
* Ake Bjorck and Victor Pereyra
* "Solution of Vandermonde Systems of Equations"
* Mathematics of Computation, v.24 n.112 pp.893-903, October 1970
*-----------------------------------------------------------------------
*                            I/O LIST
*___Name_____Type_______________In/Out___Description____________________
*   X(0:N)   Double Precision   In       Independent variable values
*   Y(0:N)   Double Precision   In       Dependent variable values
*   C(0:N)   Double Precision   Out      Coefficients of polynomial equation
*   N        Integer            In       Order of the polynomial
*-----------------------------------------------------------------------
       SUBROUTINE POLY (X, Y, C, N)
        IMPLICIT NONE
        INTEGER I, J, N
        DOUBLE PRECISION X, Y, C
        DIMENSION X(0:N), Y(0:N), C(0:N)

        DO J = 0, N
          C(J) = Y(J)
        END DO
        DO J = 0, N-1
          DO I = N, J+1, -1
            C(I) = (C(I) - C(I-1)) / (X(I) - X(I-J-1))
          END DO
        END DO
        DO J = N-1, 0, -1
          DO I = J, N-1
            C(I) = C(I) - X(J) * C(I+1)
          END DO
        END DO
        RETURN
       END  ! of poly