subroutine qc25f(f,a,b,omega,integr,nrmom,maxp1,ksave,result,
     *   abserr,neval,resabs,resasc,momcom,chebmo)
c***begin prologue  qc25f
c***date written   810101   (yymmdd)
c***revision date  211011   (yymmdd)
c***category no.  h2a2a2
c***keywords  integration rules for functions with cos or sin
c             factor, clenshaw-curtis, gauss-kronrod
c***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
c           de doncker,elise,appl. math. & progr. div. - k.u.leuven
c***purpose  to compute the integral i=integral of f(x) over (a,b)
c            where w(x) = cos(omega*x) or (wx)=sin(omega*x)
c            and to compute j=integral of abs(f) over (a,b). for small
c            value of omega or small intervals (a,b) 15-point gauss-
c            kronrod rule used. otherwise generalized clenshaw-curtis us
c***description
c
c        integration rules for functions with cos or sin factor
c        standard fortran subroutine
c        real version
c
c        parameters
c         on entry
c           f      - real
c                    function subprogram defining the integrand
c                    function f(x). the actual name for f needs to
c                    be declared e x t e r n a l in the calling 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 weight function
c
c           integr - integer
c                    indicates which weight function is to be used
c                       integr = 1   w(x) = cos(omega*x)
c                       integr = 2   w(x) = sin(omega*x)
c
c           nrmom  - integer
c                    the length of interval (a,b) is equal to the length
c                    of the original integration interval divided by
c                    2**nrmom (we suppose that the routine is used in an
c                    adaptive integration process, otherwise set
c                    nrmom = 0). nrmom must be zero at the first call.
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(bb-aa)*2**(-l),
c                    l = 0,1,2, ..., maxp1-2.
c
c           ksave  - integer
c                    key which is one when the moments for the
c                    current interval have been computed
c
c         on return
c           result - real
c                    approximation to the integral i
c
c           abserr - real
c                    estimate of the modulus of the absolute
c                    error, which should equal or exceed abs(i-result)
c
c           neval  - integer
c                    number of integrand evaluations
c
c           resabs - real
c                    approximation to the integral j
c
c           resasc - real
c                    approximation to the integral of abs(f-i/(b-a))
c
c         on entry and return
c           momcom - integer
c                    for each interval length we need to compute the
c                    chebyshev moments. momcom counts the number of
c                    intervals for which these moments have already been
c                    computed. if nrmom.lt.momcom or ksave = 1, the
c                    chebyshev moments for the interval (a,b) have
c                    already been computed and stored, otherwise we
c                    compute them and we increase momcom.
c
c           chebmo - real
c                    array of dimension at least (maxp1,25) containing
c                    the modified chebyshev moments for the first momcom
c                    momcom interval lengths
c
c***references  (none)
c***routines called  qcheb,qk15w,qwgtf,r1mach,sgtsl
c***end prologue  qc25f
c
      real a,abserr,ac,an,an2,as,asap,ass,b,centr,chebmo,
     *  cheb12,cheb24,conc,cons,cospar,d,qwgtf,
     *  d1,r1mach,d2,estc,ests,f,fval,hlgth,oflow,omega,parint,par2,
     *  par22,p2,p3,p4,resabs,resasc,resc12,resc24,ress12,ress24,
     *  result,sinpar,v,x
      integer i,iers,integr,isym,j,k,ksave,m,maxp1,momcom,neval,
     *  noequ,noeq1,nrmom
c
      dimension chebmo(maxp1,25),cheb12(13),cheb24(25),d(25),d1(25),
     *  d2(25),fval(25),v(28),x(11)
c
      external f,qwgtf
c
c           the vector x contains the values cos(k*pi/24)
c           k = 1, ...,11, to be used for the chebyshev expansion of f
c
      data x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8),x(9),
     *  x(10),x(11)/
     *     0.9914448613738104e+00,     0.9659258262890683e+00,
     *     0.9238795325112868e+00,     0.8660254037844386e+00,
     *     0.7933533402912352e+00,     0.7071067811865475e+00,
     *     0.6087614290087206e+00,     0.5000000000000000e+00,
     *     0.3826834323650898e+00,     0.2588190451025208e+00,
     *     0.1305261922200516e+00/
c
c           list of major variables
c           -----------------------
c
c           centr  - mid point of the integration interval
c           hlgth  - half-length of the integration interval
c           fval   - value of the function f at the points
c                    (b-a)*0.5*cos(k*pi/12) + (b+a)*0.5,
c                    k = 0, ..., 24
c           cheb12 - coefficients of the chebyshev series expansion
c                    of degree 12, for the function f, in the
c                    interval (a,b)
c           cheb24 - coefficients of the chebyshev series expansion
c                    of degree 24, for the function f, in the
c                    interval (a,b)
c           resc12 - approximation to the integral of
c                    cos(0.5*(b-a)*omega*x)*f(0.5*(b-a)*x+0.5*(b+a))
c                    over (-1,+1), using the chebyshev series
c                    expansion of degree 12
c           resc24 - approximation to the same integral, using the
c                    chebyshev series expansion of degree 24
c           ress12 - the analogue of resc12 for the sine
c           ress24 - the analogue of resc24 for the sine
c
c
c           machine dependent constant
c           --------------------------
c
c           oflow is the largest positive magnitude.
c
c***first executable statement  qc25f
      oflow = r1mach(2)
c
      centr = 0.5e+00*(b+a)
      hlgth = 0.5e+00*(b-a)
      parint = omega*hlgth
c
c           compute the integral using the 15-point gauss-kronrod
c           formula if the value of the parameter in the integrand
c           is small.
c
      if(abs(parint).gt.0.2e+01) go to 10
      call qk15w(f,qwgtf,omega,p2,p3,p4,integr,a,b,result,
     *  abserr,resabs,resasc)
      neval = 15
      go to 170
c
c           compute the integral using the generalized clenshaw-
c           curtis method.
c
   10 conc = hlgth*cos(centr*omega)
      cons = hlgth*sin(centr*omega)
      resasc = oflow
      neval = 25
c
c           check whether the chebyshev moments for this interval
c           have already been computed.
c
      if(nrmom.lt.momcom.or.ksave.eq.1) go to 120
c
c           compute a new set of chebyshev moments.
c
      m = momcom+1
      par2 = parint*parint
      par22 = par2+0.2e+01
      sinpar = sin(parint)
      cospar = cos(parint)
c
c           compute the chebyshev moments with respect to cosine.
c
      v(1) = 0.2e+01*sinpar/parint
      v(2) = (0.8e+01*cospar+(par2+par2-0.8e+01)*sinpar/
     *  parint)/par2
      v(3) = (0.32e+02*(par2-0.12e+02)*cospar+(0.2e+01*
     *  ((par2-0.80e+02)*par2+0.192e+03)*sinpar)/
     *  parint)/(par2*par2)
      ac = 0.8e+01*cospar
      as = 0.24e+02*parint*sinpar
      if(abs(parint).gt.0.24e+02) go to 30
c
c           compute the chebyshev moments as the
c           solutions of a boundary value problem with 1
c           initial value (v(3)) and 1 end value (computed
c           using an asymptotic formula).
c
      noequ = 25
      noeq1 = noequ-1
      an = 0.6e+01
      do 20 k = 1,noeq1
        an2 = an*an
        d(k) = -0.2e+01*(an2-0.4e+01)*(par22-an2-an2)
        d2(k) = (an-0.1e+01)*(an-0.2e+01)*par2
        d1(k+1) = (an+0.3e+01)*(an+0.4e+01)*par2
        v(k+3) = as-(an2-0.4e+01)*ac
        an = an+0.2e+01
   20 continue
      an2 = an*an
      d(noequ) = -0.2e+01*(an2-0.4e+01)*(par22-an2-an2)
      v(noequ+3) = as-(an2-0.4e+01)*ac
      v(4) = v(4)-0.56e+02*par2*v(3)
      ass = parint*sinpar
      asap = (((((0.210e+03*par2-0.1e+01)*cospar-(0.105e+03*par2
     *  -0.63e+02)*ass)/an2-(0.1e+01-0.15e+02*par2)*cospar
     *  +0.15e+02*ass)/an2-cospar+0.3e+01*ass)/an2-cospar)/an2
      v(noequ+3) = v(noequ+3)-0.2e+01*asap*par2*(an-0.1e+01)*
     *   (an-0.2e+01)
c
c           solve the tridiagonal system by means of gaussian
c           elimination with partial pivoting.
c
      call sgtsl(noequ,d1,d,d2,v(4),iers)
      go to 50
c
c           compute the chebyshev moments by means of forward
c           recursion.
c
   30 an = 0.4e+01
      do 40 i = 4,13
        an2 = an*an
        v(i) = ((an2-0.4e+01)*(0.2e+01*(par22-an2-an2)*v(i-1)-ac)
     *  +as-par2*(an+0.1e+01)*(an+0.2e+01)*v(i-2))/
     *  (par2*(an-0.1e+01)*(an-0.2e+01))
        an = an+0.2e+01
   40 continue
   50 do 60 j = 1,13
        chebmo(m,2*j-1) = v(j)
   60 continue
c
c           compute the chebyshev moments with respect to sine.
c
      v(1) = 0.2e+01*(sinpar-parint*cospar)/par2
      v(2) = (0.18e+02-0.48e+02/par2)*sinpar/par2
     *  +(-0.2e+01+0.48e+02/par2)*cospar/parint
      ac = -0.24e+02*parint*cospar
      as = -0.8e+01*sinpar
      if(abs(parint).gt.0.24e+02) go to 80
c
c           compute the chebyshev moments as the
c           solutions of a boundary value problem with 1
c           initial value (v(2)) and 1 end value (computed
c           using an asymptotic formula).
c
      an = 0.5e+01
      do 70 k = 1,noeq1
        an2 = an*an
        d(k) = -0.2e+01*(an2-0.4e+01)*(par22-an2-an2)
        d2(k) = (an-0.1e+01)*(an-0.2e+01)*par2
        d1(k+1) = (an+0.3e+01)*(an+0.4e+01)*par2
        v(k+2) = ac+(an2-0.4e+01)*as
        an = an+0.2e+01
   70 continue
      an2 = an*an
      d(noequ) = -0.2e+01*(an2-0.4e+01)*(par22-an2-an2)
      v(noequ+2) = ac+(an2-0.4e+01)*as
      v(3) = v(3)-0.42e+02*par2*v(2)
      ass = parint*cospar
      asap = (((((0.105e+03*par2-0.63e+02)*ass+(0.210e+03*par2
     *  -0.1e+01)*sinpar)/an2+(0.15e+02*par2-0.1e+01)*sinpar-
     *  0.15e+02*ass)/an2-0.3e+01*ass-sinpar)/an2-sinpar)/an2
      v(noequ+2) = v(noequ+2)-0.2e+01*asap*par2*(an-0.1e+01)
     *  *(an-0.2e+01)
c
c           solve the tridiagonal system by means of gaussian
c           elimination with partial pivoting.
c
      call sgtsl(noequ,d1,d,d2,v(3),iers)
      go to 100
c
c           compute the chebyshev moments by means of
c           forward recursion.
c
   80 an = 0.3e+01
      do 90 i = 3,12
        an2 = an*an
        v(i) = ((an2-0.4e+01)*(0.2e+01*(par22-an2-an2)*v(i-1)+as)
     *  +ac-par2*(an+0.1e+01)*(an+0.2e+01)*v(i-2))
     *  /(par2*(an-0.1e+01)*(an-0.2e+01))
        an = an+0.2e+01
   90 continue
  100 do 110 j = 1,12
        chebmo(m,2*j) = v(j)
  110 continue
  120 if (nrmom.lt.momcom) m = nrmom+1
       if (momcom.lt.maxp1-1.and.nrmom.ge.momcom) momcom = momcom+1
c
c           compute the coefficients of the chebyshev expansions
c           of degrees 12 and 24 of the function f.
c
      fval(1) = 0.5e+00*f(centr+hlgth)
      fval(13) = f(centr)
      fval(25) = 0.5e+00*f(centr-hlgth)
      do 130 i = 2,12
        isym = 26-i
        fval(i) = f(hlgth*x(i-1)+centr)
        fval(isym) = f(centr-hlgth*x(i-1))
  130 continue
      call qcheb(x,fval,cheb12,cheb24)
c
c           compute the integral and error estimates.
c
      resc12 = cheb12(13)*chebmo(m,13)
      ress12 = 0.0e+00
      k = 11
      do 140 j = 1,6
        resc12 = resc12+cheb12(k)*chebmo(m,k)
        ress12 = ress12+cheb12(k+1)*chebmo(m,k+1)
        k = k-2
  140 continue
      resc24 = cheb24(25)*chebmo(m,25)
      ress24 = 0.0e+00
      resabs = abs(cheb24(25))
      k = 23
      do 150 j = 1,12
        resc24 = resc24+cheb24(k)*chebmo(m,k)
        ress24 = ress24+cheb24(k+1)*chebmo(m,k+1)
        resabs = resabs+abs(cheb24(k))+abs(cheb24(k+1))
        k = k-2
  150 continue
      estc = abs(resc24-resc12)
      ests = abs(ress24-ress12)
      resabs = resabs*abs(hlgth)
      if(integr.eq.2) go to 160
      result = conc*resc24-cons*ress24
      abserr = abs(conc*estc)+abs(cons*ests)
      go to 170
  160 result = conc*ress24+cons*resc24
      abserr = abs(conc*ests)+abs(cons*estc)
  170 return
      end