"Mean radius" proposal

[MarkWieczorek]

Description of the desired feature:

In geodesy, there are several definitions of "radius" or "mean radius", and Boule implements some of these. The standard definitions that one finds for WGS84 are

• R_1 : mean radius of semi-axes, i.e., (a+b+c)/3.
• R_2 : radius of a sphere with the same area as the ellipsoid.
• R_3 : radius of a sphere with the same volume as the ellipsoid.

At the present time, we implement "R_1" and call this mean_radius. There is a PR open to add R_3 which is called volume_equivalent_radius.

The problem is that R_1 is not the mean radius of the ellipsoid. I was surprised to learn this when comparing the degree 0 spherical harmonic coefficient of an ellipsoid with R_1, but perhaps I shoudn't have been surprised...

My proposal is the following:

1. Rename mean_radius (R_1) to semiaxes_mean_radius (R_1)
2. Add a new definition of the true mean radius called mean_radius with the symbol R_0.
3. This would leave us with R_0, R_1, R_2, R_3, with names of mean_radius, semiaxes_mean_radius, area_equivalent_radius and volume_equivalent_radius.

As far as I can tell, there is not an analytic equation for the mean radius of an ellipsoid. Therefore I propose that we compute this numerically using Gauss-Legendre quadrature. This shouldn't take up more than about 5 lines of code, and there are two routines that we could use out of the box. The first is scipy.integrate.fixed_quad and the second is numpy.polynomial.legendre.leggauss.

The scipy method is probably one line of code (for the integration), whereas for the numpy version, you first need to compute the points where the function is evaluated and then do a sum (so two lines of code). We already use numpy, so we wouldn't need to import another module. However, we don't use scipy, so we would need to import this for this for one line of code. Since the radius is very close to being a degree-2 polynomial, using 50 integration points is probably sufficient to get machine precision accuracy (of course I will verify this): numpy claims that their routine is tested to degree 100, so that is probably good enough.

If we do this, we could probably also compute the area of the ellipsoid numerically using a similar approach.

Lastly, I note that these issues are really not important for Earth. However, for Vesta, we have the following:

a = 286.300
b = 278.600
c = 223.200

R_0 = 259.813
R_1 = 262.700
R_3 = 261.115

So there are differences of more than a km. And there are bodies that are even more triaxial than Vesta...

Are you willing to help implement and maintain this feature? Yes.

[leouieda]

@MarkWieczorek thank you for all the thought you put into this! I'm definitely on board with having more accurate property names.

I don't know if we use the mean radius anywhere that is significant so it probably wouldn't hurt to just swap these out in the next release. Or maybe we add a warning to mean_radius that it had changed.

Could you share a reference for the R0 definition you had?

[MarkWieczorek]

I don't have a reference, but this should help make it clear:

The mean value of a function $r$ in spherical coordinates (let's just call this $R_0$) is by definition

$$R_0 = \frac{1}{4\pi} \int_0^{\pi} \int_0^{2\pi} r(\theta, \phi) \sin \theta d\theta d\phi$$

Next, compare this with the definition of the spherical harmonic coefficients of a function $f$, which are defined to be

$$f_{lm} = \frac{1}{4\pi} \int_0^{\pi} \int0^{2\pi} f(\theta, \phi) Y{lm}(\theta, \phi) \sin \theta d\theta d\phi$$

For degree and order 0, $Y{00}$ is unity, and the two equations are thus identical. This just shows that the mean radius is equal to the $f{00}$ spherical harmonic coefficient.

For an ellipsoid, the equation for the geocentric radius $r(\theta, \phi)$ is somewhat complicated, and I don't believe that an analytic expression exists. The integration over longitude is trivial for a biaxial ellipsoid, but the integration over latitude needs to be done numerically.