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 files accessed as needed via the IPTA Clock Repository.
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'
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 files from the IPTA Clock Repository. 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 from the IPTA Clock Repository. Those corrections are derived from BIPM Circular T. Whether this correction is applied is controlled by the apply_gps2utc property of the Observatory object, and is defined in observatories.json or at observatory creation. The file can be found here pint.config.runtimefile("observatories.json"). (Note: For advanced users, this can be overridden by specifying the apply_gps2utc kwarg 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 IPTA Clock Repository. Whether this correction is enabled is controlled by the include_bipm argument to pint.toa.get_TOAs() and other functions.
Note: The current difference between TT(TAI) and TT(BIPM) is about 27 microseconds.
NOTE: For TOAs that are already in TDB (site is "bat" or "@"), no GPS or BIPM corrections are 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).
The clock correction flow in astropy is shown in the following figure:
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.
The SOFA code for dtdb.c, uses the Fairhead & Bretagnon (1990) model for computing the geocentric conversion from TT to TDB. Note that Tempo2 selects between this and Irwin & Fukushima (1999) by setting TIMEEPH to either FB90 or IF99.
It is this DTDB transformation that requests the PINT Time objects representing TOAs have their location property set.
The location-dependent part of the transformation is also describe in the SOFA code:
The topocentric part of the model is from Moyer (1981) and Murray (1983), with fundamental arguments adapted from Simon et al. 1994. It is an approximation to the expression ( v / c ) . ( r / c ), where v is the barycentric velocity of the Earth, r is the geocentric position of the observer and c is the speed of light. By supplying zeroes for [the position], the topocentric part of the model can be nullified, and the function will return the Fairhead & Bretagnon result alone.
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
locationto the spacecraft location when using spacecraft photon TOAs do the 'right' thing? For example, does setting the location implicitly assume that the observatory is co-rotating with the Earth and use that to compute the velocity used in the time dilation function? Answer: Yes, I believe that is the right thing to do. That term "expresses the fact that events simultaneous in the barycentric system usually are not simultaneous in the geocentric system" (see Eqn 13 in https://ui.adsabs.harvard.edu/abs/1991CeMDA..52..355S/abstract).
[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.