[f,a,b,omega,integr,epsabs,epsrel,result,abserr,neval,ier,leniw,maxp1,lenw,last,iwork,work]=dqawo(f,a,b,omega,integr,epsabs,epsrel,result,abserr,neval,ier,leniw,maxp1,lenw,last,iwork,work);

function [f,a,b,omega,integr,epsabs,epsrel,result,abserr,neval,ier,leniw,maxp1,lenw,last,iwork,work]=dqawo(f,a,b,omega,integr,epsabs,epsrel,result,abserr,neval,ier,leniw,maxp1,lenw,last,iwork,work);
%***BEGIN PROLOGUE  DQAWO
%***PURPOSE  Calculate an approximation to a given definite integral
%            I= Integral of F(X)*W(X) over (A,B), where
%                   W(X) = COS(OMEGA*X)
%               or  W(X) = SIN(OMEGA*X),
%            hopefully satisfying the following claim for accuracy
%                ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
%***LIBRARY   SLATEC (QUADPACK)
%***CATEGORY  H2A2A1
%***TYPE      doubleprecision (QAWO-S, DQAWO-D)
%***KEYWORDS  AUTOMATIC INTEGRATOR, CLENSHAW-CURTIS METHOD,
%             EXTRAPOLATION, GLOBALLY ADAPTIVE,
%             INTEGRAND WITH OSCILLATORY COS OR SIN FACTOR, QUADPACK,
%             QUADRATURE, SPECIAL-PURPOSE
%***AUTHOR  Piessens, Robert
%             Applied Mathematics and Programming Division
%             K. U. Leuven
%           de Doncker, Elise
%             Applied Mathematics and Programming Division
%             K. U. Leuven
%***DESCRIPTION
%
%        Computation of oscillatory integrals
%        Standard fortran subroutine
%        doubleprecision version
%
%        PARAMETERS
%         ON ENTRY
%            F      - doubleprecision
%                     function subprogram defining the function
%                     F(X).  The actual name for F needs to be
%                     declared E X T E R N A L in the driver program.
%
%            A      - doubleprecision
%                     Lower limit of integration
%
%            B      - doubleprecision
%                     Upper limit of integration
%
%            OMEGA  - doubleprecision
%                     Parameter in the integrand weight function
%
%            INTEGR - Integer
%                     Indicates which of the weight functions is used
%                     INTEGR = 1      W(X) = COS(OMEGA*X)
%                     INTEGR = 2      W(X) = SIN(OMEGA*X)
%                     If INTEGR.NE.1.AND.INTEGR.NE.2, the routine will
%                     end with IER = 6.
%
%            EPSABS - doubleprecision
%                     Absolute accuracy requested
%            EPSREL - doubleprecision
%                     Relative accuracy requested
%                     If EPSABS.LE.0 and
%                     EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
%                     the routine will end with IER = 6.
%
%         ON RETURN
%            RESULT - doubleprecision
%                     Approximation to the integral
%
%            ABSERR - doubleprecision
%                     Estimate of the modulus of the absolute error,
%                     which should equal or exceed ABS(I-RESULT)
%
%            NEVAL  - Integer
%                     Number of integrand evaluations
%
%            IER    - Integer
%                     IER = 0 Normal and reliable termination of the
%                             routine. It is assumed that the requested
%                             accuracy has been achieved.
%                   - IER.GT.0 Abnormal termination of the routine.
%                             The estimates for integral and error are
%                             less reliable. It is assumed that the
%                             requested accuracy has not been achieved.
%            ERROR MESSAGES
%                     IER = 1 Maximum number of subdivisions allowed
%                             has been achieved (= LENIW/2). One can
%                             allow more subdivisions by increasing the
%                             value of LENIW (and taking the according
%                             dimension adjustments into account).
%                             However, if this yields no improvement it
%                             is advised to analyze the integrand in
%                             order to determine the integration
%                             difficulties. If the position of a local
%                             difficulty can be determined (e.g.
%                             SINGULARITY, DISCONTINUITY within the
%                             interval) one will probably gain from
%                             splitting up the interval at this point
%                             and calling the integrator on the
%                             subranges. If possible, an appropriate
%                             special-purpose integrator should be used
%                             which is designed for handling the type of
%                             difficulty involved.
%                         = 2 The occurrence of roundoff error is
%                             detected, which prevents the requested
%                             tolerance from being achieved.
%                             The error may be under-estimated.
%                         = 3 Extremely bad integrand behaviour occurs
%                             at some interior points of the
%                             integration interval.
%                         = 4 The algorithm does not converge.
%                             Roundoff error is detected in the
%                             extrapolation table. It is presumed that
%                             the requested tolerance cannot be achieved
%                             due to roundoff in the extrapolation
%                             table, and that the returned result is
%                             the best which can be obtained.
%                         = 5 The integral is probably divergent, or
%                             slowly convergent. It must be noted that
%                             divergence can occur with any other value
%                             of IER.
%                         = 6 The input is invalid, because
%                             (EPSABS.LE.0 and
%                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28))
%                             or (INTEGR.NE.1 AND INTEGR.NE.2),
%                             or LENIW.LT.2 OR MAXP1.LT.1 or
%                             LENW.LT.LENIW*2+MAXP1*25.
%                             RESULT, ABSERR, NEVAL, LAST are set to
%                             zero. Except when LENIW, MAXP1 or LENW are
%                             invalid, WORK(LIMIT*2+1), WORK(LIMIT*3+1),
%                             IWORK(1), IWORK(LIMIT+1) are set to zero,
%                             WORK(1) is set to A and WORK(LIMIT+1) to
%                             B.
%
%         DIMENSIONING PARAMETERS
%            LENIW  - Integer
%                     Dimensioning parameter for IWORK.
%                     LENIW/2 equals the maximum number of subintervals
%                     allowed in the partition of the given integration
%                     interval (A,B), LENIW.GE.2.
%                     If LENIW.LT.2, the routine will end with IER = 6.
%
%            MAXP1  - Integer
%                     Gives an upper bound on the number of Chebyshev
%                     moments which can be stored, i.e. for the
%                     intervals of lengths ABS(B-A)*2**(-L),
%                     L=0,1, ..., MAXP1-2, MAXP1.GE.1
%                     If MAXP1.LT.1, the routine will end with IER = 6.
%
%            LENW   - Integer
%                     Dimensioning parameter for WORK
%                     LENW must be at least LENIW*2+MAXP1*25.
%                     If LENW.LT.(LENIW*2+MAXP1*25), the routine will
%                     end with IER = 6.
%
%            LAST   - Integer
%                     On return, LAST equals the number of subintervals
%                     produced in the subdivision process, which
%                     determines the number of significant elements
%                     actually in the WORK ARRAYS.
%
%         WORK ARRAYS
%            IWORK  - Integer
%                     Vector of dimension at least LENIW
%                     on return, the first K elements of which contain
%                     pointers to the error estimates over the
%                     subintervals, such that WORK(LIMIT*3+IWORK(1)), ..
%                     WORK(LIMIT*3+IWORK(K)) form a decreasing
%                     sequence, with LIMIT = LENW/2 , and K = LAST
%                     if LAST.LE.(LIMIT/2+2), and K = LIMIT+1-LAST
%                     otherwise.
%                     Furthermore, IWORK(LIMIT+1), ..., IWORK(LIMIT+
%                     LAST) indicate the subdivision levels of the
%                     subintervals, such that IWORK(LIMIT+I) = L means
%                     that the subinterval numbered I is of length
%                     ABS(B-A)*2**(1-L).
%
%            WORK   - doubleprecision
%                     Vector of dimension at least LENW
%                     On return
%                     WORK(1), ..., WORK(LAST) contain the left
%                      end points of the subintervals in the
%                      partition of (A,B),
%                     WORK(LIMIT+1), ..., WORK(LIMIT+LAST) contain
%                      the right end points,
%                     WORK(LIMIT*2+1), ..., WORK(LIMIT*2+LAST) contain
%                      the integral approximations over the
%                      subintervals,
%                     WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST)
%                      contain the error estimates.
%                     WORK(LIMIT*4+1), ..., WORK(LIMIT*4+MAXP1*25)
%                      Provide space for storing the Chebyshev moments.
%                     Note that LIMIT = LENW/2.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  DQAWOE, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   800101  DATE WRITTEN
%   890831  Modified array declarations.  (WRB)
%   890831  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
%***end PROLOGUE  DQAWO
%
persistent l1 l2 l3 l4 limit lvl momcom ; 

if isempty(limit), limit=0; end;
if isempty(lvl), lvl=0; end;
if isempty(l1), l1=0; end;
if isempty(l2), l2=0; end;
if isempty(l3), l3=0; end;
if isempty(l4), l4=0; end;
if isempty(momcom), momcom=0; end;
%
iwork_shape=size(iwork);iwork=reshape(iwork,1,[]);
work_shape=size(work);work=reshape(work,1,[]);
%
%
%         CHECK VALIDITY OF LENIW, MAXP1 AND LENW.
%
%***FIRST EXECUTABLE STATEMENT  DQAWO
ier = 6;
neval = 0;
last = 0;
result = 0.0d+00;
abserr = 0.0d+00;
if( leniw>=2 && maxp1>=1 && lenw>=(leniw.*2+maxp1.*25) )
%
%         PREPARE CALL FOR DQAWOE
%
limit = fix(fix(leniw./2));
l1 = fix(limit + 1);
l2 = fix(limit + l1);
l3 = fix(limit + l2);
l4 = fix(limit + l3);
[f,a,b,omega,integr,epsabs,epsrel,limit,dumvar9,maxp1,result,abserr,neval,ier,last,dumvar16,dumvar17,dumvar18,dumvar19,idumvar16,idumvar17,momcom,dumvar23]=dqawoe(f,a,b,omega,integr,epsabs,epsrel,limit,1,maxp1,result,abserr,neval,ier,last,work(sub2ind(size(work),max(1,1)):end),work(sub2ind(size(work),max(l1,1)):end),work(sub2ind(size(work),max(l2,1)):end),work(sub2ind(size(work),max(l3,1)):end),iwork(sub2ind(size(work),max(1,1)):end),iwork(sub2ind(size(work),max(l1,1)):end),momcom,work(l4:end));   dumvar16i=find((work(sub2ind(size(work),max(1,1)):end))~=(dumvar16));dumvar17i=find((work(sub2ind(size(work),max(l1,1)):end))~=(dumvar17));dumvar18i=find((work(sub2ind(size(work),max(l2,1)):end))~=(dumvar18));dumvar19i=find((work(sub2ind(size(work),max(l3,1)):end))~=(dumvar19));dumvar23i=find((work(l4:end))~=(dumvar23));   work(1-1+dumvar16i)=dumvar16(dumvar16i); work(l1-1+dumvar17i)=dumvar17(dumvar17i); work(l2-1+dumvar18i)=dumvar18(dumvar18i); work(l3-1+dumvar19i)=dumvar19(dumvar19i); work(l4-1+dumvar23i)=dumvar23(dumvar23i); 
%
%         CALL ERROR HANDLER IF NECESSARY
%
lvl = 0;
end;
if( ier==6 )
lvl = 0;
end;
if( ier~=0 )
[dumvar1,dumvar2,dumvar3,ier,lvl]=xermsg('SLATEC','DQAWO','ABNORMAL RETURN',ier,lvl);
end;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
end
%DECK DQAWSE