Lambert’s theorem

Posted on 26 October 2022 by John

At the start of this series we pointed out that a conic section has five degrees of freedom, and so claims to be able to determine an orbit from less than five numbers must be hiding information somewhere. That is the case with Lambert’s theorem which reportedly determines an orbit from two numbers.

Lambert’s theorem says that the time between two observations of an object in an elliptical orbit can be determined by the sum of the distances to the two observations, the length of the chord connecting the two observations, and the semi-major axis of the ellipse. From this one can determine the equation of the ellipse.

Let r₁ be the distance to the first observation and r₂ the distance to the second observation, and C the length of the chord between the two observations. We don’t need to know r₁ and r₂ individually but only their sum r₁ + r₂. And if we know the angle between the two observations we can use the law of cosines to find C. However we get there, we have two lengths: r₁ + r₂ and C. We also know a, the length of the semi-major axis. This gives us three lengths.

The observations are taken from the body being orbited, so we implicitly know one of the foci of the ellipse. That’s four pieces of information. The fifth is the mass of the object being orbited.

Lambert’s theorem divides into three cases: elliptic, parabolic, and hyperbolic.

Elliptical orbits

For an object in an elliptic orbit, Lambert’s theorem becomes

$\sqrt{\frac{\mu}{a^3}} (t_2 - t_1) = \alpha - \beta - (\sin \alpha - \sin\beta)$

where μ is the standard gravitational parameter, equal to GM where G is the gravitational constant and M is the mass of the body orbited. The angles α and β are given by

$\begin{align*} \sin\frac{\alpha}{2} = \frac{1}{2}\left( \frac{r_1 + r_2 + C}{a} \right)^{1/2} \\ \sin\frac{\beta}{2} = \frac{1}{2}\left( \frac{r_1 + r_2 - C}{a} \right)^{1/2} \\ \end{align*}$

Hyperbolic orbits

The hyperbolic case is very similar, replacing sine with hyperbolic sine.

$\sqrt{\frac{\mu}{a^3}} (t_2 - t_1) = \gamma - \delta - (\sin \gamma - \sin\delta)$

where

$\begin{align*} \sinh\frac{\gamma}{2} = \frac{1}{2}\left( \frac{r_1 + r_2 + C}{a} \right)^{1/2} \\ \sinh\frac{\delta}{2} = \frac{1}{2}\left( \frac{r_1 + r_2 - C}{a} \right)^{1/2} \\ \end{align*}$

Parabolic orbits

The parabolic case was worked out by Newton:

$t_2 - t_1 = \frac{1}{6\sqrt{mu}}\left( (r_1 + r_2 + C)^{3/2} - (r_1 + r_2 - C)^{3/2}\right)$

Gibbs’ method of determining an orbit

Posted on 25 October 2022 by John

Josiah Willard Gibbs

Josiah Willard Gibbs (1839–1903) was prominent American scientist at a time when America had yet to produce many prominent scientists. I first heard of him via Gibbs phenomenon and later by attending one of the Gibbs lectures at an AMS meeting.

Gibbs came up with a method of determining an orbit from three observations. As discussed earlier, a conic section has five degrees of freedom. So how can we determine a conic section, i.e. a two-body problem orbit, from three observations?

Mars Orbiter

As a warm-up to Gibbs’ method, let’s look back at the post on India’s Mars Orbiter Mission. There we determined an orbit by just two observations. What determined the missing three degrees of freedom?

The post on the Mars Orbiter Mission (MOM) didn’t use two arbitrary observations but two very special observations, namely the nearest and furthest approaches of the probe. Because the distances were known relative to Mars, we know the position of one of the foci of the orbit, i.e. the center of mass of Mars. And because the two observations are on opposite sides of the orbit, we know where the other focus is. And because these observations were along the major axis of the ellipse, we know the orientation of the ellipse. All together we know five facts about the orbit.

Gibbsian method

Gibbs method starts with three facts about the orbit, the position of a satellite at three points in its orbit. These are not simply three points on the ellipse but three measurements taken from the body that the satellite is orbiting. That means we know the position of one of the foci, i.e. the body being orbited. Gibbs’ method also assumes we know the mass of body being orbited, so that’s our fifth piece of information.

So with the Mars Orbiter Mission we had two special observations. The “specialness” of these observations provided two more pieces of data, and the fact that the measurements were taken from a focus of the orbit gave us the fifth piece of information.

With Gibbs’ method, we have an extra observation, but none of the observations are special. So we’re up one piece of information from observation but down two pieces of information related to specialness. The missing piece of information is supplied by a physical quantity, the mass of the object orbited.

Radar observation

There is a variation on Gibbs’ method that uses only one observation. It is based on radar observation, not optical observation, and includes the velocity of the object as well as its position. So two partial derivatives, rates of change in two perpendicular directions, replace two of the observations. We still need five bits of information to determine five degrees of freedom.

Preliminary determination

The methods this series of posts use the minimum amount of information algebraically necessary to determine an orbit. In practice, these methods are used for preliminary orbit determination, i.e. they give an approximate result that could be bootstrapped to obtain a more accurate solution.

Under ideal circumstances more further observations would be redundant, but in practice having more observations lets measurement errors cancel each other out to some extent. This observation was a milestone in the development of modern statistics: the problem of determining orbits from multiple observations lead to regression theory.

Notes

You can find the details of how to determine an orbit from three observations in Fundamentals of Astrodynamics by Bate, Mueller, and White. The book also covers determining an orbit from one radar observation with derivatives.

The next post in this series looks at Lambert’s theorem for determining an orbit from two observations (plus three other pieces of information).

Five points determine a conic section

Posted on 25 October 2022 by John

This post is the first in a series looking at determining an orbit. Lambert’s theorem is often summarized by saying you can determine an orbit from two observations. This statement isn’t true without further assumptions, assumptions I plan to make explicit.

A solution to the two-body problem is an orbit given by a conic section, and the general equation of a conic section in the plane is

So conic sections have five degrees of freedom: if you know five out of the six coefficients A, B, C, D, E, and F then the equation above determines the sixth coefficient. And if you know five points on a conic section, there is an elegant way to find an equation of the conic. Given points (x_i, y_i) for i = 1, …, 5, the following determinant yields an equation for the conic section.

$\begin{vmatrix} x^2 & xy & y^2 & x & y & 1 \\ x_1^2 & x_1 y_1 & y_1^2 & x_1 & y_1 & 1 \\ x_2^2 & x_2 y_2 & y_2^2 & x_2 & y_2 & 1 \\ x_3^2 & x_3 y_3 & y_3^2 & x_3 & y_3 & 1 \\ x_4^2 & x_4 y_4 & y_4^2 & x_4 & y_4 & 1 \\ x_5^2 & x_5 y_5 & y_5^2 & x_5 & y_5 & 1 \\ \end{vmatrix} = 0$

This means that five observations are enough to determine a conic section, and since Keplerian orbits are conic sections, such an orbit can be determined by five observations. How do we get from five down to two?

Astronomical observations have more context than merely points in a plane. These observations take place over time. So we have not only the positions of objects but their positions at particular times. And we know that the motion of an object in orbit is constrained by Kepler’s laws. In short, we have more data than (x, y) pairs; we have (x, y, t) triples plus physical constraints.

We also have implicit information, and future posts in this series will make this implicit information explicit. For example, suppose you’re planning a trajectory to send a probe to Mars. The path of the probe will essentially be an orbit around the sun. You can determine this orbit by two observations: the position of Earth when the probe leaves and the position of Mars when it arrives. This orbit is not simply an ellipse passing through two points. It is an ellipse with one focus at the sun. I intend to unpack this in a future post, or series of posts, making implicit data explicit.

When I write a series of blog posts, the post don’t always come out consecutively. I expect I’ll write about other topics in between posts in this series.

Update: The next post in the series considers Gibbs’ method of determining an orbit from three observations (plus two other pieces of information). The post after that is about Lambert’s theorem.

Elliptical orbit example: Mars Orbiter Mission

Posted on 22 October 2022 by John

This post will look at India’s first interplanetary mission, Mars Orbiter Mission, to illustrate points in recent posts.

Mars Orbiter Mission insignia

As suggested by the logo, the probe had a very eccentric orbit of Mars with periareion 3,812 km and apoareion 80,384 km. We can derive everything else from these numbers. [1]

Peri-what?! As mentioned in footnote 2 of this post the nearest and furthest points of a satellite’s orbit go by various names depending on what the satellite is orbiting. The terms perigee and apogee for Earth satellites are more familiar; the terms periareion and apoareion are the Martian counterparts. The names for the closest and furthest approaches of satellites of Jupiter would be perijove and apojove. The “areion” suffix refers to Ares, the Greek name for Mars.

So the Mars Orbiter Mission (MOM) was much closer to Mars on its nearest approach than on the opposite side of its orbit, 21 times closer. That means its tangential velocity was 21 times faster at periareion than at apoareion [2].

The orbit was an ellipse with major axis 3,812 km + 80,384 km = 84,196 km. The semi-major axis was a = 42,098 km. The center of Mars was a distance c = a – 3,812 km = 38,286 km from the center. The eccentricity was thus e = c/a = 0.909.

Using the equations from this post we find that the aspect ratio of the ellipse was 0.4158.

Orbit of Mars Orbiter Mission

Although the deviation of the orbit from circular is substantial, the distance of Mars from the center of the orbit is dramatic. More on this here.

Using the formulas from this post we can plot the mean anomaly, true anomaly, and eccentric anomaly for MOM.

True anomaly and eccentric anomaly of Mars Orbiter Mission

The steep lines on the far left and far right show how quickly the probe is moving at nearest approach to the planet. The true anomaly is the angle to the satellite based at Mars, and this angle changes very rapidly when the probe is near the planet. Eccentric anomaly is measured from the center of the orbit and so it does not change as quickly.

The period along with either the true anomaly or eccentric anomaly is enough to know the position of the satellite at all times. We can determine the period of a satellite from Kepler’s laws:

T = 2π √(a³/μ)

where μ = GM, with G being the gravitational constant and M being the mass of Mars. Now G = 6.674 × 10^-11 N m²/kg² and M = 6.417 × 10²³ kg. We’ll need to multiply a by 1,000 to work in units of meters. This gives us a period of 262,242 seconds or about 73 hours.

***

[1] My understanding is that while perigee and apogee are measured as altitude from the Earth’s surface, peri[planet] and apo[planet] are conventionally measured from the center of the planet for other planets. That’s how I’m using the terms here. When I first wrote this post, I used a value of periapeion from Wikipedia which was less than the radius of Mars and hence impossible. In retrospect the Wikipedia article was reporting an altitude.

[2] For objects in elliptical orbits, v_peri / v_apo = r_apo / r_peri. This follows from Kepler’s law that the vector from a planet to its satellite sweeps out equal area in equal time.

Mean anomaly, true anomaly, and eccentric anomaly

Posted on 22 October 2022 by John

Orbital mechanics has a lot of arcane terminology because it has been studied for centuries. V. I. Arnold said that orbital mechanics was one of the three main sources of modern mathematics.

Mean anomaly, true anomaly, and eccentric anomaly are three ways of describing where an object is in its orbit. All would be the same if planets moved in circular orbits. And since planets nearly do move in circular orbits [1], at least in our solar system, these three anomalies will be similar. We will demonstrate this with plots.

Mean anomaly

Suppose a planet orbiting a star has period T, i.e. it takes time T for the planet to complete one orbit. Let t₀ be the time when the planet was closest to the star, i.e. the time since the planet passed through periapsis [2]. Let t be the current time. Then the mean anomaly M is simply

M = 2π(t – t₀) / T

if you’re working in radians. Change 2π to 360° if you’re working in degrees.

True anomaly

The true anomaly f is the angle with vertex at the star and sides running from the star to periapsis and from the star to the planet. (I don’t know why true anomaly is denoted f. I suppose t and T were taken.)

Eccentric anomaly

The eccentric anomaly E is sorta like the true anomaly, but with the vertex at the center of the elliptical orbit rather than at the star. I say sorta because that’s not quite right.

The eccentric anomaly is an angle with vertex at the center of the elliptical orbit. (The star is not at the center because elliptical orbits are eccentric, literally off-center.) One side does run along the major axis as with the true anomaly, but the other side doesn’t extend from the center to the planet but rather from the center to a projection of the planet’s orbit onto a circle.

Imagine an elliptical orbit in the plane, centered at the origin, with major axis along the x-axis. Now draw a circle around the ellipse with radius equal to the semi-major axis of the ellipse. So the circle touches the ellipse on the two ends, apsis and periapsis, but otherwise extends above and below the ellipse. Take a point P in the orbit and push it out vertically to a point P‘ on the circumscribed circle. That is, keep the x coordinate of P the same, but push the y coordinate up if it’s positive, or down if it’s negative, to get a point on the circle.

Equations

Kepler’s equation relates mean anomaly and eccentric anomaly:

M = E – e sin E

where M is the mean anomaly, E is the eccentric anomaly, and e is the eccentricity of the orbit. So when e is small, M is approximately equal to E.

True anomaly f is related to eccentric anomaly by

f = E + 2 arctan( β sin E / (1 – β cos E) )

where

β = e / (1 + √(1 – e²)).

I’ve written about Kepler’s equation a couple times. Here is a post about a simple way to solve Kepler’s equation for E given M, and here is a post about more efficient ways to solve Kepler’s equation.

Examples

Let’s make some plots comparing mean anomaly, eccentric anomaly, and true anomaly. Mean anomaly is proportional to time, so we’ll use it as our independent variable. Eccentric anomaly and true anomaly are roughly equal to mean anomaly, so we’ll plot the difference between these anomaly and mean anomaly.

Here’s what we get for Earth, e = 0.01671.

Note that the three anomalies never differ by more than 2 degrees.

And here’s what we get for Pluto, e = 0.25.

Now the differences between the anomalies is about 15 times greater. The shapes of the curves are roughly the same, but the curves for Earth are more symmetric about their peak and their trough than the corresponding curves for Pluto.

And here’s an extreme example: Halley’s comet. Its orbit has very high eccentricity, e = 0.967.

Newton’s method

I used Newton’s method to solve Kepler’s equation when making the plots above, using the mean anomaly as the initial guess when solving for the eccentric anomaly. I thought I might run into problems with high eccentricity, but everything worked smoothly, even for Halley’s comet. I suppose the basin of attraction for Newton’s method applied to Kepler’s equation is large, at least large enough to include mean anomaly as a starting point.

I didn’t even supply the root-finding method with a derivative functon. I used scipy.optimize.newton with just the function and starting point as arguments, letting the routine construct its own approximate derivative.

Notes

[1] The planets in our solar system move in very nearly circular orbits as explained here. The dwarf planet Pluto has a rather eccentric orbit (e = 0.25) and yet its orbit isn’t far from circular. However, the center of its orbit is far from the sun, as explained here.

[2] This point goes by many names in different contexts. The generic term for this point is called periapsis. For a satellite orbiting the earth it’s called perigee. For an object orbiting our sun it’s called perihelion. For an object orbiting our moon it’s called perilune. …

Cryptography, hydrodynamics, and celestial mechanics

Posted on 21 October 2022 by John

V. I. Arnold

Last night I was reading a paper by the late Russian mathematician V. I. Arnold “Polymathematics: is mathematics a single science or a set of arts?” and posted a lightly edited extract of it on Twitter. It begins

All mathematics is divided into three parts: cryptography, hydrodynamics, and celestial mechanics.

Arnold is alluding to the opening line to Julius Caesar’s Gallic Wars [1] which suggests he’s being a little playful. The unedited version leaves no doubt that he’s being playful or cynical.

All mathematics is divided into three parts: cryptography (paid for by CIA, KGB and the like), hydrodynamics (supported by manufacturers of atomic submarines), and celestial mechanics (financed by military and other institutions dealing with missiles, such as NASA).

I edited out the parenthetical remarks, in part edit the sentence down to a tweet, but also because when you take out the humor/cynicism he makes a valid if hyperbolic point. He goes on to back this up.

Cryptography has generated number theory, algebraic geometry over finite fields, algebra, combinatorics and computers.

Hydrodynamics has procreated complex analysis, partial differential equations, Lie groups and algebra theory, cohomology theory and scientific computing.

Celestial mechanics is the origin of dynamical systems, linear algebra, topology, variational calculus and symplectic geometry.

Arnold adds a footnote to the comment about cryptography:

The creator of modern algebra, Viète, was the cryptographer of King Henry IV of France.

Of course not all mathematics was motivated by cryptography, hydrodynamics, and celestial mechanics, but an awful lot of it was. And his implicit argument that applied math gives birth to pure math is historically correct. Sometimes pure math gives rise to applied math, but much more often it’s the other way around.

His statements roughly match my own experience. Much of the algebra and number theory that I’ve learned has been motivated by cryptography. I dove into Knuth’s magnum opus, volume 2, because I wanted to implement the RSA algorithm in C++.

I got started in scientific computing in a computational fluid dynamics (CDF) lab. I didn’t work in the lab, but my roommate did, and I went there to use the computers. That’s where I would try my hand at numerical analysis and where I wrote simulation code for my dissertation. My dissertation in partial differential equations was related (very abstractly) to fluid flow in porous media.

I didn’t know anything about celestial mechanics until I sat in on Richard Arenstorf‘s class as a postdoc. So celestial mechanics was not my personal introduction to dynamical systems etc., but Arnold is correct that these fields came out of celestial mechanics.

[1] “Gallia est omnis divisa in partes tres.” which translates “Gaul is a whole divided into three parts.”

How eccentricity matters

Posted on 20 October 2022 by John

I wrote last week that the eccentricities of planet orbits in our solar system do not effect the shape of the orbit very much. Here’s a plot of all the orbits, shifted to have the same center and scaled to have the same minor axis.

However, the planet orbits do not have a common center. The eccentricity of an orbit doesn’t affect its shape as much as its position. Eccentric orbits are off-center. That’s literally what eccentric means.

Kepler’s laws say that each planet orbits in an ellipse with the sun at one focus of that ellipse. That means the more eccentric an orbit is, the further the center of the orbit is from the sun.

If a planetary orbit has semi-major axis a and the foci are a distance c from the center, then the eccentricity e equals c/a. If we scale all the planetary orbits so that they each have a semi-major axis equal to 1, then e = c, i.e. the eccentricity equals the distance from each focus to the center.

Imagine the sun as the center of our coordinate system and that each planet’s orbit is aligned so that its major axis is along the x axis. Each orbit has the sun, i.e. the origin, as a focus, and so the center of each orbit is at –e where e is the eccentricity of that orbit. Here’s the plot above redone to fix the sun rather than to fix the center. The black dot in the center is the sun.

Here is the same plot with only Venus and Pluto.

The orbit of Venus is essentially a circle centered at the sun. The orbit of Pluto is an ellipse with major axis about 3% longer than its minor axis. It’s hard to see that the shape of Pluto’s orbit is not a circle, but it’s easy to see that the orbit is off-center, i.e. eccentric.

This post explains why moderately large eccentricities do not have much impact on the aspect ratio of an ellipse.

The image below from this post shows an ellipse with eccentricity 0.8. The two foci are far from the center, and yet the the aspect ratio is 3 : 5. It’s obviously an ellipse, but the foci are further apart than you might imagine.

The view from a Galilean moon

Posted on 18 October 2022 by John

The Galilean moons are the four largest moons of Jupiter, first observed by Galileo, contra Stigler’s law of eponymy.

This post shows what the Jovian system look like from the perspective of each of these moons, a sort of pre-Copernican perspective in a Jovian context.

The view from Io

Here’s what Europa would look like from Io. The orbital period for Europa is very nearly twice the orbital period of Io.

Here’s what Ganymede would look like from Io. The orbital period of Ganymede is very nearly 4 times that of Io.

Here’s what Callisto would look like from Io.

Here’s a combined view of the whole system, what Europa, Ganymede, Callisto, and Jupiter would look like from Io.

The view from Europa

The orbit of Io as seen from Europa is the same as the orbit of Europa as seen from Io.

Here’s what the orbit of Ganymede would look like from Europa. It’s similar to the view of Europa from Io because it is also a 2 : 1 resonance, i.e. the orbital period of Ganymede is approximately twice that of Europa.

Here’s what the orbit of Callisto would look like from Europa.

Here’s a combined view of the whole system, what Io, Ganymede, Callisto, and Jupiter would look like from Io.

The view from Ganymede

The view of Io from Ganymede is the same as the view of Ganymede from Io above. Similarly for Europa and Ganymede.

Here’s what the orbit of Callisto would look like from Ganymede. The two moons are in a 7 : 3 resonance.

Here’s a combined view of the whole system from Ganymede.

View of Callisto

The views of each of the moons from Callisto are presented above by symmetry. Here’s what the whole Jovian system would look like from Callisto.

View from Jupiter

Finally, here’s what the orbits of each of the moons looks like from Jupiter.

What if Copernicus had been a Martian?

Posted on 18 October 2022 by John

The genius of Copernicus was to realize that this plot

view of Mars, Jupiter, and Saturn from Earth

becomes this plot

View of Earth, Mars, Jupiter, and Saturn from the sun

under a clever change of coordinates.

The first plot is the apparent motion of Mars, Jupiter, and Saturn as viewed from Earth over the course of 30 years. Of course Copernicus knew of Mercury and Venus as well, but I’ve omitted them to make the plot simpler. The second plot is the view of Earth, Mars, Jupiter, and Saturn from the sun.

If Copernicus had lived on Mars, what view of the solar system would he have started with? Here’s what 30 (Earth) years of the apparent motion of Earth, Jupiter, and Saturn would look like from Mars.

View of Earth, Jupiter, and Saturn from Mars

These plots were produced by modifying the Python script given in the previous post.

How the orbit of one planet appears from another

Posted on 17 October 2022 by John

This post shows how the orbits some planets appear from other planets. I’ve give a few of my favorite examples and include a Python program you could use to create your own plots.

We will assume the planets move in circles around the sun. They don’t exactly—they don’t exactly move in ellipses either—but their orbits are much closer to circles than most people realize. More on that here.

Venus / Earth

The ratio of Earth’s orbital period to that of Venus is almost exactly 13 to 8, so the view of Venus from Earth is nearly periodic with period 8 (Earth) years.

Venus from Earth

This is also the view of Earth from Venus. As the code below shows, the view of X from Y is the same as the view of Y from X.

Jupiter / Saturn

The ratio of Saturn’s orbital period to that of Jupiter is nearly 5 to 2. If this ratio were exactly, the figure below would repeat itself every two Saturnine years. But the following plot is over 4 Saturnine years, showing that years 3 and 4 don’t exactly retrace the path of years 1 and 2.

Jupiter from Saturn

Uranus / Earth

Here’s what the orbit of Uranus looks like from Earth.

Uranus from Earth

We see a big circle of tight little circles because the ratio of the distances from the sun is large.

Mercury / Venus

Here’s what Mercury looks like from Venus.

Mercury from Venus

Jovian moons

The same math that describes planets around the sun describes moons orbiting a planet. Here is what the orbit of Ganymede looks like from Callisto, and vice versa, as they orbit Jupiter. Ganymede completes 7 orbits in the time it takes Callisto to complete 3 orbits.

Ganymede and Callisto

Python code

Here’s the code that produced the plots.

    import matplotlib.pyplot as plt
    from numpy import linspace, sin, cos, pi
    
    period = {
        'mercury' : 87.96926,
        'venus'   : 224.7008,
        'earth'   : 365.25636,
        'mars'    : 686.97959,
        'ceres'   : 1680.22107,
        'jupiter' : 4332.8201,
        'saturn'  : 10755.699,
        'uranus'  : 20687.153,
        'neptune' : 60190.03
    }
    
    dist = lambda T : T**(2/3) # Kepler
    
    def plot_orbit(planet1, planet2, periods=10):
        T1 = period[planet1]
        T2 = period[planet2]
        d1 = dist(T1)
        d2 = dist(T2)

        theta = linspace(0, 2*pi*periods, 1000)
        x = d1*cos(T2*theta/T1) - d2*cos(theta)
        y = d1*sin(T2*theta/T1) - d2*sin(theta)

        plt.gca().set_aspect("equal")
        plt.axis('off')
        plt.plot(x, y)
        plt.show()
    
    plot_orbit("venus", "earth", 8)
    plot_orbit("jupiter", "saturn", 4)
    plot_orbit("uranus", "earth", 57)
    plot_orbit("mercury", "venus", 9)

Elliptical orbits

Hyperbolic orbits

Parabolic orbits

Other posts in this series

Mars Orbiter

Gibbsian method

Radar observation

Preliminary determination

Notes

Related posts

Mean anomaly

True anomaly

Eccentric anomaly

Equations

Examples

Newton’s method

Notes

Related posts

The view from Io

The view from Europa

The view from Ganymede

View of Callisto

View from Jupiter

Related posts

Venus / Earth

Jupiter / Saturn

Uranus / Earth

Mercury / Venus

Jovian moons

Python code