Orbital Dynamics 1

3 Orbital Dynamics

3.1 The Two-Body Problem

In this section we consider a system of two particles with masses m1 and m2 which experience
only the mutual gravitation acting between them. Although only point masses are covered
by this treatment, the application to satellites orbiting the earth will be justified in a latter
section.

FI f 12 f 21
~ 6 ~ -
~ u
u
r1 
3m :

~ 1 @ m2
  r @

 -
2
 ~ @
@@
O r = r1 − r2
~ ~ ~

Newton’s Law of Gravitation gives the force exerted on m1 by m2 as

Gm1 m2
f 12 = − r = −f 21 (1)
~ r3 ~ ~
where r = r1 − r2 . Applying Newton’s Second Law to m1 and m2 gives
~ ~ ~
Gm1 m2 Gm2
m1 r̈1 = f 12 = − 3
r ⇒ r̈1 = − 3 r (2)
~ ~ r ~ ~ r ~
Gm1 m2 Gm1
m2 r̈2 = f 21 = 3
r ⇒ r̈2 = r (3)
~ ~ r ~ ~ r3 ~
Subtracting (3) from (2) gives
µ ∆
r̈ = − 3 r , µ = G(m1 + m2 ) (4)
~ r~
The goal of this chapter is to provide a complete solution of this differential equation for r(t)
~
and v(t) = ṙ given initial conditions r(0) and v(0).
~ ~ ~ ~
Addition of (2) and (3) gives
m1 r̈1 + m2 r̈2 = 0 (5)
~ ~ ~
But
(m1 + m2 )rc = m1 r1 + m2 r2 (6)
~ ~ ~
Orbital Dynamics 2

defines the location of the centre of mass, rc . Therefore, r̈c = 0 which shows that the centre
~ ~ ~
of mass moves at constant velocity. As such, we may take the origin of our inertial frame to
be located at that point. Hence, rc = 0 and without loss in generality,
~ ~
m1 r1 + m2 r2 = 0 (7)
~ ~ ~

An Equivalent Problem
Now, define
∆ m1 m2 ∆
M = m1 + m2 , m= (8)
m1 + m2
. .
For m1 ≪ m2 , M = m2 and m = m1 . Multiplying both sides of (4) by m gives
GMm Gm1 m2
mr̈ = − 3
r=− r = f 12
~ r ~ r3 ~ ~
Hence, the motion of m1 with respect to m2 is the same as if m2 were fixed and m1 were
replaced by the reduced mass m.
Energy
The kinetic energy of the system is given by

T = 21 m1 ṙ1 · ṙ1 + 12 m2 ṙ2 · ṙ2 (9)

~ ~ ~ ~
Using (7),
r1 = −(m2 /m1 )r2 = −(m2 /m1 )(r1 − r)
~ ~ ~ ~
Therefore,
m2
r1 = r
~ m1 + m2~
m1 m1
r2 = − r1 = − r
~ m2~ m1 + m2~
Substituting into (9) gives

1 m1 m22 1 m2 m21
T = 2
ṙ · ṙ + ṙ · ṙ
(m1 + m2 )2~ ~ 2 (m1 + m2 )2~ ~
= 12 mṙ · ṙ
~ ~
Hence, the interpretation introduced above would give the correct kinetic energy.
The potential energy is given by
Gm1 m2 GMm µm
V =− =− =−
r r r
Notice that f 12 = −∇V , i.e., the gravitational force is derivable from a potential function
~ ~
since it is conservative.
Orbital Dynamics 3

The total energy (per unit of reduced mass) is given by

µ µ
E = 21 ṙ · ṙ − = 21 v 2 − , v 2 = ṙ · ṙ (10)
~ ~ r r ~ ~
The time derivative is given by
µ
Ė = ṙ · r̈ + 2 ṙ
~ ~ r
µ µ
= −ṙ · r 3 + 3 r ṙ
~ ~r r
But
d d
r · ṙ = 12 (r · r) = 12 (r 2 ) = r ṙ (11)
~ ~ dt ~ ~ dt
Therefore, Ė = 0 and the total energy is conserved.
Angular Momentum
The angular momentum of the system per unit mass is given by
∆
h = r × v , v = ṙ (12)
~ ~ ~ ~ ~
The time derivative is given by
ḣ = ṙ × ṙ + r × r̈
~ ~ ~ ~ ~
µ
= −r × r 3 = 0 (13)
~ ~r ~
Therefore h is a constant vector. Since it is perpendicular to r and v, the motion is confined
~ ~ ~
to a plane.
Solution
Since the motion is planar and can treated as solving for the motion of the reduced mass
with respect to m2 , we introduce polar coordinates:
hm
*


12
22 ~ 
6
~
AKAθ 21 r
 
A *
~
u 6θ - 11
A
m2 ~

The position and angular velocity can be expressed as

0
   
r
T T  T T 
r = F2 r = F2  0  , ω = F2 ω = F2  0  (14)
 
~ ~ ~ 0 ~ ~ ~ θ̇
Orbital Dynamics 4

The velocity and acceleration of m are given by:

ṙ = F T2 (ṙ + ω × r)
~ ~ 
ṙ


= F T2  θ̇r  (15)
 
~ 0

and

r̈ = F T2 (r̈ + 2ω × ṙ + ω̇ × r + ω × ω × r)
~ ~
r̈ − r θ̇2
 

= F T2  r θ̈ + 2ṙ θ̇  (16)
 
~ 0

The motion equation (4) becomes

−µ/r 2
 
µ T
r̈ = − 3 F 2 r = F T2  0 (17)
 
r

~ ~ ~ 0

Therefore,

r̈ − r θ̇2 = −µ/r 2 (18)

r θ̈ + 2ṙ θ̇ = 0 (19)

The second of these can be integrated to give

d 2
(r θ̇) = 2r ṙθ̇ + r 2 θ̈ = 0
dt
or
r 2 θ̇ = h (20)
where h is a constant (it is easily seen using (14) and (15) to be the magnitude of h = r × v,
~ ~ ~
the angular momentum).
In order to solve Eq. (18) we make the substitution r = 1/u. Therefore,
h
θ̇ = = hu2
r2
dr du 1 1 du du
ṙ = · = − 2 u̇ = − 2 θ̇ = −h
du dt u u dθ dθ
2 2
dṙ d u du
r̈ = θ̇ = −h 2 θ̇ = −h2 u2 2
dθ dθ dθ
Making these substitutions in (18) gives

d2 u
−h2 u2 − h2 u3 = −µu2
dθ2
Orbital Dynamics 5

or
d2 u µ
2
+u= 2 (21)
dθ h
The solution to this linear equation can be written as
µ
u= + C cos(θ − θ0 ) (22)
h2
where C and θ0 are constants of integration. Therefore,
1
r = µ
+ C cos(θ − θ0 )
h2
h2 /µ
= 2 (23)
1 + h µC cos(θ − θ0 )

When θ = θ0 + 2nπ, r takes on its minimum value. Without loss in generality, we shall take
θ0 = 0 and assign t = t0 when r = rmin.
Defining
h2 h2 C
p= , e= (24)
µ µ
Eq. (23) can be written as
p
r= (25)
1 + e cos θ
Note that
ṙ = ep θ̇ sin θ/(1 + e cos θ)2
Hence, ṙ = 0 when θ = 0. Let us evaluate the energy at this point in the trajectory. From
(15), v 2 = v · v = ṙ 2 + r 2 θ̇2 . Hence,
~ ~
2
1 2 2 1h µ
E = 2
r θ̇ − µ/r = 2 2
−
r r
1 2 2
= 2
h u − µu
µ µ
= 21 h2 ( 2 + C)2 − µ( 2 + C)
h h
2
µ 2µC µ2
= 12 h2 ( 4 + 2 + C 2 ) − 2 − µC
h h h
2
µ
= 21 h2 C 2 − 21 2 (26)
h
Solving for C gives s
2E µ2
C= +
h2 h4
Using (24), the solution for e becomes:
s s
h2 C 2Eh2 2Ep
e= = 1+ 2 = 1+ (27)
µ µ µ
Orbital Dynamics 6

Geometry of Conic Sections

Eq. (25) is the polar equation of a conic section where e is the eccentricity and p is called the
semilatus rectum. Hence, m2 occupies one of the foci of a conic section. It shall be referred
to as the occupied focus. Given the constants h and E, p and e are uniquely determined by
(24) and (27) respectively.
The solution of (25) can be classified according to the value of e:

e=0 → circle (E < 0)

0≤e<1 → ellipse (E < 0)
e=1 → parabola (E = 0)
e>1 → hyperbola (E > 0)

The orbits of interest to us are ellipses and circles wherein the trajectory of m1 forms a closed
path.

@ p/e
@
a b @a p sP r/e
@ r

@
θ
@s
AK
*

s
ae ae 
 
F′ O a(1−e) a 1−e
F e
b

directrix
 -
2a

Ellipse Geometry

Elliptical Orbits
Let us restrict ourselves to the case of an elliptical orbit (0 ≤ e < 1). The minimum value
of r occurs when θ = 0 in Eq. (25):
p
rmin = rp =
1+e
The maximum value of r occurs when θ = π in Eq. (25):
p
rmax = ra =
1−e
The point at which rmin occurs is called periapsis (perigee for an earth orbit). Correspond-
ingly, rmax occurs at apoapsis (apogee for an earth orbit).
Orbital Dynamics 7

The semimajor axis of an ellipse is given by

p h2
a = 12 (ra + rp ) = = (28)
1 − e2 µ(1 − e2 )
Therefore, the size of the orbit (a) and its shape (eccentricity = e) are completely specified
by the constants h and E which are functions of position r and velocity v.
~ ~
2
Using p = a(1 − e ), we can write
p a(1 − e2 )
rp = = = (1 − e)a (29)
1+e 1+e
p a(1 − e2 )
ra = = = (1 + e)a (30)
1−e 1−e
The semiminor axis of the ellipse is given by
√ √
b = a2 − a2 e2 = a 1 − e2 (31)

The dynamical constants E and h can be expressed in terms of the geometric constants a
and e using (24) and (27):
√ q
h = pµ = aµ(1 − e2 ) (32)
(e2 − 1)µ µ
E = =− (33)
2p 2a

The last relationship states that the energy only depends on the size of the orbit. A very
useful relationship can be derived by rearranging Eq. (10) and substituting (33):
q
v = 2(E + µ/r)
s 
2 1

= µ − (34)
r a
which is called the vis-viva equation. It yields the velocity in terms of the spacecraft position.
Note that the magnitude of v varies inversely with r.
Parabolic Orbits (e = 1, E = 0)

6 y
6
p r
I
@
θ
@- -x
p/2
Orbital Dynamics 8

Since r = p/(1 + cos θ), the periapsis distance is

rp = r|θ=0 = p/2

As θ → ±π, r → ∞. For the velocity, we have

µ
r q
v= 2(E + ) = 2µ/r
r
Therefore, as r → ∞, v → 0. Therefore, m1 escapes with vanishing velocity. As an exercise,
show that the parabola may be described by y 2 = −2px.
Hyperbolic Orbits (e > 1, E > 0)
QQ 


Q 
Q
Q Q 
Q
Q
b′

Q
Q
δ -

Q

Q 


Q
Q 
Q Q 
Q Q  θ∞ 
Q γ Q

 
Qt
Q JJ
Q
t
F  - Q F′
 Q
 a′ Q
 Q
 Q
 Q
 Q
 Q
 Q
 QQ

By rewriting the polar equation of the hyperbola in Cartesian form, it can be shown that
the semimajor axis is given by
p
a′ = −a = 2
e −1
The periapsis distance is again
p
rp = = a′ (e − 1)
1+e
As r → ∞, θ → θ∞ = cos−1 (−1/e). Hence,

γ = π − θ∞ = cos−1 (1/e)
δ = π − 2γ = 2θ∞ − π = 2 sin−1 (1/e)

We also have

b′ = (a′ + rp ) sin γ
s
1
= a′ e 1 − 2
e
′
√
= a e2 − 1
p
= √ 2 = ib
e −1
Orbital Dynamics 9

q
Using the vis-viva equation, v = (2µ/r) − µ/a, we have that as r → ∞,
q
v → v∞ = µ/a′ = hyperbolic excess speed

The value of the energy,

E = 12 vp2 − (µ/rp ) = 21 v∞
2

gives the periapsis velocity s

2 +
2µ
vp = v∞
rp

Kepler’s Laws
Johannes Kepler was the first to state the following based on an empirical study of Tycho
Brahe’s observational data:

I. The orbit of each planet is an ellipse with the sun at one focus.
II. The radius vector drawn from the sun to a planet sweeps out equal areas in
equal times.
III. The squares of the periods of the planets are proportional to the cubes of
the semimajor axes of the respective orbits.

These were subsequently derived mathematically by Newton using his Law of Gravitation
and his Laws of Motion. We have followed this route and proven Kepler’s First Law. The
Second and Third Laws are proven below.
Kepler’s Second Law
Consider the diagram below which shows the motion of m in time ∆t. The area swept out
is given by

∆A = 1
2
r · r∆θ + O(∆θ)2
1 2
= 2
r ∆θ + O(∆θ)2

Therefore,
dA dθ h
= 21 r 2 = (a constant) (35)
dt dt 2
This demonstrates Kepler’s Second Law which states that dA/dt is constant.
Orbital Dynamics 10

r(θ + ∆θ) ∆A

 
 r(θ)
 ∆θ 
@
I@
 I θ
@
 @

Kepler’s Third Law

Let T denote the period of the orbit. Since the areal rate is constant,
dA
T = A = πab (area of an ellipse)
dt
Hence,
dA
T = πab/( )
dt
2πab
=
h √
2πa2 1 − e2
= q
aµ(1 − e2 )
q
= 2π a3 /µ (36)

where we have used the fact that dA/dt = h/2, Eq. (32) for h, and Eq. (31) for b. The last
result states that T 2 is proportional to a3 which is a statement of the Third law.
The mean angular rate is defined by
2π rad q
n= = µ/a3
T

3.2 The Measurement of Time and Position

In this section, we would like to relate time to the motion of a body in an elliptical orbit.
Mathematically, what is θ(t − t0 )? Since r 2 θ̇ = h, the rate of change of the true anomaly,
θ̇ = h/r 2 , is not constant unless r is constant (circular orbit).
Eccentric Anomaly
Consider the construction of the auxiliary circle below. It is centred at the centre of the
ellipse and has radius a. The eccentric anomaly E is defined as shown in the figure.
Orbital Dynamics 11

Ac
P ′

BB Ae
 B
 AK B
P A B
b a AB
 rAB
 AKA AH
YB θ
E A ABHH
O A F
 ae -

True and Eccentric Anomaly

Mean Anomaly
The mean anomaly M is defined by

∆ t − t0
M = 2π (37)
T
Since Ṁ = 2π/T is the average angular rate, the mean anomaly represents the angle covered
if the body moved constantly at the average rate. The key to solving the time and distance
problem is to relate of θ, E, and M since M brings in time in a simple way.
From Kepler’s II Law, the area swept out in time ∆t is given by

∆A = (h/2)∆t (38)

Let Ae (t − t0 ) represent the area swept out in time (t − t0 ) and AT be the area of the ellipse.
By virtue of (37),
Ae t − t0
=
AT T
Using the definition of M, (37), we can write

t − t0 Ae
M = 2π = 2π (39)
T AT
An ellipse can be interpreted as the projection of a circle with all lines perpendicular to the
major axis equally foreshortened. Therefore, Ae is the projection of the area Ac shown in
the previous figure.
Hence,
Ae b
= (40)
Ac a
Orbital Dynamics 12

since all vertical elements of Ae are shortened by a factor of b/a. But from the figure,

Ac = [area of circular section subtended by angle E, radius a]

− [area of a triangle with base ae, height a sin E]
E
= (πa2 ) − 1 (ae)(a sin E)
2π 2
a2 E a2 e
= − sin E
2 2
a2
= (E − e sin E) (41)
2
Combining (40) and (41) gives
ba
Ae = (E − e sin E)
2
But AT = πab (the area of an ellipse) and hence we can write the mean anomaly defined by
(39) as
Ae
M = 2π = E − e sin E (42)
AT
which is called Kepler’s Equation. Given the definitions of M, (37), and T , (36), this can
also be written as
t − t0 µ
r
M = 2π = (t − t0 ) = E − e sin E (43)
T a3
which relates time to the eccentric anomaly.
We shall now try to relate θ and E. From the previous figures
r sin θ b
=
a sin E a
Therefore, r sin θ = b sin E. But from the two-body solution, r = p/(1 + e cos θ). Hence,

p sin θ a(1 − e2 )
r sin θ = = sin θ = b sin E
1 + e cos θ 1 + e cos θ
√
But, b = a 1 − e2 and therefore,
√
1 − e2 sin θ
= sin E (44)
1 + e cos θ
which is the desired result. It is also possible to show that

r cos θ = a(cos E − e)
2
a(1 − e )
⇒ cos θ = a(cos E − e)
1 + e cos θ
Therefore,
e + cos θ
= cos E (45)
1 + e cos θ
Orbital Dynamics 13

The application of trigonometric identities to (44) and (45) yields:

s
1+e
tan(θ/2) = tan(E/2) (46)
1−e
The latter relationship allows the correct quadrant to be determined. Equations (46) and
(43) give the complete solution for θ as a function of (t − t0 ). In summary:

(i) Given the elapsed time since periapsis, (t − t0 ), calculate M using (37).
(ii) Solve Kepler’s Equation, (43), for E.
(iii) Using (46), calculate θ.

The solution of Kepler’s Equation is equivalent to solving the nonlinear equation

∆
g(E) = E − e sin E − M = 0

Newton’s method provides an iterative scheme which converges quite rapidly:

Ek+1 = Ek − g(Ek )/g ′(Ek ) , k = 0, 1, 2, · · ·

where
g ′ (E) = dg/dE = 1 − e cos E
A suitable initial guess for E is E0 = M (usually). The algorithm terminates when |Ek+1 −
Ek | < ǫ or |g(Ek )| < δ where δ and ǫ are appropriately chosen small, positive numbers.
Parabolic and Hyperbolic Orbits
For a parabola, the time of flight can be written as
s " #
p3 1 θ 1 θ
t − t0 = tan + tan3
µ 2 2 6 2
In the hyperbolic case, E and M have imaginary values and we define the hyperbolic eccentric
anomaly and the hyperbolic mean anomaly by E ′ = −iE and M ′ = iM, respectively. The
relationships analogous to (43) and (46) are
µ
s
′ ′ ′
M = e sinh E − E = (t − t0 )
(a′ )3
s
θ e+1 E′
tan = tanh
2 e−1 2

3.3 The Classical Orbital Elements

The two-body problem has been written as
µ
r̈ = − 3 r
~ r~
Orbital Dynamics 14

which consists of three second order differential equations. Equivalently,

ṙ = v , r(0) = r0 (47)
~ ~ ~ ~
µ
v̇ = − 3 r , v(0) = v0 (48)
~ r~ ~ ~
which is six first order differential equations. Hence, the solution is completely specified by
the six initial conditions {r0 , v0 }. It is more convenient to have six other constants.
~ ~
Classical Orbital Elements
We have noted that the orbit lies in a plane. For an ellipse, the size and shape are specified
by a and e. A third parameter is required to temporally locate the satellite in the orbit. The
time of perigee passage, t0 , is used. Hence, {a, e, t0 } completely specify the motion within
the orbital plane. We need three parameters to specify the orientation of the orbital plane.
Two parameters are used to orient the plane of the orbit and one is needed to locate the
orbit within that plane.
Inertial Reference Frame
We require an inertial frame, F i , which is used as a reference frame for locating the orbit.
~
In the case of the earth F i = [i1 i2 i3 ]T is usually the geocentric equatorial coordinate
~ ~ ~ ~
system. The unit vectors {i1 , i2 } lie in the equatorial plane and i3 points in the direction
~ ~ ~
of the geographical north pole. The vector i1 points in the direction of Υ, Aries. Υ is the
~
direction of the vernal equinox (in the direction of the sun along a line connecting the earth
and sun on the first day of spring when the sun crosses the equator).

Heliocentric Ecliptic Coordinate System

Orbital Dynamics 15

Geocentric Equatorial Coordinate System

Orbit Orientation
Define the point in the orbit where the body crosses the equatorial plane from south to north
as the ascending node. A line from the centre of the earth through the ascending node is
the line of nodes. Denote the unit vector in this direction by n.
~
The angle between i1 and n is Ω, the right ascension of the ascending node. It takes on
~ ~
values of 0 ≤ Ω < 2π. In the case of the earth it is also called the longitude of the ascending
node.
The angle between the equatorial plane and the orbital plane measured at the ascending
node is i, the inclination. It takes on values in the range 0 ≤ i ≤ π.
The angle measured in the plane of the orbit between n and perigee is ω, the argument of
~
perigee. It takes on values of 0 ≤ ω < 2π.
The angles {Ω, i, ω} specify the orientation of the orbit with respect to an inertial frame.
Collectively, {a, e, Ω, i, ω, t0 } are the classical orbital elements.

The Classical Orbital Elements

Orbital Dynamics 16

Orbital Reference Frame

A reference frame F o = [o1 o2 o3 ]T is now attached to the orbit. It shares its origin with F i ,
~ ~ ~ ~ ~
i.e., at the occupied focus, and o1 is chosen to point in the direction of perigee. Hence, ω is
~
the angle between n and o1 . The vector o3 points in the direction of the angular momentum,
~ ~ ~
o3 = h/h, and hence is normal to the orbital plane. The vector o2 lies in the orbital plane
~ ~ ~
and completes a right-handed coordinate system.
The orientation of F o with respect to F i can be specified by the rotation matrix
~ ~
Coi = F o · F Ti
~ ~
Refering to the previous figure, we recognize Ω, i, and ω as a 3-1-3 Euler sequence describing
the orientation of F o with respect to F i . Hence,
~ ~
Coi = C3 (ω)C1 (i)C3 (Ω)
   
cω sω 0 1 0 0 cΩ sΩ 0
=  −sω cω 0   0 ci si   −sΩ cΩ 0 
   

0 0 1 0 −si ci 0 0 1
 
cΩ cω − sω ci sΩ s Ω cω + ci sω cΩ si sω
=  −cΩ sω − cω ci sΩ −sΩ sω + cω ci cΩ si cω  (49)
 

sΩ si −si cΩ ci

Therefore, {a, e, Ω, i, ω, t0 } are a complete set of orbital elements. They completely specify
the orbit and hence, are equivalent to knowing {r0 , v0 }.
~ ~

3.4 Determination of the Orbital Elements

Problem: Given r and v at any time t, what are the orbital elements?
~ ~
Solution:
Step 1. Determine the angular momentum and the components of o3 in F i :
~ ~
h = r × v , o3 = h−1 h
~ ~ ~ ~ ~
Step 2. Calculate the energy:


 r = (r · r)1/2
E= 1 2
v − µ/r , ~ ~ 1/2
2
 v = (v · v)

~ ~
Step 3. Determine: s
µ h2 a−p
a=− , p= , e=
2E µ a
Orbital Dynamics 17

Step 4. Since i is the angle between i3 and o3 ,

~ ~
cos i = i3 · o3 = [0 0 1]T hh−1 = h3 /h
~ ~
where h = F Ti h. Since 0 ≤ i ≤ π, this determines i uniquely.
~ ~
Step 5. Since n lies in the equatorial plane it has the form
~
n = cos Ω i1 + sin Ω i2 (50)
~ ~ ~
But n is a unit vector perpendicular to i3 and h. Therefore
~ ~ ~
i3 × h
n= ~ ~ (51)
~ |i3 × h|
~ ~
where h = F Ti h. Eqs. (50) and (51) allow Ω to be uniquely determined.
~ ~
Step 6. Determine the true anomaly corresponding to {r, v}. Since
~ ~
p p 1
r= ⇒ cos θ = −
1 + e cos θ re e
Since 0 ≤ θ < 2π, we need an expression for sin θ. Differentiating the above gives
dr pe sin θ he sin θ
ṙ = θ̇ = θ̇ = (52)
dθ (1 + e cos θ)2 p
But r · v = r ṙ. Hence,
~ ~
pr · v
her sin θ
r · v = r ṙ = ⇒ sin θ = ~ ~
~ ~ p her

Step 7. The angle between n and r measured in the orbital plane is ω + θ. Threfore,
~ ~
n · r = nr cos(ω + θ) = r cos(ω + θ)
~ ~
Therefore,  
ω = cos−1 (n · r)/r − θ (53)
~ ~
Note that ω + θ < π if r3 > 0 and ω + θ > π if r3 < 0.
Step 8. Knowing θ, we can solve for E using (46), and then (43) can be used to solve for
t0 .
Steps 1-8 determine {a, e, Ω, i, ω, t0 }.
Orbital Dynamics 18

3.5 Position and Velocity as a Function of Time

We can now provide a complete solution of the two-body problem. Given {r(0), v(0)},
~ ~
expressions for {r(t), v(t)} can be generated.
~ ~
From the previous analysis, we may assume that the orbital elements are known (i.e., cal-
culated using {r(0), v(0)}. The position can be expressed in F 0 as
~ ~ ~
r = (r cos θ)o1 + (r sin θ)o2 (54)
~ ~ ~

o2 r
~ 6 ~ θ
I
@
- @
o1
~

Therefore,
v = ṙ = (ṙ cos θ − r θ̇ sin θ)o1 + (ṙ sin θ + r θ̇ cos θ)o2 (55)
~ ~ ~ ~
From Eq. (52),
√
ṙ = (he/p) sin θ , h= µp
s
µ
= e sin θ (56)
p

Using this expression and r θ̇ = h/r, (55) becomes

s
µ
 
v= −(sin θ)o1 + (cos θ + e)o2
~ p ~ ~
Given the time t and the orbital elements, we can solve for E from Kepler’s Equation:
µ
r
(t − t0 ) = E − e sin E
a3
Then θ is given by Eq. (46). Hence,
 
r cos θ
p
r(t) = F Ti Cio  r sin θ  , r = (57)
 
~ ~ 1 + e cos θ
0
Orbital Dynamics 19

Also, Cio = CToi is the transpose of the matrix presented in Eq. (49). The velocity is given
by  
− sin θ s
T  µ
v(t) = F i Cio  cos θ + e  (58)

~ ~ p
0