5
$\begingroup$

I am trying to implement the conversion between TEME and GCRF frames following the same procedure as done in SPICE. After digging into CSPICE's code, this transformation is performed through the zzteme routine. Following the function calls in that code leads me to finding out that the transformation is as follows:

  • Obtain a rotation to bring coordinates in TEME frame to MEME (I guess Mean Equator Mean Equinox) frame. This is done by using 1980 IAU theory of nutation, which is achieved by CSPICE's routine zzenut80. This routine basically calculates 3 angles describing nutation (mean obliquity, nutation in obliquity and nutation in longitude), together with their time-derivatives, and then combines these 6 values to obtain a series of 3 Euler angles describing an XZX rotation together with the time derivative of each of the 3 Euler angles. Note that the resulting transformation matrix is not a 3-by-3 matrix, but instead a 6-by-6 transformation matrix which multiplies a vector of length 6 including position and velocity. But this is equivalent I believe to applying the transportation theorem to the velocity as explained in this answer.
  • Obtain a rotation to bring coordinates in MEME frame to J2000 (which I am treating for now as equivalent to GCRF). This is done with 1976 IAU theory of precession, implemented in zzeprec76. This routine outputs as well 3 angles describing precession (zeta, theta and z) together with their time-derivatives, which are used to construct another 6-by-6 state transformation matrix as is done for nutation (in this case though, rotation sequence is ZYZ).
  • The two transformations are then combined, resulting in a transformation matrix from TEME to J2000.

So in summary, nutation and precession angles are calculated together with their derivatives, these are then combined into 2 sets of rotations described each by 3 Euler angles and their derivatives, then the rotations described by Euler angles and their derivatives are converted into state transformation matrices (which differ from rotation matrices in that they also account for conversion of velocities) and then the transformations are applied.

I understand how to obtain nutation and precession angles, and it is very straightforward to obtain the corresponding Euler angles and sequence of rotations for both nutation and precession. I also have found multiple sources indicating how to convert a sequence of rotations described by 3 Euler angles into a quaternion or DCM for each of the 12 valid rotation sequences, such as this.

However, I am wondering how would one convert the Euler angle derivatives into the angular velocity vector as required to apply the transport theorem? I guess this will also need to take into account the rotation sequence?

$\endgroup$
4
  • $\begingroup$ The angular derivatives are the angular velocity. $\endgroup$ Commented Jan 16 at 19:37
  • 2
    $\begingroup$ @GregMiller Angular velocity in three dimensional space is a vector quantity in the sense of how they transform from one frame to another. If one frame is rotating with respect to another, the transport theorem also needs to be taken into account. The derivatives of an Euler angle sequence do not behave as do vectors. $\endgroup$ Commented Jan 16 at 22:43
  • 1
    $\begingroup$ they're generally called "Euler rates" or "Euler-angle rates" in my experience. Tons of literature out there about how to convert between those and angular velocity vectors $\endgroup$
    – Erin Anne
    Commented Jan 17 at 1:12
  • $\begingroup$ Thanks! Yes, I have found multiple sources describing the conversion for specific rotation sequences, such as this , providing conversions for rotation sequences ZYX and ZXZ in equations 22 and 24. But I was looking for possibly some reference covering all 12 valid rotation sequences $\endgroup$
    – Rafa
    Commented Jan 17 at 1:21

1 Answer 1

5
$\begingroup$

However, I am wondering how would one convert the Euler angle derivatives into the angular velocity vector as required to apply the transport theorem? I guess this will also need to take into account the rotation sequence?

Regarding the second question, yes you will need to take into account the rotation sequence as matrix multiplication in general is not commutative (i.e., $AB \ne BA$ in general, but see When is matrix multiplication commutative?). Regarding the first question, consider four frames of reference $A$, $B$, $C$, and $D$, each of which has three orthogonal basis vectors, ordered such that they follow the right hand rule.

Preliminaries

The transformation matrices $T_{A\to B}$, $T_{B\to C}$, and $T_{C\to D}$ respectively transform a vector from frame $A$ to frame $B$, from $B$ to $C$, and from $C$ to to $D$. As the reference frames all follow the right hand rule, each of these transformation matrices will inherently represent a proper rotation1 (one that does not involve reflections), i.e., a member of $\operatorname{SO}(3)$. A question arises: What is the transformation matrix from the original frame $A$ to the final frame $D$? The answer is

$$T_{A\to D} = T_{C\to D} T_{B\to C} T_{A\to B} \tag{1}$$

Suppose that some (and perhaps all) of these transformations vary with time. Taking the time derivative of (1) results in

$$\dot T_{A\to D} = \dot T_{C\to D} T_{B\to C} T_{A\to B} + T_{C\to D} \dot T_{B\to C} T_{A\to B} + T_{C\to D} T_{B\to C} \dot T_{A\to B}\tag{2}$$

It will be advantageous to look at the transpose of equation (1):

$${T_{A\to D}}^T = T_{D\to A} = \left(T_{C\to D} T_{B\to C} T_{A\to B}\right)^T = T_{B\to A} T_{C\to B} T_{D\to C} \tag{3}$$

Differentiating with respect to time yields

$$\dot T_{D\to A} = \dot T_{B\to A} T_{C\to B} T_{D\to C} + T_{B\to A} \dot T_{C\to B} T_{D\to C} + T_{B\to A} T_{C\to B} \dot T_{D\to C}\tag{4}$$

On the transform of a cross product -- and of a $3\times 3$ skew symmetric matrix

In three dimensional space, a skew symmetric matrix can be represented by a vector $\vec s$ (more correctly, by a pseudovector $\vec s$) that is typically denoted as $\left[\vec s\right]_x$. This is the cross product matrix generated from the vector $\vec s$: $\left[\vec s\right]_x \vec a \equiv \vec s \times \vec a$. If $\vec s = [s_1, s_2, s_3]^T$ then the cross product matrix generated from this vector is

$$\left[\vec s\right]_x = \left[\begin{matrix} \phantom{-}0& -s_3 & \phantom{-}s_2 \\ \phantom{-}s_3& \phantom{-}0 & -s_1 \\ -s_2 & \phantom{-}s_1 & \phantom{-}0 \end{matrix}\right]\tag{5}$$

Generating a pseudovector given a 3x3 skew symmetric matrix is similarly trivial: Simply read off the components of the pseudovector using equation (5).

Given pair of reference frames $X$ and $Y$ in $\mathbb R_3$, a proper (det=1) orthonormal transformation matrix between these frames $T_{X\to Y}$, and a pair of vectors $\vec a_X$ and $\vec b_X$ represented in frame $X$ coordinates, the cross product of the corresponding vectors in frame $Y$ is given by $\left(T_{X\to Y} \vec a_X\right)\times\left(T_{X\to Y} \vec b_X\right)$. This is equal to $T_{X\to Y}(\vec a_X \times \vec b_X)$ as the cross product transforms as does an ordinary vector (given a proper orthogonal transformation matrix). Since this is the case for all $\vec b$, the cross product matrix formed from $\vec a$ transforms as $$T_{X\to Y} [\vec a_X]_x T_{X\to Y}^T = [T_{X\to Y} \vec a_X]_x\tag{6}$$

On the time derivative of a $3x3$ transformation matrix

Equation (4) is a bit of a mess. To clean this up, it will be helpful to take a side track on the nature of the time derivative of a proper orthonormal transformation matrix $T$ (note thee lack of subscripts). That $TT^T=I$ means that time derivative of a real transformation matrix $T$ is antisymmetric as $\frac{dT}{dt} T^T + T\frac{dT^T}{dt} = \dot I = 0$ and $\frac{dT^T}{dt} = \left(\frac{dT}{dt}\right)^T$. This in turn means that there exists skew symmetric matrices $S_L$ and $S_R$ such that $$\dot T = S_L T = T S_R \tag{7}$$

The subscripts $_L$ and $_R$ denote whether the skew symmetric matrix is to the left or right of the transformation matrix $T$.

The transport theorem yields a clue regarding the physical nature of $S_R$: $$\frac{d \vec{q}_Y}{dt} = T_{X\to Y} \left(\frac{d \vec{q}_X}{dt} + \vec{\omega}_{Y\to X:X}\times \vec{q}_X\right)$$ Here, $\vec q$ is some vector quantity in three dimensional Cartesian space ($\mathbb R_3$), $X$ and $Y$ are two reference frames for $\mathbb R_3$ with a common origin with right-handed Cartesian bases, and $\vec{\omega}_{Y\to X:X}$ is the angular velocity of frame $X$ with respect to frame $Y$, expressed in frame $X$ coordinates. The transport theorem results from taking the time derivative of $\vec q_Y = T_{X\to Y} \vec q_X$. This in turn implies that $$\frac{d T_{X\to Y}}{dt} = T_{X\to Y} \left[ \vec{\omega}_{Y\to X:X}\right]_x\tag{8}$$ This means that the skew symmetric matrix $S_R$ in equation (7) is the skew symmetric matrix formed from the angular velocity of frame $X$ with respect to frame $Y$, expressed in frame $X$ coordinates. Given that $S_L = T S_R T^T$ from equation (7), the vector that underlies $S_L$ in equation (7) is $$T_{X\to Y} \vec{\omega}_{Y\to X:X} = \vec{\omega}_{Y\to X:Y}$$ In other words, $S_L$ in equation (7) is the skew symmetric matrix formed from the angular velocity of frame $X$ with respect to frame $Y$, expressed in frame $Y$ coordinates.

Angular velocity due to a sequence of time-varying transformations

With the above, equation (4) can be rewritten as $$\begin{aligned} \dot T_{D\to A} \equiv T_{D\to A}\left[\vec{\omega}_{A\to D:D} \right]_x =& \phantom{+\,\,\;}T_{B\to A} \left[\vec{\omega}_{A\to B:B}\right]_x T_{C\to B} T_{D\to C} \\ & + T_{B\to A} T_{C\to B} \left[\vec{\omega}_{B\to C:C}\right]_x T_{D\to C} \\ &+ T_{B\to A} T_{C\to B} T_{D\to C} \left[\vec{\omega}_{D\to C:D}\right]_x \end{aligned}\tag{9}$$ This is still a mess, but the solution is straightforward: Move the intermediate skew symmetric matrices to the right2. This yields $$\begin{aligned} T_{D\to A}\left[ \vec{\omega}_{A\to D:D} \right]_x =& \phantom{+\,\,\;}T_{B\to A} T_{C\to B} T_{D\to C} \left[T_{C\to D} T_{B\to C}\vec{\omega}_{A\to B:B}\right]_x \\ & + T_{B\to A} T_{C\to B} T_{D\to C} \left[T_{C\to D} \vec{\omega}_{B\to C:C}\right]_x \\ &+ T_{B\to A} T_{C\to B} T_{D\to C} \left[\vec{\omega}_{D\to C:D}\right]_x \end{aligned}\tag{10}$$ Or, more simply, $$\vec{\omega}_{A\to D:D} = T_{C\to D} T_{B\to C}\vec{\omega}_{A\to B:B} + T_{C\to D} \vec{\omega}_{B\to C:C} + \vec{\omega}_{D\to C:D}\tag{11}$$

Footnotes

1 A proper 3x3 orthogonal transformation matrix has several key essential properties:

  • For any two columns (or any two rows) $i$ and $j$ of the matrix, the inner product of the vectors extracted from those columns (or rows) is $\delta_{ij}$: +1 if $i=j$, 0 otherwise.
  • The cross product of any two columns $i$ and $j$ of the matrix is either the omitted column $k$ or the negation of that column, depending on whether $ijk$ is an even or odd permutation of 1,2,3.
  • The determinant of the transformation matrix is +1. An improper transformation matrix has a determinant of -1. Such improper transformations can result from transforming between frames that involve a reflection such as a between right handed and left handed coordinate systems.
  • One of the eigenvalues of a proper 3x3 orthonormal matrix will be +1. The product of the other two eigenvalues will also be +1. These other two eigenvalues can be real (in which case they are both +1 or both -1) or can be a complex conjugate pair. Only in the case of the identity matrix are the other two eigenvalues also +1.
  • Except for the identity matrix, the eigenvector corresponding to the +1 eigenvalue represents the axis of rotation. (Any axis can serve as an axis of rotation for the identity matrix as a rotation of zero degrees about any axis results in the identity matrix.)
  • All proper 3x3 orthonormal matrices can be expressed as the matrix exponential of a 3x3 skew symmetric matrix. The matrix exponential of a 3x3 skew symmetric matrix is a proper 3x3 orthonormal matrix.
  • The 3x3 skew symmetric matrices form an algebra, $\mathfrak{so}(3)$. The 3x3 proper orthonormal transformation matrices form a group, $\operatorname{SO}(3)$. The relation via the exponential map between $\mathfrak{so}(3)$ and $\operatorname{SO}(3)$ means these are respectively a Lie algebra and a Lie group.

2 Angular velocity is typically expressed in some body fixed frame because that is the frame in which an observer fixed with respect to the body observes angular velocity. For example, a rate gyro rigidly attached to a spacecraft body measures the angular velocity with respect to inertial, expressed in a body-fixed case frame.

References

Gallo, Eduardo. "The SO (3) and SE (3) Lie Algebras of Rigid Body Rotations and Motions and their Application to Discrete Integration, Gradient Descent Optimization, and State Estimation." arXiv preprint arXiv:2205.12572] (2023).
Chapter 4 is relevant to the discussion above.

Zhao, Shiyu. "Time derivative of rotation matrices: A tutorial." arXiv preprint arXiv:1609.06088 (2016).
A short but non-handwaving derivation of equation (8).

International Astronomical Union. "SOFA Tools for Earth Attitude." IAU Standards Of Fundamental Astronomy (2017).
Compared to JPL's SPICE, the IAU SOFA software is easier to use, easier to read, better documented, better maintained, and much more up-to-date.

Engø, Kenth. "On the BCH-formula in so (3)." BIT Numerical Mathematics 41 (2001): 629-632.
The Baker–Campbell–Hausdorff formula relates the exponential maps of two values $X$ and $Y$ to the exponential of a third value: $\exp(X)\exp(Y) = \exp(Z)$. Solving this for $Z$ is in general non-trivial. The referenced article provides a closed form solution for $X,Y,Z \in \mathfrak{so}(3)$. The link is to the official, paywalled version of the paper. Finding a non-paywalled version is a task left for the reader of this answer.

Iserles, Arieh, et al. "Lie-group methods." Acta numerica 9 (2000): 215-365.
Lie group methods are (IMHO) the correct way to numerically integrate transformation matrices. This huge paper takes a very deep dive down the abstract algebra rabbit hole. As with the above reference, finding a non-paywalled version is a task left for the reader of this answer.

$\endgroup$
2
  • 1
    $\begingroup$ Thanks a lot for taking the time to write this excellent (and well-referenced) answer. As with time systems, the more I read about transformation of coordinate systems, the bigger the rabbit hole seems to become! But I now see a clear path to doing the transformation through equation 11 together with transport theorem. Luckily I have not (yet) gotten into transforming accelerations too... $\endgroup$
    – Rafa
    Commented Jan 29 at 2:38
  • 2
    $\begingroup$ @Rafa It is possible, at least for quaternions as opposed to matrices; see for example the documentation regarding how the Johnson Engineering Orbital Dynamics package (JEOD) on quaternions. In particular, look at the appendices. $\endgroup$ Commented Jan 29 at 10:03

Not the answer you're looking for? Browse other questions tagged or ask your own question.