# dqawf.f

          SUBROUTINE DQAWF (F, A, OMEGA, INTEGR, EPSABS, RESULT, ABSERR,
         +   NEVAL, IER, LIMLST, LST, LENIW, MAXP1, LENW, IWORK, WORK)
    C***BEGIN PROLOGUE  DQAWF
    C***PURPOSE  The routine calculates an approximation result to a given
    C            Fourier integral I=Integral of F(X)*W(X) over (A,INFINITY)
    C            where W(X) = COS(OMEGA*X) or W(X) = SIN(OMEGA*X).
    C            Hopefully satisfying following claim for accuracy
    C            ABS(I-RESULT).LE.EPSABS.
    C***LIBRARY   SLATEC (QUADPACK)
    C***CATEGORY  H2A3A1
    C***TYPE      DOUBLE PRECISION (QAWF-S, DQAWF-D)
    C***KEYWORDS  AUTOMATIC INTEGRATOR, CONVERGENCE ACCELERATION,
    C             FOURIER INTEGRALS, INTEGRATION BETWEEN ZEROS, QUADPACK,
    C             QUADRATURE, SPECIAL-PURPOSE INTEGRAL
    C***AUTHOR  Piessens, Robert
    C             Applied Mathematics and Programming Division
    C             K. U. Leuven
    C           de Doncker, Elise
    C             Applied Mathematics and Programming Division
    C             K. U. Leuven
    C***DESCRIPTION
    C
    C        Computation of Fourier integrals
    C        Standard fortran subroutine
    C        Double precision version
    C
    C
    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 to be
    C                     declared E X T E R N A L in the driver program.
    C
    C            A      - Double precision
    C                     Lower limit of integration
    C
    C            OMEGA  - Double precision
    C                     Parameter in the integrand WEIGHT function
    C
    C            INTEGR - Integer
    C                     Indicates which of the WEIGHT functions is used
    C                     INTEGR = 1      W(X) = COS(OMEGA*X)
    C                     INTEGR = 2      W(X) = SIN(OMEGA*X)
    C                     IF INTEGR.NE.1.AND.INTEGR.NE.2, the routine
    C                     will end with IER = 6.
    C
    C            EPSABS - Double precision
    C                     Absolute accuracy requested, EPSABS.GT.0.
    C                     If EPSABS.LE.0, the routine will end with IER = 6.
    C
    C         ON RETURN
    C            RESULT - Double precision
    C                     Approximation to the integral
    C
    C            ABSERR - Double precision
    C                     Estimate of the modulus of the absolute error,
    C                     Which should equal or exceed ABS(I-RESULT)
    C
    C            NEVAL  - Integer
    C                     Number of integrand evaluations
    C
    C            IER    - Integer
    C                     IER = 0 Normal and reliable termination of the
    C                             routine. It is assumed that the requested
    C                             accuracy has been achieved.
    C                     IER.GT.0 Abnormal termination of the routine.
    C                             The estimates for integral and error are
    C                             less reliable. It is assumed that the
    C                             requested accuracy has not been achieved.
    C            ERROR MESSAGES
    C                    If OMEGA.NE.0
    C                     IER = 1 Maximum number of cycles allowed
    C                             has been achieved, i.e. of subintervals
    C                             (A+(K-1)C,A+KC) where
    C                             C = (2*INT(ABS(OMEGA))+1)*PI/ABS(OMEGA),
    C                             FOR K = 1, 2, ..., LST.
    C                             One can allow more cycles by increasing
    C                             the value of LIMLST (and taking the
    C                             according dimension adjustments into
    C                             account). Examine the array IWORK which
    C                             contains the error flags on the cycles, in
    C                             order to look for eventual local
    C                             integration difficulties.
    C                             If the position of a local difficulty
    C                             can be determined (e.g. singularity,
    C                             discontinuity within the interval) one
    C                             will probably gain from splitting up the
    C                             interval at this point and calling
    C                             appropriate integrators on the subranges.
    C                         = 4 The extrapolation table constructed for
    C                             convergence acceleration of the series
    C                             formed by the integral contributions over
    C                             the cycles, does not converge to within
    C                             the requested accuracy.
    C                             As in the case of IER = 1, it is advised
    C                             to examine the array IWORK which contains
    C                             the error flags on the cycles.
    C                         = 6 The input is invalid because
    C                             (INTEGR.NE.1 AND INTEGR.NE.2) or
    C                              EPSABS.LE.0 or LIMLST.LT.1 or
    C                              LENIW.LT.(LIMLST+2) or MAXP1.LT.1 or
    C                              LENW.LT.(LENIW*2+MAXP1*25).
    C                              RESULT, ABSERR, NEVAL, LST are set to
    C                              zero.
    C                         = 7 Bad integrand behaviour occurs within
    C                             one or more of the cycles. Location and
    C                             type of the difficulty involved can be
    C                             determined from the first LST elements of
    C                             vector IWORK.  Here LST is the number of
    C                             cycles actually needed (see below).
    C                             IWORK(K) = 1 The maximum number of
    C                                          subdivisions (=(LENIW-LIMLST)
    C                                          /2) has been achieved on the
    C                                          K th cycle.
    C                                      = 2 Occurrence of roundoff error
    C                                          is detected and prevents the
    C                                          tolerance imposed on the K th
    C                                          cycle, from being achieved
    C                                          on this cycle.
    C                                      = 3 Extremely bad integrand
    C                                          behaviour occurs at some
    C                                          points of the K th cycle.
    C                                      = 4 The integration procedure
    C                                          over the K th cycle does
    C                                          not converge (to within the
    C                                          required accuracy) due to
    C                                          roundoff in the extrapolation
    C                                          procedure invoked on this
    C                                          cycle. It is assumed that the
    C                                          result on this interval is
    C                                          the best which can be
    C                                          obtained.
    C                                      = 5 The integral over the K th
    C                                          cycle is probably divergent
    C                                          or slowly convergent. It must
    C                                          be noted that divergence can
    C                                          occur with any other value of
    C                                          IWORK(K).
    C                    If OMEGA = 0 and INTEGR = 1,
    C                    The integral is calculated by means of DQAGIE,
    C                    and IER = IWORK(1) (with meaning as described
    C                    for IWORK(K),K = 1).
    C
    C         DIMENSIONING PARAMETERS
    C            LIMLST - Integer
    C                     LIMLST gives an upper bound on the number of
    C                     cycles, LIMLST.GE.3.
    C                     If LIMLST.LT.3, the routine will end with IER = 6.
    C
    C            LST    - Integer
    C                     On return, LST indicates the number of cycles
    C                     actually needed for the integration.
    C                     If OMEGA = 0, then LST is set to 1.
    C
    C            LENIW  - Integer
    C                     Dimensioning parameter for IWORK. On entry,
    C                     (LENIW-LIMLST)/2 equals the maximum number of
    C                     subintervals allowed in the partition of each
    C                     cycle, LENIW.GE.(LIMLST+2).
    C                     If LENIW.LT.(LIMLST+2), the routine will end with
    C                     IER = 6.
    C
    C            MAXP1  - Integer
    C                     MAXP1 gives an upper bound on the number of
    C                     Chebyshev moments which can be stored, i.e. for
    C                     the intervals of lengths ABS(B-A)*2**(-L),
    C                     L = 0,1, ..., MAXP1-2, MAXP1.GE.1.
    C                     If MAXP1.LT.1, the routine will end with IER = 6.
    C            LENW   - Integer
    C                     Dimensioning parameter for WORK
    C                     LENW must be at least LENIW*2+MAXP1*25.
    C                     If LENW.LT.(LENIW*2+MAXP1*25), the routine will
    C                     end with IER = 6.
    C
    C         WORK ARRAYS
    C            IWORK  - Integer
    C                     Vector of dimension at least LENIW
    C                     On return, IWORK(K) FOR K = 1, 2, ..., LST
    C                     contain the error flags on the cycles.
    C
    C            WORK   - Double precision
    C                     Vector of dimension at least
    C                     On return,
    C                     WORK(1), ..., WORK(LST) contain the integral
    C                      approximations over the cycles,
    C                     WORK(LIMLST+1), ..., WORK(LIMLST+LST) contain
    C                      the error estimates over the cycles.
    C                     further elements of WORK have no specific
    C                     meaning for the user.
    C
    C***REFERENCES  (NONE)
    C***ROUTINES CALLED  DQAWFE, XERMSG
    C***REVISION HISTORY  (YYMMDD)
    C   800101  DATE WRITTEN
    C   890831  Modified array declarations.  (WRB)
    C   891009  Removed unreferenced variable.  (WRB)
    C   891009  REVISION DATE from Version 3.2
    C   891214  Prologue converted to Version 4.0 format.  (BAB)
    C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
    C***END PROLOGUE  DQAWF