[Numpy-discussion] composing Euler rotation matrices

[Seb]

On Tue, 31 Jan 2017 21:23:55 -0500,
Joseph Fox-Rabinovitz <jfoxrabinovitz at gmail.com> wrote:

> Could you show what you are doing to get the statement "However, I
> cannot reproduce this matrix via composition; i.e. by multiplying the
> underlying rotation matrices.". I would guess something involving the
> `*` operator instead of `@`, but guessing probably won't help you
> solve your issue.

Sure, although composition is not something I can take credit for, as
it's a well-described operation for generating linear transformations.
It is the matrix multiplication of two or more transformation matrices.
In the case of Euler transformations, it's matrices specifying rotations
around 3 orthogonal axes by 3 given angles.  I'm using `numpy.dot' to
perform matrix multiplication on 2D arrays representing matrices.

However, it's not obvious from the link I provided what particular
rotation matrices are multiplied and in what order (i.e. what
composition) is used to arrive at the Z1Y2X3 rotation matrix shown.
Perhaps I'm not understanding the conventions used therein.  This is one
of my attempts at reproducing that rotation matrix via composition:

---<--------------------cut here---------------start------------------->---
import numpy as np

angles = np.radians(np.array([30, 20, 10]))

def z1y2x3(alpha, beta, gamma):
    """Z1Y2X3 rotation matrix given Euler angles"""
    return np.array([[np.cos(alpha) * np.cos(beta),
                      np.cos(alpha) * np.sin(beta) * np.sin(gamma) -
                      np.cos(gamma) * np.sin(alpha),
                      np.sin(alpha) * np.sin(gamma) +
                      np.cos(alpha) * np.cos(gamma) * np.sin(beta)],
                     [np.cos(beta) * np.sin(alpha),
                      np.cos(alpha) * np.cos(gamma) +
                      np.sin(alpha) * np.sin(beta) * np.sin(gamma),
                      np.cos(gamma) * np.sin(alpha) * np.sin(beta) -
                      np.cos(alpha) * np.sin(gamma)],
                     [-np.sin(beta), np.cos(beta) * np.sin(gamma),
                      np.cos(beta) * np.cos(gamma)]])

euler_mat = z1y2x3(angles[0], angles[1], angles[2])

## Now via composition

def rotation_matrix(theta, axis, active=False):
    """Generate rotation matrix for a given axis

    Parameters
    ----------

    theta: numeric, optional
        The angle (degrees) by which to perform the rotation.  Default is
        0, which means return the coordinates of the vector in the rotated
        coordinate system, when rotate_vectors=False.
    axis: int, optional
        Axis around which to perform the rotation (x=0; y=1; z=2)
    active: bool, optional
        Whether to return active transformation matrix.

    Returns
    -------
    numpy.ndarray
    3x3 rotation matrix
    """
    theta = np.radians(theta)
    if axis == 0:
        R_theta = np.array([[1, 0, 0],
                            [0, np.cos(theta), -np.sin(theta)],
                            [0, np.sin(theta), np.cos(theta)]])
    elif axis == 1:
        R_theta = np.array([[np.cos(theta), 0, np.sin(theta)],
                            [0, 1, 0],
                            [-np.sin(theta), 0, np.cos(theta)]])
    else:
        R_theta = np.array([[np.cos(theta), -np.sin(theta), 0],
                            [np.sin(theta), np.cos(theta), 0],
                            [0, 0, 1]])
    if active:
        R_theta = np.transpose(R_theta)
    return R_theta

## The rotations are given as active
xmat = rotation_matrix(angles[2], 0, active=True)
ymat = rotation_matrix(angles[1], 1, active=True)
zmat = rotation_matrix(angles[0], 2, active=True)
## The operation seems to imply this composition
euler_comp_mat = np.dot(xmat, np.dot(ymat, zmat))
---<--------------------cut here---------------end--------------------->---

I believe the matrices `euler_mat' and `euler_comp_mat' should be the
same, but they aren't, so it's unclear to me what particular composition
is meant to produce the matrix specified by this Z1Y2X3 transformation.
What am I missing?

-- 
Seb