************************************************************************
*   FUNCTION:  Estimate the numerical value of a definite integral.    *
*   AUTHORS:   Lawrence Shampine, Richard Allen, Steven Pruess  for    *
*              the text  Fundamentals of Numerical Computing           *
*   DATE:      June 25, 1994                                           *
*   PURPOSE:   Array based implementation of subroutine Adapt.         *
************************************************************************
      SUBROUTINE ADAPT(F,A,B,ABSERR,RELERR,ANSWER,ERREST,FLAG)
C
C     ADAPT estimates the definite integral of F(X) from A to B
C     using an adaptive quadrature scheme based on Gauss-Kronrod
C     (3,7) formulas.
C
      INTEGER FLAG
      DOUBLE PRECISION A,B,ABSERR,RELERR,ANSWER,ERREST
      EXTERNAL F
C
C  Input parameters:
C     F      = name of the function subprogram defining F(X).
C              This program must have the form:
C                 DOUBLE PRECISION FUNCTION F(X)
C                 DOUBLE PRECISION X
C                 F = ...
C                 RETURN
C                 END
C              The function name F must appear in an EXTERNAL
C              statement in the calling program.  
C     A, B   = end points of integration interval.
C     ABSERR = absolute error tolerance desired.
C     RELERR = relative error tolerance desired.
C  Output parameters:
C     ANSWER = computed estimate of the integral.
C     ERREST = estimate of the absolute error in ANSWER.
C     FLAG   = 0  for normal return;
C            = 1  insufficient storage in queue (120);
C            = 2  too many function evaluations (3577);
C            =-1  invalid input parameters; ABSERR .LE. 0  or
C                 RELERR too small.
C  Local variables:
      INTEGER BOTTOM,LENGTH,MAXF,MAXQ,NUMF,TOP
      DOUBLE PRECISION ALPHA,BETA,E,EL,ER,H,Q,QL,QR,QUEUE,TOL,U
      DOUBLE PRECISION ZERO,HALF,TEN
C
      PARAMETER (MAXF = 3577, MAXQ = 120)
      DIMENSION QUEUE(MAXQ,4)
      DATA ZERO/0.D0/,HALF/0.5D0/,TEN/10.D0/
************************************************************************
*  Machine dependent constant:                                         *
*                                                                      *
*  Set U to the appropriate unit roundoff level for your machine;      *
*  here we assume IEEE double precision.                               *
      DATA U/2.22D-16/
************************************************************************
C
C     Test input data.
C
      IF ((ABSERR .LE. ZERO) .OR. (RELERR .LT. TEN*U)) THEN
         FLAG = -1
         RETURN
      ENDIF
C
C     Initialization.
C
      LENGTH = 0
      TOP = 1
      BOTTOM = 1
      FLAG = 0
      NUMF = 0
C
C     Form an initial approximation ANSWER to the integral over [A,B].
C     If it is not sufficiently accurate, initialize the queue and
C     begin the main loop.
C
      CALL QUAD(F,A,B,ANSWER,ERREST,NUMF)
      IF(ABS(ERREST) .GT. MAX(ABSERR,RELERR*ABS(ANSWER)))
     &   CALL ADD(QUEUE,MAXQ,ANSWER,ERREST,A,B,LENGTH,BOTTOM)
C
C     Main loop; if queue is empty then return, else subdivide
C     top entry.
C
   20 IF (LENGTH .EQ. 0) RETURN
      CALL REMOVE(QUEUE,MAXQ,Q,E,ALPHA,BETA,LENGTH,TOP)
      H = HALF*(BETA-ALPHA)
      CALL QUAD(F,ALPHA,ALPHA+H,QL,EL,NUMF)
      CALL QUAD(F,ALPHA+H,BETA,QR,ER,NUMF)
C
C     Update ANSWER and the error estimate.
C
      ANSWER = ANSWER+((QL+QR)-Q)
      ERREST = ERREST+((EL+ER)-E)
C
C     Test for failures.
C
      IF (LENGTH .GE. MAXQ-1) THEN
         FLAG = 1
         RETURN
      ENDIF
      IF (NUMF .GE. MAXF) THEN
         FLAG = 2
         RETURN
      ENDIF
C
C     Test for convergence.
C
      TOL = MAX(ABSERR,RELERR*ABS(ANSWER))
      IF (ABS(ERREST) .LE. TOL) RETURN
C
C     Add new subintervals to queue if errors are too big.
C
      TOL = TOL*H/(B-A)
      IF (ABS(EL) .GT. TOL)
     &  CALL ADD(QUEUE,MAXQ,QL,EL,ALPHA,ALPHA+H,LENGTH,BOTTOM)
      IF (ABS(ER) .GT. TOL)
     &  CALL ADD(QUEUE,MAXQ,QR,ER,ALPHA+H,BETA,LENGTH,BOTTOM)
      GOTO 20
      END
C-----------------------------------------------------------------------
      SUBROUTINE ADD(QUEUE,MAXQ,Q,E,ALPHA,BETA,LENGTH,BOTTOM)
      INTEGER MAXQ,LENGTH,BOTTOM
      DOUBLE PRECISION QUEUE(MAXQ,4),Q,E,ALPHA,BETA
C
C     Add an entry to the end of the queue.
C
      QUEUE(BOTTOM,1) = Q
      QUEUE(BOTTOM,2) = E
      QUEUE(BOTTOM,3) = ALPHA
      QUEUE(BOTTOM,4) = BETA
      LENGTH = LENGTH+1
      IF (BOTTOM .LT. MAXQ+1) BOTTOM = BOTTOM+1
      IF (BOTTOM .EQ. MAXQ+1) BOTTOM = 1
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE REMOVE(QUEUE,MAXQ,Q,E,ALPHA,BETA,LENGTH,TOP)
      INTEGER MAXQ,LENGTH,TOP
      DOUBLE PRECISION QUEUE(MAXQ,4),Q,E,ALPHA,BETA
C
C     Delete the top entry in the queue.
C
      Q = QUEUE(TOP,1)
      E = QUEUE(TOP,2)
      ALPHA = QUEUE(TOP,3)
      BETA = QUEUE(TOP,4)
      LENGTH = LENGTH-1
      IF (TOP .LT. MAXQ+1) TOP = TOP+1
      IF (TOP .EQ. MAXQ+1) TOP = 1
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE QUAD(F,ALPHA,BETA,Q,E,NUMF)
      DOUBLE PRECISION F,ALPHA,BETA,Q,E
      DOUBLE PRECISION A(4),F1,F2,F3,H,MIDPT,QKRON,X(3)
      DOUBLE PRECISION HALF,FIVE,EIGHT,NINE
      INTEGER NUMF
      EXTERNAL F
      DATA HALF/.5D0/,FIVE/5.D0/,EIGHT/8.D0/,NINE/9.D0/
C
C     Gauss-Kronrod(3,7) quadrature over (ALPHA,BETA).
C
      DATA A/0.2684880898683334D0, 0.1046562260264672D0,
     &       0.4013974147759622D0, 0.4509165386584744D0/
      DATA X/0.7745966692414834D0, 0.9604912687080202D0,
     &       0.4342437493468026D0/
      H = HALF*(BETA-ALPHA)
      MIDPT = ALPHA+H
      F1 = F(MIDPT-H*X(1))
      F2 = F(MIDPT)
      F3 = F(MIDPT+H*X(1))
      Q = H*(FIVE*(F1+F3)+EIGHT*F2)/NINE
      QKRON = H*(A(1)*(F1+F3) + A(4)*F2 + A(2)*(F(MIDPT-H*X(2)) + 
     &   F(MIDPT+H*X(2))) + A(3)*(F(MIDPT-H*X(3)) + F(MIDPT+H*X(3))))
      E = QKRON-Q
      NUMF = NUMF+7
      RETURN
      END
C=======================================================================