Add true mean radius for Sphere, Ellipsoid and TriaxialEllipsoid classes

[MarkWieczorek]

This PR adds a property that computes the "true" mean radius of the ellipsoid as described in #174. The following changes were made:

1. The old Sphere/Ellipsoid/TriaxialEllipsoid mean_radius properties were renamed to semiaxes_mean_radius ($R_1$).
2. A new property mean_radius ($R_0$) was added that computes the mean radius using Gauss Legendre quadrature. The quadrature routine is taken directly from numpy, and requires no additional imports.

Note: There is a single parameter that needs to be defined, which is the number of quadrature points to use. (For triaxial ellipsoids, I set the number of points in longitude to be twice that in latitude). For the ellipsoids that are already defined in boule, tests show that about n=10 is sufficient to get machine precision results. Tests for extreme ellipsoids with a flattening of 0.5 (which is 100 times greater than for Mars) requires n=30. Tests for a Vesta-like object where the semimajor axis is reduced by a factor of two requires n=38. In an extreme abundance of caution, I have chosen to use n=50, which is probably sufficient for any solar system object. I have cross checked these results using the degree-0 spherical harmonic coefficients calculated with the independent code pyshtools, which compare to the last digit (or two).

Relevant issues/PRs: resolves #174

[VascoSch92]

Hey Very interesting

Do you have a reference about the n=10 parameter to be enough precise?

Do you have benchmark for different n parameters?

[MarkWieczorek]

No reference, sorry, those are my calculations :) I just computed the mean radius for all values of n and saw how it converged. Then I verified that the value at large n was equal to the degree-0 spherical harmonic coefficient from an independent code (pyshtools). I could copy the results here if you really want to see them.

If desired, we could make "n" an optional parameter that the user sets, but this would turn the property into a function.

[MarkWieczorek]

I redid the benchmarks after asking the question: What are the most oblate and triaxial objects in our Solar System (and beyond)?

For bi-axial objects, I think that Wenu (which is part of the Arrokoth contact binary, see: https://en.wikipedia.org/wiki/486958_Arrokoth) is probably the most flattened. I took the dimensions from wikipedia (a=21 km, b=9 km) which gives a flattening of f=0.57. Machine precision results are obtained for n=35. The result from pyshtools is listed at the end.

code output

2 13.331032139727506
3 15.426051291782759
4 14.655648468695782
5 14.924707863783672
6 14.825403037262024
7 14.862063866071386
8 14.848315738349147
9 14.853500228581536
10 14.851531399055808
11 14.85228252022791
12 14.85199474739618
13 14.852105373675421
14 14.85206272077701
15 14.852079207403024
16 14.852072820766614
17 14.85207529962679
18 14.852074335850151
19 14.852074711138155
20 14.852074564802848
21 14.85207462193383
22 14.852074599604139
23 14.852074608340661
24 14.852074604919295
25 14.85207460626031
26 14.852074605734273
27 14.85207460594078
28 14.852074605859661
29 14.85207460589154
30 14.852074605878984
31 14.852074605883926
32 14.852074605881999
33 14.852074605882756
34 14.852074605882455
35 14.852074605882592
36 14.852074605882521
37 14.852074605882539
38 14.852074605882544
39 14.852074605882565
40 14.852074605882581
41 14.852074605882573
42 14.852074605882567
43 14.852074605882525
44 14.852074605882532
45 14.852074605882564
46 14.852074605882533
47 14.85207460588253
48 14.852074605882525
49 14.85207460588252

pyshtools: 14.852074605882546

For the most triaxial object, this is probably the interstellar object Oumuamua (see https://en.wikipedia.org/wiki/%CA%BBOumuamua). The dimensions are uncertain, but here is one of the more extreme end-members : a=230, b=35, c=35, with a 6.5:1 aspect ratio. For this case, we don't achieve machine precision for n<100. However, at n=50, we have 6 significant figures, increasing to 10 at n=100.

code output

2 42.62769309788046
3 49.076029915094736
4 47.06587210359279
5 49.4577542401064
6 48.72465260254689
7 49.68645156983228
8 49.45125371234477
9 49.89206958433329
10 49.80170283560964
11 50.01778965363778
12 49.98355983620022
13 50.09380011763802
14 50.082318858120054
15 50.14027496041555
16 50.13761005907834
17 50.168815145990536
18 50.16929404375516
19 50.18643318587273
20 50.1877901941027
21 50.19736896017173
22 50.19875279899406
23 50.20419132725212
24 50.205333736014744
25 50.20846698537194
26 50.209327754317
27 50.211157666389056
28 50.21177482455779
29 50.212857249235824
30 50.21328654759997
31 50.21393444872784
32 50.214227231516674
33 50.214619302701934
34 50.21481631965695
35 50.21505595952294
36 50.21518730179966
37 50.21533510455376
38 50.2154220902638
39 50.21551399432274
40 50.215571336750294
41 50.215628898243324
42 50.215666576975366
43 50.21570286100459
44 50.21572756428916
45 50.21575056574457
46 50.21576673847619
47 50.215781392391484
48 50.21579197107286
49 50.21580134767154
50 50.21580826418095
51 50.215814286959315
52 50.215818808548
53 50.215822690091564
54 50.21582564640475
55 50.21582815535153
56 50.215830088851234
57 50.2158317147966
58 50.215832979912705
59 50.21583403604118
60 50.215834864272445
61 50.21583555167732
62 50.21583609422373
63 50.215836542450084
64 50.21583689808667
65 50.21583719083112
66 50.21583742410865
67 50.215837615584825
68 50.21583776870865
69 50.21583789411396
70 50.21583799469616
71 50.21583807692793
72 50.21583814304252
73 50.21583819702409
74 50.21583824051248
75 50.21583827598549
76 50.21583830461116
77 50.21583832794213
78 50.21583834679657
79 50.21583836215558
80 50.21583837458184
81 50.21583838470056
82 50.215838392896174
83 50.21583839956724
84 50.215838404975656
85 50.215838409377334
86 50.21583841294831
87 50.215838415855394
88 50.215838418213416
89 50.21583842013365
90 50.215838421692986
91 50.21583842296151
92 50.215838423993006
93 50.21583842483249
94 50.21583842551519
95 50.21583842607049
96 50.21583842652272
97 50.21583842689002
98 50.21583842718916
99 50.21583842743262

pyshtools: 50.21583842850183

Based on these tests, I think that n=50 is still a good default choice.

[VascoSch92]

Hey,

thank you for the clear and clean answer :-)

Now the situation is much clear

[leouieda]

Looks great! Merging it in!