Integral of Multivariable Gaussian across a Circular Domain. Basically, trying to compute the probability, a.k.a volume under the bi-variate gaussian distribution: $$ f(x,y) = \frac{1}{2 \pi \sigma^2} \cdot \exp(\frac{-x^2 - y^2}{2 \sigma^2}) $$ over axis-aligned circles in the xy-plane. So, the circles all are of the form: $(x - c_x)^2 + y^2 = R^2$ So the centers always lie on the x-axis somewhere, and the gaussian always has mean $(x, y) = (0, 0)$
I've been using a numerical integration library which computes the following integral (setup in cartesian notation):
$$ 2 \cdot \int_{x = c_x-R}^{c_x + R}\int_{y=0}^{\sqrt{R^2 - (x - c_x)^2}} f(x, y) \, dy\, dx $$
But I also know that computing this is very easy when the circle is centered at the origin (i.e. when $c_x = 0$). Since the integral can just be written in polar coordinates and turns into:
$$ \frac{1}{2 \pi \sigma^2} \int_{\theta = 0}^{2\pi}\int_{r=0}^R re^{-r^2 / 2\sigma^2} \, dr\, d\theta = \frac{1 - e^{-R^2 / 2\sigma^2}}{2 \pi \sigma^2}$$
Buut, once I introduce the $c_x$ term and the circle is not centered at the origin, the polar integral gets a weird term added on, making things really scary lol. I don't know how I would go about finding a closed form, if one even exists? Here is what the polar integral becomes when $c_x \neq 0$:
$$ \frac{1}{2 \pi \sigma^2} \int_{\theta = 0}^{2\pi}\int_{r=0}^R r e^{-r^2 / 2\sigma^2}e^{2r\cdot c_x \cdot\cos\theta / 2\sigma^2} \, dr\, d\theta$$
Anybody have a clue if/how I can find a closed form for the last integral? or perhaps that first integral, but I'd imagine it's harder from the cartesian perspective? 
Any help appreciated,
Cheers :-)
 A: There isn't much you can do, but here are a few comments.  First, the result $f(\sigma,c_x,R)$ is invariant under scaling:
$$f(\sigma,c_x,R) = f(t\sigma,tc_x,tR)$$
for all $t > 0$.  Therefore you may solve the problem for the $\sigma = 1$ case without loss of generality:
$$f(1,c_x,R) = \frac{e^{-c_x^2/2}}{2\pi} \int_0^R \int_0^{2\pi} r e^{-r^2/2} e^{r c_x \cos \theta} \,d\theta dr.$$
The inner integral may be evaluated in terms of $I_0(x)$, a modified Bessel function of the first kind:
$$f(1,c_x,R) = e^{-c_x^2/2} \int_0^R r e^{-r^2/2} I_0(c_x r)\,dr.$$
That's about it.  This integral doesn't appear to simplify for finite $R$, nor does integrating in the other order produce anything simpler.  If you'd prefer a series to a sum you could use the Taylor series for $I_0(x)$ to perform a term-by-term integration
$$f(1,c_x,R) = e^{-c_x^2/2} \sum_{k=0}^\infty \int_0^R r e^{-r^2/2} \frac{(c_xr/2)^{2k}}{k!^2} \,dr.$$
This simplifies to
$$f(1,c_x,R) = 1 - e^{-c_x^2/2} \sum_{k=0}^\infty \frac{(c_x^2/2)^k}{k!^2} \Gamma(k+1,R^2/2),$$
where $\Gamma(a,x)$ denotes the incomplete Gamma function.
