time representation accuracy

[bgodard]

AstroTime.jl represents an epoch as a julian day value separated in 2 float64 numbers jd1 and jd2: jd1+jd2. While in theory this representation can allow a sub-nanosecond accuracy for contemporary epochs, the accuracy is only guaranteed if the operations on epochs are designed to ensure that accuracy, but this is not the case here. Taking as example a round trip light time to a probe in the far solar system (such as Voyager 1 and 2) which could be something like 1.5 days. If the reception time TR is known, we compute the transmission time TT:

julia> using AstroTime

julia> TR = TDBEpoch("2018-06-01T00:00:00.0") # reception time
2018-06-01T00:00:00.000 TDB

julia> RTLTa = 1.5 * days                        # round trip light time
1.5 days

julia> RTLTb = (1.5+1e-6/86400.)* days # 1 microsecond different from LTa
1.500000000011574 days

julia> TTa = TR - RTLTa                          # transmission time
2018-05-30T12:00:00.000 TDB

julia> TTb = TR - RTLTb 
2018-05-30T12:00:00.000 TDB

julia> TTb-TTa
0.0 days

In this example, we see that the time representation is in error by one microsecond after a trivial operation.

One microsecond accuracy is not sufficient for application like deep space probe navigation. Ideally you want nanosecond accuracy or better in absolute time representation. Some navigation software actually manage to work with about 50 nanoseconds accuracy in absolute time but this requires some additional complexity and/or introduce certain modelling limitations.

Note that time representation is not in general as accurate as (short) duration representation: indeed the light time is normally modeled to nanosecond or better accuracy while the computed explicit difference between transmission and reception time is of the accuracy of absolute time representation.

Another aspect to take into account is that to achieve accurate modelling of deep space navigation observable, it is not sufficient to have accurate epoch representation but the interpolation of trajectory states should make use of the extra epoch precision:

• for JPL DExxx and OBSPM INPOP ephemerides, this can be obtained if you decompose your epoch in JD1 and JD2 with JD1 integer and 0<=JD2<1.0 but this is not guaranteed in general for other decompositions.
• for custom numerically integrated trajectory, the independent argument could be eg a float64 storing the time (in days or seconds) since the start of the numerical integration interval and the accurate starting epoch could be stored separately.

The following considerations should be taken into account for the choice of time representation and operations on time:

• absolute accuracy better than a nanosecond for contemporary epochs (accuracy is not so important in far future and far past (since you do not have any accurate observations in those times). This is in favor of choosing a contemporary epoch as a reference instead of JD 0.0 if using floating point time representation since floating point absolute accuracy decrease with the size of the stored value. Most navigation software use reference 2000/01/01 or 1950, 1970...
• accuracy maintained when operating between epochs and small duration (up to a few days). It could also work for long duration if duration is stored in extended precision similarly to epoch but this is not as important and can introduce some additional computation costs.
• sufficiently large range of epochs (eg if you are interested in cosmological scale although there is not much use of the library for this since all available timescales are irrelevant before there was a solar system).
• provide routines which decompose accurately the time in two floating point numbers for best accuracy interpolation of JPL and OBSPM ephemerides as well as for SOFA routines.
• finally efficiency considerations.

All of the above can be achieved with the current JD1, JD2 double double storage but I am not sure this is the easiest to implement and/or the most CPU usage efficient.

Some other options: (1) seconds since reference epoch (Integer 64) and femtoseconds in second (Integer 64): representation accurate to 0.5 femtoseconds for all epochs (would also work using attoseconds but more tricky to implement since the sum of two fraction of seconds could then overflow in Int64). (2) seconds since reference epoch (Integer 64) and fraction of second (double): representation accurate to 55 attoseconds (worst case if 0.5<=second fraction<1) (3) days since reference epoch (Integer 64) and seconds in days (double): representation accurate to 7 picoseconds (worst case if seconds in day>=65536)

Options (1) and (2) using "seconds since" can store epoch more than 10 times the age of the universe from reference epoch. Option (3) using "days since" can store a million times the age of the universe...

The avantage of using integer is that representation accuracy does not get worse as you get further away from the reference epoch.

Example from Orekit AbsoluteDate.java:

   /** Reference epoch in seconds from 2000-01-01T12:00:00 TAI.
     * <p>Beware, it is not {@link #J2000_EPOCH} since it is in TAI and not in TT.</p> */
    private final long epoch;
    /** Offset from the reference epoch in seconds. */
    private final double offset;
    /** Create an instance with a default value ({@link #J2000_EPOCH}).
     */
    public AbsoluteDate() {
        epoch  = J2000_EPOCH.epoch;
        offset = J2000_EPOCH.offset;
    }

using a variant of option (3) with fraction of days instead of seconds in days.

I am not suggesting the current representation should be changed but that the necessary accuracy and efficiency should be achieved and for this, it might be easier/more efficient (or not) to use a different representation.

[helgee]

Thanks for dragging all the skeletons out of the closet 😅 I am on vacation until Thursday so I won’t be able to give this a proper look until Friday.

[bgodard]

I probably said something wrong concerning Orekit: their offset from reference epoch is not constrained to be within a day.

I think for now the easiest way to fix this is simply replace:

function (+)(ep::Epoch{S,T1}, p::Period{U,T2}) where {T1,T2,S,U<:TimeUnit}
    delta = get(days(p))
    if delta >= oneunit(T2)
        ep1 = Epoch{S}(julian1(ep) + delta, julian2(ep))
    else
        ep1 = Epoch{S}(julian1(ep), julian2(ep) + delta)
    end
end

function (-)(ep::Epoch{S,T1}, p::Period{U,T2}) where {T1,T2,S,U<:TimeUnit}
    delta = get(days(p))
    if delta >= oneunit(T2)
        ep1 = Epoch{S}(julian1(ep) - delta, julian2(ep))
    else
        ep1 = Epoch{S}(julian1(ep), julian2(ep) - delta)
    end
end

by:

function (+)(ep::Epoch{S,T1}, p::Period{U,T2}) where {T1,T2,S,U<:TimeUnit}
    delta = get(days(p))
    Epoch{S}(julian1(ep), julian2(ep) + delta)
end

function (-)(ep::Epoch{S,T1}, p::Period{U,T2}) where {T1,T2,S,U<:TimeUnit}
    delta = get(days(p))
     Epoch{S}(julian1(ep), julian2(ep) - delta)
end

Since long duration are not accurate anyway, it is not a problem degrading the accuracy of an epoch when adding a large duration such as several years.

However this is still not perfect:

t1 = TDBEpoch("2018-06-01T00:00:00.0")
years = 365.25 * days
microsecond = 1e-6 * seconds 
t2 = t1+100*years # now both jd1 and jd2 are big so the absolute accuracy for future ops is very poor
t3 = t2+1*microsecond 
t3-t2 # inaccurate (0.0?)

This can be solved if the epoch is always re-normalised to eg 0<=jd2<1.0 but of course this requires more CPU time.

[bgodard]

I had typos in the 2 previous posts where I was saying millisecond but was actually talking about microseconds. This is now fixed.

[helgee]

We might want to have a look at the approach taken by the unpublished Chrono.jl. Which takes great care to preserve precision.

[helgee]

All in all, I am in favor of approach 3 (Orekit-style). This would preserve some genericity, i.e. being able to use arbitrary floating point types for the offset from the reference epoch, and keep the implementation simple as well. I looked around a bit for prior art and found this re-normalization routine in AstroPy: https://github.com/astropy/astropy/blob/master/astropy/time/utils.py My naive Julia port benchmarks like this:

BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     21.743 ns (0.00% GC)
  median time:      22.154 ns (0.00% GC)
  mean time:        23.041 ns (0.00% GC)
  maximum time:     196.103 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000

[helgee]

Fixed by #43