We demonstrate how the equation for the location of a satellite orbiting a large primary can be solved with the help of quaternions. We also treat the issue of orbit determination, given (1) a location and a velocity of the satellite; (2) two observations of its location; and (3) three or more values of its position.

Quaternions

We define a quaternion as a symbolic sum of a scalar, say , and a three-dimensional vector , collectively denoted . We add and subtract quaternions component by component; multiplication (denoted by ) of two quaternions is done by

(1)

where and are the dot and cross products, respectively. One can easily verify that this multiplication is associative (multiplication of three or more quaternions is meaningful without parentheses), but it is clearly noncommutative (changing the order of quaternions in a product will, in general, yield a different answer). Taking a quaternion’s conjugate changes the sign of its vector part, thus: . Based on the corresponding Taylor expansion, we can easily construct functions of quaternion arguments, such as, for example,

(2)

where is the magnitude of and is the corresponding unit vector. (In general, the answer for is the same as it would be for a complex argument, with assuming the role of .)

3D Rotations

One can easily verify that rotating a vector with respect to (whose direction and magnitude represent the axis and angle of rotation, respectively) is achieved by

(3)

where both and are shorthand for the corresponding pure-vector quaternions (i.e., and , respectively). This enables us to easily find any composition of two or more rotations by simply multiplying the corresponding quaternions.

Here are routines for multiplying two quaternions, taking a conjugate, and evaluating (2), respectively.

Next, , , and are unit vectors in the direction of , , and , respectively.

The following multiplication finds the composition of a 90° rotation around followed by a 90° rotation around .

The result is a rotation around the direction.

A general 3D rotation can be also parametrized by three Euler angles , , and , thus:

(4)

and visualized as a composition of three rotations (around , , and again, by the angles , , and respectively, and in that order).

Satellite in Orbit

It is well known [1] that the equation for establishing a satellite’s location (relative to the primary’s center), as a function of time , reads

(5)

where is a the sum of the two masses (primary’s and satellite’s) multiplied by the gravitational constant and is the magnitude of r.

One can conveniently solve this equation [2] by introducing a new dependent (quaternion) variable such that

(6)

and a new independent (scalar) variable (called the eccentric anomaly) via

(7)

where is the (regular) time, is a positive constant (we choose its value to equal the orbit’s semimajor axis), and is a scalar (equal to the magnitude of ).

The resulting so-called Kustaanheimo-Stiefel equation (which we do not need to quote here) is solved by the following quaternion function of t.

The orbit’s eccentricity is . In terms of (6), this is the usual ellipse in the – plane.

The previous solution can be further rotated by three Euler angles.

The following two lines prove that we have thus solved (5).

At the same time, this solution is fully general, capable of describing any elliptical orbit.

Since , (7) can be easily integrated, yielding the so-called Kepler equation:

(8)

where is the value of at .

Orbit Determination

Given Initial Values

Suppose that we are given the location and velocity of a satellite at time (the initial conditions). The corresponding orbital elements can then be computed by the following program.

Thus, for example, having an Earth-like primary with (we use 10,000 km as our unit distance and 1 hour as a unit of time), and a satellite currently at with velocity , we find the values of its orbital elements.

Solving Kepler’s Equation

To find the location and the velocity of the satellite some time (say 20 hours) later, we first need to solve (8) for , knowing that, at , the value of was 2.14254.

This tells us (as a byproduct) that over six full orbits will be completed during that time.

The new location and velocity can then be computed.

Note that by setting to zero, the above procedure returns the original location and velocity, thus checking its correctness.

Given Two Locations

Establishing values of the orbital elements when two locations and the corresponding time interval are given (these are called boundary conditions) is slightly more difficult and not fully analytic (some numerical root-finding is necessary), but is in general more accurate.

This gives the solution in terms of , , , and , and the three Euler angles.

We can verify its correctness.

Given Three or More Positions

This is the classical problem, where we know the observer’s location and the direction toward the satellite (given by two spherical angles and —astronomers call this the satellite’s position) at several (at least three) distinct times. There are many techniques of various degrees of sophistication for computing the orbital elements based on such data; here we can demonstrate only the simplest, brute-force approach.

Corresponding to each observation, the “theoretical” direction (a unit 3D vector) to the satellite is computed as a function of its (as yet unspecified) orbital elements. This is then subtracted from the observed direction (converted into a unit vector), and the sum of squares of the three components of this difference is cumulatively totaled for all observations. Similarly, differences between the actual and “theoretical” time based on Kepler’s equation are computed at each observation. Then, the sum of squares of consecutive differences of these quantities, each also divided by the corresponding time interval for proper scaling, is added to the previous total. All we need to do then is to minimize the resulting quantity, thereby reducing the discrepancy between the observed and theoretical values to near 0.

The function cts converts two spherical angles to a unit vector.

The following routine does the rest.

The first argument of classical, d, is a list of the observed satellite’s directions; the second argument, L, is a list of the corresponding observer’s locations; and the last argument, T, is a list of times at which these observations were taken.

Here is an example.

The final sum of squares is roughly , reflecting the fact that our input values had only a two-digit accuracy. This is clearly not a typical situation (astronomical observations tend to be substantially more precise); nevertheless, it is of some interest to note that when such inaccurate data is used, the resulting orbital elements may be up to 50% off the mark. Things improve dramatically when we reach a typical “medieval” accuracy of three digits (the resulting error is usually not more than 1-2%); with five or more digits we can normally expect the answers to have a comparable (five digit or better) accuracy, and the resulting sum of squares to be smaller than .

Further Challenges

In this article, we have discussed the motion of a satellite not perturbed by any extra forces, save the gravitational pull of a perfectly round primary. We know that in the real world there are always additional forces due to

• other celestial bodies (the remaining planets, when the “satellite” is the Earth; the Sun, when the “satellite” is the Moon orbiting the Earth, etc.),
• the Earth not having a perfectly round shape (the so called “oblateness” perturbations, most relevant for artificial satellites),
• other, less important forces: atmospheric drag, tidal forces, Kepler shear (being hit by co-orbiting small objects), radiation pressure, and so on.

To incorporate such forces in our solution is more difficult and will be left to a follow-up article.

References

[1]	L. G. Taff, Celestial Mechanics: A Computational Guide for the Practitioner, New York: Wiley, 1985.
[2]	J. Vrbik, “A Novel Solution to Kepler’s Problem,” European Journal of Physics, 24(6), 2003 pp. 575-583. www.df.uba.ar/users/sgil/physics_paper_doc/papers_phys/mechan/kepler_prob.pdf.
J. Vrbik, “Solving Kepler’s Problem,” The Mathematica Journal, 2011. dx.doi.org/doi:10.3888/tmj.13-12.

About the Author

Jan Vrbik
Department of Mathematics, Brock University
500 Glenridge Ave., St. Catharines
Ontario, Canada, L2S 3A1
jvrbik@brocku.ca