- Email: [email protected]
A Direct Method of Deriving Perturbation Equation of Motion
CHINESE ASTRONOMY AND ASTROPHYSICS
ELSEVIER
Chinese Astronomy and Astrophysics 31 (2007) 205–210
A Direct Method of Deriving Perturbation Equation of Motion† FENG He-sheng
LI Yu-qiang
National Astronomical Observatories / Yunnan Observatory, Chinese Academy of Sciences, Kunming 650011
Abstract A direct method of deriving perturbation equation of motion is proposed. Namely, the variations in the orbital elements of an instantaneous ellipse in the sense of two-body problem under the action of three components of the perturbing force are directly derived by means of vector differentiation, with a sharp geometric image and explicit mechanical significance. The derivation of the two orbital elements, ω˙ and M˙ , is especially more concise and objective than in the traditional method. Key words: celestial mechanics
1. INTRODUCTION The derivation of perturbation equation in orbital elements is an old classical problem. The traditional method is as follows. On the basis of analytic mechanics, the relation between the variation in the canonical constants and the perturbation function is derived, then according to the correspondence of the canonical constants to the orbital elements the perturbation equation of the orbital elements is derived. Although this kind of method is of universal significance in mechanics, the process of derivation is very lengthy and tedious, and the existence of a force function for the perturbing forces has also to be identified. The direct derivation method which takes the orbital elements as the target appeared in the subsequent works [1−3] , and perturbation was no longer regarded as a conservative force, but the derivation method was based on the analytical conception of the variation of constants of the solution of differential equation and the inference thread was not clear enough. Burns[4] and Murray et al.[5] also adopted the direct method to derive the perturbation equation, †
Supported by Natural Science Foundation of Yunnan Province Received 2005–11–08; revised version 2006–01–05 A translation of Acta Astron. Sin. Vol. 48, No. 1, pp. 54–59, 2007
c 2007 Elsevier B . V. All rights reserved. 0275-1062/07/$-see front matter doi:10.1016/j.chinastron.2007.03.003
206
FENG He-sheng et al. / Chinese Astronomy and Astrophysics 31 (2007) 205–210
but the derivation process, especially the derivation of the last two orbital elements, is not sufficiently objective. In this article, the variation in the orbital elements of the instantaneous ellipse under the action of the three components of the perturbing force is directly inferred by taking advantage of the differential operation of vectors, which makes the derivation process much more concise and transparent, making the perturbation equation more clearly understood.
2. DERIVATION PROCESS We firstly establish a coordinate system which moves with the disturbed body (as shown in Fig.1), where R0 is the radial unit vector, T 0 the transverse unit vector and W 0 the normal unit vector of the orbital plane. The components of the disturbing force Fe in the above 3 directions are S, T and W , respectively. Through the coordinate rotation it is not difficult to calculate the components of R0 , T 0 and W 0 in the equatorial coordinate system O − XY Z, and the results are: ⎛ ⎞ cos u cos Ω − sin u cos i sin Ω (1) R0 = ⎝ cos u sin Ω + sin u cos i cos Ω ⎠ , sin u sin i ⎛
⎞ − sin u cos Ω − cos u cos i sin Ω T 0 = ⎝ − sin u sin Ω + cos u cos i cos Ω ⎠ , cos u cos i ⎛ ⎞ sin i sin Ω W 0 = ⎝ − sin i cos Ω ⎠ , cos i where u = ω + f is the argument of latitude.
Fig. 1
Geometry of the coordinate system moving with the disturbed body
(2)
(3)
FENG He-sheng et al. / Chinese Astronomy and Astrophysics 31 (2007) 205–210
207
For the convenience of exposition, the mass of the disturbed body is taken to be the unit mass, then the acceleration of the disturbed body is the force and the velocity is the momentum. The disturbed body is acted by the central gravitation F 0 and the perturbing force F e . We have d2 r μr = F0 + Fe = − 2 + Fe , (4) 2 dt r r r where − denotes the unit vector pointing to the center of the central body from the r disturbed body and μ is the gravitational constant. The moment of momentum of the unit mass is h, r × r˙ = h , (5) By differentiating the two sides of the above equation, we have dh μr dr dr d2 r = × + r × 2 = r × (− 3 + F e ) = r × F e . dt dt dt dt r
(6)
Thus we arrive at the conclusion: the derivative of the angular momentum of the disturbed body is equal to the moment of the perturbing force. The instantaneous ellipse at time t has dr dr 2 1 · = μ( − ) . dt dt r a
(7)
Differentiating the two sides, we have 2
dr d2 r 2 dr dr 1 da 2 1 da = μ(− 2 · + 2 ) = μ(− 3 r · + 2 ). dt dt2 r dt a dt r dt a dt
(8)
According to Eqs.(4) and (8) we have da dr 2a2 = Fe · . dt μ dt
(9)
The above equation may also be written as d −μ dr ( ) = Fe · . dt 2a dt
(10)
From the above formula we come to the conclusion that the rate of change of the total energy of the system is equal to the power of the work done by the perturbing force. In the new coordinate system (see Fig.1) we have F e = SR0 + T T 0 + W W 0 ,
(11)
r × F e = −rW T 0 + rT W 0 , h = μa(1 − e2 )W 0 , h = |h| = μa(1 − e2 ) .
(12) (13)
In the light of Eqs.(6) and (12) we have −rW T 0 + rT W 0 =
d dW 0 μa(1 − e2 )W 0 + μa(1 − e2 ) . dt dt
(14)
208
FENG He-sheng et al. / Chinese Astronomy and Astrophysics 31 (2007) 205–210
From the above equation we may obtain d μa(1 − e2 ) = rT , dt
(15)
dW 0 μa(1 − e2 ) = −rW T 0 , dt
(16)
dW 0 where and T 0 can be gotten from Eqs.(3) and (2). dt Thus, we have ⎛ di dΩ ⎞ cos i sin Ω + sin i cos Ω ⎛ ⎞ ⎜ dt dt ⎟ − sin u cos Ω − cos u cos i sin Ω ⎟ ⎜ di dΩ ⎟ = −rW ⎝ − sin u sin Ω + cos u cos i cos Ω ⎠ . μa(1 − e2 ) ⎜ ⎜ − cos i cos Ω + sin i sin Ω ⎟ dt dt ⎠ ⎝ cos u sin i di − sin i dt (17) Through some simple calculation with the above equation we can obtain
From Eq.(9) we have
r sin (ω + f ) dΩ W, = dt μa(1 − e2 ) sin i
(18)
r cos (ω + f ) di = W. dt μa(1 − e2 )
(19)
da 2a2 dr df = (S + rT ). dt μ dt dt
(20)
As for the instantaneous ellipse we also have μ dr = e sin f , dt P √ μP df r = , dt r
(21)
(22)
P = a(1 − e2 ) .
(23)
2 da [Se sin f + T (1 + e cos f )] . = √ dt n 1 − e2
(24)
Thus,
In accordance with Eq.(15) we obtain
da μ de dh = − 2ae = rT . (1 − e2 ) dt dt dt n μa(1 − e2 ) Substituting Eq.(24) into Eq.(25), we get √ 1 − e2 de = [S sin f + T (cos f + cos E)] . dt na
(25)
(26)
FENG He-sheng et al. / Chinese Astronomy and Astrophysics 31 (2007) 205–210
209
dω can be directly and quickly carried out by means of vector differentiaThe derivation of dt tion with the help of the above results. For this, a unit vector P 0 is defined, which is located on the line of apsides of the instantaneous ellipse and pointing to the apsis (see Fig.1). The instantaneous ellipse equation at any given time is r= Its vector form is
h2 . μ(1 + e cos f )
h2 = r · (R0 + eP 0 ) = r + reR0 · P 0 . μ
(27)
(28)
Differentiating the above equation, we have de dR0 dP 0 2h dh dr dr = + eP 0 · R0 + r R0 · P 0 + re · P 0 + reR0 · , μ dt dt dt dt dt dt
(29)
and obtain the derivative of the unit vector, perpendicular to the unit vector and with modulus equal to the angular rate of change of the unit vector. It is obvious that the df angular rate of change of R0 is , so the angular rate of change of P 0 is the rotation, dt dω , of the line of apsides itself plus the projection of the precessional motion of the line of dt dΩ cos i. Then one may directly intersection between the orbital planes on the orbital plane, dt write out dR0 df · P 0 = − sin f , (30) dt dt
dω dΩ dP 0 = sin f + cos i , (31) R0 · dt dt dt 2h dh dr dr de df dω dΩ = + e cos f + r cos f − re sin f + re sin f ( + cos i) . μ dt dt dt dt dt dt dt
(32)
For the elliptical motion one always has dr dr df + e cos f − re sin f = 0, dt dt dt
dω dΩ 2h de rT = r cos f + re sin f + cos i . μ dt dt dt Substituting into the relevant results, we have √ dΩ r 1 − e2 dω −S cos f + T (1 + ) sin f − = cos i . dt nae P dt
(33) (34)
(35)
Finally, we select the mean anomaly M as the sixth orbital element. By means of the Kepler’s equation for the instantaneous orbit M = E − e sin E .
(36)
210
FENG He-sheng et al. / Chinese Astronomy and Astrophysics 31 (2007) 205–210
and the above relevant results, we directly derive its variation, dM / dt. The eccentric anomaly E depends not only on the true anomaly f , but also on the eccentricity e. Making use of the geometric relation in an ellipse e + cos f . (37) cos E = 1 + e cos f to differentiate the above equation, we obtain dE df sin f r 1 de = √ −√ , (38) dt a 1 − e2 dt 1 − e2 1 + e cos f dt df dθ where the rate of change of the true anomaly depends not only on the rotation of dt dt the radius vector r itself, but simultaneously on the rotation of the line of apsides, which is also involved in the definition of the true anomaly. It is obvious that
df dω dΩ dθ = − + cos i , (39) dt dt dt dt dθ = h is the instantaneous angular momentum. So dt μa(1 − e2 ) df dω dΩ = + cos i) . −( dt r2 dt dt Differentiating the two sides of Eq.(36), we have
while r2
dM dE de = (1 − e cos E) − sin E . dt dt dt Substituting Eqs.(26), (35), (38) and (40) into Eq.(41), we obtain dM 1 − e2 r r =n− −S(cos f − 2e ) + T (1 + ) sin f . dt nae P P
(40)
(41)
(42)
3. CONCLUSIONS We have now derived the 6 perturbation equations of motion in quite a simple way. The application of the method of vector differentiation made the derivation simple, the geometric image sharp and the mechanical significance definite. The perturbation equation of motion is well known in the orbit determination through analytical means. The purpose of our re-derivation using a more objective method is to gain a better understanding of the perturbation equation of motion and to possibly provide some helps to both teaching and research. References 1
Dubiago A.D., Orbit Determination, Leningrad-Moscow, 1949
2
Brouwer D., Clemence G.M., Methods of Celestial Mechanics, New York, Academic Press, 1961
3
Liu L., Orbital Mechanics of Artificial Earth Satellites, Beijing: Higher Education Publishing House,
4
Burns J.A., Am. J. Physics, 1976, 44: 944
5
Murray C. D., Dermott S. F., Solar System Dynamics, Cambridge: Cambridge University Press, 1999,
1992, 79
54
ELSEVIER
Chinese Astronomy and Astrophysics 31 (2007) 205–210
A Direct Method of Deriving Perturbation Equation of Motion† FENG He-sheng
LI Yu-qiang
National Astronomical Observatories / Yunnan Observatory, Chinese Academy of Sciences, Kunming 650011
Abstract A direct method of deriving perturbation equation of motion is proposed. Namely, the variations in the orbital elements of an instantaneous ellipse in the sense of two-body problem under the action of three components of the perturbing force are directly derived by means of vector differentiation, with a sharp geometric image and explicit mechanical significance. The derivation of the two orbital elements, ω˙ and M˙ , is especially more concise and objective than in the traditional method. Key words: celestial mechanics
1. INTRODUCTION The derivation of perturbation equation in orbital elements is an old classical problem. The traditional method is as follows. On the basis of analytic mechanics, the relation between the variation in the canonical constants and the perturbation function is derived, then according to the correspondence of the canonical constants to the orbital elements the perturbation equation of the orbital elements is derived. Although this kind of method is of universal significance in mechanics, the process of derivation is very lengthy and tedious, and the existence of a force function for the perturbing forces has also to be identified. The direct derivation method which takes the orbital elements as the target appeared in the subsequent works [1−3] , and perturbation was no longer regarded as a conservative force, but the derivation method was based on the analytical conception of the variation of constants of the solution of differential equation and the inference thread was not clear enough. Burns[4] and Murray et al.[5] also adopted the direct method to derive the perturbation equation, †
Supported by Natural Science Foundation of Yunnan Province Received 2005–11–08; revised version 2006–01–05 A translation of Acta Astron. Sin. Vol. 48, No. 1, pp. 54–59, 2007
c 2007 Elsevier B . V. All rights reserved. 0275-1062/07/$-see front matter doi:10.1016/j.chinastron.2007.03.003
206
FENG He-sheng et al. / Chinese Astronomy and Astrophysics 31 (2007) 205–210
but the derivation process, especially the derivation of the last two orbital elements, is not sufficiently objective. In this article, the variation in the orbital elements of the instantaneous ellipse under the action of the three components of the perturbing force is directly inferred by taking advantage of the differential operation of vectors, which makes the derivation process much more concise and transparent, making the perturbation equation more clearly understood.
2. DERIVATION PROCESS We firstly establish a coordinate system which moves with the disturbed body (as shown in Fig.1), where R0 is the radial unit vector, T 0 the transverse unit vector and W 0 the normal unit vector of the orbital plane. The components of the disturbing force Fe in the above 3 directions are S, T and W , respectively. Through the coordinate rotation it is not difficult to calculate the components of R0 , T 0 and W 0 in the equatorial coordinate system O − XY Z, and the results are: ⎛ ⎞ cos u cos Ω − sin u cos i sin Ω (1) R0 = ⎝ cos u sin Ω + sin u cos i cos Ω ⎠ , sin u sin i ⎛
⎞ − sin u cos Ω − cos u cos i sin Ω T 0 = ⎝ − sin u sin Ω + cos u cos i cos Ω ⎠ , cos u cos i ⎛ ⎞ sin i sin Ω W 0 = ⎝ − sin i cos Ω ⎠ , cos i where u = ω + f is the argument of latitude.
Fig. 1
Geometry of the coordinate system moving with the disturbed body
(2)
(3)
FENG He-sheng et al. / Chinese Astronomy and Astrophysics 31 (2007) 205–210
207
For the convenience of exposition, the mass of the disturbed body is taken to be the unit mass, then the acceleration of the disturbed body is the force and the velocity is the momentum. The disturbed body is acted by the central gravitation F 0 and the perturbing force F e . We have d2 r μr = F0 + Fe = − 2 + Fe , (4) 2 dt r r r where − denotes the unit vector pointing to the center of the central body from the r disturbed body and μ is the gravitational constant. The moment of momentum of the unit mass is h, r × r˙ = h , (5) By differentiating the two sides of the above equation, we have dh μr dr dr d2 r = × + r × 2 = r × (− 3 + F e ) = r × F e . dt dt dt dt r
(6)
Thus we arrive at the conclusion: the derivative of the angular momentum of the disturbed body is equal to the moment of the perturbing force. The instantaneous ellipse at time t has dr dr 2 1 · = μ( − ) . dt dt r a
(7)
Differentiating the two sides, we have 2
dr d2 r 2 dr dr 1 da 2 1 da = μ(− 2 · + 2 ) = μ(− 3 r · + 2 ). dt dt2 r dt a dt r dt a dt
(8)
According to Eqs.(4) and (8) we have da dr 2a2 = Fe · . dt μ dt
(9)
The above equation may also be written as d −μ dr ( ) = Fe · . dt 2a dt
(10)
From the above formula we come to the conclusion that the rate of change of the total energy of the system is equal to the power of the work done by the perturbing force. In the new coordinate system (see Fig.1) we have F e = SR0 + T T 0 + W W 0 ,
(11)
r × F e = −rW T 0 + rT W 0 , h = μa(1 − e2 )W 0 , h = |h| = μa(1 − e2 ) .
(12) (13)
In the light of Eqs.(6) and (12) we have −rW T 0 + rT W 0 =
d dW 0 μa(1 − e2 )W 0 + μa(1 − e2 ) . dt dt
(14)
208
FENG He-sheng et al. / Chinese Astronomy and Astrophysics 31 (2007) 205–210
From the above equation we may obtain d μa(1 − e2 ) = rT , dt
(15)
dW 0 μa(1 − e2 ) = −rW T 0 , dt
(16)
dW 0 where and T 0 can be gotten from Eqs.(3) and (2). dt Thus, we have ⎛ di dΩ ⎞ cos i sin Ω + sin i cos Ω ⎛ ⎞ ⎜ dt dt ⎟ − sin u cos Ω − cos u cos i sin Ω ⎟ ⎜ di dΩ ⎟ = −rW ⎝ − sin u sin Ω + cos u cos i cos Ω ⎠ . μa(1 − e2 ) ⎜ ⎜ − cos i cos Ω + sin i sin Ω ⎟ dt dt ⎠ ⎝ cos u sin i di − sin i dt (17) Through some simple calculation with the above equation we can obtain
From Eq.(9) we have
r sin (ω + f ) dΩ W, = dt μa(1 − e2 ) sin i
(18)
r cos (ω + f ) di = W. dt μa(1 − e2 )
(19)
da 2a2 dr df = (S + rT ). dt μ dt dt
(20)
As for the instantaneous ellipse we also have μ dr = e sin f , dt P √ μP df r = , dt r
(21)
(22)
P = a(1 − e2 ) .
(23)
2 da [Se sin f + T (1 + e cos f )] . = √ dt n 1 − e2
(24)
Thus,
In accordance with Eq.(15) we obtain
da μ de dh = − 2ae = rT . (1 − e2 ) dt dt dt n μa(1 − e2 ) Substituting Eq.(24) into Eq.(25), we get √ 1 − e2 de = [S sin f + T (cos f + cos E)] . dt na
(25)
(26)
FENG He-sheng et al. / Chinese Astronomy and Astrophysics 31 (2007) 205–210
209
dω can be directly and quickly carried out by means of vector differentiaThe derivation of dt tion with the help of the above results. For this, a unit vector P 0 is defined, which is located on the line of apsides of the instantaneous ellipse and pointing to the apsis (see Fig.1). The instantaneous ellipse equation at any given time is r= Its vector form is
h2 . μ(1 + e cos f )
h2 = r · (R0 + eP 0 ) = r + reR0 · P 0 . μ
(27)
(28)
Differentiating the above equation, we have de dR0 dP 0 2h dh dr dr = + eP 0 · R0 + r R0 · P 0 + re · P 0 + reR0 · , μ dt dt dt dt dt dt
(29)
and obtain the derivative of the unit vector, perpendicular to the unit vector and with modulus equal to the angular rate of change of the unit vector. It is obvious that the df angular rate of change of R0 is , so the angular rate of change of P 0 is the rotation, dt dω , of the line of apsides itself plus the projection of the precessional motion of the line of dt dΩ cos i. Then one may directly intersection between the orbital planes on the orbital plane, dt write out dR0 df · P 0 = − sin f , (30) dt dt
dω dΩ dP 0 = sin f + cos i , (31) R0 · dt dt dt 2h dh dr dr de df dω dΩ = + e cos f + r cos f − re sin f + re sin f ( + cos i) . μ dt dt dt dt dt dt dt
(32)
For the elliptical motion one always has dr dr df + e cos f − re sin f = 0, dt dt dt
dω dΩ 2h de rT = r cos f + re sin f + cos i . μ dt dt dt Substituting into the relevant results, we have √ dΩ r 1 − e2 dω −S cos f + T (1 + ) sin f − = cos i . dt nae P dt
(33) (34)
(35)
Finally, we select the mean anomaly M as the sixth orbital element. By means of the Kepler’s equation for the instantaneous orbit M = E − e sin E .
(36)
210
FENG He-sheng et al. / Chinese Astronomy and Astrophysics 31 (2007) 205–210
and the above relevant results, we directly derive its variation, dM / dt. The eccentric anomaly E depends not only on the true anomaly f , but also on the eccentricity e. Making use of the geometric relation in an ellipse e + cos f . (37) cos E = 1 + e cos f to differentiate the above equation, we obtain dE df sin f r 1 de = √ −√ , (38) dt a 1 − e2 dt 1 − e2 1 + e cos f dt df dθ where the rate of change of the true anomaly depends not only on the rotation of dt dt the radius vector r itself, but simultaneously on the rotation of the line of apsides, which is also involved in the definition of the true anomaly. It is obvious that
df dω dΩ dθ = − + cos i , (39) dt dt dt dt dθ = h is the instantaneous angular momentum. So dt μa(1 − e2 ) df dω dΩ = + cos i) . −( dt r2 dt dt Differentiating the two sides of Eq.(36), we have
while r2
dM dE de = (1 − e cos E) − sin E . dt dt dt Substituting Eqs.(26), (35), (38) and (40) into Eq.(41), we obtain dM 1 − e2 r r =n− −S(cos f − 2e ) + T (1 + ) sin f . dt nae P P
(40)
(41)
(42)
3. CONCLUSIONS We have now derived the 6 perturbation equations of motion in quite a simple way. The application of the method of vector differentiation made the derivation simple, the geometric image sharp and the mechanical significance definite. The perturbation equation of motion is well known in the orbit determination through analytical means. The purpose of our re-derivation using a more objective method is to gain a better understanding of the perturbation equation of motion and to possibly provide some helps to both teaching and research. References 1
Dubiago A.D., Orbit Determination, Leningrad-Moscow, 1949
2
Brouwer D., Clemence G.M., Methods of Celestial Mechanics, New York, Academic Press, 1961
3
Liu L., Orbital Mechanics of Artificial Earth Satellites, Beijing: Higher Education Publishing House,
4
Burns J.A., Am. J. Physics, 1976, 44: 944
5
Murray C. D., Dermott S. F., Solar System Dynamics, Cambridge: Cambridge University Press, 1999,
1992, 79
54