- Email: [email protected]
Decay of satellite orbits using K-S elements in an oblate diurnally varying atmosphere with scale height dependent on altitude
Available
Pergamon
online SCIENCE
www.elsevier.com/locate/asr
at www.sciencedirect.com DIRECT*
doi: 10,1016/SO273-1177(03)00153-4
DECAY OF SATELLITE ORBITS USING K-S ELEMENTS IN AN OBLATE DJURNALLY VARYING ATMOSPHERE WITH SCALE HEIGHT DEPENDENT ON ALTITUDE L. S. Nair’, and R.K. Sharma’ ‘Department 2Applied
of Mathematics, Mathematics
H H M S P B N S S College, Neeramonkara ,Thiruvananthapuram Division, Vikram Sarabhai Space Centre, Thiruvananfhapuram
- 695040, India. - 695022, India.
ABSTRACT
A non-singular analytical theory for the motion of near-Earth satellite orbits with the air drag effect is developed in terms of the K-S elements utilizing an analytical, oblate diurnally varying atmospheric model with varying scale height, dependent on altitude, The series expansions include up to fourth-order terms in eccentricity (e) and c, a small parameter dependent on the flattening of the atmosphere. Only two of the nine equations are solved analytically to compute the state vector at the end of each revolution due to symmetry in the K-S element equations. Numerical studies are done over a wide range of orbital parameters. A numerical comparison with numerically integrated values of the change in the orbital parameters: semi-major axis and eccentricity is made with the extended theory up to fourthorder terms of Swinerd and Boulton. It is observed that the theory in terms of the K-S elements provides better accuracy than the extended theory of Swinerd and Boulton, when compared with the numerically integrated values. 0 2003 COSPAR. Published
by Elsevier
Science Ltd. All rights reserved.
INTRODUCTION
It is well known that the Earth’s atmosphere is non-stationary and does not possess spherical symmetry and varies, in general, considerably during the lifetime of a satellite. The change in the distribution of atmospheric density with time is primarily caused by the variability of the ultraviolet and corpuscular fluxes from the Sun. As these fluxes change relatively slowly with time, the atmosphere can be considered stationary at time intervals comparable to the satellite orbital period around the Earth. The non-spherical nature of the atmosphere is primarily due to Earth’s gravitational field and the heating of the illuminated part of the Earth’s atmosphere by solar radiation. The atmospheric density at a given distance from the center of the Earth will depend upon the geographic latitude and longitude, and on the Sun’s location in the celestial sphere. Heating of the illuminated part of the Earth’s atmosphere by the Sun results in a deformation in the upper atmosphere and an acquired pear shape directed in an acute ring approximately towards the Sun. Thus the selection of an appropriate atmospheric model has great significance in the study of the perturbations of satellite motion caused by air drag. In order to account for the drag perturbations, it is necessary to accurately model the density of the Earth’s atmosphere. Some of the important atmospheric models used frequently in the literature are CIRA (1972), Jacchia (1977), Barlier et al. (1978), MSIS-86 (Hedin, 1987) and MSIS-90 (Hedin, 1991). MSIS-90 uses analytic models to model the lower altitudes to account for disturbances such as solar activity, magnetic storms, and daily variations, as well as latitude, longitude, and monthly variations. Analytical solutions of artificial satellite motion in an atmosphere, which considers the non-sphericity of the distribution of the atmdspheric density in which diurnal and atmospheric oblateness effects are treated in combination
Adv. Space Res. Vol. 31, No. 8, pp. 201 l-2017,2003
Q 2003 COSPAR. Published by Elsevier Printed in Great Britain 0273-l 177/03 $30.00 + 0.00
Science
Ltd. All rights
reserved
2012
L. S. Nair and R. K. Sharma
were carried out by Santora (1975), Swinerd and Boulton (1982) and a few others. Swinerd-Boulton’s theory is up to third-order terms in e and c. In both, Swinerd-Boulton and Boulton (1983) theories, the variation of density scale height with altitude was included. Sharma (1997) developed an analytical solution in terms of the K-S elements up to third-order terms in e and C, by including the effects of the atmospheric oblateness and diurnal bulge. However, the density scale height was assumed to be constant. In this paper, we extend Sharma’s theory up to fourth-order terms in e and c and include the variation of the density scale-height with altitude. Only two of the nine equations need to be solved analytically to compute the state vector at the end of each revolution due to symmetry in the K-S element equations. Swinerd and Boulton, and Boulton solutions have also been extended up to fourth-order terms in e and c for comparison purpose. These solutions are compared numerically with the present solution obtained in terms of the K-S elements. Numerical comparisons are made with four test caseshaving perigee altitudes of 175,240,300 and 400 km, with each case having low and high inclination of 5’ and 80’ and with a high value of 0.2 for the rate of increase of the density scale height H with the distance r. At each perigee height, two values of e are chosen as 0.1 and 0.2. The series expansions and simplifications are done using the software MATHEMATJCA (Wolfram, 1996), available at V&ram Sarabhai Space Centre, Thiruvananthapumrn. As our aim is to compare the numerical and analytical solutions, we have utilized a relatively simple, Jacchia (1977) atmospheric model to compute the density and density scale height values. Two distinct advantages of the present solution are seen over the theory of Swinerd and Boulton, and Boulton. In the present theory only two equations are handled analytically against three by Swinerd and Boulton, and Boulton. The present solution is non-singular, whereas in Boulton’s expression of argument of perigee variation (Ao), e is in the denominator. EQUATIONS
OF MOTION
The K-S element equations of motion of a satellite under the effect of additional perturbing force F are (Stiefel and Scheifele, 1971) dw -=--&z*,LTP,, dE
(1) (2)
u’ = (u1,u2,u3,u4)
z=t+&ii*), W
=iEcos(f)+gsin(f),
u”
j;=(x1,x2,x3)=L(ii)ii, r=
B,
-
a3
P,
+ a2 ru*
and
L(G) =
B3 -
-u2
al -u3
u*
u, -uq
u3
u4
Lu4
-u3
Ul u2
P4
(4)
= 0, u4
-u3
1
I
(5)
u2 -%j
where Eq. (4), is the bilinear relation and E, w, t, r and K* are, respectively, the eccentric anomaly, angular frequency, physical time, radial distance and the gravitational constant. Knowing the position and velocity vector j; and 1 at the instant t = 0, the values of r, w, z, t& and h* can be computed and by adopting E = 0 as the initial value of the eccentric anomaly; we obtain
Jf 13 is the aerodynamic drag force per unit mass on a satellite of mass m, (King-Hele, D.G, 1987) F =[1-rp,,Acosi,lv,,]2, F = -ipSliJlG, S=FAC,fm,
(6)
Decay
of Satellite
Orbits
Using
K-S Elements
2013
where p is the atmospheric density, A is the rotational rate of the atmosphere about the Earth’s axis, rpo is the initial perigee radius, &,is the initial inclination, vpois the velocity at the initial perigee, Cu and A are, respectively, the drag coefficient and the effective area of the satellite. If the density, p is assumed to vary sinusoidally with 4, where 4 is the geocentric angular distance from the direction of the density maximum, F= LI&=”
P=Po(l+Fcosp)exp{-p(r-a)},
p, ,B=$, (7) 0 !lU” where p o is the average density on the reference spheroid when 9=90° and H,, is the average density scale height. Following (Cook et al, 1961) we choose o with 8, geocentric latitude, as (l-&sin* 0) CT=rpo (l-&sin’ QPO) ’ The ellipticity of each of the oblate spheroidal atmospheric surfaces is assumed to be the same as that of the Earth’s ellipticity, 0.00335. The density p, in Eq. (7) can be written as 1
p=k
{ 1 -I
F A(cosE-e)
(1-ecosE) x[l+ccos(2o+2E)
+ FB(l-e*)
:sinE
(1-ecosE)
~}xexp(zcosE)
+ce(cos(2o+3E)-cos(2w+E)}
+ce2{cos2w-4cos(20+2E)+3cos(20+4E)}/4 -cos(2@+5E)}/2+ -cos(4w+3E)}/2
-ce3{cos(2w+3E)
c2 {l+cos(4w+4E)}/4+c”
e{cos(40+5E)
+c2e2 {3cos(4w+2E)-8cos(4w+4E)
+5cos(4o+6E)}/8
+c3{3~~~(2~+2E)+~~~(6~+6E)}/24
-c3e{cos(2w+E)-cos(2u+3E)+cos(6w+5E)-cos(6w+7E)}/8 + c4 (3+4cos(4w+4E)+cos(8~+8E)}/192 k=Ppoexp{P(rP,-a)-ccos2w, A=sinS,
sinisinw
+cos~,{cos(S2-a,
+O(ce4,c2
e3,c3e2,c4 e,ce4,c5)1, 1 }, c=&$$,sin’ -i, 2 )cosw-cosisin(G-oc, )sinw},
(8)
B=sinS, sinicoso-cos6,{cos(SZ-ol,)sinw+cosisin(SZ-~:,)cosw], (a,, 6,) and (a~ , &,.,), where oM = a, +h and & = 6, are the right ascension and declination of the Sun and the centre of the day time bulge, i.e. the bulge lags behind in right ascension by an angle h. w and Q are argument of perigee and right ascension of ascending node of the orbit. In reality the scale height is known to vary with altitude, and a satellite will sample the atmosphere over a significant height range, even when the orbit is of low eccentricity. Similar to that of Cook and King-Hele (1963), we adopt an expression, H = Ho +p (r-o), (9) in which the scale height H varies linearly with height where p, the altitude gradient of H,, , is assumed to be constant in the range lpl< 0.2. Now, by definition, H and p are related through the equation - 1 = --- 1 dP(T) p(r)dr’ . H(r) Hence, if we substitute Eq. (9) in to Eq. (lo), the density is p(r)
=pl (l+~(r-a)l~, HO
(10)
(11)
2014
L. S. Nair and R. K. Sharma
Ignoring the terms of order pz and replacing p , by ‘p O (1 + F cos cp),the density equation p(r)= { l+p (r-o)‘/2H,” )p, (l+Fcos#)exp{-(r-a)/Ho }, 02) allows an approximately linear altitude variation in the density scale height. Using r = a (1 - e cosE) and z = ae / II,,, Eq. (12) can be expressed in terms of E as p(r) =~p{l++(l-
2cosE+cos2E)+O(~c)).
(13)
Following the terminology, notations and procedure as in Sharma (1991, 1992) s=-pc!TK2
[{fo+5e2(2fo+fl)/8+27e4(3fo+2f,)}+{e(2fo+fl)/64
+3e3 (3f, + 2f,)/4}cosE+(e+3
e3/4) f, sinE+{J
+27e4 (Sf, +7f, )/128}cos2E+(f, +3e3(3f,+2fO)/8}cos3E+
+5e2 /4+135e4
+5e2 (f, +fi )I4 /128)sin2E+{ef,
(e+9e3/8)f2sin3E+{5e2f,/8+27e4/64
+(f,+2f,))cos4E+(5e2/8+27e4/32)f2sin4E+3e3f,cos5E/8 3
+3e3 f2sin5E/8+27e4 ANALYTIC
fIcos6E/128+27e4
f2sin6E1128]/16a2
.
(14)
INTEGRATION
The explicit form of the perturbations Aw is obtained by substituting the value of p from Eq. (8) in to Eq. (14) and noting that the non-zero integrals in the resulting expression are of the form of linear combinations of the Bessel functions 2II
Z,(z) =&
j exp (z cos E) cos nEdE . 0
After some algebra, we obtain the change in w after one revolution as: Aw = ( C$/ a ) [A,,+ eA, +c{coQo A2 +sin2o A3)} +e2A4 +ce{cos2 o A5 + sin2 o A,} + c* {A, +cos4 o Aa +sin4 o AP} +e3Alo + ce2 (cos2 o A,, +sin2 o A12} + c2e { Al3 +cos40 Al4 + sin4o Ar5} +c3 {COS~OAt6+sin20 At7 + coda Ata + sin6e.rA19}+e4AU,+ ce3 { cos20 A2i+ sin2o AZ} +c2e2{ Aa +cos4u1Au +sin4o AZ} +c3e {COS~OA+ sin2o A+ coda Aa +sin6o A29 }+ c4 {Aso + COS~WA31 +sin4o A~~+COS~OI A33 +sin8w A34 }+ O(r) +FA [A35+ eA36 +c ( cos20 A3, +sin2o A3s)} +e2Ajs ice{ cos2k AN + sin2o A41} + c2 { h2 +cos4w A43 +sin4o Au} +e3AU + ce2 { cos2o A& +sin2o &T) + c2e {AM +cos4w A49 + sin40 Aso} +c3 (cosiw ASI +sin2o A52 + coda A& sin6o AN} +e4Ass+ ce3 { cos2w As6+ sin2w AZ7) +c2e2 {Ass +cos40 A59+sin4u) AW}+c3e { cos20 A61 +sin2o A,~~+cos~u A63}+sin60 A&c’{ A&cos40 AM +sin4o A67+cos80 Abs+sin8w A69 )I +FB [ATo+ eA,, +c{ cos20 AT2+sin2o AT3)} +e2AT4+ce{ cos2o AT5+ sin2o AX} + c2 (A,, +cos40 AT8+sin4w ATg} +e3Aso + ce2 (~0~20 Asr +sin2o Aaz}+ c2e (As3 +cos40 As4 + sin4w Aas} +c3 {costi Aw+sin2o Aa7 + cos60 Aas+ sin6w Asq} +e4Aw +ce3{cos2ur Ag,+ sin2w Ag2}+c2e2{Ag3+cos40 Ag4+sin4ccAg,}+c3e (~0~20 Ag6+sin2cuA97 +cos60 AgS}+sin6o Agg+c4{AlW+ cos40 Alol +sin4@ Alo*+ cos80 Aro3 +sin8 o A104 III , = (@la)A,
(13
2015
Decay of Satellite Orbits Using K-S Elements
where Ap ,q=O,l,....,, 104, are linear combinations of the modified Bessel functions IO(z),n = O,l,. . . ...16. The change in the state variables given by Eq. (3) after one revolution is A& =(@/w)h,
i=l.....
8,
(17)
where the coefficients A, in Eq. (16) change in Eq. (17). The changed coefficients are also linear combinations of I,(z), n =O,l,.., 16, and q~=-nkK6
exp(-pa)/80i.
When we consider the scale height variation, p is replaced by Eq. (13). Then the change in w and & are given by Aw, = Aw + (l/2) Joz2 N(Aw) +0(r). where v denotes the assumption of the variability in H, and N is an operator defined by N(1,) = I(O,n) - 2 I(l,n) + I(2,n). Here I(p,n) = 2’Z PC,In+pZr,r = 0, p; p = 0,1,2, where pC, are the usual binomial coefficients. By a similar process, we have A& = A& +( l/2) p Z2 N(A&) +0(r). NUMJZRICAL RESULTS We have programmed Eqs. (1) and (3) in double precision arithmetic to compute the K-S elements Q and pi. Then the position, velocity and the orbital elements can be easily computed. Four test caseswith perigee heights of 175, 240, 300 and 400 km, where the effects of the atmospheric oblateness and the diurnal bulge are significant, are chosen for detailed numerical study. At each perigee height two different inclinations of 5’ and 80’ are chosen to cover a wide range of inclination, and also to vary the parameter c, which is used in series expansions, from a very small value of I
Table 1. Decrease in semi-major axis (km) after 250 revolutions with percentage error orbit Parameters u = 0.2 I Method HP1 e I a I -I (km) NUM ( KSANAL4 [ % error ( ANAL4 1 % error 7279.78 93.9459 1 90.6000 1 3.56 1 89.1875 1 5.07 175 0.1 8183.50 95. 2295 I 94.1912 I 1.09 0.2 I 91.2925 I bTii 1 15.7828 -__- 1 4.46 ---240 0.1 7352.00 16.5201 1 16.0057 ( 3% 8271.00 0.2 0.47 16.0945 ~1 s.05 16.9510 1 16.8712
I
l
5
] 101.2190 1 34.94 1 106.6745 1 105.0527 1 80
3.93
1 99.6059
1.52
-_--
[ 5.5114 1 7525.50 1 0.9361
1 5.4770 0.9025
1
0.62
1 5.1622 I
0.8
1
.--5.46
I
]
1
2016
L. S. Nair and R. K. Sharma
0.002 to a.large value of 0.39, At each perigee height we have computed the change in semi-major axis, by which the energy of an orbit can be calculated, after 250 and 500 revolutions given in Tables 1 and 2, with reasonably high values of e (0.1, 0.2). e is a very important parameter in the present solution, as lot of series expansions are involved with it. In the entire test cases reported here, the values of O, &Z and mean anomaly are 60, 30 and 0 degrees, respectively. The density scale height variation with altitude is included and a high value of p = 0.2 is taken to test the analytical solution more rigorously. The effect of the flattening of the Earth has been included in computing the perigee and apogee heights. Tables 1 and 2 provide the change in semi-major axis, obtained by the present theory (KSANAU), Swinerd and Boulton method (ANAL4) and by numerical integration (NUM) of the K-S element Eqs. (1) and (3) with a fixed step size fourth order Runga-Kutta-Gill method having a very small step size of half degree in the eccentric anomaly by using the expressions for p given in Eq. (12). Jacchia (1977) atmospheric density model, which is relatively easier to use, is employed to compute the value of p. and I&, being used in the Eq. (13). Arbitrarily 22 August 2002 is chosen as the initial epoch. The values of solar flux (Fr0.r) and averaged geomagnetic index (A,) Table 2. Decrease in semi-major axis (km) after 500 revolutions with percentage error Orbit Parameters u=o.2
are taken as 150 and 10, respectively, which approximately represents an average density model and results in exospheric temperatures between 1000 and 1100 K, in the present study. The values of A, E and B, = (m /CnA) utilized during the computations are 1.2,0.0035 and 50.0, respectively. It can be noticed that for the epoch under consideration the values of FA and FB in Eq. (16) slightly decrease as the number of revolution increases. It is observed that these values increase as perigee height increases. Percentage error = 100 x (NUM - ANAL) / NUM is calculated with respect to the present solution as well as with Swinerd and Boulton’s theory. It is seen from Tables 1 and 2 that at lower altitudes small increase in percentage error is noticed with KSANAL4, when inclination increases from 5’ to SO’. Also, the percentage error decreases when e increases from 0.1 to 0.2. It is noticed from these tables that the drag perturbed orbital parameter, a, match quite well with NUM values with the present solution. On altitudes above 175 km, the percentage error seems to be less than 3.The maximum percentage error after 250 and 500 revolutions for the test cases under study, is less than 4.0 and 6.3, respectively, whereas with ANAL4 these values are 8.1 and 13.4. The accuracies of the numerical computations examined with the help of the bilinear relation in the K-S elements are found to be very satisfactory. The left-hand side of the bilinear equation, Eq. (4), gave the values of order lo”‘- lo-” with respect to KSANAL4 and NUM computations after 500 revolutions in the present study.
Decay of Satellite Orbits Using K-S Elements
2017
CONCLUSION The K-S element equations are integrated analytically by a series expansion method by assuming an oblate diurnal atmosphere with density scale height varying with altitude and a non-singular solution containing up to fourth-order terms in eccentricity and c (a small parameter dependent on the flattening of the atmosphere) is obtained. Only two of the nine equations are solved analytically to compute the state vector at the end of each revolution due to symmetry in the K-S element equations. Numerical comparison is made with numerically integrated values and with the modified solution of Swinerd and Boulton up to fourth order terms in e and c, with different perigee altitudes for reasonably high eccentric orbits, having small and high orbital inclination, after 250 and 500 revolutions. The present solution is found to provide better accuracies. ACKNOWLEDGMENTS We gratefully acknowledge the referees for their valuable comments and suggestions, which helped in bringing the paper to its present form. REFERENCES Barlier, F., C. Berger, J. L.. Falin, et al., A Thermospheric Model based on satellite drag data, Ann. Geophys., 34,9-24, 1978.
Boulton, W.J., The advance of a satellite in a rotating oblate atmosphere with a diurnal variation in density, Proc. R. Sot. Lord, A 389,433444,1983. Cook, G.E., D.G. King-Hele, and D.M.C. Walker, The contraction of satellite orbits under the influence of air drag. Part II. With oblate tmosphere, Proc. R.. Soc.Lond. A 264,88-121,1961. Cook, G.E., D.G. King-Hele, The contraction of satellite orbits under the influence of air drag. Part IV. With scale height dependent on altitude. Proc. R. Sot. Land. A 215,357-390,1963. Hedin, A.E., MSIS-86 Thermospheric model , J. Geophys. Res., 92(A5), 4649-4662, 1987. Hedin, A.E., Extension of the MSIS Thermosphere model into the middle and lower atmosphere, J. Geophys.Res.,96 (A2), 1159-l 172, 1991. Jacchia, L.G. Smithsonian Astrophysics Observatory Special report, No. 375, 1977. King-Heie, D.G. Satellite Orbits in an Atmosphere: Theory and Applications, Blackie, Glasgow and London, 1987. Santora, F. A., Satellite drag perturbations in an oblate diurnal atmosphere, AZAA J., 13, 1212-1216, 1975. Sharma, R.K., Analytical approach using K-S elements to near-Earth orbit predictions including drag. Proc. R. Sot. Lo&. A 433, 121-130,199l. Sharma, R.K., A third - order theory for the effect of drag on Earth satellite orbits, Proc. R.. Sot. Land. A 438, 467475,1992. Sharma, R.K., Contraction of satellite orbits using K-S elements in an oblate diurnally varying atmosphere. Proc. R. Sot. Lmd. A 453, 2353-2368, 1997. Stiefel, E.L., and G. Scheifele, Linear and Regular Celestial Mechanics, Springer - Veralg, Berlin, Heidelberg, New York, 1971. Swinerd, G.G., and W.J. Boulton, Contraction of satellite orbits in an oblate atmosphere with a diurnal density variation, Proc. R. Sot. Land. A 383,127-145, 1982. Swinerd, G.G., and W.J. Boulton, Near circular satellite orbits in an oblate diurnally varying atmosphere, Proc. R. Sot. Land. A 389, 153-170, 1983. Wolfram, S, Mathematics Version 3, Wolfram Media, Cambridge University Press, USA, 1996. E-mail address of Lila S Nair [email protected] Manuscript received 7 November 2002; revised 9 March 2003, accepted 14 March 2003
Pergamon
online SCIENCE
www.elsevier.com/locate/asr
at www.sciencedirect.com DIRECT*
doi: 10,1016/SO273-1177(03)00153-4
DECAY OF SATELLITE ORBITS USING K-S ELEMENTS IN AN OBLATE DJURNALLY VARYING ATMOSPHERE WITH SCALE HEIGHT DEPENDENT ON ALTITUDE L. S. Nair’, and R.K. Sharma’ ‘Department 2Applied
of Mathematics, Mathematics
H H M S P B N S S College, Neeramonkara ,Thiruvananthapuram Division, Vikram Sarabhai Space Centre, Thiruvananfhapuram
- 695040, India. - 695022, India.
ABSTRACT
A non-singular analytical theory for the motion of near-Earth satellite orbits with the air drag effect is developed in terms of the K-S elements utilizing an analytical, oblate diurnally varying atmospheric model with varying scale height, dependent on altitude, The series expansions include up to fourth-order terms in eccentricity (e) and c, a small parameter dependent on the flattening of the atmosphere. Only two of the nine equations are solved analytically to compute the state vector at the end of each revolution due to symmetry in the K-S element equations. Numerical studies are done over a wide range of orbital parameters. A numerical comparison with numerically integrated values of the change in the orbital parameters: semi-major axis and eccentricity is made with the extended theory up to fourthorder terms of Swinerd and Boulton. It is observed that the theory in terms of the K-S elements provides better accuracy than the extended theory of Swinerd and Boulton, when compared with the numerically integrated values. 0 2003 COSPAR. Published
by Elsevier
Science Ltd. All rights reserved.
INTRODUCTION
It is well known that the Earth’s atmosphere is non-stationary and does not possess spherical symmetry and varies, in general, considerably during the lifetime of a satellite. The change in the distribution of atmospheric density with time is primarily caused by the variability of the ultraviolet and corpuscular fluxes from the Sun. As these fluxes change relatively slowly with time, the atmosphere can be considered stationary at time intervals comparable to the satellite orbital period around the Earth. The non-spherical nature of the atmosphere is primarily due to Earth’s gravitational field and the heating of the illuminated part of the Earth’s atmosphere by solar radiation. The atmospheric density at a given distance from the center of the Earth will depend upon the geographic latitude and longitude, and on the Sun’s location in the celestial sphere. Heating of the illuminated part of the Earth’s atmosphere by the Sun results in a deformation in the upper atmosphere and an acquired pear shape directed in an acute ring approximately towards the Sun. Thus the selection of an appropriate atmospheric model has great significance in the study of the perturbations of satellite motion caused by air drag. In order to account for the drag perturbations, it is necessary to accurately model the density of the Earth’s atmosphere. Some of the important atmospheric models used frequently in the literature are CIRA (1972), Jacchia (1977), Barlier et al. (1978), MSIS-86 (Hedin, 1987) and MSIS-90 (Hedin, 1991). MSIS-90 uses analytic models to model the lower altitudes to account for disturbances such as solar activity, magnetic storms, and daily variations, as well as latitude, longitude, and monthly variations. Analytical solutions of artificial satellite motion in an atmosphere, which considers the non-sphericity of the distribution of the atmdspheric density in which diurnal and atmospheric oblateness effects are treated in combination
Adv. Space Res. Vol. 31, No. 8, pp. 201 l-2017,2003
Q 2003 COSPAR. Published by Elsevier Printed in Great Britain 0273-l 177/03 $30.00 + 0.00
Science
Ltd. All rights
reserved
2012
L. S. Nair and R. K. Sharma
were carried out by Santora (1975), Swinerd and Boulton (1982) and a few others. Swinerd-Boulton’s theory is up to third-order terms in e and c. In both, Swinerd-Boulton and Boulton (1983) theories, the variation of density scale height with altitude was included. Sharma (1997) developed an analytical solution in terms of the K-S elements up to third-order terms in e and C, by including the effects of the atmospheric oblateness and diurnal bulge. However, the density scale height was assumed to be constant. In this paper, we extend Sharma’s theory up to fourth-order terms in e and c and include the variation of the density scale-height with altitude. Only two of the nine equations need to be solved analytically to compute the state vector at the end of each revolution due to symmetry in the K-S element equations. Swinerd and Boulton, and Boulton solutions have also been extended up to fourth-order terms in e and c for comparison purpose. These solutions are compared numerically with the present solution obtained in terms of the K-S elements. Numerical comparisons are made with four test caseshaving perigee altitudes of 175,240,300 and 400 km, with each case having low and high inclination of 5’ and 80’ and with a high value of 0.2 for the rate of increase of the density scale height H with the distance r. At each perigee height, two values of e are chosen as 0.1 and 0.2. The series expansions and simplifications are done using the software MATHEMATJCA (Wolfram, 1996), available at V&ram Sarabhai Space Centre, Thiruvananthapumrn. As our aim is to compare the numerical and analytical solutions, we have utilized a relatively simple, Jacchia (1977) atmospheric model to compute the density and density scale height values. Two distinct advantages of the present solution are seen over the theory of Swinerd and Boulton, and Boulton. In the present theory only two equations are handled analytically against three by Swinerd and Boulton, and Boulton. The present solution is non-singular, whereas in Boulton’s expression of argument of perigee variation (Ao), e is in the denominator. EQUATIONS
OF MOTION
The K-S element equations of motion of a satellite under the effect of additional perturbing force F are (Stiefel and Scheifele, 1971) dw -=--&z*,LTP,, dE
(1) (2)
u’ = (u1,u2,u3,u4)
z=t+&ii*), W
=iEcos(f)+gsin(f),
u”
j;=(x1,x2,x3)=L(ii)ii, r=
B,
-
a3
P,
+ a2 ru*
and
L(G) =
B3 -
-u2
al -u3
u*
u, -uq
u3
u4
Lu4
-u3
Ul u2
P4
(4)
= 0, u4
-u3
1
I
(5)
u2 -%j
where Eq. (4), is the bilinear relation and E, w, t, r and K* are, respectively, the eccentric anomaly, angular frequency, physical time, radial distance and the gravitational constant. Knowing the position and velocity vector j; and 1 at the instant t = 0, the values of r, w, z, t& and h* can be computed and by adopting E = 0 as the initial value of the eccentric anomaly; we obtain
Jf 13 is the aerodynamic drag force per unit mass on a satellite of mass m, (King-Hele, D.G, 1987) F =[1-rp,,Acosi,lv,,]2, F = -ipSliJlG, S=FAC,fm,
(6)
Decay
of Satellite
Orbits
Using
K-S Elements
2013
where p is the atmospheric density, A is the rotational rate of the atmosphere about the Earth’s axis, rpo is the initial perigee radius, &,is the initial inclination, vpois the velocity at the initial perigee, Cu and A are, respectively, the drag coefficient and the effective area of the satellite. If the density, p is assumed to vary sinusoidally with 4, where 4 is the geocentric angular distance from the direction of the density maximum, F= LI&=”
P=Po(l+Fcosp)exp{-p(r-a)},
p, ,B=$, (7) 0 !lU” where p o is the average density on the reference spheroid when 9=90° and H,, is the average density scale height. Following (Cook et al, 1961) we choose o with 8, geocentric latitude, as (l-&sin* 0) CT=rpo (l-&sin’ QPO) ’ The ellipticity of each of the oblate spheroidal atmospheric surfaces is assumed to be the same as that of the Earth’s ellipticity, 0.00335. The density p, in Eq. (7) can be written as 1
p=k
{ 1 -I
F A(cosE-e)
(1-ecosE) x[l+ccos(2o+2E)
+ FB(l-e*)
:sinE
(1-ecosE)
~}xexp(zcosE)
+ce(cos(2o+3E)-cos(2w+E)}
+ce2{cos2w-4cos(20+2E)+3cos(20+4E)}/4 -cos(2@+5E)}/2+ -cos(4w+3E)}/2
-ce3{cos(2w+3E)
c2 {l+cos(4w+4E)}/4+c”
e{cos(40+5E)
+c2e2 {3cos(4w+2E)-8cos(4w+4E)
+5cos(4o+6E)}/8
+c3{3~~~(2~+2E)+~~~(6~+6E)}/24
-c3e{cos(2w+E)-cos(2u+3E)+cos(6w+5E)-cos(6w+7E)}/8 + c4 (3+4cos(4w+4E)+cos(8~+8E)}/192 k=Ppoexp{P(rP,-a)-ccos2w, A=sinS,
sinisinw
+cos~,{cos(S2-a,
+O(ce4,c2
e3,c3e2,c4 e,ce4,c5)1, 1 }, c=&$$,sin’ -i, 2 )cosw-cosisin(G-oc, )sinw},
(8)
B=sinS, sinicoso-cos6,{cos(SZ-ol,)sinw+cosisin(SZ-~:,)cosw], (a,, 6,) and (a~ , &,.,), where oM = a, +h and & = 6, are the right ascension and declination of the Sun and the centre of the day time bulge, i.e. the bulge lags behind in right ascension by an angle h. w and Q are argument of perigee and right ascension of ascending node of the orbit. In reality the scale height is known to vary with altitude, and a satellite will sample the atmosphere over a significant height range, even when the orbit is of low eccentricity. Similar to that of Cook and King-Hele (1963), we adopt an expression, H = Ho +p (r-o), (9) in which the scale height H varies linearly with height where p, the altitude gradient of H,, , is assumed to be constant in the range lpl< 0.2. Now, by definition, H and p are related through the equation - 1 = --- 1 dP(T) p(r)dr’ . H(r) Hence, if we substitute Eq. (9) in to Eq. (lo), the density is p(r)
=pl (l+~(r-a)l~, HO
(10)
(11)
2014
L. S. Nair and R. K. Sharma
Ignoring the terms of order pz and replacing p , by ‘p O (1 + F cos cp),the density equation p(r)= { l+p (r-o)‘/2H,” )p, (l+Fcos#)exp{-(r-a)/Ho }, 02) allows an approximately linear altitude variation in the density scale height. Using r = a (1 - e cosE) and z = ae / II,,, Eq. (12) can be expressed in terms of E as p(r) =~p{l++(l-
2cosE+cos2E)+O(~c)).
(13)
Following the terminology, notations and procedure as in Sharma (1991, 1992) s=-pc!TK2
[{fo+5e2(2fo+fl)/8+27e4(3fo+2f,)}+{e(2fo+fl)/64
+3e3 (3f, + 2f,)/4}cosE+(e+3
e3/4) f, sinE+{J
+27e4 (Sf, +7f, )/128}cos2E+(f, +3e3(3f,+2fO)/8}cos3E+
+5e2 /4+135e4
+5e2 (f, +fi )I4 /128)sin2E+{ef,
(e+9e3/8)f2sin3E+{5e2f,/8+27e4/64
+(f,+2f,))cos4E+(5e2/8+27e4/32)f2sin4E+3e3f,cos5E/8 3
+3e3 f2sin5E/8+27e4 ANALYTIC
fIcos6E/128+27e4
f2sin6E1128]/16a2
.
(14)
INTEGRATION
The explicit form of the perturbations Aw is obtained by substituting the value of p from Eq. (8) in to Eq. (14) and noting that the non-zero integrals in the resulting expression are of the form of linear combinations of the Bessel functions 2II
Z,(z) =&
j exp (z cos E) cos nEdE . 0
After some algebra, we obtain the change in w after one revolution as: Aw = ( C$/ a ) [A,,+ eA, +c{coQo A2 +sin2o A3)} +e2A4 +ce{cos2 o A5 + sin2 o A,} + c* {A, +cos4 o Aa +sin4 o AP} +e3Alo + ce2 (cos2 o A,, +sin2 o A12} + c2e { Al3 +cos40 Al4 + sin4o Ar5} +c3 {COS~OAt6+sin20 At7 + coda Ata + sin6e.rA19}+e4AU,+ ce3 { cos20 A2i+ sin2o AZ} +c2e2{ Aa +cos4u1Au +sin4o AZ} +c3e {COS~OA+ sin2o A+ coda Aa +sin6o A29 }+ c4 {Aso + COS~WA31 +sin4o A~~+COS~OI A33 +sin8w A34 }+ O(r) +FA [A35+ eA36 +c ( cos20 A3, +sin2o A3s)} +e2Ajs ice{ cos2k AN + sin2o A41} + c2 { h2 +cos4w A43 +sin4o Au} +e3AU + ce2 { cos2o A& +sin2o &T) + c2e {AM +cos4w A49 + sin40 Aso} +c3 (cosiw ASI +sin2o A52 + coda A& sin6o AN} +e4Ass+ ce3 { cos2w As6+ sin2w AZ7) +c2e2 {Ass +cos40 A59+sin4u) AW}+c3e { cos20 A61 +sin2o A,~~+cos~u A63}+sin60 A&c’{ A&cos40 AM +sin4o A67+cos80 Abs+sin8w A69 )I +FB [ATo+ eA,, +c{ cos20 AT2+sin2o AT3)} +e2AT4+ce{ cos2o AT5+ sin2o AX} + c2 (A,, +cos40 AT8+sin4w ATg} +e3Aso + ce2 (~0~20 Asr +sin2o Aaz}+ c2e (As3 +cos40 As4 + sin4w Aas} +c3 {costi Aw+sin2o Aa7 + cos60 Aas+ sin6w Asq} +e4Aw +ce3{cos2ur Ag,+ sin2w Ag2}+c2e2{Ag3+cos40 Ag4+sin4ccAg,}+c3e (~0~20 Ag6+sin2cuA97 +cos60 AgS}+sin6o Agg+c4{AlW+ cos40 Alol +sin4@ Alo*+ cos80 Aro3 +sin8 o A104 III , = (@la)A,
(13
2015
Decay of Satellite Orbits Using K-S Elements
where Ap ,q=O,l,....,, 104, are linear combinations of the modified Bessel functions IO(z),n = O,l,. . . ...16. The change in the state variables given by Eq. (3) after one revolution is A& =(@/w)h,
i=l.....
8,
(17)
where the coefficients A, in Eq. (16) change in Eq. (17). The changed coefficients are also linear combinations of I,(z), n =O,l,.., 16, and q~=-nkK6
exp(-pa)/80i.
When we consider the scale height variation, p is replaced by Eq. (13). Then the change in w and & are given by Aw, = Aw + (l/2) Joz2 N(Aw) +0(r). where v denotes the assumption of the variability in H, and N is an operator defined by N(1,) = I(O,n) - 2 I(l,n) + I(2,n). Here I(p,n) = 2’Z PC,In+pZr,r = 0, p; p = 0,1,2, where pC, are the usual binomial coefficients. By a similar process, we have A& = A& +( l/2) p Z2 N(A&) +0(r). NUMJZRICAL RESULTS We have programmed Eqs. (1) and (3) in double precision arithmetic to compute the K-S elements Q and pi. Then the position, velocity and the orbital elements can be easily computed. Four test caseswith perigee heights of 175, 240, 300 and 400 km, where the effects of the atmospheric oblateness and the diurnal bulge are significant, are chosen for detailed numerical study. At each perigee height two different inclinations of 5’ and 80’ are chosen to cover a wide range of inclination, and also to vary the parameter c, which is used in series expansions, from a very small value of I
Table 1. Decrease in semi-major axis (km) after 250 revolutions with percentage error orbit Parameters u = 0.2 I Method HP1 e I a I -I (km) NUM ( KSANAL4 [ % error ( ANAL4 1 % error 7279.78 93.9459 1 90.6000 1 3.56 1 89.1875 1 5.07 175 0.1 8183.50 95. 2295 I 94.1912 I 1.09 0.2 I 91.2925 I bTii 1 15.7828 -__- 1 4.46 ---240 0.1 7352.00 16.5201 1 16.0057 ( 3% 8271.00 0.2 0.47 16.0945 ~1 s.05 16.9510 1 16.8712
I
l
5
] 101.2190 1 34.94 1 106.6745 1 105.0527 1 80
3.93
1 99.6059
1.52
-_--
[ 5.5114 1 7525.50 1 0.9361
1 5.4770 0.9025
1
0.62
1 5.1622 I
0.8
1
.--5.46
I
]
1
2016
L. S. Nair and R. K. Sharma
0.002 to a.large value of 0.39, At each perigee height we have computed the change in semi-major axis, by which the energy of an orbit can be calculated, after 250 and 500 revolutions given in Tables 1 and 2, with reasonably high values of e (0.1, 0.2). e is a very important parameter in the present solution, as lot of series expansions are involved with it. In the entire test cases reported here, the values of O, &Z and mean anomaly are 60, 30 and 0 degrees, respectively. The density scale height variation with altitude is included and a high value of p = 0.2 is taken to test the analytical solution more rigorously. The effect of the flattening of the Earth has been included in computing the perigee and apogee heights. Tables 1 and 2 provide the change in semi-major axis, obtained by the present theory (KSANAU), Swinerd and Boulton method (ANAL4) and by numerical integration (NUM) of the K-S element Eqs. (1) and (3) with a fixed step size fourth order Runga-Kutta-Gill method having a very small step size of half degree in the eccentric anomaly by using the expressions for p given in Eq. (12). Jacchia (1977) atmospheric density model, which is relatively easier to use, is employed to compute the value of p. and I&, being used in the Eq. (13). Arbitrarily 22 August 2002 is chosen as the initial epoch. The values of solar flux (Fr0.r) and averaged geomagnetic index (A,) Table 2. Decrease in semi-major axis (km) after 500 revolutions with percentage error Orbit Parameters u=o.2
are taken as 150 and 10, respectively, which approximately represents an average density model and results in exospheric temperatures between 1000 and 1100 K, in the present study. The values of A, E and B, = (m /CnA) utilized during the computations are 1.2,0.0035 and 50.0, respectively. It can be noticed that for the epoch under consideration the values of FA and FB in Eq. (16) slightly decrease as the number of revolution increases. It is observed that these values increase as perigee height increases. Percentage error = 100 x (NUM - ANAL) / NUM is calculated with respect to the present solution as well as with Swinerd and Boulton’s theory. It is seen from Tables 1 and 2 that at lower altitudes small increase in percentage error is noticed with KSANAL4, when inclination increases from 5’ to SO’. Also, the percentage error decreases when e increases from 0.1 to 0.2. It is noticed from these tables that the drag perturbed orbital parameter, a, match quite well with NUM values with the present solution. On altitudes above 175 km, the percentage error seems to be less than 3.The maximum percentage error after 250 and 500 revolutions for the test cases under study, is less than 4.0 and 6.3, respectively, whereas with ANAL4 these values are 8.1 and 13.4. The accuracies of the numerical computations examined with the help of the bilinear relation in the K-S elements are found to be very satisfactory. The left-hand side of the bilinear equation, Eq. (4), gave the values of order lo”‘- lo-” with respect to KSANAL4 and NUM computations after 500 revolutions in the present study.
Decay of Satellite Orbits Using K-S Elements
2017
CONCLUSION The K-S element equations are integrated analytically by a series expansion method by assuming an oblate diurnal atmosphere with density scale height varying with altitude and a non-singular solution containing up to fourth-order terms in eccentricity and c (a small parameter dependent on the flattening of the atmosphere) is obtained. Only two of the nine equations are solved analytically to compute the state vector at the end of each revolution due to symmetry in the K-S element equations. Numerical comparison is made with numerically integrated values and with the modified solution of Swinerd and Boulton up to fourth order terms in e and c, with different perigee altitudes for reasonably high eccentric orbits, having small and high orbital inclination, after 250 and 500 revolutions. The present solution is found to provide better accuracies. ACKNOWLEDGMENTS We gratefully acknowledge the referees for their valuable comments and suggestions, which helped in bringing the paper to its present form. REFERENCES Barlier, F., C. Berger, J. L.. Falin, et al., A Thermospheric Model based on satellite drag data, Ann. Geophys., 34,9-24, 1978.
Boulton, W.J., The advance of a satellite in a rotating oblate atmosphere with a diurnal variation in density, Proc. R. Sot. Lord, A 389,433444,1983. Cook, G.E., D.G. King-Hele, and D.M.C. Walker, The contraction of satellite orbits under the influence of air drag. Part II. With oblate tmosphere, Proc. R.. Soc.Lond. A 264,88-121,1961. Cook, G.E., D.G. King-Hele, The contraction of satellite orbits under the influence of air drag. Part IV. With scale height dependent on altitude. Proc. R. Sot. Land. A 215,357-390,1963. Hedin, A.E., MSIS-86 Thermospheric model , J. Geophys. Res., 92(A5), 4649-4662, 1987. Hedin, A.E., Extension of the MSIS Thermosphere model into the middle and lower atmosphere, J. Geophys.Res.,96 (A2), 1159-l 172, 1991. Jacchia, L.G. Smithsonian Astrophysics Observatory Special report, No. 375, 1977. King-Heie, D.G. Satellite Orbits in an Atmosphere: Theory and Applications, Blackie, Glasgow and London, 1987. Santora, F. A., Satellite drag perturbations in an oblate diurnal atmosphere, AZAA J., 13, 1212-1216, 1975. Sharma, R.K., Analytical approach using K-S elements to near-Earth orbit predictions including drag. Proc. R. Sot. Lo&. A 433, 121-130,199l. Sharma, R.K., A third - order theory for the effect of drag on Earth satellite orbits, Proc. R.. Sot. Land. A 438, 467475,1992. Sharma, R.K., Contraction of satellite orbits using K-S elements in an oblate diurnally varying atmosphere. Proc. R. Sot. Lmd. A 453, 2353-2368, 1997. Stiefel, E.L., and G. Scheifele, Linear and Regular Celestial Mechanics, Springer - Veralg, Berlin, Heidelberg, New York, 1971. Swinerd, G.G., and W.J. Boulton, Contraction of satellite orbits in an oblate atmosphere with a diurnal density variation, Proc. R. Sot. Land. A 383,127-145, 1982. Swinerd, G.G., and W.J. Boulton, Near circular satellite orbits in an oblate diurnally varying atmosphere, Proc. R. Sot. Land. A 389, 153-170, 1983. Wolfram, S, Mathematics Version 3, Wolfram Media, Cambridge University Press, USA, 1996. E-mail address of Lila S Nair [email protected] Manuscript received 7 November 2002; revised 9 March 2003, accepted 14 March 2003