-
Notifications
You must be signed in to change notification settings - Fork 101
Clock Corrections and Timescales in PINT
Clock corrections and timescale conversions in PINT are handled by a combination of astropy Time objects (and its built in timescale transformation capabilities) and clock correction files stored in pint/datafiles
.
Here we walk through the transformations and how they happen for a typical TOA read from a .tim file.
TOA times read from TOA (.tim) files are assumed to be in the UTC timescale, referenced to an observatory clock (often a GPS-disciplined rubidium clock). This timescale is denoted UTC(obs), where obs is the name of the observatory [2]. The MJD representation used in these files is the pulsar_mjd
, as described below. When reading from a file, precision is maintained by reading the integer and fractional parts into separate values (arg1
and arg2
). Then the Time
object is instantiated like this:
time.Time(arg1, arg2, scale=scale,
location=site.earth_location,
format='pulsar_mjd', precision=9)
People from the precise timing community know well that using MJD format for UTC times is a bad idea. There is not a unique way to assign MJDs to times during days with leap seconds, and MJD1 - MJD2 does not correctly give the time interval between two times, because of possible leap seconds between MJD1 and MJD2 (this has been the source of incorrect results in the literature). Nevertheless, MJDs are commonly used for UTC times in many places. In pulsar timing, the MJD times in TOA (.tim) files are in the UTC time scale. PINT follows TEMPO and Tempo2 in defining a pulsar MJD (in pint called pulsar_mjd
) format, in which the integer part is the normal integer MJD and the fractional part is seconds_of_day/86400 [1]. The means that MJDs 'tick' at a constant rate, but there is no representation for a time during a leap second, so there is no way to have a TOA during that time.
This example illustrates the issue. On a leap second day (MJD 56108), the MJD 56108.25 corresponds to exactly 6:00 AM in the pulsar_mjd
format, but is 0.25 seconds later in the mjd
format that astropy uses. On any other day, both will be exactly 6:00 AM.
In [5]: from astropy.time import Time
In [6]: import pint.pulsar_mjd
In [7]: t1 = Time(56108.25,scale='utc',format='mjd',precision=9)
In [8]: t2 = Time(56108.25, scale='utc',format='pulsar_mjd',precision=9)
In [9]: t1.isot
Out[9]: '2012-06-30T06:00:00.250000000'
In [10]: t2.isot
Out[10]: '2012-06-30T06:00:00.000000000'
WARNING: I (@paulray) think we need to look through PINT for Time inits with
format='mjd'
and make sure they are correct. Those may break on leap second days. They should probably all beformat='pulsar_mjd'
Given a UTC(obs), PINT applies clock corrections to convert to UTC(GPS). This can be done using either TEMPO or Tempo2 format clock files. The default is to use the set of TEMPO clock files distributed with PINT in pint/datafiles
. These correct UTC as measured by the observatory clock to the UTC(GPS) timescale from GPS.
NOTE: I (@paulray) really don't have a good understanding of the TEMPO clock correction files. What is the deal with the two columns? If you apply the clock correction, exactly what timescale are you on? Is it always UTC(GPS)? Could someone add a better explanation of them here? In most comments and docs, it says things like "apply the clock correction" without clearly defining what timescale the correction goes from and to. Tempo2 is generally clearer in its clock correction strategy.
A further correction can be applied to convert UTC(GPS) to "UTC" (which I think means UTC(USNO), but this should be confirmed). This is done by default using the Tempo2-format gps2utc.clk file in pint/datafiles
. Those corrections are derived from BIPM Circular T. Whether this correction is applied is controlled by the include_gps
argument to pint.observatory.observatory.get_observatory()
.
UTC is converted to TAI by adding an integer number of seconds. The number changes to account for leap seconds, so that TAI is a continuous timescale. The current value (since January 2017) of TAI-UTC is 37.0 seconds. In PINT, this conversion is handled by astropy, which in turn uses ERFA (the free version of IAU SOFA), which has the leap seconds hard coded. Versions of astropy 1.3 or later include the 2017 January 1 leap second. Note the 37 second difference in this example:
In [10]: t = Time("2017-01-15T06:00:00.000",scale='utc',format='isot',precision=9)
In [11]: t.utc.isot
Out[11]: '2017-01-15T06:00:00.000000000'
In [12]: t.tai.isot
Out[12]: '2017-01-15T06:00:37.000000000'
TT (also known as TDT) ticks at the same rate as TT and UTC, but for reasons of continuity has an offset. TT has a unit of duration 86400 SI seconds on the geoid and is the independent argument of apparent geocentric ephemerides.
The most common realization of TT is TT(TAI), which is defined as:
TT(TAI) = TAI + 32.184 seconds
In [13]: t.tt.isot
Out[13]: '2017-01-15T06:01:09.184000000'
PINT can also use TT(BIPM), which is a realization of TT published by the BIPM on their web page. In PINT, this clock correction is read from the Tempo2-style clock file pint/datafiles/tai2tt_bipm2015.clk
. Whether this correction is enabled is controlled by the include_bipm
argument to pint.observatory.observatory.get_observatory()
.
WARNING: I think the current implementation of TT(BIPM) corrections has issues. It is applied in
TopoObs.clock_correction
. These corrections are supposed to be part of the UTC(obs) -> UTC transformation. Currently, the SpecialObservatory classes all have clock_correction = 0.0. If PINT is trying to use TT(BIPM), then probably observatories that have TT TOAs (like Geocenter, NICER, Fermi) want to have the BIPM corrections applied. I guess that for Barycenter TOAs, they can be assumed to be TDB times regardless of whether PINT is set to TT(TAI) or TT(BIPM) since the TT -> TDB correction has already been applied.
TDB is the independent argument of barycentric solar system ephemerides. TDB varies from TT only by periodic variations. The difference TDB−TT is quasi-periodic, dominated by an annual term of amplitude 1.7 ms.
The conversion from TT to TDB is done by astropy using the delta_tdb_tt()
attribute. The way this works is fully described in the SOFA Documentation (PDF).
From the SOFA Docs:
The SOFA routine DTDB provides a detailed model (nearly 800 terms) for TDB−TT, accurate to a few nanoseconds in the present era when compared with the latest work. Its arguments are TDB (TT can be used instead without materially affecting the results), UT1, and the observer’s coordinates, expressed as longitude and the distances from the rotation axis and equator. If the application can tolerate errors of up to 2 μs, zeroes can be used for the final four arguments, giving a geocentric result. The example in the next Section demonstrates use of the DTDB routine, including what is involved in transforming the observer’s terrestrial coordinates into the required form.
I am not actually sure whether the SOFA routine is identical to either Irwin & Fukushima (1999) or Fairhead & Bretagnon (1990). Note that Tempo2 selects between these two published models by setting TIMEEPH to either IF99 or FB90.
It is this DTDB transformation that requests the PINT Time objects representing TOAs have their location
property set.
Following the examples above, for this time that has no location
set, we get this (an offset of 361 µs from TT):
In [14]: t.tdb.isot
Out[14]: '2017-01-15T06:01:09.184361443'
Setting the location
to a spot on the Earth's surface (e.g. longitude +36 deg, latitude +21 deg) changes the TDB value by 0.68 µs.
Question: Does setting the
location
to the spacecraft location when using spacecraft photon TOAs do the 'right' thing?
[1] Note that some other codes (notably IAU SOFA and thus astropy) use seconds_of_day/86400 for non-leap-second days and seconds_of_day/86401 for days with a leap second. This makes MJDs 'tick' at a different rate on leap second days, but allows all times to have a representation. PINT code that writes TOA files needs to ensure it is using the pulsar_mjd
format for compatibility with the PINT, TEMPO and Tempo2 tim file readers. There is discussion of this issue here
[2] Note that this is not sufficient to distinguish between multiple instruments or backends at an observatory. In the past, this has been dealt with by making two different observatory names (e.g. ncy and ncyobs), which have the same position, but different clock correction files. PINT really should do this a better way.
[3] PINT really needs to support TCB as well as TDB, but this is not currently implemented.