[AstroPy] Numerical issues with rotation_matrix

[Maik Riechert]

Hi,

I discovered today that astropy has a rotation_matrix function (in
coordinates.angles) and wanted to use it instead of the one from
Christoph Gohlke:
http://www.lfd.uci.edu/~gohlke/code/transformations.py.html

While doing some testing, I noticed that the values differ between

rotation_matrix(angle, [1,0,0])
and
rotation_matrix(angle, 'x')

The latter version agrees with Gohlke's function. The difference between
both versions in astropy is around the order of e-13 in my case, which
leads to severe differences in further processing.

Looking at the source code it seems as if the general array version uses
sqrt and other things which may introduce numerical issues.

I think a test case should be added which checks that 'x'/'y'/'z'
produces the same values as [1,0,0] etc.

My test case:

>>> import numpy as np
>>> from astropy.coordinates.angles import rotation_matrix as rot_astro
>>> from gohlke.transformations import rotation_matrix as rot_gohlke
>>> m1 = rot_astro(-0.07687069267192707, [1,0,0])
>>> m2 = rot_astro(-0.07687069267192707, 'x')
>>> m2-m1
matrix([[  0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   2.22044605e-16,   1.30859299e-13],
        [  0.00000000e+00,  -1.30859299e-13,   2.22044605e-16]])

>>> m3 = rot_gohlke(np.deg2rad(-0.07687069267192707), [-1,0,0])[:3,:3]
>>> m3-m2
matrix([[ 0.,  0.,  0.],
        [ 0.,  0.,  0.],
        [ 0.,  0.,  0.]])

I'm not completely sure which one is the more accurate one. This has to
be validated. But if I should take a guess, it's m2/m3.

Cheers
Maik