Add true mean radius for Sphere, Ellipsoid and TriaxialEllipsoid classes #177
Conversation
MarkWieczorek
changed the title
|
Hey Do you have a reference about the n=10 parameter to be enough precise? Do you have benchmark for different n parameters? |
|
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. |
|
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 output2 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 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 output2 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 Based on these tests, I think that n=50 is still a good default choice. |
|
Hey, thank you for the clear and clean answer :-) Now the situation is much clear |
There was a problem hiding this comment.
@MarkWieczorek very nice! The implementation is elegant. Left some minor comments but can merge as soon as they're handled.
I'm also assuming that the test values came from spherical harmonic models right?
|
Looks great! Merging it in! |
This PR adds a property that computes the "true" mean radius of the ellipsoid as described in #174. The following changes were made:
mean_radiusproperties were renamed tosemiaxes_mean_radius(mean_radius(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