Thanks for tracking down this issue @demorest!

Could you provide your test case with the correct transformation difference?

One thing I want to do is understand why our test case was insufficient. We verify that the transforms match what is shown in the SOFA documentation. The test case is essentially running through all the transforms. Starting from a TT time of '2006-01-15 21:25:42.684000000' we get:

In [18]: t_tt.iso
Out[18]: '2006-01-15 21:25:42.684000000'

In [19]: t_tt.tdb.iso
Out[19]: '2006-01-15 21:25:42.683798890'

In [20]: t_tt.lon, t_tt.lat
Out[20]: (<Longitude -155.933222 deg>, <Latitude 19.48125 deg>)

This matches exactly the SOFA result (page 25 of http://www.iausofa.org/sofa_ts_c.pdf), so maybe there is a bug in SOFA / ERFA?

Yes, it looks like the SOFA example that you copied your test case from has the same bug: iauGd2gc() returns meters, while iauDtdb() wants km. I don't have any reason so far to think there is a bug in the actual SOFA/ERFA code, just this error in the example usage in the manual.

I'll see if I can come up with a test case you can use for this. However, it's not completely trivial to separate out just this calculation from all the other things TEMPO2 does. Also I'm not yet sure how closely to expect the two TT-TDB versions to agree, as they are using different methods for some parts of the calculation.

@taldcroft: I just filed a bug report to the SOFA folks (which I guess includes me!) about their example code.

OK, thanks for filing the documentation bug report @scottransom. I should be able to patch this here pretty soon. It would be really useful to get some independent transformation test values.

Looks like Pat Wallace has already fixed the SOFA documentation, @taldcroft. Here are the new correct values for the results of the SOFA tests on p25 of the document:

TDB 2006/01/15 21:25:42.684372
TCB 2006/01/15 21:25:56.893951

@scottransom @demorest - patch now attached, please have a look.

cc: @mhvk - same issue will need fixing in #1928 for 0.4.

Yikes, that is pretty bad. Guess we should take more seriously what we teach our students - check everything independently! Will rebase #1928 on this once you merge.

Hey Marten: the cool thing is that this was found with an incredibly independent test. Paul was checking the timing residuals as calculated with PINT for a high-precision MSP, and that showed a big annual term. So a "pulsar time standard" is what uncovered the error! That's pretty cool, IMO....

@scottransom, @demorest - That's great indeed! More generally, great that this is being helpful. Do also not hesitate to let us know missing stuff. E.g., some of the calculations now made in spice and pint should probably be covered by astropy (barycentric correction being a prime example). Feel free to open issues on those, ideally with links to code that implements it.

It will certainly mean a lot if astropy.time gets the millisecond pulsar stamp of approval!