Chapter-3-Orbital-Position-as-a-Funct 2014 Orbital-Mechanics-for-Engineeri PDF

Preview

Summary

Orbital mechanics mathematical background for solving 2 body equations....

Description

CHAPTER

Orbital Position as a Function of Time

3

CHAPTER OUTLINE 3.1 Introduction ..................................................................................................................................145 3.2 Time since periapsis .....................................................................................................................145 3.3 Circular orbits (e [ 0)...................................................................................................................147 3.4 Elliptical orbits (e < 1)...................................................................................................................147 3.5 Parabolic trajectories (e [ 1)........................................................................................................163 3.6 Hyperbolic trajectories (e > 1)........................................................................................................165 3.7 Universal variables .......................................................................................................................173 Problems.............................................................................................................................................183 Section 3.2 .........................................................................................................................................183 Section 3.4 .........................................................................................................................................184 Section 3.5 .........................................................................................................................................186 Section 3.6 .........................................................................................................................................186 Section 3.7 .........................................................................................................................................186

3.1 Introduction In Chapter 2, we found the relationship between position and true anomaly for the two-body problem. The only place time appeared explicitly was in the expression for the period of an ellipse. Obtaining position as a function of time is a simple matter for circular orbits. For elliptical, parabolic, and hyperbolic paths, we are led to the various forms of Kepler’s equation relating position to time. These transcendental equations must be solved iteratively using a procedure like Newton’s method, which is presented and illustrated in this chapter. The different forms of Kepler’s equation are combined into a single universal Kepler’s equation by introducing universal variables. Implementation of this appealing notion is accompanied by the introduction of an unfamiliar class of functions known as Stumpff functions. The universal variable formulation is required for the Lambert and Gauss orbit determination algorithms in Chapter 5. The road map of Appendix B may aid in grasping how the material presented here depends on that of Chapter 2.

3.2 Time since periapsis The orbit formula, r ¼ ðh2 =mÞ=ð1 þ ecos qÞ, gives the position of body m2 in its orbit around m1 as a function of the true anomaly. For many practical reasons, we need to be able to determine the position Orbital Mechanics for Engineering Students. http://dx.doi.org/10.1016/B978-0-08-097747-8.00003-7 Copyright Ó 2014 Elsevier Ltd. All rights reserved.

145

146

CHAPTER 3 Orbital Position as a Function of Time

of m2 as a function of time. For elliptical orbits, we have a formula for the period T (Eqn (2.82)), but we cannot yet calculate the time required to fly between any two true anomalies. The purpose of this section is to come up with the formulas that allow us to do that calculation. _ which The one equation we have that relates true anomaly directly to time is Eqn (2.47), h ¼ 2rq, can be written as dq h ¼ 2 r dt Substituting r ¼ ðh2 =mÞ=ð1 þ ecos qÞ we find, after separating variables, dq m2 dt ¼ h3 ð1 þ ecos qÞ2 Integrating both sides of this equation yields  m2  t  tp ¼ 3 h

Zq 0

dw ð1 þ e cos wÞ2

(3.1)

where the constant of integration pt is the time at periapsis passage, where by definition q ¼ 0. tp is the sixth constant of the motion that was missing in Chapter 2. The origin of time is arbitrary. It is convenient to measure time from periapsis passage, so we will usually set pt ¼ 0. In that case we have m2 t¼ h3

Zq 0

dw ð1 þ e cos wÞ2

(3.2)

The integral on the right may be found in any standard mathematical handbook, such as Beyer (1991), in which we find Z

! pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ffiffiffiffiffiffiffiffiffiffiffi 2  b2 sin x 1 x a  b b dx a ¼ 2atan1 ðb < aÞ tan  aþb a þ bcos x 2 ða þ bcos xÞ2 ða2  b2 Þ3=2

(3.3)

Z

  1 1 dx x 1 3x ¼ ðb ¼ aÞ þ tan tan 2 2 6 ða þ bcos xÞ2 a2 2

(3.4)

Z

2

3 ! pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p ffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffi x 2 2 b þ a b  a þ tan 7 dx 1 6b b  a sin x  aln ¼ p ffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffi 2x 5 ðb > aÞ (3.5) 4 a þ bcos x b þ a  b  a tan 2 ða þ bcos xÞ2 ðb2  a2 Þ3=2

3.4 Elliptical orbits (e < 1)

3.3 Circular orbits (e [ 0) If e ¼ 0, the integral in Eqn (3.2) is simply

Z

147

q

dw, which means

0

t¼

h3 q m2 3

3

Recall that for a circle (Eqn (2.62)), r ¼ h2/m. Therefore h3 ¼ r 2 m 2 , so that 3

r2 t ¼ p ffiffiffi q m

3p ffiffiffi Finally, substituting the formula (Eqn (2.64)) for the period T of a circular orbit, T ¼ 2pr2 = m, yields

t¼

q T 2p

or q¼

2p t T

The reason that t is directly proportional to q in a circular orbit is simply that the angular velocity 2p/T is constant. Therefore, the time Dt to fly through a true anomaly of Dq is (Dq/2p)T. Because the circle is symmetric about any diameter, the apse linedand therefore the periapsisd can be chosen arbitrarily. D r t=0 C, F

P

Apse line

FIGURE 3.1 Time since periapsis is directly proportional to true anomaly in a circular orbit.

3.4 Elliptical orbits (e < 1) Set a ¼ 1 and b ¼ e in Eqn (3.3) to obtain Zq 0

dw ð1 þ ecos wÞ2

¼

1 3

ð1  e2 Þ 2

"

2tan

1

! # pffiffiffiffiffiffiffiffiffiffiffiffiffi r ffiffiffiffiffiffiffiffiffiffiffi q 1e e 1  e2 sin q tan  1þe 2 1 þ ecos q

148

CHAPTER 3 Orbital Position as a Function of Time

Therefore, Eqn (3.2) in this case becomes ! # " pffiffiffiffiffiffiffiffiffiffiffiffiffi r ffiffiffiffiffiffiffiffiffiffiffi e 1  e2 sin q q 1e m2 1 1  tan t¼ 2tan 3 2 1 þ ecos q h3 1þe ð1  e2 Þ 2 or Me ¼ 2tan

1

! pffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffi q 1e e 1  e2sin q tan  1 þ ecos q 1þe 2

(3.6)

where Me ¼

3 m2  1  e2 2 t 3 h

(3.7)

2p t T

(3.8)

Me is called the mean anomaly. The subscript e reminds us that this is for an ellipse and not for parabolas and hyperbolas, which have their own “mean anomaly” formulas. Equation (3.6) is plotted in Figure 3.2. Observe that for all values of the eccentricity e, Me is a monotonically increasing function of the true anomaly q. From 3Eqn (2.82), the formula for the period T of an elliptical orbit, we have m2 ð1  e2 Þ2 =h3 ¼ 2p=T, so that the mean anomaly can be written much more simply as Me ¼

Mean anomaly, Me

e=0 e = 0.2 e = 0.5 e = 0.8 e = 0.9

0

FIGURE 3.2 Mean anomaly vs true anomaly for ellipses of various eccentricities.

3.4 Elliptical orbits (e < 1)

149

The angular velocity of the position vector of an elliptical orbit is not constant, but since 2p radians are swept out per period T, the ratio 2p/T is the average angular velocity, which is given the symbol n and called the mean motion, n¼

2p T

(3.9)

In terms of the mean motion, Eqn (3.5) can be written simpler still, Me ¼ nt The mean anomaly is the azimuth position (in radians) of a fictitious body moving around the ellipse at the constant angular speed n. For a circular orbit, the mean anomaly M e and the true anomaly q are identical. It is convenient to simplify Eqn (3.6) by introducing an auxiliary angle E called the eccentric anomaly, which is shown in Figure 3.3. This is done by circumscribing the ellipse with a concentric auxiliary circle having a radius equal to the semimajor axis a of the ellipse. Let S be that point on the ellipse whose true anomaly is q. Through point S we pass a perpendicular to the apse line, intersecting the auxiliary circle at point Q and the apse line at point V. The angle between the apse line and the radius drawn from the center of the circle to Q on its circumference is the eccentric anomaly E. Observe that E lags q from periapsis P to apoapsis Að0  q < 180  Þ, whereas it leads q from A to P ð180  q < 360 Þ. To find E as a function of q, we first observe from Figure 3.3 that, in terms of the eccentric anomaly, OV ¼ acos E, whereas in terms of the true anomaly, OV ¼ ae þ rcos q. Thus, acos E ¼ ae þ rcos q B Q S b a

A

a

O

D

FIGURE 3.3 Ellipse and the circumscribed auxiliary circle.

E ae

r

F

V

P

150

CHAPTER 3 Orbital Position as a Function of Time

Using Eqn (2.72), r ¼ að1  e2 Þ=ð1 þ ecos qÞ, we can write this as   a 1  e2 cos q acos E ¼ ae þ 1 þ ecos q Simplifying the right-hand side, we get cos E ¼

e þ cos q 1 þ ecos q

(3.10a)

Solving this for cos q we obtain the inverse relation, cos q ¼

e  cos E ecos E  1

(3.10b)

2 2 Substituting Eqn (3.10a) into the trigonometric identity sin E þ cos E ¼ 1 and solving for sin E yields pffiffiffiffiffiffiffiffiffiffiffiffiffi 1  e2sin q sin E ¼ (3.11) 1 þ ecos q

Equation (3.10a) would be fine for obtaining E from q, except that, given a value of cos E between 1 and 1, there are two values of E between 0 and 360 , as illustrated in Figure 3.4. The same comments hold for Eqn (3.11). To resolve this quadrant ambiguity, we use the following trigonometric identity: 1  cos E sin2 E=2 1  cos E 2 tan ¼ ¼ ¼ 2 1 þ cos E cos E=2 1 þ cos E 2 2 2E

(3.12)

From Eqn (3.10a) 1  cos E ¼

1  cos q ð1  eÞ and 1 þ ecos q

1 þ cos E ¼

1 þ cos q ð1 þ eÞ 1 þ ecos q

1 cos E

0

E EI

EIV 2

2

1

FIGURE 3.4 For 0 < cos E < 1, E can lie in the first or fourth quadrant. For 1 < cos E < 0, E can lie in the second or third quadrant.

3.4 Elliptical orbits (e < 1)

151

Therefore, tan2

E 1  e 1  cos q 1  e q ¼ ¼ tan2 2 2 1 þ e 1 þ cos q 1 þ e

where the last step required applying the trig identity in Eqn (3.12) to the term (1  cos q)/(1 þ cos q). Finally, therefore, we obtain rffiffiffiffiffiffiffiffiffiffiffi q 1e E tan (3.13a) tan ¼ 2 1þe 2 or

E ¼ 2tan

1

! r ffiffiffiffiffiffiffiffiffiffiffi q 1e tan 2 1þe

(3.13b)

Observe from Figure 3.5 that for any value of tan (E/2), there is only one value of E between 0 and 360 . There is no quadrant ambiguity. Substituting Eqns (3.11) and (3.13b) into Eqn (3.6) yields Kepler’s equation, Me ¼ E  esin E

(3.14)

This monotonically increasing relationship between mean anomaly and eccentric anomaly is plotted for several values of eccentricity in Figure 3.6. Given the true anomaly q, we calculate the eccentric anomaly E using Eqn (3.13). Substituting E into Kepler’s formula, Eqn (3.14), yields the mean anomaly directly. From the mean anomaly and the period T we find the time (since periapsis) from Eqn (3.5), t¼

tan

Me T 2p

(3.15)

E 2

0

FIGURE 3.5 For any value of tan (E/2), there is a corresponding unique value of e in the range 0 to 2p.

E

152

CHAPTER 3 Orbital Position as a Function of Time

e = 1.0 e = 0.8 e = 0.6

Mean anomaly, Me

e = 0.4 e = 0.2 e=0

0 Eccentric anomaly, E

FIGURE 3.6 Plot of Kepler’s equation for an elliptical orbit.

On the other hand, if we are given the time, then Eqn (3.15) yields the mean anomaly M e. Substituting Me into Kepler’s equation, we get the following expression for the eccentric anomaly: E  esin E ¼ Me We cannot solve this transcendental equation directly for E. A rough value of E might be read from Figure 3.6. However, an accurate solution requires an iterative, “trial and error” procedure. Newton’s method, or one of its variants, is one of the more common and efficient ways of finding the root of a well-behaved function. To find a root of the equation f(x) ¼ 0 in Figure 3.7, we estimate it to be xi and evaluate the function f(x) and its first derivative f 0 (x) at that point. We then extend the tangent to the curve at f(xi) until it intersects the x-axis at xiþ1 , which becomes our updated estimate of the root. The intercept xiþ1 is found by setting the slope of the tangent line equal to the slope of the curve at xi, that is, f 0 ðxi Þ ¼

0  f ðxi Þ xiþ1  xi

from which we obtain xiþ1 ¼ xi 

f ðxi Þ f 0 ðxi Þ

(3.16)

The process is repeated, using xiþ1 to estimate xiþ2, and so on, until the root has been found to the desired level of precision. To apply Newton’s method to the solution of Kepler’s equation, we form the function f ðEÞ ¼ E  esin E  Me

3.4 Elliptical orbits (e < 1)

153

f

Root xi

x xi+1

Slope = f(xi)

df dx

x= xi

FIGURE 3.7 Newton’s method for finding a root of f(x) ¼ 0.

and seek the value of eccentric anomaly that makes f(E) ¼ 0. Since f 0 ðEÞ ¼ 1  ecos E for this problem, Eqn (3.16) becomes Eiþ1 ¼ Ei 

Ei  esin Ei  Me 1  ecos Ei

(3.17)

n ALGORITHM 3.1 Solve Kepler’s equation for the eccentric anomaly E given the eccentricity e and the mean anomaly M e. See Appendix D.11 for the implementation of this algorithm in MATLAB. 1. Choose an initial estimate of the root E as follows (Prussing and Conway, 1993). If Me < p, then E ¼ M e þ e/2. If Me > p, then E ¼ M e  e/2. Remember that the angles E and Me are in radians. (When using a hand held calculator, be sure that it is in the radian mode.) 2. At any given step, having obtained Ei from the previous step, calculate f ðEi Þ ¼ Ei  esin Ei  Me and f 0 ðEi Þ ¼ 1  ecos Ei . 3. Calculate ratioi ¼ f ðEi Þ=f 0 ðEi Þ. 8 4. If jratioi j exceeds the chosen tolerance (e.g., 10 ), then calculate an updated value of E, Eiþ1 ¼ Ei  ratioi Return to Step 2. 5.

If jratioi j is less than the tolerance, then accept Ei as the solution to within the chosen accuracy.

n

154

CHAPTER 3 Orbital Position as a Function of Time

EXAMPLE 3.1 A geocentric elliptical orbit has a perigee radius of 9600 km and an apogee radius of 21,000 km, as shown in Figure 3.8. Calculate the time to fly from perigee P to a true anomaly of 120 .

B

120°

A

P

Earth 21,000 km 9600 km

FIGURE 3.8 Geocentric elliptical orbit.

Solution Before anything else, let us find the primary orbital parameters e and h. The eccentricity is readily obtained from the perigee and apogee radii by means of Eqn (2.84), e¼

ra  rp 21;000  9600 ¼ 0:37255 ¼ ra þ rp 21;000 þ 9600

(a)

We find the angular momentum using the orbit equation, 9600 ¼

h2 1 0 h ¼ 72;472 km2 =s 398;600 1 þ 0:37255 cosð0Þ

With h and e, the period of the orbit is obtained from Eqn (2.82),

T ¼

3 3   72; 472 h 2p 2p p ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ 18;834 s pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ 2 2 2 2 398;600 m 1  0:37255 1e

Equation (3.13a) yields the eccentric anomaly from the true anomaly, tan

E ¼ 2

r ffiffiffiffiffiffiffiffiffiffiffiffi r ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1e 1  0:37255 q 120 tan ¼ tan ¼ 1:1711 0 E ¼ 1:7281 rad 2 2 1þe 1 þ 0:37255

Then Kepler’s equation, Eqn (3.14), is used to find the mean anomaly, Me ¼ 1:7281  0:37255 sin 1:7281 ¼ 1:3601 rad

(b)

3.4 Elliptical orbits (e < 1)

155

Finally, the time follows from Eqn (3.12), t¼

1:3601 Me 18;834 ¼ 4077s ð1:132 hÞ T ¼ 2p 2p

EXAMPLE 3.2 In the previous example, find the true anomaly at three hours after perigee passage.

Solution Since the time (10,800 s) is greater than one-half the period, the true anomaly must be greater than 180 . First, we use Eqn (3.12) to calculate the mean anomaly for t ¼ 10,800 s. Me ¼ 2p

10;800 t ¼ 2p ¼ 3:6029 rad T 18;830

(a)

Kepler’s equation, E  esin (E) ¼ Me (with all angles in radians), is then employed to find the eccentric anomaly. This transcendental equation will be solved using Algorithm 3.1 with an error tolerance of 106. Since M e > p, a good starting value for the iteration is E 0 ¼ M e  e/2 ¼ 3.4166. Executing the algorithm yields the following steps: Step 1: E0 ¼ 3:4166 f ðE0 Þ ¼ 0:085124 f 0 ðE0 Þ ¼ 1:3585 ratio ¼

0:085124 ¼ 0:062658 1:3585

jratioj > 106 ; so repeat: Step 2: E1 ¼ 3:4166  ð0:062658Þ ¼ 3:4793 f ðE1 Þ ¼ 0:0002134 f 0 ðE1 Þ ¼ 1:3515 ratio ¼

0:0002134 ¼ 1:5778  104 1:3515

jratioj > 106 ; so repeat:

156

CHAPTER 3 Orbital Position as a Function of Time

Step 3:   E2 ¼ 3:4793  1:5778  104 ¼ 3:4794 f ðE2 Þ ¼ 1:5366  109 f 0 ðE2 Þ ¼ 1:3515 ratio ¼

1:5366  109  1:137  109 1:3515

jratioj < 106 ; so accept E ¼ 3:4794 as the solution: Convergence to even more than the desired accuracy occurred after just two iterations. With this value of the eccentric anomaly, the true anomaly is found from Eqn (3.13a) to be

tan

q ¼ 2

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ffiffiffiffiffiffiffiffiffiffiffiffi 1 þ 0:37255 1þe 3:4794 E ¼ 8:67210 q ¼ 193:2 tan tan ¼ 2 2 1  0:37255 1e

EXAMPLE 3.3 Let a satellite be in a 500 km by 5000 km orbit with its apse line parallel to the line from the earth to the sun, as shown in Figure 3.9. Find the time that the satellite is in the earth’s shadow if: (a) the apogee is toward the sun and (b) the perigee is toward the sun.

Solution We start by using the given data to find the primary orbital parameters, e and h. The eccentricity is obtained from Eqn (2.84), e¼

ra  rp ð6378 þ 5000Þ  ð6378 þ 500Þ ¼ ¼ 0:24649 ra þ rp ð6378 þ 5000Þ þ ð6378 þ 500Þ

b

c r

RE To the sun

A

P Earth

d

FIGURE 3.9 Satellite passing in and out of the earth’s shadow.

a

(a)

3.4 Elliptical orbits (e < 1)

157

The orbit equation can then be used to find the angular momentum

rp ¼

h2

1 1 h2 0 h ¼ 58;458km2 =s 0 6878 ¼ 1 þ ecosð0Þ 398;600 1 þ 0:24649 m

(b)

The semimajor axis may be found from Eqn (2.71),

a¼

h2 1 1 58;4582 ¼ 9128 km ¼ m 1  e2 398;600 1  0:246492

(c)

or from the fact that a ¼ ðrp þ ra Þ=2. The period of the orbit follows from Eqn (2.83), 2p 2p T ¼ pffiffiffi a 3=2 ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 91283=2 ¼ 8679:1 s ð2:4109 hÞ m 398;600

(d)

(a) If the apogee is toward the sun, as in Figure 3.9, then th...