Consider a point with spherical coordinates $\vec{r}_0=(r_0, \theta_0, 0)$. The spherical gaussian distribution centered at $\vec{r}_0$ is $f(\vec{r})=Ne^{|\vec{r}-\vec{r}_0|^2/A}$, where $N$ is the normalizing factor and $A$ is a measure of the spread.
I am interested in the projection of $f$ onto the $\theta$ coordinate, $f_\theta(\theta)=\int_0^\infty\int_0^{2\pi} f(\vec{r})r\sin\theta\,d\phi\,dr$.
Using the euclidean distance formula $\sqrt{(x-x_0)^2 + (y-y_0)^2+(z-z_0)^2}$ and substituting in the cartesian coordinates for $\vec{r}=(r,\theta,\phi)$ with $x=r\cos\phi\sin\theta\;;\;y=r\sin\phi\sin\theta\;;\;z=r\cos\theta$ and similarly for $\vec{r}_0$ with
$x_0=r_0\sin\theta_0\; ; \; y_0=0\;;\;z_0=r\cos\theta.$ Then $f(\vec{r})=Ne^{(-r^2+2rr_0(\cos\phi\sin\theta\sin\theta_0+\cos\theta\cos\theta_0)-r_0^2)/A},$ so $f_\theta(\theta)=\int_0^\infty\int_0^{2\pi} Ne^{(-r^2+2rr_0(\cos\phi\sin\theta\sin\theta_0+\cos\theta\cos\theta_0)-r_0^2)/A}r\sin\theta\,d\phi\,dr.$
Using the fact that $\int_0^{2\pi}e^{x\cos\phi}=2\pi I_0(x)$, where $I_0(x)$ is the modified Bessel function of the first kind, I can simplify $f_\theta(\theta)$ to $f_\theta(\theta)=N\sin\theta2\pi e^{-r_0^2/A}\int_0^\infty e^{(-r^2+2rr_0\cos\theta\cos\theta_0)/A}I_0(2rr_0\sin\theta\sin\theta_0/A)r\,dr$
Is there a way to evaluate this integral further?
I will ultimately be using this as a kernel for kernel density estimation. If it is not possible to further evaluate the integral, I would be happy with an approximate evaluation as long as I will be able to implement it and it was not prohibitively inefficient to evaluate.