I'm trying to get the probability density function of the range from Gaussian distributed $\mathrm{\bar{y}} = [x,y,z]$ coordinates (covariance matrix identity). I know,
\begin{align} x &= r.cos(\alpha)cos(\theta)\\ y &= r.cos(\alpha)sin(\theta)\\ z &= r.sin(\alpha) \end{align}
where $r$ is the range, $\theta$ azimuth and $\alpha$ is the elevation.
I transformed the Gaussian distribution and arrived at the following joint distribution in spherical coordinates, \begin{align} p(r,\theta,\alpha) = \frac{r^2cos(\alpha)}{(2\pi)^{3/2}} \exp(-1/2 (r^2 + \mu^T\mu))\exp(r\mu^T\beta) \end{align}
where $\mu$ is the mean of Gaussian distributed variables, and $\beta$ is, \begin{align} \beta = \begin{bmatrix} cos(\alpha)cos(\theta)\\ cos(\alpha)sin(\theta)\\ sin(\alpha) \end{bmatrix} \end{align}
Then I marginalized $\theta$ and $\alpha$ to get the integral, \begin{align} p(r) = \frac{r^2}{(2\pi)^{3/2}} \exp(-1/2 (r^2 + \mu^T\mu)) \int_{-\pi/2}^{\pi/2} cos(\alpha)\int_0^{2\pi} \exp(r\mu^T\beta)d\theta d\alpha \end{align}
Using $\mu^T = [a,b,c]$, the inner integral can be solved. Hence \begin{align} p(r) = \frac{r^2cos(\alpha)}{(2\pi)^{3/2}} \exp(-1/2 (r^2 + \mu^T\mu)) \int_{-\pi/2}^{\pi/2} cos(\alpha) \exp(rcsin(\alpha)) \mathrm{I}_0(r\sqrt{a^2+b^2}cos(\alpha)) d\alpha \end{align} Where $\mathrm{I}_0(.)$ is the modified Bessel function of the first kind and zeroth order. There shouldn't be any mistakes thus far but I cannot seem to make any progress now. Is there a general form of such integral, \begin{align} \int_{-\pi/2}^{\pi/2} cos(\alpha) \exp(rcsin(\alpha)) \mathrm{I}_0(r\sqrt{a^2+b^2}cos(\alpha)) d\alpha \end{align}
Later on I would probably also have to calculate the integral, \begin{align} \int_{-\pi/2}^{\pi/2} cos^2(\alpha) \exp(rcsin(\alpha)) \mathrm{I}_1(r\sqrt{a^2+b^2}cos(\alpha)) d\alpha \end{align}
Notice the order of modified Bessel function has increased.
Thanks.