Numerical approximation of the modified Bessel function $I_0$ with radical argument for integration purposes I have to numerically calculate the following definite integral

$$\int_{\alpha}^{\beta}I_0(a\sqrt{1-x^2})dx$$

for different values of $\alpha$ and $\beta$, where $a$ has a value of, say, $30$. I'm not using quadrature rules, for the sake of calculation speed. Conversely, I'm trying to approximate the integrand so to obtain a closed form primitive.

According to 

J.M. Blair, "Rational Chebyshev approximations for the modified Bessel functions I_0(x) and I_1(x)", Math. Comput., vol. 28, n. 126, pp. 581-583, Apr. 1974.

when the argument of $I_0$ is less than $15$, then $I_0$ can be approximated as a rational function involving only even powers of $x$. In this case, the presence of $\sqrt{1-x^2}$ introduces no difficulty, being $I_0(a\sqrt{1-x^2})$ still a rational function, so that the integration can be done in closed form for the interval $(\alpha,\sqrt{1-(15/a)^2})$.

The problem arises for the interval $(\sqrt{1-(15/a)^2},\beta)$. In this case, according to Blair's paper,

$$I_0(a\sqrt{1-x^2})\simeq \frac{e^{-a\sqrt{1-x^2}}}{(1-x^2)^{1/4}}p(a\sqrt{1-x^2})$$

where $p$ is a Chebyshev polynomial. In this case, closed form integration is not possible anymore (myself and Mathematica are not able to find a closed form primitive). I have also tried with Mathematica to find a rational function approximation for the range $(\sqrt{1-(15/a)^2},\beta)$. Unfortunately, I was not successful to achieve a satisfactory solution (of the order of machine double precision accuracy), perhaps due to my bad skills in using Mathematica.

Is anyone aware of any other approximation of $I_0(a\sqrt{1-x^2})$ in $(\sqrt{1-(15/a)^2},\beta)$ (or even of a better usage of Mathematica) which could be integrated in a closed form?
There is a nice expression for the modified Bessel function of the first kind of order zero:

$$I_0(z)=\frac1{\pi}\int_0^\pi \exp(z\cos\,u)\mathrm du$$

The nice thing about this expression is that one can use the simple-minded trapezoidal rule or midpoint rule to evaluate the integral numerically, with good accuracy; for an explanation of why this is so, see e.g. this article or this article.

As an aside, the substitution $x=\sin t$ should improve the numerical behavior of your integral.
In fact the best approach should be first expand the bessel function and then integrate it term by term.

$\int_\alpha^\beta I_0(a\sqrt{1-x^2})~dx$

$=\int_\alpha^\beta\sum\limits_{n=0}^\infty\dfrac{a^{2n}(1-x^2)^n}{4^n(n!)^2}dx$

$=\int_\alpha^\beta\sum\limits_{n=0}^\infty\sum\limits_{k=0}^n\dfrac{a^{2n}C_k^n(-1)^kx^{2k}}{4^n(n!)^2}dx$

$=\int_\alpha^\beta\sum\limits_{n=0}^\infty\sum\limits_{k=0}^n\dfrac{(-1)^ka^{2n}x^{2k}}{4^nn!k!(n-k)!}dx$

$=\left[\sum\limits_{n=0}^\infty\sum\limits_{k=0}^n\dfrac{(-1)^ka^{2n}x^{2k+1}}{4^nn!k!(n-k)!(2k+1)}\right]_\alpha^\beta$

$=\sum\limits_{n=0}^\infty\sum\limits_{k=0}^n\dfrac{(-1)^ka^{2n}(\beta^{2k+1}-\alpha^{2k+1})}{4^nn!k!(n-k)!(2k+1)}$