*DECK DQAWF 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 C DOUBLE PRECISION A,ABSERR,EPSABS,F,OMEGA,RESULT,WORK INTEGER IER,INTEGR,IWORK,LENIW,LENW,LIMIT,LIMLST,LL2,LVL, 1 LST,L1,L2,L3,L4,L5,L6,MAXP1,NEVAL C DIMENSION IWORK(*),WORK(*) C EXTERNAL F C C CHECK VALIDITY OF LIMLST, LENIW, MAXP1 AND LENW. C C***FIRST EXECUTABLE STATEMENT DQAWF IER = 6 NEVAL = 0 RESULT = 0.0D+00 ABSERR = 0.0D+00 IF(LIMLST.LT.3.OR.LENIW.LT.(LIMLST+2).OR.MAXP1.LT.1.OR.LENW.LT. 1 (LENIW*2+MAXP1*25)) GO TO 10 C C PREPARE CALL FOR DQAWFE C LIMIT = (LENIW-LIMLST)/2 L1 = LIMLST+1 L2 = LIMLST+L1 L3 = LIMIT+L2 L4 = LIMIT+L3 L5 = LIMIT+L4 L6 = LIMIT+L5 LL2 = LIMIT+L1 CALL DQAWFE(F,A,OMEGA,INTEGR,EPSABS,LIMLST,LIMIT,MAXP1,RESULT, 1 ABSERR,NEVAL,IER,WORK(1),WORK(L1),IWORK(1),LST,WORK(L2), 2 WORK(L3),WORK(L4),WORK(L5),IWORK(L1),IWORK(LL2),WORK(L6)) C C CALL ERROR HANDLER IF NECESSARY C LVL = 0 10 IF(IER.EQ.6) LVL = 1 IF (IER .NE. 0) CALL XERMSG ('SLATEC', 'DQAWF', + 'ABNORMAL RETURN', IER, LVL) RETURN END