#### 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" 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 (

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:

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].

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 -

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.

However, one can also carry out the transformation using

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

Quaternions are generally represented in the form:

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:

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:

and then one carries out the quaternion multiplication:

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:

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.

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:

such that:

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

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

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

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

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:

and then carrying out the quaternion multiplication:

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

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:

and also:

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:

to construct the 'rotation quaternion' as:

and its conjugate as:

and the position of the ISS in perifocal coordinates as:

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

And from this, we can immediately read off the rotated 3-vector -

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.

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.

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]

[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]

[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: