Double integral of product of Bessel functions in Python I would like to perform a numerical integral (in Python) of the form 
$$\displaystyle F_\nu(k) = \int_0^\infty dk' \ k' \int_0^\infty dr\ r \ f(k', r) \ J_\nu(k' r) \ J_\nu(k r),$$
where $J_\nu$ is the $\nu$-th order Bessel function of the first kind. 
Moreover, note that the Bessel functions form an orthogonal basis with respect to the weighting factor $r$, i.e.
$$ \int_0^\infty dr \ r \ J_\nu(k'r) \ J_\nu(kr) = \frac{\delta(k' - k)}{k}, \quad k,\ k'>0,$$
where $\delta$ is the Dirac delta. 
Although this is not impossible by naively performing this integral in Python, it takes a long time to evaluate. However, given the fact that the Bessel functions form an orthogonal basis, and since since this form somewhat resembles a Hankel transform (for which packages are available that perform an order of magnitude better than their naive counterparts) I was hoping there might be faster way to evaluate this integral. 
Any ideas on how to tackle this integral?
 A: $$
\int_0^\infty k' dk' \int_0^\infty r dr f(k',r) J_\nu (k' r) J_\nu (kr)\\
f(k',r) = \int_0^\infty k'' dk'' F_\mu (k',k'') J_\mu (k'' r)\\
\int_0^\infty k' dk' \int_0^\infty r dr \int_0^\infty k'' dk'' F_\mu (k',k'') J_\mu (k'' r) J_\nu (k' r) J_\nu (kr)\\
$$
Where we have substituted $F_\mu(k',k'')$ which is a Hankel transformation of $f(k',r)$
The integral
$$
D((m,\gamma),(m',\gamma'),(m'',\gamma'')) \equiv \int_0^\infty r dr J_m (\gamma r) J_{m'} (\gamma' r) J_{m''} (\gamma '' r)
$$
was studied in Auluck for integer orders.
So assuming I can Fubini the $r$ and $k''$ integrals
$$
\int_0^\infty k' dk' \int_0^\infty k'' dk'' F_\mu (k',k'') \int_0^\infty r dr J_\mu (k'' r) J_\nu (k' r) J_\nu (kr)\\
$$
I'm assuming $\nu$ is an integer now.
$$
\int_0^\infty k' dk' \int_0^\infty k'' dk'' F_\mu (k',k'') D((\mu ,k''),(\nu,k'),(\nu,k))\\
$$
Some examples of $D$ have exact solutions. For example,
$$
D((0,\gamma),(m',\gamma'),(m',\gamma''))
$$
is at the bottom of page 6. But that gets you some expression in Legendre functions and $\frac{k^2+k'^2-k''^2}{2kk'}$, so I don't know how useful that will be.
You could use the approximation on page 24 which gives you an approximation to $D$ as a distribution on $k''$ for smooth weight functions which would hopefully include $F_\mu (k',k'')$. Note however, that formula is not proven there.
