Determine source of difference between Kapteyn and PyAST for FK4/FK5

[astrofrog]

For the FK4 <-> FK5 conversion, there is a disagreement of 4mas between PyAST and Kapteyn. Pyephem is worse, but for now I'm interested in understanding the difference between PyAST and Kapteyn, because I can locally get astropy to agree with kapteyn to machine precision (almost) but it therefore disagrees with PyAST by ~4mas.

I think both PyAST and Kapteyn are correctly assuming the epoch of observation is J2000. What else could be the source of the difference? @timj and @dsberry - do you have any ideas? Kapteyn is including the e-terms for FK4 - is that also the case for PyAST?

I wanted to check if it's to do with the calculation of the e-terms - how do we not include the e-terms in PyAST if they are included by default?

[timj]

I think @dsberry mentioned that the difference is probably in the nutation/precession models. AST uses PAL which uses the most recent IAU models from SOFA.

[astrofrog]

@timj @dsberry - I just noticed that the disagreement appears to become larger towards the poles - is that an effect of the precession/nutation models? For example:

import numpy as np
import starlink.Ast as Ast
from astropy import units as u
from astropy import coordinates as coord
from astropy.time import Time
from kapteyn import celestial

frame_fk4 = Ast.SkyFrame('System=FK4-NO-E,Equinox=B1950,Epoch=J2000')
frame_fk5 = Ast.SkyFrame('System=FK5,Equinox=J2000')
frameset = frame_fk4.convert(frame_fk5)

ra_in = 4
dec_in = 89.9

coords = np.degrees(frameset.tran([[np.radians(ra_in)], [np.radians(dec_in)]]))
ra, dec = coords[0,0], coords[1,0]
print "PyAST:   %15.10f %15.10f" % (ra, dec)

c1 = coord.FK4NoETermCoordinates(ra_in, dec_in, 
                                 unit=(u.degree, u.degree),
                                 obstime=Time("J2000", scale='utc'),
                                 equinox=Time("B1950", scale='utc'))
c2 = c1.transform_to(coord.FK5Coordinates)
print "Astropy: %15.10f %15.10f" % (c2.ra.degrees, c2.dec.degrees)

before = (celestial.equatorial, celestial.fk4_no_e, "B1950", "J2000_OBS")
after = (celestial.equatorial, celestial.fk5, "J2000")

coords = celestial.sky2sky(before, after, ra_in, dec_in)
ra, dec = coords[0,0], coords[0,1]
print "Kapteyn: %15.10f %15.10f" % (ra, dec)

gives:

PyAST:    177.9056601073   89.8212115670
Astropy:  177.9052563909   89.8212116654
Kapteyn:  177.9052563869   89.8212116657

which is a difference of >1" (the difference is smaller away from the poles). Or is this more likely due to the way in which the transformation is done and numerical inaccuracies? The declinations agree better, it's mostly the RAs that are very different.

[astrofrog]

Ignore me for now, I think I was using the routine to find the distance between two coordinates wrongly. I will post an update once I've checked the values.

[astrofrog]

Ignore the issue I was reporting above - I was using the spherical separation routines from astropy wrongly (see astropy/astropy#562 for why). The differences are now less than 0.1" which is fine for the 0.2 release of Astropy. It looks like the difference are now dependent most on the observation epoch, so there must be a subtle difference in the algorithms.

[eteq]

So should this be closed @astrofrog , or do you want to leave it open for those smaller differences?

[astrofrog]

Let me update it with the latest version then I'll merge.

[astrofrog]

Ah wait, wrong issue ;-)

[astrofrog]

I haven't determined exactly what the difference is due to, so this should remain open.

[cdeil]

@sargas asked an (I think different but possibly related) question about Galactic to FK5 coordinate conversion differences here: http://mail.scipy.org/pipermail/astropy/2014-November/003510.html

To make it easier for us to discuss differences in the implementation I've put a copy of the latest Kapteyn here: https://github.com/cdeil/kapteyn-mirror/blob/master/kapteyn/celestial.py (and I've asked them if they could put a mirror of their repo on Github via email).

[cdeil]

Is this resolved or is there something left to understand / fix? (see https://github.com/astropy/coordinates-benchmark/issues/17#issuecomment-70984965)

[astrofrog]

The issue reported here is still true and the 4mas precision difference. What areas of astronomy would likely need the highest precision here?