      subroutine qsimp(func,a,b,s)
*
* Numerical integration via Simpson's rule. (From Numerical Recipies)
*
      implicit double precision(a-h,o-z)

      external func
      parameter(eps = 1.d-6, jmax=25) 

      ost = -1.e-30
      os =  -1.e-30
      do 11 j = 1,jmax
        call trapzd(func,a,b,st,j)
        s = (4.*st-ost)/3. 
        if ( dabs(s-os) .lt. eps*dabs(os)) return
        os = s
        ost = st 
11    continue
      stop
      end

************************************************************************

