Merge pull request astropy#2619 from eteq/fix-rot-matrix

Fix quirks in rotation_matrix

Author: eteq

CHANGES.rst

@@ -41,6 +41,10 @@ New Features
   - The deprecated functions for pre-0.3 coordinate object names like
     ``ICRSCoordinates`` have been removed. [#2422]
 
+  - The ``rotation_matrix`` and ``angle_axis`` functions in
+    ``astropy.coordinates.angles`` were made more numerically consistent and
+    are now tested explicitly [#2619]
+
 - ``astropy.cosmology``
 
   - Added ``z_at_value`` function to find the redshift at which a cosmology

astropy/coordinates/angles.py

@@ -701,54 +701,41 @@ def rotation_matrix(angle, axis='z', unit=None):
     rmat: `numpy.matrix`
         A unitary rotation matrix.
     """
-    # TODO: This doesn't handle arrays of angles
-
     if unit is None:
         unit = u.degree
 
     angle = Angle(angle, unit=unit)
 
+    s = np.sin(angle)
+    c = np.cos(angle)
+
+    # use optimized implementations for x/y/z
     if axis == 'z':
-        s = np.sin(angle)
-        c = np.cos(angle)
         return np.matrix(((c, s, 0),
                           (-s, c, 0),
                           (0, 0, 1)))
     elif axis == 'y':
-        s = np.sin(angle)
-        c = np.cos(angle)
         return np.matrix(((c, 0, -s),
                           (0, 1, 0),
                           (s, 0, c)))
     elif axis == 'x':
-        s = np.sin(angle)
-        c = np.cos(angle)
         return np.matrix(((1, 0, 0),
                           (0, c, s),
                           (0, -s, c)))
     else:
-        x, y, z = axis
-        w = np.cos(angle / 2)
-
-        # normalize
-        if w == 1:
-            x = y = z = 0
-        else:
-            l = np.sqrt((x * x + y * y + z * z) / (1 - w * w))
-            x /= l
-            y /= l
-            z /= l
+        axis = np.asarray(axis)
+        axis = axis / np.sqrt((axis * axis).sum())
 
-        wsq = w * w
-        xsq = x * x
-        ysq = y * y
-        zsq = z * z
-        return np.matrix(((wsq + xsq - ysq - zsq, 2 * x * y - 2 * w * z, 2 * x * z + 2 * w * y),
-                          (2 * x * y + 2 * w * z, wsq - xsq + ysq - zsq, 2 * y * z - 2 * w * x),
-                          (2 * x * z - 2 * w * y, 2 * y * z + 2 * w * x, wsq - xsq - ysq + zsq)))
+        R = np.diag((c, c, c))
+        R += np.outer(axis, axis) * (1. - c)
+        axis *= s
+        R += np.array([[0., axis[2], -axis[1]],
+                       [-axis[2], 0., axis[0]],
+                       [axis[1], -axis[0], 0.]])
+        return R.view(np.matrix)
 
 
-def angle_axis(matrix, unit=None):
+def angle_axis(matrix):
     """
     Computes the angle of rotation and the rotation axis for a given rotation
     matrix.
@@ -758,29 +745,20 @@ def angle_axis(matrix, unit=None):
     matrix : array-like
         A 3 x 3 unitary rotation matrix.
 
-    unit : UnitBase
-        The output unit.  If `None`, the output unit is degrees.
-
     Returns
     -------
     angle : `Angle`
         The angle of rotation for this matrix.
 
     axis : array (length 3)
-        The axis of rotation for this matrix.
+        The (normalized) axis of rotation for this matrix.
     """
-    # TODO: This doesn't handle arrays of angles
-
     m = np.asmatrix(matrix)
     if m.shape != (3, 3):
         raise ValueError('matrix is not 3x3')
 
-    angle = np.acos((m[0, 0] + m[1, 1] + m[2, 2] - 1) / 2)
-    denom = np.sqrt(2 * ((m[2, 1] - m[1, 2]) + (m[0, 2] - m[2, 0]) + (m[1, 0] - m[0, 1])))
-    axis = np.array((m[2, 1] - m[1, 2], m[0, 2] - m[2, 0], m[1, 0] - m[0, 1])) / denom
-    axis /= np.sqrt(np.sum(axis ** 2))
+    axis = np.array((m[2, 1] - m[1, 2], m[0, 2] - m[2, 0], m[1, 0] - m[0, 1]))
+    r = np.sqrt((axis * axis).sum())
+    angle = np.arctan2(r, np.trace(m) - 1)
 
-    angle = Angle(angle, u.radian)
-    if unit is None:
-        unit = u.degree
-    return angle.to(unit), axis
+    return Angle(angle, u.radian), -axis / r

astropy/coordinates/tests/test_angles.py

@@ -753,3 +753,47 @@ def test_mixed_string_and_quantity():
     a2 = Angle(['1d', 1 * u.rad * np.pi, '3d'])
     assert_array_equal(a2.value, [1., 180., 3.])
     assert a2.unit == u.deg
+
+def test_rotation_matrix():
+    from ..angles import rotation_matrix
+
+    assert_array_equal(rotation_matrix(0*u.deg, 'x'), np.eye(3))
+
+    assert_allclose(rotation_matrix(90*u.deg, 'y'), [[ 0, 0,-1],
+                                                     [ 0, 1, 0],
+                                                     [ 1, 0, 0]], atol=1e-12)
+
+    assert_allclose(rotation_matrix(-90*u.deg, 'z'), [[ 0,-1, 0],
+                                                      [ 1, 0, 0],
+                                                      [ 0, 0, 1]], atol=1e-12)
+
+    assert_allclose(rotation_matrix(45*u.deg, 'x'),
+                    rotation_matrix(45*u.deg, [1, 0, 0]))
+    assert_allclose(rotation_matrix(125*u.deg, 'y'),
+                    rotation_matrix(125*u.deg, [0, 1, 0]))
+    assert_allclose(rotation_matrix(-30*u.deg, 'z'),
+                    rotation_matrix(-30*u.deg, [0, 0, 1]))
+
+    assert_allclose(np.dot(rotation_matrix(180*u.deg, [1, 1, 0]).A, [1, 0, 0]),
+                    [0, 1, 0], atol=1e-12)
+
+    #make sure it also works for very small angles
+    assert_allclose(rotation_matrix(0.000001*u.deg, 'x'),
+                    rotation_matrix(0.000001*u.deg, [1, 0, 0]))
+
+def test_angle_axis():
+    from ..angles import rotation_matrix, angle_axis
+
+    m1 = rotation_matrix(35*u.deg, 'x')
+    an1, ax1 = angle_axis(m1)
+
+    assert an1 - 35*u.deg < 1e-10*u.deg
+    assert_allclose(ax1, [1, 0, 0])
+
+
+    m2 = rotation_matrix(-89*u.deg, [1, 1, 0])
+    an2, ax2 = angle_axis(m2)
+
+    assert an2 - 89*u.deg < 1e-10*u.deg
+    assert_allclose(ax2, [-2**-0.5, -2**-0.5, 0])
+