Why does FK4 <-> Galactic disagree between Astropy/Kapteyn and PyAST/PySLALIB

[astrofrog]

@timj and @dsberry: Astropy and Kapteyn currently agree well with PyAST/PySLALIB for the FK5 <-> Galactic conversion, with differences around 1mas. On the other hand, for FK4 <-> Galactic the differences can be up to 0.2 arcseconds which is quite large. This is a bit weird since the FK4 <-> FK5 conversion agrees well between all packages:

http://www.astropy.org/coordinates-benchmark/summary_matrix.html#astropy

I think that both Astropy and Kapteyn do:

FK4 -> FK4 (no e terms) -> Galactic

Do PyAST and PySLALIB do something different? What would cause such a large difference? Why do we agree for FK5 -> Galactic but not FK4 -> Galactic? (the relevant scripts are in e.g. tools/kapteyn/convert.py, tools/pyast/convert.py, etc.

[astrofrog]

Just to make sure we are all on the same page, here is what the benchmarks should be testing:

https://github.com/astropy/coordinates-benchmark#the-benchmarks

[dsberry]

If you ensure the epoch is the same in all systems (which I see you do), then pyast will do fk4->fk5, followed by fk5->galactic. So since those two conversion both agree with kapteyn/astropy, I'm at a bit of a loss to know why the total fk4->galactic conversions don't agree. I'll hopefully be able to look further into this tomorrow.

[astrofrog]

@dsberry - thanks for offering to look into this!

[dsberry]

The pyast and pyslalib tests implement FK4->galactic as FK4->FK5, then FK5-> galactic. The first part (FK4->FK5) depends on the epoch (which the coordinates tests set at Julian 2000), but the second part (FK5->galactic) does not. SLALIB contains a direct FK4 to galactic conversion routine (sla_eg50) which does not go via FK5 and does not allow the epoch to be specified. If I change the pyslalib tests to use this direct route rather than the FK5 route, then pyslalib agrees with kapteyn and astropy. So the question is, which is the best approach?

I've talked a bit to Pat Wallace about this. I gather that the position of the galactic north pole is very nearly fixed (i.e. does not vary with time) in FK5 coordinates. It's only "very nearly" because FK5 is not quite an inertial frame. But that means that the galactic north pole does move (significantly) in FK4 coordinates, because FK4 moves with respect to FK5. It seems from Pat that the direct FK4->galactic routine in slalib (sla_eg50) gets round this by simply assuming that the input FK4 position is epoch (and equinox) B1950 and using the known position of the galactic pole in FK4 at that time. That's why sla_eg50 has no epoch parameter. And so I assume therefore that the astropy/kapteyn fk4->galactic conversion also assumes epoch B1950 inputs, since their results agree with sla_eg50. But in the case of this test suite, the supplied positions are specified as epoch J2000. The direct route used by the astropy and kapteyn tests will ignore this J2000 epoch and assume the values are actually epoch B1950. But the pyslalib and pyast tests honour the J2000 epoch and so give different results.

That's my best stab at explaining the differences at the moment. If you modify the astropy/kapteyn tests to first transform the epoch J2000 test FK4 positions to the epoch B1950 equivalents, before then converting to galactic, I suspect the differences would go away.

[astrofrog]

@dsberry - thanks for the clarification! I'll keep this issue open to remind myself to look into it. I think @eteq and I will just need to go through exactly how the conversion is being done in astropy to see how we can improve things.

[eteq]

@dsberry - thanks for looking into this! One quick bit I'm not clear on:

I gather that the position of the galactic north pole is very nearly fixed (i.e. does not vary with time) in FK5 coordinates.
But that means that the galactic north pole does move (significantly) in FK4 coordinates, because FK4 moves with respect to FK5.

Why does he say the NGP is fixed to FK5? Our reading at the time we implemented this in astropy was that the most authoritative definition of Galactic coordinates is the Blaauw et al. 1958 IAU definition, which is based on a specified FK4 NGP. So by that reading it should be fixed in FK4 rather than FK5. Was there some update to that definition that we missed?

[dsberry]

My understanding is that the NGP was indeed defined in FK4, but at a specified epoch (B1950). The point on the sky that corresponds to FK4 RA=12:49:02 Dec=27:22.7 at equinox=B1950 and epoch=B1950 will have different FK4 RA and Dec values at some other epoch due to the rotation of the FK4 system with time.

[dsberry]

See also section 7) of Murray C.A., 1989, A&A, 218, 325. These are the matrices used by the slalib routines sla_eg50 (FK4 to galactic) and sla_eqgal (FK5 to galactic).

Also see the introduction to Jia-Cheng Liu et al (http://arxiv.org/pdf/1010.3773.pdf) in which the authors propose a redefinition of galactic coordinates so that it is fixed with respect to ICRS rather than FK5.

[eteq]

Ah, I think I get what you mean now - you mean it moves insofar as FK4 has a signficant obstime/epoch component. So as long as we compare at a fixed obstime/epoch, it should still be the same.

For the record, I am 100% with Liu et al. in that it really ought to be tied to ICRS... But it's not clear there's much point to it now given how there's already >50 yrs of literature referencing the "old" Galactic system. My take (which admittedly may be due to my field) is that it's safer to just work in ICRS or another better-defined system, and use Galactic for convenience rather than as a real "working" coordinate system.

[astrofrog]

@eteq - in general I do think that since the old system is ill-defined, there's still no reason why one can't redefine it relative to ICRS such that it 'best fits' the old system. Of course, we can't just arbitrarily re-define it though, so it's not up to us...

[cdeil]

So now that astropy/astropy#3107 is merged, does that mean the definition of Galactic coordinates in Astropy is fixed and OK?
Or is there a need to have any of the different definitions discussed here implemented and available to users as well?

In any case, it would be great if one of the experts here could write a few sentences about "Galactic" coordinates for the Astropy docs. If I'm not mistaken there currently isn't any info:
http://astropy.readthedocs.org/en/latest/search.html?q=galactic&check_keywords=yes&area=default
I don't think much text is needed, just:

• Give a definition and reference which Galactic coordinates Astropy has implemented.
• Maybe mention a few other definitions / codes for Galactic coordinates and roughly by how much they differ from the ones implemented in Astropy (milli-arcsec? arcsec? arcmin?). E.g. in practice it matters to astronomers whether Astropy gives the same results as the NED coordinate calculater and / or matches SIMBAD?

[astrofrog]

@cdeil - regarding your first question, this is not yet resolved because as we were discussing above we don't always get the same as PyAST.