JuliaSpaceMissionDesign / Ephemerides.jl

A Modern Binary Ephemeris Reader for Julia, in Julia.
https://juliaspacemissiondesign.github.io/Ephemerides.jl/
MIT License
22 stars 1 forks source link

Loss of precision in time argument #10

Closed MicheleCeresoli closed 11 months ago

MicheleCeresoli commented 11 months ago

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 commented 11 months ago

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.