************************************************************************
************************************************************************
************************************************************************
      SUBROUTINE POLEX(NV,XEST,YEST,NCALL,NUSE,  YOUT,DY)
C                      ^^^^^^^^^^^^^^^^^^^^^^^   ^^^^^^^
C                                 in               out


C     extrapolates functions to X=0.

C     NV functions are evaluated at X = 0 by polynomial extrapolation
C     to a progressively smaller sequence of estimates XEST and the
C     function vectors YEST.  This call is number NCALL.  Use at most
C     NUSE estimates.
C
C     Extrapolated values are output as YOUT and their estimated error
C     as DY.
C
C     Based on ideas from NUMERICAL RECIPES.

C     Code written by Phillip Helbig, Hamburg Observatory, Sep 1995.


************************************************************************
************************************************************************
C     declarations


************************************************************************
C     named constants as array bounds

C     maximum expected values of NCALL and NV

      INTEGER IMAX, NMAX

      PARAMETER (IMAX=13, NMAX=50)


************************************************************************
C     input variables from calling routine (parameter list)


C     number of functions, number of call, maximum number of previous
C     estimates to be used

      INTEGER NV, NCALL, NUSE

C     independent and dependent variables

      DOUBLE PRECISION XEST, YEST(NV)



************************************************************************
C     output variables to calling routine (parameter list)


C     output variable and estimated uncertainty

      DOUBLE PRECISION YOUT(NV), DY(NV)



************************************************************************
C     internally used quantities


C     variables

C     DO loops, calculated loop limit

      INTEGER J, K, MM

C     independent variable, table, work variables

      DOUBLE PRECISION X(IMAX), QCOL(NMAX,IMAX), D(NMAX),
     $                 DELTA, F1, F2, Q

C     save until the next call

      SAVE QCOL, X



************************************************************************
************************************************************************
C     calculations


      X(NCALL) = XEST

      DO 100 J = 1, NV

         DY(J) = YEST(J)

         YOUT(J) = YEST(J)

100   CONTINUE

C     store first estimate in the first column

      IF(NCALL.EQ.1) THEN

         DO 200 J = 1, NV

            QCOL(J,1) = YEST(J)

200      CONTINUE

      ELSE

C     use at most NUSE previous estimates

        MM = MIN(NCALL,NUSE)

         DO 300 J = 1, NV

            D(J) = YEST(J)

300     CONTINUE

        DO 500 K = 1, (MM - 1)

           DELTA = 1.0D0/(X(NCALL - K) - XEST)

           F1 = XEST*DELTA

           F2 = X(NCALL - K)*DELTA

C          propagate tableau one diagonal more

           DO 400 J = 1, NV

              Q = QCOL(J, K)

              QCOL(J,K) = DY(J)

              DELTA = D(J) - Q

              DY(J) = F1*DELTA

              D(J) = F2*DELTA

              YOUT(J) = YOUT(J) + DY(J)

400        CONTINUE

500     CONTINUE

        DO 600 J = 1, NV

           QCOL(J,MM) = DY(J)

600     CONTINUE

      ENDIF



************************************************************************
************************************************************************
      END