Integral over three Bessel functions $J_{0} ( \alpha\ x ) K_{iP} ( \beta x ) K_{iP} ( \xi x )$ I'm considering the following integral:
$I(\alpha,\beta,\xi) = \int_{0}^{\infty} dx\ x\ J_{0} ( \alpha\ x )\ K_{iP} ( \beta x )\ K_{iP} ( \xi x )$
Where $P >0$, and $\alpha > 0$ and $\beta, \xi \in \mathbb{R}$. Is there any way to do an integral over three Bessel functions? I have seen plenty of examples for two Bessel functions up until this point.
Is there any way to evaluate the above?
EDIT: I found this link which is basically the same as the above, except with $K_{iP} \mapsto J_{iP}$.
 A: Eq. 6.578.10 of Gradshteyn and Ryzhik quting Ederlyi TII p.67 (30) expresses the result with an associated Legendre function: $$
\frac{\sqrt{\pi}}{2^{3/2}}\frac{\Gamma(1+iP)\Gamma(1-iP)}{\beta\xi (u^2-1)^{1/4}}P_{iP-1/2}^{-1/2}(u)
$$
with $2\beta\xi u=\alpha^2+\beta^2+\xi^2$ when $\Re(\xi)>0,\Re(\beta)>0,\alpha>0$. 
This result can be obtained by representing (DLMF)
\begin{equation}
 \ K_{iP} ( \beta x )\ K_{iP} ( \xi x )=\frac{1}{2}\int_0^\infty \exp\left( -\frac{t}{2}-\frac{\beta^2+\xi^2}{2t}x^2 \right)K_{iP}\left(\frac{\beta\xi}{t}  x^2\right)\frac{dt}{t}
\end{equation} 
With $t=x^2v$, and changing the order of integration it follows
\begin{equation}
  I(\alpha,\beta,\xi) = \frac{1}{2}\int_0^\infty \exp\left(-\frac{\beta^2+\xi^2}{2v} \right)K_{iP}\left(\frac{\beta\xi}{v}  \right)\frac{dv}{v}\int_{0}^{\infty} dx\ x\ J_{0} ( \alpha\ x )\exp\left( -\frac{x^2v}{2} \right)
\end{equation} 
Inner integral is a classical one: (DLMF)
\begin{equation}
 I(\alpha,\beta,\xi) = \frac{1}{2}\int_0^\infty \exp\left(-\frac{\alpha^2+\beta^2+\xi^2}{2v} \right)K_{iP}\left(\frac{\beta\xi}{v}  \right)\frac{dv}{v^2}
\end{equation} 
With $v=\frac{\beta\xi}{w}$ and $u$ defined as above,
\begin{equation}
 I(\alpha,\beta,\xi) = \frac{1}{2\beta\xi}\int_0^\infty \exp\left(-uw \right)K_{iP}\left(w \right)\,dw
\end{equation} 
which is a Laplace transform of the Bessel function. Using an integral representation for this function (DLMF)
\begin{equation}
 I(\alpha,\beta,\xi) = \frac{1}{2\beta\xi}\int_0^\infty \,dw\exp\left(-uw \right)\int_0^\infty \exp\left( -w\cosh t \right)\cosh\left( iPt \right)\,dt
\end{equation} 
and thus
\begin{equation}
 I(\alpha,\beta,\xi) = \frac{1}{2\beta\xi}\int_0^\infty \,dt\frac{\cosh\left( iPt \right)}{u+\cosh t}\,dt
\end{equation} 
This integral may be expressed using an integral representation of an associated Legendre function (DLMF) giving the above result.
EDIT 29/04/2017: Obtained result can be simplified by using the Gamma reflection formula and, noting that, if $\alpha,\beta,\xi$ are real, then $u>1$. Denoting $u=\cosh z$, the given Legendre associated function can be simplified (DLMF) to express
$$ 
 I(\alpha,\beta,\xi) = \frac{\pi}{2\beta\xi}\frac{\sin Pz}{\sinh \pi P\sinh z}
$$
