SUBROUTINE qromb(func,a,b,ss,q) ! one parameter
      INTEGER JMAX,JMAXP,K,KM
      REAL*8 a,b,func,ss,q,EPS
      EXTERNAL func
      PARAMETER (EPS=1.d-6, JMAX=20, JMAXP=JMAX+1, K=5, KM=K-1)
CU    USES polint,trapzd
      INTEGER j
      REAL*8 dss,h(JMAXP),s(JMAXP)
      h(1)=1d0
      do 11 j=1,JMAX
        call trapzd(func,a,b,s(j),j,q)
        if (j.ge.K) then
          call polint(h(j-KM),s(j-KM),K,0d0,ss,dss)
          if (dabs(dss).le.EPS*dabs(ss)) return
        endif
        s(j+1)=s(j)
        h(j+1)=0.25d0*h(j)
11    continue
      pause 'too many steps in qromb'
      END
C  (C) Copr. 1986-92 Numerical Recipes Software D041&0(9p#3.
      SUBROUTINE trapzd(func,a,b,s,n,q) ! one paramete
      INTEGER n
      REAL*8 a,b,s,func,q
      EXTERNAL func
      INTEGER it,j
      REAL*8 del,sum,tnm,x
      if (n.eq.1) then
        s=0.5d0*(b-a)*(func(a,q)+func(b,q))
      else
        it=2**(n-2)
        tnm=dble(it)
        del=(b-a)/tnm
        x=a+0.5d0*del
        sum=0d0
        do 11 j=1,it
          sum=sum+func(x,q)
          x=x+del
11      continue
        s=0.5d0*(s+(b-a)*sum/tnm)
      endif
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software D041&0(9p#3.
      SUBROUTINE polint(xa,ya,n,x,y,dy)
      INTEGER n,NMAX
      REAL*8 dy,x,y,xa(n),ya(n)
      PARAMETER (NMAX=50)
      INTEGER i,m,ns
      REAL*8 den,dif,dift,ho,hp,w,c(NMAX),d(NMAX)
      ns=1
      dif=dabs(x-xa(1))
      do 11 i=1,n
        dift=dabs(x-xa(i))
        if (dift.lt.dif) then
          ns=i
          dif=dift
        endif
        c(i)=ya(i)
        d(i)=ya(i)
11    continue
      y=ya(ns)
      ns=ns-1
      do 13 m=1,n-1
        do 12 i=1,n-m
          ho=xa(i)-x
          hp=xa(i+m)-x
          w=c(i+1)-d(i)
          den=ho-hp
          if(den.eq.0d0)pause 'failure in polint'
          den=w/den
          d(i)=hp*den
          c(i)=ho*den
12      continue
        if (2*ns.lt.n-m)then
          dy=c(ns+1)
        else
          dy=d(ns)
          ns=ns-1
        endif
        y=y+dy
13    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software D041&0(9p#3.