SUBROUTINE DQ1DA(A,B,EPS,R,E,KF,IFLAG)
C***BEGIN PROLOGUE  DQ1DA
C***DATE WRITTEN  870525
C***REVISION DATE  870525
C***CATEGORY NO.  H2A1A1
C***KEYWORDS   ADAPTIVE  QUADRATUE, AUTOMATIC  QUADRATURE
C***AUTHOR  KAHANER, DAVID K., SCIENTIFIC COMPUTING DIVISION, NBS.
C***PURPOSE  APPROXIMATES ONE DIMENSIONAL INTEGRALS OF USER DEFINED
C            FUNCTIONS, EASY TO USE. DOUBLE PRECISION VERSION OF Q1DA.
C
C***DESCRIPTION
C       DQ1DA IS A SUBROUTINE FOR THE AUTOMATIC EVALUATION
C            OF THE DEFINITE INTEGRAL OF A USER DEFINED FUNCTION
C            OF ONE VARIABLE. DOUBLE PRECISION VERSION OF Q1DA
C
C       A R G U M E N T S   I N   T H E   C A L L   S E Q U E N C E
C
C       A
C       B     (INPUT) THE ENDPOINTS OF THE INTEGRATION INTERVAL
C       EPS   (INPUT) THE ACCURACY TO WHICH YOU WANT THE INTEGRAL
C                COMPUTED.  IF YOU WANT 2 DIGITS OF ACCURACY SET
C                EPS=.01, FOR 3 DIGITS SET EPS=.001, ETC.
C                EPS MUST BE POSITIVE.
C       R     (OUTPUT) DQ1DA'S BEST ESTIMATE OF YOUR INTEGRAL
C       E     (OUTPUT) AN ESTIMATE OF ABS(INTEGRAL-R)
C       KF    (OUTPUT) THE COST OF THE INTEGRATION, MEASURED IN
C                   NUMBER OF EVALUATIONS OF YOUR INTEGRAND.
C                   KF WILL ALWAYS BE AT LEAST 30.
C       IFLAG (OUTPUT) TERMINATION FLAG...POSSIBLE VALUES ARE
C               0   NORMAL COMPLETION, E SATISFIES
C                        E<EPS  AND  E<EPS*ABS(R)
C               1   NORMAL COMPLETION, E SATISFIES
C                        E<EPS, BUT E>EPS*ABS(R)
C               2   NORMAL COMPLETION, E SATISFIES
C                        E<EPS*ABS(R), BUT E>EPS
C               3   NORMAL COMPLETION BUT EPS WAS TOO SMALL TO
C                     SATISFY ABSOLUTE OR RELATIVE ERROR REQUEST.
C
C               4   ABORTED CALCULATION BECAUSE OF SERIOUS ROUNDING
C                     ERROR.  PROBABLY E AND R ARE CONSISTENT.
C               5   ABORTED CALCULATION BECAUSE OF INSUFFICIENT STORAGE.
C                     R AND E ARE CONSISTENT.
C               6   ABORTED CALCULATION BECAUSE OF SERIOUS DIFFICULTIES
C                     MEETING YOUR ERROR REQUEST.
C               7   ABORTED CALCULATION BECAUSE EPS WAS SET <= 0.0
C
C            NOTE...IF IFLAG=3, 4, 5 OR 6 CONSIDER USING DQ1DAX INSTEAD.
C
C       NOTE... A,B,EPS, R AND E MUST BE DECLARED DOUBLE PRECSION IN
C                YOUR CALLING PROGRAM.
C    W H E R E   I S   Y O U R   I N T E G R A N D ?
C
C        YOU MUST WRITE A FORTRAN FUNCTION, CALLED F, TO EVALUATE
C        THE INTEGRAND.  USUALLY THIS LOOKS LIKE...
C                 DOUBLE PRECISION FUNCTION F(X)
C                    DOUBLE PRECISION X
C                    F=(EVALUATE THE INTEGRAND AT THE POINT X)
C                    RETURN
C                 END
C
C
C    T Y P I C A L   P R O B L E M   S E T U P
C
C          DOUBLE PRECISION A,B,EPS,R,E
C
C          A=0.0
C          B=1.0          (SET INTERVAL ENDPOINTS TO [0,1])
C          EPS=0.001       (SET ACCURACY REQUEST FOR 3 DIGITS)
C          CALL DQ1DA(A,B,EPS,R,E,KF,IFLAG)
C          END
C          DOUBLE PRECISION FUNCTION F(X)
C            DOUBLE PRECISION X
C            F=SIN(2.*X)-SQRT(X)     (FOR EXAMPLE)
C            RETURN
C          END
C      FOR THIS SAMPLE PROBLEM, THE OUTPUT IS
C  0.0    1.0     .001    .041406750    .69077D-07    30    0
C
C    R E M A R K   I.
C
C           A SMALL AMOUT OF RANDOMIZATION IS BUILT INTO THIS PROGRAM.
C           CALLING DQ1DA A FEW TIMES IN SUCCESSION WILL GIVE DIFFERENT
C           BUT HOPEFULLY CONSISTENT RESULTS.
C
C   R E M A R K   II.
C
C           THIS ROUTINE IS DESIGNED FOR INTEGRATION OVER A FINITE
C           INTERVAL.  THUS THE INPUT ARGUMENTS A AND B MUST BE
C           VALID REAL NUMBERS ON YOUR COMPUTER.  IF YOU WANT TO DO
C           AN INTEGRAL OVER AN INFINITE INTERVAL SET A OR B OR BOTH
C           LARGE ENOUGH SO THAT THE INTERVAL [A,B] CONTAINS MOST OF
C           THE INTEGRAND.  CARE IS NECESSARY, HOWEVER.  FOR EXAMPLE,
C           TO INTEGRATE EXP(-X*X) ON THE ENTIRE REAL LINE ONE COULD
C           TAKE A=-20., B=20. OR SIMILAR VALUES TO GET GOOD RESULTS.
C           IF YOU TOOK A=-1.E10 AND B=+1.E10 TWO BAD THINGS WOULD
C           OCCUR. FIRST, YOU WILL CERTAINLY GET AN ERROR MESSAGE FROM
C           THE EXP ROUTINE, AS ITS ARGUMENT IS TOO SMALL.  OTHER
C           THINGS COULD HAPPEN TOO, FOR EXAMPLE AN UNDERFLOW.
C           SECOND, EVEN IF THE ARITHMETIC WORKED PROPERLY DQ1DA WILL
C           SURELY GIVE AN INCORRECT ANSWER, BECAUSE ITS FIRST TRY
C           AT SAMPLING THE INTEGRAND IS BASED ON YOUR SCALING AND
C           IT IS VERY UNLIKELY TO SELECT EVALUATION POINTS IN THE
C           INFINITESMALLY SMALL INTERVAL [-20,20] WHERE ALL THE
C           INTEGRAND IS CONCENTRATED, WHEN A, B ARE SO LARGE.
C
C    M O R E   F L E X I B I L I T Y
C
C           DQ1DA IS AN EASY TO USE DRIVER FOR ANOTHER PROGRAM, DQ1DAX.
C           DQ1DAX PROVIDES SEVERAL OPTIONS WHICH ARE NOT AVAILABLE
C                WITH DQ1DA.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED (DQ1DAX)
C    E N D   O F   D O C U M E N T A T I O N
C***END PROLOGUE  DQ1DA
C
      DOUBLE PRECISION A,B,E,EPS,F,FMAX,FMIN,R,W(50,6)
      LOGICAL RST
      EXTERNAL F
C
C***FIRST EXECUTUTABLE STATEMENT DQ1DA
      NINT=1
      RST = .FALSE.
      NMAX=50
      CALL DQ1DAX(F,A,B,EPS,R,E,NINT,RST,W,NMAX,FMIN,FMAX,KF,IFLAG)
      RETURN
      END
      SUBROUTINE DQ1DAX(
     *      F,A,B,EPS,R,E,NINT,RST,W,NMAX,FMIN,FMAX,KF,IFLAG)
C***BEGIN PROLOGUE DQ1DAX
C***DATE WRITTEN 870525
C***REVISION DATE 870525
C***CATEGORY NO. H2A1A1
C***KEYWORDS  ADAPTIVE QUADRATURE, AUTOMATIC QUADRATURE
C***AUTHOR  KAHANER, DAVID K., SCIENTIFIC COMPUING DIVISION, NBS.
C***PURPOSE  APPROXIMATES ONE DIMENSIONAL INTEGRALS OF USER DEFINED
C            FUNCTIONS, FLEXIBLE USAGE. DOUBLE PRECISION VERSION OF
C            Q1DAX.
C
C***DESCRIPTION
C
C       DQ1DAX IS A FLEXIBLE SUBROUTINE FOR THE AUTOMATIC EVALUATION
C             OF DEFINITE INTEGRALS OF A USER DEFINED FUNCTION
C             OF ONE VARIABLE. DOUBLE PRECISION VERSION OF Q1DAX.
C
C             FOR AN EASIER TO USE ROUTINE SEE DQ1DA.
C
C             CAPABILITIES OF DQ1DAX (IN ADDITION TO THOSE OF DQ1DA)
C             INCLUDE:
C                ABILITY TO RESTART A CALCULATION TO GREATER
C                  ACCURACY WITHOUT PENALTY...
C                ABILITY TO SPECIFY AN INITIAL PARTITION OF
C                  THE INTEGRATION INTERVAL...
C                ABILITY TO INCREASE THE WORK SPACE TO HANDLE
C                  MORE DIFFICULT PROBLEMS...
C                OUTPUT OF LARGEST/SMALLEST INTEGRAND VALUE FOR
C                  APPLICATIONS SUCH AS SCALING GRAPHS...
C
C       A R G U M E N T S   I N   T H E   C A L L   S E Q U E N C E
C
C       F      (INPUT) THE NAME OF YOUR INTEGRAND FUNCTION.
C                      THIS NAME MUST APPEAR IN AN EXTERNAL STATEMENT
C                      IN ANY PROGRAM WHICH CALLS DQ1DAX.
C                      YOU MUST WRITE F IN THE FORM
C                           DOUBLE PRECISION FUNCTION F(X)
C                              DOUBLE PRECISION X
C                              F=(EVALUATE INTEGRAND AT THE POINT X)
C                              RETURN
C                           END
C       A
C       B       (INPUT) ENDPOINTS OF INTEGRATION INTERVAL
C       EPS     (INPUT)  ACCURACY TO WHICH THE INTEGRAL IS TO BE CALCULATED.
C                          DQ1DAX WILL TRY TO ACHIEVE RELATIVE ACCURACY,
C                          E.G. SET EPS=.01 FOR 2 DIGITS, .001 FOR 3, ETC.
C       R       (OUTPUT) THE ESTIMATE OF THE INTEGRAL
C       E       (OUTPUT) THE ESTIMATE OF THE ABSOLUTE ERROR IN R.
C       NINT    (INPUT
C                OUTPUT)
C                        AS AN INPUT QUANTITY, NINT MUST BE SET TO
C                          THE NUMBER OF SUBINTERVALS IN THE INITIAL
C                          PARTITION OF [A,B].  FOR MOST PROBLEMS
C                          THIS IS JUST 1, THE INTERVAL [A,B] ITSELF.
C                          NINT MUST BE LESS THAN NMAX, SEE BELOW.
C                          NINT IS USEFUL IF YOU WOULD LIKE TO HELP
C                          DQ1DAX LOCATE A DIFFICULT SPOT ON [A,B].
C                          IN THIS REGARD NINT IS USED ALONG
C                          WITH THE ARRAY W (SEE BELOW).  IF YOU SET
C                          NINT=1 IT IS NOT NECESSARY TO BE CONCERNED
C                          WITH W, EXCEPT THAT IT MUST BE DIMENSIONED...
C                          AS AN EXAMPLE OF MORE GENERAL APPLICATIONS,
C                          IF [A,B]=[0,1] BUT THE INTEGRAND JUMPS AT 0.3,
C                          IT WOULD BE WISE TO SET NINT=2 AND THEN SET
C                                  W(1,1)=0.0  (LEFT ENDPOINT)
C                                  W(2,1)=0.3  (SINGULAR POINT)
C                                  W(3,1)=1.0  (RIGHT ENDPOINT)
C                          IF YOU SET NINT GREATER THAN 1, BE SURE TO
C                          CHECK THAT YOU HAVE ALSO SET
C                                W(1,1)=A  AND  W(NINT+1,1)=B
C                        AS AN OUTPUT QUANTITY, NINT GIVES THE
C                          NUMBER OF SUBINTERVALS IN THE FINAL
C                          PARTITION OF [A,B].
C       RST     (INPUT) A LOGICAL VARIABLE (E.G. TRUE OR FALSE)
C                       SET RST=.FALSE. FOR INITIAL CALL TO DQ1DAX
C                       SET RST=.TRUE. FOR A SUBSEQUENT CALL,
C                           E.G. ONE FOR WHICH MORE ACCURACY IS
C                           DESIRED (SMALLER EPS).  A RESTART ONLY
C                           MAKES SENSE IF THE PRECEDING CALL RETURNED
C                           WITH A VALUE OF IFLAG (SEE BELOW) LESS THAN 3.
C                           ON A RESTART YOU MAY NOT CHANGE THE VALUES OF ANY
C                           OTHER ARGUMENTS IN THE CALL SEQUENCE, EXCEPT EPS.
C       W(NMAX,6) (INPUT &
C                  SCRATCH)
C                       W IS AN ARRAY USED BY DQ1DAX.
C                        YOU  M U S T  INCLUDE A DIMENSION STATEMENT IN
C                        YOUR CALLING PROGRAM TO ALLOCATE THIS STORAGE.
C                        THIS SHOULD BE OF THE FORM
C                                   DIMENSION W(NMAX,6)
C                        WHERE NMAX IS AN INTEGER. AN ADEQUATE VALUE OF
C                        NMAX IS 50.  IF YOU SET NINT>1 YOU MUST ALSO
C                        INITIALIZE W, SEE NINT ABOVE.
C       NMAX    (INPUT) AN INTEGER EQUAL TO THE FIRST SUBSCRIPT IN THE
C                        DIMENSION STATEMENT FOR THE ARRAY W.  THIS IS
C                        ALSO EQUAL TO THE MAXIMUM NUMBER OF SUBINTERVALS
C                        PERMITTED IN THE INTERNAL PARTITION OF [A,B].
C                        A VALUE OF 50 IS AMPLE FOR MOST PROBLEMS.
C       FMIN
C       FMAX    (OUTPUT) THE SMALLEST AND LARGEST VALUES OF THE INTEGRAND
C                          WHICH OCCURRED DURING THE CALCULATION.  THE
C                          ACTUAL INTEGRAND RANGE ON [A,B] MAY, OF COURSE,
C                          BE GREATER BUT PROBABLY NOT BY MORE THAN 10%.
C       KF      (OUTPUT) THE ACTUAL NUMBER OF INTEGRAND EVALUATIONS USED
C                          BY DQ1DAX TO APPROXIMATE THIS INTEGRAL.  KF
C                          WILL ALWAYS BE AT LEAST 30.
C       IFLAG   (OUTPUT) TERMINATION FLAG...POSSIBLE VALUES ARE
C                  0   NORMAL COMPLETION, E SATISFIES
C                           E<EPS  AND  E<EPS*ABS(R)
C                  1   NORMAL COMPLETION, E SATISFIES
C                           E<EPS, BUT E>EPS*ABS(R)
C                  2   NORMAL COMPLETION, E SATISFIES
C                           E<EPS*ABS(R), BUT E>EPS
C                  3   NORMAL COMPLETION BUT EPS WAS TOO SMALL TO
C                        SATISFY ABSOLUTE OR RELATIVE ERROR REQUEST.
C                  4   ABORTED CALCULATION BECAUSE OF SERIOUS ROUNDING
C                        ERROR.  PROBABLY E AND R ARE CONSISTENT.
C                  5   ABORTED CALCULATION BECAUSE OF INSUFFICIENT STORAGE.
C                        R AND E ARE CONSISTENT.  PERHAPS INCREASING NMAX
C                        WILL PRODUCE BETTER RESULTS.
C                  6   ABORTED CALCULATION BECAUSE OF SERIOUS DIFFICULTIES
C                        MEETING YOUR ERROR REQUEST.
C                  7   ABORTED CALCULATION BECAUSE EITHER EPS, NINT OR NMAX
C                        HAS BEEN SET TO AN ILLEGAL VALUE.
C                  8   ABORTED CALCULATION BECAUSE YOU SET NINT>1 BUT FORGOT
C                        TO SET W(1,1)=A  AND  W(NINT+1,1)=B
C       NOTE... A,B,EPS,R,E,W,FMIN AND FMAX MUST BE DECLARED DOUBLE
C                  PRECISION IN YOUR CALLING PROGRAM.
C
C     T Y P I C A L   P R O B L E M   S E T   U P
C
C      DOUBLE PRECISION A,B,EPS,W(50,6),R,E,FMIN,FMAX
C      LOGICAL RST
C      EXTERNAL F
C      A=0.0
C      B=1.0
C      W(1,1)=A
C      W(2,1)=.3      [SET INTERNAL PARTITION POINT AT .3]
C      W(3,1)=B
C      NINT=2         [INITIAL PARTITION HAS 2 INTERVALS]
C      RST=.FALSE.
C      EPS=.001
C      NMAX=50
C
C    1 CALL DQ1DAX(F,A,B,EPS,R,E,NINT,RST,W,NMAX,FMIN,FMAX,KF,IFLAG)
C
C      IF(EPS.EQ. .0001 .OR. IFLAG.GE.3)STOP
C      RST=.TRUE
C      EPS=.0001      [ASK FOR ANOTHER DIGIT]
C      GO TO 1
C      END
C      DOUBLE PRECISION FUNCTION F(X)
C      DOUBLE PRECISION X
C      IF(X.LT. .3)
C     1  THEN
C          F=X**(0.2)*ALOG(X)
C        ELSE
C          F=SIN(X)
C      ENDIF
C      RETURN
C      END
C
C
C            R E M A R K
C
C               WHEN YOU USE DQ1ADX WITH NINT=1, WE HAVE BUILT A SMALL
C            AMOUNT OF RANDOMIZATION INTO IT.  REPEATED CALLS DURING
C            THE SAME RUN WILL PRODUCE DIFFERENT, BUT HOPEFULLY
C            CONSISTENT, RESULTS.
C
C    E N D   O F   D O C U M E N T A T I O N
C***REFERENCES  (NONE)
C***ROUTINES CALLED (D1MACH,UNI,DGL15T,IDAMAX)
C***END PROLOGUE DQ1DAX
C
      DOUBLE PRECISION A,B,E,EB,EPMACH,EPS,F,FMAX,FMAXL,FMAXR,FMIN,
     *  FMINL,FMINR,FMN,FMX,R,D1MACH,RAB,RABS,RAV,T,TE,TE1,TE2,TR,
     *  TR1,TR2,UFLOW,W(NMAX,6),XM
      INTEGER C
      EXTERNAL F
      LOGICAL RST
C
C***FIRST EXECUTABLE STATEMENT  DQ1DAX
      EPMACH = D1MACH(4)
      UFLOW = D1MACH(1)
      MXTRY=NMAX/2
C          In case there is no more room, we can toss out easy intervals,
C             at most MXTRY times.
      IF(A.EQ.B) THEN
          R=0.
          E=0.
          NINT=0
          IFLAG=0
          KF=1
          FMIN=F(A)
          FMAX=FMIN
          GO TO 20
      ENDIF
      IF(RST) THEN
         IF(IFLAG.LT.3) THEN
           EB=MAX(100.*UFLOW,MAX(EPS,50.*EPMACH)*ABS(R))
           DO 19 I=1,NINT
               IF(ABS(W(I,3)).GT.(EB*(W(I,2)-W(I,1))/(B-A)))THEN
                                 W(I,3)=ABS(W(I,3))
               ELSE
                                 W(I,3)=-ABS(W(I,3))
               ENDIF
   19      CONTINUE
           GOTO 15
         ELSE
           GOTO 20
         ENDIF
      ENDIF
      KF=0
      IF(EPS .LE. 0. .OR. NINT .LE. 0 .OR. NINT .GE. NMAX) THEN
          IFLAG=7
          GO TO 20
      ENDIF
      IF(NINT.EQ.1)
     1  THEN
          W(1,1)=A
          W(2,2)=B
          W(1,5)=A
          W(1,6)=B
          W(2,5)=A
          W(2,6)=B
C          SELECT FIRST SUBDIVISION RANDOMLY
          W(1,2)=A+(B-A)/2.*(2*UNI()+7.)/8.
          W(2,1)=W(1,2)
          NINT=2
        ELSE
          IF(W(1,1).NE.A .OR. W(NINT+1,1).NE.B) THEN
               IFLAG=8
               GO TO 20
          ENDIF
          W(1,5)=A
          DO 89 I=1,NINT
             W(I,2)=W(I+1,1)
             W(I,5)=W(I,1)
             W(I,6)=W(I,2)
   89     CONTINUE
      ENDIF
C
      IFLAG = 0
      IROFF=0
      RABS=0.0
      DO 3 I=1,NINT
          CALL DGL15T(F,W(I,1),W(I,2),W(I,5),W(I,6),
     1          W(I,4),W(I,3),RAB,RAV,FMN,FMX)
          KF=KF+15
          IF(I.EQ.1)
     1       THEN
               R=W(I,4)
               E=W(I,3)
               RABS=RABS+RAB
               FMIN=FMN
               FMAX=FMX
             ELSE
               R=R+W(I,4)
               E=E+W(I,3)
               RABS=RABS+RAB
               FMAX=MAX(FMAX,FMX)
               FMIN=MIN(FMIN,FMN)
          ENDIF
    3 CONTINUE
      DO 10 I=NINT+1,NMAX
          W(I,3) = 0.
   10 CONTINUE
   15 CONTINUE
C
C   MAIN SUBPROGRAM LOOP
C
      IF(100.*EPMACH*RABS.GE.ABS(R) .AND. E.LT.EPS)GO TO 20
      EB=MAX(100.*UFLOW,MAX(EPS,50.*EPMACH)*ABS(R))
      IF(E.LE.EB) GO TO 20
      IF (NINT.LT.NMAX)
     1 THEN
        NINT = NINT+1
        C = NINT
       ELSE
        C=0
   16   IF(C.EQ.NMAX .OR. MXTRY.LE.0) THEN
            IFLAG=5
            GO TO 20
        ENDIF
        C=C+1
        IF(W(C,3).GT.0.0) GO TO 16
C            Found an interval to throw out
        MXTRY=MXTRY-1
      END IF
      LOC=IDAMAX(NINT,W(1,3),1)
      XM = W(LOC,1)+(W(LOC,2)-W(LOC,1))/2.
      IF ((MAX(ABS(W(LOC,1)),ABS(W(LOC,2)))).GT.
     1   ((1.+100.*EPMACH)*(ABS(XM)+0.1D+04*UFLOW)))
     2    THEN
            CALL DGL15T(F,W(LOC,1),XM,W(LOC,5),W(LOC,6),
     1                    TR1,TE1,RAB,RAV,FMINL,FMAXL)
            KF=KF+15
            IF (TE1.LT.(EB*(XM-W(LOC,1))/(B-A))) TE1=-TE1
            CALL DGL15T(F,XM,W(LOC,2),W(LOC,5),W(LOC,6),
     1                    TR2,TE2,RAB,RAV,FMINR,FMAXR)
            KF=KF+15
            FMIN=MIN(FMIN,FMINL,FMINR)
            FMAX=MAX(FMAX,FMAXL,FMAXR)
            IF (TE2.LT.(EB*(W(LOC,2)-XM)/(B-A))) TE2=-TE2
            TE = ABS(W(LOC,3))
            TR = W(LOC,4)
            W(C,3) = TE2
            W(C,4) = TR2
            W(C,1) = XM
            W(C,2) = W(LOC,2)
            W(C,5) = W(LOC,5)
            W(C,6) = W(LOC,6)
            W(LOC,3) = TE1
            W(LOC,4) = TR1
            W(LOC,2) = XM
            E = E-TE+(ABS(TE1)+ABS(TE2))
            R = R-TR+(TR1+TR2)
            IF(ABS(ABS(TE1)+ABS(TE2)-TE).LT.0.001*TE) THEN
                IROFF=IROFF+1
                IF(IROFF.GE.10) THEN
                     IFLAG=4
                     GO TO 20
                ENDIF
            ENDIF
          ELSE
            IF (EB.GT.W(LOC,3))
     1         THEN
                   W(LOC,3) = 0.
               ELSE
                   IFLAG=6
                   GO TO 20
            END IF
      END IF
      GO TO 15
C
C        ALL EXITS FROM HERE
C
   20 CONTINUE
      IF(IFLAG.GE.4)RETURN
      IFLAG=3
      T=EPS*ABS(R)
      IF(E.GT.EPS .AND. E.GT.T)RETURN
      IFLAG=2
      IF(E.GT.EPS .AND. E.LT.T)RETURN
      IFLAG=1
      IF(E.LT.EPS .AND. E.GT.T)RETURN
      IFLAG=0
      RETURN
      END
      SUBROUTINE DGL15T(F,A,B,XL,XR,R,AE,RA,
     1                 RASC,FMIN,FMAX)
C
C***AUTHORS          ROBERT PIESSENS AND ELISE DE DONCKER
C                    APPL. MATH. AND PROGR. DIV.- K.U.LEUVEN
C                    DAVID KAHANER, NBS WASHINGTON
C
C***PURPOSE
C              TO COMPUTE I = INTEGRAL OF G(X) OVER (A,B),
C                             WITH ERROR ESTIMATE
C                         J = INTEGRAL OF ABS(G) OVER (A,B)
C              DOUBLE PRECISION VERSION OF GL15T
C
C***DESCRIPTION
C           PARAMETERS
C            ON ENTRY
C              F      - DOUBLE PRECISION
C                       FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
C                       FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS
C                       TO BE DECLARED E X T E R N A L IN THE
C                       CALLING PROGRAM.
C                       THE FUNCTION G(X) IS DEFINED TO BE
C                       G(X)=F(PHI(X))*PHIP(X)
C                       WHERE PHI(X) IS THE CUBIC GIVEN BY
C                       THE ARITHMETIC STATEMENT FUNCTION BELOW.
C                       PHIP(X) IS ITS DERIVATIVE.  THE VARIABLES
C                       XL AND XR ARE THE LEFT AND RIGHT ENDPOINTS
C                       OF A PARENT INTERVAL OF WHICH (A,B) IS A PART.
C
C              A      - DOUBLE PRECISION
C                       LOWER LIMIT OF INTEGRATION
C
C              B      - DOUBLE PRECISION
C                       UPPER LIMIT OF INTEGRATION
C
C              XL     - DOUBLE PRECISION
C              XR     - DOUBLE PRECISION
C                       LOWER AND UPPER LIMITS OF PARENT INTERVAL
C                       OF WHICH [A,B] IS A PART.
C
C            ON RETURN
C              R - DOUBLE PRECISION
C                       APPROXIMATION TO THE INTEGRAL I
C                       R IS COMPUTED BY APPLYING THE 15-POINT
C                       KRONROD RULE (RESK) OBTAINED BY OPTIMAL
C                       ADDITION OF ABSCISSAE TO THE 7-POINT GAUSS
C                       RULE (RESG).
C
C              AE - DOUBLE PRECISION
C                       ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
C                       WHICH SHOULD NOT EXCEED ABS(I-R)
C
C              RA - DOUBLE PRECISION
C                       APPROXIMATION TO THE INTEGRAL J
C
C              RASC - DOUBLE PRECISION
C                       APPROXIMATION TO THE INTEGRAL OF ABS(G-I/(B-A))
C                       OVER (A,B)
C
C              FMAX, FMIN - DOUBLE PRECISION
C                       MAX AND MIN VALUES OF THE FUNCTION F ON (A,B)
C           SUBROUTINES OR FUNCTIONS NEEDED
C                 - F (USER-PROVIDED FUNCTION)
C                 - DOUBLE PRECISION
C                 - FORTRAN ABS, MAX, MIN
C
C     .................................................................
C***END PROLOGUE  DGL15T
C
      SAVE EPMACH,UFLOW
      DOUBLE PRECISION A,AE,B,DHLGTH,EPMACH,F,FC,FMAX,FMIN,FSUM,FV1,FV2,
     *  FVAL1,FVAL2,
     *  HLGTH,PHI,PHIP,PHIU,R,D1MACH,RA,RASC,RESG,RESK,RESKH,SL,SR,
     *  UFLOW,WG,WGK,XGK
      DOUBLE PRECISION XL,XR,CENTR,ABSC,U
      INTEGER J,JTW,JTWM1
C
      DIMENSION FV1(7),FV2(7),WG(4),WGK(8),XGK(8)
C
C           THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1)
C           BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR
C           CORRESPONDING WEIGHTS ARE GIVEN.
C
C           XGK    - ABSCISSAE OF THE 15-POINT KRONROD RULE
C                    XGK(2), XGK(4), ...  ABSCISSAE OF THE 7-POINT
C                    GAUSS RULE
C                    XGK(1), XGK(3), ...  ABSCISSAE WHICH ARE OPTIMALLY
C                    ADDED TO THE 7-POINT GAUSS RULE
C
C           WGK    - WEIGHTS OF THE 15-POINT KRONROD RULE
C
C           WG     - WEIGHTS OF THE 7-POINT GAUSS RULE
C
      DATA WG  (  1) / 0.1294849661 6886969327 0611432679 082 D0 /
      DATA WG  (  2) / 0.2797053914 8927666790 1467771423 780 D0 /
      DATA WG  (  3) / 0.3818300505 0511894495 0369775488 975 D0 /
      DATA WG  (  4) / 0.4179591836 7346938775 5102040816 327 D0 /
C
      DATA XGK (  1) / 0.9914553711 2081263920 6854697526 329 D0 /
      DATA XGK (  2) / 0.9491079123 4275852452 6189684047 851 D0 /
      DATA XGK (  3) / 0.8648644233 5976907278 9712788640 926 D0 /
      DATA XGK (  4) / 0.7415311855 9939443986 3864773280 788 D0 /
      DATA XGK (  5) / 0.5860872354 6769113029 4144838258 730 D0 /
      DATA XGK (  6) / 0.4058451513 7739716690 6606412076 961 D0 /
      DATA XGK (  7) / 0.2077849550 0789846760 0689403773 245 D0 /
      DATA XGK (  8) / 0.0000000000 0000000000 0000000000 000 D0 /
C
      DATA WGK (  1) / 0.0229353220 1052922496 3732008058 970 D0 /
      DATA WGK (  2) / 0.0630920926 2997855329 0700663189 204 D0 /
      DATA WGK (  3) / 0.1047900103 2225018383 9876322541 518 D0 /
      DATA WGK (  4) / 0.1406532597 1552591874 5189590510 238 D0 /
      DATA WGK (  5) / 0.1690047266 3926790282 6583426598 550 D0 /
      DATA WGK (  6) / 0.1903505780 6478540991 3256402421 014 D0 /
      DATA WGK (  7) / 0.2044329400 7529889241 4161999234 649 D0 /
      DATA WGK (  8) / 0.2094821410 8472782801 2999174891 714 D0 /
C
C
      PHI(U)=XR-(XR-XL)*U*U*(2.*U+3.)
      PHIP(U)=-6.*U*(U+1.)
C
C           LIST OF MAJOR VARIABLES
C           -----------------------
C
C           CENTR  - MID POINT OF THE INTERVAL
C           HLGTH  - HALF-LENGTH OF THE INTERVAL
C           ABSC   - ABSCISSA
C           FVAL*  - FUNCTION VALUE
C           RESG   - R OF THE 7-POINT GAUSS FORMULA
C           RESK   - R OF THE 15-POINT KRONROD FORMULA
C           RESKH  - APPROXIMATION TO THE MEAN VALUE OF F OVER (A,B),
C                    I.E. TO I/(B-A)
C
C           MACHINE DEPENDENT CONSTANTS
C           ---------------------------
C
C           EPMACH IS THE LARGEST RELATIVE SPACING.
C           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
C
C***FIRST EXECUTABLE STATEMENT
      DATA EPMACH,UFLOW/0.0,0.0/
      IF(EPMACH.EQ.0.0) THEN
         EPMACH=D1MACH(4)
         UFLOW=D1MACH(1)
      ENDIF
C
      IF(XL.LT.XR)THEN
         SL=(XL)
         SR=(XR)
        ELSE
         SL=(XR)
         SR=(XL)
      ENDIF
      HLGTH = 0.5D+00*(B-A)
      CENTR = A+HLGTH
      DHLGTH = ABS(HLGTH)
C
C           COMPUTE THE 15-POINT KRONROD APPROXIMATION TO
C           THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR.
C
      U=(CENTR-XR)/(XR-XL)
      PHIU=PHI(U)
      IF(PHIU.LE.SL .OR. PHIU.GE.SR) PHIU=CENTR
      FMIN=F(PHIU)
      FMAX=FMIN
      FC=FMIN*PHIP(U)
      RESG = FC*WG(4)
      RESK = FC*WGK(8)
      RA = ABS(RESK)
      DO 10 J=1,3
        JTW = J*2
        ABSC = HLGTH*XGK(JTW)
        U=(CENTR-ABSC-XR)/(XR-XL)
        PHIU=PHI(U)
        IF(PHIU.LE.SL .OR. PHIU.GE.SR) PHIU=CENTR
        FVAL1=F(PHIU)
        FMAX=MAX(FMAX,FVAL1)
        FMIN=MIN(FMIN,FVAL1)
        FVAL1=FVAL1*PHIP(U)
        U=(CENTR+ABSC-XR)/(XR-XL)
        PHIU=PHI(U)
        IF(PHIU.LE.SL .OR. PHIU.GE.SR) PHIU=CENTR
        FVAL2=F(PHIU)
        FMAX=MAX(FMAX,FVAL2)
        FMIN=MIN(FMIN,FVAL2)
        FVAL2=FVAL2*PHIP(U)
        FV1(JTW) = FVAL1
        FV2(JTW) = FVAL2
        FSUM = FVAL1+FVAL2
        RESG = RESG+WG(J)*FSUM
        RESK = RESK+WGK(JTW)*FSUM
        RA = RA+WGK(JTW)*(ABS(FVAL1)+ABS(FVAL2))
   10 CONTINUE
      DO 15 J = 1,4
        JTWM1 = J*2-1
        ABSC = HLGTH*XGK(JTWM1)
        U=(CENTR-ABSC-XR)/(XR-XL)
        PHIU=PHI(U)
        IF(PHIU.LE.SL .OR. PHIU.GE.SR) PHIU=CENTR
        FVAL1=F(PHIU)
        FMAX=MAX(FMAX,FVAL1)
        FMIN=MIN(FMIN,FVAL1)
        FVAL1=FVAL1*PHIP(U)
        U=(CENTR+ABSC-XR)/(XR-XL)
        PHIU=PHI(U)
        IF(PHIU.LE.SL .OR. PHIU.GE.SR) PHIU=CENTR
        FVAL2=F(PHIU)
        FMAX=MAX(FMAX,FVAL2)
        FMIN=MIN(FMIN,FVAL2)
        FVAL2=FVAL2*PHIP(U)
        FV1(JTWM1) = FVAL1
        FV2(JTWM1) = FVAL2
        FSUM = FVAL1+FVAL2
        RESK = RESK+WGK(JTWM1)*FSUM
        RA = RA+WGK(JTWM1)*(ABS(FVAL1)+ABS(FVAL2))
   15 CONTINUE
      RESKH = RESK*0.5D+00
      RASC = WGK(8)*ABS(FC-RESKH)
      DO 20 J=1,7
        RASC = RASC+WGK(J)*(ABS(FV1(J)-RESKH)+ABS(FV2(J)-RESKH))
   20 CONTINUE
      R = RESK*HLGTH
      RA = RA*DHLGTH
      RASC = RASC*DHLGTH
      AE = ABS((RESK-RESG)*HLGTH)
      IF(RASC.NE.0.0D+00.AND.AE.NE.0.0D+00)
     *  AE = RASC*MIN(0.1D+01,
     *  (0.2D+03*AE/RASC)**1.5D+00)
      IF(RA.GT.UFLOW/(0.5D+02*EPMACH)) AE = MAX(
     *  (EPMACH*0.5D+02)*RA,AE)
      RETURN
      END
      SUBROUTINE DQK15(F,A,B,RESULT,ABSERR,RESABS,RESASC)
C***BEGIN PROLOGUE  DQK15
C***DATE WRITTEN   800101   (YYMMDD)
C***REVISION DATE  870530   (YYMMDD)
C***CATEGORY NO.  H2A1A2
C***KEYWORDS  15-POINT GAUSS-KRONROD RULES
C***AUTHOR  PIESSENS, ROBERT, AND DE DONCKER, ELISE,
C             APPLIED MATH. AND PROGR. DIV. - K. U. LEUVEN
C***PURPOSE  To compute I = Integral of F over (A,B), with error estimate,
C               and     J = integral of ABS(F) over (A,B)
C***DESCRIPTION
C           Double precision version
C
C           PARAMETERS ON ENTRY
C              F      - Double precision
C                       Function subprogram defining the integrand
C                       FUNCTION F(X). The actual name for F needs to be
C                       Declared E X T E R N A L in the calling program.
C
C              A      - Double precision: Lower limit of integration
C
C              B      - Double precision: Upper limit of integration
C
C            PARAMETERS ON RETURN
C              RESULT - Double precision: Approximation to the integral I
C                       Result is computed by applying the 15-POINT
C                       KRONROD RULE (RESK) obtained by optimal addition
C                       of abscissae to the7-POINT GAUSS RULE(RESG).
C
C              ABSERR - Double precision: Estimate of the modulus of the
C                       absolute error, which should not exceed ABS(I-RESULT)
C
C              RESABS - Double precision: Approximation to the integral J
C
C              RESASC - Double precision: Approximation to the integral of
C                       ABS(F-I/(B-A))  over (A,B)
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH
C***END PROLOGUE  DQK15
C
      DOUBLE PRECISION A,ABSC,ABSERR,B,CENTR,DABS,DHLGTH,DMAX1,DMIN1,
     1  D1MACH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,RESABS,RESASC,
     2  RESG,RESK,RESKH,RESULT,UFLOW,WG,WGK,XGK
      INTEGER J,JTW,JTWM1
      EXTERNAL F
C
      DIMENSION FV1(7),FV2(7),WG(4),WGK(8),XGK(8)
C
C           THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1).
C           BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR
C           CORRESPONDING WEIGHTS ARE GIVEN.
C
C           XGK    - ABSCISSAE OF THE 15-POINT KRONROD RULE
C                    XGK(2), XGK(4), ...  ABSCISSAE OF THE 7-POINT
C                    GAUSS RULE
C                    XGK(1), XGK(3), ...  ABSCISSAE WHICH ARE OPTIMALLY
C                    ADDED TO THE 7-POINT GAUSS RULE
C
C           WGK    - WEIGHTS OF THE 15-POINT KRONROD RULE
C
C           WG     - WEIGHTS OF THE 7-POINT GAUSS RULE
C
C
C GAUSS QUADRATURE WEIGHTS AND KRONRON QUADRATURE ABSCISSAE AND WEIGHTS
C AS EVALUATED WITH 80 DECIMAL DIGIT ARITHMETIC BY L. W. FULLERTON,
C BELL LABS, NOV. 1981.
C
      DATA WG  (  1) / 0.1294849661 6886969327 0611432679 082 D0 /
      DATA WG  (  2) / 0.2797053914 8927666790 1467771423 780 D0 /
      DATA WG  (  3) / 0.3818300505 0511894495 0369775488 975 D0 /
      DATA WG  (  4) / 0.4179591836 7346938775 5102040816 327 D0 /
C
      DATA XGK (  1) / 0.9914553711 2081263920 6854697526 329 D0 /
      DATA XGK (  2) / 0.9491079123 4275852452 6189684047 851 D0 /
      DATA XGK (  3) / 0.8648644233 5976907278 9712788640 926 D0 /
      DATA XGK (  4) / 0.7415311855 9939443986 3864773280 788 D0 /
      DATA XGK (  5) / 0.5860872354 6769113029 4144838258 730 D0 /
      DATA XGK (  6) / 0.4058451513 7739716690 6606412076 961 D0 /
      DATA XGK (  7) / 0.2077849550 0789846760 0689403773 245 D0 /
      DATA XGK (  8) / 0.0000000000 0000000000 0000000000 000 D0 /
C
      DATA WGK (  1) / 0.0229353220 1052922496 3732008058 970 D0 /
      DATA WGK (  2) / 0.0630920926 2997855329 0700663189 204 D0 /
      DATA WGK (  3) / 0.1047900103 2225018383 9876322541 518 D0 /
      DATA WGK (  4) / 0.1406532597 1552591874 5189590510 238 D0 /
      DATA WGK (  5) / 0.1690047266 3926790282 6583426598 550 D0 /
      DATA WGK (  6) / 0.1903505780 6478540991 3256402421 014 D0 /
      DATA WGK (  7) / 0.2044329400 7529889241 4161999234 649 D0 /
      DATA WGK (  8) / 0.2094821410 8472782801 2999174891 714 D0 /
C
C
C           LIST OF MAJOR VARIABLES
C           -----------------------
C
C           CENTR  - MID POINT OF THE INTERVAL
C           HLGTH  - HALF-LENGTH OF THE INTERVAL
C           ABSC   - ABSCISSA
C           FVAL*  - FUNCTION VALUE
C           RESG   - RESULT OF THE 7-POINT GAUSS FORMULA
C           RESK   - RESULT OF THE 15-POINT KRONROD FORMULA
C           RESKH  - APPROXIMATION TO THE MEAN VALUE OF F OVER (A,B),
C                    I.E. TO I/(B-A)
C
C           MACHINE DEPENDENT CONSTANTS
C           ---------------------------
C
C           EPMACH IS THE LARGEST RELATIVE SPACING.
C           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
C
C***FIRST EXECUTABLE STATEMENT  DQK15
      EPMACH = D1MACH(4)
      UFLOW = D1MACH(1)
C
      CENTR = 0.5D+00*(A+B)
      HLGTH = 0.5D+00*(B-A)
      DHLGTH = DABS(HLGTH)
C
C           COMPUTE THE 15-POINT KRONROD APPROXIMATION TO
C           THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR.
C
      FC = F(CENTR)
      RESG = FC*WG(4)
      RESK = FC*WGK(8)
      RESABS = DABS(RESK)
      DO 10 J=1,3
        JTW = J*2
        ABSC = HLGTH*XGK(JTW)
        FVAL1 = F(CENTR-ABSC)
        FVAL2 = F(CENTR+ABSC)
        FV1(JTW) = FVAL1
        FV2(JTW) = FVAL2
        FSUM = FVAL1+FVAL2
        RESG = RESG+WG(J)*FSUM
        RESK = RESK+WGK(JTW)*FSUM
        RESABS = RESABS+WGK(JTW)*(DABS(FVAL1)+DABS(FVAL2))
   10 CONTINUE
      DO 15 J = 1,4
        JTWM1 = J*2-1
        ABSC = HLGTH*XGK(JTWM1)
        FVAL1 = F(CENTR-ABSC)
        FVAL2 = F(CENTR+ABSC)
        FV1(JTWM1) = FVAL1
        FV2(JTWM1) = FVAL2
        FSUM = FVAL1+FVAL2
        RESK = RESK+WGK(JTWM1)*FSUM
        RESABS = RESABS+WGK(JTWM1)*(DABS(FVAL1)+DABS(FVAL2))
   15 CONTINUE
      RESKH = RESK*0.5D+00
      RESASC = WGK(8)*DABS(FC-RESKH)
      DO 20 J=1,7
        RESASC = RESASC+WGK(J)*(DABS(FV1(J)-RESKH)+DABS(FV2(J)-RESKH))
   20 CONTINUE
      RESULT = RESK*HLGTH
      RESABS = RESABS*DHLGTH
      RESASC = RESASC*DHLGTH
      ABSERR = DABS((RESK-RESG)*HLGTH)
      IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.0D+00)
     1  ABSERR = RESASC*DMIN1(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00)
      IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = DMAX1
     1  ((EPMACH*0.5D+02)*RESABS,ABSERR)
      RETURN
      END