## Two-body orbits: example from initial vectors

Edit: 2011 Apr 29, mu = GM not Gm!

## Introduction & Review

I am going to work a numerical example of the equations I worked out in the previous post. To be specific, I am going to assume we are given the position and velocity vectors for an object in orbit, at one instant of time.

Let me begin by doing what I should have done at the end of the previous post: here’s a summary of the previous post.

We took the vector differential equation

$\ddot{\vec{r}} + \frac{\mu}{r^3}\ \vec{r} = 0$

as our starting point for the two-body problem, where G is the universal gravitation constant in Newton’s equation for the gravitational force between two bodies of masses M and m:

$F = G m M / r^2$

and

$\mu = G\ (m+M)$

— although we will almost always, perhaps always, assume that mass m is negligibly small compared with mass M, and we will take

$\mu = G\ M\$.

(Edit: That’s a correction; I can’t believe I had LC m instead of UC M!)
First, we dotted that vector differential equation with the velocity vector $\dot{\vec{r}}\$

$\dot{\vec{r}} \cdot \ddot{\vec{r}} + \frac{\mu}{r^3}\ \dot{\vec{r}} \cdot \vec{r}$

and got conservation of energy:

$\frac{v^2}{2} - \frac{\mu}{r} = \cal{E}\$, a constant.

Second, we crossed that vector differential equation with the position vector $\vec{r}\$

$\vec{r} \times \ddot{\vec{r}} + \frac{\mu}{r^3} \vec{r} \times \vec{r} = 0\$.

and got conservation of angular momentum:

$\vec{r} \times \vec{v} = \vec{h}\$, a constant.

Third, we crossed the vector differential equation with the angular momentum vector $\vec{h}\$

$\ddot{\vec{r}} \times \vec{h} + \frac{\mu}{r^3}\ \vec{r} \times \vec{h} = 0\$

and got a vector orbit equation for the position and velocity on the orbit:

$\dot{\vec{r}} \times \vec{h} - \frac{\mu\ \vec{r}}{r} = \vec{B}\$.

Fourth, we dotted the position vector $\vec{r}\$ into that vector orbit equation

$\vec{r} \cdot \dot{\vec{r}} \times \vec{h} - \frac{\mu}{r}\ \vec{r} \cdot \vec{r} = \vec{r} \cdot \vec{B}\$

and got a scalar orbit equation for the position as a function of angle:

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

where the angle $\nu\$ is defined as the angle between the vector constant of integration $\vec{B}\$ and the position vector $\vec{r}\$:

$\vec{r} \cdot \vec{B} = r\ B\ cos(\nu)\$

but we can show that the angle $\nu\$ is also the angle between periapse (closest approach) and r, so we infer that the vector $\vec{B}\$ points from the primary (the origin) to periapse. (When $\nu = 0\$, we have $\vec{r}\$ pointing to periapse, and $\vec{B}\$ parallel to it, so $\vec{B}\$, which is constant, points toward periapse, too.)

Oh, and

$e = B/\mu\$ and $p = h^2 / \mu \$.

## The example

This example comes from p. 18 of “Fundamentals of Astrodynamics” by Bate, Mueller, and White (see the bibliography page).

Given position and velocity vectors at one time in some inertial coordinate system,
find specific energy & momentum. (They don’t say this is the “earth”, but it is. And they don’t say the origin is at the center of the earth, but it is.)

Incidentally, let me call this the “observer coordinate system”.

They also want to
find the flight path angle (whatever that is; I’ll you what it is shortly).

First, let’s just go to town. We are given position and velocity vectors:

$\vec{r} = (4.1852, 6.2778, 10.463)\ 10^7\ Feet$

$\vec{v} =(2.5936, 5.1872, 0)\ 10^4\ Feet/Second$

Angular momentum (actually, specific angular momentum or angular momentum per unit mass) is just the cross product

$\vec{h} = \vec{r} \times \vec{v}$

$\vec{h} = (-5.42737*10^{12},\ 2.71368*10^{12},\ 5.42737*10^{11} )\ Feet^2/Second$

For the energy, we need $\mu\$ of earth. Mathematica will happily give me the gravitational constant…

The product is… $\mu\$ = G M = 3.98735*10^14 Meter^2 Newton/Kilogram

which simplifies to… $\mu\$ = 3.98735*10^14 Meter^3/Second^2

and I convert to English units… $\mu\$ = 1.40812*10^16 Feet^3/Second^2

Then the specific energy is

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

$\cal{E}\$ = 1.57253*10^9 Feet^2/Second^2

Note that the energy is positive, so this is a hyperbolic orbit.

Last part of the question. What is this thing called the flight path angle? Here’s a picture:

This a redrawing of their figure 1.4-1, on p 17. $\gamma\$ is the angle between r and v; $\phi\$ is the flight path angle.

Because $\phi + \gamma = 90{}^{\circ}\$, we are replacing the scalar equation (for the magnitude of the cross product)

$h = r\ v\ sin(\gamma)$

by

$h = r\ v\ cos(\phi)\$.

(To be more specific, $\gamma\$ is the angle between r and v; by using its complement $\phi\$, we get to write the magnitude of h as though it were some kind of dot product — but the fact remains that $\phi\$ is not the angle between r and v, and h is not the dot product of r and v.)

They say, plausibly, that the sign of $\phi\$ will be the sign of r.v, which is also the sign of $cos(\phi)\$. I infer that if $\gamma > 90{}^{\circ}\$, then $\phi\$ will be chosen negative, so that we still have $\phi + \gamma = 90{}^{\circ}\$.

OK, we need magnitudes.

r = 1.28997*10^8 Feet

v = 57994.7 Feet/Second

h = 6.0922*10^12 Feet^2/Second

(They did show all of those magnitudes, and we have the same answers. And, just in case you’re very new at this, the magnitude of a vector is just the square root of its dot product with itself: e.g. $r = \sqrt {\vec{r} \cdot \vec{r}}\$, and that’s how I computed the magnitudes.)

Then

$cos(\phi)\$ = h / r v = 0.814345

and $\phi\ = 35.4773{}^{\circ}\$

(Once again, as we would hope, the book and I agree.)

That’s the end of the example.

At the very least, however, I want to convert the distance and speed into units I have a feeling for. Mathematica will happily convert.

Just how far from the center of the earth is this thing? Convert r to miles…

24431.2 Mile

(That’s about where a geosynchronous orbit would be.)

… and how fast is it moving? Convert v to mph…

39541.8 Mile/Hour

(Smoking!)

## the scalar orbit

But we can get so much more. In particular, I write the vector equation for the orbit in the form

$\vec{v} \times \vec{h} - \frac{\mu}{r}\ \vec{r} = \vec{B}$

and solve for B. The vector is…

$\vec{B} = (2.35843*10^{16} ,-2.09292*10^{16} , 3.40489*10^{17}) Feet^3/Second^2$

and its magnitude is

B = 3.41946*10^17 Feet^3/Second^2.

Then $p = h^2/\mu\$, $e = B/\mu\$

p = 2.63578*10^9 Feet

e = 24.2839

(That’s good: e > 1, which is another way of seeing that the orbit is hyperbolic.)

We can now write the equation of the orbit:

$r = \frac{2.63578*10^9 Feet}{1+24.2839\ Cos(\nu)}$

Now, where was this thing when we saw it? Since $\vec{r} \cdot \vec{B} = r\ B\ cos(\nu)\$, we have

$\frac{4.41099\times 10^{25} \text{Feet}^4\ \cos (\nu)} {\text{Second}^2} = \frac{3.52985\times 10^{25} \text{Feet}^4}{\text{Second}^2}$

which gives us two values for $\nu\$: ±0.643099 (radians).

Hmm. We have the magnitude of $\nu\$, but we don’t know whether we’re approaching perigee (closest approach) or leaving it.

But we have the position vector. Can’t we do a change-of-basis to get from the observer coordinate system to the orbit (local) coordinate system?

Sure we can!

But first, let’s just try drawing a picture. I want to show you what the scalar orbit equation gives us.

In the orbit coordinate system, as we saw, the orbit has polar equation

$r = \frac{2.63578*10^9 Feet}{1+24.2839\ Cos(\nu)}$

What’s the radius of the earth? (I want to show it on the drawing.)

Now here’s the drawing. It has only two pieces: the scalar equation for the orbit, and a circle whose radius is that of the earth.

We’d be home free if we could transform our initial r and v into that coordinate system.

But that’s easy.

In the oberver coordinate system, we define a new x axis using a unit vector along B; just divide the vector $\vec{B}\$ by its magnitude.

$\hat{i} = (0.0689708, -0.0612062, 0.995739)$

(In the orbit coordinate system, the x-axis always points toward periapse… but the B vector does that, so it becomes our new x axis, hence defines a new unit vector.)

We define a new z axis using a unit vector along h:

$\hat{k} = (-0.890871, 0.445435, 0.0890871)$

(The angular momentum is perpendicular to r and v, i.e. to the plane of the orbit… and we take that plane to be the xy-plane in the orbit coordinate system… so h gives us the z-axis.)

To get the third basis vector in a right-handed coordinate system, we define the new y axis using a cross-product, namely

$\hat{j} = \hat{k} \times \hat{i}\$.

Exactly what are these vectors? They are the old (observer) components of the new basis vectors.

Therefore the attitude matrix has these vectors for rows…

$att = \left(\begin{array}{ccc} 0.0689708 & -0.0612062 & 0.995739 \\ 0.44899 & 0.89322 & 0.0238048 \\ -0.890871 & 0.445435 & 0.0890871\end{array}\right)$

The transition matrix T (mapping new to old, orbit to observer) is the transpose of the attitude matrix:

$T = \left(\begin{array}{ccc} 0.0689708 & 0.44899 & -0.890871 \\ -0.0612062 & 0.89322 & 0.445435 \\ 0.995739 & 0.0238048 & 0.0890871\end{array}\right)$

and the inverse transition matrix, mapping old to new, observer to orbit) is

$T^{-1} = \left(\begin{array}{ccc} 0.0689708 & -0.0612062 & 0.995739 \\ 0.44899 & 0.89322 & 0.0238048 \\ -0.890871 & 0.445435 & 0.0890871\end{array}\right)$

Notice something?

The inverse transition matrix is the same as the attitude matrix, which is the transpose of the transition matrix. That is, the inverse of the transition matrix is also its transpose.

Oh, the transition matrix is a rotation. (We didn’t change the physical scale for the orbit coordinate system, so all we’re doing is rotating the coordinates.)

That shouldn’t be surprising, and I could have simply said so up front — but computing power is cheap…. It’s probably a good thing to see it, whether before or after computing it.

What we really want to know is the orbit coordinates for our initial position and velocity vectors $\vec{R}\$ and $\vec{V}\$. We apply the inverse trasition matrix to $\vec{R}\$ and to $\vec{V}\$

$T^{-1}\ \vec{R} = (1.03228*10^8,\ 7.73564*10^7,\ 7.45058*10^{-9}) Feet$

$T^{-1}\ \vec{V} = (-1386.06,\ 57978.1,\ 3.63798*10^{-12}) Feet/Second$

and we see that the z-components are, in fact, zero. (Please note the terms 10^-9 and 10^-12 for the z-components.) That’s reassuring, because our orbit coordinate system was chosen so that the orbit was in the xy-plane!

Let me just get 2D vectors, instead of 3D, for our 2D drawing… and drop the units:

$\vec{r} = (1.03228*10^8,\ 7.73564*10^7)$

$\vec{v} = (-1386.06,\ 57978.1)$

NOTE that v is much smaller than r. I’m going to multiply v by 1000 for the graph.

So, take the previous pictutre of the orbit… draw the r vector from the origin… and then draw (magnified 1000 times) the v vector from r.

Whatever it is, it is outbound, moving away from the earth. (That is, the correct choice for $\nu\$ was positive.)

## Summary

Given the position and velocity vectors at one time — not all that easy to get in the real world — we…

• computed the specific angular momentum $h = \vec{r} \times \vec{v}\$ as the cross product
• computed the specific energy using the velocity, distance, and $\mu\$, $\frac{v^2}{2} - \frac{\mu}{r} = \cal{E}\$
• defined the flight path angle, and computed it from $cos(\phi)\ = \frac{h}{r\ v}$
• computed $\vec{B}\$ and used it to write out the equation for the scalar orbit
• then we could draw the orbit (and the primary)
• then we constructed the transition matrix (etc.) relating the (old ) initial coordinate system and the (new) orbit coordinate system
• and we transformed the given r and v to the orbit coordinate system
and added them to the drawing.

I have a few additional remarks.

The flight path angle may seem a little odd. In the elementary theory of curves, we define a really nice coordinate system (the Frenet frame) by starting with the tangent vector — which can often be interpreted as the velocity. From it we would get the binormal and normal vectors.

By contrast, our starting point for the flight path angle is the “local vertical” — the position vector — and then the local horizontal; the velocity vector is not generally one of those two directions.

And that brings me to a second remark. The flight path angle is zero exactly when the position and velocity are perpendicular… i.e. exactly at periapse and apoapse (closest and furthest distances from the primary; yes, apoapse is infinite for parabolic and hyperbolic orbits)… so we might suspect that it has one sign as we go from peripase to apoapse, and the other sign coming the other way.

That is, we should expect to see that the flight path angle has one sign “inbound” and the other sign “outbound”.

Go ahead, draw both possibilities for an ellipse. The point is that the positive value of flight path angle implies that the object was on its way to apoapse.

I could have told that it was outbound without drawing a picture.

In addition, although I never wrote out the explicit vector equation for the orbit, the computation of $\vec{B}\$ was all we needed — I just headed straight for the scalar orbit.