module numinteg
	integer :: numcalls
	interface
		subroutine integrate(f, a, b, acc, eps, integral, err)
			real*8, external :: f
			real*8, intent(in) :: a, b, acc, eps
			real*8, intent(out) :: integral, err
		end subroutine integrate
	
		recursive subroutine adapt(f, a, b, acc, eps, f2, f3, integral, err)
			real*8, external :: f
			real*8, intent(in) :: a, b, acc, eps, f2, f3
			real*8, intent(out) :: integral, err
		end subroutine adapt
	end interface
end module

subroutine integrate(f, a, b, acc, eps, integral, err)
	use numinteg, only : adapt
	implicit none
	real*8, external :: f
	real*8, intent(in) :: a, b, acc, eps
	real*8, intent(out) :: integral, err
	real*8 :: h, f2, f3
	
	! Evaluate the first function points to start the recursive routine:
	h = b-a
	f2 = f(a+h*1.d0/3.d0)
	f3 = f(a+h*2.d0/3.d0)
	! Start the recursive routine:
	call adapt(f, a, b, acc, eps, f2, f3, integral, err)
end subroutine

recursive subroutine adapt(f, a, b, acc, eps, f2, f3, integral, err)
	implicit none
	real*8, external :: f
	real*8, intent(in) :: a, b, acc, eps, f2, f3
	real*8, intent(out) :: integral, err
	real*8 :: tol, int1, int2, err1, err2
	real*8 :: h, f1, f4, integral2
	
	! Evaluate the function-points:
	h = b-a
!	print *, h
!	print *, f2,f3
!	read(*,*)
	f1 = f(a+h*1.d0/6.d0)
	f4 = f(a+h*5.d0/6.d0)
	
	! Calculate the integrals:
	integral2 = (h/2.d0)*(f2+f3)
	integral = (h/6.d0)*(2.d0*f1+f2+f3+2.d0*f4)
	
	! Estimate of the error:
	err = abs(integral-integral2)
	tol = acc+eps*abs(integral)
	
	if (err > sqrt(2.d0)*tol) then
		call adapt(f, a, (a+b)/2.d0, acc/sqrt(2.d0), eps, f1, f2, int1, err1)
		call adapt(f, (a+b)/2.d0, b, acc/sqrt(2.d0), eps, f3, f4, int2, err2)

		integral = int1+int2
		err = sqrt(err1**2 + err2**2)
	endif
end subroutine adapt