subroutine dqawf(f,a,omega,integr,epsabs,result,abserr,neval,ier, * limlst,lst,leniw,maxp1,lenw,iwork,work) c***begin prologue dqawf c***date written 800101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a3a1 c***keywords automatic integrator, special-purpose,fourier c integral, integration between zeros with dqawoe, c convergence acceleration with dqelg c***author piessens,robert ,appl. math. & progr. div. - k.u.leuven c de doncker,elise,appl. math & progr. div. - k.u.leuven 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***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 accelaration 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 extimates over the cycles. c further elements of work have no specific c meaning for the user. c c***references (none) c***routines called dqawfe,xerror c***end prologue dqawf c double precision a,abserr,epsabs,f,omega,result,work integer ier,integr,iwork,last,leniw,lenw,limit,limlst,ll2,lvl, * lst,l1,l2,l3,l4,l5,l6,maxp1,neval c dimension iwork(leniw),work(lenw) c external f c c check validity of limlst, leniw, maxp1 and lenw. c c***first executable statement dqawf ier = 6 neval = 0 last = 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. * (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, * abserr,neval,ier,work(1),work(l1),iwork(1),lst,work(l2), * 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 xerror(26habnormal return from dqawf,26,ier,lvl) return end