HEOSAT: A mean elements orbit propagator program for Highly Elliptical Orbits
CCEAS Space Journal manuscript No. (will be inserted by the editor)
HEOSAT: A mean elements orbit propagator program for HighlyElliptical Orbits
Martin Lara · Juan F. San-Juan · Denis Hautesserres
CEAS Space Journal (ISSN: 1868-2502, ESSN: 1868-2510) (2018) doi:10.1007/s12567-017-0152-x (Pre-print version)
Abstract
The algorithms used in the construction of a semi-analytical propagator for the long-term propagation of High-ly Elliptical Orbits (HEO) are described. The software prop-agates mean elements and include the main gravitational andnon-gravitational effects that may affect common HEO or-bits, as, for instance, geostationary transfer orbits or Mol-niya orbits. Comparisons with numerical integration showthat it provides good results even in extreme orbital config-urations, as the case of SymbolX.
Keywords
HEO · Geopotential · third-body perturbation · tesseral resonances · SRP · atmospheric drag · meanelements · semi-analytic propagation The subject of analytical or semi-analytical propagation isvery old. Since the first analytical orbit propagators based on
Funded by CNES contract Ref. DAJ-AR-EO-2015-8181. Preliminaryresults were presented at the 6th ICATT, Darmstadt, Germany, March14 – 17, 2016M. LaraGRUCACI – University of La RiojaC/ Madre de Dios 53, Edificio CCT,26006 Logro˜no, SpainTel.: +34-941-299440Fax: +34-941-299460E-mail: [email protected]. San-JuanGRUCACI – University of La RiojaE-mail: [email protected]. HautesserresCCT – Centre National d’ ´Etudes Spatiales18 Av. Edouard Belin,31401 Toulouse CEDEX 4, FranceE-mail: [email protected] intermediary solutions to the J problem [34,11], the contin-uous increase in the accuracy of observations demanded theuse of more complex dynamical models to achieve a sim-ilar precision in the orbit predictions. In particular higherdegrees in the Legendre polynomials expansion of the third-body disturbing function are commonly required (see [17,24], for instance). Useful analytical theories needed to dealwith a growing number of effects, a fact that made that thetrigonometric series evaluated by the theory comprised tensof thousands of terms [5].In an epoch of computational plenty, the vast possibili-ties offered by special perturbation methods clearly surpassthose of general perturbation methods in their traditionalapplication to orbit propagation. Apparently by this reasonanalytical perturbations have these days been cornered to adowngraded role of providing some insight into the problemunder investigation, a task for which a first order averag-ing is usually considered to be enough, yet the computationsof higher orders may provide important details on the dy-namics [33,21]. However, analytical theories like the pop-ular SGP4 [16] still enjoy a wide number of users mainlyinvolved in catalog propagation duties, a role in which othertools like the Draper Semi-Analytic Satellite Theory [29,6]or the numeric-analytic theory THEONA [12] can competeto numerical integration up to a limited accuracy. On the other hand, new needs in satellite propagation,like the challenges derived of compliance with Space Law,motivate the development of software tools based on analyti-cal or semi-analytical methods, as, for instance, STELA . Inaddition, design of end of life disposal strategies may requirethe long term propagation of thousands of trajectories to findan optimal solution; accurate ephemeris are not needed in A partial list of orbit propagators can be found in http://facul-ty.nps.edu/bneta/papers/list.pdf, accessed September 29, 2016. https://logiciels.cnes.fr/content/stela a r X i v : . [ a s t r o - ph . E P ] J u l Martin Lara et al. the preliminary design and using semi-analytic propagationmakes the approach quite feasible [1].Current needs for long-term propagation at the CentreNational d’ ´Etudes Spatiales motivate the present research.HEOSAT, a semi-analytical orbit propagator to study thelong-term evolution of spacecraft in Highly Elliptical Orbits(HEO) is presented. The perturbation model used includesthe gravitational effects produced by the more relevant zonalharmonics as well as the main tesseral harmonics affectingto the 2:1 resonance of earth’s gravitational potential, whichhas an impact on Molniya-type orbits; the third body per-turbations in the mass-point approximation, which only in-clude the Legendre polynomial of second order for the sunand the polynomials from second order to sixth order in thecase of the moon; solar radiation pressure, in the cannonballapproximation, and atmospheric drag.The forces of gravitational origin are modeled taking ad-vantage of the Hamiltonian formalism. Besides, the prob-lem is formulated in the extended phase space in order toavoid time-dependence issues. The solar radiation pressureand the atmospheric drag are added as generalized forces.The semi-analytical theory is developed using perturbationtechniques based on Lie transforms. Deprit’s perturbationalgorithm [7] is applied up to the second order of the secondzonal harmonics, J . In order to avoid as far as possible thelost of long-period effects from the mean elements Hamilto-nian, the theory is corrected by the inclusion of long-periodterms of the Kozai-type [20,10]. The transformation is de-veloped in closed-form of the eccentricity except for tesseralresonances, and the coupling between J and the moon’s dis-turbing effects are neglected.This paper describes the semi-analytical theory and pres-ents relevant examples of the numerical validation. An ex-tensive description of the tests performed in the validationof the HEOSAT software is given in [22]. Satellites in earth’s orbits are affected by a variety of per-turbations of different nature. A full account of this can befound in textbooks on orbital mechanics (see, for instance,[32,35]). All known perturbations must be taken into ac-count in orbit determination problems. However, for orbitprediction the accuracy requirements are notably relaxed,and hence some of the disturbing effects may be consideredof higher order in the perturbation model.Furthermore, for the purpose of long-term predictions itis customary to ignore short-period effects, which occur ontime-scales comparable to the orbital period. Thus, in thecase of the gravitational potential the focus is on the ef-fect of even-degree zonal harmonics, which are known tocause secular effects. Odd-degree zonal harmonics may also be important because they give rise to long-period effects,whereas the effects of tesseral harmonics in general aver-age out to zero. The latter, however, can have an importantimpact on resonant orbits, as in the case of geostationarysatellites (1 to 1 resonance) or GPS and Molniya orbits (2 to1 resonance).The importance of each perturbation acting on an earth’ssatellite depends on the orbit characteristics, and fundamen-tally on the altitude of the satellite. Thus, for instance theatmospheric drag, which can have an important impact inthe lower orbits, may be taken as a higher order effect foraltitudes above, say, 800 km over the earth’s surface, and isalmost negligible above 2000 km.A sketch of the order of different perturbations whencompared to the Keplerian attraction is presented in Fig. 1based on approximate formulas borrowed from [32, p. 114].As illustrated in the figure, the non-centralities of the Geopo-tential have the most important effect in those parts of the or-bit that are below the geosynchronous distance, where the J contribution is a first order perturbation and other harmonicscause second order effects. To the contrary, in those parts ofthe orbit that are farther than the geosynchronous distancethe gravitational pull of the moon is the most important per-turbation, whereas that of the sun is of second order whencompared to the disturbing effect of the moon, and pertur-bations due to J and solar radiation pressure (SRP) are ofthird order. J SRP (cid:72) A (cid:144) M (cid:61) m (cid:144) kg (cid:76) Moon Sun J J (cid:45) (cid:45) (cid:45) (cid:72) thousands of km (cid:76) O r d e r o f m a gn it ud e Fig. 1
Perturbation order relative to the Keplerian attraction (after [1]).An area to mass ratio of 0 . / kg was taken for estimating the SRP.The horizontal, dashed line corresponds to 3 . × − times the Kep-lerian attraction. Finally, we recall that integration of osculating elementsis properly done only in the True of Date system [9]. For thisreason Chapront’s solar and lunar ephemeris are used [4,3],which are directly referred to the mean of date in this wayincluding the effect of equinoctial precession.In the case of a highly elliptical orbit the distance fromthe satellite to the earth’s center of mass varies notably alongthe orbit, a fact that makes particularly difficult to establish
EOSAT: A mean elements orbit propagator program for Highly Elliptical Orbits 3 the main perturbation over which to hinge the correct pertur-bation arrangement. This issue is aggravated by the impor-tance of the gravitational pull of the moon on high altitudeorbits, which requires taking higher degrees of the Legendrepolynomial expansion of the third-body disturbing functioninto account. Because this expansion converges slowly, thesize of the multivariate Fourier series representing the moonperturbation will soon grow enormously. Besides, for oper-ational reasons, different HEOs may be synchronized withthe earth rotation, and hence, being affected by resonanceeffects. Also, a HEO satellite will spend most of the time inthe apogee region, where, depending on the physical charac-teristics of the satellite, the solar radiation pressure may havean observable effect in the long-term. Finally, the perigee ofusual geostationary transfer orbits (GTO) will enter repeat-edly the atmosphere, yet only for short periods.Therefore, a long-term orbit propagator for HEO aimingat describing at least qualitatively the orbit evolution mustconsider the following perturbations: – the effects of the main zonal harmonics of the Geopoten-tial ( J — J in our case), as well of second order effectsof J ; – the effects of the tesseral harmonics of the Geopotentialthat affect the most important resonances, and, in partic-ular, the 2:1 resonance that impacts Molniya orbits – lunisolar perturbations (mass-point approximation). Thefirst few terms of the Legendre polynomials expansionof the third-body perturbation are needed in the case ofthe moon, whereas the effect of the polynomial of thesecond degree is enough for the sun; – solar radiation pressure effects; – and the circularizing effects of atmospheric drag affect-ing the orbit semi-major axis and eccentricity.In what respects to the ephemeris of the sun and moon,needed for evaluating the corresponding disturbing poten-tials, the precision provided by simplified analytical formu-las is enough for the accuracy requirements of a perturbationtheory. For this reason, truncated series from Chapront’s so-lar and lunar ephemeris are taken from [31]. The precisionof these short truncations is enough for mean elements prop-agation and notably speed the computations. The aim is to find a canonical transformation ( x , X ) → ( x (cid:48) , X (cid:48) ; ε ) , where x ∈ R m are coordinates, X ∈ R m their conjugate mo-menta, and ε is a small parameter, that converts the Hamil-tonian H = ∑ m ≥ ε m m ! H m , (1) which is written in terms of the x and X , into a new Hamil-tonian (cid:102) H = ∑ m ≥ ε m m ! H , m , (2)in the new, prime variables, which is free from short-periodterms.The transformation is derived from a generating function W = ∑ m ≥ ε m m ! W m + , (3)each term of which is computed as a solution of the homo-logical equation L ( W m ) + (cid:102) H , m = H , m , (4)as follows: – terms (cid:102) H , m are known. They are obtained from previouscomputations based on Deprit’s recurrence H n , q + = H n + , q + n ∑ m ≥ (cid:18) nm (cid:19) { H n − m , q ; W m + } , (5)where { ; } means the Poisson bracket operator. – the new term H , m is chosen as the part of (cid:102) H , m that isfree from short-period terms – the term W m is then solved from Eq. (4), where L ( W m ) = { H , ; W m } + ∂ W m ∂ t , (6)is customarily known as the Lie derivative of W m .After solving the homological equation up to the desiredorder, the new Hamiltonian in Eq. (2) is obtained by chang-ing the old variables in which the terms H , m are written bythe new ones.Once the generating function is known, the transforma-tion from old to new variables is itself computed by applyingDeprit’s recurrence in Eq. (5) to the vector functions x = ∑ m ≥ ε m m ! x m , ( x , X ) , X = ∑ m ≥ ε m m ! X m , ( x , X ) , where x , = x , x m , = m > X , = X , X m , = m > The Hamiltonian Eq. (1) is arranged as follows H , = − µ a , H , = − µ r R ⊕ r C , P ( sin ϕ ) , H , = ( Z + T + V ( | + V (cid:12) ) , H m , = , m ≥ , Martin Lara et al. where µ is the earth’s gravitational parameter, ( r , ϕ , λ ) arethe spherical coordinates of the satellite and a is the semi-major axis of the satellite’s orbit, R ⊕ is the earth’s equatorialradius, P m are the usual Legendre polynomials, and C m , arezonal harmonic coefficients.4.1 Geopotential second order perturbationsThe second order Hamiltonian term H , comprises the sec-ond order perturbations of gravitational origin. Thus, thehigher order terms of the zonal potential are Z = − µ r ∑ m ≥ R m ⊕ r m C m , P m ( sin ϕ ) . The tesseral disturbing function is T = − µ r ∑ m ≥ R m ⊕ r m (7) × m ∑ n = [ C m , n cos n λ + S m , n sin n λ ] P m , n ( sin ϕ ) where P m , n are associated Legendre polynomials, and C n , m and S n , m , m (cid:54) =
0, are the non-zonal harmonic coefficients.Besides, for orbital applications, instead of using sphericalcoordinates the satellite’s radius vector is expressed in or-bital elements. This can be done by first replacing the spher-ical coordinates by Cartesian onessin ϕ = zr , sin λ = y (cid:112) x + y , cos λ = x (cid:112) x + y . Then, the orbital and inertial frame are related by means ofsimple rotations, to give xyz = R ( − Ω ) R ( − I ) R ( − θ ) r (8)where R and R are usual rotation matrices, θ = f + ω , r = p + e cos f , ( a , e , I , Ω , ω , M ) are traditional orbital elements and p = a ( − e ) , is the semi-latus rectum of the osculating ellipse.4.2 Third-body perturbationsUnder the assumption of point masses, the third-body dis-turbing potential is V (cid:63) = − µ (cid:63) r (cid:63) (cid:18) r (cid:63) || r − r (cid:63) || − r · r (cid:63) r (cid:63) (cid:19) (9)where V (cid:63) ≡ V ( | for the moon, and V (cid:63) ≡ V (cid:12) in the case of thesun. If the disturbing body is far away from the perturbed body, which will always be the case when dealing with per-turbed Keplerian motion, Eq. (9) can be expanded in powerseries of the ratio r / r (cid:63) V (cid:63) = − β n (cid:63) a (cid:63) r (cid:63) ∑ m ≥ r m r m (cid:63) P m ( cos ψ (cid:63) ) , (10)where β = m (cid:63) m (cid:63) + m , (11)andcos ψ (cid:63) = xx (cid:63) + yy (cid:63) + zz (cid:63) rr (cid:63) . (12)The semi-analytic theory only considers P in Eq. (10)for the sun disturbing potential, whereas P – P are taken forthe moon. Both solar and lunar ephemerides are taken fromthe low precision formulas given by [31] (see also [30]).The Lie transforms averaging is only developed up tothe second order in the small parameter, so there is no cou-pling between the different terms of the disturbing function.Hence, the generating function of the averaged Hamiltoniancan be split into different terms which are simply added atthe end.4.3 Time dependencyOrbits with large semi-major axis will experience high per-turbations due to the moon’s gravitational attraction, andmay be better described in the restricted three-body prob-lem model approximation than as a perturbed Keplerian or-bit. However, since our approach is based on perturbed Ke-plerian motion, we still take the Keplerian term as the zeroorder Hamiltonian.Because sun and moon ephemeris are known functionsof time, the perturbed problem remains of three degrees offreedom, but the Lie derivative in Eq. (6) must take the timedependency into account, viz. L ( W ) ≡ { H , ; W } + ∂ W ∂ t = n ∂ W ∂ M + ∂ W ∂ t , where n stands for mean motion.Dealing explicitly with time can be avoided by movingto the extended phase space. Then, assuming that the semi-major axis, eccentricity, and inclination of the third-body or-bits, remain constant L ( W ) = n ∂ W ∂ M + n ( | ∂ W ∂ M ( | + ∂ W ∂ ω ( | d ω ( | d t + ∂ W ∂ Ω ( | d Ω ( | d t + n (cid:12) ∂ W ∂ M (cid:12) + ∂ W ∂ ω (cid:12) d ω (cid:12) d t + ∂ W ∂ Ω (cid:12) d Ω (cid:12) d t , which, in view of the period of the lunar perigee motion of8.85 years (direct motion) and the period of the lunar nodemotion of 18.6 years (retrograde motion), and the higher pe-riods in the case of the sun relative orbit, can be safely ap-proximated by L ( W ) ≈ n ∂ W ∂ M + n ( | ∂ W ∂ M ( | + n (cid:12) ∂ W ∂ M (cid:12) . (13) EOSAT: A mean elements orbit propagator program for Highly Elliptical Orbits 5 W . In view of the form of theLie derivative in Eq. (13), the directions of both the sun andthe moon must be expressed as functions of the sun andmoon mean anomaly, respectively. In the case of the sun,we use x (cid:12) y (cid:12) z (cid:12) = R ( − ε ) R ( − λ (cid:12) ) R ( β (cid:12) ) r (cid:12) , (14)where ε is the mean obliquity of the ecliptic, r (cid:12) is the ra-dius of the sun, and β (cid:12) and λ (cid:12) are the ecliptic latitude andlongitude of the sun, respectively. The sun’s latitude can beneglected because it never exceeds 1 . x ( | y ( | z ( | = R ( − ε ) R ( − Ω ( | ) R ( − J ) R ( − θ ( | ) r ( | , where θ ( | is the argument of the latitude of the moon (whichis also expressed as a function of the moon’s mean anomalyby standard formulas), J ≈ . ◦
15 is the inclination of themoon orbit over the ecliptic, which is affected of periodicoscillations whose period slightly shorter than half a yearbecause the retrograde motion of the moon’s line of nodes,and Ω ( | is the longitude of the ascending node of the moonorbit with respect to the ecliptic measured from the meanequinox of date.Note, however, that in the present stage, the semi-analyt-ical theory only deals with mean elements and correspond-ing homological equations of the type in Eq. (13) do notneed to be solved at the second order. Future evolutions ofthe theory will take the short-period corrections due to third-body perturbations into account. The zonal part of the Hamiltonian is first transformed. Be-cause of the actual values of the zonal coefficients of theearth, the old Hamiltonian is arranged in the form H , = − µ a (15) H , = − µ r R ⊕ r C , P ( s sin θ ) (16) H , = − µ r ∑ m ≥ R m ⊕ r m C m , P m ( s sin θ ) (17) where we abbreviate s ≡ sin I , and all the symbols, viz. a , r , I , and θ , are assumed to be functions of some set of canoni-cal variables. In particular, the averaging is carried out basedon the canonical set of Delaunay variables, which is madeof the coordinates (cid:96) = M , g = ω , h = Ω , and their conjugatemomenta L = √ µ a , G = L η , H = G cos I , respectively, where η = (cid:112) − e , is customarily known as the “eccentricity function”. The De-launay variables are action-angle variables for the Keplerianmotion, and fit particularly well in the construction of per-turbation theories for perturbed Keplerian motion.The homological equation of the zonal problem is ob-tained from Eq. (4) as − n ∂ W m ∂ (cid:96) + (cid:102) H , m = H , m , (18)In consequence, the terms of the generating function arecomputed from quadratures of the form W m = n (cid:90) ( (cid:102) H , m − H , m ) d (cid:96). (19)5.1 Elimination of the parallaxIt is known that the removal of short-period terms from theHamiltonian is facilitated to a large extent by a preprocess-ing of the original Hamiltonian in order to remove paral-lactic terms [8]. That is, by first applying the “parallacticidentity”1 r m = r r m − = r ( + e cos f ) m − p m − , m > , (20)and next selecting H , m by removing all the trigonometricterms that explicitly contain the true anomaly in the expan-sion of H m , as a Fourier series [26,25]. At the first step of Deprit’s recurrence in Eq. (5), the knownterms are simply (cid:102) H , ≡ H , , where H , = C , (cid:110) ( − s )( + e cos f ) + s (cid:104) e cos ( f + ω )+ ( f + ω ) + e cos ( f + ω ) (cid:105)(cid:111) R ⊕ r n p η Then, the new Hamiltonian term H , is selected as, H , = µ p C , R ⊕ p (cid:18) − s (cid:19) p r (21) Martin Lara et al. and the generating function term W must be computed fromEq. (19) with m =
1. That is, W = n (cid:90) C , R ⊕ r n p η (cid:104) (cid:0) − s (cid:1) e cos f (22) + s ∑ j = jE , j cos ( j f + ω ) (cid:105) d (cid:96), in which E , = e , E , = , E , = e . (23)Equation (22) is solved in closed form by recalling thedifferential relationd M = ( r / p ) η d f , (24)based on the preservation of the angular momentum of theKeplerian motion. We get W = k + nR ⊕ C , η (cid:104) ( − s ) e sin f + s ∑ j = E , j sin ( j f + ω ) (cid:105) , where k is an arbitrary function independent of (cid:96) .To avoid the appearance of hidden long-period terms in W , we choose k in such a way that12 π (cid:90) π W d M = . In this way it is guaranteed that W only comprises short-period terms. That is k = − π nR ⊕ C , η (cid:90) π (cid:104) ( − s ) e sin f + s ∑ j = E , j sin ( j f + ω ) (cid:105) d M . Using, again, Eq. (24) to compute this quadrature, we get k = nR ⊕ C , η − η + η ( + η ) s sin 2 ω . Therefore, calling E , = ( + η ) − η + η = + η ( + η ) e , the first term of the generating function of the elimination ofthe parallax simplification is written W = nR ⊕ C , η (cid:104) ( − s ) e sin f + s ∑ j = E , j sin ( j f + ω ) (cid:105) which is now free from hidden long-period terms. Deprit’s recurrence now gives H , = { H , , W } + H , , H , = { H , , W } + { H , , W } + H , . Hence, the computable terms of the homological equation(18) are (cid:102) H , = { H , , W } + { H , , W } + H , . (25)The computation of the Poisson brackets is straightforwardeven for the partial derivatives of the true anomaly, whichcannot be explicitly written in terms of the Delaunay vari-ables. For the reader’s convenience, a table of partial deriva-tives is given in Appendix A, where all of them are ex-pressed as functions of classical orbital elements and relatedfunctions rather than in Delaunay variables; this way of pro-ceeding eases computations in the whole process by avoid-ing the appearance of square roots, on the one hand, andretains the engineering insight provided by the orbital ele-ments, on the other.After carrying out the required operations, parallacticterms are removed from (cid:102) H , using Eq. (20). Then, (cid:102) H , is written in the form of a Poisson series, from which H , is chosen by removing the trigonometric terms of (cid:102) H , thatexplicitly depend on the true anomaly, to give: H , = µ p p r ∑ i = J ∗ i − µ p p r C , R ⊕ p (cid:26) − s (26) + s + (cid:18) c − s (cid:19) e − e s cos 2 ω × (cid:20) − s − (cid:0) − s (cid:1) η ( + η ) (cid:21)(cid:27) with c ≡ cos I and J ∗ i = C i , R i ⊕ p i (cid:98) i / (cid:99)− ∑ j = e l Q i , l s l B i , l ˙ıı m exp ( ˙ıı l ω ) (27)where l = j + m , (cid:98) (cid:99) notes an integer division, m ≡ i mod 2,and ˙ıı = ( − ) / . Note that Eq. (27) applies also to i = H , = µ p p r J ∗ . (28)The eccentricity polynomials Q i , j in Eq. (27) are given inTable 1 of Appendix B, while the inclination polynomials B i , j are given in Tables 2–3.Finally, the simplified Hamiltonian is obtained by re-placing the old variables by the new ones in all the H , m terms. That is, the new Hamiltonian is obtained by assum-ing that all symbols that appear in Eqs. (21) and (26) arefunctions of the new (prime) Delaunay variables. EOSAT: A mean elements orbit propagator program for Highly Elliptical Orbits 7 H (cid:48) = ∑ ( ε m / m ! ) H (cid:48) m , where H (cid:48) , is the same as H , , H (cid:48) , is the same as H , in Eq. (21), and, H (cid:48) , isthe same as H , in Eq. (26), but all of them expressed inthe new, prime variables.The new Hamiltonian term H (cid:48) , is chosen by removingthe short-period terms from H (cid:48) , , namely H (cid:48) , = π (cid:90) π H (cid:48) , d M = µ p C , R ⊕ p η (cid:18) − s (cid:19) (29)which was trivially solved using Eq. (24).The first term of the new generating function is solvedby quadrature from the homological equation, from which W (cid:48) = n (cid:20) − H (cid:48) , (cid:96) + (cid:90) H (cid:48) , r p η d f (cid:21) = φ n H (cid:48) , , (30)where φ = f − M , is the equation of the center. Since φ is made only of short-period terms, there is no need of introducing additional in-tegration constants in the solution of W (cid:48) .Analogously to Eq. (25), the computable terms of thesecond order homological equation are (cid:102) H (cid:48) , = { H (cid:48) , , W (cid:48) } + { H (cid:48) , , W (cid:48) } + H (cid:48) , , from which expression the new Hamiltonian term H (cid:48) , ischosen by removing the short-period terms.After performing the required operations and replacingprime variables by new, double prime variables in all the H (cid:48) , m terms, we get the averaged Hamiltonian H (cid:48)(cid:48) = H (cid:48) , + H (cid:48) , + ( / ) H (cid:48) , , viz. H (cid:48)(cid:48) = − µ a + µ p η ∑ i = J ∗ i + µ p η C , R ⊕ p (cid:40) c (31) × ( − c ) − (cid:18) + s − s (cid:19) e − η ( − c ) − (cid:20) ( − c ) − ( − c ) η ( + η ) (cid:21) e s cos 2 ω (cid:27) . and the orbit evolution is obtained from Hamilton equationsd ( (cid:96) (cid:48)(cid:48) , g (cid:48)(cid:48) , h (cid:48)(cid:48) ) d t = ∂ H (cid:48)(cid:48) ∂ ( L (cid:48)(cid:48) , G (cid:48)(cid:48) , H (cid:48)(cid:48) ) , (32)d ( L (cid:48)(cid:48) , G (cid:48)(cid:48) , H (cid:48)(cid:48) ) d t = − ∂ H (cid:48)(cid:48) ∂ ( (cid:96) (cid:48)(cid:48) , g (cid:48)(cid:48) , h (cid:48)(cid:48) ) , (33) Now, the perturbation Hamiltonian is arranged as H , = − µ a , H , = , H , = ( V ( | + V (cid:12) ) where the sun and moon potentials are computed from theLegendre series expansion in Eq. (10).Because the maximum power of χ in P m ( χ ) is χ m , inview of Eqs. (10) and (12), we check that the satellite’s ra-dius r appears now in numerators, contrary to the Geopo-tential case where r always appear in denominators. For thatreason, the closed form theory for third-body perturbationsis approached using the eccentric anomaly u instead of thetrue one f . Hence, r = a ( − e cos u ) , and the Cartesian coordinates of the satellite are expressedin terms of the eccentric anomaly replacing in Eq. (8) theknown relations from the ellipse geometry r sin f = a η sin u , r cos f = a ( cos u − e ) . (34)In view of H , ≡
0, we choose H , = W =
0. Then, the second order of the homological equation(4) is L ( W ) + (cid:102) H , = H , , where (cid:102) H , = H , from De-prit’s recurrence (5). The new Hamiltonian term H , is cho-sen by removing the short-period terms from H , . Again,this is done in closed-form by computing the average H , = π (cid:90) π H , ra d u , (35)where we used the differential relationd M = ( − e cos u ) d u = ( r / a ) d u (36)which is obtained from Kepler equation.After computing Eq. (35) we obtain the long-term Ham-iltonian H , = − µ / ( a ) , H , =
0, and H , = ( na ) β ∗ a (cid:63) r (cid:63) n (cid:63) n ∑ m ≥ a m − r m − (cid:63) Γ m , (37)with the non-dimensional coefficients Γ m = (cid:98) m / (cid:99) ∑ j = A m , j m ∑ l = − m P m , j , l ( S (cid:63) m , l cos α + T (cid:63) m , l sin α ) , (38)where α = ( j + k ) ω + l Ω , and k = m mod 2, cf. [27].The eccentricity coefficients A m , j ≡ A m , j ( e ) , the inclina-tion ones P m , j , l ≡ P m , j , l ( I ) , and the third-body direction coef-ficients T (cid:63) m , l ≡ T (cid:63) m , l ( u (cid:63) , v (cid:63) , w (cid:63) ) , S (cid:63) m , l ≡ S (cid:63) m , l ( u (cid:63) , v (cid:63) , w (cid:63) ) , where u (cid:63) = x (cid:63) r (cid:63) , v (cid:63) = y (cid:63) r (cid:63) , w (cid:63) = z (cid:63) r (cid:63) , are given in Tables 4, 5–7, and 8 of Appendix B, respec-tively. They are valid for both the moon ( (cid:63) ≡ ( | ) and thesun ( (cid:63) ≡ (cid:12) ) by using the proper third-body direction vec-tor ( u (cid:63) , v (cid:63) , w (cid:63) ) .The contribution of lunisolar perturbations to the meanelements equations is by adding to Eqs. (32)–(33) the termsof Eq. (37) derived from corresponding Hamilton equations. Martin Lara et al.
The tesseral potential is no longer symmetric with respectto the earth’s rotation axis. Therefore, longitude dependentterms will explicitly depend on time when referred to theinertial frame.To avoid the explicit appearance of time in the Hamilto-nian, we move to a rotating frame with the same frequencyas the earth’s rotation rate n ⊕ . The argument of the node inthe rotating frame is h = Ω − n ⊕ t , and, in order to preserve the symplectic character, we furtherintroduce the Coriolis term − n ⊕ H into the Hamiltonian. Itis then simple to check that H = Θ cos I still remains as theconjugate momentum to h .Then, the tesseral Hamiltonian is arranged as a perturba-tion problem in which H , = − µ a − n ⊕ Θ cos I , H , = , H , = T , where, now, H , is the Keplerian in the rotating frame, andthe tesseral potential is given in Eq. (7). Now, the Lie deriva-tive in Eq. (6) reads { H , ; W } = − n ∂ W ∂ (cid:96) + n ⊕ ∂ W ∂ h , and the solution of the homological equation (4) will intro-duce denominators of the type ( in − jn ⊕ ) , with i and j inte-gers. Therefore, resonances n / n ⊕ = j / i between the rotationrate of the node in the rotating frame and the mean motionof the satellite introduce the problem of small divisors.In fact, resonant tesseral terms introduce long-period ef-fects in the semi-major axis that may be not negligible evenat the limited precision of a long-term propagation. There-fore, these terms must remain in the long-term Hamiltonian.Furthermore, these terms must be traced directly in the meananomaly, contrary to true anomaly, to avoid leaving short-period terms in the Hamiltonian, which will destroy the per-formance of the semi-analytical integration. Hence, trigono-metric functions of the true anomaly must be expanded asFourier series in the mean anomaly whose coefficients are(truncated) power series in the eccentricity.After the short-period terms have been removed from thetesseral Hamiltonian, we return to the inertial frame by drop-ping the Coriolis term and replacing h by the right ascensionof the ascending node (RAAN), in this way the time explic-itly appears into resonant terms of the long-period Hamilto-nian.From Kaula expansions [18], we find that the main termsof the Geopotential that are affected by the 2:1 tesseral res-onance are R = − µ a R ⊕ a (cid:110) F , , G , , − × [ C , cos ( α + ω ) + S , sin ( α + ω )]+ F , , G , , ( C , cos α + S , sin α ) + F , , × G , , [ C , cos ( α + ω ) + S , sin ( α − ω )] (cid:111) , in which α = ( Ω − n ⊕ t ) + M , is the (slowly evolving) resonant angle, F , , = ( + c ) , F , , = s , F , , = ( − c ) (39)and, up to O ( e ) G , , − = − e + e − e − e − e − e − e − e G , , = e + e + e + e + e + e + e + e G , , = e + e + e + e + e + e + e Other 2:1-resonant terms can be found in [23].For the 1:1 tesseral resonance, we find R = − µ a R ⊕ a (cid:110) F , , G , , (cid:104) C , cos ( α + ω )+ S , sin ( α + ω ) (cid:105) + F , , G , , × [ C , cos 2 α + S , sin 2 α ] (cid:111) − µ a R ⊕ a (cid:110) F , , G , , − (cid:104) C , sin ( α + ω ) − S , cos ( α + ω ) (cid:105) + F , , G , , × ( C , sin α − S , cos α ) + F , , G , , × [ C , sin ( α − ω ) − S , cos ( α − ω )] (cid:111) , where, now, α = Ω − n ⊕ t + M . In the particular case of the earth, C , = O ( − ) and S , = O ( − ) . Due to the smallness of these values, corre-sponding terms are commonly neglected from the resonanttesseral potential R . Therefore, the only needed inclina-tion polynomials are F , , and F , , , which were alreadygiven in Eq. (39), whereas the required eccentricity func-tions, up to O ( e ) , are G , , = − e + e − e − e − e − e − e − e G , , = e + e + e + e + e + e + e + e Terms of the Hamilton equations derived from the dis-turbing functions R or R , will be added to the evolu-tion equations (32)–(33), only for the propagation of thoseorbits which experience the corresponding resonance. EOSAT: A mean elements orbit propagator program for Highly Elliptical Orbits 9
The evolution equations must be completed by adding to theright hand side of Hamilton equations, Eqs. (32)–(33), theaveraged effects of the generalized forces. Because these ef-fects are derived from Gauss equations, we recall thatd L d t = na d a d t d G d t = d L d t η − na e η d e d t d H d t = d G d t c − na η s d I d t α srp = − F srp i (cid:12) , is always in the opposite direction of the unit vector of thesun i (cid:12) . If, besides, it is assumed that [19], – the parallax of the sun is negligible – the solar flux is constant along the satellite’s orbit – there is no re-radiation from the earth’s surfacethe magnitude of the SRP acceleration is F srp = ( + β ) P (cid:12) a (cid:12) r (cid:12) Am , where β is the index of reflection (0 < β < A / m is thearea-to-mass ratio of the spacecraft, a (cid:12) is the semi-majoraxis of the sun’s orbit around earth, r (cid:12) is the radius of thesun’s orbit around earth, and the solar radiation pressureconstant at one AU is P (cid:12) ≈ . × − N / m [32, p. 77],The components of i (cid:12) in the radial, tangent, and normaldirections, respectively, are obtained by simple rotations i (cid:12) = R ( θ ) R ( I ) R ( Ω ) R ( − ε ) R ( − λ (cid:12) ) . Then, calling F = − F srp / µ Kozai’s analytical expres-sions for perturbations due to SRP [19] are easily recoveredfrom the usual Gauss equations. After averaging over themean anomaly, which is done in closed form based on thedifferential relation (36), we getd a d t = e d t = na (cid:110) sin ω (cid:104) ( cos ε − ) cos ( λ (cid:12) + Ω ) (41) − ( cos ε + ) cos ( λ (cid:12) − Ω ) (cid:105) + cos ω (cid:104) s sin ε × sin λ (cid:12) + c ( cos ε + ) sin ( λ (cid:12) − Ω )+ c ( cos ε − ) sin ( λ (cid:12) + Ω ) (cid:105)(cid:111) η F d I d t = na e η F cos ω (cid:104) s ( cos ε + ) sin ( λ (cid:12) − Ω ) (42) − c sin ε sin λ (cid:12) + s ( cos ε − ) sin ( λ (cid:12) + Ω ) (cid:105) d Ω d t = na e η s F sin ω (cid:104) s ( cos ε + ) sin ( λ (cid:12) − Ω ) (43) − c sin ε sin λ (cid:12) + s ( cos ε − ) sin ( λ (cid:12) + Ω ) (cid:105) d ω d t = − na Fe η (cid:40) sin ω (cid:104) ( cos ε + ) sin ( λ (cid:12) − Ω ) (44) × c − (cid:18) e s − s (cid:19) sin ε sin λ (cid:12) + c ( cos ε − ) × sin ( λ (cid:12) + Ω ) (cid:105) + η cos ω (cid:104) ( cos ε + ) × cos ( λ (cid:12) − Ω ) + ( − cos ε ) cos ( λ (cid:12) + Ω ) (cid:105)(cid:41) d M d t = n + na e + e F (cid:110) sin ω (cid:104) c ( cos ε + ) (45) × sin ( λ (cid:12) − Ω ) + c ( cos ε − ) sin ( λ (cid:12) + Ω )+ s sin ε sin λ (cid:12) (cid:105) + cos ω (cid:104) ( cos ε + ) × cos ( λ (cid:12) − Ω ) + ( − cos ε ) cos ( λ (cid:12) + Ω ) (cid:105)(cid:111) ρ and the cross-sectional area A of the spacecraft in the direction of motion. The drag forceper unit of mass m is α drag = − ( / ) n d V , where V is velocity of the spacecraft relative to the atmo-sphere, of modulus V , we abbreviated n d = ρ BV > , (46)and B = ( A / m ) C drag , is the so-called ballistic coefficient, inwhich the dimensionless drag coefficient C drag ranges from1.5–3.0 for a typical satellite. Note that n d ≡ n d ( t ) .A reasonable approximation of the relative velocity isobtained with the assumption that the atmosphere co-rotates with the earth. Then, from the derivative of a vector in arotating frame, V = d r d t − ω ⊕ × r . We further take ω ⊕ = n ⊕ k , in the direction of the earth’s ro-tation axis, and compute its projections in the radial, normal,and bi-normal directions as ω ⊕ = R ( θ ) R ( I ) n ⊕ . Then, the velocity components in the radial, normal, andbi-normal direction relative to a rotating atmosphere are V = R ( Θ / r ) − rn ⊕ cos Irn ⊕ cos θ sin I , where R = d r d t = Θ p e sin f , Θ = r d θ d t = √ µ p . Models giving the atmospheric density are usually com-plex. Furthermore, since the atmospheric density depends onthe solar flux which is not easily predictable, reliable pre-dictions of the disturbing effects caused by the atmosphericdrag are not expected for long-term propagation. Hence, theaim is rather to show the effect that the atmospheric dragmight have in the orbit, as opposite from a drag-free model.Therefore, to speed evaluation of the semi-analytical prop-agator, we take advantage of the simplicity of the Harris-Priester atmospheric density model [13], which is imple-mented with the modifications of [28].After replacing α drag into Gauss planetary equations, thelong-term effects are computed by averaging the equationsover the mean anomaly, viz.d a d t = − a η π (cid:90) π n d (47) × (cid:16) + e cos f + e − n ⊕ n η c (cid:17) d M d e d t = − π (cid:90) π n d (48) × (cid:104) e + cos f − δ c (cid:16) e + cos f − e f (cid:17)(cid:105) d M d I d t = − s π (cid:90) π n d δ cos θ d M (49)d Ω d t = −
12 12 π (cid:90) π n d δ sin θ cos θ d M (50)d ω d t = − c d Ω d t (51) − π (cid:90) π n d e sin f (cid:104) − δ c (cid:16) + e f (cid:17)(cid:105) d M d M d t = n + π (cid:90) π n d e η ra sin f d M (52) + π (cid:90) π n d η e sin f (cid:104) − δ c (cid:16) + e f (cid:17)(cid:105) d M in which δ = ( n ⊕ / n )( r / p ) η .Both the relative velocity with respect to the rotating at-mosphere V , and the atmospheric density ρ are naturally ex-pressed as a function of the true anomaly [28], then it hap-pens that n d ≡ n d ( f ) from the definition of n d in Eq. (46).Hence, the quadratures above are conveniently integrated in f rather than in M using the differential relation in Eq. (24).Besides, due to the complex representation of the atmos-pheric density, these quadratures are evaluated numerically. To illustrate the performance of the mean elements theorywe describe two test cases: one for a Molniya-type orbit,and the other for a SymbolX-type orbit, in which the propa-gations are extended to 100 years.A full account of the different tests that have been car-ried out in the development of the HEOSAT software can beconsulted in [22]. The test cases include GTO, super GTO,and SSTO orbits, as well as Tundra orbits, and orbits of thetelescope satellites’ missions Integral and XMM-Newton.The semi-analytical theory generally runs one or two ordersof magnitude faster than the Cowell integration, althoughthese ratios notably reduce when the atmospheric drag has anon-negligible effect. In spite of that, in all the tested casesHEOSAT runs more than 5 times faster than the numericalintegration.9.1 Molniya orbitThe first test presented is for a Molniya type orbit. The ini-tial conditions used in the test correspond to the osculatingelements a = . e = . I = . Ω = . ω =
280 deg M = EOSAT: A mean elements orbit propagator program for Highly Elliptical Orbits 11
Fig. 2, the numerical reference and the mean elements prop-agation fit quite well, with a slight shift to the right of themean elements orbit. This shift is due to the difference be-tween the osculating and mean elements used in both kind ofpropagations. Indeed, as far as the conversion from osculat-ing to mean elements is not implemented in the current ver-sion HEOSAT, the initial mean elements used in launchingthe semi-analytical propagation do not correspond exactlyto the initial conditions used in the propagation of the os-culating reference orbit, even though the semi-major axis ofthe mean elements propagation has been manually adjustedto the mean value of 26653 . ∼ . a = . e = . I = . Ω = .
351 deg ω = − .
992 deg M = (cid:72) years (cid:76) S e m i (cid:45) m a j o r a x i s (cid:72) k m (cid:76) (cid:72) years (cid:76) E cce n t r i c it y (cid:72) years (cid:76) I n c li n a ti on (cid:72) d e g (cid:76) (cid:72) years (cid:76) R AAN (cid:72) d e g (cid:76) (cid:72) years (cid:76) A r gu m e n t o f t h e p e r i g ee (cid:72) d e g (cid:76) Fig. 2
Time history of the orbit elements of the Molniya orbit. Dots:mean elements propagation; gray line: numerical reference about 1000 km. These irregular variations are caused by themoon third-body perturbation, that make the SymbolX orbitto experience different resonances, which include a 1:7 reso-nance of the Laplace type, due to the orbital period of 4 days,as well as secular resonances of the Kozai type, cf. [14].Because of the irregularities in the time history of theosculating semi-major axis, we did not perform any adjust-ment in the computation of the HEOSAT mean semi-majoraxis, and the osculating initial elements are directly usedas mean initial elements for launching the semi-analyticalpropagation. In spite of that, the time history of the mean or- (cid:45) (cid:45)
500 Time (cid:72) years (cid:76) (cid:68) a (cid:72) k m (cid:76) (cid:45) (cid:45) (cid:45) (cid:72) years (cid:76) (cid:68) e (cid:45) (cid:45) (cid:72) years (cid:76) (cid:68) I (cid:72) d e g (cid:76) (cid:72) years (cid:76) (cid:68) (cid:87) (cid:72) d e g (cid:76) (cid:45) (cid:45) (cid:72) years (cid:76) (cid:68) Ω (cid:72) d e g (cid:76) Fig. 3
Errors between the numerical reference and the mean elementspropagation: Molniya case. bital elements shows very good agreement with the osculat-ing elements provided by the numerical reference. The detailon Fig. 5 shows that this agreement extends for more than70 years, and the discrepancies become important passed 95years. Again, notable improvements are expected when themean elements theory is completed with the analytic trans-formation from osculating to mean elements.
10 Conclusions
HEO propagation is a challenging problem due to the differ-ent perturbations that have an effect in highly elliptical or-bits, the relative influence of which may notably vary alongthe orbit. However, modern tools and methods allow to ap-proach the problem by means of analytical methods. Indeed,using perturbation theory we succeeded in the implementa-tion of a fast and efficient semi-analytical propagator whichis able to capture the main frequencies of the HEO motionover long time spans, even in extreme cases, as corroboratedwith the tests performed on the SymbolX orbit. In particu-lar, we used the Lie transforms method, which is standard (cid:72) years (cid:76) S e m i (cid:45) m a j o r a x i s (cid:72) k m (cid:76) (cid:72) years (cid:76) E cce n t r i c it y (cid:72) years (cid:76) I n c li n a ti on (cid:72) d e g (cid:76) (cid:72) years (cid:76) R AAN (cid:72) d e g (cid:76) (cid:72) years (cid:76) A r gu m e n t o f t h e p e r i g ee (cid:72) d e g (cid:76) Fig. 4
Time history of the orbit elements of the SymbolX. Dots: ana-lytical propagation; gray line: numerical reference these days in the construction of perturbation theories. Thismethod is specifically designed for automatic computationby machine, and is easily implemented with modern, com-mercial, general purpose software.Future evolutions of the semi-analytical theory shouldincorporate the transformation from osculating to mean ele-ments, in this way enhancing the precision of the mean el-ements predictions based on it. Also, in spite of commonHEO orbits are not affected by singularities, a reformula-tion in non-singular variables will make the orbit propagator
EOSAT: A mean elements orbit propagator program for Highly Elliptical Orbits 13 (cid:45) (cid:45) (cid:72) years (cid:76) (cid:68) a (cid:72) k m (cid:76) (cid:45) (cid:45) (cid:45) (cid:72) years (cid:76) (cid:68) e (cid:45) (cid:45) (cid:45) (cid:72) years (cid:76) (cid:68) I (cid:72) d e g (cid:76) (cid:45) (cid:45) (cid:45) (cid:72) years (cid:76) (cid:68) (cid:87) (cid:72) d e g (cid:76) (cid:45) (cid:45) (cid:45) (cid:72) years (cid:76) (cid:68) Ω (cid:72) d e g (cid:76) Fig. 5
Errors between the numerical reference and the mean elementspropagation: SymbolX case. software more versatile, widening its scope to the propaga-tion of the majority of objects in a catalogue of earth satelliteand debris orbits.
Acknowledgements
Partial support is acknowledged from the Min-istry of Economic Affairs and Competitiveness of Spain, via ProjectsESP2013-41634-P (M.L.) and ESP2014-57071-R (M.L. and J.F.S.).
A Some useful partial derivatives in elliptic motion
When dealing with automatic manipulation of literal expressions, itresults practical to limit the symbolic algebra to the basic arithmeticoperations, to wit, addition, subtraction, multiplication and division—integer powers being a particular case of multiplication. However,square roots and trigonometric functions appear naturally in the for-mulation of the perturbation in canonical variables. To avoid dealingexplicitly with square roots, it is wise to be equipped with a battery ofpartial derivatives that ease handling and simplifying symbolic expres-sions.Since our perturbation approach relies on the use of Delaunay andpolar-nodal canonical variables, the partial derivatives of the classicalKeplerian variables ( a , e , I , Ω , ω , M ) , standing for semi-major axis, ec-centricity, inclination, right ascension of the ascending node, argumentof the periapsis and mean anomaly, respectively, as well as the usualfunctions of the Keplerian variables – the mean motion n = (cid:112) µ / a – the parameter (or semilatus rectum ) p = a ( − e ) – the eccentricity function η = √ − e – the cosine of the inclination c = cos I – the sine of the inclination s = √ − c – the eccentric anomaly u , given by M = u − e sin u – the true anomaly f , given by ( − e ) tan f = η tan u – the radial distance r = a ( − e cos u ) = p / ( + e cos f ) – the radial velocity R = an ( a / r ) e sin u = an ( e / η ) sin f – the argument of the latitude θ = f + ω – the projections of the eccentricity vector in the orbital frame: k = e cos f = − + p / r , and q = e sin f = R η / ( na ) are provided both in the Delaunay chart ( (cid:96), g , h , L , G , H ) and in thepolar-nodal chart ( r , θ , ν , R , Θ , N ) .We found convenient to express all the partial derivatives by meansof the Keplerian functions: p , n , e , η , k , q , s , c from whose definition it is obtained L = n p / η G = L η = n p / η = Θ H = Gc = n p c / η = NR = n pq / η r = p / ( + k ) Recall also that, from the ellipse geometry, the following relations ap-ply, cf. Eq. (34),sin f = ( a / r ) η sin u , cos f = ( a / r )( cos u − e ) . Finally, it worths to mention that the use of logarithmic deriva-tives is helpful in finding the differentials that eased the computationof the partial derivatives. Note that only non-vanishing derivatives arepresented.
A.1 With respect to Delaunay variables
A.1.1 Orbital elements and related functions – Semi-major axis a :d a d L = L µ = η n p (55) – Mean motion n :d n d L = − η p , (56) – Parameter p :d p d G = G µ = η n p . (57) – Eccentricity function η :d η d L = η (cid:18) − L (cid:19) = − η n p (58)d η d G = η (cid:18) G (cid:19) = η n p (59) – Eccentricity e :d e d L = η e L = η en p (60)d e d G = − η e G = − η en p (61)4 Martin Lara et al. – Cosine of inclination c :d c d G = − cG = − c η n p (62)d c d H = cH = G = η n p (63) – Sine of inclination s :d s d G = c s η n p (64)d s d H = − cs η n p (65) – Eccentric anomaly u :d u d (cid:96) = pr η = + k η (66)d u d L = R η e n p = q η e n p (67)d u d G = − R η e n p = − q η e n p (68) A.1.2 Polar-nodal variables and related functions – Radial distance r :d r d (cid:96) = Rn = pq η (69)d r d L = η e n (cid:18) e rp + p − r (cid:19) = η n p (cid:18) + k − ke (cid:19) (70)d r d G = η e n (cid:18) r − p (cid:19) = η ke n p (71) – Radial velocity R :d R d (cid:96) = pn η p r (cid:16) pr − (cid:17) = pn η k ( + k ) (72)d R d L = η Rnr (cid:18) e − r p (cid:19) = η qp (cid:20) ( + k ) e − (cid:21) (73)d R d G = − η e Rnr = − ( + k ) e qp (74) – True anomaly f :d f d (cid:96) = η ran ( p − r ) d R d (cid:96) = p r η = ( + k ) η (75)d f d L = R η e n p (cid:18) p + r (cid:19) = q η e n p ( + k ) (76)d f d G = − R η e n p (cid:18) p + r (cid:19) = − q η e n p ( + k ) (77) – Argument of the latitude θ :d θ d (cid:96) = d f d (cid:96) (78)d θ d g = θ d L = d f d L (80)d θ d G = d f d G (81) – Modulus of the angular momentum Θ :d Θ d G = – Argument of the node ν :d ν d h = – Polar component of the angular momentum N :d N d H = – Eccentricity vector k :d k d (cid:96) = − pRr n = − q η ( + k ) (85)d k d L = η nr (cid:18) ke − k + (cid:19) (86)d k d G = − η nr (cid:18) ke − k + (cid:19) (87) – Eccentricity vector q :d q d (cid:96) = η p r (cid:16) pr − (cid:17) = k η p r = k η ( + k ) (88)d q d L = q η nr (cid:18) e − r p (cid:19) = q η n p (cid:20) ( + k ) e − (cid:21) (89)d q d G = − q η nr (cid:18) e − r p (cid:19) = − q η n p (cid:20) ( + k ) e − (cid:21) (90) A.2 With respect to polar-nodal variables
A.2.1 Delaunay variables – Delaunay action L :d L d r = − Θ p k η ( + k ) = − np k η ( + k ) (91)d L d R = Θ n qp = p q η (92)d L d Θ = p η r = ( + k ) η (93) – Modulus of the angular momentum Θ :d G d Θ = – Polar component of the angular momentum N :d H d N = – Mean anomaly (cid:96) :d (cid:96) d r = η qr (cid:18) + ke − + k (cid:19) (96)d (cid:96) d R = η qR (cid:18) ke − + k (cid:19) (97)d (cid:96) d Θ = − η q Θ + ke (98) – Argument of the perigee g :d g d r = − qr + ke (99)d g d θ = g d R = − qR ke (101)d g d Θ = q Θ + ke (102) – Right ascension of the ascending node h :d h d ν = A.2.2 Orbital elements and related functions – Parameter p :d p d Θ = p Θ = η n p (104) – Eccentricity vector k :d k d r = − pr = − ( + k ) p (105)d k d Θ = pr Θ = qr R = ( + k ) η n p (106) – Eccentricity vector q :d q d R = qR = p Θ = η n p (107)d q d Θ = q Θ = q R p = q η n p (108) – Eccentricity e :d e d r = − ke pr = − ke ( + k ) p (109)d e d R = q eR = q η en p (110)d e d Θ = η k ( + k ) + q en p (111) – Eccentricity function η :d η d r = k η pr = k η ( + k ) p (112)d η d R = − q η R = − q η n p (113)d η d Θ = = − η k ( + k ) + q n p (114) – Semi-major axis a :d a d r = − k η ( + k ) (115)d a d R = pq R η = qn η (116)d a d Θ = p η Θ r = ( + k ) n p η (117) – Mean motion n :d n d r = nk η ( + k ) p (118)d n d R = − nq R η = − q η p (119)d n d Θ = − n p η Θ r = − η r (120) – True anomaly f :d f d r = qr + ke (121)d f d R = qR ke (122)d f d Θ = − q Θ + ke (123) – Eccentric anomaly u :d u d r = η pq (cid:20) + e − ke k η ( + k ) (cid:21) ( + k ) (124)d u d R = anR (cid:18) − r − ar d ee d R − d aa d R (cid:19) (125)d u d Θ = anR (cid:18) − r − ar d ee d Θ − d aa d Θ (cid:19) (126) – Equation of the center φ :d φ d r = qr (cid:18) + k + η + η + k (cid:19) (127)d φ d R = qR (cid:18) k + η + η + k (cid:19) (128)d φ d Θ = − q Θ + k + η (129) – Cosine of inclination c :d c d Θ = − c Θ = − cqR p = − c η n p (130)d c d N = Θ = qR p = η n p (131) – Sine of inclination s :d s d Θ = c η sn p (132)d s d N = − c η sn p (133) B Tables of coefficients
The coefficients of the trigonometric series used by HEOSAT are pro-vided in following tables
References
1. Armellin, R., San-Juan, J.F., Lara, M.: End-of-life disposal of highelliptical orbit missions: The case of INTEGRAL. Advances inSpace Research (3), 479–493 (2015). DOI 10.1016/j.asr.2015.03.020. Advances in Asteroid and Space Debris Science and Tech-nology - Part 12. Breakwell, J.V., Vagners, J.: On Error Bounds and Initialization inSatellite Orbit Theories. Celestial Mechanics , 253–264 (1970).DOI 10.1007/BF012294993. Chapront, J., Francou, G.: The lunar theory ELP revisited. Intro-duction of new planetary perturbations. Astronomy and Astro-physics , 735–742 (2003). DOI 10.1051/0004-6361:200305294. Chapront-Touze, M., Chapront, J.: ELP 2000-85 - A semi-analytical lunar ephemeris adequate for historical times. Astron-omy and Astrophysics , 342–352 (1988)5. Coffey, S.L., Neal, H.L., Segerman, A.M., Travisano, J.J.: Ananalytic orbit propagation program for satellite catalog mainte-nance. In: K.T. Alfriend, I.M. Ross, A.K. Misra, C.F. Peters (eds.)AAS/AIAA Astrodynamics Conference 1995, Advances in the As-tronautical Sciences , vol. 90, pp. 1869–1892. American Astronau-tical Society, Univelt, Inc., USA (1996)6. Danielson, D.A., Neta, B., Early, L.W.: Semianalytic Satel-lite Theory (SST): Mathematical algorithms. Technical ReportNPS-MA-94-001, Naval Postgraduate School, Naval Postgradu-ate School, Monterey, CA. Dept. of Mathematics. (1994)7. Deprit, A.: Canonical transformations depending on a small pa-rameter. Celestial Mechanics (1), 12–30 (1969). DOI 10.1007/BF012306298. Deprit, A.: The elimination of the parallax in satellite theory.Celestial Mechanics (2), 111–153 (1981). DOI 10.1007/BF012291929. Efroimsky, M.: Gauge Freedom in Orbital Mechanics. Annals ofthe New York Academy of Sciences , 346–374 (2005). DOI10.1196/annals.1370.01610. Exertier, A.: Orbitographie des satellites artificiels sur de grandesperiodes de temps. Possibilites d’applications. PhD. Thesis. Ob-servatoire de Paris, Paris (1988)6 Martin Lara et al. Table 1
Eccentricity polynomials Q m , k + m mod 2 in Eq. (27); Q m , m − = Q m , m − = m − + e and Q , = + e + e . k m = m = m = m =
100 3 ( + e + e ) ( + e + e + e ) ( + e + e + e ) ( + e + e + e + e ) ( + e + e ) + e + e ( + e + e + e ) ( + e + e ) Table 2
Even inclination polynomials B m , k in Eq. (27). k m = m = m = m = (cid:0) c − (cid:1) − (cid:0) c − c + (cid:1) (cid:0) c − c + c − (cid:1) − (cid:0) c − c + c − c + (cid:1) − (cid:0) c − (cid:1) (cid:0) c − c + (cid:1) − (cid:0) c − c + c − (cid:1) (cid:0) c − (cid:1) − (cid:0) c − c + (cid:1) − (cid:0) c − (cid:1) m = (cid:0) c − c + c − c + c − (cid:1) (cid:0) c − c + c − c + (cid:1) (cid:0) c − c + c − (cid:1) (cid:0) c − c + (cid:1) (cid:0) c − (cid:1) Table 3
Odd inclination polynomials B m + , k + in Eq. (27). k m = m = m = m = − (cid:0) c − (cid:1) (cid:0) c − c + (cid:1) − (cid:0) c − c + c − (cid:1) (cid:0) c − c + c − c + (cid:1) (cid:0) c − (cid:1) − (cid:0) c − c + (cid:1) (cid:0) c − c + c − (cid:1) − (cid:0) c − (cid:1) (cid:0) c − c + (cid:1) (cid:0) c − (cid:1)
11. Garfinkel, B.: On the motion of a satellite of an oblate planet.The Astronomical Journal (1257), 88–96 (1958). DOI 10.1086/10769712. Golikov, A.R.: THEONA—a numerical-analytical theory of mo-tion of artificial satellites of celestial bodies. Cosmic Research (6), 449–458 (2012). DOI 10.1134/S001095251206002013. Harris, I., Priester, W.: Time-Dependent Structure of the Upper At-mosphere. Journal of Atmospheric Sciences , 286–301 (1962).DOI 10.1175/1520-0469(1962)019 (cid:104) (cid:105) , in press (2016). DOI 10.1007/s10569-016-9736-6. URL http://arxiv.org/pdf/1605.00525.pdf
16. Hoots, F.R., Roehrich, R.L.: Models for Propagation of the NO-RAD Element Sets. Project SPACETRACK, Rept. 3, U.S. AirForce Aerospace Defense Command, Colorado Springs, CO(1980)17. Kaufman, B.: First order semianalytic satellite theory with recov-ery of the short period terms due to third body and zonal per-turbations. Acta Astronautica (5–6), 611 – 623 (1981). DOI10.1016/0094-5765(81)90108-918. Kaula, W.M.: Theory of satellite geodesy. Applications of satel-lites to geodesy. Blaisdell, Waltham, Massachusetts (1966)19. Kozai, Y.: Effects of Solar Radiation Pressure on the Motion of anArtificial Satellite. SAO Special Report , 25–34 (1961) 20. Kozai, Y.: Second-Order Solution of Artificial Satellite Theorywithout Air Drag. The Astronomical Journal (7), 446–461(1962)21. Lara, M.: Simplified Equations for Computing Science OrbitsAround Planetary Satellites. Journal of Guidance Control Dynam-ics (1), 172–181 (2008). DOI 10.2514/1.3110722. Lara, M., San-Juan, J., Hautesserres, D.: Semi-analytical propaga-tor of high eccentricity orbits. Technical Report R-S15/BS-0005-024, Centre National d’ ´Etudes Spatiales, 18, avenue Edouard Be-lin - 31401 Toulouse Cedex 9, France (2016)23. Lara, M., San-Juan, J.F., Folcik, Z.J., Cefola, P.: Deep ResonantGPS-Dynamics Due to the Geopotential. The Journal of theAstronautical Sciences (4), 661–676 (2011). DOI 10.1007/BF0332153624. Lara, M., San-Juan, J.F., L´opez, L.M., Cefola, P.J.: On the third-body perturbations of high-altitude orbits. Celestial Mechanicsand Dynamical Astronomy , 435–452 (2012). DOI 10.1007/s10569-012-9433-z25. Lara, M., San-Juan, J.F., L´opez-Ochoa, L.M.: Delaunay variablesapproach to the elimination of the perigee in Artificial SatelliteTheory. Celestial Mechanics and Dynamical Astronomy (1),39–56 (2014). DOI 10.1007/s10569-014-9559-226. Lara, M., San-Juan, J.F., L´opez-Ochoa, L.M.: Proper AveragingVia Parallax Elimination (AAS 13-722). In: Astrodynamics 2013, Advances in the Astronautical Sciences , vol. 150, pp. 315–331.American Astronautical Society, Univelt, Inc., USA (2014)27. Lara, M., Vilhena de Moraes, R., Sanchez, D.M., Prado, A.F.B.A.:Efficient computation of short-period analytical corrections dueto third-body effects (AAS 15-295). In: Proceedings of the 25thAAS/AIAA Space Flight Mechanics Meeting, Williamsburg, VA,January 11 – 15, 2015,
Advances in the Astronautical Sciences ,EOSAT: A mean elements orbit propagator program for Highly Elliptical Orbits 17
Table 4
Eccentricity polynomials A m , j in Eq. (38). m j = ( + e ) − e e ( + e ) e ( + e + e ) e ( + e ) e e ( + e + e ) e ( + e ) e + e + e + e e + e + e e + e e Table 5
Inclination polynomials P m , j , l in Eq. (38) ( χ = c ± P , j , l P , j , l P , j , l l j = j = j = j = j = j = j = ( c − ) − s − ( c − ) s − s − ( c − c + ) − ( c − ) s − s ± cs χ s − χ ( c ∓ c − ) − χ s c ( − c ) s χ ( c ∓ c − ) s χ s ± − s χ χ ( c ∓ ) s − χ s ( c − ) s χ ( c ∓ c + ) χ s ± χ s χ cs χ ( c ∓ ) s − χ s ± − s − χ s − χ Table 6
Inclination polynomials P , j , l in Eq. (38) ( χ = c ± l j = j = j = (cid:0) c − c + (cid:1) s ( c − ) s s ± χ ( c ∓ c − c ± c + ) χ ( c ∓ c − ) s χ s ± − χ ( c ∓ c − c ± ) s χ ( c ∓ c + ) s χ s ± χ ( c ∓ c − ) s χ ( c ∓ )( c ∓ ) χ s ± − χ ( c ∓ ) s − χ ( c ∓ ) s χ s ± χ s χ s χ vol. 155, pp. 437–455. American Astronautical Society, Univelt,Inc., USA (2015)28. Long, A.C., Cappellari, J.O., Velez, C.E., Fluchs, A.J.: Mathe-matical Theory of the Goddard Trajectory Determination Sys-tem. Technical Report FDD/552-89/001, National Aeronauticsand Space Administration, Goddard Space Flight Center, Green-belt, MD. (1989)29. McClain, W.D.: A Recursively Formulated First-Order Semiana-lytic Artificial Satellite Theory Based on the Generalized Methodof Averaging, Volume 1: The Generalized Method of AveragingApplied to the Artificial Satellite Problem, 2nd edn. NASA CR-156782. NASA, Greenbelt, Maryland (1977)30. Meeus, J.: Mathematical astronomy morsels. Willmann-Bell,Richmond, VA (1997)31. Meeus, J.: Astronomical algorithms, 2nd edn. Willmann-Bell,Richmond, VA (1998)32. Montenbruck, O., Gill, E.: Satellite Orbits. Models, Methods andApplications. Physics and Astronomy. Springer-Verlag, Berlin,Heidelberg, New York (2001)33. San-Juan, J.F., Lara, M., Ferrer, S.: Phase Space Structure AroundOblate Planetary Satellites. Journal of Guidance Control Dynam-ics , 113–120 (2006). DOI 10.2514/1.1338534. Sterne, T.E.: The gravitational orbit of a satellite of an oblateplanet. The Astronomical Journal , 28–40 (1958). DOI10.1086/10767335. Vallado, D.A.: Fundamentals of Astrodynamics and Applications,2nd edn. Microcosm (2001)8 Martin Lara et al. Table 7
Inclination polynomials P , j , l in Eq. (38) ( χ = c ± l j = j = j = j = − ( c − c + c − ) − ( c − c + ) s − ( c − ) s − s ± − c ( c − c + ) s χ ( c ∓ c − c ± c + ) s χ ( c ∓ c − ) s χ s ± − ( c − c + ) s − χ ( c ∓ c + c ± c − ) − χ ( c ∓ c + ) s − χ s ± − c ( c − ) s − χ ( c ∓ c + c ± ) s χ ( c ∓ c + ) s χ s ± ( − c ) s − χ ( c ∓ c + ) s − χ ( c ∓ c + ) − χ s ± − cs ( ∓ c ) χ s − χ ( c ∓ ) s χ s ± s χ s χ s χ Table 8
Third-body direction polynomials in Eq. (38); u ≡ u (cid:63) , v ≡ v (cid:63) , w ≡ w (cid:63) . m l S m , l T m , l − + w ± − vw ± uw ± u − v ± uv w ( w − ) ± ± u ( w − ) v ( w − ) ± ± uvw w ( v − u ) ± ± u ( u − v ) − v ( v − u ) − w + w ± vw ( − w ) ± uw ( − + w ) ± ( u − v )( − + w ) ± uv ( − + w ) ± v ( − u + v ) w ± u ( u − v ) w ± ( u − u v + v ) ± uv ( u − v ) w ( − w + w ) ± ± u ( − w + w ) v ( − w + w ) ± ± uvw ( − + w ) ( u − v ) w ( − w ) ± ± u ( u − v )( − w ) v ( − u + v )( − + w ) ± ± uv ( − u + v ) w ( u − u v + v ) w ± ± u ( u − u v + v ) v ( u − u v + v ) − w + w − w ± v ( − w + w ) w ∓ u ( − w + w ) w ± ( u − v )( − w + w ) ± uv ( − w + w ) ± v ( u − v )( − w ) w ± u ( u − v )( − + w ) w ± ( − w )( u − u v + v ) ± ( − w ) uv ( u − v ) ± v ( u − u v + v ) w ∓ u ( u − u v + v ) w ± ( v − u )( u − u v + v ) ∓ ( u − u v + uv ))