[AstroPy] astropy.coordinates vs kapteyn.celestial Coordinate Transformations

[vogelaar at astro.rug.nl]

Hello,

There are two articles which support the decisions I made for the
conversion of positions from fk5 to galactic in module celestial of the
Kapteyn Package. The first article documents the rotation matrix that I
used and the second article (written after the current version of
celestial) documents the same matrix derived in a slightly different way.

1) Murray, C.A., 1989. The transformation of coordinates between systems
of B1950.0 and J2000.0 and the principal galactic axes referred to
J2000.0, Astron. Astrophys, 218, p.325-329,
2) Jia-Cheng Liu, Zi Zhu, Hong Zhang, Reconsidering the galactic coordinate
system, arXiv:1010.3773 [astro-ph.GA]

All the procedures implemented in celestial are documented in
https://www.astro.rug.nl/software/kapteyn/celestialbackground.html

The problem conversion mentioned by Joseph is:

from kapteyn import celestial
print celestial.sky2sky(celestial.galactic, celestial.fk5,
[292.03306305555554], [1.7592747222222223])
[[ 171.158163857372  -59.263193188529]]

which differs from the result by pyregion by a third of an arcsec.

Murray proposed that the axes of the IAU1958 galactic coordinate system
should be considered as absolute directions, unaffected by the E-term of
aberration.
This system is (in FITS terms) called FK4-NO-E.
This approach has been approved by Adriaan Blaauw (private communications,
2008), so in fact I use:

m1 = celestial.skymatrix(celestial.fk5, celestial.fk4_no_e)
print m1[0]
[[  9.999256794957e-01   1.118148323917e-02   4.859003772314e-03]
 [ -1.118148322047e-02   9.999374848933e-01  -2.717029374405e-05]
 [ -4.859003815359e-03  -2.716259471425e-05   9.999881946024e-01]]

This is the same matrix as Xo in Liu et al. in equation 9

m2 = celestial.skymatrix(celestial.fk4_no_e, celestial.galactic)

In [27]: print m2[0]
[[-0.066988739415 -0.872755765852 -0.483538914632]
 [ 0.492728466075 -0.45034695802   0.744584633283]
 [-0.867600811151 -0.188374601723  0.460199784784]]

This is the same matrix as Nb in Liu et al. in equation 9
and the same as in Slalib
(http://stsdas.stsci.edu/cgi-bin/gethelp.cgi?eg50.src)

print m2[0]*m1[0]
[[-0.054875539396 -0.873437104728 -0.48383499177 ]
 [ 0.494109453628 -0.444829594298  0.7469822487  ]
 [-0.867666135683 -0.198076389613  0.455983794521]]

This is the same matrix as Nj in Liu et al. in equation 9
and Murray equation 33.

The same matrix is constructed for:
m3 = celestial.skymatrix(celestial.fk5, celestial.galactic)

print m3[0]
[[-0.054875539396 -0.873437104728 -0.48383499177 ]
 [ 0.494109453628 -0.444829594298  0.7469822487  ]
 [-0.867666135683 -0.198076389613  0.455983794521]]


The inverse transformation is:
lonlat= celestial.sky2sky(celestial.galactic, celestial.fk5, [0], [0])

With this expression we find the galactic center in equatorial coordinates:
celestial.lon2hms(lonlat[0,0],prec=4)
'17h45m37.1991s'

celestial.lat2dms(lonlat[0,1],prec=4)
'-28d56m10.2207s'

which are the same values as in Liu et al. equation 11.

For the pole:
lonlat= celestial.sky2sky(celestial.galactic, celestial.fk5, [0], [90])

celestial.lon2hms(lonlat[0,0],prec=4)
'12h51m26.2755s'

celestial.lat2dms(lonlat[0,1],prec=4)
' 27d07m41.7043s'

which are the same as the position given in Murray section 7.1
and Liu equation 10.

So the transformation requires that the e-terms are removed.
and it uses a rotation matrix which is documented in several articles and
which guarantees that the galactic coordinate system is orthogonal.

The suggestion of Thomas Robitaille seems therefore worthwhile to discuss.

Martin




> Hello,
>
> I've been porting pyregion to use astropy instead of kapteyn, and tests
> with coordinate system conversions are slightly off.
>
> I think I've narrowed down the problem to my expectation that this should
> be nearly zero:
>
> In [21]: from astropy.coordinates import SkyCoord
>
> In [22]: from kapteyn import celestial
>
> In [23]: a = SkyCoord('292.03306305555554d 1.7592747222222223d',
> frame='galactic').transform_to('fk5'); print(a)
> <SkyCoord (FK5: equinox=J2000.000): ra=171.158093022 deg,
> dec=-59.2630875829 deg>
>
> In [24]: celestial.sky2sky(celestial.galactic, celestial.fk5,
> [292.03306305555554], [1.7592747222222223])
> Out[24]: matrix([[ 171.15816386,  -59.26319319]])
>
> In [25]: SkyCoord('171.15816386d -59.26319319d',
> frame='fk5').separation(a).to('arcsecond')
> Out[25]: <Angle 0.4019071919711007 arcsec>
>
>
> My question is: am I misunderstanding something about these coordinate
> transformations to make them not equivalent? A third of an arcsecond is
> significantly big deviation, particularly for HST or interferometry. AFAIK
> fk5 is J2000 in both libraries and galactic coordinates have no concept of
> epoch or equinox time.
>
> Thanks,
> Joseph Booker
> _______________________________________________
> AstroPy mailing list
> AstroPy at scipy.org
> http://mail.scipy.org/mailman/listinfo/astropy
>