real function ccquad(f,a,b,tolerr,limit,esterr,used,
     *  csxfrm)

      implicit none
c
c  input arguments-
c
      real f,a,b,tolerr
      integer limit
c  output arguments-
      real esterr,csxfrm(limit)
      integer used
c
c  using clenshaw-curtis quadrature, this function sub-
c  program attempts to integrate the function f from a to b
c  to at least the requested relative accuracy tolerr, while
c  using no more than limit function evaluations.  if this
c  can be done, ccquad returns the value of the integral,
c  esterr returns an estimate of the absolute error actually
c  committed, used returns the number of function values
c  actually used, and csxfrm(1),...,csxfrm(used) contains
c  n=used-1 times the discrete cosine transform, as usually
c  defined, of the integrand in the interval.  if the
c  requested accuracy cannot be obtained with the number of
c  function evaluations permitted, the last (and presumably
c  best) answer obtained is returned.
c
      real pi,rt3,centre,width,shift,fund,angle,c,s
      real oldint,newint
      real t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12
c  insert the following statement to trace program flow.
c     real sclint,sclerr
      integer n,n2,n3,n_less_1,n_less_3,max,m_max,j,step
      integer l(8),l1,l2,l3,l4,l5,l6,l7,l8
      integer j1,j2,j3,j4,j5,j6,j7,j8,j_rev
      equivalence (l(1),l1),(l(2),l2),(l(3),l3),(l(4),l4),
     *  (l(5),l5),(l(6),l6),(l(7),l7),(l(8),l8),(j8,j_rev)
      data pi,rt3/3.141592653589e0, 1.732050807568e0 /
      data m_max / 8 /
c
c  initialization.
c
      centre=(a+b)*.5e0
      width=(b-a)*.5e0
      max=min0(limit,2*3**(m_max+1))
      do j=1,m_max
        l(j)=1
      end do
c
c  cosine transform.
c  compute double the cosine transform with n=6.
c
      n=6
c
c  sample the function.
c
      csxfrm(1)=f(a)
      csxfrm(7)=f(b)
      shift=width*rt3*.5e0
      csxfrm(2)=f(centre-shift)
      csxfrm(6)=f(centre+shift)
      shift=width*.5e0
      csxfrm(3)=f(centre-shift)
      csxfrm(5)=f(centre+shift)
      csxfrm(4)=f(centre)
c  evaluate the factored n=6 cosine transform.
      t1=csxfrm(1)+csxfrm(7)
      t2=csxfrm(1)-csxfrm(7)
      t3=2.e0*csxfrm(4)
      t4=csxfrm(2)+csxfrm(6)
      t5=(csxfrm(2)-csxfrm(6))*rt3
      t6=csxfrm(3)+csxfrm(5)
      t7=csxfrm(3)-csxfrm(5)
      t8=t1+2.e0*t6
      t9=2.e0*t4+t3
      t10=t2+t7
      t11=t1-t6
      t12=t4-t3
      csxfrm(1)=t8+t9
      csxfrm(2)=t10+t5
      csxfrm(3)=t11+t12
      csxfrm(4)=t2-2.e0*t7
      csxfrm(5)=t11-t12
      csxfrm(6)=t10-t5
      csxfrm(7)=t8-t9
      used=7
c  go to integral computation, but first compute integral for n=2.
      go to 200
c
c  compute refined approximation.
c  sample function at intermediate points in digit reversed
c  order.  as the sequence is generated, compute the first
c  (radix four transform) pass of the fast fourier transform.
c
100   do j=2,m_max
        l(j-1)=l(j)
      end do
      l(m_max)=3*l(m_max-1)
      j=used
      fund=pi/float(3*n)

      do j1=1,l1,1
        do j2=j1,l2,l1
          do j3=j2,l3,l2
            do j4=j3,l4,l3
              do j5=j4,l5,l4
                do j6=j5,l6,l5
                  do j7=j6,l7,l6
                    do j8=j7,l8,l7
                      angle=fund*float(3*j_rev-2)
                      shift=width*cos(angle)
                      t1=f(centre-shift)
                      t3=f(centre+shift)
                      shift=width*sin(angle)
                      t2=f(centre+shift)
                      t4=f(centre-shift)
                      t5=t1+t3
                      t6=t2+t4
                      csxfrm(j+1)=t5+t6
                      csxfrm(j+2)=t1-t3
                      csxfrm(j+3)=t5-t6
                      csxfrm(j+4)=t2-t4
                      j=j+4
                    end do
                  end do
                end do
              end do
            end do
          end do
        end do
      end do
c
c  do radix 3 passes of fast fourier transform.
      n2=2*n
      step=4
150   j1=used+step
      j2=used+2*step
      call r3pass(n2,step,n2-2*step,csxfrm(used+1),
     *  csxfrm(j1+1),csxfrm(j2+1))
      step=3*step
      if (step .lt. n) go to 150
c
c  combine results.
c
c  first do j=0 and j=n.
c
      t1=csxfrm(1)
      t2=csxfrm(used+1)
      csxfrm(1)=t1+2.e0*t2
      csxfrm(used+1)=t1-t2
      t1=csxfrm(n+1)
      t2=csxfrm(n2+2)
      csxfrm(n+1)=t1+t2
      csxfrm(n2+2)=t1-2.e0*t2
c  now do remaining values of j.
      n3=3*n
      n_less_1=n-1
      do j=1,n_less_1
        j1=n+j
        j2=n3-j
        angle=fund*float(j)
        c=cos(angle)
        s=sin(angle)
        t1=c*csxfrm(j1+2)-s*csxfrm(j2+2)
        t2=(s*csxfrm(j1+2)+c*csxfrm(j2+2))*rt3
        csxfrm(j1+2)=csxfrm(j+1)-t1-t2
        csxfrm(j2+2)=csxfrm(j+1)-t1+t2
        csxfrm(j+1)=csxfrm(j+1)+2.e0*t1
      end do
c  now unscramble.
      t1=csxfrm(n2+1)
      t2=csxfrm(n2+2)
      do j=1,n_less_1
        j1=used+j
        j2=n2+j
        csxfrm(j2)=csxfrm(j1)
        csxfrm(j1)=csxfrm(j2+2)
      end do
      csxfrm(n3)=t1
      csxfrm(n3+1)=t2
      n=n3
      used=n+1
c
c  go to integral computation.
c
      go to 210
c
c  integral evaluation.
c  integral estimates are not scaled by width*n/2
c  until function ccquad returns.
c
c  when n=6, evaluate integral for n=2.
c
200   oldint=(t1+2.e0*t3)/3.e0
c
c  evaluate new estimate of integral.
c
210   n_less_3=n-3
      newint=.5e0*csxfrm(used)/float(1-n**2)
      do j=1,n_less_3,2
        j_rev=n-j
        newint=newint+csxfrm(j_rev)/float(j_rev*(2-j_rev))
      end do
      newint=newint+.5e0*csxfrm(1)
c
c  test if done.
c  test if estimated error adequate.
      esterr=abs(oldint*3.e0-newint)
c
c  insert the following four statements to trace program flow.
c     sclint=width*newint/float(n/2)
c     sclerr=width*(oldint*3.e0-newint)/float(n/2)
c     write(6,900) n,sclint,sclerr
c900  format ( 3h n=,i5,23h integral estimated as ,e15.8,
c    *  7h error ,e15.8 )
      if ( abs(newint)*tolerr .ge. esterr ) go to 400
c  if estimated error too large, refine sampling if permitted.
      oldint=newint
      if (3*n+1 .le. max ) go to 100
c  if refinement not permitted, or if estimated error
c  satisfactory, rescale answers and return.
c  insert the following two statements to trace program flow.
c     write (6,910)
c 910 format ( 25h refinement not permitted )
c
400   ccquad=width*newint/float(n/2)
      esterr=width*esterr/float(n/2)
      return
      end
      subroutine r3pass(n2,m,length,x0,x1,x2)

      implicit none

c  radix 3 pass for fast fourier transform of real sequence
c  of length n2.
      integer n2,m,length
      real x0(length),x1(length),x2(length)
c  the notation of references 2 and 3 is used in this
c  subroutine.
c  m is the length of the transform already accomplished.
c  i.e., the number of distinct values of the frequency index
c  c hat of these transforms, and the spacing of the
c  sequences to be transformed.  explicit use is made of the
c  fact that m is even and not less than 4.
      integer half_m,m3,k,k0,k1,j,j0,j1
      real twopi,hafrt3,rsum,rdiff,rsum2,isum,idiff,idiff2
      real fund,angle,c1,s1,c2,s2,r0,r1,r2,i0,i1,i2
      data twopi, hafrt3 / 6.283185307e0, .866025403e0 /
      half_m=(m-1)/2
      m3=m*3
      fund=twopi/float(m3)
c  do all transforms for c hat = 0, i.e., twiddle factor unity.
      do k=1,n2,m3
        rsum=(x1(k)+x2(k))
        rdiff=(x1(k)-x2(k))*hafrt3
        x1(k)=x0(k)-rsum*.5e0
        x2(k)=rdiff
        x0(k)=x0(k)+rsum
      end do
c  do all transforms for c hat = cap c/2, i.e., twiddle factor.
c  e(b/6)
      j=m/2+1
      do k=j,n2,m3
        rsum=(x1(k)+x2(k))*hafrt3
        rdiff=(x1(k)-x2(k))
        x1(k)=x0(k)-rdiff
        x2(k)=rsum
        x0(k)=x0(k)+rdiff*.5e0
      end do
c  do all transforms for remaining values of c hat.  observe
c  that c hat and cap c-c hat must be paired.
c  choose a frequency index.
      do j=1,half_m
        j0=j+1
        j1=m-j+1
c  compute the twiddle factor.
        angle=fund*float(j)
        c1=cos(angle)
        s1=sin(angle)
        c2=c1**2-s1**2
        s2=2.e0*s1*c1
c  choose the replication.
        do k0=j0,n2,m3
          k1=k0-j0+j1
c  obtain twiddled values.
          r0=x0(k0)
          i0=x0(k1)
          r1=c1*x1(k0)-s1*x1(k1)
          i1=s1*x1(k0)+c1*x1(k1)
          r2=c2*x2(k0)-s2*x2(k1)
          i2=s2*x2(k0)+c2*x2(k1)
c  compute transforms and return in place.
          rsum=r1+r2
          rdiff=(r1-r2)*hafrt3
          rsum2=r0-.5e0*rsum
          isum=i1+i2
          idiff=(i1-i2)*hafrt3
          idiff2=i0-.5e0*isum
          x0(k0)=r0+rsum
          x0(k1)=rsum2+idiff
          x1(k0)=rsum2-idiff
          x1(k1)=rdiff+idiff2
          x2(k0)=rdiff-idiff2
          x2(k1)=i0+isum
        end do
      end do

      return
      end