[AstroPy] Numerical issues with rotation_matrix

[Maik Riechert]

Hi Thomas,

>>>>> 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]])
> 
> These differences are close to what might be expected for machine
> precision. However, I don't see why we couldn't try and recognize
> special cases (for example [1,0,0] -> x). In your case, why not simply
> use 'x' since it gives the expected result?

I did some more testing and used sympy (symbolic execution/arbitrary
precision floats) to give me the "ground truth". For example, using an
angle of 0.000004deg the difference between sympy and gohlke/astropy's
special cases is at most e-18. When an array is used with astropy, the
difference is e-7, which I think is significant as it may add up when
multiple rotations are combined.

Tell me if you like to see the test code, then I prettify it a bit.

Cheers
Maik