- Email: [email protected]
Computer algebra methods for developing an artificial satellite motion theory within elliptic functions
Computer Physics Communications ELSEVIER
Computer Physics Communications 118 (1999) 17-20
Computer algebra methods for developing an artificial satellite motion theory within elliptic functions A k m a l A. V a k h i d o v , N i c k o l a y N. Vasiliev Institute of Theoretical Astronomy of Russian Academy of Sciences, Nab. Kutuzova, 10, 191187 Saint Petersburg, Russian Federation Received 16 May 1998; revised 18 September 1998
Abstract The problem of developing a semianalytical theory of motion for artificial Earth satellites with large eccentricities is considered using elliptic function techniques. Computer-algebraic methods for symbolic computation of the coefficients of expansion for the disturbing function in multiples of the elliptic anomaly are described. Possibilities of computing analytically the Hamiltonian of a satellite perturbed motion by means of computer algebra are shown. The numerical technique to compute elliptic Hansen coefficients with the help of the Chebyshev approximation is presented. (g) 1999 Published by Elsevier Science B.V. PACS: 95.10.Ce; 95.40.+s; 02.30.Mv Keywords: Celestial mechanics; Artificial Earth satellites; Approximations and expansions
1. Introduction One of the important problems of modern celestial mechanics is the development of efficient theories of motion for artificial Earth satellites with large eccentricities. For solving this problem it is very convenient to use the expansions of the Hamiltonian of the satellite disturbed motion into series in multiples of the elliptic anomaly [ 1 ]. The coefficients of these expansions (called elliptic Hansen coefficients) can be determined in a closed form via complete elliptic integrals of the first kind with the eccentricity as a modulus [2]. The convergence of expansions in multiples of the elliptic anomaly is very fast even for large eccentricities, which gives a possibility of using them even for highly eccentric orbits, when the eccentricity is larger than the Laplace limit [3,4]. (It is impossible for the
classical expansions of satellite disturbing function.) In this paper we consider some problems which appear in constructing a satellite motion theory within elliptic functions. Solving these problems we obtain a possibility of developing a highly efficient theory of motion for artificial Earth satellites with large eccentricities.
2. Some problems of computation of elliptic Hansen coefficients In Refs. [ 5,6] the authors have obtained the analytical expressions for the Hamiltonian presenting the motion of an artificial Earth satellite under the influence of the Earth gravitational field and lunar-solar gravitational perturbations in the form of trigonometric series in multiples of the elliptic anomaly. These expansions contain some special functions of celestial mechanics
0010-4655/99/$ - see front matter (g) 1999 Published by Elsevier Science B.V. All rights reserved. PII S 0 0 1 0 - 4 6 5 5 ( 9 8 ) 0 0 1 7 7 - 5
18
A.A. Vakhidov, N.N. Vasiliev /Computer Physics Communications 118 (1999) 17-20
(inclination functions, Hansen coefficients, and others), and there are no difficulties to compute them as many efficient algorithms are elaborated for their computations (see also Refs. [7-9] ). Conversely, the problem of evaluation of elliptic Hansen coefficients (which are included into this Hamiltonian, too) faces a lot of difficulties. Let us consider some of them and show some possible ways for solving the problems presented here. We notice that for the prediction of motion of highly elliptic Earth satellites we need to take into consideration, as a rule, perturbations from many harmonics of the geopotential (near to the perigee of the orbit) and lunar-solar perturbations including terms of high order of disturbing function expansion with respect to the parallax (near to the orbit apogee). As a result, in the process of prediction of motion of a satellite the number of elliptic Hansen coefficients which should be computed is equal to many thousands. Therefore, it is necessary to develop special techniques for saving computer memory and computing time in constructing the motion theory. Another problem in computing the elliptic Hansen coefficients is the usage of exact analytical expressions for their definition. Simple computer-algebraic experiments show that the exact analytical formulae for these coefficients obtained by means of recurrence relations [2] are very cumbersome even for geopotential harmonics of the fourth or fifth order. It leads to some inconveniences to use the analytical expressions for the elliptic coefficients in the satellite theory. Moreover, these experiments show that numerical values of derivatives of elliptic Hansen coefficients with respect to the eccentricity are often larger than the numerical values of the coefficients themselves. Therefore, a small variation of eccentricity causes substantial variations of elliptic Hansen coefficients, and we need to evaluate over and over again all coefficients at each step of integration of the satellite motion equations. We have also experienced that computing the derivatives of elliptic coefficients is more complicated than evaluating the coefficients themselves, because the formulae for these derivatives are substantially more cumbersome [ 10]. For solving the problems described above we have prepared the computing techniques to overcome those difficulties.
3. Chebyshev approximation of elliptic Hansen coefficients In order to solve the problem of efficient computation of elliptic Hansen coefficients let us consider now some of their properties. From computer experiments we can see that all the elliptic coefficients and their first and second derivatives show no change of sign in the domain of large eccentricities. This means that these coefficients are monotonic functions of eccentricity. Therefore, it is very convenient to accomplish the "Chebyshev approximation" of elliptic Hansen coefficients in the interval of possible values of the eccentricity. Then for computation of coefficients one can use the compact approximating polynomials instead of very cumbersome exact analytical formulae and keep in computer memory the numerical values of approximating coefficients of the "Chebyshev expansions" only. For approximation it is very efficient to use the standard MAPLE-package ORTHOPOLY [ 11 ]. The accuracy of approximation is here one of the input parameters of the MAPLE-function CHEBYSHEV. Moreover, the optimal number of approximation points is defined in this case by means of the system MAPLE. We notice that the mean computing time of approximation of a coefficient using MAPLE on an IBM PC AT-486 (66 MHz) is not larger than 1 s. Furthermore, the coefficients of approximating polynomials are saved in order to be used in the process of computation of perturbations in the satellite motion. The recurrence process of computation of elliptic Hansen coefficients can be realized in such a way that each coefficient is computed symbolically using the "Chebyshev approximations" of other coefficients included into the corresponding recurrence formula and determined at the previous steps of the recurrences. After this we approximate the coefficient computed symbolically by means of Chebyshev polynomials and at the next steps of recurrences we compute this coefficient using its "Chebyshev approximation" only. As a result, we obtain the "Chebyshev expansions" for all elliptic Hansen coefficients, which should be evaluated in the satellite motion theory. In such a way we solve two problems simultaneously. First, we obtain the compact analytical expressions for elliptic Hansen coefficients instead of very cumbersome exact analytical formulae for them. Sec-
A.A. Vakhidov, N.N. Vasiliev/ Computer Physics Communications 118 (1999) 17-20
ond, we have now a possibility to do recurrences prior to the integration of differential equations of satellite motion. After this, in the satellite motion theory we evaluate the elliptic Hansen coefficients using their "Chebyshev approximations". It allows us to save computer memory and to decrease the time of computations. The above described scheme of computation of elliptic Hansen coefficients can be used also for computation of derivatives of these coefficients which are included into equations of motion, too.
4. Some important remarks Let us make now some remarks concerning the accuracy of computing the elliptic Hansen coefficients in the satellite motion theory. For constructing analytical and semianalytical theories of motion for artificial Earth satellites we have to evaluate the coefficients with an accuracy higher than the order of terms eliminated from consideration. This means that for theories of the first order we should evaluate the elliptic Hansen coefficients with the relative accuracy 10 -4, because the perturbations of the second order are approximately 103 times smaller (in their absolute values) than the perturbations of the first order [ 12]. For theories of the second order we should compute the first order perturbations with errors smaller than 10 -7 and the second order perturbations with errors smaller than 10 -4 . Then the errors of the theory will be smaller than the perturbations of the third order. Because the analytical and semianalytical theories are most efficient for prediction of motion of satellites with not a very high accuracy but on a large interval of time, then it is usually enough to take into consideration the perturbations of the first order only [ 13] and only in special cases should we consider the perturbations of the second order. In order to reach the accuracy 10 -7 the approximation of elliptic Hansen coefficients should be made with the help of the Chebyshev approximation up to and including the polynomials of the sixth order, if the interval of approximation is [e0 - 0.05; e0 + 0.05] and the initial value of eccentricity e0 is placed in the domain 0.45 < e0 < 0.85. The interval of eccentricity [e0 - 0.05; e0 + 0.05] is usually enough, if we predict the motion on such intervals of time as one or two
19
months or less, because the eccentricity has no secular perturbations. But if we would like to study the evolution of the orbit within, for example, several years, then we can divide the interval of all possible values of eccentricity by separate fragments and realize the Chebyshev approximation scheme on each interval independently. The summary of our experiments is as follows: ( 1) for computing elliptic Hansen coefficients it is enough to use the Chebyshev approximation of the sixth (maximal, for the emergency, seventh) order; (2) the interval of eccentricity [ e0 - 0.05; e0 + 0.05 ] is enough for doing approximations from the point of view of practical purposes; (3) under conditions (1), (2) we evaluate the elliptic Hansen coefficients with relative errors of 10 -7 which is quite acceptable for all analytical and semianalytical theories of the first and the second order; (4) the above described approach is valid up to the values of eccentricity e < 0.9; there are no limitations concerning critical inclinations, resonances, etc.; (5) this approach gives a possibility to improve the existing analytical theories substantially from the point of view of their compactness (not accuracy). The computer experiments show that the computation of elliptic Hansen coefficients by means of their "Chebyshev approximations" instead of the exact analytical formulae allows one to speed up the process of computing the Hamiltonian (and respectively the right parts of differential equations of motion) by approximately ten times for low harmonics of the geopotential. For high harmonics of the Earth potential the effect of decrease of computing time is more substantial because the exact formulae for elliptic Hansen coefficients are more complicated and more cumbersome for high harmonics but "Chebyshev expansions" do not suffer any complications with the increase of the order of harmonics. Saving computer memory in evaluating elliptic Hansen coefficients and their derivatives gives a possibility of extending the number of geopotential harmonics and the number of terms of the lunar-solar disturbing function which can be taken into consideration.
20
A.A. Vakhidov, N.N. Vasiliev/Computer Physics Communications 118 (1999) 17-20
All our computer investigations have been made on an IBM PC AT computer with an Intel 80486 processor in MS DOS 6.0 operating system. The computer experiments have shown a high efficiency of the described approach and possibility to use it for computation of perturbations in the satellite motion.
Acknowledgements The authors are very grateful to Prof. Andre Deprit for his kind attention to our work and for useful suggestions for improving the text of the manuscript. The authors thank also Prof. V.A. Brumberg for useful discussions.
5. Conclusion References Special computing techniques presented in this paper give the possibility of solving some problems in evaluating elliptic Hansen coefficients and their derivatives from the point of view of saving computer memory and increasing the rapidity of satellite motion computation. On the basis of the described approach one can develop a very efficient satellite motion theory within elliptic functions. At present such a theory is constructed by the authors for an artificial Earth satellite with the eccentricity e = 0.74 [ 14]. The accuracy of this theory is quite acceptable for practical purposes. The results of comparison of this theory with the direct numerical integration of equations of motion in Cartesian coordinates are also presented in Ref. [ 14]. In the future we also plan to apply this theory to other artificial Earth satellites with large eccentricities.
[ 1] E.V. Brumberg, in: Proc. 25th Symp. on Celest. Mech. (Tokyo, 1992) p. 139. [2] E.V. Brumberg, T. Fukushima, Celest. Mech. Dynam. Astron. 60 (1994) 69. 13] N.N. Vasiliev, A.A. Vakhidov, S.A. Klioner, Preprint Institute of Theoretical Astronomy of the Russian Academy of Sciences No. 60 (1996) [in Russian]. [4] S.A. Klioner, A.A. Vakhidov, N.N. Vasiliev, Celest. Mech. Dynam. Astron. 68 (1998) 257. [5] N.N. Vasiliev, A.A. Vakhidov, A.G. Sokolsky, Preprint Institute of Theoretical Astronomy of the Russian Academy of Sciences No. 38 (1994) [in Russian]. [6] A.A. Vakhidov, N.N. Vasiliev, J. Symbolic Comput. 24 (1997) 705. [7] G.E.O. Giacaglia, Celest. Mech. 14 (1976) 515. [8] R.H. Gooding, Celest. Mech. 4 (1971) 91. [9] D.R. Cok, Celest. Mech. 16 (1977) 459. [10] A.G. Sokolsky, A.A. Vakhidov, N.N. Vasiliev, Celest. Mech. Dynam. Astron. 63 (1996) 357. I l l ] B.W. Char, K.O. Geddes, G.H. Gonnet, B.L. Leong, M.B. Monagan, S.M. Watt, MAPLE V Library Reference Manual (Springer, New York, 1993). [12] A.M. Fominov, Bull. ITA 14 (1980) 621 [in Russianl. [13] E.P. Aksenov, Theory of Motion of Artificial Earth Satellites (Nauka, Moscow, 1977) lin Russian]. [ 14] A.A. Vakhidov, N.N. Vasiliev, Astron. J. 112 (1996) 2330.
Computer Physics Communications 118 (1999) 17-20
Computer algebra methods for developing an artificial satellite motion theory within elliptic functions A k m a l A. V a k h i d o v , N i c k o l a y N. Vasiliev Institute of Theoretical Astronomy of Russian Academy of Sciences, Nab. Kutuzova, 10, 191187 Saint Petersburg, Russian Federation Received 16 May 1998; revised 18 September 1998
Abstract The problem of developing a semianalytical theory of motion for artificial Earth satellites with large eccentricities is considered using elliptic function techniques. Computer-algebraic methods for symbolic computation of the coefficients of expansion for the disturbing function in multiples of the elliptic anomaly are described. Possibilities of computing analytically the Hamiltonian of a satellite perturbed motion by means of computer algebra are shown. The numerical technique to compute elliptic Hansen coefficients with the help of the Chebyshev approximation is presented. (g) 1999 Published by Elsevier Science B.V. PACS: 95.10.Ce; 95.40.+s; 02.30.Mv Keywords: Celestial mechanics; Artificial Earth satellites; Approximations and expansions
1. Introduction One of the important problems of modern celestial mechanics is the development of efficient theories of motion for artificial Earth satellites with large eccentricities. For solving this problem it is very convenient to use the expansions of the Hamiltonian of the satellite disturbed motion into series in multiples of the elliptic anomaly [ 1 ]. The coefficients of these expansions (called elliptic Hansen coefficients) can be determined in a closed form via complete elliptic integrals of the first kind with the eccentricity as a modulus [2]. The convergence of expansions in multiples of the elliptic anomaly is very fast even for large eccentricities, which gives a possibility of using them even for highly eccentric orbits, when the eccentricity is larger than the Laplace limit [3,4]. (It is impossible for the
classical expansions of satellite disturbing function.) In this paper we consider some problems which appear in constructing a satellite motion theory within elliptic functions. Solving these problems we obtain a possibility of developing a highly efficient theory of motion for artificial Earth satellites with large eccentricities.
2. Some problems of computation of elliptic Hansen coefficients In Refs. [ 5,6] the authors have obtained the analytical expressions for the Hamiltonian presenting the motion of an artificial Earth satellite under the influence of the Earth gravitational field and lunar-solar gravitational perturbations in the form of trigonometric series in multiples of the elliptic anomaly. These expansions contain some special functions of celestial mechanics
0010-4655/99/$ - see front matter (g) 1999 Published by Elsevier Science B.V. All rights reserved. PII S 0 0 1 0 - 4 6 5 5 ( 9 8 ) 0 0 1 7 7 - 5
18
A.A. Vakhidov, N.N. Vasiliev /Computer Physics Communications 118 (1999) 17-20
(inclination functions, Hansen coefficients, and others), and there are no difficulties to compute them as many efficient algorithms are elaborated for their computations (see also Refs. [7-9] ). Conversely, the problem of evaluation of elliptic Hansen coefficients (which are included into this Hamiltonian, too) faces a lot of difficulties. Let us consider some of them and show some possible ways for solving the problems presented here. We notice that for the prediction of motion of highly elliptic Earth satellites we need to take into consideration, as a rule, perturbations from many harmonics of the geopotential (near to the perigee of the orbit) and lunar-solar perturbations including terms of high order of disturbing function expansion with respect to the parallax (near to the orbit apogee). As a result, in the process of prediction of motion of a satellite the number of elliptic Hansen coefficients which should be computed is equal to many thousands. Therefore, it is necessary to develop special techniques for saving computer memory and computing time in constructing the motion theory. Another problem in computing the elliptic Hansen coefficients is the usage of exact analytical expressions for their definition. Simple computer-algebraic experiments show that the exact analytical formulae for these coefficients obtained by means of recurrence relations [2] are very cumbersome even for geopotential harmonics of the fourth or fifth order. It leads to some inconveniences to use the analytical expressions for the elliptic coefficients in the satellite theory. Moreover, these experiments show that numerical values of derivatives of elliptic Hansen coefficients with respect to the eccentricity are often larger than the numerical values of the coefficients themselves. Therefore, a small variation of eccentricity causes substantial variations of elliptic Hansen coefficients, and we need to evaluate over and over again all coefficients at each step of integration of the satellite motion equations. We have also experienced that computing the derivatives of elliptic coefficients is more complicated than evaluating the coefficients themselves, because the formulae for these derivatives are substantially more cumbersome [ 10]. For solving the problems described above we have prepared the computing techniques to overcome those difficulties.
3. Chebyshev approximation of elliptic Hansen coefficients In order to solve the problem of efficient computation of elliptic Hansen coefficients let us consider now some of their properties. From computer experiments we can see that all the elliptic coefficients and their first and second derivatives show no change of sign in the domain of large eccentricities. This means that these coefficients are monotonic functions of eccentricity. Therefore, it is very convenient to accomplish the "Chebyshev approximation" of elliptic Hansen coefficients in the interval of possible values of the eccentricity. Then for computation of coefficients one can use the compact approximating polynomials instead of very cumbersome exact analytical formulae and keep in computer memory the numerical values of approximating coefficients of the "Chebyshev expansions" only. For approximation it is very efficient to use the standard MAPLE-package ORTHOPOLY [ 11 ]. The accuracy of approximation is here one of the input parameters of the MAPLE-function CHEBYSHEV. Moreover, the optimal number of approximation points is defined in this case by means of the system MAPLE. We notice that the mean computing time of approximation of a coefficient using MAPLE on an IBM PC AT-486 (66 MHz) is not larger than 1 s. Furthermore, the coefficients of approximating polynomials are saved in order to be used in the process of computation of perturbations in the satellite motion. The recurrence process of computation of elliptic Hansen coefficients can be realized in such a way that each coefficient is computed symbolically using the "Chebyshev approximations" of other coefficients included into the corresponding recurrence formula and determined at the previous steps of the recurrences. After this we approximate the coefficient computed symbolically by means of Chebyshev polynomials and at the next steps of recurrences we compute this coefficient using its "Chebyshev approximation" only. As a result, we obtain the "Chebyshev expansions" for all elliptic Hansen coefficients, which should be evaluated in the satellite motion theory. In such a way we solve two problems simultaneously. First, we obtain the compact analytical expressions for elliptic Hansen coefficients instead of very cumbersome exact analytical formulae for them. Sec-
A.A. Vakhidov, N.N. Vasiliev/ Computer Physics Communications 118 (1999) 17-20
ond, we have now a possibility to do recurrences prior to the integration of differential equations of satellite motion. After this, in the satellite motion theory we evaluate the elliptic Hansen coefficients using their "Chebyshev approximations". It allows us to save computer memory and to decrease the time of computations. The above described scheme of computation of elliptic Hansen coefficients can be used also for computation of derivatives of these coefficients which are included into equations of motion, too.
4. Some important remarks Let us make now some remarks concerning the accuracy of computing the elliptic Hansen coefficients in the satellite motion theory. For constructing analytical and semianalytical theories of motion for artificial Earth satellites we have to evaluate the coefficients with an accuracy higher than the order of terms eliminated from consideration. This means that for theories of the first order we should evaluate the elliptic Hansen coefficients with the relative accuracy 10 -4, because the perturbations of the second order are approximately 103 times smaller (in their absolute values) than the perturbations of the first order [ 12]. For theories of the second order we should compute the first order perturbations with errors smaller than 10 -7 and the second order perturbations with errors smaller than 10 -4 . Then the errors of the theory will be smaller than the perturbations of the third order. Because the analytical and semianalytical theories are most efficient for prediction of motion of satellites with not a very high accuracy but on a large interval of time, then it is usually enough to take into consideration the perturbations of the first order only [ 13] and only in special cases should we consider the perturbations of the second order. In order to reach the accuracy 10 -7 the approximation of elliptic Hansen coefficients should be made with the help of the Chebyshev approximation up to and including the polynomials of the sixth order, if the interval of approximation is [e0 - 0.05; e0 + 0.05] and the initial value of eccentricity e0 is placed in the domain 0.45 < e0 < 0.85. The interval of eccentricity [e0 - 0.05; e0 + 0.05] is usually enough, if we predict the motion on such intervals of time as one or two
19
months or less, because the eccentricity has no secular perturbations. But if we would like to study the evolution of the orbit within, for example, several years, then we can divide the interval of all possible values of eccentricity by separate fragments and realize the Chebyshev approximation scheme on each interval independently. The summary of our experiments is as follows: ( 1) for computing elliptic Hansen coefficients it is enough to use the Chebyshev approximation of the sixth (maximal, for the emergency, seventh) order; (2) the interval of eccentricity [ e0 - 0.05; e0 + 0.05 ] is enough for doing approximations from the point of view of practical purposes; (3) under conditions (1), (2) we evaluate the elliptic Hansen coefficients with relative errors of 10 -7 which is quite acceptable for all analytical and semianalytical theories of the first and the second order; (4) the above described approach is valid up to the values of eccentricity e < 0.9; there are no limitations concerning critical inclinations, resonances, etc.; (5) this approach gives a possibility to improve the existing analytical theories substantially from the point of view of their compactness (not accuracy). The computer experiments show that the computation of elliptic Hansen coefficients by means of their "Chebyshev approximations" instead of the exact analytical formulae allows one to speed up the process of computing the Hamiltonian (and respectively the right parts of differential equations of motion) by approximately ten times for low harmonics of the geopotential. For high harmonics of the Earth potential the effect of decrease of computing time is more substantial because the exact formulae for elliptic Hansen coefficients are more complicated and more cumbersome for high harmonics but "Chebyshev expansions" do not suffer any complications with the increase of the order of harmonics. Saving computer memory in evaluating elliptic Hansen coefficients and their derivatives gives a possibility of extending the number of geopotential harmonics and the number of terms of the lunar-solar disturbing function which can be taken into consideration.
20
A.A. Vakhidov, N.N. Vasiliev/Computer Physics Communications 118 (1999) 17-20
All our computer investigations have been made on an IBM PC AT computer with an Intel 80486 processor in MS DOS 6.0 operating system. The computer experiments have shown a high efficiency of the described approach and possibility to use it for computation of perturbations in the satellite motion.
Acknowledgements The authors are very grateful to Prof. Andre Deprit for his kind attention to our work and for useful suggestions for improving the text of the manuscript. The authors thank also Prof. V.A. Brumberg for useful discussions.
5. Conclusion References Special computing techniques presented in this paper give the possibility of solving some problems in evaluating elliptic Hansen coefficients and their derivatives from the point of view of saving computer memory and increasing the rapidity of satellite motion computation. On the basis of the described approach one can develop a very efficient satellite motion theory within elliptic functions. At present such a theory is constructed by the authors for an artificial Earth satellite with the eccentricity e = 0.74 [ 14]. The accuracy of this theory is quite acceptable for practical purposes. The results of comparison of this theory with the direct numerical integration of equations of motion in Cartesian coordinates are also presented in Ref. [ 14]. In the future we also plan to apply this theory to other artificial Earth satellites with large eccentricities.
[ 1] E.V. Brumberg, in: Proc. 25th Symp. on Celest. Mech. (Tokyo, 1992) p. 139. [2] E.V. Brumberg, T. Fukushima, Celest. Mech. Dynam. Astron. 60 (1994) 69. 13] N.N. Vasiliev, A.A. Vakhidov, S.A. Klioner, Preprint Institute of Theoretical Astronomy of the Russian Academy of Sciences No. 60 (1996) [in Russian]. [4] S.A. Klioner, A.A. Vakhidov, N.N. Vasiliev, Celest. Mech. Dynam. Astron. 68 (1998) 257. [5] N.N. Vasiliev, A.A. Vakhidov, A.G. Sokolsky, Preprint Institute of Theoretical Astronomy of the Russian Academy of Sciences No. 38 (1994) [in Russian]. [6] A.A. Vakhidov, N.N. Vasiliev, J. Symbolic Comput. 24 (1997) 705. [7] G.E.O. Giacaglia, Celest. Mech. 14 (1976) 515. [8] R.H. Gooding, Celest. Mech. 4 (1971) 91. [9] D.R. Cok, Celest. Mech. 16 (1977) 459. [10] A.G. Sokolsky, A.A. Vakhidov, N.N. Vasiliev, Celest. Mech. Dynam. Astron. 63 (1996) 357. I l l ] B.W. Char, K.O. Geddes, G.H. Gonnet, B.L. Leong, M.B. Monagan, S.M. Watt, MAPLE V Library Reference Manual (Springer, New York, 1993). [12] A.M. Fominov, Bull. ITA 14 (1980) 621 [in Russianl. [13] E.P. Aksenov, Theory of Motion of Artificial Earth Satellites (Nauka, Moscow, 1977) lin Russian]. [ 14] A.A. Vakhidov, N.N. Vasiliev, Astron. J. 112 (1996) 2330.