Orbital mechanics

Solar Day vs Sidereal Day

Posted on by John

How long does it take the earth to complete one rotation on its axis? The answer depends on your frame of reference. A solar day is the time it takes for the sun to appear at the same position in the sky. A sidereal day is the time it takes for a distant star to appear in the same position. These are not the same.

This post will illustrate the difference between a solar day and a sidereal day. To make things a little simpler, assume the earth has a perfectly circular orbit around the sun, and that the earth’s axis of rotation is perpendicular to the orbital plane, i.e. there is no axial tilt.

Sidereal day

Imagine an astronomer Alice observing our solar system from the vantage point of a distant star. How distant? We’ll show below that it doesn’t make much difference, but we’ll assume for now that she is “infinitely” far away, which means “far enough away that we don’t have to worry about exactly how far.”

The time it takes for Alice to observe one rotation of the earth on its axis is 1 sidereal day = 24 sidereal hours. We will suppose that from Alice’s perspective it takes 360 (sidereal) days for the earth to orbit the sun.

Here comes the sun

Let’s set the origin of our coordinate system at the sun. Assume that at time t = 0 the earth is located at (1, 0) and that an observer Bob is at the bottom of a deep well on the equator looking up at the sun. The time between Bob’s observations of the sun is one solar day.

Twenty four (sidereal) hours later, the earth is located at (cos 1°, sin 1°) and Bob’s well is parallel to the x-axis, but not looking directly at the sun. He will be looking at the sun a few minutes later when the earth’s rotation brings the sun into view.

Just how long will Bob have to wait to see the sun again? About 4 minutes, because the sun is 1° away from his line of sight, and he needs the earth to turn 1°, which takes 24/360 hours, or 4 minutes. This is not exactly correct though, because the earth has moved in its orbit during that 4 minutes.

So a solar day is about 4 minutes longer than the time it takes for the earth to rotate on its axis (from Alice’s perspective). Even though we’ve made several simplifying assumptions, our estimate only differs from the exact value by about 4 seconds.

Equations

Let the distance from the center of the earth to the sun be 1 and the radius of the earth be ε. Let t be time in (sidereal) days. The position of the center of the earth as a function of time will be

(cos t°, sin t°)

and Bob’s position is

(cos t° − ε cos 360t° , sin t° − ε sin 360t°).

Bob sees the sun at time t = 0. When will he see the sun next? When the slope of the line from the center of the earth to his position equals the slope of the line from the sun to the center of the earth, i.e.

tan t° = tan 360t°

The value of ε doesn’t matter.

There are two solutions to this equation, one when Bob is facing the sun and another when he is on the opposite side of the earth from the sun. We know t ≈ 1 and the solution near 1 is on the correct side of the sun, i.e.

t = 1.002786 days = 1 day + 4.011142 minutes.

Distant stars

Now let’s imagine at that time t = 0, while Bob is looking up at the sun, on the opposite side of the earth Charlie at the bottom of another well looking up a star a distance R away, located at (R + 1, 0).

Charlie’s position as a function of t is

(cos t° + ε cos 360t° , sin t° + ε sin 360t°).

When will Charlie see his star next? When the line from the center of the earth through his position has the same slope as the line from Charlie to the star. This happens when

tan 360t° = sin t° / R,

the right side above being the tangent of the angle at R of a right triangle with base on the x-axis and hypotenuse running from the star to Charlie. As R goes to infinity, the right side goes to 0, and the solutions for t are integer numbers of days.

Charlie will see the star at time slightly less than 1, so let t = 1 − x. So we need to solve

tan 360(1 − x)° = − tan x° = − sin (1 − x)°/R.

Using the approximation sin θ ≈ θ ≈ tan θ for small angles θ (in radians) we have

−2πx ≈ −2π(1 − x)/(360R)

and so

x = 1/(360R + 1).

Now R is very large. For the nearest star, Proxima Centauri, R is about 270,000 AU. So x is on the order of 10−8 or smaller. This is why it doesn’t matter which star Alice is located near: a sidereal day is essentially the same whatever distant star you use as your reference point.

Sphere of influence

Posted on by John

Apollo 11 orbit diagram via NASA

Suppose a spaceship is headed from the earth to the moon. At some point we say that the ship has left the earth’s sphere of influence is now in the moon’s sphere of influence (SOI). What does that mean exactly?

Wrong explanation #1

One way you’ll hear it described is that the moon’s sphere of influence is the point at which the earth is no longer pulling on the spaceship, but that’s nonsense. Everything has some pull on everything else, so how do you objectively say the earth’s pull is small enough that we’re now going to call it zero? And as we’ll see below, the earth’s pull is still significant even when the spaceship leaves earth’s SOI.

Wrong explanation #2

Another explanation you’ll hear is the moon’s sphere of influence is the point at which the moon is pulling on the spaceship harder than the earth is. That’s a better explanation, but still not right.

The distance from the earth to the moon is about 240,000 miles, and the radius of the moon’s SOI is about 40,000 miles. So when a spaceship first enters the moon’s SOI, it is five times closer to the moon than to the earth.

Newton’s law of gravity says gravitational force between two bodies is proportional to the product of their masses and inversely proportional to the square of the distance. The mass of the earth is about 80 times that of the moon. So at the moon’s SOI boundary, the pull of the earth is 80/25 times as great as that of the moon, about three times greater.

Correct explanation

So what does sphere of influence mean? The details are a little complicated, but essentially the moon’s sphere of influence is the point at which it’s more accurate to say the ship is orbiting the moon than to say it is orbiting the earth.

How can we say it’s better to think of the ship orbiting the moon than the earth when the earth is pulling on the ship three times as hard as the moon is? What matters is not so much the force of earth’s gravity as the effect of that force on the equations of motion.

The motion of an object between the earth and the moon could be viewed as an orbit around earth, with the moon exerting a perturbing influence, or as an orbit around the moon, with the earth exerting a perturbing influence.

At the boundary of the moon’s SOI the effect of the earth perturbing the ship’s orbit around the moon is equal to the effect of the moon perturbing its orbit around the earth. It’s a point at which it is convenient to switch perspectives. It’s not a physical boundary [1]. Also, the “sphere” of influence is not exactly a sphere but an approximately spherical region.

The moon has an effect on the ship’s motion when it’s on our side of the moon’s SOI, and the earth still has an effect on its motion after it has crossed into the moon’s SOI.

Calculating the SOI radius

As a rough approximation, the SOI boundary is where the ratio of the distances to the two bodies, e.g. moon and earth, equals the ratio of their masses to the exponent 2/5:

r/R = (m/M)2/5.

This approximation is better when the mass M is much larger than the mass m. For the earth and the moon, the equation is good enough for back-of-the-envelope equations but not accurate enough for planning a mission to the moon. Using the round numbers in this post, the left side of the equation is 1/5 = 0.2 and the right side is (1/80)0.4 = 0.17.

Context

Everything above has been in the context of the earth-moon system. Sphere of influence is defined relative to two bodies. When we spoke of a spaceship leaving the earth’s sphere of influence, we implicitly meant that it was leaving the earth’s sphere of influence relative to the moon.

Relative to the sun, the earth’s sphere of influence reaches roughly 600,000 miles. You could calculate this distance using the equation above. A spaceship like Artemis leaves the earth’s sphere of influence relative to the moon at some point, but never leaves the earth’s sphere of influence relative to the sun.

Related posts

[1] The sphere of influence sounds analogous to a continental divide, where rain falling on one side of the line ends up in one ocean and rain falling on the other side ends up in another ocean. But it’s not that way. I suppose you could devise an experiment to determine which side of the SOI you’re on, but it would not be a simple experiment. An object placed between the earth and the moon at the SOI boundary would fall to the earth unless it had sufficient momentum toward the moon.

Lagrange’s quintic and Descartes’ rule

Posted on by John

Do fifth degree polynomial equations come up in applications? Yes, and this post will give an example.

In general the three-body problem, describing the motion of three objects interacting under gravity, does not have a closed-form solution. However, Euler and Lagrange discovered a few special cases that do have closed-form solutions. We will look at Lagrange’s quintic, an equation that came out of Lagrange’s elaboration on Euler’s solution involving three masses moving so that they remain colinear.

Lagrange’s quintic equation is

\begin{align*} (m_1 + m_2)x^5 &+ (3m_1 + 2m_2)x^4 + (3m_1 + m_2)x^3 \\ &- (m_2 + 3m_3)x^2 - (2m_2 + 3m_3)x - (m_2+m_3) = 0 \end{align*}

where m1, m2, and m3 are the masses of the three bodies and x represents the ratio of the distances from the second body to the other two bodies.

Descartes’ rule of signs says that the equation above has exactly one positive root. This is because the signs of the coefficients only change once: the first three coefficients are positive and the next three are negative. Only positive roots are physically meaningful since x represents a ratio of (unsigned) distances, so Laplace’s quintic has a unique meaningful solution. Note that this argument places no restrictions on the relative masses of the three objects.

We can set x = −w and apply Descartes’ rule to the polynomial equation in w. This tells us that the equation has 0, 2, or 4 solutions for positive w, i.e. for negative x. But this doesn’t tell us anything we wouldn’t know from more general principles. Fifth degree polynomials have five roots, and complex roots to equations with real coefficients come in conjugate pairs, so Lagrange equation either has two complex roots (and thus two negative real roots) or four complex roots (and no negative real roots).

Related posts

Artemis lunar orbit

Posted on by John

I haven’t been able to find technical details of the orbit of Artemis I, and some of what I’ve found has been contradictory, but here are some back-of-the-envelope calculations based on what I’ve pieced together. If someone sends me better information I can update this post.

Artemis is in a highly eccentric orbit around the moon, coming within 130 km (80 miles) of the moon’s surface at closest pass, and this orbit will take 14 days to complete. The weak link in this data is “14 days.” Surely this number has been rounded for public consumption.

If we assume Artemis is in a Keplerian orbit, i.e. we can ignore the effect of the Earth, then we can calculate the shape of the orbit using the information above. This assumption is questionable because as I understand it the reason for such an eccentric orbit has something to do with Lagrange points, which means the Earth’s gravity matters. Still, I image the effect of Earth’s gravity is a smaller source of error than the lack of accuracy in knowng the period.

Solving for axes

Artemis is orbiting the moon similarly to how the Mars Orbiter Mission orbited Mars. We can use Kepler’s equation for period T to solve for the semi-major axis a of the orbit.

T = 2π √(a³/μ)

Here μ = GM, with G being the gravitational constant and M being the mass of the moon. Now

G = 6.674 × 10−11 N m²/kg²

and

M = 7.3459 × 1022 kg.

If we assume T is 14 × 24 × 3600 seconds, then we get

a = 56,640 km

or 35,200 miles. The value of a is rough since the value of T is rough.

Assuming a Keplerian orbit, the moon is at one focus of the orbit, located a distance c from the center of the ellipse. If Artemis is 130 km from the surface of the moon at perilune, and the radius of the moon is 1737 km, then

c = a − (130 + 1737) km = 54,770 km

or 34,000 miles. The semi-minor axis b satisfies

b² = a² − c²

and so

b = 14,422 km

or 8962 miles.

Orbit shape

The eccentricity is c/a = 0.967. As I’ve written about before, eccentricity is hard to interpret intuitively. Aspect ratio is much easier to imaging than eccentricity, and the relation between the two is highly nonlinear.

Assuming everything above, here’s what the orbit would look like. The distances on the axes are in kilometers.

Artemis moon orbit

The orbit is highly eccentric: the center of the orbit is far from the foci of the orbit. But the aspect ratio is about 1/4. The orbit is only about 4 times wider in one direction than the other. It’s obviously an ellipse, but it’s not an extremely thin ellipse.

Lagrange points

In an earlier post I showed how to compute the Lagrange points for the Sun-Earth system. We can use the same equations for the Earth-Moon system.

The equations for the distance r from the Lagrange points L1 and L2 to the moon are

\frac{M_1}{(R\pm r)^2} \pm \frac{M_2}{r^2}=\left(\frac{M_1}{M_1+M_2}R \pm r\right)\frac{M_1+M_2}{R^3}

The equation for L1 corresponds to taking ± as − and the equation for L2 corresponds to taking ± as +. Here M1 and M2 are the masses of the Earth and Moon respectively, and R is the distance between the two bodies.

If we modify the code from the earlier post on Lagrange points we get

L1 = 54784 km
L2 = 60917 km

where L1 is on the near side of the moon and L2 on the far side. We estimated the semi-major axis a to be 56,640 km. This is about 3% larger than the distance from the moon to L1. So the orbit of Artemis passes near or through L1. This assumes the axis of the Artemis orbit is aligned with a line from the moon to Earth, which I believe is at least approximately correct.

Python code to solve Kepler’s equation

Posted on by John

The previous post looked at solving Kepler’s equation using Newton’s method. The problem with using Newton’s method is that it may not converge when the eccentricity e is large unless you start very close to the solution. As discussed at the end of that post, John Machin came up with a clever way to start close. His starting point is defined as follow.

\begin{align*} n &= \sqrt{5 + \sqrt{16 + \frac{9}{e}}} \\ M &= n \left((1-e)s + \frac{e(n^2 - 1) + 1}{6}s^3 \right) \\ x_0 &= n \arcsin s \end{align*}

The variable s is implicitly defined as the root of a cubic polynomial. This could be messy. Maybe there are three real roots and we have to decide which one to use. Fortunately this isn’t the case.

The discriminant of our cubic equation is negative, so there is only one real root. And because our cubic equation for s has no s² term the expression for the root isn’t too complicated.

Here’s Python code to solve Kepler’s equation using Newton’s method with Machin’s starting point.

    from numpy import sqrt, cbrt, pi, sin, cos, arcsin, random
    
    # This will solve the special form of the cubic we need.
    def solve_cubic(a, c, d):
        assert(a > 0 and c > 0)
        p = c/a
        q = d/a
        k = sqrt( q**2/4 + p**3/27 )
        return cbrt(-q/2 - k) + cbrt(-q/2 + k)
    
    # Machin's starting point for Newton's method
    # See johndcook.com/blog/2022/11/01/kepler-newton/
    def machin(e, M):
        n = sqrt(5 + sqrt(16 + 9/e))
        a = n*(e*(n**2 - 1)+1)/6
        c = n*(1-e)
        d = -M
        s = solve_cubic(a, c, d)
        return n*arcsin(s)    
    
    def solve_kepler(e, M):
        "Find E such that M = E - e sin E."
        assert(0 <= e < 1)
        assert(0 <= M <= pi) 
        f = lambda E: E - e*sin(E) - M 
        E = machin(e, M) 
        tolerance = 1e-10 

        # Newton's method 
        while (abs(f(E)) > tolerance):
            E -= f(E)/(1 - e*cos(E))
        return E

To test this code, we’ll generate a million random values of e and M, solve for the corresponding value of E, and verify that the solution satisfies Kepler’s equation.

    random.seed(20221102)
    N = 1_000_000
    e = random.random(N)
    M = random.random(N)*pi
    for i in range(N):
        E = solve_kepler(e[i], M[i])
        k = E - e[i]*sin(E) - M[i]
        assert(abs(k) < 1e-10)
    print("Done")

All tests pass.

Machin’s starting point is very good, and could make an adequate solution on its own if e is not very large and if you don’t need a great deal of accuracy. Let’s illustrate by solving Kepler’s equation for the orbit of Mars with eccentricity e = 0.09341.

Error in Machin's starting guess as a function of M

Here the maximum error is 0.01675 radians and the average error is 0.002486 radians. The error is especially small for small values of M. When M = 1, the error is only 1.302 × 10−5 radians.

Solving Kepler’s equation with Newton’s method

Posted on by John

Postage stamps featuring Kepler and Newton

In the introduction to his book Solving Kepler’s Equation Over Three Centuries, Peter Colwell says

In virtually every decade from 1650 to the present there have appeared papers devoted to the Kepler problem and its solution.

This is remarkable because Kepler’s equation isn’t that hard to solve. It cannot be solved in closed form using elementary functions, but it can be solved in many other ways, enough ways for Peter Colwell to write a survey about. One way to find a solution is simply to guess a solution, stick it back in, and iterate. More on that here.

Researchers keep writing about Kepler’s equation, not because it’s hard, but because it’s important. It’s so important that a slightly more efficient solution is significant. Even today with enormous computing resources at our disposal, people are still looking for more efficient solutions. Here’s one that was published last year.

Kepler’s equation

What is Kepler’s equation, and why is it so important?

Kepler’s problem is to solve

M = E - e \sin E

for E, given M and e, assuming 0 ≤ M ≤ π and 0 < e < 1.

This equation is important because it essentially tells us how to locate an object in an elliptical orbit. M is mean anomaly, e is eccentricity, and E is eccentric anomaly. Mean anomaly is essentially time. Eccentric anomaly is not exactly the position of the orbiting object, but the position can be easily derived from E. See these notes on mean anomaly and eccentric anomaly. This is because we’re using our increased computing power to track more objects, such as debris in low earth orbit or things that might impact Earth some day.

Newton’s method

A fairly obvious approach to solving Kepler’s equation is to use Newton’s method. I think Newton himself applied his eponymous method to this equation.

Newton’s method is very efficient when it works. As it starts converging, the number of correct decimal places doubles on each iteration. The problem is, however, that it may not converge. When I taught numerical analysis at Vanderbilt, I used a textbook that quoted this nursery rhyme at the beginning of the chapter on Newton’s method.

There was a little girl who had a little curl
Right in the middle of her forehead.
When she was good she was very, very good
But when she was bad she was horrid.

To this day I think of that rhyme every time I use Newton’s method. When Newton’s method is good, it is very, very good, converging quadratically. When it is bad, it can be horrid, pushing you far from the root, and pushing you further away with each iteration. Finding exactly where Newton’s method converges or diverges can be difficult, and the result can be quite complicated. Some fractals are made precisely by separating converging and diverging points.

Newton’s method solves f(x) = 0 by starting with an initial guess x0 an iteratively applies

x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}

Notice that

A sufficient condition for Newton’s method to converge is for x0 to belong to a disk around the root where

\left| \frac{f(x)f''(x)}{f'(x)^2}\right| < 1

throughout the disk.

In our case f is the function

f(x; e, M) = x - e \sin(x) - M

where I’ve used x instead of E because Mathematica reserves E for the base of natural logarithms. To see whether the sufficient condition for convergence given above applies, let’s define

g(x; e, M) = \left| \frac{e(x - e \sin x - M) \sin x}{(1 - e\cos x)^2} \right|

Notice that the denominator goes to 0 as e approaches 1, so we should expect difficulty as e increases. That is, we expect Newton’s method might fail for objects in a highly elliptical orbit. However, we are looking at a sufficient condition, not a necessary condition. As I noted here, I used Newton’s method to solve Kepler’s equation for a highly eccentric orbit, hoping for the best, and it worked.

Starting guess

Newton’s method requires a starting guess. Suppose we start by setting E = M. How bad can that guess be? We can find out using Lagrange multipliers. We want to maximize

(E-M)^2

subject to the constraint that E and E satisfy Kepler’s equation. (We square the difference to get a differentiable objective function to maximize. Minimizing the squared difference minimizes the absolute difference.)

Lagrange’s theorem tells us

\begin{pmatrix} 2x \\ -2M\end{pmatrix} = \lambda \begin{pmatrix} 1 - e \cos x \\ -1\end{pmatrix}

and so λ = 2M and

2x = 2M(1 - e\cos x)

We can conclude that

|x - M| \leq \frac{|e \cos x|}{2} \leq \frac{e}{2} \leq \frac{1}{2}

This says that an initial guess of M is never further than a distance of 1/2 from the solution x, and its even closer when the eccentricity e is small.

If e is less than 0.55 then Newton’s method will converge. We can verify this in Mathematica with

    NMaximize[{g[x, e, M], 0 < e < 0.55, 0 <= M <= Pi, 
        Abs[M - x] < 0.5}, {x, e, M}]

which returns a maximum value of 0.93. The maximum value is an increasing function of the upper bound on e, so converting for e = 0.55 implies converging for e < 0.55. On the other hand, we get a maximum of 1.18 when we let e be as large as 0.60. This doesn’t mean Newton’s method won’t converge, but it means our sufficient condition cannot guarantee that it will converge.

A better starting point

John Machin (1686–1751) came up with a very clever, though mysterious, starting point for solving Kepler’s equation with Newton’s method. Machin didn’t publish his derivation, though later someone was able to reverse engineer how Machin must have been thinking. His starting point is as follows.

\begin{align*} n &= \sqrt{5 + \sqrt{16 + \frac{9}{e}}} \\ M &= n \left((1-e)s + \frac{e(n^2 - 1) + 1}{6}s^3 \right) \\ x_0 &= n \arcsin s \end{align*}

This produces an adequate starting point for Newton’s method even for values of e very close to 1.

Notice that you have to solve a cubic equation to find s. That’s messy in general, but it works out cleanly in our case. See the next post for a Python implementation of Newton’s method starting with Machin’s starting point.

There are simpler starting points that are better than starting with M but not as good as Machin’s method. It may be more efficient to spend less time on the best starting point and more time iterating Newton’s method. On the other hand, if you don’t need much accuracy, and e is not too large, you could use Machin’s starting point as your final value and not use Newton’s method at all. If e < 0.3, as it is for every planet in our solar system, then Machin’s starting point is accurate to 4 decimal places (See Appendix C of Colwell’s book).

Lambert’s theorem

Posted on 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 r1 be the distance to the first observation and r2 the distance to the second observation, and C the length of the chord between the two observations. We don’t need to know r1 and r2 individually but only their sum r1 + r2. 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: r1 + r2 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)

Other posts in this series

Source: Adventures in Celestial Mechanics by Victor Szebehely

Gibbs’ method of determining an orbit

Posted on 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 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 (xi, yi) 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.

Related posts

Elliptical orbit example: Mars Orbiter Mission

Posted on 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 × 1023 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, vperi / vapo = rapo / rperi. This follows from Kepler’s law that the vector from a planet to its satellite sweeps out equal area in equal time.