# qawfe.f

          SUBROUTINE QAWFE (F, A, OMEGA, INTEGR, EPSABS, LIMLST, LIMIT,
         +   MAXP1, RESULT, ABSERR, NEVAL, IER, RSLST, ERLST, IERLST, LST,
         +   ALIST, BLIST, RLIST, ELIST, IORD, NNLOG, CHEBMO)
    C***BEGIN PROLOGUE  QAWFE
    C***PURPOSE  The routine calculates an approximation result to a
    C            given Fourier integral
    C            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      SINGLE PRECISION (QAWFE-S, DQAWFE-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        Real version
    C
    C        PARAMETERS
    C         ON ENTRY
    C            F      - Real
    C                     Function subprogram defining the integrand
    C                     Function F(X). The actual name for F needs to
    C                     be declared E X T E R N A L in the driver program.
    C
    C            A      - Real
    C                     Lower limit of integration
    C
    C            OMEGA  - Real
    C                     Parameter in the WEIGHT function
    C
    C            INTEGR - Integer
    C                     Indicates which WEIGHT function 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 will
    C                     end with IER = 6.
    C
    C            EPSABS - Real
    C                     absolute accuracy requested, EPSABS.GT.0
    C                     If EPSABS.LE.0, the routine will end with IER = 6.
    C
    C            LIMLST - Integer
    C                     LIMLST gives an upper bound on the number of
    C                     cycles, LIMLST.GE.1.
    C                     If LIMLST.LT.3, the routine will end with IER = 6.
    C
    C            LIMIT  - Integer
    C                     Gives an upper bound on the number of subintervals
    C                     allowed in the partition of each cycle, LIMIT.GE.1
    C                     each cycle, LIMIT.GE.1.
    C
    C            MAXP1  - Integer
    C                     Gives an upper bound on the number of
    C                     Chebyshev moments which can be stored, I.E.
    C                     for the intervals of lengths ABS(B-A)*2**(-L),
    C                     L=0,1, ..., MAXP1-2, MAXP1.GE.1
    C
    C         ON RETURN
    C            RESULT - Real
    C                     Approximation to the integral X
    C
    C            ABSERR - Real
    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    - IER = 0 Normal and reliable termination of
    C                             the routine. It is assumed that the
    C                             requested accuracy has been achieved.
    C                     IER.GT.0 Abnormal termination of the routine. The
    C                             estimates for integral and error are less
    C                             reliable. It is assumed that the requested
    C                             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).
    C                             Examine the array IWORK which contains
    C                             the error flags on the cycles, in order to
    C                             look for eventual local integration
    C                             difficulties. If the position of a local
    C                             difficulty can be determined (e.g.
    C                             SINGULARITY, DISCONTINUITY within the
    C                             interval) one will probably gain from
    C                             splitting up the interval at this point
    C                             and calling appropriate integrators on
    C                             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. As in the case of
    C                             IER = 1, it is advised to examine the
    C                             array IWORK which contains the error
    C                             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.3.
    C                              RESULT, ABSERR, NEVAL, LST are set
    C                              to zero.
    C                         = 7 Bad integrand behaviour occurs within one
    C                             or more of the cycles. Location and type
    C                             of the difficulty involved can be
    C                             determined from the vector IERLST. Here
    C                             LST is the number of cycles actually
    C                             needed (see below).
    C                             IERLST(K) = 1 The maximum number of
    C                                           subdivisions (= LIMIT) has
    C                                           been achieved on the K th
    C                                           cycle.
    C                                       = 2 Occurrence of roundoff error
    C                                           is detected and prevents the
    C                                           tolerance imposed on the
    C                                           K th cycle, from being
    C                                           achieved.
    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
    C                                           extrapolation procedure
    C                                           invoked on this cycle. It
    C                                           is assumed that the result
    C                                           on this interval is the
    C                                           best which can be obtained.
    C                                       = 5 The integral over the K th
    C                                           cycle is probably divergent
    C                                           or slowly convergent. It
    C                                           must be noted that
    C                                           divergence can occur with
    C                                           any other value of
    C                                           IERLST(K).
    C                    If OMEGA = 0 and INTEGR = 1,
    C                    The integral is calculated by means of DQAGIE
    C                    and IER = IERLST(1) (with meaning as described
    C                    for IERLST(K), K = 1).
    C
    C            RSLST  - Real
    C                     Vector of dimension at least LIMLST
    C                     RSLST(K) contains the integral contribution
    C                     over the interval (A+(K-1)C,A+KC) where
    C                     C = (2*INT(ABS(OMEGA))+1)*PI/ABS(OMEGA),
    C                     K = 1, 2, ..., LST.
    C                     Note that, if OMEGA = 0, RSLST(1) contains
    C                     the value of the integral over (A,INFINITY).
    C
    C            ERLST  - Real
    C                     Vector of dimension at least LIMLST
    C                     ERLST(K) contains the error estimate corresponding
    C                     with RSLST(K).
    C
    C            IERLST - Integer
    C                     Vector of dimension at least LIMLST
    C                     IERLST(K) contains the error flag corresponding
    C                     with RSLST(K). For the meaning of the local error
    C                     flags see description of output parameter IER.
    C
    C            LST    - Integer
    C                     Number of subintervals needed for the integration
    C                     If OMEGA = 0 then LST is set to 1.
    C
    C            ALIST, BLIST, RLIST, ELIST - Real
    C                     vector of dimension at least LIMIT,
    C
    C            IORD, NNLOG - Integer
    C                     Vector of dimension at least LIMIT, providing
    C                     space for the quantities needed in the subdivision
    C                     process of each cycle
    C
    C            CHEBMO - Real
    C                     Array of dimension at least (MAXP1,25), providing
    C                     space for the Chebyshev moments needed within the
    C                     cycles
    C
    C***REFERENCES  (NONE)
    C***ROUTINES CALLED  QAGIE, QAWOE, QELG, R1MACH
    C***REVISION HISTORY  (YYMMDD)
    C   800101  DATE WRITTEN
    C   890531  Changed all specific intrinsics to generic.  (WRB)
    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***END PROLOGUE  QAWFE