I'm trying to calculate hit probabilities on a dart board if the dart thrower has some Gaussian angle distribution function with width $\Delta$ and some systematic angle offsets $\phi_0, \theta_0$. This is a problem related to my work.

The distribution function is

$p(\theta, \phi) = \frac{1}{2 \pi \Delta^2} e^{-\frac{(\phi-\phi_0)^2}{2\Delta^2}}e^{-\frac{(\theta-\theta_0)^2}{2\Delta^2}}$

To get the hit probability on the dart boad, this distribution function must be integrated over a circle centered around 0 with some radius $R$ at a far away distance $D$.

So far I turned in into an integral in circle coordinates and now I'm stuck with this bad boy:

$p_{tot} = e^{-\rho_0^2} \int_0^{\frac{R}{\sqrt{2} D \Delta}} d\rho \rho e^{-\rho^2} I_0(2\rho \rho_0)$

where $\rho = \sqrt{\phi^2 + \theta^2}$ and $\rho_0$ equivalently. $I_0$ is the 0th modified Bessel function.

Any ideas on how to properly evaluate this? I'd be happy about a reasonable approximation. Taylor expansions don't really work because the function in the integral has a peak that shifts to the right with increasing $\rho_0$. The area under the peak will grow like $e^{\rho_0^2}$ which can lead to some numerical issues but I guess I can approximate it with a delta distribution for really large $\rho_0$.