      SUBROUTINE QUANC8(FUN,A,B,ABSERR,RELERR,RESULT,ERREST,NOFUN,FLAG)
C
      DOUBLE PRECISION FUN, A, B, ABSERR, RELERR, RESULT, ERREST, FLAG
      INTEGER NOFUN
C
C   ESTIMATE THE INTEGRAL OF FUN(X) FROM A TO B
C   TO A USER PROVIDED TOLERANCE.
C   AN AUTOMATIC ADAPTIVE ROUTINE BASED ON
C   THE 8-PANEL NEWTON-COTES RULE.
C
C   INPUT ..
C
C   FUN     THE NAME OF THE INTEGRAND FUNCTION SUBPROGRAM FUN(X).
C   A       THE LOWER LIMIT OF INTEGRATION.
C   B       THE UPPER LIMIT OF INTEGRATION.(B MAY BE LESS THAN A.)
C   RELERR  A RELATIVE ERROR TOLERANCE. (SHOULD BE NON-NEGATIVE)
C   ABSERR  AN ABSOLUTE ERROR TOLERANCE. (SHOULD BE NON-NEGATIVE)
C
C   OUTPUT ..
C
C   RESULT  AN APPROXIMATION TO THE INTEGRAL HOPEFULLY SATISFYING THE
C           LEAST STRINGENT OF THE TWO ERROR TOLERANCES.
C   ERREST  AN ESTIMATE OF THE MAGNITUDE OF THE ACTUAL ERROR.
C   NOFUN   THE NUMBER OF FUNCTION VALUES USED IN CALCULATION OF RESULT.
C   FLAG    A RELIABILITY INDICATOR.  IF FLAG IS ZERO, THEN RESULT
C           PROBABLY SATISFIES THE ERROR TOLERANCE.  IF FLAG IS
C           XXX.YYY , THEN  XXX = THE NUMBER OF INTERVALS WHICH HAVE
C           NOT CONVERGED AND  0.YYY = THE FRACTION OF THE INTERVAL
C           LEFT TO DO WHEN THE LIMIT ON  NOFUN  WAS APPROACHED.
C
      DOUBLE PRECISION W0,W1,W2,W3,W4,AREA,X0,F0,STONE,STEP,COR11,TEMP
      DOUBLE PRECISION QPREV,QNOW,QDIFF,QLEFT,ESTERR,TOLERR
      DOUBLE PRECISION QRIGHT(31),F(16),X(16),FSAVE(8,30),XSAVE(8,30)
      DOUBLE PRECISION DABS,DMAX1
      INTEGER LEVMIN,LEVMAX,LEVOUT,NOMAX,NOFIN,LEV,NIM,I,J
C
C   ***   STAGE 1 ***   GENERAL INITIALIZATION
C   SET CONSTANTS.
C
      LEVMIN = 1
      LEVMAX = 30
      LEVOUT = 6
      NOMAX = 5000
      NOFIN = NOMAX - 8*(LEVMAX-LEVOUT+2**(LEVOUT+1))
C
C   TROUBLE WHEN NOFUN REACHES NOFIN
C
      W0 =   3956.0D0 / 14175.0D0
      W1 =  23552.0D0 / 14175.0D0
      W2 =  -3712.0D0 / 14175.0D0
      W3 =  41984.0D0 / 14175.0D0
      W4 = -18160.0D0 / 14175.0D0
C
C   INITIALIZE RUNNING SUMS TO ZERO.
C
      FLAG = 0.0D0
      RESULT = 0.0D0
      COR11  = 0.0D0
      ERREST = 0.0D0
      AREA   = 0.0D0
      NOFUN = 0
      IF (A .EQ. B) RETURN
C
C   ***   STAGE 2 ***   INITIALIZATION FOR FIRST INTERVAL
C
      LEV = 0
      NIM = 1
      X0 = A
      X(16) = B
      QPREV  = 0.0D0
      F0 = FUN(X0)
      STONE = (B - A) / 16.0D0
      X(8)  =  (X0  + X(16)) / 2.0D0
      X(4)  =  (X0  + X(8))  / 2.0D0
      X(12) =  (X(8)  + X(16)) / 2.0D0
      X(2)  =  (X0  + X(4))  / 2.0D0
      X(6)  =  (X(4)  + X(8))  / 2.0D0
      X(10) =  (X(8)  + X(12)) / 2.0D0
      X(14) =  (X(12) + X(16)) / 2.0D0
      DO 25 J = 2, 16, 2
         F(J) = FUN(X(J))
   25 CONTINUE
      NOFUN = 9
C
C   ***   STAGE 3 ***   CENTRAL CALCULATION
C   REQUIRES QPREV,X0,X2,X4,...,X16,F0,F2,F4,...,F16.
C   CALCULATES X1,X3,...X15, F1,F3,...F15,QLEFT,QRIGHT,QNOW,QDIFF,AREA.
C
   30 X(1) = (X0 + X(2)) / 2.0D0
      F(1) = FUN(X(1))
      DO 35 J = 3, 15, 2
         X(J) = (X(J-1) + X(J+1)) / 2.0D0
         F(J) = FUN(X(J))
   35 CONTINUE
      NOFUN = NOFUN + 8
      STEP = (X(16) - X0) / 16.0D0
      QLEFT  =  (W0*(F0 + F(8))  + W1*(F(1)+F(7))  + W2*(F(2)+F(6))
     1  + W3*(F(3)+F(5))  +  W4*F(4)) * STEP
      QRIGHT(LEV+1)=(W0*(F(8)+F(16))+W1*(F(9)+F(15))+W2*(F(10)+F(14))
     1  + W3*(F(11)+F(13)) + W4*F(12)) * STEP
      QNOW = QLEFT + QRIGHT(LEV+1)
      QDIFF = QNOW - QPREV
      AREA = AREA + QDIFF
C
C   ***   STAGE 4 *** INTERVAL CONVERGENCE TEST
C
      ESTERR = DABS(QDIFF) / 1023.0D0
      TOLERR = DMAX1(ABSERR,RELERR*DABS(AREA)) * (STEP/STONE)
      IF (LEV .LT. LEVMIN) GO TO 50
      IF (LEV .GE. LEVMAX) GO TO 62
      IF (NOFUN .GT. NOFIN) GO TO 60
      IF (ESTERR .LE. TOLERR) GO TO 70
C
C   ***   STAGE 5   ***   NO CONVERGENCE
C   LOCATE NEXT INTERVAL.
C
   50 NIM = 2*NIM
      LEV = LEV+1
C
C   STORE RIGHT HAND ELEMENTS FOR FUTURE USE.
C
      DO 52 I = 1, 8
         FSAVE(I,LEV) = F(I+8)
         XSAVE(I,LEV) = X(I+8)
   52 CONTINUE
C
C   ASSEMBLE LEFT HAND ELEMENTS FOR IMMEDIATE USE.
C
      QPREV = QLEFT
      DO 55 I = 1, 8
         J = -I
         F(2*J+18) = F(J+9)
         X(2*J+18) = X(J+9)
   55 CONTINUE
      GO TO 30
C
C   ***   STAGE 6   ***   TROUBLE SECTION
C   NUMBER OF FUNCTION VALUES IS ABOUT TO EXCEED LIMIT.
C
   60 NOFIN = 2*NOFIN
      LEVMAX = LEVOUT
      FLAG = FLAG + (B - X0) / (B - A)
      GO TO 70
C
C   CURRENT LEVEL IS LEVMAX.
C
   62 FLAG = FLAG + 1.0D0
C
C   ***   STAGE 7   ***   INTERVAL CONVERGED
C   ADD CONTRIBUTIONS INTO RUNNING SUMS.
C
   70 RESULT = RESULT + QNOW
      ERREST = ERREST + ESTERR
      COR11  = COR11  + QDIFF / 1023.0D0
C
C   LOCATE NEXT INTERVAL.
C
   72 IF (NIM .EQ. 2*(NIM/2)) GO TO 75
      NIM = NIM/2
      LEV = LEV-1
      GO TO 72
   75 NIM = NIM + 1
      IF (LEV .LE. 0) GO TO 80
C
C   ASSEMBLE ELEMENTS REQUIRED FOR THE NEXT INTERVAL.
C
      QPREV = QRIGHT(LEV)
      X0 = X(16)
      F0 = F(16)
      DO 78 I = 1, 8
         F(2*I) = FSAVE(I,LEV)
         X(2*I) = XSAVE(I,LEV)
   78 CONTINUE
      GO TO 30
C
C   ***   STAGE 8   ***   FINALIZE AND RETURN
C
   80 RESULT = RESULT + COR11
C
C   MAKE SURE ERREST NOT LESS THAN ROUNDOFF LEVEL.
C
      IF (ERREST .EQ. 0.0D0) RETURN
   82 TEMP = DABS(RESULT) + ERREST
      IF (TEMP .NE. DABS(RESULT)) RETURN
      ERREST = 2.0D0*ERREST
      GO TO 82
      END
