      subroutine dqawo(f,a,b,omega,integr,epsabs,epsrel,result,abserr,
     *   neval,ier,leniw,maxp1,lenw,last,iwork,work)
c***begin prologue  dqawo
c***date written   800101   (yymmdd)
c***revision date  830518   (yymmdd)
c***category no.  h2a2a1
c***keywords  automatic integrator, special-purpose,
c             integrand with oscillatory cos or sin factor,
c             clenshaw-curtis method, (end point) singularities,
c             extrapolation, globally adaptive
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            definite integral i=integral of f(x)*w(x) over (a,b)
c            where w(x) = cos(omega*x)
c            or w(x) = sin(omega*x),
c            hopefully satisfying following claim for accuracy
c            abs(i-result).le.max(epsabs,epsrel*abs(i)).
c***description
c
c        computation of oscillatory integrals
c        standard fortran subroutine
c        double precision version
c
c        parameters
c         on entry
c            f      - double precision
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      - double precision
c                     lower limit of integration
c
c            b      - double precision
c                     upper 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 will
c                     end with ier = 6.
c
c            epsabs - double precision
c                     absolute accuracy requested
c            epsrel - double precision
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 - 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                     ier = 1 maximum number of subdivisions allowed
c                             (= leniw/2) has been achieved. 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   - double precision
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  dqawoe,xerror
c***end prologue  dqawo
c
       double precision a,abserr,b,epsabs,epsrel,f,omega,result,work
       integer ier,integr,iwork,last,limit,lenw,leniw,lvl,l1,l2,l3,l4,
     *  maxp1,momcom,neval
c
       dimension iwork(leniw),work(lenw)
c
       external f
c
c         check validity of leniw, maxp1 and lenw.
c
c***first executable statement  dqawo
      ier = 6
      neval = 0
      last = 0
      result = 0.0d+00
      abserr = 0.0d+00
      if(leniw.lt.2.or.maxp1.lt.1.or.lenw.lt.(leniw*2+maxp1*25))
     *  go to 10
c
c         prepare call for dqawoe
c
      limit = leniw/2
      l1 = limit+1
      l2 = limit+l1
      l3 = limit+l2
      l4 = limit+l3
      call dqawoe(f,a,b,omega,integr,epsabs,epsrel,limit,1,maxp1,result,
     *   abserr,neval,ier,last,work(1),work(l1),work(l2),work(l3),
     *   iwork(1),iwork(l1),momcom,work(l4))
c
c         call error handler if necessary
c
      lvl = 0
10    if(ier.eq.6) lvl = 0
      if(ier.ne.0) call xerror(26habnormal return from dqawo,26,ier,lvl)
      return
      end
