Quaternions, rotations and orbital elements

MontBlanc2012

Active member
Joined
Aug 13, 2017
Messages
138
Reaction score
26
Points
28
This is a short (but somewhat technical note) on the rather arcane subject of converting a point in the perifocal reference frame to a more general x-y-z reference frame using quaternions.

To get the ball rolling, let's introduce the perifocal reference frame.


The perifocal reference frame

The "perifocal reference frame" is the natural (and traditional) reference frame in which to calculate the position of a satellite in a standard Keplerian orbit. Specifically, in the perifocal reference frame, the periapsis is located at the 3-vector point [imath]\left(a\,(1-e),\,0,\,0\right)[/imath] where, as usual, [imath]a[/imath] denotes the orbital semi-major axis, and [imath]e[/imath] the orbital eccentricity. Similarly, the point of apoapsis is [imath]\left(-a\,(1+e),\,0,\,0\right)[/imath]. More generally, all points of an elliptical orbit in the perifocal reference frame lie in the x-y plane of that reference frame; such that the orbital periapsis lies on the positive [imath]x[/imath]-axis; and such that the focus of the ellipse (i.e., where the gravitating body is) is located at the origin of the coordinate system - hence the name 'peri'-'focal'.

In the perifocal reference frame, any point on the ellipse can then be written as [imath]\left(r\,\cos(\nu),\, r\,\sin(\nu),\, 0\right)[/imath] where [imath]r[/imath] is the orbital radius; and [imath]\nu[/imath] is the true anomaly. For a Keplerian orbit, the orbital radius is given by the expression:

[math]r = \frac{a\,(1-e^2)}{1 + e\,\cos(\nu)}[/math]​

and where the true anomaly, [imath]\nu[/imath] is usually determined from the mean anomaly, [imath]M[/imath], (a measure of ordinary clock time) via the application of Kepler's equation. In other words, we need just three of the six orbital elements to specify the location of a point on a Keplerian orbit in the perifocal reference frame - the semi-major axis, [imath]a[/imath]; the orbital eccentricity, [imath]e[/imath]; and the mean anomaly, [imath]M[/imath].


Rotation from the perifocal reference frame

All of this is pretty standard stuff. However, one often wants to convert a point in the perifocal reference to a more general reference frame - e.g., Orbiter's Earth-centric equatorial or eccliptic reference frames. This transformation is defined in terms of the three remaining orbital elements: [imath]\Omega[/imath], the longitude of the ascending node; [imath]\iota[/imath], the orbital inclination; and [imath]\omega[/imath], the argument of periapsis. These three orbital elements are, of course, all angles and together form an Euler angle triplet.

From these orbital three elements, the normal procedure for converting from the perifocal reference frame to the equatorial or elliptic frames is to construct a rotation matrix. The procedures for building this rotation matrix from the Euler Angles are well-known and I'm not going to restate them here.


Quaternions

However, one can also carry out the transformation using quaternions. If nothing else, totating 3-vectors using quaternions rather than rotation matrices tends to be more elegant than using rotation matrices. So, in this section, I'll show how to construct (and apply) the quaternions needed to carry out the rotation from the perifocal reference frame.

But first, a few comments about quaternions. Accoridng to Wikipedia,

In mathematics, the quaternions are a number system that extends the complex numbers. They were first described by Irish mathematician William Rowan Hamilton in 1843 and applied to mechanics in three-dimensional space.

Quaternions are generally represented in the form:

[math]\mathcal{Q} = q_0 + q_1\,\mathbf{i} + q_2\,\mathbf{j} + q_3\,\mathbf{k}[/math]​

where [imath]q_0[/imath], [imath]q_1[/imath], [imath]q_2[/imath], and [imath]q_3[/imath] are real numbers, and [imath]\mathbf{i}[/imath], [imath]\mathbf{j}[/imath], and [imath]\mathbf{k}[/imath] are the fundamental unit quaternions. The general rules for adding, subtracting, multiplying and dividing quaternions is treated well in innumerable books and articles, so I'm going to assume that you either know them or can Google them. It is worth pointing out that the conjugate [imath]\mathcal{Q}^\dagger[/imath] of the quaternion [imath]\mathcal{Q}[/imath] is given by:

[math]\mathcal{Q}^\dagger = q_0 - q_1\,\mathbf{i} - q_2\,\mathbf{j} - q_3\,\mathbf{k}[/math]​

Given all of this, a rotation of a 3-vector [imath](x,\,y,\,z)[/imath], say, is carried out using quaternions by first constructing a new quaternion, [imath]\mathcal{P}[/imath], such that:

[math]\mathcal{P} = 0 + x\,\mathbf{i} + y\,\mathbf{j} + z\,\mathbf{k}[/math]​

and then one carries out the quaternion multiplication:

[math]\mathcal{P}' = \mathcal{Q}\,\mathcal{P}\,\mathcal{Q}^\dagger[/math]​

such that [imath]q_0^2 + q_1^2 + q_2^2 + q_3^2 = 1[/imath]. In this case, the new quaternion, [imath]\mathcal{P}'[/imath] is of the form:

[math]\mathcal{P}' = 0 + X\,\mathbf{i} + Y\,\mathbf{j} + Z\,\mathbf{k}[/math]​

so that the initial 3-vector [imath](x,\,y,\,z)[/imath] after rotation can simply be read off as [imath](X,\,Y,\,Z)[/imath]. All of this is actually rather elegant - and all one needs to know is how to multiply quaternions.


Application to orbital coordinate transformations

Of course we need to know how to construct the quaternion [imath]Q[/imath] so that we can carry out the rotation. As it turns out, this is pretty simple. In terms of the Euler angles, [imath]\Omega[/imath] (the longitude of the ascending node), [imath]\iota[/imath] (the orbital inclination); and [imath]\omega[/imath] (the argument of periapsis), we construct the quaternions components as:

[math]\begin{aligned} q_0 &= \cos\left(\frac{\iota}{2}\right)\,\cos\left(\frac{\Omega + \omega}{2}\right) \\ q_1 &= \sin\left(\frac{\iota}{2}\right)\,\cos\left(\frac{\Omega - \omega}{2}\right) \end{aligned}[/math]
[math]\begin{aligned} q_2 &= \sin\left(\frac{\iota}{2}\right)\,\sin\left(\frac{\Omega - \omega}{2}\right) \\ q_3 &= \cos\left(\frac{\iota}{2}\right)\,\sin\left(\frac{\Omega + \omega}{2}\right) \\ \end{aligned}[/math]​

such that:

[math]\mathcal{Q} = q_0 + q_1\,\mathbf{i} + q_2\,\mathbf{j} + q_3\,\mathbf{k}[/math]​

We then construct the conjugate of [imath]\mathcal{Q}[/imath]:

[math]\mathcal{Q}^\dagger = q_0 - q_1\,\mathbf{i} - q_2\,\mathbf{j} - q_3\,\mathbf{k}[/math]​

Then for a 3-vector (in the perifocal reference frame) that we want to rotate, [imath](x,\,y,\,z)[/imath], we construct:

[math]\mathcal{P} = 0 + x\,\mathbf{i} + y\,\mathbf{j} + z\,\mathbf{k}[/math]​

And, finally, we calculate the rotated form of the 3-vector by calculating:

[math]\mathcal{P}' = 0 + X\,\mathbf{i} + Y\,\mathbf{j} + Z\,\mathbf{k} = \mathcal{Q}\,\mathcal{P}\,\mathcal{Q}^\dagger[/math]​

And so we read off the rotated form of the 3-vector as [imath](X,\,Y,\,Z)[/imath].


Inverting the transformation

If we want to transform from an equatorial or ecliptic reference frame back to the perifocal reference frame, this can easily be achieved by constructing [imath]P'[/imath] as:

[math]\mathcal{P}' = 0 + X\,\mathbf{i} + Y\,\mathbf{j} + Z\,\mathbf{k}[/math]​

and then carrying out the quaternion multiplication:

[math]\mathcal{P} = \mathcal{Q}^\dagger\,\mathcal{P}'\,\mathcal{Q}[/math]​

And then reading off the [imath](x,\,y,\,z)[/imath] coordinates from the result. Simple!


A worked example

OK, let's work through a more or less realistic example of how this all works: Let's take a specific example from the "Docked at the ISS" XR2 Ravenstar stock scenario. With the scenario paused at MJD = 51984.64994, and using the Scenarion Editor, we can read off directly the orbital elements of the ISS' orbit:



Clearly, in the Earth-centric ecliptic reference frame we have:

[math]\begin{aligned} a &= 6735949.639\\ e &= 0.00100408 \\ \nu &= 256.384^\circ \end{aligned}[/math]​

and also:

[math]\begin{aligned} \Omega &= 162.194^\circ \\ \iota &= 73.681^\circ \\ \omega &= 112.48^\circ \end{aligned}[/math]​

where we have used the identity that the longitude of periapsis ([imath]274.674^\circ[/imath]) is just the sum of the argument of periapsis ([imath]\omega[/imath]) and the longitude of the ascending node ([imath]162.194^\circ[/imath]).

From this, we can use:

[math]\begin{aligned} q_0 &= \cos\left(\frac{\iota}{2}\right)\,\cos\left(\frac{\Omega + \omega}{2}\right) = -0.5885082\\ q_1 &= \sin\left(\frac{\iota}{2}\right)\,\cos\left(\frac{\Omega - \omega}{2}\right) = 0.5440433 \end{aligned}[/math]
[math]\begin{aligned} q_2 &= \sin\left(\frac{\iota}{2}\right)\,\sin\left(\frac{\Omega - \omega}{2}\right) = 0.2520404\\ q_3 &= \cos\left(\frac{\iota}{2}\right)\,\sin\left(\frac{\Omega + \omega}{2}\right) = 0.5423565\\ \end{aligned}[/math]​

to construct the 'rotation quaternion' as:

[math]\mathcal{Q} = -0.5885082 + 0.5440433\,\mathbf{i} + 0.2520404\,\mathbf{j} + 0.5423565\,\mathbf{k}[/math]​

and its conjugate as:

[math]\mathcal{Q}^\dagger = -0.5885082 - 0.5440433\,\mathbf{i} - 0.2520404\,\mathbf{j} - 0.5423565\,\mathbf{k}[/math]​

and the position of the ISS in perifocal coordinates as:

[math]\mathcal{P} = \frac{a\,(1-e^2)}{1 + e\,\cos(\nu)}\,\left(0 + \cos(\nu)\,\mathbf{i} + \sin(\nu)\,\mathbf{j} + 0\,\mathbf{k}\right) = 0 -1586106.976\,\mathbf{i} -6548179.005\,\mathbf{j} + 0\,\mathbf{k}[/math]​

Finally, we can carry out the rotation to the ecliptic frame of this perifocal position vector by applying the quaternion rotation:

[math]\begin{aligned} \mathcal{P}' &= \mathcal{Q}\,\mathcal{P}\,\mathcal{Q}^\dagger \\ &= 0 - 6427381.91\,\mathbf{i} + 1757957.97\,\mathbf{j} + 996357.96\,\mathbf{k} \end{aligned}[/math]​

And from this, we can immediately read off the rotated 3-vector - i.e., the x-y-z coordinates of the ISS' position vector in the Earth-centric ecliptic reference frame as [imath](-6427381.91,\, 1757957.97,\, 996357.96)[/imath]. Finally, to convert from a conventional right-handed coordinate system to Orbiter's somewhat quirky left-handed coordinate system, we swap the 'y' and 'z' components to end up with an Orbiter compatible representation [imath](-6427381.91,\, 996357.96,\, 1757957.97)[/imath].

How does this result compare with the state vector components calculated by Orbiter? Again using the Scenario Orbiter, we can read off the state vector components directly:



as [imath](-6427366.6,\, 996435.8,\, 175975.5)[/imath]. This isn't exactly the same as the above but it is close - with 100 meters or so. Why the difference? Well, the Scenario Editor only reports the orbital elements to a relatively low precision (about five or six significant figures). If we were to use the full machine precision, obtained by interrogating Orbiter's system via the API, the two methods would yield the same results with very high precision.


In summary

Quaternions do the same job as rotation matrices in mapping points in one coordinate system to another. Arguably, however, they are more elegant - even if they may not be hugely more computationally efficient. In any event, one comes across them periodically - so understanding what they are and how they can be used is no bad thing.
 
Last edited:

GLS

Well-known member
Orbiter Contributor
Addon Developer
Joined
Mar 22, 2008
Messages
5,877
Reaction score
2,870
Points
188
Website
github.com
Arguably, however, they are more elegant - even if they may not be hugely more computationally efficient.

Not an expert, but just now during lunch I was reading about quaternions in the shuttle, and they started using them because it was more efficient. :shrug:
 

MontBlanc2012

Active member
Joined
Aug 13, 2017
Messages
138
Reaction score
26
Points
28
Not an expert, but just now during lunch I was reading about quaternions in the shuttle, and they started using them because it was more efficient.

That may indeed be the case. However, I just think quaternions are interesting.
 

martins

Orbiter Founder
Orbiter Founder
Joined
Mar 31, 2008
Messages
2,448
Reaction score
462
Points
83
Website
orbit.medphys.ucl.ac.uk
One of the nice things about quaternions is that they describe the rotation more concisely. A rotation matrix consists of 9 scalar values, but there is some redundancy since the conditions of orthogonality and normalisation introduce constraints that don't allow for just any combination of 9 values. Indeed, rotation matrices that are continuously updated tend to need periodic re-orthogonalisation and re-normalisation, since numerical errors build up over time. Quaternions don't have that problem.
 

KnightOwl

New member
Joined
May 27, 2022
Messages
1
Reaction score
1
Points
0
Location
Earth
Sorry to give this a revive, but I thought it might be interesting to know that this very implementation has a huge performance improvement when implementing orbital mechanics in a simulation game. We happened to implement this very concept in our game, LithoBreak, not knowing if it had ever been done before. We went from being able to simulate dozens of orbits to thousands at 60 FPS.

By the way, does anyone know of further references that might help calculate transfer orbits?
 

MontBlanc2012

Active member
Joined
Aug 13, 2017
Messages
138
Reaction score
26
Points
28
Sorry to give this a revive, but I thought it might be interesting to know that this very implementation has a huge performance improvement when implementing orbital mechanics in a simulation game. We happened to implement this very concept in our game, LithoBreak, not knowing if it had ever been done before. We went from being able to simulate dozens of orbits to thousands at 60 FPS.

By the way, does anyone know of further references that might help calculate transfer orbits?
What sort of transfer orbits do you want to calculate?
 
Top