FUNCTION qtrap(func,a,b, eps_in)
	USE nrtype; USE nrutil, ONLY : nrerror
	USE nr, ONLY : trapzd
	IMPLICIT NONE
	REAL(dp), INTENT(IN) :: a,b
	REAL(dp) :: qtrap, eps
  real(dp), optional :: eps_in
	INTERFACE
		FUNCTION func(x)
		USE nrtype
		REAL(dp), DIMENSION(:), INTENT(IN) :: x
		REAL(dp), DIMENSION(size(x)) :: func
		END FUNCTION func
	END INTERFACE
	INTEGER(I4B), PARAMETER :: JMAX=20
	REAL(dp), PARAMETER :: EPS_std=1.0e-6_dp
	REAL(dp) :: olds
	INTEGER(I4B) :: j
  if (present(eps_in)) then
    eps = eps_in
  else
    eps = eps_std
  endif
	olds=0.0
	do j=1,JMAX
		call trapzd(func,a,b,qtrap,j)
		if (j > 5) then
			if (abs(qtrap-olds) <= EPS*abs(olds) .or. &
				(qtrap == 0.0 .and. olds == 0.0)) RETURN
		end if
		olds=qtrap
	end do
	call nrerror('qtrap: too many steps')
	END FUNCTION qtrap