juliaastro / astrotime.jl Goto Github PK

This package helpfully defines GPS_EPOCH here


However, I think it's incorrect to define it in TT. It should be 1980-01-06 UTC according to wikipedia and also section 11.4.3.3 of the manual from the US Coast Guard. So the definition should be equivalent to from_utc(1980,1,6) or TAIEpoch(1980, 1, 6, 0, 0, 19.0).

Based on @bgodard's comments on #45 it seems we still do not have an optimal design. The current implementation probably works for anybody not doing interplanetary flight dynamics (modulo bugs). Addressing this topic will probably entail another redesign and this time I want to get it right. Thus, I would like to collect some user stories and test cases here.

From #45

julia> t0=TAIEpoch(2000,1,1,12,00,32.0)
2000-01-01T12:00:32.000 TAI

julia> t1=UTCEpoch(2000,1,1,12,00,32.0)
2000-01-01T12:00:32.000 UTC

julia> t2=UTCEpoch(0.0,t0)
2000-01-01T12:00:32.000 UTC

julia> t2-t1
0.0 seconds

julia> t2==t1
false

julia> t2>t1
false

julia> t2<t1
false

julia> t2,t1
(2000-01-01T12:00:32.000 UTC, 2000-01-01T12:00:32.000 UTC)

I have also decided to restrict the input to TAIEpochs when overriding S-TAI.

This is not convenient at all because transformation between TDB and spacecraft proper time in a deep space trajectory has no meaningful relation to TAI. User would have to invent a fake intermediate TAI transformation.

julia> t2=TDBEpoch(2000,1,1)
2000-01-01T00:00:00.000 TDB

julia> t1=TDBEpoch(2000,1,1,12)
2000-01-01T12:00:00.000 TDB

julia> t2-t1
-43199.99998551589 seconds

But I would like to get the number of TDB seconds not TAI (TAI difference is not very useful for solar system dynamics though this is a good approximation for most applications).

Maybe we could have difference of Epoch in same timescale return seconds in that timescale.
For mixed timescales Epoch difference, raise an Error or convert the second epoch argument to the timescale of the first epoch argument (could be explicitly enabled by "using UnsafeEpochConversions")

At the moment I see it as a big problem that all timescales are directly related to TAI. I would like to have possibility for epochs which do not have any valid default conversion to TAI (timescale graph is disconnected), but possibly to other custom timescales.

The function is straightforward to implement. Leapseconds makes it quite easier.

Regarding this test
leap = UTCEpoch(2016, 12, 31, 23, 59, 60). I'm a bit unsure, here's what I think, should I subtract the drift from the seconds part because 60 is going out of range.
@helgee

#63 (review)

E.g.

# Three days from J2000
AstroTime.UTCEpoch(3days)

# Twenty seconds from J2000
AstroTime.UTCEpoch(20seconds)

When installing the AstroTime, a warning is flagged saying that downloaded files are corrupted and to execute EarthOrientation.update(force = true). This still does not fix the issue as this error appears:

julia> EarthOrientation.update(force = true)
[ Info: Downloading file 'finals.csv' from 'https://datacenter.iers.org/data/csv/finals.all.csv' with cURL.
[ Info: Downloading file 'finals2000A.csv' from 'https://datacenter.iers.org/data/csv/finals2000A.all.csv' with cURL.
ERROR: BoundsError: attempt to access 19011-element Vector{Float64} at index [0]
Stacktrace:
[1] getindex(A::Vector{Float64}, i1::Int64)
@ Base ./essentials.jl:13
[2] EOParams(iau1980file::String, iau2000file::String)
@ EarthOrientation ~/.julia/packages/EarthOrientation/1ITdn/src/EarthOrientation.jl:161
[3] push!
@ ~/.julia/packages/OptionalData/sGxJ1/src/OptionalData.jl:98 [inlined]
[4] push!
@ ~/.julia/packages/OptionalData/sGxJ1/src/OptionalData.jl:101 [inlined]
[5] update(; force::Bool)
@ EarthOrientation ~/.julia/packages/EarthOrientation/1ITdn/src/EarthOrientation.jl:62
[6] top-level scope
@ REPL[8]:1

Any help on this would be appreciated as it means that we cannot use the AstroTime package for satellite propagation.

To be tackled when the porting from ERFA is done, e.g.

function shift(jd1, jd2, offset)
    if jd1 > jd2
        jd1 = jd1
        jd2 += offset
    else
        jd1 += offset
        jd2 = jd2
    end
    jd1, jd2
end

# TAI -> TT
taitt(jd1, jd2) = shift(jd1, jd2, OFFSET_TT_TAI/SECONDS_PER_DAY)
# TT -> TAI
tttai(jd1, jd2) = shift(jd1, jd2, -OFFSET_TT_TAI/SECONDS_PER_DAY)

The following expression generates a vector:

julia> (1:3)*seconds
3-element Vector{AstroPeriod{AstroTime.Periods.Second, Float64}}:
1.0 seconds
2.0 seconds
3.0 seconds

I argue that it should create a new range:

julia> (1:3)*seconds
1 seconds: 1 seconds: 3 seconds

The reason is that the expression "(1:typemax(Int32))*seconds" is inefficient in terms of time and memory, because it generates a huge vector, instead of a range type.

If you want to generate a vector, then use the dot notation, i.e., "(1:3).*seconds".

This change just requires the addition of two functions, one each for the implied and explicit versions of step.

Also, the expression "(1:typemax(Int64))*seconds" generates an incorrect result, because the maximum value is converted to a float and then back to an integer. This corner case should be handled properly. Int32 works as expected.

If I have 2 observers on the surface of the Earth take some measurements at a given TT = tt0 and I want to model a combined observation in the barycentric reference system, I can compute TDB1 (for observer 1) and TDB2 (for observer 2) using an accurate position dependent TDB-TT difference model.

tdb1 = TDBEpoch(tai_offset_observer1,tt0)
tdb2 = TDBEpoch(tai_offset_observer2,tt0)

but then:

julia> tdb2-tdb1
0.0 seconds

julia> tdb2==tdb1
true

but both of these are not correct. The assumption that time can be compared independently of timescale is incorrect (simultaneity is a relative concept).

CC #65

The following code currently fails because AstroPeriod is constructed implicitly using Float64

julia> typemax(Int64)*seconds
ERROR: InexactError: trunc(Int64, 9.223372036854776e18)
Stacktrace:
 [1] trunc
   @ ./float.jl:716 [inlined]
 [2] floor
   @ ./float.jl:294 [inlined]
 [3] AstroPeriod
   @ ~/dev/julia/AstroTime/src/Periods.jl:97 [inlined]
 [4] *(dt::Int64, unit::AstroTime.Periods.Second)
   @ AstroTime.Periods ~/dev/julia/AstroTime/src/Periods.jl:152
 [5] top-level scope
   @ REPL[27]:1

I'm sure there's probably a simple fix using oftype, I don't have time at the moment to dig into it though.

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

The ports of taiutc and utctai are too simplistic and give an incorrect result during a leap second, see https://github.com/liberfa/erfa/blob/master/src/utctai.c, https://github.com/liberfa/erfa/blob/master/src/taiutc.c.

I have rolled the changes back and added two broken tests that should pass once the issue is fixed: https://github.com/JuliaAstro/AstronomicalTime.jl/blob/master/test/runtests.jl#L181

The following range expression produces an incorrect value when using TT time.

TTEpoch("2021-05-09T07:01:10"):TTEpoch("2021-05-09T08:02:00")
2021-05-09T07:01:10.000 TT:86400.0 seconds:2021-05-09T07:01:10.000 TT

The conversion between TDB and TT in AstroTime.jl taken from ERFA/SOFA is time consuming because it is computing a large series.

The accuracy is probably of the order of a few nanoseconds but for what purpose: the difference between TDB and TT is a quantity which is different for each observers and dependent on their trajectory (their position and velocities in GCRS). The quantity can vary by a few microseconds from one ground station at the surface of the Earth to another.

In the routine dtdb(jd1, jd2, ut, elong, u, v), arguments ut, elong, u and v are only used to compute the observer position correction (observer position is actually defined by elong, u and v but the correction needs additionally the value of ut) neglecting radial velocity correction (which is only important for observer not on the surface such as an Earth satellite), but they are set to zeros (correction disabled) in the rest of the code because there is no observer trajectory information. Thus the quantity computed is the TT minus TDB for an observer at the geocenter.

While it is good to have the routine dtdb available for high accuracy computation, for most use cases where the observer location is unspecified a faster and lower accuracy routine would probably be better.

Additionally the difference between TDB and TT can be computed faster and more accurately using state of the art numerically integrated geocenter time ephemeris. Example:

using BenchmarkTools

using AstroTime
using CALCEPH

astrotime_TTminusTDB(tdb1,tdb2)= - AstroTime.Epochs.dtdb(tdb1,tdb2,0.0,0.0,0.0,0.0)

eph1=Ephem("inpop13c_TDB_m100_p100_tt.dat")
calceph_TTminusTDB(tdb1,tdb2)=compute(eph1,tdb1,tdb2,naifId.id[:ttmtdb],naifId.id[:timecenter],unitSec+useNaifId,0)[1]

eph2=Ephem("inpop13c_TDB_m100_p100_tt.dat")
prefetch(eph2)
calceph_TTminusTDB_prefetched(tdb1,tdb2)=compute(eph2,tdb1,tdb2,naifId.id[:ttmtdb],naifId.id[:timecenter],unitSec+useNaifId,0)[1]


tdb = 2458257.0:2458257.0+366*10

difference(tdb) = astrotime_TTminusTDB(tdb,0.0)-calceph_TTminusTDB(tdb,0.0)

using Plots
pyplot()

plot(tdb,difference.(tdb))
xlabel!("JD")
ylabel!("seconds")

tdb = 2414105.5:0.49:2488984.5

function test1(fun,tab)
   for x in tab
      fun(floor(x),x-floor(x))
   end
end


b1 = @benchmark test1(astrotime_TTminusTDB,$tdb)
display(b1)
println()
b2 = @benchmark test1(calceph_TTminusTDB,$tdb)
display(b2)
println()
b3 = @benchmark test1(calceph_TTminusTDB_prefetched,$tdb)
display(b3)
println()

Result:

BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     2.024 s (0.00% GC)
  median time:      2.026 s (0.00% GC)
  mean time:        2.026 s (0.00% GC)
  maximum time:     2.029 s (0.00% GC)
  --------------
  samples:          3
  evals/sample:     1

BenchmarkTools.Trial: 
  memory estimate:  30.31 MiB
  allocs estimate:  916890
  --------------
  minimum time:     45.561 ms (3.81% GC)
  median time:      46.358 ms (3.98% GC)
  mean time:        46.627 ms (5.28% GC)
  maximum time:     50.704 ms (3.63% GC)
  --------------
  samples:          108
  evals/sample:     1

BenchmarkTools.Trial: 
  memory estimate:  30.31 MiB
  allocs estimate:  916890
  --------------
  minimum time:     39.932 ms (4.32% GC)
  median time:      41.296 ms (4.55% GC)
  mean time:        41.379 ms (6.10% GC)
  maximum time:     46.952 ms (9.21% GC)
  --------------
  samples:          121
  evals/sample:     1

The computed values differ by only a few nanoseconds but the runtime is much better.
If AstroTime could target microsecond TDB-TT conversion accuracy (better accuracy doesn't make sense without a location specification) while being much faster than interpolating the time ephemeris that would be very useful.

EDIT: I mixed up nanoseconds and picoseconds. This is fixed now!

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.

Replace the calls to ERFA by pure Julia functions:

• eraD2dtf
• eraDat #31
• eraDtdb #14
• eraDtf2d #36
• eraTaitt #4
• eraTaiut1 #8
• eraTaiutc #13 #39
• eraTcbtdb #25
• eraTcgtt #9
• eraTdbtcb #24
• eraTdbtt #14
• eraTttai #6
• eraTttcg #10
• eraTttdb #20
• eraTtut1 #21
• eraUt1tai #7
• eraUt1tt #22
• eraUt1utc #40
• eraUtctai #15 #37
• eraUtcut1 #38

Printing relies on Dates functionality which does not know about leap seconds.

Hi,
I'm too dumb to use git or I would have created an actual pull request (if that's the proper term).

AstroTime.now() returns the computer's local time as it were UTC.
In the file Epochs.jl, line 352ish, I've changed

• now() = UTCEpoch(Dates.now())
  into

• now() = UTCEpoch(Dates.now(Dates.UTC))
  And that change makes AstroTime.now() return the UTC time independently of your local time.

We currently download an LSK kernel with leap seconds data once and the parse it on every import. This is not incredibly useful since the download URL will change everytime NAIF publishes a new kernel. In consequence we will need to publish a new release everytime a new leap second is introduced which means that we could also (semi-)hard code leap seconds instead.

I propose to use a Julia script to generate a file with the required constants from an LSK kernel similar to the approach taken by Chrono.jl. This has the added benefit that AstroTime then does not need to depend on RemoteFiles.jl and OptionalData.jl anymore.

Hi, I am a new user of AstroTime.jl, and I found this package very useful. I encountered an issue as below. Note that a leap second was inserted before 2017-1-1 00:00:00 UTC, but this is the case around 2016-12-31 00:00:00 UTC. The same happens for 2015-6-30.

julia> ep = UTCEpoch(2016, 12, 31, 0, 0, 0.0)
2016-12-30T23:59:59.000 UTC

julia> ep = UTCEpoch(2016, 12, 31, 0, 0, 0.1)
Error showing value of type Epoch{CoordinatedUniversalTime,Float64}:
ERROR: ArgumentError: Seconds are out of range
Stacktrace:
(followed by a long message)

julia> ep = UTCEpoch(2016, 12, 31, 0, 1, 0.0)
2016-12-31T00:00:60.000 UTC

Thanks,
Toshi

Issue JuliaLang/julia#21905 contains a request for a date of year format, almost identical to the 2nd format in date_conv.pro at IDL AstroLib.
Such a feature seems highly suitable to be present in AstronomicalTime.jl

Hi

I was wondering if you were planning on bumping the version of RemoteFiles to 0.4. I am having some dependency issues with RemoteFiles and HTTP that require this.

I've done the port and it's working and providing the julian of the particular day. But the DateTime is failing to parse the julian date as seconds is going outside the range of DateTime.

I know d2dtf is required to get the julian dates to Gregorian format but still, DateTime would fail.

since the Date type only resolves to the precision of a single date (i.e. no hours, minutes, or seconds), normal considerations for time zones, daylight savings/summer time, and leap seconds are unnecessary and avoided.

in the Julia docs.

Do I have to design a new show function?
Any suggestions

If I define a custom time scale inside a package, conversion will not work:

module FooBar

using AstroTime

@timescale GPSTime TAI
AstroTime.Epochs.getoffset(
    ::GPSTimeScale, 
    ::InternationalAtomicTime, 
    second, 
    fraction) = 19
AstroTime.Epochs.getoffset(
    ::InternationalAtomicTime, 
    ::GPSTimeScale, 
    second, 
    fraction) = -19

end

Note that when copying this into the REPL, it will actually work. It does not work when it is a real package (not a module inside Main).
I get the following error:

julia> using FooBar, AstroTime
julia> TAIEpoch(FooBar.GPSTimeEpoch(2021,5,31))
ERROR: No conversion path between 'GPSTime' and 'TAI' available.
Stacktrace:
 [1] apply_offset
   @ ~/.julia/packages/AstroTime/944cN/src/Epochs/offsets.jl:154 [inlined]
 [2] (TAIEpoch{T} where T)(ep::Epoch{IGS.GPSTimeScale, Float64})
   @ AstroTime.Epochs ~/.julia/packages/AstroTime/944cN/src/Epochs/offsets.jl:39
 [3] top-level scope
   @ REPL[3]:1