Question: Numerical calculation of a Hankel transform

I want to calculate the numerical value of the following integral at a given (r,z) where r and z are real numbers:

g(r,z)=int(exp(-sqrt(k^2 - k0^2)*z)*BesselJ(0, k*r)*k/sqrt(k^2 - k0^2), k = 0 .. infinity) where k0 is a real number.

In other words, it is the hankel transform g(r,z) of the function f(k):

f(k)=exp(-sqrt(k^2 - k0^2)*z)/sqrt(k^2 - k0^2)

The theoretical result is exp(-I*k0*R)/R with R=sqrt(r^2+z^2)

I tried to use the commande evalf(Int(exp(-sqrt(k^2 - k0^2)*z)*BesselJ(0, k*r)*k/sqrt(k^2 - k0^2), k = 0 .. infinity)) but it takes hours to evaluate and I did not get any results yet. I tried to calculate this integral with another software which uses a vectorized adaptive quadrature algorithm and I got the result in a few seconds.

How could I calculate this integral with Maple ?

Thanks.