Development of a New Form of Equations of Disturbed Motion of a Satellite in Nearly Circular Orbits

The use of simple physical reasoning instead of the method of varying constants has made it possible to elaborate a short scheme of deriving equations of disturbed motion of a satellite in nearly circular orbits. The use of a circular Keplerian orbit as a reference orbit has ensured the nondegeneracy of equations and their simple relation with time. All this taken together has made it possible to propose a form of equations convenient for carrying out numerical and analytical studies with its variables having a simple physical meaning.<br><br>A relationship between the introduced variables and the Keplerian elements of the orbits was described for undisturbed motion. It was shown that the variables describing the deviation of the orbit radius from the radius of the reference orbit are proportional to eccentricity and deviation of the focal parameter is proportional to the square of the eccentricity.<br><br>Relationships were constructed that describe the connection between the introduced variables and the Cartesian coordinates of position and velocity in the inertial coordinate system as well as arguments for choosing the radius of the reference orbit. From the condition of equality of energies of motion along circular reference orbit and in elliptical Keplerian orbit, it is expedient to take the radius of the reference orbit equal to the semi-major axis of the Keplerian ellipse.<br><br>Approaches to the possible development of the proposed equations were presented. They make it possible to describe changes in the argument of the orbit perigee. The proposed change of variables makes it possible to avoid degeneracy of equations at very small eccentricities when studying the change in the orbit perigee.<br><br>The advantages of using the proposed equations for numerical and analytical studies of satellite motion in nearly circular orbits were shown on concrete calculation examples. It was shown that the results of numerical integration in the proposed variables give almost five orders of magnitude less error than the results of the integration of equations in Cartesian coordinates.


Introduction
Space Earth missions often use nearly circular orbits for which changes in radius per revolution do not exceed tenths of a percent. In particular, this applies to satellites of remote Earth sensing for which proximity of trajectories to a circular orbit is important for the observation quality [1][2][3].
The disturbed motion equations include knowledge of basic motion laws. These laws are incorporated into equations of disturbed motion by fixing expressions describing laws of undisturbed motion [4]. Simplicity and clarity of the equations which largely determine the simplicity and success of the study depend on the choice of the expressions of undisturbed motion being used. Formulas of elliptical motion in Keplerian orbits are conventionally used in celestial mechanics as the main expressions describing regularities of undisturbed motion [5,6]. This leads to rather cumbersome and inconvenient calculations that determine the relationship between time and parameters of motion. As will be shown below, the use of a circular Keplerian orbit as a reference orbit and choice of "good" variables describing motion deviation from the circular orbit makes it possible to avoid these problems.
Derivation of the proposed equations in one or another form was previously described in [4,7,8]. In [4], this conclusion had rather the nature of theoretical reasoning about the possibility of using new approaches to constructing equations of disturbed Keplerian motion. A similar but somewhat different form of equations adapted to the problem was used in [7]. It is clear that equations, as a research tool, are desirable to be selected for a specific problem. Practical work related to calculations and analysis of the motion of low-orbit Earth satellites has shown the convenience of just such a system of equations that is presented in this article. Development of the proposed form of equations required derivation of formulas describing the relationship between the introduced variables and the vectors of position and velocity of the satellite, determination of connection of the introduced variables with the Keplerian elements of orbits. The problems of studying changes in the orbit ellipticity and changes in the perigee argument are associated with further development of the proposed form of equations.
To date, there are a large number of various forms of equations of satellite motion and methods of their study. Despite this, the development of forms of equations describing satellite motion continues so far. This is connected with new challenges emerging at the present stage of space exploration. Derivation of planetary Lagrange equations and the Gaussian form of equations of disturbed motion which are most common among researchers is based on the method of variation of arbitrary constants. This conclusion is unnecessarily cumbersome and complicates the introduction of new variables. Therefore, it is urgent to develop an approach to deriving equations of disturbed motion which would make it possible to eliminate this drawback. For nearly circular satellite orbits, it is important to use circular orbits as reference orbits because they provide a simple relationship of variables with time.

Literature review and problem statement
There are many sets of variables and their corresponding systems of equations that describe satellite orbiting. Among their variety, the following groups can be distinguished [5,[9][10][11]: equations of change of Cartesian coordinates of a satellite (projections of the radius-vector R  of the satellite and its velocity R   on the selected Cartesian coordinate system are considered as variables); equations in osculating elements (different sets of variables that describe the change in Keplerian orbital elements are used). For example, a, the semimajor axis of the orbit (focal parameter p is used sometimes); e, eccentricity; Ω, the longitude of the ascending node; ω, perigee argument; I, the orbit inclination; ν, the true anomaly (or another variable determining the position of the satellite in the orbit at each moment of time) are considered as variables as well as their various combinations; equations of motion of the satellite in a canonical form with construction of various Hamilton functions.
In addition to these variables and equations of their change, there are also equations of motion in cylindrical coordinates, spherical coordinates, polar Hansen coordinates, various canonical forms of equations [11], etc. Such a multitude of possible forms of equations is apparently explained by complexity and variety of the problems the researchers are facing and attempts to find a form of equations suitable for a specific problem using a rich mathematical apparatus.
The need for practical application of equations of disturbed motion encounters difficulties associated with the study of existing models of motion and the choice of suitable equations for specific studies of satellite dynamics [12].
Apparently, numerical integration of equations of motion in rectangular Cartesian coordinates is the simplest method for constructing a solution to the disturbed satellite motion. In celestial mechanics, this method is known as Cowell's method [14]. However, a series of features of equations of mo-tion in rectangular coordinates reduces the efficiency of this method. These features, in particular, include the Lyapunov instability of equations and the need to carry out computational operations with large quantities which are known to be burdened with large errors of rounding-off.
To improve the numerical construction of solutions in the problem of the disturbed motion of two bodies, various methods were developed. They imply introducing new variables describing motion [15]. A new method based on the introduction of new variables of motion was developed in [12,16,17] for the numerical solution of the problem of the disturbed motion of two bodies. There is also a detailed review of the history of methods for solving the problem of the disturbed motion of two bodies. New variables were introduced in [12,16,17] according to the same scheme as suggested in [18]. Equations of change in orientation of the orbital coordinate system (OCS) were derived on the basis of the theorem on change in kinetic momentum. The expressions obtained for the projections of angular velocities make it possible to introduce arbitrary kinematic parameters describing the OCS orientation. The introduction of variables describing a change in the orbit radius is based on the scalar equation of change in the orbit radius and the integral of energy. As stated in these publications, variables introduced in [12,16,17] significantly improve the efficiency of the numerical construction of problem solutions. However, the obtained equations of motion and their derivation are rather cumbersome and the variables themselves do not have physical clarity.
The applied form of equations of disturbed motion is largely determined by the problem. For example, in [19], the possibility of effective application of vector orbital elements of Milankovitch type for constructing a century-long motion of a satellite under the influence of solar pressure was shown.
In many problems of studying satellite motion, high accuracy of calculations is not required. For example, in the problems of developing systems for removal of objects of space debris [20][21][22], assessment of the effectiveness of the systems being developed involves assessment of the time of satellite removal. In problems of satellite attitude control [23-25], a high-precision long-term forecast of motion of its center of mass is not required as well. Understanding of physics of the process of orbital motion control is required above all in the problems of maintaining target orbits by a satellite [26,27]. In such problems, simplicity of deriving equations of motion of the satellite itself including physical clarity (simplicity of physical interpretation) of the variables describing motion is an important factor. From this point of view, planetary Lagrange and Gaussian equations of disturbed motion are the most common forms of equations [5]. However, classical derivation of these equations based on the variation of constants is very cumbersome. The application of the inference scheme proposed in [18] makes it possible to simplify it significantly.
The Lagrange and Gauss equations of disturbed motion of the satellite degenerate for orbits with an eccentricity close to zero. There are conventional methods of solving the problem of the singularity of these equations based on the change of variables but they lead to complications of equations. It is argued in [28] that since the inverse transformation has never been carried out, the problem of the singularity of the Lagrange and Gauss equations has not been solved. A solution to this problem was proposed in [28]. At the same time, it is obvious from a physical point of view that even a small perturbation can lead to a very rapid change in the perigee argument at a very small orbit eccentricity.
Another drawback of the Lagrange and Gauss equations for nearly circular orbits is the presence in them of a fast variable, the mean anomaly, and the need to calculate complex expressions giving its connection with the true anomaly. These disadvantages are easily overcome when a circular Keplerian orbit is used as an undisturbed orbit.
Thus, it seems important to have a simple form of equations of motion for many problems of studying satellite motion in nearly circular orbits. This form of equations should enable a sufficiently efficient numerical construction of solutions on the one hand and have physically clear variables for analyzing the motion laws on the other hand. Taking into account the variety of problems in the study of satellite motion, the derivation of equations of disturbed motion should be simple and enable the introduction of new variables most suitable for the problem under study.

The aim and objectives of the study
The study objective is to develop a new form of equations of disturbed motion of a satellite in nearly circular orbits convenient for calculations and analytical motion studies.
To achieve the objective, the following tasks were set: -to present a short scheme of the derivation of equations of disturbed motion with the introduction of the variables ensuring nondegeneracy of the equations and their simple relationship with time. The introduced variables should have a simple physical meaning and be convenient for numerical and analytical studies; to determine the relationship between the considered osculating elements and components of the radius-vector R  and the satellite velocity R   ; -to describe the relationship of the introduced variables with the Keplerian elements of orbits; to show the ways of possible development (transformation) of equations depending on the study task; to evaluate the convenience of using the proposed equations for numerical and analytical studies of satellite motion using concrete examples.

1. Scheme of deriving the equations of disturbed motion
The proposed equations of disturbed motion of a satellite in nearly circular orbits are based on the theorem of change in angular momentum and fixing the Kepler circular orbit as an undisturbed reference orbit. Two right-hand coordinate systems are used: the inertial coordinate system (ICS) OXYZ (the OXY plane lies in the equatorial plane; the OX axis is directed to the spring point; the OZ axis is directed along the Earth's rotation axis to the North Pole); the orbital coordinate system (OCS) Oxyz (the Ox axis is directed along R  ; the Oy axis is in the plane of the instantaneous orbit and lying in the direction of the satellite motion; the Oz axis complements the system to the right rectangular one).
The orientation of Oxyz in OXYZ is described conventionally: by the Euler angles i, Ω, u, that is an inclination, longitude of the ascending node, and argument of latitude, respectively.
Derivation of equations follows the conclusion of [8] and is given here for completeness of presentation.
The following equation is original: where R  is the radius-vector of the satellite center of mass relative to the Earth center; μ is the gravitational parameter of the Earth; F  is the disturbing acceleration; the dots denote time derivatives.
Multiply the equation of motion vectorially by R  from the left and proceed to the differentiation of the vector of specific angular momentum L  ( )  in the OCS to obtain the following: where the prime denotes the relative derivative L  in the OCS; ω  is the angular velocity of the OCS rotation relative to the is the specific moment of disturbing forces. (2) on the OCS axis gives the following: 2 1 , Using the known kinematic relations for the Euler angles, the following is obtained: 3), and the variable p (focal parameter of the orbit) is introduced from equality = µ .
L p To construct a group of equations describing the change in the orbit shape, a circular orbit is fixed as an undisturbed motion, that is, deviations of R and p from a constant magnitude are considered in the disturbed motion (radius is a constant value in the circular Keplerian orbit). To describe these deviations, three variables b 1 , b 2 , γ are introduced which are small in the disturbed motion along nearly circular orbits.
The variables b 1 , b 2 , γ are introduced by the following relations: where R 0 is the radius of the undisturbed circular orbit; b 1 , γ are deviations of the current radius and focal parameter of the disturbed orbit from the radius of the undisturbed orbit, respectively, referred to R 0 ; b 2 is the radial velocity in the disturbed orbit referred to as the velocity of motion in the undisturbed circular orbit. Differentiation of equalities (5) gives the following: Taking into account the equation of change in R obtained by projecting (1) onto the axis Ox, it is easy to obtain from equalities (5) where the following designations are introduced: z=1+b; s=1+γ; Δu=u-u 0 ; notes the derivative with respect to u 0 . Note that 2 0 R µ is the acceleration of gravity for a given height (R 0 ). The quantity z is the dimensionless radius of the orbit equal to the ratio of the orbit radius to the radius of the reference orbit. The quantity s is a dimensionless focal parameter of the orbit equal to the ratio of the focal parameter of the orbit to the focal parameter of the reference orbit (since the reference orbit is circular, its focal parameter is equal to R 0 ). Note

2. Determination of the values of osculating elements from velocity and radius. Selection of the reference orbit radius
The numerical study includes verification of the results obtained by comparison with previously known results or the results of observation of satellites that are publicly available on the Internet. Return of satellite coordinates and velocities in the ICS is often the most widely used form of displaying results in such resources.
Let us consider the problem of determining the osculating elements of the orbit i, Ω, γ, u, b 1 , b 2 from known R  and , R   given by projections in the ICS.
The vector L R R = ×     and the unit vector 3 / , | | e L L L L = =    set the basis vectors of the orbital coordinate system: The focal parameter p was introduced from equality , L p = µ therefore, p=L 2 /μ. Then  There remains arbitrariness with the choice of R 0 as a radius of the circular reference orbit. Generally speaking, equations (7) are true for any R 0 but the convenience of using these equations presupposes smallness of γ, b 1 , b 2 . Therefore, the only formal requirement for the choice of R 0 is the requirement that the initial values of γ, b 1 , b 2 be small.
In the case when R  and R   are such that the corresponding Keplerian orbit has a very small eccentricity (e<0.001). Then changes in motion in a low near-earth orbit caused by the difference between the Keplerian orbit from the circular one (eccentric oscillations) are comparable in magnitude with the changes caused by the disturbing effects of the second zonal harmonic of the Earth's gravitational potential. In this case, the ellipticity of the disturbed orbit is of conditional nature and R 0 can (conveniently) be taken equal to the focal parameter of the undisturbed orbit, 2 0 н / . R p L = = µ When the eccentricity of the Keplerian orbit is significant (e>0.001), it is advisable to take a value of the semi-major axis a of the Keplerian ellipse as R 0 . In this case, average motion along a circular orbit of radius R 0 will be equal to the average motion in the Kepler orbit. Another way to say: R 0 is chosen so that energy of motion in the circular reference orbit is equal to the energy of motion in an elliptical Keplerian orbit: The transition from osculating elements to R  and R   in the ICS presents no difficulties and is not presented here.
The relationship between the proposed variables and components of vectors , R  R   makes it possible to compare the calculation results with the results obtained using known software products and sets of other variables describing the satellite orbiting dynamics.

3. Relationship with Keplerian elements in undisturbed motion
Since In undisturbed motion, change in the elements is described by equations: Write the change b 1 in the form where γ=const. Equation (13) has a first integral (energy integral).
where h b is a constant. The following is found from the last equality: For the considered motion, h b <0, z oscillates periodically between roots of the equation There is a single root z=1 at h b =-1/2(s=1). Equation (15) is solved in quadrature where z 1 , z 2 are the roots of equation (16), z 1 <z 2 . The integral in (17) is a tabular integral. It can be calculated but the resulting expression is too cumbersome and it seems very difficult to express z explicitly as a function of u 0 . But this does not really matter since b 1 <<1 for the considered case of motion and equation (13) is close to the equation of a linear oscillator.

Development of forms of equations of disturbed motion
In a series of problems, it is necessary to monitor the position of the orbital perigee (apogee). The variables introduced above do not directly describe the position of the orbit perigee. Approaches to the possible development of a system of equations (7) are proposed below. They make it possible to describe changes in the perigee argument.
Deviation of the orbit radius from the reference orbit radius R 0 is described by the following equations or the equation Since, γ, b 1 , b 2 <<1, then the following is obtained by referring quantities of the second order of smallness to the disturbing accelerations: where ( )  Let us introduce variables A, α conventional for the averaging method as follows: Here, roughly speaking, A is the eccentricity of the orbit and α is the apogee argument. Then where ( ) (22) have a singularity when A is close to zero. To avoid this singularity, introduce variables τ 1 , τ 2 : τ 1 =Acosα, τ 2 =Asinα conventional for celestial mechanics.
and the equations for τ 1 , τ 2 will take this form: The application of formulas (21) to (23) assumes that only the periodic component is retained in the changes of b 1 .

5. Examples of calculations
Let us compare the results of calculations of the satellite's orbital elements obtained using the proposed equations and equations of motion in Cartesian coordinates. Fig. 1 shows the change in radius of a nearly circular (e=0.0001) undisturbed orbit with an altitude of 300 km relative to the mean radius of the Earth in a time interval of 1 day calculated using equations of motion in Cartesian coordinates (R xyz ), the system of equations (7) (R oe ) and analytical ex- The change in the true anomaly ν over time was calculated as follows [5]: where ( ) The calculations were carried out with a step of return equal to 1 second.
As seen in Fig. 1-3, the solutions are very close. In such an integration interval (about 16 orbital turns), the error δR xyz does not exceed 1.8·10 -5 %, however, the error δR oe is five orders of magnitude less and does not exceed 4.3·10 -10 %. This difference in errors is explained by the fact that the right-hand sides of the equations in Cartesian coordinates contain quantities that are orders of magnitude larger than the right-hand sides of the proposed equations. The small-ness of the right-hand sides in (7) makes it possible to use larger values of the integration steps in comparison with the use of equations in Cartesian coordinates. Fig. 4-6 show the results of calculations for initial conditions the same as in the previous case taking into account the effect of the second zonal harmonic of decomposition of the geopotential among spherical functions [5]. Fig. 4, similar to Fig. 1, shows the change in R xyz and R oe and a linear approximation of the change in the radius R lin obtained on the basis of (7). The solution of ( ) was obtained according to the procedure presented in [8]. In the case under consideration, the initial value of u 0st =0, then, integrating the equation for the change in b 1 , the following is obtained: 2  It can be seen from Fig. 5, 6 that the solutions are in good agreement and the line ar analytical approximation (25) in the integration interval of about 15 orbital turns gives an approximation no worse than 4·10 -3 %.

Discussion of the results obtained in the development of equations of disturbed satellite motion in nearly circular orbits
The use of a circular Keplerian orbit as a reference orbit has made it possible to develop a simple and clear form of equations of motion of satellites in nearly circular orbits (7). In contrast to the conventional equations of disturbed Keplerian motion [5,6], the proposed system of equations (7) contains only slow variables, and the relationship between their independent variable with time is described by the simplest expression. The use of simple physical considerations, rather than the method of variation of constants, has made it possible to introduce physically clear variables describing the deviation of the trajectory from the reference orbit. The above derivation of equations of disturbed motion leaves many opportunities for introducing new variables convenient for a concrete problem. The proposed form of equations is also valid for orbits with large eccentricities. However, the efficiency of its use significantly decreases with an increase in orbit eccentricity as soon as the right-hand sides in (7) cease to be small quantities.
The fixing of the circular reference orbit has allowed us to obtain the simplest relationship between the independent variable of the equations and time. However, this fixation of the radius of the circular reference orbit causes difficulties in studying the satellite motion in cases when the radius of its orbit changes significantly. For example, such a case occurs when solving problems of satellite removal. It seems that introduction of intermediate reference orbits will make it possible to avoid these difficulties without a significant loss of efficiency. However, the algorithms of introducing intermediate reference orbits and algorithms of recalculating the initial data of motions require additional studies.
The proposed form of equations has shown the possibility of its effective use in elaboration of approximate analytical solutions for the disturbed motion of satellites in low near-earth orbits [8]. It seems promising to use it for constructing an approximate analytical model of such a motion of satellites. This will require the construction of a second approximation in terms of small parameters for solving the equations. Taking into account the aerodynamic resistance as an essential factor in the motion of satellites in low nearearth orbits can cause a particular difficulty.

Conclusions
1. The use of simple physical reason instead of the method of variation of constants has made it possible to develop a short scheme of deriving equations of disturbed motion. The use of a circular Keplerian orbit as a reference orbit has ensured the nondegeneracy of equations and their simple connection with time. All this taken together has made it possible to propose a form of equations convenient for nu-merical and analytical studies with its variables having a simple physical meaning.
2. The solution to the problem of introducing new variables presupposes a description of their connection with the known variables. Connections between the introduced variables and the Keplerian elements of the orbit was described for undisturbed motion. It was shown that the variables describing the deviation of the orbit radius from the radius of the reference orbit are proportional to eccentricity and deviations of the focal parameter are proportional to the square of eccentricity.
3. The use of new forms of equations presupposes the connection of their variables with the Cartesian coordinates of position and velocity in the inertial coordinate system. Relationships describing this connection are given and arguments for choosing the radius of the reference orbit are given. Proceeding from the condition of equality of energies of motion in a circular reference orbit and the elliptical Keplerian orbit, it is expedient to take the radius of the reference orbit equal to the semi-major axis of the Keplerian ellipse.
4. The introduced new variables do not directly describe perigee position in the orbit. Approaches to the possible development of the proposed equations were presented. They make it possible to describe changes in the argument of the orbit perigee. The proposed change of variables makes it possible to avoid degeneracy of equations at very small eccentricities when studying the change in the orbit perigee.
5. The advantages of using the proposed equations for numerical studies were demonstrated on specific design examples in comparison with the integration of the equations of motion in Cartesian coordinates. It was shown that the results of numerical integration in the proposed variables give almost five orders of magnitude less error compared to the results of the integration of equations in Cartesian coordinates. The advantages of the proposed equations for analytical studies of satellite motion in nearly circular orbits were also shown. It has appeared that the linear analytical approximation gives an approximation no worse than 4-10 -3 % when taking into account the influence of the second zonal harmonic of geopotential in an integration interval of about 15 orbital turns.