SLATEC Routines --- QAWO ---
*DECK QAWO
SUBROUTINE QAWO (F, A, B, OMEGA, INTEGR, EPSABS, EPSREL, RESULT,
+ ABSERR, NEVAL, IER, LENIW, MAXP1, LENW, LAST, IWORK, WORK)
C***BEGIN PROLOGUE QAWO
C***PURPOSE Calculate an approximation to a given definite integral
C I = Integral of F(X)*W(X) over (A,B), where
C W(X) = COS(OMEGA*X)
C or W(X) = SIN(OMEGA*X),
C hopefully satisfying the following claim for accuracy
C ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
C***LIBRARY SLATEC (QUADPACK)
C***CATEGORY H2A2A1
C***TYPE SINGLE PRECISION (QAWO-S, DQAWO-D)
C***KEYWORDS AUTOMATIC INTEGRATOR, CLENSHAW-CURTIS METHOD,
C EXTRAPOLATION, GLOBALLY ADAPTIVE,
C INTEGRAND WITH OSCILLATORY COS OR SIN FACTOR, QUADPACK,
C QUADRATURE, SPECIAL-PURPOSE
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 oscillatory integrals
C Standard fortran subroutine
C Real version
C
C PARAMETERS
C ON ENTRY
C F - Real
C Function subprogram defining the function
C 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 - Real
C Lower limit of integration
C
C B - Real
C Upper limit of integration
C
C OMEGA - Real
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 will
C end with IER = 6.
C
C EPSABS - Real
C Absolute accuracy requested
C EPSREL - Real
C Relative accuracy requested
C If EPSABS.LE.0 and
C EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
C the routine will end with IER = 6.
C
C ON RETURN
C RESULT - Real
C Approximation to the integral
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 - 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 IER = 1 Maximum number of subdivisions allowed
C has been achieved(= LENIW/2). One can
C allow more subdivisions by increasing the
C value of LENIW (and taking the according
C dimension adjustments into account).
C However, if this yields no improvement it
C is advised to analyze the integrand in
C order to determine the 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 the integrator on the
C subranges. If possible, an appropriate
C special-purpose integrator should be used
C which is designed for handling the type of
C difficulty involved.
C = 2 The occurrence of roundoff error is
C detected, which prevents the requested
C tolerance from being achieved.
C The error may be under-estimated.
C = 3 Extremely bad integrand behaviour occurs
C at some interior points of the
C integration interval.
C = 4 The algorithm does not converge.
C Roundoff error is detected in the
C extrapolation table. It is presumed that
C the requested tolerance cannot be achieved
C due to roundoff in the extrapolation
C table, and that the returned result is
C the best which can be obtained.
C = 5 The integral is probably divergent, or
C slowly convergent. It must be noted that
C divergence can occur with any other value
C of IER.
C = 6 The input is invalid, because
C (EPSABS.LE.0 and
C EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28))
C or (INTEGR.NE.1 AND INTEGR.NE.2),
C or LENIW.LT.2 OR MAXP1.LT.1 or
C LENW.LT.LENIW*2+MAXP1*25.
C RESULT, ABSERR, NEVAL, LAST are set to
C zero. Except when LENIW, MAXP1 or LENW are
C invalid, WORK(LIMIT*2+1), WORK(LIMIT*3+1),
C IWORK(1), IWORK(LIMIT+1) are set to zero,
C WORK(1) is set to A and WORK(LIMIT+1) to
C B.
C
C DIMENSIONING PARAMETERS
C LENIW - Integer
C Dimensioning parameter for IWORK.
C LENIW/2 equals the maximum number of subintervals
C allowed in the partition of the given integration
C interval (A,B), LENIW.GE.2.
C If LENIW.LT.2, the routine will end with IER = 6.
C
C MAXP1 - Integer
C Gives an upper bound on the number of Chebyshev
C moments which can be stored, i.e. for the
C 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
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 LAST - Integer
C On return, LAST equals the number of subintervals
C produced in the subdivision process, which
C determines the number of significant elements
C actually in the WORK ARRAYS.
C
C WORK ARRAYS
C IWORK - Integer
C Vector of dimension at least LENIW
C on return, the first K elements of which contain
C pointers to the error estimates over the
C subintervals, such that WORK(LIMIT*3+IWORK(1)), ..
C WORK(LIMIT*3+IWORK(K)) form a decreasing
C sequence, with LIMIT = LENW/2 , and K = LAST
C if LAST.LE.(LIMIT/2+2), and K = LIMIT+1-LAST
C otherwise.
C Furthermore, IWORK(LIMIT+1), ..., IWORK(LIMIT+
C LAST) indicate the subdivision levels of the
C subintervals, such that IWORK(LIMIT+I) = L means
C that the subinterval numbered I is of length
C ABS(B-A)*2**(1-L).
C
C WORK - Real
C Vector of dimension at least LENW
C On return
C WORK(1), ..., WORK(LAST) contain the left
C end points of the subintervals in the
C partition of (A,B),
C WORK(LIMIT+1), ..., WORK(LIMIT+LAST) contain
C the right end points,
C WORK(LIMIT*2+1), ..., WORK(LIMIT*2+LAST) contain
C the integral approximations over the
C subintervals,
C WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST)
C contain the error estimates.
C WORK(LIMIT*4+1), ..., WORK(LIMIT*4+MAXP1*25)
C Provide space for storing the Chebyshev moments.
C Note that LIMIT = LENW/2.
C
C***REFERENCES (NONE)
C***ROUTINES CALLED QAWOE, XERMSG
C***REVISION HISTORY (YYMMDD)
C 800101 DATE WRITTEN
C 890831 Modified array declarations. (WRB)
C 890831 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 QAWO