IntBesselK Function

FUNCTION IntBesselK(a,b)
	REAL(rp), INTENT(IN) :: a
	REAL(rp), INTENT(IN) :: b
	REAL(rp) :: IntBesselK
	REAL(rp) :: Iold
	REAL(rp) :: Inew
	REAL(rp) :: rerr
	REAL(rp) :: sum_f
	REAL(rp) :: v,h,z
	INTEGER :: ii,jj,npoints
	LOGICAL :: flag

	v = 5.0_rp/3.0_rp
	h = b - a
	sum_f = 0.5*(besselk(v,a) + besselk(v,b))

	Iold = 0.0_rp
	Inew = sum_f*h

	ii = 1_idef
	flag = .TRUE.
	do while (flag)
		Iold = Inew

		ii = ii + 1_idef
		npoints = 2_idef**(ii-2_idef)
		h = 0.5_rp*(b-a)/REAL(npoints,rp)
		sum_f = 0.0_rp
		do jj=1_idef,npoints
			z = a + h + 2.0_rp*(REAL(jj,rp) - 1.0_rp)*h
			sum_f = sum_f + besselk(v,z)
		end do

		Inew = 0.5_rp*Iold + sum_f*h
		rerr = ABS((Inew - Iold)/Iold)
		flag = .NOT.(rerr.LT.Tol)
	end do
	IntBesselK = Inew
END FUNCTION IntBesselK