Computing the volume of a region on the unit $n$-sphere

I would like to compute the surface volume of a region on the unit $n-1$-sphere:

$$\sum_{i=1}^n x_i^2 = 1,$$

bounded by an ellipsoid

$$\sum_{i=1}^n a_ix_i^2 \leq a_2,$$

where $1=a_1 < a_2 < \dots < a_n$. A couple nice properties are that

1. The region is symmetric across each axis, so I need only worry about computing the volume in one orthant and then multiply by $2^n$, and
2. the limits of integration are somewhat easy to compute.

My approach so far has been to use Cartesian coordinates ($n$-spherical seemed a bit too daunting with all the powers of sines):

$$\int_0^{\ell_2} dx_2 \int_0^{\ell_3(x_2)} dx_3 \dots \int_0^{\ell_n(x_2,\dots,x_{n-1})} \frac{dx_n}{\sqrt{(1-x_2^2-\dots - x_n^2)}},$$

where $\ell_2 = 1$, $\ell_3(x_2) = \sqrt{\frac{(a_2-a_1) - (a_2-a_1)x_2^2}{a_3-a_1}}$, $\ell_4(x_2, x_3) = \sqrt{\frac{(a_2-a_1) - (a_2-a_1)x_2^2 - (a_3-a_1)x_3^2}{a_4-a_1}}$, and so on, so $\ell_i(x_2, \dots, x_{i-1}) = \sqrt{\frac{1}{a_i-a_1}\left((a_2-a_1) - \sum_{j=2}^{i-1}(a_j-a_1)x_j^2\right)}$.

The problem for me is that as soon as I compute the most internal integral, I get a rather nasty formula involving arcsin:

$$\sin^{-1}\left(\sqrt{\frac{(a_2-a_1) - \sum_{j=2}^{n-1}(a_j-a_1)x_j^2}{(a_n-a_1)(1-\sum_{j=2}^{n-1}x_j^2)}}\right)$$

Since this is the innermost formula, I don't even know where to go from here. Playing tricks with the trig doesn't seem to help either.

My questions are:

1. Is there a better way to attack this? Is it even possible to get a solution in terms of the $a_i$?
2. Is there an approximation to the value of the integral that I can work out in terms of the $a_i$ (perhaps with $O()$-notation)?

I should also note that $n$ could easily be extremely large (on the order of thousands or millions), so I'd need a very general way to attack this problem.

You can write your area as a probability times the area of $S^n$, by writing $X_i = Z_i/\sqrt{\sum_{j=1}^n Z_j^2}$ for i.i.d. $N(0,1)$ gaussian random variables. Then your desired area is $2\pi^{n/2}/\Gamma(n/2)$ times the probability that $\sum_j a_j Z_j^2 \le a_2\sum_j Z_j^2$, which might be estimated (say) by a form of the central limit theorem.