Two-body orbits: example from initial scalars


I want to work out a problem similar to the previous one — but different.

Instead of being given the position and velocity vectors, we will be given their magnitudes, i.e. the distance (from the center of mass of the primary) and the speed — and the flight path angle.

These are enough to let us determine the scalar orbit; and that is what I want to show you how to do.

The problem has one new feature: it is measured in what are called “canonical units”. I will explain them before the end of the post, but I want to focus on the problem itself first. At this point, we need to know only two things about canonical units, in addition to their names.

DU stands for “distance unit” and TU stands for “time unit”. The two facts we need for this problem are that

  1. mu = 1;
  2. DU = the radius of the primary.

(Oh, okay… TU is a time unit such that an object in a circular orbit at 1 DU would have a speed of 1 DU/TU. And for the sun, we take DU = 1 astronomical unit, rather than the radius of the sun; mu = 1 still follows, because it only depends on circular orbit speed = 1 at 1 DU.)

The example

Here’s the problem (#1.19 from pp. 47-48 of BMW, i.e. “Fundamentals of Astrodynamics” by Bate, Mueller, and White).

BMEWS (the ballistic missile early warning system) detects an unidentified object with the following parameters:

altitude = .5 DU
speed = \sqrt{2/3}\ DU/TU
flight-path angle = 30{}^{\circ}

Is it possible that this object is a space probe intended to escape the Earth, an Earth satellite or a ballistic missile?

Alright then, let’s see what kind of orbit we have.

Oh, we are given altitude instead of distance.

But 1 DU = radius of primary, so we just add 1 to the altitude to get the distance… therefore the distance is

r = 1.5.

We now have all the numbers we are given:

ps=\{v\to 0.816497,r\to 1.5,\mu \to 1,\phi \to 30 {}^{\circ}\}

(Let me remark that those arrows specify rules, so that the symbol “r”, for example, remains just a symbol, while r/.ps would be a number, namely 1.5 . When I had trouble with the ship versus fort ballistics problem, I was using the wrong rule in one place.)

We can compute the energy immediately, because it involves only the scalars speed and distance. Here’s the equation…

{\cal{E}}=\frac{v^2}{2}-\frac{\mu }{r}

and we plug in, getting…

{\cal{E}} =-0.333333\ .

The energy is negative, so it’s a bound orbit. More to the point, this is not a space probe leaving the earth.

Next, the flight path angle is exactly what we need for computing h, and therefore p. We have the scalar equation for h, where \phi\ is the flight path angle (defined in this post, but we could very well take this equation as its definition):

h = r v \cos (\phi )\

and we plug in, getting…

h = 1.06066 .

We know that p depends only on h:

p = \frac{h^2}{\mu }

and we plug in, getting…

p = 1.125 .

The equation for the scalar orbit is

r=\frac{p}{e \cos (\nu )+1}

so we still need to find the eccentricity e.

What to do?

The key is that the energy is constant, so we can get additional information by computing the energy at a special point on the orbit, namely at periapse (closest approach). It is customary to introduce the semi-major axis of the orbit (the usual for an ellipse, infinite for a parabola, negative for a hyperbola), but we don’t need it. If we did get it, we would find one relationship between energy \cal{E}\ and semi-major axis a, and then another relationship between a, p, and e. Then, having a and p, we would compute the eccentricity e.

Rather than flip open a shopping list of equations and grab a couple, however, I’m going to fall back on the principle — the following calculation is effectively how we would derive the pair of equations involving a. In a sense, what I’m doing is eliminating a — but I’m doing it up front, rather than afterwards. I’m avoiding introducing it at all.

Another way to look at this is: I am using the tool I’ve got. Mathematica® will do this for me, rather easily. I know that there is a relationship between energy \cal{E}\ , p, and e, and I am going to find it.

Here are the scalar equations for energy, (the magnitude of) the angular momentum, and the orbit (I have chosen not to set \mu = 1\ ):

\left\{\mathcal{E}=\frac{v^2}{2}-\frac{\mu }{r},h=r\ v \sin (\gamma), r=\frac{p}{e \cos (\nu )+1}\right\}\

(Yes, this time I used the angle \gamma\ between the r and v vectors, rather than the flight path angle \phi\ .)

At periapse, r is rp, \nu\ is 0, and the angle \gamma\ between r and v is 90{}^{\circ}: we get

\left\{\mathcal{E}=\frac{v^2}{2}-\frac{\mu }{\text{rp}},h=\text{rp}\  v,\text{rp}=\frac{p}{e+1}\right\}\ .

Mathematica is quite happy to eliminate speed v…

h^2=\text{rp}\ (2\ \text{rp}\ \mathcal{E}+2 \mu )\land p=(e+1) \text{rp}

(that \land\ (“and”) separates two equations) and then to eliminate rp…

(2 e+2) p\ \mu +2 p^2 \mathcal{E}=\left(e^2+2 e+1\right) h^2

and then to replace h^2 by p…

(2 e+2) p\ \mu +2\ p^2 \mathcal{E}=\left(e^2+2\ e+1\right) p

and that’s all we really need: we have just gotten an equation relating the energy \cal{E}\ , the parameter p, and the eccentricity e. Given any two, we can find the third.

Yes, I know it’s not pretty. I could clean it up, but why bother? I’m not trying to display formulas, I’m trying to compute e. (In retrospect, I probably should have done some minimal cleaning up; instead, I checked it another way.)

Given \cal{E}\ and p, we can use that equation to find e. Set \mu = 1\ , plug in \cal{E}\ and p:

1.125 (2 e+2)-0.84375=1.125 \left(e^2+2 e+1\right)\ .

Take the non-negative solution:

{e -> 0.5}

With an eccentricity less than 1, this is an ellipse. From the general scalar orbit…

r=\frac{p}{e \cos (\nu )+1}

we have gotten the particular case:

r=\frac{1.125}{0.5 \cos (\nu )+1}

Finally, we should get the true anomaly \nu\ from the given value of r:

\{\{\nu \to -2.0944\},\{\nu \to 2.0944\}\}\

Here it is (the gray disk is the earth); incidentally, you can confirm something I didn’t show you: the semi-major axis a = 1.5, i.e. the major axis itself is 3.

The object had a positive flight path angle (30{}^{\circ}), and that’s the point shown in red. Whatever it is, the object will have a maximum altitude (apoapse) about 1.2 DU (hence about 4800 miles altitude). That’s pretty far up there — but I have no idea how high or how fast an intercontinental missle goes.

And if it isn’t a missile, it had better change orbit when it gets to apoapse. I suppose if it doesn’t, that’s when we commit to destroying it.

(The authors were professors at the U.S. Air Force Academy; I suspect that their students were expected to know whether this was a missile or a satellite on its way to another orbit.)

But I’ve gone as far as I can. We’re done with the problem itself.

Let’s return to the ugly but computationally useful equation

(2 e+2) p \mu +2 p^2 \mathcal{E}=\left(e^2+2 e+1\right) p\ .

As I said, this is not the way the equation is usually presented. What I did was follow the usual derivation, but omit an auxiliary variable, namely the semi-major axis of the ellipse. (The following equations are also valid for hyperbolic orbits. For parabolic orbits, they reduce to \mathcal{E} = 0\ and p = 0 \cdot \infty\ .)

Usually one would see the following two equations:

{\cal{E}} = - \mu / 2a

p = a (1 - e^2)\ .

(Now I show them to you! Sheesh! And if it makes you more comfortable — it did me — I confirmed my answer for e, and got a = 1.5, using this pair of equations!)

They are certainly useful, and a whole lot prettier, and they may provide insight, or they may be easier to pick up again, and I have certainly been known to use them… but I didn’t want to get too far afield. The simple fact is that there is a relationship between \cal{E}\ , p, and e, and that is what I derived. It is more commonly expressed as these two relationships: one between \cal{E}\ and a, and one between a, p, and e.

Canonical units

DU is simple enough: the radius of the primary (or, for the sun as primary, 1 astronomical unit = the radius of the earth’s orbit around the sun — whether the object is earth, Mars, Jupiter, Voyager…).

To get TU, however, we must figure out the speed in a circular orbit; once we do that, we take a circular orbit of radius 1 DU, and find TU so that the speed is 1 DU/TU. (I can’t tell you how many times I have written DU/DT right there. Calculus is deeply embedded in my brain.)

Let me remark that canonical units are convenient because…

  • \mu = 1\ ;
  • numbers are small;
  • converting between altitude & distance is trivial;
  • code can be written independently of the primary.

Speed in a circular orbit.

Let’s return to that set of 3 equations from which we obtained the \cal{E}\ , p, e equation.

\left\{\mathcal{E}=\frac{v^2}{2}-\frac{\mu }{\text{rp}},h=\text{rp}\  v,\text{rp}=\frac{p}{e+1}\right\}

In fact, since we’re now considering a circular orbit, we have e = 0. (And, yes, rp = r is constant… let me change it to R.)

\left\{\mathcal{E}=\frac{v^2}{2}-\frac{\mu }{R},h=R\ v,R=p\right\}

Let me replace p by h^2/mu…

\left\{\mathcal{E}=\frac{v^2}{2}-\frac{\mu }{R},h=R v,R=\frac{h^2}{\mu}\right\}

Now ask Mathematica to eliminate energy… (edit: the \land\ is a logical “and”, joining the two equations being fed to Mathematica.)

h=R\ v\land \mu =R\ v^2

That second equation will do nicely… solve it for v… take the positive solution:

\left\{v\to \frac{\sqrt{\mu }}{\sqrt{R}}\right\}\ .

Nice. We see that if R = 1 and v = 1, then \mu = 1\ . That was one of the two properties I asked you to take on faith.

But what is TU?

If I simply write…

DU/TU = speed


TU = DU/speed

and then…

TU = DU / (\sqrt{\mu} / \sqrt{r}) = r \sqrt{r}/ \sqrt{\mu}\ .

Let’s try it for the earth. Here’s the gravitational constant G…

(So that’s \mu\ .)

So that’s TU for an orbit around the earth. (And BMW agrees pretty much with that number — what we disagree on, a little bit, is the radius of the earth!)

Note that I used “ct” and “cd” for those numerical values of TU and DU.

We know that the object was spotted at an altitude of about 2000 miles (DU = .5), but what is its speed? We have speed in canonical units… then miles per second… and finally miles per hour.

(Incidentally, this picture shows the use of rules in Mathematica.)

While we’re at it, let’s explicitly convert \mu\ to canonical units. We have…

\mu = (1.40812*10^{16}\ Feet^3)/Second^2

Convert feet to miles…

tm1 = (95661.6\ Mile^3)/Second^2

As we expected, \mu = 1\ .


Given distance r, speed v, and flight path angle… we found the scalar orbit as follows:

  • Compute energy \cal{E}\ from r and v;
  • Compute h and p from r, v, and flight path angle;
  • Compute e from p and \cal{E}\ ;

Once having the orbit, we used the given distance to find the two angles \nu\ corresponding to that distance — and the flight path angle tells us which of the two angles \nu\ is right because a positive flight path angle means the object is outbound from periapse.

Finally, we saw what canonical units actually were… and converted the object’s speed to miles/hour… and confirmed numerically that \mu = 1\ for the earth.

OK, I’ve earned the Memorial Day barbecue I’m heading off to.


Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )


Connecting to %s

%d bloggers like this: