[AstroPy] Numerical issues with rotation_matrix

[Erik Tollerud]

Hello Maik,

Just to clarify something: `rotation_matrix` is technically not part
of the public API for astropy.  It's mainly used internally in
coordinates. Furthermore, it's not actually used internally for
anything other than the 'x', 'y', and 'z' use cases. So just be aware
that this function comes with no "guarantees" like we do for the
public API.

If you actually look at the code, you'll see the matricies for the
'x', 'y', and 'z' cases are constructed quite differently from the
generic angle, such that we do *not* expect them to come out exactly
the same to machine precision.  So I think special-casing [1,0,0],
[0,1,0], and [0,0,1] to be the same as 'x', 'y', and 'z' is bad,
because it would cause discontinuous jumps when you're using the
axis-specified form.  So I would say the fact that they are not the
same to machine precision is a feature, not a bug (because the
'x'/'y'/'z' cases are quite a bit faster to run).

I am concerned about the "When an array is used with astropy, the
difference is e-7" bit, though: can you clarify a bit more what you
mean by "when an array is used with astropy"?




On Thu, Jun 5, 2014 at 3:13 PM, Maik Riechert <maik.riechert at arcor.de> wrote:
> 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
> _______________________________________________
> AstroPy mailing list
> AstroPy at scipy.org
> http://mail.scipy.org/mailman/listinfo/astropy



-- 
Erik T