Projection of Gaussian in Spherical Coordinates 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.  
 A: This integral has a quite simple evaluation in the limit $A\rightarrow 0$, very small spreading. Indeed, in the distribution sense,
$$\lim_{A\rightarrow 0}\frac{1}{\sqrt{\pi A}}e^{-\frac{(r-r_0\cos\theta\cos\theta_0)^2}{A}}=\delta(r-r_0\cos\theta\cos\theta_0)$$
and so
$$f_\theta(\theta)=N'\sin\theta 2\pi e^{-\frac{r_0^2}{A}(1-\cos^2\theta\cos^2\theta_0)}\int_0^\infty\delta(r-r_0\cos\theta\cos\theta_0)I_0\left(\frac{2}{A}rr_0\sin\theta\sin\theta_0\right)rdr$$
with $N'$ again a proportionality constant and I have assumed it will depend on $A$ from the start due to normalization (this should be given explicitly in the question). At this point the evaluation of the integral should be a simple matter.
One can also extend this simple argument to a more general asymptotic expansion in $\frac{1}{A}$ with standard techniques. The idea is this. Consider the integral
$$I(\sigma)=\frac{1}{\sqrt{2\pi}\sigma}\int_0^r e^{-\frac{(r-r_0)^2}{2\sigma^2}}\phi(r)r dr$$
and your aim is to compute the asymptotic series for $\sigma\rightarrow 0$. Now, change the variable to $y=\frac{r-r_0}{\sqrt{2}\sigma}$ and you get
$$I(\sigma)=\frac{1}{\sqrt{\pi}}\int_{-\frac{r_0}{\sqrt{2}\sigma}}^\infty e^{-y^2}\phi\left(\sqrt{2}\sigma y+r_0\right) \left(\sqrt{2}\sigma y+r_0\right)dy$$
and you can expand $\phi\left(\sqrt{2}\sigma y+r_0\right)$ around $\sigma=0$ obtaining
$$I(\sigma)=\frac{1}{\sqrt{\pi}}\sum_{n=0}^\infty\frac{\phi^{(n)}(r_0)}{n!}(\sqrt{2}\sigma)^n\int_{-\frac{r_0}{\sqrt{2}\sigma}}^\infty e^{-y^2}y^n\left(\sqrt{2}\sigma y+r_0\right)dy.$$
Finally, extend the lower end of the integration interval to $-\infty$ and you are left with the evaluation of the integrals
$$\int_{-\infty}^\infty e^{-y^2}y^k dy$$
that are given by the $\Gamma$ function and so are well known. Note that the first term of this series is just the one I have obtained with the Dirac distribution. The other I have obtained now are all the corrections in powers of $A$.
