Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Results from scripts/integrate_ephemeris.jl change from Julia v1.6 to Julia v1.8 #10

Open
LuEdRaMo opened this issue Mar 8, 2023 · 2 comments

Comments

@LuEdRaMo
Copy link
Collaborator

LuEdRaMo commented Mar 8, 2023

The results from scripts/integrate_ephemeris.jl change from Julia v1.6 to Julia v1.8 due to differences in rounding. For instance, consider the Moon's radius R_moon = 1.1617812418502559e-5. We need R_moon ^ 5 to compute the Moon's inertia tensor, however:

  • In Julia 1.6.7:
julia> R_moon ^ 5
2.1165171911845607e-25
  • In Julia 1.8.5:
julia> R_moon ^ 5
2.116517191184561e-25

Fortunately, this does not have a huge impact on the ephemeris. The following figure plots the logarithm of the absolute difference between integrations in Julia v1.6 and Julia v1.8. The blue (red) curves represent the position [au] (velocity [au/day]) with highest and lowest error; i.e. all positions (velocites) lie within the two curves. The third Euler angle [rad] for both the Moon's core and mantle (green lines in the figure) are the leading sources of error. All other Euler angles lie within the blue/red hatched region.

PE_v16_vsv18

@lbenet
Copy link
Collaborator

lbenet commented Mar 9, 2023

Thanks for reporting this subtlety. I guess this is related to some internals of Julia, which involve the power function, and that may have changed in Julia 1.8.

Are you using here the current master of NEOs.jl?

Aside from that, I'd like to note that (in Juliav1.8):

julia> big(R_moon)^5
2.116517191184560855255820794719949721474183242017690151038850992707508088212283e-25

julia> R_moon*R_moon*R_moon*R_moon*R_moon
2.1165171911845607e-25

It seems then that the nearest Float64 is the one produced in Julia v1.6 (and v1.7)...

@PerezHz
Copy link
Owner

PerezHz commented Apr 2, 2023

Thanks for spotting this! In a 10-year integration, the (absolute) difference between PlanetaryEphemeris.jl and DE430 is already around the 0.1 arcsec (i.e., about 1e-6 rad) level for the lunar mantle angles, so I agree that this more reflective of a difference between julia 1.6 and 1.8, while not having much impact in the ephemeris. In extended-body acceleration terms, maybe it would be better to do e.g. (R_M/r)^5 rather than (R_M^5)/(r^5) as it is done now, but I doubt this would make a huge difference wrt how things are now... still, it's good to keep track of this numerical subtleties!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants