- Email: [email protected]
RETRACTED: Analytical orbit predictions for low and high eccentricity orbits using uniformly regular KS canonical elements in an oblate atmosphere
Acta Astronautica 110 (2015) 214–224
Contents lists available at ScienceDirect
Acta Astronautica journal homepage: www.elsevier.com/locate/actaastro
Analytical orbit predictions for low and high eccentricity orbits using uniformly regular KS canonical elements in an oblate atmosphere P.S. Rajendran a,c, M. Xavier James Raj b,n a b c
Research Scholar, Manonmaniam Sundernar University, Tirunelveli, Tamil Nadu, India Orbital Analytics Section, Applied Mathematics Division, Vikram Sarabhai Space Centre, Trivandrum, India Department of Mathematics, College of Engineering, Trivandrum, India
a r t i c l e in f o
abstract
Article history: Received 24 April 2014 Received in revised form 24 December 2014 Accepted 30 January 2015 Available online 10 February 2015
Orbit prediction of satellites is an important requirement for mission planning, spacecraft navigation, re-entry, orbital lifetime estimates, etc. The prediction becomes complicated in near Earth environment as Earth's gravitational field and atmosphere influences the orbit. The effects of atmosphere are difficult to determine, since the atmospheric density undergoes large fluctuations. Method of KS elements was employed to generate analytical solutions with different atmospheric models. Particular form of KS elements known as uniformly regular KS canonical elements are used to generate number of analytical solutions with different atmospheric models, which are found to be better than KS theory. But, all are for low eccentric (e o 0.2) orbits only. In this paper, a non singular analytical theory for the motion of Earth satellite for high eccentricity orbits (e 4 0.2) in an oblate exponential atmosphere using uniformly regular KS canonical elements is developed, and integrated with the low eccentricity theory. Due to symmetry, only two of the ten equations need to be solved analytically to get complete solutions. Comparison of the important orbital parameters semi-major axis and eccentricity, which define the size and shape of orbit with KS theory, shows the superiority of the present theory. & 2015 IAA. Published by Elsevier Ltd. All rights reserved.
Keywords: Equations of motion Canonical elements Orbit prediction Oblate atmosphere
1. Introduction Accurate orbit prediction of an artificial satellite under the influence of air drag is one of the most difficult and untraceable problem in orbital dynamics. Major perturbations experienced by low Earth satellites, which are within 1000 km of the Earth surface, are due to the atmospheric drag and the in-homogeneity of the geo-potential. It is well known that the in-homogeneity of the Earth gravity field induces small periodic changes in right ascension of acceding node (Ω) and argument of perigee (ω). As such, this perturbation by itself does not influence the lifetime of a satellite [1,2]. The satellites on low Earth orbits, experiencing the influence of air drag can have perturbing effect of many times than the magnitude of other perturbing forces. Due to atmospheric resistance, there will be a continuous loss of kinetic energy from the satellite to the atmosphere, causing continuous altitude drop. The selection of an appropriate atmospheric model has
n
Corresponding author. Tel.: þ91 4712564486. E-mail address: [email protected] (M.X. James Raj).
http://dx.doi.org/10.1016/j.actaastro.2015.01.021 0094-5765/& 2015 IAA. Published by Elsevier Ltd. All rights reserved.
P.S. Rajendran, M.X. James Raj / Acta Astronautica 110 (2015) 214–224
215
great significance in the study of the perturbations of artificial Earth satellite motion caused by air resistance. Even with a good atmospheric model, describing variations with time, seasons, latitude and altitude, a complete specification of orbital decay is not possible, because of the uncertainties in the predictions of satellite's altitude and geo-magnetic indices. To predict the orbit precisely, a mathematical model representing the forces must be selected properly for integrating in the resultant differential equations of motion. The options for the mathematical solutions can be classified as analytical, semianalytical and numerical. For a first estimation of the orbit of a satellite and for a correlation between the semi-major axis and eccentricity of the orbit, which undergoes a contraction due to the perturbing effect of atmospheric drag, an analytical theory is adequate. The classical theory was documented in a monograph prepared by King-Hele [3]. Cook et al. [4] developed an analytical theory for the contraction of satellite orbits with air drag using oblate atmospheric model with constant density scale height. Cook and King-Hele [5] developed an analytical theory for an atmosphere in which the density scale height varied linearly with the distance from the Earth centre. But, all these theories were developed for low eccentricity orbits. Most of the satellites orbits are elliptical. To transfer one orbit to another, a satellite must go into an elliptical transfer orbit. Geostationary transfer orbits are intermediate orbits used for transferring the satellites into geostationary orbits. They are highly elliptical orbits with an apogee altitude of about 35,700 km and perigee altitude of a few hundred kilometers, which passes through high atmospheric region near perigee. These orbits have large variations in velocity due to perturbations during a revolution. During the planning of such missions, the analytical solutions are preferred as they represent a manifold of solutions for a large domain of initial conditions and find indispensable application to mission planning and qualitative analysis. An analytical theory for high eccentricity orbits in a spherical atmosphere was developed by King-Hele [6]. The classical Newtonian equations of motion, which is non linear is not suitable for long-term integration. Many transformations have emerged in the literature to stabilize the equations of motion, either to reduce the accumulation of local numerical errors or allowing the use of large integration step sizes, or both. One such transformation is known as KS transformation by Kustaanheimo and Stiefel [7], who regularized the nonlinear Kepler equation of motion and reduced it into linear differential equations of a harmonic oscillator of constant frequency. The method of KS total energy element equations [8] has been found to be a very powerful method for obtaining numerical solution with respect to any type of perturbing forces, as the equations are less sensitive to round off and truncation errors [9]. The equations are everywhere regular comparing to the classical Newtonian equations, which are singular at the collision of two bodies. In literatures many papers are available for orbit predictions using KS theory with oblateness and atmospheric drag as perturbing forces. The uniformly regular KS canonical (KSC) equations are a particular canonical form of the KS differential equations where all the ten KS Canonical elements αi and βi are constants for unperturbed motion. These equations permit the uniform formulation of the basic laws of elliptic, parabolic and hyperbolic motion [8]. The equations were found to provide accurate long term orbit predictions numerically with Earth's oblateness [10]. Xavier James Raj and Sharma [11] utilized these equations to generate an analytical solution for short term orbit predictions with Earth's zonal harmonic terms J2–J4. Further, these equations were utilized to include the canonical forces and analytical theories with air drag were developed for the low eccentricity orbits (eo0.2) with different atmospheric models [12]. But all these theories are for low eccentricity orbits. For high eccentricity orbits, Sharma and Xavier James Raj [13,14] utilized KS theory for analytical solutions using oblate atmosphere with variation in density scale height with altitude. In this paper a non-singular analytical solution is developed for the motion of high eccentricity satellite orbits under the influence of air drag in terms of the uniformly regular KS canonical elements using oblate atmospheric model. The analytical solutions are generated up to fourth-order terms using the independent variable λ introduced by Sterne [2]. The series expansions include upto fourth-order terms in λ and c (a small parameter dependent on the flattening of the atmosphere). The theory is developed on the assumption that the density is constant on the surfaces of spheroids of fixed ellipticity ε (equal to the Earth's ellipticity, 0.00335) whose axes coincide with the Earth's axis. The developed theory is integrated with the low eccentricity theory [12] and generated a theory for both low and high eccentricity orbits using uniformly regular KS canonical elements. Numerical experimentation with the present analytical solution for all eccentricities (0oeo1) has been carried out upto 1000 revolutions. Comparisons are made with numerically integrated values as well as KS theory. The numerical results obtained from the present analytical solutions match quite well with the numerically integrated values and found to be better than the existing KS theory, which shows the superiority of the present solution.
2. Equations of motion If Ak and Bk are the canonical forces attached to the system of equations with the uniformly regular KS canonical variables uk and wk (k ¼0–4), then the canonical equations of motion become [8, p. 250] dβ k ∂H ¼ Bk ; ∂αk ds with Hamltonian H¼
rV μ ; 4 4
d αk ∂H ¼ þ Ak ; ds ∂βk
ð1Þ
216
P.S. Rajendran, M.X. James Raj / Acta Astronautica 110 (2015) 214–224
where s is the ficitus time, V is perturbation due to oblatness, and Ak ¼
4 X
Wj
j¼0
W0 ¼
∂uj ; ∂β k
Bk ¼
4 X
Wj
j¼0
! Þ;
∂uj ; ∂ αk
! ! W j ¼ 2r LT ð u Þ: D ;
! 2r ð v : D
μ is the gravitational constant
ðk ¼ 0; 1; 2; 3; 4Þ ð2Þ ðj ¼ 1; 2; 3; 4Þ
! ! r ¼ j u j2 ¼ u21 þ u22 þ u23 þu24 ; j w j2 ¼ w21 þ w22 þ w23 þ w24 ; 2 3 u1 u2 u3 u4 6 7 u4 u3 7 6 u2 u1 ! 7: Lð u Þ ¼ 6 6 u3 u4 u1 u2 7 4 5 u4 u3 u2 u1 ! ! If D is the aerodynamic drag force per unit mass acting on a satellite of mass m with velocity vector v [1, p.29] then ! !! ð3Þ D ¼ 12 ρδ v v ; where δ ¼ FAC D =m; F ¼ ½1 r p0 Λ cos i0 =vp0 2 ; ρ is the atmospheric density, Λ is the rotational rate of the atmosphere about Earth's axis, r p0 is the initial perigee radius, i0 is the initial inclination, vp0 is the velocity at the initial perigee, CD the drag coefficient, and A and m are respectively, effective area and mass of the satellite. Since our aim is to develop an analytical solution only with air drag force, we consider only the canonical forces in the equations of motion (1). Then the equations of motion become dβ k ¼ Bk ; ds Once
d αk ¼ Ak ; ds
ðk ¼ 0; 1; 2; 3; 4Þ:
ð4Þ
αk and βk are known, they will be converted into the position and velocity vectors using the relation -x ¼ L u u;
-
where
2 -x ¼ L u w; r
_
2 4 6 X
pffiffiffiffiffiffi 1 ffiffiffiffi p α0 α2k β2k sin 2 α0 s þ
3
7 7 pffiffiffiffiffiffi 5; 2 s 2 k¼1 α0 f cos 2 α0 s 1g þ α0 βk þ α0 αk pffiffiffiffiffiffi pffiffiffiffiffiffi βk uk ¼ pffiffiffiffi α0 s αk cos α0 s ; α0 sin
u0 ¼ β0 þ 14
2
α
6 4 αk βk
3 0
w 0 ¼ α0 ;
pffiffiffiffiffiffi pffiffiffiffiffiffi pffiffiffiffiffiffi α0 s þ αk α0 sin α0 s ; ! Substituting the expression of D from Eq. (3) in Eq. (2) and substituting the resultant expression in the equations of motion (4) then [13] d α0 1 ! ¼ pffiffiffiffiffiffi ρδ v μ½1 þ e cos E; 8 α0 dE d αi 1 ! ¼ pffiffiffiffiffiffiρδ v f 0i þ f 1i cos E þ f 2i cos 2E þ f 3i sin E þ f 4i sin 2E þ f 5i E þf 6i E cos E ; ð5Þ dE 8 α0 wk ¼ βk cos
with αi þ 4 ¼ βi for i¼1–4. where f 0j ¼ mP 2j þ aQ 2j þ eðmP 0j aQ oj Þ=2; f 1j ¼ mP 0j þ aQ 0j þ eðmP 2j aQ 2j Þ; f 2j ¼ eðmP 0j aQ oj Þ=2; f 3j ¼ mP 1j þ aQ 1j; f 4j ¼ eðmP 1j aQ 1j Þ=2; f 5j ¼ mP 3j; f 6j ¼ emP 3j : For j¼ 1–4, For Aj Bj
Poj
αj/(4α0) βj/(4α0)
P1j
P2j 1.5
βj/(4α0 ) 0.5 0 )
αj/(4α
P3j
Q0j
Q1j
1.5 0 )
αj
0.5 0 )
βj/(α0 )
βj
αj (α0 )
αj/(4α0)
βj/(4α
βj/(4α0)
αj/(4α
Q2j 0.5
αj
0.5
βj
P.S. Rajendran, M.X. James Raj / Acta Astronautica 110 (2015) 214–224
217
3. Oblate exponential atmosphere In an atmosphere of constant scale height H and ellipticity ε (taken as equal to the Earth's ellipticity, 0.00335), the density ρ may be written as r σ ; ð6Þ ρ ¼ ρP exp H where
h
n
o
i
σ ¼ r P 1 ε sin 2 i sin 2 ω þ θ sin 2 ω þO ε2 ;
ð7Þ
is the radius vector of the surface of constant density that passes through the perigee point and ρP is the density at perigee, θ is the true anomaly of the satellite. In view of (6), we get
2 2 ε rp r σ r r p εr p sin i f sin 2 ðω þ θÞ sin 2 ωg þO ¼ þ : H H H H Now introduce a parameter z ¼ ae=H; so that ðr r P Þ=H ¼ zð1 cos EÞ: Writing c ¼ ðεr p sin 2 iÞ=2H; we get n o ρ ¼ ρP exp½ zð1 cos EÞ 2c sin 2 ω þ θ sin 2 ω þOðcεÞ; ¼ ρP exp½ zð1 cos EÞ þ c cos 2 ω þ θ c cos 2ω þ OðcεÞ:
ð8Þ
4. Analytical integration for low eccentricity For low eccentricity orbits, expanding expfc cos 2ðω þ θÞg up to fourth power in c, and converting θ in terms of eccentric anomaly E, substituting in (8) and retaining terms up to fourth order in e and c, we get ρ ¼ k expf βða x cos EÞg 1 þ c cos 2ðω þ EÞ 2ce sin 2ðω þEÞ sin E þ 14 c2 f1 þ cos 4ðω þ EÞg 12 ce2 f cos 2ω þ 2 cos 2ðω þ EÞ 3 cos 2ðω þ 2EÞ cos 2Eg c2 e sin 4ðω þ EÞ sin E ce3 f cos ð2ω þ3EÞ
c2 e2 c3 f3 cos ð4ω þ 2EÞ 8 cos ð4ω þ4EÞ þ 5 cos ð4ω þ6EÞg þ f3 cos ð2ω þ 2EÞ 8 24 c3 e c4 f cos ð2ω þ EÞ cos ð2ω þ 3EÞ þ cos ð6ω þ 5EÞ cos ð6ω þ 7EÞg þ f3 þ4 cos ð4ω þ 4EÞ þ cos ð6ω þ 6EÞg 8 192 4 2 3 3 2 4 4 5 þ cos ð8ω þ 8EÞg þ Oðce ; c e ; c e ; c e; ce ; c Þ; cos ð2ω þ 5EÞg=2 þ
ð9Þ
where k ¼ ρp0 expðβ r p0 c cos 2ω0 Þ and β ¼ 1=H. The explicit form of the perturbations, Δαk (k¼0, 1, 2, …., 9) is obtained by substituting the value of ρ obtained from Eq. (9) into the equations of motion (5) and noting that the non-zero integrals in the resulting expression are of the form of linear combinations of the modified Bessel functions Z 2π 1 I n ðzÞ ¼ expðz cos EÞ cos nE dE: ð10Þ 2π 0 After some algebra the resulting equations are integrated from 0 to 2π using Eq. (10). We obtain the changes in Δα0 and after one revolution as " # 3 6 X X Δα0 ¼ μT 1 Bmn cos 2nωI m ; "
Δαi ¼ T 1
m¼0 n¼0 4 9 X X
#
Δαi
ð11Þ
ðGmni cos 2nω þ H mni sin 2nωÞI n ; with ai þ 4 ¼ bi for i ¼ 1; 2; 3; 4:
m¼0 n¼0
qffiffiffi μ where T 1 ¼ 4p1ffiffiffiffi α0 ρp δμ a and the coefficients Bmn, Gmni and Hmni are provided in Appendix-A. 5. Analytical integration for high eccentricity For high eccentricity orbits, introduce a new independent variable λ defined in [2] as ð12Þ cos E ¼ 1 λ =z; where z ¼ ae=H , pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Expanding v ¼ ðμð1 þ e cos EÞÞ=ðað1 e cos EÞÞ appearing in Eq. (5) up to fourth-order terms in λ with the help of Eq. (12) 2
218
P.S. Rajendran, M.X. James Raj / Acta Astronautica 110 (2015) 214–224
we get sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4 μð1 þ eÞ X ! 2i V iλ ; v ¼ að1 eÞ
ð13Þ
i¼0
where V 0 ¼ 1;
V 1 ¼ e=½2s2 ð1 þeÞ2 ;
V 2 ¼ ð2e þ 1Þe2 =½8s4 ð1 þ eÞ4 ; V 3 ¼ ð2e2 þ eþ 1Þe3 =½16s6 ð1 þ eÞ6 ; V 4 ¼ 3 þ12eþ 12e2 þ 8e3 e4 =½128s8 ð1 þeÞ8 ; with s2 ¼ zð1 eÞ=f2ð1 þeÞg Substituting Eqs. (12) and (13) into the equations of motion (5) and expanding the resultant integrals in powers of Kð ¼ λ=sÞ, we get qffiffiffiffiffi X 4 d α0 1 ¼ pffiffiffiffiffiffi ρδμ 2zaμ W j K 2j ; 8 α0 dλ j¼0 dαi dλ
ð14Þ
qffiffiffiffiffi X 8 2μ ¼ 8p1ffiffiffiffi Sj K j ; with αi þ 4 ¼ β i for i ¼ 1; 2; 3; 4: α0 ρδ za j¼o
the coefficients Wi and Si are provided in Appendix B. Substituting Eq. (12) in Eq. (9) and expanding up to 8th power in K we get
ρ ¼ ρp expð λ2 Þ½1 þ R1 K þ R2 K 2 þR3 K 3 þ R4 K 4 þ R54 K 5 þR6 K 6 þ R7 K 7 þ R8 K 8 ;
ð15Þ
where R1 ¼ 2c sin 2ω; R2 ¼ 2c½c sin 2 2ω cos 2ω; h i 2 sin 2ω R3 ¼ c12ð1 þ eÞ 16cð1 þeÞfc sin 2ω 3 cos 2ωg 3ð5 þ 7eÞ ; " # 4cð1 þ eÞfc sin 2 2ωðc sin 2 2ω 6 cos 2ωÞ þ 3 cos 2 2ωg c ; R4 ¼ 6ð1 þeÞ 3c sin 2 2ωð7eþ 5Þ þ cos 2ωð9e þ 3Þ 2 R5 ¼
3
7 c sin 2ω 6 6 96c cos 2ωfð7 þ 20e þ13e2 Þ þ 8cð1 þeÞ2 cos 2ωg 7; 5 192ð1 þ eÞ2 4 3ð7 þ74eþ 79e2 Þ 2
R6 ¼
f512c3 cos 2ωð1 þeÞ2 þ 96c2 ð7e2 þ 12eþ 5Þg sin 3 2ω
c sin 2 2ωf2c2 sin 2 2ωð7e2 þ 12e þ 5Þ
3
6 7 6 24c2 cos 2 2ωð1 þeÞ2 12c cos 2ωð5e2 þ 8e þ 3Þ 7 6 7; 6 6ð1 þeÞ2 4 3ð8e2 þ9e þ2Þg þ cos 2ωf8c2 cos 2 2ωð1 þ e2 Þ 7 5 þ 6c cos 2ωð3e2 þ4e þ1Þ þ 3eð1 þ 2eÞg c
2
16c2 sin 3 2ωf32c2 cos 2ωð27e3 þ71e2 þ61eþ 17Þ
3
6 7 6 þ 3ð177e3 þ 391e2 þ 271e þ57Þg 7 7 c sin 2ω 6 6 2 2 3 27 R7 ¼ ω f256c ð1 þ eÞ cos 2 ω þ 48c cos 2 ω ð19e þ47e 16c cos 2 6 7; 36 7 1536ð1 þeÞ 6 7 þ 37e þ9Þ 4 5 þ 3ð227e3 þ 421e2 þ 221e þ27Þg 3ð407e3 þ431e2 þ 61e 3Þ 2
c sin 2 2ωfc2 sin 2 2ωð113e3 þ255e2 þ 183eþ 41Þ
3
6 7 6 24c2 cos 2 2ωð13e3 þ 33e2 þ 27e þ7Þ 7 6 7 6 7 3 2 3 6 6c cos 2ωð61e þ 123e þ 75e þ 13Þ 3ð30e 7 6 7 c 2 6 þ42e þ15eþ 1Þg 7: R8 ¼ 7 36 24ð1 þ eÞ 6 7 6 þ cos 2ωf16c3 ð1 þeÞ3 cos 3 2ω þ24 cos 2 2ωð3e3 7 6 7 2 6 þ7e þ 5e þ 1Þ 7 4 5 þ3c cos 2ωð25e3 þ 39e2 þ15eþ 1Þ þ 3e2 ð5eþ 3Þg
P.S. Rajendran, M.X. James Raj / Acta Astronautica 110 (2015) 214–224
By integrating the equations of motion (14), the changes in
2n ¼4 R pffiffiffiffi nX λ Δα0 ¼ μΔ1 0 2z ρ Wn dλ; s n¼0 ¼ 8 n R pffiffiffiffi nX λ Δαi ¼ Δ1 0 2z ρ Sn dλ s n¼0
219
α0 and αi during one revolution is,
ð16Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffi α0 ÞÞρp δ ð2μ=zaÞ: pffiffiffiffiffi with Δ1 ¼ pðð1Þ=ð8 ffiffiffiffiffi Since 2z 46 for e40.2 the limit 2z in Eq. (16) may be replaced by 1, with relative error less than 10 16 [2] Also if we write Z 1 1 λ2n expð λ2 Þdλ ¼ Gn ; then Gn ¼ ð2n 1ÞGn 1 for n Z 1: 2 0 pffiffiffiffi Since G0 ¼ ð1=2Þ π ; we have pffiffiffiffi pffiffiffiffi pffiffiffiffi pffiffiffiffi π and G4 ¼ 105 π: G1 ¼ 14 π ; G2 ¼ 38 π ; G3 ¼ 15 16 32 Employing the expression for ρ in Eq. (15) into the equations of motion (16) and simplifying, using the above integrals, we get hpffiffiffi i 1 Δα0 ¼ Δ1 μ 16π 8T 0 þ 4T 2 þ 6T 4 þ 15T 6 þ 105 2 T 8 þ 2ðT 1 þ T 3 þ2T 5 þ6T 7 Þ ; 2 pffiffiffi 3 π 105 ð17Þ 16 8Z k0 þ 4Z k2 þ 6Z k4 þ 15Z k6 þ 2 Z k8 5 4 ; Δαk ¼ Δ1 þ 12ðZ k1 þ Z k3 þ2Z k5 þ6Z k7 Þ with αk þ 4 ¼ βk for k ¼1–4. where Z ki ¼
i X Sj Ri j j¼0
T 2i ¼
si
;
for i ¼ 0; 1; 2; …:; 7
i X W j R2ði jÞ j¼0
s2i
and T 2i þ 1 ¼
i X W j R2ði jÞ þ 1 j¼0
s2i þ 1
for i ¼ 0; 1; 2; 3; 4
with R0 ¼1. 6. Numerical results The four equations (11) and (17) are programmed in double precision arithmetic, to compute the variations in αi and βi (i¼0–4) at the end of each revolution. Once αi and βi are known, they are converted into the position and velocity vectors, which are further converted into the osculating orbital elements. The KS elements equations of motions are numerically integrated with the analytical density model by the fourth order-Runge–Kutta–Gill method with a small step size of half degree in the eccentric anomaly. Due to air drag the maximum changes occur in the important orbital parameters semimajor axis (a) and eccentricity (e). Hence, we study the changes in these parameters. The values of the gravitational constant, Earth's equatorial radius, rotational rate of the atmosphere and the ballistic coefficient (bn ¼m/cDA) are chosen as 398600.8 km3/s2, 6378.135 km, 1.2 and 50 kg/m2, respectively. Table 1 Decrease in semi-major axis after 100 revolutions for orbits having Hp ¼ 200 km and i ¼15o. e
0.05 0.1 0.15 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
NUM (km)
20.053 16.656 16.267 16.967 20.380 26.969 38.678 60.980 110.193 252.880 1029.191
Analytical (km)
Percentage error
KSC
KS
KSC
KS
20.066 16.664 16.274 16.973 20.380 26.969 38.678 60.981 110.195 252.887 1029.248
19.876 16.542 16.173 16.879 20.387 26.976 38.687 60.996 110.228 253.014 1030.959
0.068 0.048 0.040 0.037 0.0004 0.0008 0.0012 0.0015 0.0020 0.0029 0.0056
0.088 0.068 0.058 0.052 0.032 0.026 0.024 0.026 0.032 0.053 0.172
220
P.S. Rajendran, M.X. James Raj / Acta Astronautica 110 (2015) 214–224
Table 2 Decrease in eccentricity after 100 revolutions for orbits having Hp ¼200 km and i ¼15o. e
0.05 0.1 0.15 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
NUM 108
Analytical 108
26.312 20.137 17.687 16.406 15.148 14.758 14.719 14.873 15.144 15.492 15.893
Percentage error
KSC
KS
KSC
KS
26.329 20.147 17.694 16.412 15.148 14.758 14.720 14.873 15.144 15.492 15.894
26.043 19.996 17.583 16.321 15.153 14.761 14.723 14.877 15.149 15.500 15.921
0.065 0.047 0.041 0.037 0.0007 0.0014 0.0018 0.0022 0.0027 0.0036 0.0063
1.019 0.699 0.584 0.519 0.032 0.027 0.025 0.026 0.033 0.054 0.175
1.8
Difference in semi-major axis(km)
1.6 KSC KS
1.4 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
Initial eccentricity
Fig. 1. Difference between numerically and analytically computed values of semi-major axis (km) after 100 revolutions for orbits having Hp ¼ 200 km and i¼ 15o.
KSC KS
Difference in eccentricity x 10000
0.25 0.2 0.15 0.1 0.05 0 -0.05 0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
Initial eccentricity Fig. 2. Difference between numerically and analytically computed values of eccentricity 104 after 100 revolutions for orbits having Hp ¼200 km and i¼ 15o.
To show the effectiveness of the present solution, simulations are carried out for the orbits with variations in eccentricity (0 o e o 1). The comparison of the results for the important orbital parameter semi-major axis and eccentricity are made with the existing KS theory. Tables 1 and 2 provide the decay in semi-major axis and eccentricity,
P.S. Rajendran, M.X. James Raj / Acta Astronautica 110 (2015) 214–224
221
Difference in semi-amjor axis (km)
5
4
3
2
KSC KS
1
0
-1 10
20
30
40
50
60
70
80
90
inclination (degrees)
Fig. 3. Difference between numerically and analytically computed values of semi-major axis (km) after 1000 revolutions for the orbit with the initial perigee height of 175 km, eccentricity 0.8 with variation in inclination.
respectively for the numerically integrated values (NUM), analytical solutions of the present theory (KSC) and KS theory (KS) after 100 revolutions for the orbits with the perigee height of 200 km, inclination of 15o and eccentricity varying from 0.01 to 0.9. From Table 1, it is noticed that the present solution provides better results for all values of eccentricities with the maximum deviation of only 58 m after 100 revolutions in semi-major axis over a decay of 1029.2 km, where the KS theory provides the difference of 1.8 km. It is also observed that the percentage errors are also less for the present solution comparing to the KS theory. It is true that the errors for small eccentricity as well as high eccentricity will be more, because small eccentricity orbit is near circular orbit and its decay due to drag is uniform throughout the orbits. But in eccentricity orbit the drag attraction is more at perigee, which results in past decay in apogee. Hence, when eccentricity is more, the decay is high in apogee, but perigee will remain same. This will result in accuracy of the prediction also. In Table 2, the same pattern is noticed for eccentricity also. Figs. 1 and 2 provide the differences between the present solution (KSC) and KS solutions (KS) with numerically integrated values for semi-major axis and eccentricity, respectively after 100 revolutions for the same orbits. The continuous lines indicate the present solutions and dotted lines indicate KS solutions. The figures clearly show the superiority of the present solution over the KS theory for all values of eccentricity. To know effect of inclinations, simulations are carried out up to 1000 revolutions for the orbits with the initial perigee height of 175 km, eccentricity 0.7 and variations in inclinations. Fig. 3 provides the differences between the present analytical solutions (KSC) and KS solutions (KS) with numerically integrated values for semi-major axis, which is a measure of energy of the orbit, after 1000 revolutions with inclination varies from 11 to 901. From the figure, it is understood that the difference increases for higher inclinations and the differences are very less for the present analytical theory (KSC) comparing to KS theory (KS) for all inclinations. The same pattern is observed for all inclinations in eccentricity also. From the above two tables as well as the figures, it is clear that the present analytical solution is matching quite well with the numerical values for all eccentricities and found to be better than KS theory, which shows the superiority of the present analytical theory for low as well as high eccentricity orbits.
7. Conclusion The uniformly regular KS canonical element equations are integrated analytically by the series expansion method with air drag, using an oblate exponential atmospheric model. A fourth-order non-singular analytical theory for the contraction of satellite orbit is developed in terms of λ and c. Due to symmetry in the equations of motion, only two of the nine equations are solved analytically to compute the state vector and change in energy at the end of each revolution. Numerical comparison of the present solution with the numerically integrated values for the important orbital parameters semi-major axis and eccentricity shows the accuracy of the present solution. Comparisons with the KS theory up to 1000 revolutions shows the superiority of the present solution, which establishes that the present analytical theory is one of the best theories available in literature for analytical orbit predictions. This theory can be utilized for predicting the orbit or life time of space object in low Earth orbit.
222
P.S. Rajendran, M.X. James Raj / Acta Astronautica 110 (2015) 214–224
Appendix A 9 8 3e2 þ4c2 þ 16 e2 c4 16c2 64 aQ 2i > > > > = < 2 2 2 4 2 1 21e þ 12 c þ 4 e þ c þ 16c þ64 mP ; G00i ¼ 64 2i > > > > ; : 8 3e2 þ 2c2 þ 8 emP 0i ( 2 1 G01i ¼ 64 5e þ 6c2 þ 24 e2 c4 16c2 64 aQ 0i 16 3e3 þ 2c2 þ 8 e mP 2i
G02i
½ 35e2 þ 18c2 þ 72 e2 þ c4 þ 16c2 þ 64mP 0i " 2 # e þ c2 þ 4 e2 aQ 2i 1 ¼ 16 ; 7e2 þ3c2 þ 12 e2 mP 2i 4 2e2 þc2 þ4 emP 0i 4
e ½ðaQ 2i 7mP 2i Þe 8mP 0i ; G04i ¼ 64
4
e G05i ¼ 128 ðaQ 0i 7mP 0i Þ;
ce 2aeQ 2i þ 2e2 c2 8 Q 0i a 6emP 2i 6ce2 mP 0i ; 16
1 2 2 4e c þ 8c 2eaQ 2i 8e2 c2 8 caQ 0i 12ce3 mP 2i 8e2 c2 8 cmP 0i ¼ 16 ( ) 2 2 3 1 2 10e þ c þ8 caQ 2i 4ac Q 0i 2 ; ¼ 16 2 2e c2 8 cmP 2i þ 3ð3c 2Þe2 2 c2 þ 8 c emP 0i 2 " 2 # 2 2 c 2 5e c 8 aeQ 2i þ 5e c 8 aQ 0i 2 ; ¼ 16 4 c þ 8 emP 2i 15e2 þ c2 8 mP 0i 2 " # 2 ce 10aeQ 2i þ c þ 8 2e aQ 0i þ 2ð2eþ 17ÞemP 2i ; ¼ 16 þ2 5e2 þc2 þ8 mP 0i
G10i ¼ G11i G12i G13i G14i
ce2 ð6aeQ 2i þ 5aQ 0i þ 28emP 2i þ 17mP 0i Þ; 16 3 ce G16i ¼ ðQ 3a þ 14mP 0i Þ; 16 0i c2 e2 G21i ¼ ð11Q 0i a mP 0i Þ; 64 2 c e G22i ¼ ð11eaQ 2i 8Q 0i a emP 2i 4mP 0i Þ; 32 1 96c2 eQ 2i a þ 69e2 2c2 24c2 Q 0i a þ48c2 emP 2i þ 15e2 þ 96e 2c2 24 c2 mP 0i G23i ¼ 192 1 51e2 c2 12 c2 Q 2i a þ 63e2 24e c2 12 mP 2i 12c2 emP 0i G24i ¼ 48 1 96c2 eaQ 2i þ 45e2 2c2 þ 24 c2 Q 0i a þ ð48ðeþ 2Þc2 emP 2i þ 39e2 þ 2c2 þ24 c2 mP 0i G25i ¼ 192 c2 e ð19aeQ 2i þ8aQ 0i þ 39emP 2i þ 12emP 0i Þ; G26i ¼ 32 2 2 c e G27i ¼ ð19aQ 0i þ 39mP 0i Þ; 64 3 c e G34i ¼ ð3aQ 0i þ2mP 0i Þ; 16 3 c G35i ¼ ð6eaQ 2i aQ 0i þ 4emP 2i mP 0i Þ; 48 c3 G36i ¼ ðaQ 2i þ mP 2i þ emP 0i Þ; 24 c3 G37i ¼ ð6aeQ 2i þ aQ 0i þ 8emP 2i þ mP 0i Þ; 48 c3 e ½3aQ 0i þ 4mP 0i ; G38i ¼ 48 c4 4 G47i ¼ c ðaQ 0i þ mP 0i Þ 384 1 c4 ðaQ 0i þmP 0i Þ; G48i ¼ ðaQ 2i þ mP 0i Þ; G49i ¼ 192 384 G15i ¼
P.S. Rajendran, M.X. James Raj / Acta Astronautica 110 (2015) 214–224
H 10i ¼
ce 16
c2 þ 8 aQ 1i þ 2c2 þ 16 11e2 mP 1i
H 11i ¼
c 10e2 c2 8 aQ 1i ð14e2 þc2 þ8ÞmP 1i 16
H 12i ¼
ce 2 4e þ3c2 þ 24 m þ 2 5e2 þ2c2 þ 8 a P 1i 16
H 13i ¼
c 15e2 c2 8 aQ 1i þðce2 c2 8ÞmP 1i 16
H 14i ¼
ce 2 8e c2 8 aQ 1i þ 2 7e2 c2 8 mP 1i 16
H 15i ¼
ce2 ð7aQ 1i þ mP 1i Þ; 16
H 22i ¼
c2 e ð2aQ 1i þ 3mP 1i Þ; 4
H 25i ¼
H 27i ¼
H 16i ¼
H 23i ¼
c3 e ½aQ 1i þ mP 1i ; 48
H 47i ¼
c4 ðaQ 1i þmP 1i Þ; 384
e2 2 7e þ 3 c2 þ4 ; 16
B04 ¼
7e4 ; 64
B12 ¼
ce 2 5e c2 8 ; 4
B15 ¼
ce3 ; 8
B24 ¼
1 2 c e; 2
B36 ¼
1 3 c e; 12
B22 ¼
3 2 2 c e ; 32
B25 ¼
ce2 ð3e 19Þ; 8
B26 ¼
c2 ðaQ 1i þmP 1i Þ; 48
c3 e ð3aQ 1i þ4mP 1i Þ; 8
B00 ¼
e 2 3e þ 2 c2 þ 4 ; 4
B03 ¼
e3 ; 4
ce2 ðe 13Þ; 8
B23 ¼
35 2 2 c e ; 32
H 35i ¼
H 38i ¼
e2 2 7e þ 3 c2 þ 4 ; 16
B11 ¼
B13 ¼
3c2 e2 ðaQ 1i þ mP 1i Þ; 16
c2 e ð2aQ 1i þ 3mP 1i Þ; 8
c3 e ð3aQ 1i þ2mP 1i Þ; 8
c4 ðaQ 1i þ mP 1i Þ; 384
B02 ¼
ce2 ðeþ 3Þ; 8
H 26i ¼
c3 ½aQ 1i þ mP 1i ; 48
H 49i ¼
B01 ¼
B10 ¼
H 34i ¼
H 37i ¼
H 21i ¼
c2 36ae2 Q 1i þ ð135e2 2c2 24ÞmP 1i 16
c2 60ae2 Q 1i þ 195e2 2c2 24 mP 1i ; 16
c2 e2 ð19aQ 1i þ39mP 1i Þ; 64
H 36i ¼
ce3 ð2aQ 1i þ13mP 1i Þ; 16
B14 ¼
3 2 ce ð5eþ 1Þ; 8
29 2 2 c e ; 32
3 2 2 c e ; 32
All other coefficients are zero:
Appendix B 3 2
W 0 ¼ ð1 þ eÞ1 ; ð1 eÞ2
W 1 ¼ 3e 8e þ11; 8ð1 e2 Þ2 2
W2 ¼
ð5e4 16e3 50e2 þ 16e 3Þ
W4 ¼
ð45e8 160e7 4228e6 4128e5 418e4 þ 32e3 228e2 þ 160e 35Þ
1
128ð1 þ eÞ2 ð1 e2 Þ2
;
W 3 ¼ ð7e
6
1
32768ð1 þ eÞ6 ð1 e2 Þ2
24e5 243e4 80e3 þ 45e2 24e þ 5Þ 1
1024ð1 þ eÞ4 ð1 e2 Þ2
:
;
223
224
P.S. Rajendran, M.X. James Raj / Acta Astronautica 110 (2015) 214–224
So ¼ ðR0 þ R2 þ R4 ÞK; S1 ¼ eR1 þR3 þR5 þ R6 ; h i ð 15e2 þ 4e þ 15ÞR0 þ ð 3e2 þ 4e þ 3ÞR2 þ ðe2 þ 4e 1ÞR4
S2 ¼
8ð1 þ eÞ2
ð3e
K;
3e 3eÞR1 3eR3 þ ð e 3e þ 1ÞR5 þ ð2e 3e 2ÞR6 ; 6ð1 þ eÞ2 3 2 4 3 2 4 35e 88e 54e þ 120e þ35 R0 þ 5e þ 8e3 þ 26e2 þ 24e 5 R2 6 7 þ 3e4 þ 40e3 þ 10e2 8e þ3 R4 6 7 7K; S4 ¼ 6 4 6 7 128ð1 þ eÞ 4 5
S3 ¼
3
2
2 6 6 S5 ¼ 6 6 4
2
2
3 15e3 30e2 R1 þ 30e3 15e2 R3 þ 4e4 þ 40e3 7e2 þ 10e 4 R5 4 7 þ 6e 10e3 27e2 20eþ 6 R6 7 7; 4 7 120ð1 þ eÞ 5
2 21e6 212e5 175e4 þ 264e3 þ 303e2 þ 140e 21R 3 0 6 þ 7e6 þ12e5 þ101e4 þ 200e3 þ27e2 20e þ 7R2 7 7 6 6 7 6 7 þ 5e6 þ 172e5 þ 129e4 þ8e3 e2 þ 12e 5 R4 6 7 S6 ¼ 6 7K; 1024ð1 þ eÞ6 7 6 7 6 6 7 4 5 2 6 6 6 6 6 S7 ¼ 6 6 6 6 4
3 105e5 315e4 105e3 R1 þ 210e5 210e4 105e3 R3 7 þ 12e6 308e5 209e4 þ 21e3 e2 28e þ12 R5 7 7 7 þ ð16e6 28e5 188e4 329e3 22e2 42e 16ÞR6 7 7; 6 1680ð1 þ eÞ 7 7 7 5
2 99e8 3056e7 3884e6 þ 1424e5 þ 6802e4 þ 7088e3 þ 724e2 336eþ 99R 3 0 6 þ 45e8 þ 80e7 þ 1556e6 þ 4816e5 þ 2354e4 þ112e3 þ 20e2 þ 112e 45R2 7 6 7 6 7 6 þ 35e8 þ2832e7 þ 3540e6 þ 2192e5 þ 274e4 þ 176e3 44e2 80e þ 35 R4 7 6 7 S8 ¼ 6 7K; 32768ð1 þ eÞ8 6 7 6 7 7 6 4 5 1 2
with K ¼ ð1 eÞ1 ð1 þ eÞ2
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]
D.G. King-Hele, Satellite Orbits in an Atmosphere: Theory and Applications, Blackie, 1987. T.E. Sterne, An Introduction to Celestial Mechanics, Interscience, New York, 1960. D.G. King-Hele, Theory of Satellite Orbits in an Atmosphere, Butter Worths, London, 1964. G.E. Cook, D.G. King-Hele, D.M.C. Walker, The contraction of satellite orbits under the influence of air drag. II with oblate atmosphere, Proc. R. Soc. A 264 (1316) (1961) 88–121. G.E. Cook, D.G. King-Hele, The contraction of satellite orbits under the influence of air drag IV with scale height depending on altitude, Proc. R. Soc. 275 (1362) (1963) 357–390. D.G. King-Hele, The contraction of satellite orbits under the influence of air drag III. High-eccentricity orbits (0.2 re o1), Proc. R. Soc. A 267 (1331) (1962) 541–557. P. Kustaanheimo, E.L. Stiefel, Perturbation theory of Kepler motion based on Spinor regularization, J. Reine Angew. Math. 218 (1965) 204–219. E.L. Stiefel, G. Scheifele, Linear and Regular Celestial Mechanics, Springer, Berlin, 1971. O.F. Graf Jr., A. Mueller, S. Starke, The Method of Averages Applied to the KS Differential Equations, NASA-CR-151607, 1977. R.K. Sharma, M. Xavier James Raj, Long-term computations with ks uniformly regular canonical elements with oblateness, Earth Moon Planets 32 (1988) 63–178. M. Xavier James Raj, R.K. Sharma, Analytical short-term orbit predictions with J2, J3, J4 in terms of K–S uniformly regular canonical elements, Adv. Space Res. 31 (2003) 2019–2025. M. Xavier James Raj, R.K. Sharma, Contraction of near-Earth satellite orbits using uniformly regular KS canonical elements in an oblate atmosphere with density scale height variation with altitude, Planet. Space Sci. 57 (2009) 34–41. R.K. Sharma, M. Xavier James Raj, Analytical approach using KS elements to high eccentricity orbit prediction including drag, Acta Astronaut. 66 (2010) 864–870. R.K. Sharma, M. Xavier James Raj, Fourth-order theories for orbit predictions for low-and high-eccentricity orbits in an oblate atmosphere with scale height dependent on altitude, Acta Astronaut. 65 (2009) 1336–1344.
Contents lists available at ScienceDirect
Acta Astronautica journal homepage: www.elsevier.com/locate/actaastro
Analytical orbit predictions for low and high eccentricity orbits using uniformly regular KS canonical elements in an oblate atmosphere P.S. Rajendran a,c, M. Xavier James Raj b,n a b c
Research Scholar, Manonmaniam Sundernar University, Tirunelveli, Tamil Nadu, India Orbital Analytics Section, Applied Mathematics Division, Vikram Sarabhai Space Centre, Trivandrum, India Department of Mathematics, College of Engineering, Trivandrum, India
a r t i c l e in f o
abstract
Article history: Received 24 April 2014 Received in revised form 24 December 2014 Accepted 30 January 2015 Available online 10 February 2015
Orbit prediction of satellites is an important requirement for mission planning, spacecraft navigation, re-entry, orbital lifetime estimates, etc. The prediction becomes complicated in near Earth environment as Earth's gravitational field and atmosphere influences the orbit. The effects of atmosphere are difficult to determine, since the atmospheric density undergoes large fluctuations. Method of KS elements was employed to generate analytical solutions with different atmospheric models. Particular form of KS elements known as uniformly regular KS canonical elements are used to generate number of analytical solutions with different atmospheric models, which are found to be better than KS theory. But, all are for low eccentric (e o 0.2) orbits only. In this paper, a non singular analytical theory for the motion of Earth satellite for high eccentricity orbits (e 4 0.2) in an oblate exponential atmosphere using uniformly regular KS canonical elements is developed, and integrated with the low eccentricity theory. Due to symmetry, only two of the ten equations need to be solved analytically to get complete solutions. Comparison of the important orbital parameters semi-major axis and eccentricity, which define the size and shape of orbit with KS theory, shows the superiority of the present theory. & 2015 IAA. Published by Elsevier Ltd. All rights reserved.
Keywords: Equations of motion Canonical elements Orbit prediction Oblate atmosphere
1. Introduction Accurate orbit prediction of an artificial satellite under the influence of air drag is one of the most difficult and untraceable problem in orbital dynamics. Major perturbations experienced by low Earth satellites, which are within 1000 km of the Earth surface, are due to the atmospheric drag and the in-homogeneity of the geo-potential. It is well known that the in-homogeneity of the Earth gravity field induces small periodic changes in right ascension of acceding node (Ω) and argument of perigee (ω). As such, this perturbation by itself does not influence the lifetime of a satellite [1,2]. The satellites on low Earth orbits, experiencing the influence of air drag can have perturbing effect of many times than the magnitude of other perturbing forces. Due to atmospheric resistance, there will be a continuous loss of kinetic energy from the satellite to the atmosphere, causing continuous altitude drop. The selection of an appropriate atmospheric model has
n
Corresponding author. Tel.: þ91 4712564486. E-mail address: [email protected] (M.X. James Raj).
http://dx.doi.org/10.1016/j.actaastro.2015.01.021 0094-5765/& 2015 IAA. Published by Elsevier Ltd. All rights reserved.
P.S. Rajendran, M.X. James Raj / Acta Astronautica 110 (2015) 214–224
215
great significance in the study of the perturbations of artificial Earth satellite motion caused by air resistance. Even with a good atmospheric model, describing variations with time, seasons, latitude and altitude, a complete specification of orbital decay is not possible, because of the uncertainties in the predictions of satellite's altitude and geo-magnetic indices. To predict the orbit precisely, a mathematical model representing the forces must be selected properly for integrating in the resultant differential equations of motion. The options for the mathematical solutions can be classified as analytical, semianalytical and numerical. For a first estimation of the orbit of a satellite and for a correlation between the semi-major axis and eccentricity of the orbit, which undergoes a contraction due to the perturbing effect of atmospheric drag, an analytical theory is adequate. The classical theory was documented in a monograph prepared by King-Hele [3]. Cook et al. [4] developed an analytical theory for the contraction of satellite orbits with air drag using oblate atmospheric model with constant density scale height. Cook and King-Hele [5] developed an analytical theory for an atmosphere in which the density scale height varied linearly with the distance from the Earth centre. But, all these theories were developed for low eccentricity orbits. Most of the satellites orbits are elliptical. To transfer one orbit to another, a satellite must go into an elliptical transfer orbit. Geostationary transfer orbits are intermediate orbits used for transferring the satellites into geostationary orbits. They are highly elliptical orbits with an apogee altitude of about 35,700 km and perigee altitude of a few hundred kilometers, which passes through high atmospheric region near perigee. These orbits have large variations in velocity due to perturbations during a revolution. During the planning of such missions, the analytical solutions are preferred as they represent a manifold of solutions for a large domain of initial conditions and find indispensable application to mission planning and qualitative analysis. An analytical theory for high eccentricity orbits in a spherical atmosphere was developed by King-Hele [6]. The classical Newtonian equations of motion, which is non linear is not suitable for long-term integration. Many transformations have emerged in the literature to stabilize the equations of motion, either to reduce the accumulation of local numerical errors or allowing the use of large integration step sizes, or both. One such transformation is known as KS transformation by Kustaanheimo and Stiefel [7], who regularized the nonlinear Kepler equation of motion and reduced it into linear differential equations of a harmonic oscillator of constant frequency. The method of KS total energy element equations [8] has been found to be a very powerful method for obtaining numerical solution with respect to any type of perturbing forces, as the equations are less sensitive to round off and truncation errors [9]. The equations are everywhere regular comparing to the classical Newtonian equations, which are singular at the collision of two bodies. In literatures many papers are available for orbit predictions using KS theory with oblateness and atmospheric drag as perturbing forces. The uniformly regular KS canonical (KSC) equations are a particular canonical form of the KS differential equations where all the ten KS Canonical elements αi and βi are constants for unperturbed motion. These equations permit the uniform formulation of the basic laws of elliptic, parabolic and hyperbolic motion [8]. The equations were found to provide accurate long term orbit predictions numerically with Earth's oblateness [10]. Xavier James Raj and Sharma [11] utilized these equations to generate an analytical solution for short term orbit predictions with Earth's zonal harmonic terms J2–J4. Further, these equations were utilized to include the canonical forces and analytical theories with air drag were developed for the low eccentricity orbits (eo0.2) with different atmospheric models [12]. But all these theories are for low eccentricity orbits. For high eccentricity orbits, Sharma and Xavier James Raj [13,14] utilized KS theory for analytical solutions using oblate atmosphere with variation in density scale height with altitude. In this paper a non-singular analytical solution is developed for the motion of high eccentricity satellite orbits under the influence of air drag in terms of the uniformly regular KS canonical elements using oblate atmospheric model. The analytical solutions are generated up to fourth-order terms using the independent variable λ introduced by Sterne [2]. The series expansions include upto fourth-order terms in λ and c (a small parameter dependent on the flattening of the atmosphere). The theory is developed on the assumption that the density is constant on the surfaces of spheroids of fixed ellipticity ε (equal to the Earth's ellipticity, 0.00335) whose axes coincide with the Earth's axis. The developed theory is integrated with the low eccentricity theory [12] and generated a theory for both low and high eccentricity orbits using uniformly regular KS canonical elements. Numerical experimentation with the present analytical solution for all eccentricities (0oeo1) has been carried out upto 1000 revolutions. Comparisons are made with numerically integrated values as well as KS theory. The numerical results obtained from the present analytical solutions match quite well with the numerically integrated values and found to be better than the existing KS theory, which shows the superiority of the present solution.
2. Equations of motion If Ak and Bk are the canonical forces attached to the system of equations with the uniformly regular KS canonical variables uk and wk (k ¼0–4), then the canonical equations of motion become [8, p. 250] dβ k ∂H ¼ Bk ; ∂αk ds with Hamltonian H¼
rV μ ; 4 4
d αk ∂H ¼ þ Ak ; ds ∂βk
ð1Þ
216
P.S. Rajendran, M.X. James Raj / Acta Astronautica 110 (2015) 214–224
where s is the ficitus time, V is perturbation due to oblatness, and Ak ¼
4 X
Wj
j¼0
W0 ¼
∂uj ; ∂β k
Bk ¼
4 X
Wj
j¼0
! Þ;
∂uj ; ∂ αk
! ! W j ¼ 2r LT ð u Þ: D ;
! 2r ð v : D
μ is the gravitational constant
ðk ¼ 0; 1; 2; 3; 4Þ ð2Þ ðj ¼ 1; 2; 3; 4Þ
! ! r ¼ j u j2 ¼ u21 þ u22 þ u23 þu24 ; j w j2 ¼ w21 þ w22 þ w23 þ w24 ; 2 3 u1 u2 u3 u4 6 7 u4 u3 7 6 u2 u1 ! 7: Lð u Þ ¼ 6 6 u3 u4 u1 u2 7 4 5 u4 u3 u2 u1 ! ! If D is the aerodynamic drag force per unit mass acting on a satellite of mass m with velocity vector v [1, p.29] then ! !! ð3Þ D ¼ 12 ρδ v v ; where δ ¼ FAC D =m; F ¼ ½1 r p0 Λ cos i0 =vp0 2 ; ρ is the atmospheric density, Λ is the rotational rate of the atmosphere about Earth's axis, r p0 is the initial perigee radius, i0 is the initial inclination, vp0 is the velocity at the initial perigee, CD the drag coefficient, and A and m are respectively, effective area and mass of the satellite. Since our aim is to develop an analytical solution only with air drag force, we consider only the canonical forces in the equations of motion (1). Then the equations of motion become dβ k ¼ Bk ; ds Once
d αk ¼ Ak ; ds
ðk ¼ 0; 1; 2; 3; 4Þ:
ð4Þ
αk and βk are known, they will be converted into the position and velocity vectors using the relation -x ¼ L u u;
-
where
2 -x ¼ L u w; r
_
2 4 6 X
pffiffiffiffiffiffi 1 ffiffiffiffi p α0 α2k β2k sin 2 α0 s þ
3
7 7 pffiffiffiffiffiffi 5; 2 s 2 k¼1 α0 f cos 2 α0 s 1g þ α0 βk þ α0 αk pffiffiffiffiffiffi pffiffiffiffiffiffi βk uk ¼ pffiffiffiffi α0 s αk cos α0 s ; α0 sin
u0 ¼ β0 þ 14
2
α
6 4 αk βk
3 0
w 0 ¼ α0 ;
pffiffiffiffiffiffi pffiffiffiffiffiffi pffiffiffiffiffiffi α0 s þ αk α0 sin α0 s ; ! Substituting the expression of D from Eq. (3) in Eq. (2) and substituting the resultant expression in the equations of motion (4) then [13] d α0 1 ! ¼ pffiffiffiffiffiffi ρδ v μ½1 þ e cos E; 8 α0 dE d αi 1 ! ¼ pffiffiffiffiffiffiρδ v f 0i þ f 1i cos E þ f 2i cos 2E þ f 3i sin E þ f 4i sin 2E þ f 5i E þf 6i E cos E ; ð5Þ dE 8 α0 wk ¼ βk cos
with αi þ 4 ¼ βi for i¼1–4. where f 0j ¼ mP 2j þ aQ 2j þ eðmP 0j aQ oj Þ=2; f 1j ¼ mP 0j þ aQ 0j þ eðmP 2j aQ 2j Þ; f 2j ¼ eðmP 0j aQ oj Þ=2; f 3j ¼ mP 1j þ aQ 1j; f 4j ¼ eðmP 1j aQ 1j Þ=2; f 5j ¼ mP 3j; f 6j ¼ emP 3j : For j¼ 1–4, For Aj Bj
Poj
αj/(4α0) βj/(4α0)
P1j
P2j 1.5
βj/(4α0 ) 0.5 0 )
αj/(4α
P3j
Q0j
Q1j
1.5 0 )
αj
0.5 0 )
βj/(α0 )
βj
αj (α0 )
αj/(4α0)
βj/(4α
βj/(4α0)
αj/(4α
Q2j 0.5
αj
0.5
βj
P.S. Rajendran, M.X. James Raj / Acta Astronautica 110 (2015) 214–224
217
3. Oblate exponential atmosphere In an atmosphere of constant scale height H and ellipticity ε (taken as equal to the Earth's ellipticity, 0.00335), the density ρ may be written as r σ ; ð6Þ ρ ¼ ρP exp H where
h
n
o
i
σ ¼ r P 1 ε sin 2 i sin 2 ω þ θ sin 2 ω þO ε2 ;
ð7Þ
is the radius vector of the surface of constant density that passes through the perigee point and ρP is the density at perigee, θ is the true anomaly of the satellite. In view of (6), we get
2 2 ε rp r σ r r p εr p sin i f sin 2 ðω þ θÞ sin 2 ωg þO ¼ þ : H H H H Now introduce a parameter z ¼ ae=H; so that ðr r P Þ=H ¼ zð1 cos EÞ: Writing c ¼ ðεr p sin 2 iÞ=2H; we get n o ρ ¼ ρP exp½ zð1 cos EÞ 2c sin 2 ω þ θ sin 2 ω þOðcεÞ; ¼ ρP exp½ zð1 cos EÞ þ c cos 2 ω þ θ c cos 2ω þ OðcεÞ:
ð8Þ
4. Analytical integration for low eccentricity For low eccentricity orbits, expanding expfc cos 2ðω þ θÞg up to fourth power in c, and converting θ in terms of eccentric anomaly E, substituting in (8) and retaining terms up to fourth order in e and c, we get ρ ¼ k expf βða x cos EÞg 1 þ c cos 2ðω þ EÞ 2ce sin 2ðω þEÞ sin E þ 14 c2 f1 þ cos 4ðω þ EÞg 12 ce2 f cos 2ω þ 2 cos 2ðω þ EÞ 3 cos 2ðω þ 2EÞ cos 2Eg c2 e sin 4ðω þ EÞ sin E ce3 f cos ð2ω þ3EÞ
c2 e2 c3 f3 cos ð4ω þ 2EÞ 8 cos ð4ω þ4EÞ þ 5 cos ð4ω þ6EÞg þ f3 cos ð2ω þ 2EÞ 8 24 c3 e c4 f cos ð2ω þ EÞ cos ð2ω þ 3EÞ þ cos ð6ω þ 5EÞ cos ð6ω þ 7EÞg þ f3 þ4 cos ð4ω þ 4EÞ þ cos ð6ω þ 6EÞg 8 192 4 2 3 3 2 4 4 5 þ cos ð8ω þ 8EÞg þ Oðce ; c e ; c e ; c e; ce ; c Þ; cos ð2ω þ 5EÞg=2 þ
ð9Þ
where k ¼ ρp0 expðβ r p0 c cos 2ω0 Þ and β ¼ 1=H. The explicit form of the perturbations, Δαk (k¼0, 1, 2, …., 9) is obtained by substituting the value of ρ obtained from Eq. (9) into the equations of motion (5) and noting that the non-zero integrals in the resulting expression are of the form of linear combinations of the modified Bessel functions Z 2π 1 I n ðzÞ ¼ expðz cos EÞ cos nE dE: ð10Þ 2π 0 After some algebra the resulting equations are integrated from 0 to 2π using Eq. (10). We obtain the changes in Δα0 and after one revolution as " # 3 6 X X Δα0 ¼ μT 1 Bmn cos 2nωI m ; "
Δαi ¼ T 1
m¼0 n¼0 4 9 X X
#
Δαi
ð11Þ
ðGmni cos 2nω þ H mni sin 2nωÞI n ; with ai þ 4 ¼ bi for i ¼ 1; 2; 3; 4:
m¼0 n¼0
qffiffiffi μ where T 1 ¼ 4p1ffiffiffiffi α0 ρp δμ a and the coefficients Bmn, Gmni and Hmni are provided in Appendix-A. 5. Analytical integration for high eccentricity For high eccentricity orbits, introduce a new independent variable λ defined in [2] as ð12Þ cos E ¼ 1 λ =z; where z ¼ ae=H , pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Expanding v ¼ ðμð1 þ e cos EÞÞ=ðað1 e cos EÞÞ appearing in Eq. (5) up to fourth-order terms in λ with the help of Eq. (12) 2
218
P.S. Rajendran, M.X. James Raj / Acta Astronautica 110 (2015) 214–224
we get sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4 μð1 þ eÞ X ! 2i V iλ ; v ¼ að1 eÞ
ð13Þ
i¼0
where V 0 ¼ 1;
V 1 ¼ e=½2s2 ð1 þeÞ2 ;
V 2 ¼ ð2e þ 1Þe2 =½8s4 ð1 þ eÞ4 ; V 3 ¼ ð2e2 þ eþ 1Þe3 =½16s6 ð1 þ eÞ6 ; V 4 ¼ 3 þ12eþ 12e2 þ 8e3 e4 =½128s8 ð1 þeÞ8 ; with s2 ¼ zð1 eÞ=f2ð1 þeÞg Substituting Eqs. (12) and (13) into the equations of motion (5) and expanding the resultant integrals in powers of Kð ¼ λ=sÞ, we get qffiffiffiffiffi X 4 d α0 1 ¼ pffiffiffiffiffiffi ρδμ 2zaμ W j K 2j ; 8 α0 dλ j¼0 dαi dλ
ð14Þ
qffiffiffiffiffi X 8 2μ ¼ 8p1ffiffiffiffi Sj K j ; with αi þ 4 ¼ β i for i ¼ 1; 2; 3; 4: α0 ρδ za j¼o
the coefficients Wi and Si are provided in Appendix B. Substituting Eq. (12) in Eq. (9) and expanding up to 8th power in K we get
ρ ¼ ρp expð λ2 Þ½1 þ R1 K þ R2 K 2 þR3 K 3 þ R4 K 4 þ R54 K 5 þR6 K 6 þ R7 K 7 þ R8 K 8 ;
ð15Þ
where R1 ¼ 2c sin 2ω; R2 ¼ 2c½c sin 2 2ω cos 2ω; h i 2 sin 2ω R3 ¼ c12ð1 þ eÞ 16cð1 þeÞfc sin 2ω 3 cos 2ωg 3ð5 þ 7eÞ ; " # 4cð1 þ eÞfc sin 2 2ωðc sin 2 2ω 6 cos 2ωÞ þ 3 cos 2 2ωg c ; R4 ¼ 6ð1 þeÞ 3c sin 2 2ωð7eþ 5Þ þ cos 2ωð9e þ 3Þ 2 R5 ¼
3
7 c sin 2ω 6 6 96c cos 2ωfð7 þ 20e þ13e2 Þ þ 8cð1 þeÞ2 cos 2ωg 7; 5 192ð1 þ eÞ2 4 3ð7 þ74eþ 79e2 Þ 2
R6 ¼
f512c3 cos 2ωð1 þeÞ2 þ 96c2 ð7e2 þ 12eþ 5Þg sin 3 2ω
c sin 2 2ωf2c2 sin 2 2ωð7e2 þ 12e þ 5Þ
3
6 7 6 24c2 cos 2 2ωð1 þeÞ2 12c cos 2ωð5e2 þ 8e þ 3Þ 7 6 7; 6 6ð1 þeÞ2 4 3ð8e2 þ9e þ2Þg þ cos 2ωf8c2 cos 2 2ωð1 þ e2 Þ 7 5 þ 6c cos 2ωð3e2 þ4e þ1Þ þ 3eð1 þ 2eÞg c
2
16c2 sin 3 2ωf32c2 cos 2ωð27e3 þ71e2 þ61eþ 17Þ
3
6 7 6 þ 3ð177e3 þ 391e2 þ 271e þ57Þg 7 7 c sin 2ω 6 6 2 2 3 27 R7 ¼ ω f256c ð1 þ eÞ cos 2 ω þ 48c cos 2 ω ð19e þ47e 16c cos 2 6 7; 36 7 1536ð1 þeÞ 6 7 þ 37e þ9Þ 4 5 þ 3ð227e3 þ 421e2 þ 221e þ27Þg 3ð407e3 þ431e2 þ 61e 3Þ 2
c sin 2 2ωfc2 sin 2 2ωð113e3 þ255e2 þ 183eþ 41Þ
3
6 7 6 24c2 cos 2 2ωð13e3 þ 33e2 þ 27e þ7Þ 7 6 7 6 7 3 2 3 6 6c cos 2ωð61e þ 123e þ 75e þ 13Þ 3ð30e 7 6 7 c 2 6 þ42e þ15eþ 1Þg 7: R8 ¼ 7 36 24ð1 þ eÞ 6 7 6 þ cos 2ωf16c3 ð1 þeÞ3 cos 3 2ω þ24 cos 2 2ωð3e3 7 6 7 2 6 þ7e þ 5e þ 1Þ 7 4 5 þ3c cos 2ωð25e3 þ 39e2 þ15eþ 1Þ þ 3e2 ð5eþ 3Þg
P.S. Rajendran, M.X. James Raj / Acta Astronautica 110 (2015) 214–224
By integrating the equations of motion (14), the changes in
2n ¼4 R pffiffiffiffi nX λ Δα0 ¼ μΔ1 0 2z ρ Wn dλ; s n¼0 ¼ 8 n R pffiffiffiffi nX λ Δαi ¼ Δ1 0 2z ρ Sn dλ s n¼0
219
α0 and αi during one revolution is,
ð16Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffi α0 ÞÞρp δ ð2μ=zaÞ: pffiffiffiffiffi with Δ1 ¼ pðð1Þ=ð8 ffiffiffiffiffi Since 2z 46 for e40.2 the limit 2z in Eq. (16) may be replaced by 1, with relative error less than 10 16 [2] Also if we write Z 1 1 λ2n expð λ2 Þdλ ¼ Gn ; then Gn ¼ ð2n 1ÞGn 1 for n Z 1: 2 0 pffiffiffiffi Since G0 ¼ ð1=2Þ π ; we have pffiffiffiffi pffiffiffiffi pffiffiffiffi pffiffiffiffi π and G4 ¼ 105 π: G1 ¼ 14 π ; G2 ¼ 38 π ; G3 ¼ 15 16 32 Employing the expression for ρ in Eq. (15) into the equations of motion (16) and simplifying, using the above integrals, we get hpffiffiffi i 1 Δα0 ¼ Δ1 μ 16π 8T 0 þ 4T 2 þ 6T 4 þ 15T 6 þ 105 2 T 8 þ 2ðT 1 þ T 3 þ2T 5 þ6T 7 Þ ; 2 pffiffiffi 3 π 105 ð17Þ 16 8Z k0 þ 4Z k2 þ 6Z k4 þ 15Z k6 þ 2 Z k8 5 4 ; Δαk ¼ Δ1 þ 12ðZ k1 þ Z k3 þ2Z k5 þ6Z k7 Þ with αk þ 4 ¼ βk for k ¼1–4. where Z ki ¼
i X Sj Ri j j¼0
T 2i ¼
si
;
for i ¼ 0; 1; 2; …:; 7
i X W j R2ði jÞ j¼0
s2i
and T 2i þ 1 ¼
i X W j R2ði jÞ þ 1 j¼0
s2i þ 1
for i ¼ 0; 1; 2; 3; 4
with R0 ¼1. 6. Numerical results The four equations (11) and (17) are programmed in double precision arithmetic, to compute the variations in αi and βi (i¼0–4) at the end of each revolution. Once αi and βi are known, they are converted into the position and velocity vectors, which are further converted into the osculating orbital elements. The KS elements equations of motions are numerically integrated with the analytical density model by the fourth order-Runge–Kutta–Gill method with a small step size of half degree in the eccentric anomaly. Due to air drag the maximum changes occur in the important orbital parameters semimajor axis (a) and eccentricity (e). Hence, we study the changes in these parameters. The values of the gravitational constant, Earth's equatorial radius, rotational rate of the atmosphere and the ballistic coefficient (bn ¼m/cDA) are chosen as 398600.8 km3/s2, 6378.135 km, 1.2 and 50 kg/m2, respectively. Table 1 Decrease in semi-major axis after 100 revolutions for orbits having Hp ¼ 200 km and i ¼15o. e
0.05 0.1 0.15 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
NUM (km)
20.053 16.656 16.267 16.967 20.380 26.969 38.678 60.980 110.193 252.880 1029.191
Analytical (km)
Percentage error
KSC
KS
KSC
KS
20.066 16.664 16.274 16.973 20.380 26.969 38.678 60.981 110.195 252.887 1029.248
19.876 16.542 16.173 16.879 20.387 26.976 38.687 60.996 110.228 253.014 1030.959
0.068 0.048 0.040 0.037 0.0004 0.0008 0.0012 0.0015 0.0020 0.0029 0.0056
0.088 0.068 0.058 0.052 0.032 0.026 0.024 0.026 0.032 0.053 0.172
220
P.S. Rajendran, M.X. James Raj / Acta Astronautica 110 (2015) 214–224
Table 2 Decrease in eccentricity after 100 revolutions for orbits having Hp ¼200 km and i ¼15o. e
0.05 0.1 0.15 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
NUM 108
Analytical 108
26.312 20.137 17.687 16.406 15.148 14.758 14.719 14.873 15.144 15.492 15.893
Percentage error
KSC
KS
KSC
KS
26.329 20.147 17.694 16.412 15.148 14.758 14.720 14.873 15.144 15.492 15.894
26.043 19.996 17.583 16.321 15.153 14.761 14.723 14.877 15.149 15.500 15.921
0.065 0.047 0.041 0.037 0.0007 0.0014 0.0018 0.0022 0.0027 0.0036 0.0063
1.019 0.699 0.584 0.519 0.032 0.027 0.025 0.026 0.033 0.054 0.175
1.8
Difference in semi-major axis(km)
1.6 KSC KS
1.4 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
Initial eccentricity
Fig. 1. Difference between numerically and analytically computed values of semi-major axis (km) after 100 revolutions for orbits having Hp ¼ 200 km and i¼ 15o.
KSC KS
Difference in eccentricity x 10000
0.25 0.2 0.15 0.1 0.05 0 -0.05 0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
Initial eccentricity Fig. 2. Difference between numerically and analytically computed values of eccentricity 104 after 100 revolutions for orbits having Hp ¼200 km and i¼ 15o.
To show the effectiveness of the present solution, simulations are carried out for the orbits with variations in eccentricity (0 o e o 1). The comparison of the results for the important orbital parameter semi-major axis and eccentricity are made with the existing KS theory. Tables 1 and 2 provide the decay in semi-major axis and eccentricity,
P.S. Rajendran, M.X. James Raj / Acta Astronautica 110 (2015) 214–224
221
Difference in semi-amjor axis (km)
5
4
3
2
KSC KS
1
0
-1 10
20
30
40
50
60
70
80
90
inclination (degrees)
Fig. 3. Difference between numerically and analytically computed values of semi-major axis (km) after 1000 revolutions for the orbit with the initial perigee height of 175 km, eccentricity 0.8 with variation in inclination.
respectively for the numerically integrated values (NUM), analytical solutions of the present theory (KSC) and KS theory (KS) after 100 revolutions for the orbits with the perigee height of 200 km, inclination of 15o and eccentricity varying from 0.01 to 0.9. From Table 1, it is noticed that the present solution provides better results for all values of eccentricities with the maximum deviation of only 58 m after 100 revolutions in semi-major axis over a decay of 1029.2 km, where the KS theory provides the difference of 1.8 km. It is also observed that the percentage errors are also less for the present solution comparing to the KS theory. It is true that the errors for small eccentricity as well as high eccentricity will be more, because small eccentricity orbit is near circular orbit and its decay due to drag is uniform throughout the orbits. But in eccentricity orbit the drag attraction is more at perigee, which results in past decay in apogee. Hence, when eccentricity is more, the decay is high in apogee, but perigee will remain same. This will result in accuracy of the prediction also. In Table 2, the same pattern is noticed for eccentricity also. Figs. 1 and 2 provide the differences between the present solution (KSC) and KS solutions (KS) with numerically integrated values for semi-major axis and eccentricity, respectively after 100 revolutions for the same orbits. The continuous lines indicate the present solutions and dotted lines indicate KS solutions. The figures clearly show the superiority of the present solution over the KS theory for all values of eccentricity. To know effect of inclinations, simulations are carried out up to 1000 revolutions for the orbits with the initial perigee height of 175 km, eccentricity 0.7 and variations in inclinations. Fig. 3 provides the differences between the present analytical solutions (KSC) and KS solutions (KS) with numerically integrated values for semi-major axis, which is a measure of energy of the orbit, after 1000 revolutions with inclination varies from 11 to 901. From the figure, it is understood that the difference increases for higher inclinations and the differences are very less for the present analytical theory (KSC) comparing to KS theory (KS) for all inclinations. The same pattern is observed for all inclinations in eccentricity also. From the above two tables as well as the figures, it is clear that the present analytical solution is matching quite well with the numerical values for all eccentricities and found to be better than KS theory, which shows the superiority of the present analytical theory for low as well as high eccentricity orbits.
7. Conclusion The uniformly regular KS canonical element equations are integrated analytically by the series expansion method with air drag, using an oblate exponential atmospheric model. A fourth-order non-singular analytical theory for the contraction of satellite orbit is developed in terms of λ and c. Due to symmetry in the equations of motion, only two of the nine equations are solved analytically to compute the state vector and change in energy at the end of each revolution. Numerical comparison of the present solution with the numerically integrated values for the important orbital parameters semi-major axis and eccentricity shows the accuracy of the present solution. Comparisons with the KS theory up to 1000 revolutions shows the superiority of the present solution, which establishes that the present analytical theory is one of the best theories available in literature for analytical orbit predictions. This theory can be utilized for predicting the orbit or life time of space object in low Earth orbit.
222
P.S. Rajendran, M.X. James Raj / Acta Astronautica 110 (2015) 214–224
Appendix A 9 8 3e2 þ4c2 þ 16 e2 c4 16c2 64 aQ 2i > > > > = < 2 2 2 4 2 1 21e þ 12 c þ 4 e þ c þ 16c þ64 mP ; G00i ¼ 64 2i > > > > ; : 8 3e2 þ 2c2 þ 8 emP 0i ( 2 1 G01i ¼ 64 5e þ 6c2 þ 24 e2 c4 16c2 64 aQ 0i 16 3e3 þ 2c2 þ 8 e mP 2i
G02i
½ 35e2 þ 18c2 þ 72 e2 þ c4 þ 16c2 þ 64mP 0i " 2 # e þ c2 þ 4 e2 aQ 2i 1 ¼ 16 ; 7e2 þ3c2 þ 12 e2 mP 2i 4 2e2 þc2 þ4 emP 0i 4
e ½ðaQ 2i 7mP 2i Þe 8mP 0i ; G04i ¼ 64
4
e G05i ¼ 128 ðaQ 0i 7mP 0i Þ;
ce 2aeQ 2i þ 2e2 c2 8 Q 0i a 6emP 2i 6ce2 mP 0i ; 16
1 2 2 4e c þ 8c 2eaQ 2i 8e2 c2 8 caQ 0i 12ce3 mP 2i 8e2 c2 8 cmP 0i ¼ 16 ( ) 2 2 3 1 2 10e þ c þ8 caQ 2i 4ac Q 0i 2 ; ¼ 16 2 2e c2 8 cmP 2i þ 3ð3c 2Þe2 2 c2 þ 8 c emP 0i 2 " 2 # 2 2 c 2 5e c 8 aeQ 2i þ 5e c 8 aQ 0i 2 ; ¼ 16 4 c þ 8 emP 2i 15e2 þ c2 8 mP 0i 2 " # 2 ce 10aeQ 2i þ c þ 8 2e aQ 0i þ 2ð2eþ 17ÞemP 2i ; ¼ 16 þ2 5e2 þc2 þ8 mP 0i
G10i ¼ G11i G12i G13i G14i
ce2 ð6aeQ 2i þ 5aQ 0i þ 28emP 2i þ 17mP 0i Þ; 16 3 ce G16i ¼ ðQ 3a þ 14mP 0i Þ; 16 0i c2 e2 G21i ¼ ð11Q 0i a mP 0i Þ; 64 2 c e G22i ¼ ð11eaQ 2i 8Q 0i a emP 2i 4mP 0i Þ; 32 1 96c2 eQ 2i a þ 69e2 2c2 24c2 Q 0i a þ48c2 emP 2i þ 15e2 þ 96e 2c2 24 c2 mP 0i G23i ¼ 192 1 51e2 c2 12 c2 Q 2i a þ 63e2 24e c2 12 mP 2i 12c2 emP 0i G24i ¼ 48 1 96c2 eaQ 2i þ 45e2 2c2 þ 24 c2 Q 0i a þ ð48ðeþ 2Þc2 emP 2i þ 39e2 þ 2c2 þ24 c2 mP 0i G25i ¼ 192 c2 e ð19aeQ 2i þ8aQ 0i þ 39emP 2i þ 12emP 0i Þ; G26i ¼ 32 2 2 c e G27i ¼ ð19aQ 0i þ 39mP 0i Þ; 64 3 c e G34i ¼ ð3aQ 0i þ2mP 0i Þ; 16 3 c G35i ¼ ð6eaQ 2i aQ 0i þ 4emP 2i mP 0i Þ; 48 c3 G36i ¼ ðaQ 2i þ mP 2i þ emP 0i Þ; 24 c3 G37i ¼ ð6aeQ 2i þ aQ 0i þ 8emP 2i þ mP 0i Þ; 48 c3 e ½3aQ 0i þ 4mP 0i ; G38i ¼ 48 c4 4 G47i ¼ c ðaQ 0i þ mP 0i Þ 384 1 c4 ðaQ 0i þmP 0i Þ; G48i ¼ ðaQ 2i þ mP 0i Þ; G49i ¼ 192 384 G15i ¼
P.S. Rajendran, M.X. James Raj / Acta Astronautica 110 (2015) 214–224
H 10i ¼
ce 16
c2 þ 8 aQ 1i þ 2c2 þ 16 11e2 mP 1i
H 11i ¼
c 10e2 c2 8 aQ 1i ð14e2 þc2 þ8ÞmP 1i 16
H 12i ¼
ce 2 4e þ3c2 þ 24 m þ 2 5e2 þ2c2 þ 8 a P 1i 16
H 13i ¼
c 15e2 c2 8 aQ 1i þðce2 c2 8ÞmP 1i 16
H 14i ¼
ce 2 8e c2 8 aQ 1i þ 2 7e2 c2 8 mP 1i 16
H 15i ¼
ce2 ð7aQ 1i þ mP 1i Þ; 16
H 22i ¼
c2 e ð2aQ 1i þ 3mP 1i Þ; 4
H 25i ¼
H 27i ¼
H 16i ¼
H 23i ¼
c3 e ½aQ 1i þ mP 1i ; 48
H 47i ¼
c4 ðaQ 1i þmP 1i Þ; 384
e2 2 7e þ 3 c2 þ4 ; 16
B04 ¼
7e4 ; 64
B12 ¼
ce 2 5e c2 8 ; 4
B15 ¼
ce3 ; 8
B24 ¼
1 2 c e; 2
B36 ¼
1 3 c e; 12
B22 ¼
3 2 2 c e ; 32
B25 ¼
ce2 ð3e 19Þ; 8
B26 ¼
c2 ðaQ 1i þmP 1i Þ; 48
c3 e ð3aQ 1i þ4mP 1i Þ; 8
B00 ¼
e 2 3e þ 2 c2 þ 4 ; 4
B03 ¼
e3 ; 4
ce2 ðe 13Þ; 8
B23 ¼
35 2 2 c e ; 32
H 35i ¼
H 38i ¼
e2 2 7e þ 3 c2 þ 4 ; 16
B11 ¼
B13 ¼
3c2 e2 ðaQ 1i þ mP 1i Þ; 16
c2 e ð2aQ 1i þ 3mP 1i Þ; 8
c3 e ð3aQ 1i þ2mP 1i Þ; 8
c4 ðaQ 1i þ mP 1i Þ; 384
B02 ¼
ce2 ðeþ 3Þ; 8
H 26i ¼
c3 ½aQ 1i þ mP 1i ; 48
H 49i ¼
B01 ¼
B10 ¼
H 34i ¼
H 37i ¼
H 21i ¼
c2 36ae2 Q 1i þ ð135e2 2c2 24ÞmP 1i 16
c2 60ae2 Q 1i þ 195e2 2c2 24 mP 1i ; 16
c2 e2 ð19aQ 1i þ39mP 1i Þ; 64
H 36i ¼
ce3 ð2aQ 1i þ13mP 1i Þ; 16
B14 ¼
3 2 ce ð5eþ 1Þ; 8
29 2 2 c e ; 32
3 2 2 c e ; 32
All other coefficients are zero:
Appendix B 3 2
W 0 ¼ ð1 þ eÞ1 ; ð1 eÞ2
W 1 ¼ 3e 8e þ11; 8ð1 e2 Þ2 2
W2 ¼
ð5e4 16e3 50e2 þ 16e 3Þ
W4 ¼
ð45e8 160e7 4228e6 4128e5 418e4 þ 32e3 228e2 þ 160e 35Þ
1
128ð1 þ eÞ2 ð1 e2 Þ2
;
W 3 ¼ ð7e
6
1
32768ð1 þ eÞ6 ð1 e2 Þ2
24e5 243e4 80e3 þ 45e2 24e þ 5Þ 1
1024ð1 þ eÞ4 ð1 e2 Þ2
:
;
223
224
P.S. Rajendran, M.X. James Raj / Acta Astronautica 110 (2015) 214–224
So ¼ ðR0 þ R2 þ R4 ÞK; S1 ¼ eR1 þR3 þR5 þ R6 ; h i ð 15e2 þ 4e þ 15ÞR0 þ ð 3e2 þ 4e þ 3ÞR2 þ ðe2 þ 4e 1ÞR4
S2 ¼
8ð1 þ eÞ2
ð3e
K;
3e 3eÞR1 3eR3 þ ð e 3e þ 1ÞR5 þ ð2e 3e 2ÞR6 ; 6ð1 þ eÞ2 3 2 4 3 2 4 35e 88e 54e þ 120e þ35 R0 þ 5e þ 8e3 þ 26e2 þ 24e 5 R2 6 7 þ 3e4 þ 40e3 þ 10e2 8e þ3 R4 6 7 7K; S4 ¼ 6 4 6 7 128ð1 þ eÞ 4 5
S3 ¼
3
2
2 6 6 S5 ¼ 6 6 4
2
2
3 15e3 30e2 R1 þ 30e3 15e2 R3 þ 4e4 þ 40e3 7e2 þ 10e 4 R5 4 7 þ 6e 10e3 27e2 20eþ 6 R6 7 7; 4 7 120ð1 þ eÞ 5
2 21e6 212e5 175e4 þ 264e3 þ 303e2 þ 140e 21R 3 0 6 þ 7e6 þ12e5 þ101e4 þ 200e3 þ27e2 20e þ 7R2 7 7 6 6 7 6 7 þ 5e6 þ 172e5 þ 129e4 þ8e3 e2 þ 12e 5 R4 6 7 S6 ¼ 6 7K; 1024ð1 þ eÞ6 7 6 7 6 6 7 4 5 2 6 6 6 6 6 S7 ¼ 6 6 6 6 4
3 105e5 315e4 105e3 R1 þ 210e5 210e4 105e3 R3 7 þ 12e6 308e5 209e4 þ 21e3 e2 28e þ12 R5 7 7 7 þ ð16e6 28e5 188e4 329e3 22e2 42e 16ÞR6 7 7; 6 1680ð1 þ eÞ 7 7 7 5
2 99e8 3056e7 3884e6 þ 1424e5 þ 6802e4 þ 7088e3 þ 724e2 336eþ 99R 3 0 6 þ 45e8 þ 80e7 þ 1556e6 þ 4816e5 þ 2354e4 þ112e3 þ 20e2 þ 112e 45R2 7 6 7 6 7 6 þ 35e8 þ2832e7 þ 3540e6 þ 2192e5 þ 274e4 þ 176e3 44e2 80e þ 35 R4 7 6 7 S8 ¼ 6 7K; 32768ð1 þ eÞ8 6 7 6 7 7 6 4 5 1 2
with K ¼ ð1 eÞ1 ð1 þ eÞ2
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]
D.G. King-Hele, Satellite Orbits in an Atmosphere: Theory and Applications, Blackie, 1987. T.E. Sterne, An Introduction to Celestial Mechanics, Interscience, New York, 1960. D.G. King-Hele, Theory of Satellite Orbits in an Atmosphere, Butter Worths, London, 1964. G.E. Cook, D.G. King-Hele, D.M.C. Walker, The contraction of satellite orbits under the influence of air drag. II with oblate atmosphere, Proc. R. Soc. A 264 (1316) (1961) 88–121. G.E. Cook, D.G. King-Hele, The contraction of satellite orbits under the influence of air drag IV with scale height depending on altitude, Proc. R. Soc. 275 (1362) (1963) 357–390. D.G. King-Hele, The contraction of satellite orbits under the influence of air drag III. High-eccentricity orbits (0.2 re o1), Proc. R. Soc. A 267 (1331) (1962) 541–557. P. Kustaanheimo, E.L. Stiefel, Perturbation theory of Kepler motion based on Spinor regularization, J. Reine Angew. Math. 218 (1965) 204–219. E.L. Stiefel, G. Scheifele, Linear and Regular Celestial Mechanics, Springer, Berlin, 1971. O.F. Graf Jr., A. Mueller, S. Starke, The Method of Averages Applied to the KS Differential Equations, NASA-CR-151607, 1977. R.K. Sharma, M. Xavier James Raj, Long-term computations with ks uniformly regular canonical elements with oblateness, Earth Moon Planets 32 (1988) 63–178. M. Xavier James Raj, R.K. Sharma, Analytical short-term orbit predictions with J2, J3, J4 in terms of K–S uniformly regular canonical elements, Adv. Space Res. 31 (2003) 2019–2025. M. Xavier James Raj, R.K. Sharma, Contraction of near-Earth satellite orbits using uniformly regular KS canonical elements in an oblate atmosphere with density scale height variation with altitude, Planet. Space Sci. 57 (2009) 34–41. R.K. Sharma, M. Xavier James Raj, Analytical approach using KS elements to high eccentricity orbit prediction including drag, Acta Astronaut. 66 (2010) 864–870. R.K. Sharma, M. Xavier James Raj, Fourth-order theories for orbit predictions for low-and high-eccentricity orbits in an oblate atmosphere with scale height dependent on altitude, Acta Astronaut. 65 (2009) 1336–1344.