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

Loss of precision in time argument #10

Closed
MicheleCeresoli opened this issue Oct 24, 2023 · 1 comment · Fixed by #11
Closed

Loss of precision in time argument #10

MicheleCeresoli opened this issue Oct 24, 2023 · 1 comment · Fixed by #11
Assignees
Labels
bug Something isn't working

Comments

@MicheleCeresoli
Copy link
Member

When computing the MOON_PA orientation, a loss of precision in the time arguments arises when converting from TDB seconds since J2000 to TDB days since J2000.

CALCEPH and SPICE instead have exact results. MWE:

using Ephemerides 
using CalcephEphemeris
using JSMDInterfaces.Ephemeris
using Tempo 
using SPICE
using ReferenceFrameRotations

kernel =  "test/assets/moon_pa_de440_200625.bpc"

ephj = EphemerisProvider(kernel);
ephc = CalcephProvider(kernel);

et = rand(0.0:1e8)
et = 1.0

yj, yc = zeros(3), zeros(3)
ephem_orient!(yj, ephj, DJ2000, et/86400.0, 31008, 1, 0)
ephem_orient!(yc, ephc, DJ2000, et/86400.0, 31008, 1, 0)

yj2 = ephem_rotation3(ephj, 1, 31008, et)

A = DCM(pxform("J2000", "MOON_PA", et))
angs = dcm_to_angle(A, :ZXZ)

yj = mod.(yj, π)
yj2 = mod.(yj2, π)
yc = mod.(ys, π)
ys = mod.([angs.a1, angs.a2, angs.a3], π)

yj-yc
yj2-yc 

yj-ys
yj2-ys

yc-ys

# yj2 is more precise than yj!
yj2 - ys 
yj2 - yj
@MicheleCeresoli MicheleCeresoli added the bug Something isn't working label Oct 24, 2023
@MicheleCeresoli MicheleCeresoli self-assigned this Oct 24, 2023
@MicheleCeresoli MicheleCeresoli changed the title Evaluate loss of precision in time argument Loss of precision in time argument Oct 24, 2023
@MicheleCeresoli
Copy link
Member Author

MicheleCeresoli commented Oct 24, 2023

Converting time as et = (jd0 - DJ2000) + time greatly improves the accuracy (yj and yj2 become equal) but does not completely solve the differences against SPICE and CALCEPH (i.e., yj2 - ys). The relative difference is about 1e-16 and seems to be related to the magnitude of the involved values.

EDIT: I've also tested the same algorithm that SPICE uses to interpolate the Chebyshev polynomials and the results are basically the same (just very slightly better in terms of accuracy, maybe because it starts from the highest-order coefficients?)

This issue is being closed with #11 because the differences produce a rotation error that is smaller than 1 micro arcsecond.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

Successfully merging a pull request may close this issue.

1 participant