Quaternion methods and models of regular celestial mechanics and astrodynamics

This paper is a review, which focuses on our work, while including an analysis of many works of other researchers in the field of quaternionic regularization. The regular quaternion models of celestial mechanics and astrodynamics in the Kustaanheimo-Stiefel (KS) variables and Euler (Rodrigues-Hamilton) parameters are analyzed. These models are derived by the quaternion methods of mechanics and are based on the differential equations of the perturbed spatial two-body problem and the perturbed spatial central motion of a point particle. This paper also covers some applications of these models. Stiefel and Scheifele are known to have doubted that quaternions and quaternion matrices can be used efficiently to regularize the equations of celestial mechanics. However, the author of this paper and other researchers refuted this point of view and showed that the quaternion approach actually leads to efficient solutions for regularizing the equations of celestial mechanics and astrodynamics. This paper presents convenient geometric and kinematic interpretations of the KS transformation and the KS bilinear relation proposed by the present author. More general (compared with the KS equations) quaternion regular equations of the perturbed spatial two-body problem in the KS variables are presented. These equations are derived with the assumption that the KS bilinear relation was not satisfied. The main stages of the quaternion theory of regularizing the vector differential equation of the perturbed central motion of a point particle are presented, together with regular equations in the KS variables and Euler parameters, derived by the aforementioned theory. We also present the derivation of regular quaternion equations of the perturbed spatial two-body problem in the Levi-Civita variables and the Euler parameters, developed by the ideal rectangular Hansen coordinates and the orientation quaternion of the ideal coordinate frame. This paper also gives new results using quaternionic methods in the perturbed spatial restricted three-body problem.


Introduction
Celestial mechanics and astrodynamics are based on Newtonian differential equations of the perturbed spatial two-body problem and the perturbed spatial restricted three-body problem. Newtonian equations of the perturbed spatial two-body problem degenerate when the second body (body of interest) collides with the first (central) body, i.e., when the distance between the bodies is zero. This makes these equations inconvenient for studying the motion of the second body in the near neighborhood of the central body or its motion along highly eccentric orbits. Singularity in the coordinate origin creates not only theoretical, but also practical (computational) difficulties in the two-body problem. Analogously, Newtonian equations of the perturbed spatial restricted three-body problem degenerate when the body of interest (with negligible mass) collides with one of the other two bodies, which have finite masses. This makes these equations inconvenient for studying the motion of the negligible mass body in the neighborhood of the first or second gravitating body, and creates not only theoretical but also practical (computational) difficulties in the restricted three-body problem. The elimination of singularities (division by zero) induced by gravitational forces in classical (Newtonian) equations of celestial mechanics and astrodynamics is known as regularization. Equations without such singularities are called regular. Among regularization methods and regular models of celestial mechanics and astrodynamics, the use of quaternion methods and models based on hypercomplex variables, i.e., the Hamilton quaternions, which use Kustaanheimo-Stiefel (KS) variables or Euler parameters as their components (elements), has recently become widespread. These methods and models offer a number of analytical and computational advantages.
Stiefel and Scheifele [37] , Bordovitsyna [38] , Bordovitsyna and Avdyushev [39] , Fukushima [40][41] , Pelaez et al. [42] , Bau et al. [43] , Amato et al. [44] , and Bau and Roa [45] demonstrated the results of comparing the numerical solutions to the equations of orbital motion of celestial and cosmic bodies in the KS variables, Euler parameters, and other variables. These results prove the efficiency of the KS variables and Euler parameters in celestial mechanics and astrodynamics.
Loginov and Chelnokov compared the accuracy of numerical integration of the classical Newtonian differential equations of the spatial restricted three-body problem (the Earth, the Moon, and a spacecraft) in the Cartesian coordinates and the regular quaternion differential equations in the four-dimensional KS variables. The equations were integrated by using the standard Runge-Kutta method of fourth-order accuracy. The regular quaternion equations in the KS variables demonstrated much higher accuracy than those in Cartesian coordinates. The accuracy was higher by 2 orders of magnitude for the circular orbit, by 4 orders for the perturbed elliptic orbits with medium eccentricity, and by 7 orders for the perturbed elliptic orbit with high eccentricity.
This work analyzes the regular quaternion methods and models of celestial mechanics and astrodynamics in the KS variables and Euler parameters. These methods and models are based on the differential equations of the perturbed spatial two-body problem. This paper also covers some applications of these methods and models.
In Section 2, preliminary works by Euler [46] and Levi-Civita [47][48][49][50] on regularizing the twobody problem equations are introduced, in which the solutions are provided for one-dimensional and two-dimensional problems involving the collision of two bodies. The KS regular equations of the perturbed spatial two-body problem in four-dimensional KS variables in scalar and matrix forms were presented by Kustaanheimo [51] , Kustaanheimo and Stiefel [52] , and Stiefel and Scheifele [37] . Primary analytical and computational advantages of these equations were described by Stiefel and Scheifele [37] , Brumberg [53] , Bordovitsyna [38] , Bordovitsyna and Avdyushev [39] , and Fukushima [40][41] . Section 2 also describes regular equations for the twobody problem in Cartesian coordinates by Bohlin [54] , in which the advantages and disadvantages are pointed out, and the development of Bohlin's idea in the works by Burdet [55][56] is also mentioned.
Sections 3 and 4 describe and analyze the quaternion methods and models proposed by the present author in the beginning of the 1980s and 1990s [20][21][22][23][24][25][26][27] , which, based on the citations by publications in this field, are mostly unknown to English-speaking researchers.
Section 3 describes and analyzes the results obtained by the present author in the field of quaternion regularization for the perturbed spatial two-body problem in the 1980s [20][21][22][23] . Soon after the discovery of the KS transformation, the usage of Hamilton quaternions (fourdimensional hypercomplex numbers) and four-dimensional quaternion matrices for regularizing the equations of the spatial two-body problem had been examined. However, Stiefel and Scheifele totally rejected this idea and wrote in their book [37] that "Any attempt to substitute the theory of the KS matrix by the more popular theory of the quaternion matrices leads to failure or at least to a very unwieldy formalism. " We refuted this statement, apparently for the first time, in Refs. [20] and [21], by demonstrating that actually, the quaternion approach to regularization reveals geometric and kinematic essence of the KS transformation. This allows the regular equations in the KS variables to be derived directly and conveniently, makes the fundamental principles of the KS regularization more intuitive and elegant, and supports the derivation of the more general quaternion regular equations of the perturbed spatial two-body problem.
Regarding the above-cited statement by Stiefel and Scheifele on quaternion matrices, i.e., having no chance of success in the regularization theory, Waldvogel [10] wrote "This statement was first refuted by Chelnokov (1981) who presented a regularization theory of the spatial Kepler problem using geometrical considerations in a rotating coordinate system and quaternion matrices. In a series of papers, including Chelnokov (1992Chelnokov ( , 1999, the same author extended the theory of quaternion regularization and also presented practical applications." Section 3 provides convenient geometric and kinematic interpretations of the KS transformation and the KS bilinear relation. This relation, according to Stiefel and Scheifele [37] , is essential to their formulation of regular celestial mechanics. These interpretations were offered by Chelnokov [20][21] , who showed that the regularizing KS transformation of the coordinates consists of transitioning from the Cartesian coordinates of the second body to the normed Euler parameters (components of the Hamilton rotation quaternion). These parameters, which in Russia are usually called Rodrigues-Hamilton parameters, characterize the orientation of some coordinate frame rotating in the inertial space. Therefore, the KS regular equations can be derived from the original Newtonian two-body problem equations by rewriting them in the rotating coordinate frame and using the Euler parameters for defining the orientation of that coordinate frame. Such a derivation, which uses quaternion matrices and Hamilton quaternions, was proposed by Chelnokov [20][21] . More general quaternion regular equations (compared with the KS equations) of the perturbed spatial two-body problem in the KS variables were derived in matrix and quaternion forms, with the assumption that the KS bilinear relation is not satisfied. These equations are presented in Section 3. Section 3 describes the regular equations of the perturbed spatial two-body problem in quaternion osculating (i.e., slowly changing) elements, corresponding to the KS variables and their first derivatives with respect to fictitious time. These equations were derived by Chelnokov [28,30] from quaternion regular equations in the KS variables by the method of variation of arbitrary constants of integration. This section also covers regular equations of the perturbed spatial two-body problem derived by Chelnokov [21] for the case when the essence of transformations of the original equations of this problem stays the same, but the KS variables are not used. These regular equations, which contain the quaternion kinematic equation in the Euler parameters (of the first order), are not linear in the case of Keplerian motion, but have a simple and convenient form.
Chelnokov [22,[26][27] used the ideas on quaternion regularization of two-body problem equations to derive the quaternion theory for regularizing the vector differential equation of the perturbed motion of a point particle in the central force field with the potential Π(r), which is assumed to be a random differentiable function of the distance r from the point particle to the center of that force field, under the influence of the perturbing force with the potential Π * (t, r) and the perturbing acceleration p(t, r, dr/dt) (r is the radius vector of a point particle, which directs from the center of the force field).
The results derived in these works are presented in Section 4. This section presents the following aspects: the derived general quaternion differential equations of perturbed central motion of a point particle with three regularizing functions (Subsection 4.1); conditions for transforming the quaternion equations of perturbed central motion to oscillator form (to the form of equations of a four-dimensional perturbed oscillator; in the case of unperturbed central motion, this performs harmonic oscillations with constant frequency), which is convenient for analytical and numerical study (Subsection 4.2); quaternion equations of perturbed central motion in normal form (Subsection 4.3); systems of regular quaternion equations of perturbed central motion, which use the KS variables or Euler parameters (Subsection 4.4); and systems of quaternion equations of perturbed motion, which contain the generalized Binet equation (Subsection 4.5).
The presented quaternion systems of equations of perturbed central motion differ in structure, order, and of dependent and independent variables. We provide the comparison of these systems and show their properties and fields of use.
In particular, the main advantage of the differential equation systems of perturbed central motion, which are derived using the Euler parameters and an independent variable τ (dτ = r −2 dt) or an independent variable ϕ (dϕ = cr −2 dt), is that in these systems, every second-order quaternion differential equation in Euler parameters is regular for perturbed motion of a point particle in a central force field with any arbitrary potential Π(r), where c is the modulus of moment of momentum of a point particle. Furthermore, in the case of unperturbed central motion, each of these quaternion equations in the Euler parameters becomes equivalent to the equation of motion of a four-dimensional single-frequency harmonic oscillator with the frequency c/2 or 1/2, which is convenient for solving certain problems. In these equation systems, the equations for the total energy and the modulus of moment of momentum vector of a point particle are also regular for any potential Π(r). Equations for the distance r, on the other hand, are regular only if the potential Π(r) has the fourth-order item relative to the value r −1 , which is a reciprocal for the distance to the center of attraction, i.e., if Π(r) = −a 1 r −1 − a 2 r −2 − a 3 r −3 − a 4 r −4 , where a i = const.
These systems of regular equations can be used, in particular, to study point particle motion in the curved spacetime described by the Schwarzschild metric. Point particle trajectories in such a space match the trajectories of motion in a central force field with the potential Π(r) = −a 1 r −1 − a 3 r −3 , where a i = const. Therefore, these regular equations can be used to predict the motion of planets considering the effects of the general theory of relativity (GTR). We have used them to derive regular quaternion equations of perturbed orbital motion of a solid body in the Earth's gravitational field, considering its zonal, tesseral, and sectorial harmonics (with regularization of equation terms which contain negative power of distance r up to and including fourth order) [57] .
In Subsection 4.6, the efficiency of using the modified four-dimensional regular variables, proposed by Chelnokov [22,25] , instead of the KS variables, is demonstrated using the example of equations of motion of a satellite in the Earth's gravitational field considering its zonal harmonics.
In Subsection 4.7, the equations for unperturbed central motion and for relating the Euler parameters and the KS variables to the elements of an orbit [23,29] are presented. In Subsection 4.8, the solution to the spatial problem of unperturbed central motion in a uniformized form is presented, which eliminates the necessity to consider branching of solutions around critical points, such as poles [23,29] .
In Subsection 4.9, the quaternion solution to the orbit orientation problem is provided [25] . These relations form the basis for determining the osculating motion if the perturbed Keplerian motion is defined by differential equations in Euler parameters and variables r and c.
In Subsection 4.10, the equations of perturbed orbital motion of a point particle are presented, which include the equations in quaternion osculating elements corresponding to Euler parameters and their first derivatives. These equations were derived by Chelnokov [25] using the method of variation of constants.
In Subsection 4.11, the solution to the problem of unambiguous estimation of Euler parameters and KS variables using the Cartesian coordinates and their derivatives [23,25] is presented. Note that the equations for calculating the KS variables and their derivatives given in the book by Stiefel and Scheifele [37] are ambiguous.
Section 5 analyzes the works by other authors on the KS and quaternion regularization of differential equations of the two-body problem.
Section 6 presents the derivation of regular quaternion equations of the perturbed spatial two-body problem in Levi-Civita variables and Euler parameters, derived using the ideal rectangular Hansen coordinates and the orientation quaternion of the ideal coordinate frame. In addition to the known advantages of KS regular equations, these equations have their own additional advantages.
Note that Levi-Civita, regarding his attempts to generalize his regularization of equations of a two-dimensional two-body problem to the spatial problem, admitted later [49] "The problem in space has long resisted my efforts, as I tried to approach it by similar coordinate changes." Stiefel and Scheifele [37] wrote in their book that Levi-Civita tried hard to generalize his method of regularizing differential equations of two-dimensional motion in the two-body problem to the general spatial two-body problem, but without success.
Aarseth and Zare [58] stated that "Because of fundamental difficulties originally clarified by Hopf (1931) [59] and Hurwitz (1933) [60] , it is impossible to generalize the Levi-Civita transformation to an equivalent set of three-dimensional variables." Nevertheless, it has been demonstrated by Chelnokov [31] that the Levi-Civita regularization can be applied successfully to derive regular equations of the perturbed spatial two-body problem. We accomplished that by rewriting the differential equations of the perturbed spatial two-body problem in the ideal coordinate frame (after which the spatial motion equations become the equations of two-dimensional motion), using the ideal rectangular Hansen coordinates, the regular Levi-Civita variables, the orientation quaternion of the ideal coordinate frame in the inertial coordinate frame, the Keplerian energy as an additional variable, and a new independent variable. We describe this quaternion regularization of the perturbed spatial two-body problem equations in Section 6. Section 7 briefly describes the development of quaternion regularization of the perturbed spatial two-body problem equations in our previous work, and introduces the applications of quaternion regular models of astrodynamics to the development and research of regular equations of motion of the Earth's satellite in the Earth's gravitational field, considering its zonal, tesseral, and sectorial harmonics; to the solutions to problems of optimal control of orbital motion of a spacecraft using the principle of maximum; and to the solutions to problems of autonomous inertial navigation in space.
Section 8 presents the results we obtained on the recently received quaternion regularization of the perturbed spatial restricted three-body problem equations.
Note that our work described regular quaternion differential equations of the perturbed spatial two-body problem, which are free from singularities (division by zero) induced by grav-itational forces. Other singularities are caused by using angular variables, in particular, Euler angles, in classical equations of celestial mechanics and astrodynamics. We also obtained new regular quaternion equations of the perturbed spatial restricted three-body problem, written in various rotating coordinate frames with the Euler parameters and free from singularities caused by angular variables. In addition, we obtained regular quaternion equations of the perturbed spatial restricted three-body problem, free from singularities induced by gravitational forces obtained from the mentioned equations in the Euler parameters.
The conclusion sums up our research in the field of quaternion regularization of the perturbed spatial two-body problem equations.

KS regularization of the perturbed spatial two-body problem equations
The differential equation of the perturbed spatial two-body problem in a vector form is where r is the radius vector of the second body's (body of interest) center of mass, directing from the first (central) body's center of mass, r = |r| is the distance between the bodies, µ is a product of gravitational constant by the sum of masses of the bodies, and p is the perturbing acceleration (a function of time t, the radius vector r, and the velocity vector v = dr/dt of the second body in the coordinate frame OXY Z, the origin of which is in the first body's center of mass, and the coordinate axes of which are parallel to the axes of inertial coordinate frame).
The problem of eliminating the singularities in Newtonian differential equations of the perturbed spatial two-body problem is known in celestial mechanics and astrodynamics as the problem of regularizing the perturbed two-body problem differential equations. This problem dates back to Euler [46] and Levi-Civita [49] , who solved one-and two-dimensional problems of collision of two bodies (for the cases of straight-line motion and two-dimensional motion). An efficient regularization of spatial two-body problem equations, the so-called spinor or KS regularization, was proposed by Kustaanheimo [51] and Kustaanheimo and Stiefel [52] . It is a generalization of the Levi-Civita regularization of two-dimensional motion equations, and it was described to the fullest extent in a widely known monograph by Stiefel and Scheifele [37] . The KS regularization takes advantage of the spinor theory, that is, instead of one complex variable from the Levi-Civita theory, it uses two complex numbers. The generalized Levi-Civita matrix, called the KS matrix, was also used. We denote it as L(u KS ). It is a four-dimensional square matrix, and its upper left corner contains the two-dimensional square Levi-Civita matrix as follows: where u j (j = 0, 1, 2, 3) are the KS variables. Instead of denoting one of the KS variables as u 4 , we denote it as u 0 here and further on. The matrix form of the KS transformation is as follows: where x, y, and z are the Cartesian coordinates of the center of mass of the body in interest in the coordinate frame OXY Z.
The scalar form of the relations of the Cartesian coordinates with the KS variables is as follows: It matches the Hopf map [59] except for permutation of subscripts. The KS equations in the scalar form are as follows: where τ is the new independent variable, called the fictitious time, related to time t by the differential relation dt = rdτ (transformation of time proposed by Sundman [61][62] ), h is an additional variable, which has the meaning of Keplerian energy and is defined by the relation where f is the gravitational constant, m and M are masses of the second and first bodies, respectively, and p x , p y , and p z are the projections of the perturbing acceleration p of the second body's center of mass on the axes of coordinate frame OXY Z, which match its projections on the axes of inertial coordinate frame. These equations form a system of ten ordinary nonlinear, nonstationary in a general case, differential equations for four KS variables u j (j = 0, 1, 2, 3), their first derivatives du j /dτ with respect to the new independent variable τ , and for the Keplerian energy h and time t.
They are equivalent to the matrix equation where u KS is the four-dimensional column vector of KS variables, u KS = (u 1 , u 2 , u 3 , u 0 ), and P KS is the four-dimensional column vector of projections p x , p y , and p z of the vector of perturbing acceleration p, P KS = (p x , p y , p z , 0). Let us list the main well-known advantages of the KS equations. Unlike the Newtonian equations, they are regular at the center of attraction, they are linear for the unperturbed Keplerian motion (for the elliptic Keplerian motion, when the Keplerian energy h = const. (h < 0), these equations in variables u j are equivalent to the motion equations of the four-dimensional single-frequency harmonic oscillator in time τ , the squared frequency of which equals the halved Keplerian energy taken with a negative sign), they allow to develop a unified approach to studying all three types of Keplerian motion, they are close to linear equations for the perturbed Keplerian motion, and they allow to present the right-hand sides of differential equations for the motion of celestial and cosmic bodies in the polynomial form, which is convenient for solving them with computers.
These circumstances made it possible to develop the efficient methods for finding analytical or numerical solutions to such problems as the study of motion in the vicinity of attracting masses or motion along high-eccentric orbits. These problems are difficult for classical methods, since in the classical models of astrodynamics, there are singularities, e.g., division by zero when passing through the center of attraction and small denominators when a particle passes near the attracting body. In these models, there are significant nonlinearities. However, equations in the KS variables do not have such singularities. Moreover, they are close (in the presence of perturbations) to linear equations. In particular, Stiefel and Scheifele [37] , Bordovitsyna [38] , and Bordovitsyna and Avdyushev [39] showed that regular equations in the KS variables allow to increase the accuracy of numeric solutions for a number of problems in celestial mechanics and astrodynamics (e.g., the problem of the perturbed motion of an Earth's artificial satellite along high-eccentricity orbits) by three to five orders compared with the solutions obtained with the classical (Newtonian) equations. Fukushima [40][41] also wrote, "The KS regularization resulted in a very efficient integration scheme, improving the accuracy and speed of numerical integration. This comes not only from the structure of the equations, but also from the use of several techniques that bring important advantages to the numerical scheme." Bohlin [54] discovered that the equations of the unperturbed spatial two-body problem (for the Keplerian motion) can be transformed to the following form: where h and A are the energy and the Laplace vector, respectively, which in this case are constants, and the constant α is the length scale, which is introduced to give the new variable τ a physical time dimension.
Therefore, Bohlin proposed the two-body problem equations in the Cartesian coordinates x, y, and z, which for the elliptic Keplerian motion have a form of three-dimensional singlefrequency harmonic oscillator under the influence of acceleration that is constant in modulus and direction in an inertial coordinate frame. The nonhomogeneous part in the analytical solution to the Bohlin equations makes it difficult to derive the equations of the perturbed two-body problem in slowly-changing variables with the method of variation of integration constants.
3 Quaternion regular equations of the perturbed spatial two-body problem The KS regularization is based on the nonlinear ambiguous transformation of Cartesian coordinates of the body in interest (a so-called KS transformation, which generalizes the Levi-Civita transformation). This transformation converts the three-dimensional space of Cartesian coordinates into the four-dimensional space of new coordinates (KS variables). Because of that, Stiefel and Scheifele suggested that direct derivation of regular equations in the threedimensional (i.e., spatial) case is impossible. In their book, Stiefel and Scheifele [37] postulated the matrix regular equation of the spatial two-body problem, which they constructed using the analogy to the matrix regular Levi-Civita equation of two-dimensional motion, and proved with the help of several theorems that the old vector Newtonian equation is satisfied.
Soon after the discovery of the KS transformation, the usage of Hamilton quaternions and four-dimensional quaternion matrices for regularizing the equations of the spatial two-body problem has been examined. However, in their book, Stiefel and Scheifele [37] dismissed that idea. In Chapter 11 of that book, devoted to the geometry of the KS transformation, they wrote, "Any attempt to substitute the theory of the KS matrix by the more popular theory of the quaternion matrices leads to failure or at least to a very unwieldy formalism." This statement has been refuted, apparently for the first time, by the author of this paper [20][21][22][23] that the quaternion approach to regularization actually makes it possible to unfold the geometric and kinematic essence of the KS transformation (including the geometric meaning of its ambiguity), and to derive directly and conveniently the regular equations in the KS variables, makes the fundamental principles of KS regularization more intuitive and elegant, and allows to derive the theory, which generalizes the KS regularization, and to propose the new, more general, quaternion regular equations of the perturbed spatial two-body problem.
Chelnokov [20] derived the quaternion regular equations of the perturbed spatial two-body problem using the four-dimensional quaternion matrix as follows: which consists of the KS variables u j . This matrix is different from the KS matrix. If we interpret the elements u j of this matrix as the Euler parameters λ j , it becomes the well-known quaternion rotation matrix, which describes the rotation in a three-dimensional space.
In the work of Chelnokov [21] , the Hamilton quaternion u proposed by Hamilton [63][64] was used for this purpose. The elements of that quaternion are the KS variables u j , Here, i, j, and k are imaginary units of Hamilton's vector.
The quaternion form of the relations between the Cartesian coordinates x, y, and z and the KS variables is as follows [37] : where an overline denotes the quaternion conjugate, and • stands for quaternion product.
Chelnokov [20][21] showed that regularizing KS transformation of coordinates converts the Cartesian coordinates of the second body's center of mass in the inertial coordinate frame to the new variables that are the specifically normed components of the conjugate quaternion of rotation λ = λ 0 − λ 1 i − λ 2 j − λ 3 k, which characterizes the orientation of a coordinate frame η rotating in the inertial coordinate frame. The axis η 1 of that coordinate frame is directed along the radius vector r of the second body's center of mass. The origin of that coordinate frame is at the center of mass of that body. The norming constant equals the square root of the distance r from the center of mass of the second body to the center of attraction. Therefore, the variables u and λ are connected by the following relation: The bilinear KS relation is which correlates the KS variables u j and their first derivatives with respect to the new independent variable τ , and as Stiefel and Scheifele [37] wrote, plays a central role in their formulation of regular celestial mechanics and puts an additional non-holonomic constraint on the motion of the coordinate system η. This constraint is that the projection ω 1 of the absolute angular velocity vector of the coordinate system η on the direction of radius vector r (axis η 1 ) should be zero, where a dot above a symbol means differentiation with respect to time t. Therefore, transitioning from the Cartesian coordinates of the second body's center of mass to the KS variables in the equations of the spatial two-body problem actually means rewriting these equations in the rotating coordinate frame η, using the four-dimensional Euler parameters λ j , which are components of the rotation quaternion of this coordinate frame, as orientation parameters of the said coordinate frame. Further transformations of these equations have to do with norming the Euler parameters (the rotation quaternion) in the above-mentioned way and transitioning to the KS variables u j , and also with introducing the Keplerian energy h and time t as additional dependent variables, and transitioning to the new independent variable τ .
Note that Chelnokov [20][21] derived the more general (compared with the KS equations) matrix (using quaternion matrices) and quaternion regular equations for the perturbed spatial two-body problem in the KS variables for the case when the above-mentioned KS bilinear relation is not satisfied. These equations are more complex and contain additional terms, with the projections ω 1 and ε 1 of vectors of angular velocity and angular acceleration of the moving coordinate system η on the direction of a radius vector r of the second body's center of mass (one of these projections is an arbitrary parameter).
These equations were derived by using the eight-dimensional screw-motion parameters λ j and λ 0 j (j = 0, 1, 2, 3) of the introduced coordinate frame η, which moves translationally and rotates, and also the matrix [20] and quaternion [21] differential equations of the perturbed motion of a point particle (second body) in the Newtonian gravitational field, derived in these variables. In these works, it was noted that the variables λ j and λ 0 j are the components of the Rodrigues-Hamilton (Euler) dual parameters Λ j = λ j + sλ 0 j (s is the Clifford symbol, s 2 = 0), which, in turn, are the components of the Clifford finite displacement biquaternion [65] , which describes the motion of a moving coordinate frame η, in which the motion equations of a point particle are written. The used equations of point particle motion contain the differential equations in variables λ j and λ 0 j , written in the matrix form using the quaternion matrices [20] , or in the quaternion form [21] using the quaternion variables λ and λ 0 0 . In the work of Chelnokov [21] , it was noted that the used quaternion differential equations in variables λ and λ 0 0 are equivalent to one biquaternion (dual quaternion) kinematic equation [66] , where ω is the angular velocity vector of a moving coordinate frame η, v is the vector of linear velocity of the origin of that coordinate frame, which coincides with the point particle, and ω k and v k are the projections of these vectors on the axes of the moving coordinate frame η. The Cartesian coordinates of the center of mass of a body in interest x, y, and z in the coordinate frame OXY Z are correlated with the variables λ j and λ 0 j by scalar relations, which have the following quaternion form: If in these scalar relations, we assume we can derive the KS formulas as follows: Therefore, it is shown that the KS transformation is a particular case of a more general transformation, which describes the spatial screw motion.
The described scalar relations, which correlate the Cartesian coordinates x, y, and z with the screw motion parameters λ j and λ 0 j , were first derived algebraically by Study [67] as formulas for finding the origin of the new coordinate frame, which was obtained by applying to the old coordinate frame a certain screw motion that is determined by biquaternion operation Λ•(·)•Λ. The derivation of these relations was also described by Kotelnikov [68] with the theory of referencing biquaternions and bivectors to a chosen reference point, and by Chelnokov [66,69] with the theory of finite displacements of a solid body.
Note that the projections η 1 , η 2 , and η 3 of the radius vector r of the second body's center of mass on the axes of a moving coordinate frame η are correlated with the variables λ j and λ 0 j by the following quaternion relation: The regular equations of the perturbed spatial two-body problem in the KS variables of a matrix form are as follows [20] :  Here and elsewhere, a prime symbol (·) means differentiating with respect to the independent variable τ . q j and the variables x, y, and z are defined by the formulas given after the KS equations in a scalar form and after the matrix form of the KS transformation. The regular quaternion equations for this problem in the KS variables are as follows [21] (in this work, these equations were written in a different form, which corresponds to the abovementioned matrix form of these equations): where u, h, and t are unknown functions of the independent variable τ , u is the conjugate for the quaternion u, the projections p x , p y , and p z of the vector p of the perturbing acceleration on the axes of the inertial coordinate frame are specified functions of time t and variables x, y, z,ẋ,ẏ,ż (u j , u j ), ω 1 and ε 1 are the projections of the angular velocity vector ω and the angular acceleration vector ε of the trihedral η on the axis η 1 , respectively, and are specified functions of variables t, u j , and u j , and the central dot stands for scalar product.
The projections of the radius vector r and the velocity vector v =ṙ of the center of mass of the body in interest on the axes of the inertial coordinate frame (the coordinates x, y, z and their time derivativesẋ,ẏ,ż with respect to time t) can be found from the variables u j and u j according to the following relations: Therefore, more general regular equations of the perturbed spatial two-body problem were derived than those by the KS. However, these equations are more complex than the KS equations. Above all, they represent theoretical interest by demonstrating that regularization can be achieved even when ω 1 = 0, i.e., when a bilinear relation, which is central to the KS regularization theory, is not satisfied. It is possible that these equations will be efficient for numerical computations with high-accuracy, because they can be integrated without satisfying the bilinear relation, which is invariably going to be violated during numerical integration due to systematic and computational errors.
By assigning ω 1 = 0 and ε 1 = 0, from Eq. (1), we obtain the quaternion form of the KS regular equations as follows: where scal(u • q) is the scalar part of quaternion u • q. Chelnokov [28,30] also proposed the regular equations of the perturbed spatial two-body problem in quaternion osculating (i.e., slowly changing) elements, derived from the quaternion regular equations in the KS variables by the method of variation of integration constants. These equations are as follows: Here, α and β are the quaternion osculating elements, related to the quaternion variables u and du/dτ * by ν and τ e are the scalar osculating elements, (b, c) is the scalar product of four-dimensional vectors b and c, µ = f (m + M ), and r = (u, u).
The element ν has the meaning of frequency and is related to the energy h by ν = (−h/2) 1/2 , where h < 0, and the time element τ e is related to time t and variables u j , du j /dτ * (j = 0, 3) by where τ * is the new independent variable, which is related to the variable τ , used in the initial equations, by the differential equation, The Cartesian coordinates of the body in interest in the coordinate frame OXY Z and the projections of its velocity can be found using For the unperturbed Keplerian motion, These equations are the quaternion equivalents for the spatial two-body problem equations in regular elements, derived by Stiefel and Scheifele [37] . These equations have the following advantages. First, they are regular (they do not have a singularity at the coordinate origin). Second, in the case of the perturbed elliptic motion, their right parts are slowly and evenly changing functions, and in the case of unperturbed Keplerian motion, the equations are integrated without systematic errors. The disadvantage of these equations is that they are limited to the elliptic motion when h < 0.
Note that Chelnokov [21] showed that regularization of differential equations of the perturbed spatial two-body problem could be achieved when the meaning of the transformations of the initial equations of this problem stayed the same, but the KS variables were not used, that is, the equations were written in the rotating coordinate frame η, the Keplerian energy h was used as an additional variable, and the new time τ , related to real time t by the differential equation dt = rdτ , was used as an independent variable.
The first-order regular equations derived are not linear in the case of Keplerian motion, but have a simple and convenient form as follows [21] : where p 1 , p 2 , and p 3 are the projections of the perturbing acceleration p of the center of mass of the second body on the axes of the rotating coordinate frame η, and the variables are the distance r between the bodies, the projection v 1 = dr/dt of the velocity vector of the body in interest on the radial direction (on the direction of the radius vector r), two projections ω 2 and ω 3 of the vector of angular velocity of the rotating coordinate frame on its own coordinate axes (its first projection ω 1 can be assigned arbitrarily and, in particular, can be prescribed to zero, in which case the equations become significantly simpler), the orientation quaternion λ (the equations include the corresponding quaternion differential equation of orientation of the rotating coordinate frame), the energy h, and time t.
Also note that the algebraic relations (2), which are used to find the projections of the velocity of the second body on the axes of the inertial coordinate frame from the KS variables u j and their first derivatives u j with respect to the new independent variable τ , are singular when r = 0, but the algebraic relations for finding the mentioned velocity projections from variables r, v 1 , ω 2 , ω 3 , and λ j , which are used in the regularization method discussed here, are free from that singularity.
By introducing new variables v 2 = rω 3 and v 3 = −rω 2 into the mentioned regular equations, which have the meaning of the projections of the velocity vector of the second body on the axes of the rotating coordinate frame η, and assuming ω 1 = 0, we can obtain the regular equations of the two-body problem in the different, simpler, and convenient forms [30] ,

Quaternion regular equations of the perturbed central motion
Chelnokov [22,[26][27]29] used the idea of quaternion regularization of the two-body problem equations to develop the quaternion theory of regularizing the vector differential equation of perturbed motion of the point particle M in the central force field, which is as follows: Here, r is the radius vector of the point particle M , directed from the center O of the force field. m is the mass of a point particle. Π is the potential of the central force field, which is assumed to be an arbitrary differentiable function of the distance r from the point particle M to the center O. Π * is the perturbing potential, which is assumed to be an arbitrary function of time t and the position coordinates x, y, and z of the point particle M in the coordinate frame OXY Z, moving translationally relative to the inertial coordinate frame. The perturbing acceleration p is assumed to be an arbitrary function of time t, the radius vector r, and the velocity vector v = dr/dt of the point particle M in the coordinate frame OXY Z, and x, y, and z are the unit vectors of the axes OX, OY , and OZ, respectively. The differential equation of unperturbed central motion of a point particle in the vector form is derived from Eq. (3) when Π * = 0 and p = 0.
In these works, the general quaternion differential equations of perturbed central motion of a point particle with three regularizing functions are derived, and the necessary and sufficient conditions are stated for transforming these equations to the oscillator form (equations of a fourdimensional perturbed oscillator, which, in the case of unperturbed central motion, performs harmonic oscillations with constant frequency), which is convenient for analytical and numerical studies. Various new (including new regular) systems of the quaternion differential equations of perturbed central motion of a point particle in normal and oscillator forms are derived, which differ by structure, dimensionality, and dependent and independent variables used. The comparison of these systems is given, and their properties and fields of use are described.
To derive the regular equations of perturbed central motion of a point particle, the abovestated vector equation was written in the rotating coordinate frame η, the axis η 1 of which was directed along the radius vector r of a point particle. Four-dimensional Euler parameters λ j were used as orientation parameters for that coordinate frame. Equations of motion of a point particle, written in a rotating coordinate frame, include the following differential equations: the scalar equation of the second order for the distance r, two scalar equations of the first order for the projections ω 2 and ω 3 of the angular velocity vector of the coordinate frame η on its own coordinate axes η 2 and η 3 , and the quaternion equation of the first order for the quaternion λ, which describes the orientation of the coordinate frame η in the inertial space and is equivalent to four scalar equations for the Euler parameters λ j . The projection ω 1 of the angular velocity of the coordinate frame η on the direction of the radius vector r of a point particle is an arbitrary parameter, and it is assumed to be zero. After that, the projections c 2 and c 3 of the moment of velocity vector c of a point particle on the axes of the coordinate frame η, defined by c 2 = r 2 ω 2 and c 3 = r 2 ω 3 (projection c 1 = 0), were introduced instead of the variables ω 2 and ω 3 . To derive the regular equations, the differential equations for the Keplerian energy h of a point particle, its total energy h * , and the modulus c of vector c of the moment of velocity of a point particle, or the square of that modulus, were used.
Further regularizing transformations of the mentioned equations, which we have performed, include the following stages: (i) Transition from the differential equations of the first order for variables c 2 , c 3 , and λ to the quaternion differential equation of the second order for the variable λ = λ 0 +λ 1 i+λ 2 j +λ 3 k. To achieve this, the initial equation of the first order for the variable λ is differentiated with respect to time t, and the differential equations for the variables c 2 and c 3 are taken into account.
(ii) Complement the derived quaternion equation for the variable λ with the scalar differential equations of the second order for the distance r and the scalar differential equations of the first order for the energies h, h * , and modulus c of the vector of velocity moment. The right parts of the equations for h, h * , and c in this case are written in terms of the Euler parameters λ j and the distance r.
(iii) Replace the variable λ in the derived quaternion equation of the second order with the new four-dimensional variable u = u 0 +u 1 i+u 2 j +u 3 k using the formula λ = κ(r)u, where κ(r) is the regularizing function (which is a twice differentiable function with respect to the variable r). When κ(r) = r −1/2 , the new scalar variables u j (the components of the quaternion variable u) are the KS variables. As a result, we obtain the main differential quaternion equation of the second order for the variable u with the regularizing function κ(r).
On this same stage, the variables λ j are replaced with the new variables u j in the equations for the variables r, h, h * , and c using the following formulas: (iv) Apply the regularizing transformation of real time t. In order to do that, we transition from time t in the derived differential quaternion equation for the variable u to the new independent variable τ using the formula dt = ν(r)dτ , where ν(r) is the second regularizing function of the distance r. As a result, we obtain the main differential quaternion equation of the second order for the variable u with two regularizing functions κ(r) and ν(r).
A similar transition from time t to the new variable τ 1 using the formula dt = ν 1 (r)dτ 1 , where ν 1 (r) is the third regularizing function of the distance r, is performed for the equation of the second order for the distance r. As a result, we derive a differential equation of the second order for the distance r, which includes the regularizing function ν 1 (r).

Quaternion equations of the perturbed central motion with regularizing func-
tions As a result of the described transformations, we obtain the following set of quaternion equations and relations for the problem of perturbed central motion of a point particle with three regularizing functions κ, ν, and ν 1 .
The relation between the generalized KS variables (u-variables) and the Euler (Rodrigues-Hamilton) parameters is

The generalized KS transformation (u-transformation) is as follows.
For the distance, For the Cartesian coordinates, For the projections of velocity, For the perturbing forces,

The main differential quaternion equation for the generalized KS variables (u-variables) is
The differential equation for the distance is The differential equations for the variables h, h * , c, and c 2 arė The equations for time are The above equations form the determined system of differential equations for the perturbed motion of a point particle, where the unknowns are time t, the parameters u j (the generalized KS variables), the distance r, the energy h of a point particle if Π * = 0 or the energy h * of a point particle if Π * = 0, and the modulus c of vector of the moment of velocity of a point particle, or the squared modulus c 2 . Either τ or τ 1 can be chosen as an independent variable.
The equation for distance (5) can be excluded from the mentioned system of differential equations in the cases when the equation (κ(r)) 2 (u 2 0 + u 2 1 + u 2 2 + u 2 3 ) = 1 can be resolved (with a specified form of function κ = κ(r)) for the distance r, i.e., when the distance r can be expressed from the variables u j .
Note that after the transition from time t to the new independent variable τ or τ 1 in the equations for h, c, and c 2 , the form of these equations does not change. Also note that in some cases, it is efficient to use the equation for the total energy h * of a point particle instead of the equation for the variable h.
To find the coordinates x, y, and z and the projections of the velocityẋ,ẏ, andż of a point particle on the axes of the coordinate frame OXY Z from the variables u j , the distance r, and their derivatives, it is necessary to use the above-stated relations.

Conditions for transforming quaternion equations of the perturbed central
motion to the oscillator form The described equations of perturbed motion of a point particle include the regularizing functions κ, ν, and ν 1 as arbitrary functions of the distance r. They were chosen [22,[26][27] in such a way that the quaternion equation (4) for the variable u and the scalar equation (5) for the distance r, or at least one of them, can be equivalent to the motion equations for harmonic oscillators for unperturbed central motion, when Π * = 0 and p = 0. Therefore, when Q = 0, P 1 = 0, h = const. and c = const.
To make the main quaternion equation (4) for u equivalent in the said case to the motion equations of four-dimensional harmonic oscillator, it is necessary that it does not include the first derivative du/dτ . For this reason, it is required that the functions κ and ν satisfy the following condition: When the condition is met, the said equation does not contain the derivative du/dτ . The general solution to Eq. (6) has the following form: Without lose of generality, we assume b = 1 and further consider such regularizing functions κ and ν interrelated by This means that after choosing the specific function κ, the function ν will also be unambiguously defined according to the last equation, and vice versa. Therefore, instead of two arbitrary functions κ and ν, only one of them is arbitrary, κ or ν. In this case, the choice of the function ν, as well as the function κ, is limited to the class of continuous differentiable functions.
The quaternion equation (4) for u, considering the relation rκ = ν 1/2 , takes the form of motion equation of the four-dimensional nonlinear perturbed oscillator, To make the equation equivalent to the motion equation of the four-dimensional harmonic oscillator for the case of unperturbed central motion of a point particle, it is necessary to require the satisfaction of the condition as follows: where a is some constant value (for example, we can set a = h).
The expanded form of this condition is where we should set h = const. and c = const. The last relation can be viewed as a differential equation of the second order for finding the regularizing function κ(r) when the form of potential Π(r) is specified. After finding the function κ from that equation, the function ν can be unambiguously determined from the condition rκ = ν 1/2 . The condition (8) can also be viewed as the differential equation of the first order for finding the potential Π(r) when the form of regularizing function κ(r) is specified.
Therefore, the relation rκ = ν 1/2 and the relation (rκ) 3 α = a = const., the expansion of which has the form (8), are the necessary and sufficient conditions for transforming the main quaternion equation (4) of the perturbed central motion of a point particle to the oscillator form. They interrelate the regularizing functions κ, ν and the potential Π(r). If the form of Π(r) is specified and the regularizing functions κ and ν are chosen in such a way that these conditions are met (recall that we should set h = const. and c = const. in the relation for (rκ) 3 α = a), the main quaternion equation (4) for u becomes equivalent to the motion equation of singlefrequency four-dimensional harmonic oscillator with the frequency k = a/2 for the case of unperturbed central motion of a point particle with the potential Π(r).
In the case of unperturbed central motion, Eq. (5) for the distance r can also be transformed to the motion equation of harmonic oscillator by choosing an appropriate regularizing function ν 1 (r). Indeed, in this case, it takes the form of the equation investigated by Belen'kiy [70] which implies the known condition [70] , for the regularizing function ν 1 (r) and the above-mentioned potential Π 1 (r). When this condition is met, the equation for r is equivalent to the motion equation of the harmonic oscillator. Unfortunately, this condition and rκ = ν 1/2 and (rκ) 3 α = const. when ν(r) = ν 1 (r) are incompatible for the general case, i.e., for an arbitrary form of the potential Π(r).
Thus, it has been determined by Chelnokov [22,[26][27] that the vector differential equation (3) of the perturbed central motion of a point particle can be transformed, when the regularizing functions are chosen in the above-described way, to the quaternion equation (7), which has the form of an oscillator. In a general case, this quaternion equation must be complemented with the scalar equations for the distance r and the variables h or h * and c. With the appropriate choice of regularizing functions, one or two of these equations (for example, the equations for r and c in the case of KS) will fall out of consideration.

Quaternion equations of the perturbed central motion in a normal form
The above-described equations of perturbed central motion of a point particle include the second-order differential quaternion equation for u as the main equation, which, with the appropriate choice of regularizing functions, takes the form of the perturbed four-dimensional oscillator motion equation. In some problems of the perturbed central motion of a point particle, it is efficient to use the perturbed central motion quaternion equations in a normal form, where the variables are the quaternion λ which describes the orientation of the above-introduced coordinate frame η rotating in the inertial space, and the projections c i of the vector c of moment of velocity of a point particle on the axes of that coordinate frame or its projections c x , c y , c z on the axes of coordinate frame OXY Z which moves translationally in the inertial coordinate frame. These motion equations have the following form [22,[26][27] .
The equations in mappings on the rotating coordinate frame η are The equations in mappings on the inertial coordinate frame (match the equations in mappings on the main coordinate frame OXY Z) are The relation for the quaternion Q is Each of these equation sets must be complemented, in the general case, with the equation for the distance r, The mentioned equations form the determined first-order system of differential equations (9) and (11) for the perturbed central motion of a point particle relative to the unknown projections c 2 and c 3 of the velocity moment vector on the axes of coordinate frame η, the Euler parameters λ j , and the distance r, and also the determined system of Eqs. (10) and (11) relative to the projections c x , c y , and c z of the velocity moment vector on the axes of coordinate frame OXY Z, the parameters λ j , and the distance r.
The coordinates x, y, z and the velocity projectionsẋ,ẏ,ż of a point particle on the axes of coordinate frame OXY Z can be found from the mentioned variables using the following relations: Note that the order of the system of Eqs. (9) and (11), in which the mappings of vectors c andċ on the coordinate frame η are used, is less by one than that of the system of Eqs. (10) and (11), in which the mappings of these vectors on the coordinate frame OXY Z (therefore, also on the inertial coordinate frame) are used. Besides, the quaternion equation for λ in the first of these systems is simpler than that for the variable in the second of these systems, since the projection c 1 of vector c on the axis η 1 (on the direction of radius vector r) is zero. Also note that the first-order equations for c 2 , c 3 , c x , c y , c z , and λ j , as expected, do not contain Π(r) which is only present in the second-order equation for r.
In some important cases, the equation systems for the perturbed central motion of a point particle (9), (11) or (10), (11), which include the first-order quaternion differential equations, can be transformed (by replacing t and r with new variables) to the systems in which all differential equations (or most of them) are linear with constant coefficients.
We have used the above-described perturbed central motion quaternion differential equations to solve a number of problems, i.e., for regularizing the perturbed motion equations (eliminating the pole-type singularity, which exists when a central body is present); for deriving the solution of the unperturbed central motion spatial problem for any arbitrary potential Π(r) in the uniformized representation, which eliminates the necessity to consider branching of solutions around the critical points; for analytical and numerical investigation of the perturbed motion in the problems of celestial mechanics and astrodynamics; and also in inertial navigation. These equations and relations also allowed us to derive the differential equations of perturbed motion of a point particle, which use quaternion elements instead of angular osculating elements.

Systems of quaternion regular equations of the perturbed central motion
Chelnokov [22,27,29] discussed the regularization problem, i.e., elimination of singularity in the differential equations of the perturbed central motion of a point particle when it is close to the center O of the central force field, i.e., when the distance r is close to zero. The above-stated quaternion equations of the perturbed central motion with regularizing functions were used, in which these functions were chosen in a certain way.
In particular, for the regularizing function κ(r), which is defined by the relation κ(r) = r −1/2 , the condition rκ = ν 1/2 is derived that the regularizing function ν(r) = r, and the relation for the coefficient (rκ) 3 α in the main quaternion equation (4) and Eq. (7) for the variable u has the following form: It is evident that if the regularizing functions are chosen as the conditions for transforming the main quaternion equation and Eq. (5) for the distance r to the oscillator forms are satisfied only for the case of the unperturbed Keplerian motion when By substituting the mentioned expression for regularizing functions κ(r), ν(r), ν 1 (r) and the coefficient (rκ) 3 α in the perturbed central motion equations with regularizing functions, we have derived the following two equation systems for the perturbed central motion in the KS variables for any arbitrary potential Π(r) of the central force field.
The system of equations for any arbitrary potential Π(r), which contains the central motion where The system of equations for any arbitrary potential Π(r), which contains the total energy of motion The differential quaternion equations of the perturbed Keplerian motion can be derived [22,25,27] from Eqs. (12) and (14) with Π = −mµr −1 .
The system of regular equations for the perturbed Keplerian motion when Π(r) = −mµr −1 , which contains the Keplerian energy h, is The system of regular equations for the perturbed Keplerian motion, which contains the total energy h * , is The systems of Eqs. (14)- (16) are complemented with Eq. (13) for the quaternion Q, which describes acting perturbances.
The systems of Eqs. (15) and (16) are the quaternion forms of the perturbed spatial two-body problem regular equations in the KS variables, taking into account the effect of the perturbing force with the potential Π * = Π * (t, r in ) = Π * (t, x, y, z), and contain the Keplerian energy h or the total energy h * as one of the variables. The differential equation for the distance r can be excluded from these systems because of the relation r = u • u = u 2 0 + u 2 1 + u 2 2 + u 2 3 . The systems of regular equations use the Euler parameters. Other choice of regularizing functions κ(r), ν(r), ν 1 (r), than for the case of KS, and also introducing the new variable c 2 (squared modulus c of the velocity moment vector) along with the variable h * allow to derive more general regular equations of the perturbed central motion [22,[26][27] .
Let us set κ = 1 and ν 1 (r) = ν(r). Then, according to the condition rκ = ν 1/2 for transforming the main quaternion equation (4) for u to the oscillator form, we get ν = ν 1 = r 2 , and the condition (8) for the coefficient (rκ) 3 α in this equation takes the following form: Since c = const. for the unperturbed central motion, from the last relation, it follows that if the regularizing functions are chosen as κ = 1 and ν = r 2 , the condition (rκ) 3 α = const. for transforming the main quaternion equation to the oscillator form is satisfied for any arbitrary potential Π(r). However, the condition for transforming Eq. (5) for the distance r to the oscillator form is not satisfied when Taking into account the last two equations and ν 1 = r 2 , from the above-given equations with regularizing functions, we derive the following quaternion system of differential equations of the perturbed central motion of a point particle in the oscillator form for any arbitrary potential Π(r) [22,24,27] : where in which the unknowns are the Euler parameters λ j , the distance r, time t, and the variables h * and c 2 . The variable τ , which is defined by the differential relation dτ = r −2 dt, is an independent variable. For the unperturbed central motion, there exists the area integral where ϕ is a polar coordinate. Compare this relation with dt = r 2 dτ which interrelates the new independent variable τ . τ corresponds to the regularizing function ν = r 2 with time t. It follows dϕ = cdτ.
By introducing a new independent variable ϕ, defined by the differential relation dϕ = cr −2 dt in the perturbed central motion equations (17)- (20) with the independent variable τ defined by the differential relation dτ = r −2 dt, and considering that in the case of the perturbed motion, c = const., we derive the perturbed central motion equations of the following form [22,24,27] : where the quaternions Q and q describe the acting perturbances, defined by the relations in Eq. (21), in which Note that the case of linear motion of a point particle (when c = 0) should be excluded from consideration when using Eqs. (22)-(25) with the independent variable ϕ, because these equations are not defined for that case.
The main advantage of the differential equation systems of the perturbed central motion derived using the Euler parameters λ j is that each one of the second-order quaternion differential equations (17) and (22) in the Euler parameters is regular for the perturbed motion of a point particle in a central force field with any arbitrary potential Π(r) if the equation terms generated by the perturbing potential Π * and the perturbing acceleration p are finite. Besides, in the case of unperturbed central motion, each one of these quaternion equations becomes equivalent to the motion equation of four-dimensional single-frequency harmonic oscillator. The frequency of oscillator, which for this case corresponds to the quaternion equation (17), equals c/2, and its value depends on the type of motion. However, the frequency of oscillator, which corresponds to the quaternion equation (22), has the constant value of 1/2 for all types of motion, which is convenient for solving some problems. The quaternion λ in the case of unperturbed central motion characterizes the orientation of orbital plane of a point particle. It follows that each one of Eqs. (17) and (22) for the quaternion λ is a regular quaternion equation of instantaneous orientation of the perturbed orbit plane.
Each one of Eqs. (19) and (24) for the total energy h * and each one of the first equations from Eqs. (20) and (25) for c 2 are also regular for any arbitrary potential Π(r), if that same condition of finiteness of perturbing forces is satisfied. However, Eqs. (18) and (23) for r are regular only if Π(r) has the form of for which Eq. (18) becomes For this reason, the described systems of Eqs. (17)- (20) and (22)- (25), derived using the Euler parameters, are generally regular for the perturbed central motion of a point particle in a force field with the potential (26), which has the fourth order relative to the inversed distance to the center of attraction r −1 . These equations are more complex than the quaternion equations in the KS variables, because they include an additional equation for r (as we have mentioned, the equations in the KS variables can be processed independently from the equation for the distance r) and the equation for the variable c 2 . Besides, the equation for the distance is not linear for unperturbed Keplerian motion. However, the equations in the KS variables are regular only for the perturbed motion of a point particle in a central force field with the Newtonian potential Π(r) = −a 1 r −1 (if that same condition of finiteness of the perturbing forces is satisfied). Besides, the main quaternion equation (12) or (14) for the KS variables is nonlinear for unperturbed central motion with any potential Π(r), except the potential Π(r) = −a 1 r −1 , unlike the quaternion equations (17) and (22).
Note that if the perturbing acceleration p and the perturbing potential Π * are not explicit functions of time t, Eqs. (17)-(20) and (22)-(25) can be processed independently from the time equation. Besides, if the perturbing acceleration from the non-conservative forces p = 0 and the potential Π * of perturbing forces is not an explicit function of time t, the total energy h * = const. and Eqs. (19) and (24) for that variable fall out of consideration.
Also note that nowadays from comparing the obtained results of computations for orbit elements offsets for planets to the ephemerides data (Simon, Bretagnon, Chapront, Chapront-Touze, Francon, Laskar 1994) for J2000 epoch, it is determined how the difference in these data conforms to relativistic effects of the GTR. Additional offset of the perihelion of planets orbit for one revolution around the Sun, according to the GTR, is determined by the known Einsteins formula (Explanation of the perihelion motion of Mercury's from the GTR). This formula follows from the solution to the problem of motion of a point particle in the curved spacetime described by the Schwarzschild metric. Point particle trajectories in such a space match the trajectories of motion in the central force field with the potential Π(r) = −a 1 r −1 − a 3 r −3 , a 1 , a 3 = const., which is a particular case of the potential (26). In such a field, the perturbing acceleration, which consists only of the radial component proportional to r −4 , is added to the Newtonian acceleration. From the perturbed motion equations in osculating elements, it follows that such perturbation only affects the behavior of the longitude of perihelion and the eccentricity of an orbit. It is determined that the relativistic effect for the offset of Mercury's perihelion is confirmed. For perihelions of other planets, the effect predicted by the GTR is within the margin of calculation error. The described regular equations (17)-(20) and (22)-(25) can be used, in particular, to predict motion of the planets, taking into account the effects of the GTR. Also, we have used them [57] to derive regular quaternion equations of the perturbed orbital motion of a solid body in the Earth's gravitational field, considering its zonal, tesseral, and sectorial harmonics (with regularization of the equation terms, which contain negative power of the distance r up to and including the fourth order).
Regular equations of perturbed central motion were also derived [27,30] in two normal quaternion forms from Eq. (9) in mappings on the rotating coordinate frame η by introducing the new independent variable τ according to the formula dt = r 2 dτ , and the new independent variable ϕ in these equations.
The equations with the independent variable τ are as follows: These differential equations of the first order for quaternion variables c η and λ are complemented with differential equations (18) and (19) for the distance r and the total energy h * , and also with the relations (21) for the quaternion Q, which describes the acting perturbances. The equations with the independent variable ϕ are as follows: These equations are complemented with the differential equations (23) and (24) for the distance r and the total energy h * , and also with the relations (21) for the quaternion Q. From Eq. (10), one can derive the equivalent equations in mappings on the inertial coordinate frame, which are more complex and have the order which is greater by one.
Systems of Eqs. (27), (28), (18), (19) and (29), (30), (23), (24), as well as systems of Eqs. (17)-(20) and (22)-(25), are regular for the perturbed motion in the central force fields with the potential (26). However, the order of these systems is 10 (the equation systems in the KS variables have the same order), which is less by 3 than the order of systems of Eqs. (17)-(20) and (22)-(25), which contain oscillator type quaternion equations. (17)-(20) and (22)-(25), as we have noted, have a common disadvantage, i.e., the included equation for r is not linear for the general case of unperturbed central motion nor for the important particular case of the unperturbed Keplerian motion. This makes it difficult to directly use these equations in the analytical research. This disadvantage for the Keplerian motion can be eliminated by replacing the distance r with the new variable ρ = 1/r, which implies that we do not consider motion of a point particle when it collides with the central mass.

Perturbed motion quaternion equation systems with the generalized Binet equation Equation systems
By introducing a new variable ρ = 1/r in Eqs. (22)-(25), we derive the equation system [22,24,27] , where This equation system with the generalized Binet equation (32) is simpler than the system (22)-(25) because it does not include the equation for the total energy h * . Besides, Eq. (32) for the case of the unperturbed Keplerian motion takes the form of harmonic oscillator motion equation, while Eq. (23) for r for this case stays highly nonlinear. For these reasons, in some cases of analytical and numerical research of the perturbed motion, the system (31)-(33) is preferable to the system (22)- (25).
Also note the following advantages of the system (31)- (33). In the case of unperturbed Keplerian motion, the universal solution to this equation system (i.e., solution, the form of which does not depend on the type of motion) can be written in elementary trigonometric functions, while the universal solution to the KS equation system (15) for this case requires special Stumpff functions [37] . This circumstance is important for analytical and numerical research of the perturbed motion, because the type of orbit can change under the effect of perturbing forces acting over finite time spans. Besides, the system (31)-(33) allows to derive the perturbed motion equations in quaternion osculating elements, which do not contain trigonometric function of slow variables, are free from singularities associated with them, and are suitable for researching any arbitrary perturbed central motion, while the perturbed motion equations in quaternion osculating elements derived from Eq. (15) can only be used for the perturbed elliptic Keplerian motion.
The described properties of Eqs. (31)-(33) are due to the fact that the quaternion equation (31) for Euler parameters and Eq. (32) for the variable that is the inverse of the distance are equivalent to the harmonic oscillator equations for all types of unperturbed Keplerian motion (elliptic, parabolic, and hyperbolic). Besides, for the case of arbitrary unperturbed central motion, Eq. (31) is equivalent to the four-dimensional single-frequency harmonic oscillator equation. However, the KS equations (the first and second equations in system (15)) are equivalent to harmonic oscillator equations only for elliptic motion when h < 0, for parabolic motion when h = 0, while for hyperbolic motion when h > 0. One disadvantage of system (31)-(33) is that it is not regular and is not suitable for research of motion in the neighborhood of coordinate origin, unlike the systems described before. It should be noted that such nonregularity is caused by Eq. (32) for the variable ρ. However, Eqs. (31) and (33) are regular for any arbitrary perturbed central motion.
Note that introducing a new variable ρ = 1/r in the system of Eqs. (29) and (30) leads to the system of quaternion equations in the normal form of which has to be complemented with the generalized Binet equation (32) for the variable ρ and the above-stated relations for the quaternion Q, which describes acting perturbances. This system also has its advantages and disadvantages in comparison to system (29), (30), (18), and (19), as well as the system (31)- (33) in comparison to the system (22)-(25).

Problem of motion of an artificial satellite in the Earth's gravitational field
In the works [22,25] , the efficiency of using the modified four-dimensional regular variables instead of the KS variables was exemplified by the motion equations of a satellite in the Earth's gravitational field.
The vector form of the differential equations of motion of an artificial satellite in the Earth's gravitational field, considering its central and zonal components, is as follows: r is the geocentric radius vector of the satellite's center of mass, f is the gravitational constant, m E is the mass of the Earth, Π E = Π + Π * z is the potential of the Earth's gravitational field, Π = Π(r) is its central component, and Π * z = Π * z (r) is its zonal component, which results from noncentrality of the Earth's gravitational field (this component contains zonal harmonics of the Earth's gravitational field).
Note that the right-hand side of the above differential motion equation of an artificial satellite, in contrast to Eq. (3) and other equations of Section 4 (except for Subsection 4.6), does not contain the factor 1/m, since the below-given potentials Π and Π * z do not include masses. The coordinate frame OXY Z, in which the orbital motion of a satellite is considered, is introduced as follows: its origin O is at the center of the Earth, the axis OZ is directed to the northern pole of the Earth, and the axis OX is directed to the point of spring equinox. The Cartesian coordinates of a satellite in this coordinate frame are x, y, and z. Therefore, the potential Π * z is the function of these coordinates Π * z = Π * z (x, y, z). The potential which describes zonal harmonics of the Earth's gravitational field is where R is the mean equatorial radius of the Earth, J n are non-dimensional constants, which characterize the shape of the Earth, P n is the Legendre polynomial of the nth order, ϑ is the angle between the axis OZ and the vector r, and ϕ is the geocentric latitude.
In the case of the KS axis, η 1 of the above-introduced rotating coordinate frame η is directed along the radius vector r. The orientation of coordinate frame η in the coordinate frame OXY Z is characterized by the Euler parameters λ j . The coordinates x, y, and z of a satellite in the coordinate frame OXY Z are interrelated with r and the Euler parameters, and also with the KS variables u j , by the following relations: where As a result, the function γ, which is present in the potential Π * z , can be expressed in terms of the distance r, the Euler parameters, and the KS variables, in the following way: Let us define Let us use the regular equation system for the perturbed Keplerian motion (14), which contains the total energy h * . In this case, in the first equation of that system for the regular quaternion variable u, we should define the perturbing potential as Π * = Π * z and the perturbing acceleration from non-conservative forces as p = 0. From the third equation of that system, it follows that in the case of the total energy of a satellite, h * = const. As a result, from the system (14) we derive, taking into account the above-described notations and assumptions, the quaternion and scalar forms of differential equations for satellite motion are as follows: where These equations take into account central and zonal harmonics of the Earth's gravitational field, and h * is the total energy per unit mass of a satellite (earlier h * meant the total energy of the whole satellite). The right parts of Eqs. (37) and (38) do not depend on time t. Time t can be calculated by additionally integrating dt/dτ = r.
The differential equation for r of a satellite to the Earth's center of mass, taking into account the equation scal (u Let us direct the axis η 3 of the coordinate frame η, instead of the axis η 1 , along the radius vector r. In this case, all the above-stated quaternion equations preserve the same form, but we should replace the unit vector i with the unit vector k (this, by the way, demonstrates the convenience of quaternion models in astrodynamics). New variables u j defined from the Euler parameters, as in the case of KS, by Eq. (35), are interrelated with the coordinates x, y, and z by which have the following quaternion form: where the new quaternion variable u = u 0 + u 1 i + u 2 j + u 3 k has the meaning different from the quaternion variable u used earlier.
The distance r, as before, can be calculated from new variables u j according to r = u 2 0 + u 2 1 + u 2 2 + u 2 3 . The variable γ in the potential of the Earth's gravitational field can be expressed from the KS variables by the relations (36), From this, taking into account the relation (35) and equation Comparing Eqs. (36) and (40), it is clear that the variable γ, on which the potential Π * z depends, can be expressed from the new variables u j in two different forms, and these expressions have a simpler and more symmetrical structure, which allows to derive simpler and more symmetrical satellite motion equations, than those for the case of using the KS variables.
Quaternion versions of the motion equations of the Earth's satellite in the KS variables and modified variables have the same form (37). Scalar version of the satellite's motion equations in the Earth's gravitational field in new (modified) four-dimensional variables u j has the following form [22,25] : where These equations must be complemented with the time equation dt/dτ = r. The derived motion equations (41) and (42) of new (modified) variables u j have a simpler and more symmetrical structure compared with Eq. (38) of the KS variables, which simplifies analytical and numerical research of these equations. The equations produce the determined system of nonlinear differential equations of the eighth order in relation to the new variables u j (j = 0, 1, 2, 3) and their first derivatives du j /dτ . Note that Eq. (41) for u 0 and u 3 of that system, complemented with the equation for r, produces the determined system of differential equations of the sixth order. Similarly, Eq. (42) for u 1 and u 2 of that system, complemented with Eq. (43) for the distance r, also produces the determined system of differential equations of the sixth order. This simplifies the analytical research of satellite motion, because if we use the KS variables, we have to deal with the determined system of differential equations of the eighth order in the considered case. Also note that the modified variables (denoted by u * j ) are interrelated with the KS variables u j by the following relations:

Equations of unperturbed central motion
Equations of spatial unperturbed central motion of a point particle are derived from the perturbed motion equations (7), (8), and (5), taking into account the relations p = 0 and Π * = 0, and have the following form [23,29] : Equations (44)-(46) must be complemented with the following relations: which allow to calculate the coordinates x, y, and z and the projections of velocityẋ,ẏ, andż of a point particle from u j , r, and their derivatives. Equation (45) matches the equation derived by Belen'kiy [70] while considering the problems of central motion in the polar coordinates r and ϕ.
Assuming that the regularizing function κ = 1, we derive Equations (44) and (46), according to this relation, are as follows: If we replace the variable τ with the polar angle ϕ, which is related to τ by dϕ = cdτ , then Eqs. (49)-(50) become Let us put the origin of the arc coordinate ϕ to the pericenter. Then, the variable ϕ in Eqs. Relations between λ and u-variables and the elements of an orbit are as follows.
The quaternion λ, involved in motion equations for a point particle M , characterizes the orientation of rotating the coordinate frame η relative to the coordinate frame OXY Z. The coordinate frame OXY Z has its origin at the center of attraction O, and the directions of its axes stay the same over time in the inertial coordinate frame. The origin of the coordinate frame η is located at the point M , and its axis M η 1 is directed along the radius vector r of a point particle. Besides, we assume that the absolute angular velocity projection ω 1 of the coordinate frame η on the axis M η 1 is zero during the whole timespan of motion.
In the case of unperturbed central motion, it is efficient to orient the axis M η 2 of the coordinate frame η in the plane of an orbit. Then, the angular position of the coordinate frame η relative to the coordinate frame OXY Z can be defined by three angles, i.e., the longitude of ascending node Ω, the inclination of an orbit I, and the argument of latitude σ = ω * + ϕ, where ω * is the angular distance of the pericenter from the node, and ϕ is the true anomaly. The plane of orbit in this case belongs to the coordinate plane M η 1 η 2 , and the projection ω 1 =İ cos σ +Ω sin I sin σ = 0.
Let us assign intrinsic quaternions of rotation (intrinsic quaternion is the quaternion defined by its projections in the coordinate frame which is transformed by that quaternion) to the finite rotations of coordinate frame η by angles Ω, I, and σ. These quaternions are as follows: The quaternion λ, associated with the resulting rotation of coordinate frame η relative to the coordinate frame OXY Z, can be found from the quaternions λ i (i = 1, 2, 3) of the component rotations using the quaternion formula of finite rotations addition, From the last relations, we derive which, together with the relation interrelates the Euler parameters λ j with the classical orbital elements Ω, I, and ω * . Taking the relation (55) into account, we rewrite Eq. (53) in the following form: where the quaternion λ * defines the orientation of the coordinate frame η * relative to the coordinate frame OXY Z. The axis M η * 1 of this coordinate frame passes through the pericenter, and the axis M η * 2 belongs to the plane of an orbit. Therefore, the quaternion λ * defines the orientation of an orbit in space.
Besides, the quaternion λ can be represented as follows: From this relation, we can find the derivative which can also be represented as Equations (54) and (55) and the relations interrelate the variables u j and the classical orbital elements.
In the quaternion form, we have Note that the two quaternions +λ and −λ correspond to any given angular position of the coordinate frame η relative to the coordinate frame OXY Z. Each of these two quaternions define the same rotation. Therefore, the relations (54) and (55) can be obtained from the specified angles Ω, I, and σ (or ω * and ϕ) and two sets of Euler parameters, which differ only by sign, +λ j and −λ j . Each of these sets defines the same angular position of the coordinate frame η relative to the coordinate frame OXY Z. Similarly, Eqs. (61), (54), and (55) can be obtained from the specified angles Ω, I, and σ (or ω * and ϕ), and two sets of values of variables u j , +u j and −u j . Assuming κ(r) = r −1/2 in Eq. (61) and taking the Euler parameters found from the angles Ω, I, and σ according to Eq. (54), with negative signs, we obtain the formulas which interrelate the KS variables with orbital elements. These formulas match those derived by Stiefel and Scheifele [37] using a different method.

Uniformized solution to the spatial problem of unperturbed central motion
To derive the solution to the spatial problem of unperturbed central motion, we can use either Eqs. (45), (46), and (50), or Eqs. (45), (51), and (52). Let us use the latter.
By integrating Eq. (51), we obtain where α and β are constant quaternions defined by in which λ 0 and (dλ/dϕ) 0 are the values of quaternions λ and dλ/dϕ at the pericenter (when ϕ = 0). The quaternion λ 0 in the general case defines the orientation of coordinate frame η 0 , the axis M η 0 1 of which directs along the radius vector r 0 = r(0), and the axes M η 0 2 and M η 0 3 of which rotate in the space around the axis M η 0 1 in an arbitrary way (i.e., in the general case, the axis M η 0 2 may or may not belong to the plane of an orbit). However, if we define the quaternion λ 0 in such a way that the axis M η 0 2 belongs to the plane of an orbit (the axis M η 0 3 in this case will be orthogonal to that plane), then the general solution to Eq. (51) can be represented in the form of Eqs. (58), (59), or (56), (57), (60), derived from geometric considerations. By comparing Eqs. (62) and (63) with Eqs. (58) and (59), we can obtain where the quaternion λ * , which defines the orientation of the coordinate frame Y * relative to X (i.e., the orientation of an orbit), can be found from the orbital elements Ω, I, and ω * according to the second relation from Eq. (56). From Eqs. (64) and (65), it follows that the quaternions α and β are unit quaternions (i.e., with norms equal to one), and the components α j and β j of quaternions α and β satisfy the condition of orthogonality, Let us consider Eq. (45). By defining the regularizing function ν 1 (r) for the specified form of the above-mentioned potential Π 1 (r) from the relation and from Eqs. (45) and (66), we derive the linear equation for calculating r, which has a form of the known harmonic oscillator motion equation [70] , The general solution to this equation is as follows: where A and ε are the constants of integration. The relations (62), (63), and (67), or (56), (60), and (67), complemented by following from the relations (52), (47), and (48), are the general solution to the spatial problem of unperturbed central motion [23,29] . The variable τ 1 is independent here. The variable ϕ involved in Eqs.
An equivalent solution to the spatial problem of unperturbed central motion can be derived using the quaternion equations (49) and (50). In that case, the variable ϕ is replaced by τ .
Let us represent the derived solution in the uniformized form, which eliminates the necessity to consider branching of solutions around the critical points, such as poles. By using the solutions for r, ϕ, and t, which describe the motion of a point particle in the uniformized form derived by Belen'kiy [70] in the two-dimensional case, Chelnokov [23,29] derived the uniformized solution to the spatial problem of unperturbed central motion as the relations (62), (63), and r = a(γ(z; g 2 , g 3 ) + µ * ), Here, γ(z; g 2 , g 3 ) is the elliptic Weierstrass function, z is the uniformizing variable, and g 2 and g 3 are invariants. The values of constants a, µ * , l 1 , l 2 , and invariants g 2 , g 3 are determined depending on the signs of values C * i by the following relations [70] : where e is the quazi-eccentricity [70] . ϕ, involved in the relations (62), (63), (74), and (75), must be calculated in advance according to the relation (73).
Equation (75) can be represented as

Quaternion solution to the orbit orientation problem
In celestial mechanics and astrodynamics, the classical orbital elements have been widely used for studying the Keplerian motion. These orbital elements are the longitude of ascending node Ω, the inclination of an orbit I, and the angular distance of the pericenter from node ω * . The relations between the λ and u-variables, which we use in our theory, and the orbital elements, have the form of relations (54), (55), and (61).
From these relations, one can derive, as particular cases, the formulas obtained by Stiefel and Scheifele [37] , which relate the KS variables to the orbital elements. Note that applying the introduced coordinate frame η and the quaternion formalism makes it trivial to derive the interrelations of λ and u-variables with the orbital elements. The method for deriving the relations between the KS variables and orbital elements, proposed by Stiefel and Scheifele, is much more complicated.
Let us consider the calculation of orientation of orbital plane in space and the orientation of an orbit in that plane for the unperturbed Keplerian motion. We assume that we know the initial values of quaternions λ and dλ/dϕ and the variables ρ = 1/r, dρ/dϕ, and c, where ϕ, as before, is the true anomaly.
Let us define the orientation of an orbit in space by two orthogonal unit vectors a and b. Let us direct the vector a from the center of attraction to the pericenter, and the vector b along the vector v(0), where v(0) is the velocity vector of a point particle at the pericenter. To find a and b, we use where r in (0) and v in (0) are the values of quaternions r in and v in when a point particle is at the pericenter, and r(0) and v(0) are the moduli of r and v for that same moment.
From Eqs. (47) and (48), we get Taking Eq. (64) and v(0) = c/r(0) into account, from the above-described relations, we find By solving Eqs. (62) and (63) for α and β, we get To calculate the value of true anomaly ϕ, which corresponds to the specified initial position of a point particle, it is necessary to use the formulas of Duboshin [71][72] , Equations (76)-(78) proposed by Chelnokov [25] provide the quaternion solution for the orbit orientation problem for the unperturbed Keplerian motion. These relations are the basis for determining the osculating motion if the perturbed Keplerian motion is described by the differential equations in variables λ, ρ, and c.

Equations in quaternion osculating elements
In the work of Chelnokov [25,29] , the following equations of the perturbed orbital motion of a satellite (a point particle) in osculating elements α, β, and c 2 were derived from Eqs. (31) and (33) by a variation method of arbitrary constants of integration as follows: and the quaternion Q is defined by the relations given after Eq. (33). Quaternion osculating elements α and β are defined from λ, dλ/dϕ, and ϕ, which describe the rotation of an auxiliary coordinate frame η in the inertial space, by Eq. (77), the scalar osculating element c 2 is the squared modulus of the velocity moment vector c of a point particle (constant value for unperturbed central motion), and the quaternion Q characterizes acting perturbances.
Note that the quaternion osculating elements α and β introduced here are different from the quaternion osculating elements α and β, which were introduced in Section 3 and are interrelated with the regular quaternion variables u and du/dτ * .
The equations in elements α, β, and c 2 are complemented with the generalized Binet equation (32), which for the Newtonian central field of gravity (when the potential Π(r) = −mµr −1 = −mµρ) takes the form of the perturbed oscillator, or with the equations in osculating elements, derived from the generalized Binet equation and related to the variable ρ = 1/r, and also the time element equation.
We have also derived the equations for the perturbed orbital motion of a satellite in other quaternion osculating variables [29] , one of which is the quaternion variable λ * , which is determined by the second relation from Eq. (56) and characterizes the osculating orientation of an orbit of a satellite in its perturbed motion.

Unambiguous calculation of λ and u-variables from Cartesian coordinates
and their derivatives Consider the problem with the initial conditions in λ and u-variables. At t = 0, the values of coordinates x, y, and z and projections of velocityẋ,ẏ, andż of a point particle are set. We need to determine the corresponding values of λ j , dλ j /dτ , or u j , du j /dτ for the moment τ = 0, which are involved in the initial conditions for integrating the quaternion equations of the perturbed motion in λ or u-variables. For the sake of simplicity, let us denote the initial values of variables by the same letters as their current values.
The Euler parameters λ j can be found from the Cartesian coordinates x, y, and z according to the following equation [23,25] : where r = (x 2 + y 2 + x 2 ) 1/2 . By using Eq. (48), we can obtain the derivativesλ j as follows: To find the derivatives dλ j /dτ and dλ j /dϕ, it is necessary to use After we find λ j andλ j using Eq. (61), it is easy to find u j and their derivatives with respect to different variables. From the obtained formulas, one can derive, as particular cases, the formulas for the KS variables and their derivatives, obtained by Stiefel and Scheifele [37] .
There is an ambiguity in estimating the values of λ j and u j and their derivatives from specified values of x, y, z andẋ,ẏ,ż according to the above-described formulas. This ambiguity is caused by the fact that for finding two variables (λ 0 and λ 1 for the case of x 0, or λ 2 and λ 3 for the case of x < 0), there is only one equation. As a result, the variables λ j and u j and their derivatives are defined ambiguously by x, y, z andẋ,ẏ,ż according to the above-stated relations. This ambiguity also exists for finding the KS variables and their derivatives according to the formulas described by Stiefel and Scheifele [37] , since these formulas, as we have noted, are the particular case of the above-stated formulas.
The described method for finding the variables λ j and u j is inconvenient not only because there is an ambiguity in it, but also because in the case of unperturbed central motion, if we use that method for estimating the variables λ j and u j , the corresponding coordinate plane M η 1 η 2 may not match the plane of an orbit, and the axis M η 3 may not match the normal to that plane.
Let us propose another method for finding the values of λ j and u j and their derivatives from specified values of variables x, y, z andẋ,ẏ,ż, which is free from the mentioned disadvantages of the above-described method and was proposed in works of Chelnokov [23,25] . In order to do that, we need an additional assumption about the orientation of trihedral η at the considered (initial) moment of time. Let us require that for that moment of time, the axis M η 3 is perpendicular to the plane, to which r and v = dr/dt belong, by directing it along the vector c = r × v, where c, as before, is a vector of velocity moment of a moving point particle M relative to the attracting center. The axis M η 1 , as before, is assumed to be directed along the radius vector r. Then, we have where the projections c x , c y , and c z of the vector c on the coordinate frame OXY Z are defined in terms of the variables x, y, z andẋ,ẏ,ż by From the relations for the quaternions r in and c in , we get the following vector equations: where λ v = vect λ = λ 1 i + λ 2 j + λ 3 k is the vector part of the quaternion λ, the unit vectors i, j, and k of the hypercomplex space from now on are equated to the unit vectors of some orthogonal three-dimensional coordinate frame, i.e., they conform to common rules of vector algebra.
Let us introduce ϑ of finite rotation of coordinate frame η relative to the coordinate frame OXY Z as follows: Considering that λ v = λ 0 ϑ and assuming that λ 0 = 0 (thus eliminating the half-turns of coordinate frame η relative to the coordinate frame OXY Z from consideration), from the vector equations, we derive For estimating the vector ϑ , from Eq. (82), we can use either the second and third equations, or the first and fourth equations of this equation set. As a result, we get Note that the condition (84) means, in particular, that r and v should not belong at the considered moment of time to the fixed plane OXY .
After we find the projections ϑ i of vector ϑ from the specified values of x, y, z andẋ,ẏ,ż according to the relations (83), (84), and (79)-(81), from the following equations: we can estimate the values of the Euler parameters λ j , and after that, we can estimate the values of their derivatives according tȯ The values of u j and their derivatives can be found using Eq. (61). The proposed method for finding the values of λ j and u j and their derivatives from the specified values of variables x, y, z andẋ,ẏ,ż at the considered moment of time, as opposed to the method described at the beginning of this subsection, not only allows to unambiguously calculate the initial conditions for integrating the perturbed motion quaternion equations in λ and u-variables, but also allows to find for the case of unperturbed central motion such values of quaternions λ and u, for which the coordinate plane M η 1 η 2 matches the plane of an orbit of a point particle, and the axis M η 3 is orthogonal to that plane, which is important in some cases, for example, for defining the osculating motion in the variables λ j or u j from the specified values of x, y, z andẋ,ẏ,ż.

Works on the KS and quaternion regularization of the two-body problem equations by other authors
By means of quaternions in vector notation, the work of Velte [1] in 1978 gave a new derivation of the KS transformation acting from a four-dimensional parameter space into the threedimensional physical space. Using the quaternions in vector notation allowed to assign an immediate geometrical interpretation to each step in the derivation. In particular, the KS transformation appears as the Levi-Civita transformation, formulated in a rotated coordinate system. As a simple application, the explicit formulas are given for the KS transformation in the case of an elliptic Kepler motion in space.
In the introduction part, Velte noted that Volk and Waldvogel pointed out a geometrical interpretation of the KS transformation involving the Eulerian angles of a rotated coordinate system [71] .
Velte also noted that in this paper, Volk pointed out that the KS transformation already occurred in a letter of Euler to Christian Goldbach by way of a special case of the so-called Euler identity [73] .
In the note of Vivarelli [2] , the KS transformation introduced by Kustaanheimo and Stiefel into celestial mechanics was formulated in terms of hypercomplex numbers as the product of a quaternion and its antiinvolute. Therefore, it represents a particular morphism of the real algebra of quaternions for image of a three-dimensional real linear subspace, and also a natural generalization of the Levi-Civita transformation. It is shown that the quaternion matrix of the product leads to the KS matrix, and the bilinear relation and the two identities which play a central role in the KS theory are easily derived using quaternions. A suitable quaternion gauge-transformation is given, which leads to the well-known fibration of the four-dimensional space. In addition, several geometrical interpretations are brought out.
Vivarelli, in particular, noted, "the product of two quaternions suggests the interpretation of the KS transformation as a transformation of the Euclidean space R 4 , i.e., a rotation and an expansion about the origin, the ratio of expansion being the norm of q", where q is a quaternion composed of the KS variables. Note that the work of Vivarelli [2] ("The KS transformation in hypercomplex form") did not introduce a rotating coordinate system and the Euler parameters characterizing its rotation in the inertial space and associated with the KS variables, as done by Chelnokov [20] . Also, Vivarelli did not give a kinematic interpretation of the KS bilinear relation proposed by Chelnokov [20] . We also note that the work of Vivarelli was received by the editors of the journal on December 21, 1981, and the work of the present author was received on March 15, 1979. We should note that the quaternion q used by Vivarelli has the following form: where, as before, i, j, and k are the imaginary units of Hamilton's vector.
The KS transformation was presented by Vivarelli in the following quaternion form: where q * is the quaternion antiinvolute to the quaternion q. Stiefel and Scheifele [37] studied the perturbed Keplerian motion not only using the regular equations in the oscillator form and the methods of the oscillation theory, but also using the regular equations in the canonical form, for which they developed the theory of canonical KS transformation. The canonical approach to the regularization problem, which uses the KS transformation, was developed by Lidov [74][75][76] and has been widely used. In a more recent work [77] , Poleshchikov considered the generalized KS matrix and associated transformations in the theory of regularizing the canonical equations of the two-body problem.
We should also note the work by Shagov [5] , in which the differential motion equations of a spacecraft were derived for the perturbed spatial two-body problem. These equations were written in the orbital coordinate frame. For describing the motion in the inertial space, the rotation quaternion normed by means of a multiplier was used, which equals the square root of modulus c of the velocity moment vector of a spacecraft. In these equations, τ , related to time t by the differential relation dτ = (1/2)cr −2 dt, was used as the independent variable. The equations of τ are linear for the unperturbed Keplerian motion of a satellite.
In the works of Deprit et al. [6] , the matrix KS theory was presented in terms of quaternions. The authors considered the KS transformation independently from its possible application to Keplerian systems and refined several theorems formulated by Stiefel.
The abstract reads, "Concerning Stiefel's own contribution to the question, on the one hand, we abandon the formalism of the matrix theory to proceed exclusively in the context of quaternion algebra; on the other hand, we explain how, in the hierarchy of hypercomplex systems, both the KS transformation and the classical projective decomposition emanate by doubling from the Levi-Civita transformation." The introduction notes, "Free as we are now of computational servitudes, we even put ourselves to the task of resetting the whole KS theory in terms of quaternions. By the time Stiefel and Scheifele had completed their monograph, they became aware of the close connection between their matrix formalism and the theory of quaternions. It was suggested that they take advantage of it; they reacted to the suggestion in excessive terms (Stiefel and Scheifele [37] , page 286). Did they really believe that a transfer from matrices to quaternions would lead to "failure or at least to a very unwieldy formalism"? Stiefel's dire predictions notwithstanding, we accepted the challenge. Did we fail? The reader is our jury. Building the KS transformation as the emanation of an alternate bilinear form over the algebra of quaternions costs no more in complications than the matrix formalism of Stiefel and Scheifele. Besides, we find rewards in the exercise: theorems are sharpened, some to a significant extent; proofs are shortened; the overall design of Stiefel comes out much enhanced as to its global and intrinsic meaning, not to mention as a collateral a programming style for manipulating quaternions through general purpose Symbol Processors." Deprit et al. [6] also considered the linearization of motion equations in cylindrical, spherical, and orbital coordinate frames by introducing new variables and a new independent variable instead of time t.
Deprit et al. [6] proposed their own transformation, that is, the DEF-transformation, as an alternative to the KS transformation. About this transformation and the known Burdet-Ferrandiz (BF) transformation, which they also considered, the following is stated.
"About alternatives to the KS transformation, Stiefel and Scheifele (1971, page 288) issued a warning little short of an injunction: 'the authors are convinced that the search for other transformations is not very promising.' To many readers, their omen conveyed the impression that the KS technique is unique in achieving jointly regularization, linearization, and dimensionraising for three-dimensional Keplerian systems. The facts disallow the claim. Kustaanheimo, Stiefel, and Scheifele never mentioned that the decisive step Fock (1935;1936) had taken in that direction thirty five years earlier, not even in their brief allusion to the same step taken, but independently, by Moser (1970). We have room here for the latest entry in the competition, the BF transformation due to Burdet (1969) for the coordinate part and completed by Ferrfindiz (1986a;1986b;1987;1988) for the moment part. We supplement it with a transformation of our own, the DEF-transformation, which we claim achieves equally well all the objectives of the KS transformation-linearization, regularization, and canonicity-although, we are inclined to believe, in a simpler and more intuitive way (Subsection 4.1). Admittedly, the construction involves heavy algebraic manipulations, no more however than is the case with the KS or the BF transformation. Besides, we pass the chore to the symbol processor. Yet, lest we create misunderstandings we hasten to emphasize the obvious, that a symbol processor is no more than a mathematical accountant. However efficient its accounting methods are, it cannot relieve users from the responsibility of creating simplifications toward clamping final results in their most significant form." The eight-dimensional DEF-transformation of the Cartesian coordinates and the projections of velocity of a point particle are as follows (using the notation of the authors of that transformation): where x = r is the radius vector of a point particle, X = v =ṙ is its velocity vector, and u 0 , U 0 and u, U are the new scalar and three-dimensional vector variables. The independent variable (time t) is replaced with the generalized true anomaly f such that where β is a parameter. The DEF-transformation has much in common with the BF-transformation, which is as follows: After transitioning to the generalized true anomaly f as the independent variable, the motion equations of a point particle in the Newtonian gravitational field (43) [6] in the DEF variables are as follows: where µ is the gravitational constant, Q = const. is the angular momentum, and h is the Keplerian energy (constant value). By differentiating the third equation of the described equation system relative to f , Deprit et al. [6] derived the linear differential equation of the second order relative to the vector variable u, For the new scalar variable σ = Q 2 /(µu 0 ), which is introduced instead of u 0 , considering the first and second equations of the above system, we obtain the following linear differential equation relative to the scalar variable σ [6] : After derivation of two last equations, Deprit et al. [6] stated that this completes the linearization of Keplerian systems in three dimensions by means of the DEF-transformation. There are no comments in Ref. [6] about the main fourth dynamic equation of motion of a point particle (equation of the first order for U ) of the derived equation system, even though it contains in the right part of the vector term, which, after substituting into the analytical solutions u 0 = u 0 (f ) and u = u(f ), derived as a result of integrating the last two equations, becomes a nonstationary expression which is an explicit function of time f . Moreover, the equation for U becomes a differential nonhomogeneous linear vector equation, the homogeneous part of which has constant coefficients and the nonhomogeneous part is an explicit function of time f , which does not allow to derive an analytical solution to this equation in a simple form. The quaternion equations of motion of a point particle in the Newtonian gravitational field in the KS variables are as follows: The first quaternion equation is equivalent to the system of four independent linear homogeneous differential equations of the second order relative to the components u j of quaternion variable u (the KS variables). These equations have constant coefficients equal to the halved constant Keplerian energy h of a point particle, taken with a negative sign, and are easily integrated in elementary functions for the case of an elliptic Keplerian motion, when h < 0, or, for any sign of energy h, in Stumpff functions. The initial conditions for integrating these equations are expressed from the initial values of the Cartesian coordinates of a point particle in the inertial coordinate frame and the projection of its velocity vector on the axes of that coordinate frame. The equation for time t can be integrated independently from these equations.
In our opinion, comparing the described equations shows that the equations in the KS variables are advantageous to the equations in the DEF variables. Note that the equations in the KS variables [28] lead to the equations in quaternion osculating elements (in slow quaternion variables), which are convenient for researching perturbed elliptic motion of a point particle.
Among other early works, we should note the papers by Vrbik [7][8] which showed the efficiency of using quaternions to solve the perturbed Kepler problem.
In particular, the paper by Vrbik [8] supplied a proof of formulas for constructing a perturbative solution to the perturbed Kepler problem by utilizing the quaternion algebra of the KS formulation of equations. The main advantage of this approach is a removal, from the corresponding solution, of fast oscillations (in the case of conservative forces) and small divisors (in the case of time-dependent forces).
Among more recent works, we should note the papers by Waldvogel [9][10] , who has works in collaboration with Stiefel and Waldvogel [78] and Stiefel et al. [79] . Waldvogel [10] in his work "Quaternions for regularizing celestial mechanics: the right way" stated that quaternions have been found to be the ideal tool for describing and developing the theory of spatial regularization in celestial mechanics.
Waldvogel [10] wrote, "This article corroborates the above statement. Beginning with a summary of quaternion algebra, we will describe the regularization procedure and its consequences in an elegant way. Also, an alternative derivation of the theory of Kepler motion based on regularization will be given. Furthermore, we will consider the regularization of the spatial restricted three-body problem, i.e., the spatial generalization of the Birkhoff transformation. Finally, the perturbed Kepler motion will be described in terms of regularized variables." Waldvogel acknowledged the authority of the author of this paper in the field of quaternion regularization of differential equations of the perturbed spatial two-body problem and, talking about the above-cited statement by Stiefel and Scheifele on the futility of using quaternion matrices in the regularization theory, wrote, "This statement was first refuted by Chelnokov (1981) who presented a regularization theory of the spatial Kepler problem using geometrical considerations in a rotating coordinate system and quaternion matrices. In a series of papers, including Chelnokov (1992Chelnokov ( , 1999, the same author extended the theory of quaternion regularization and also presented practical applications." Let us note the main characteristics of the quaternion regularization method of Waldvogel [10] . For regularization, he proposed to use "star conjugate of the quaternion", where u = u 0 + iu 1 + ju 2 + ku 3 , and also the mapping This mapping uses nonconventional representation of the three-dimensional vector x by the quaternion x = x 0 + ix 1 + jx 2 with zero k-component (note that Waldvogel did not use the special symbol • for quaternion multiplication). x is a formal generalization (extension) of the complex variable x = x 0 + ix 1 , used by Levi-Civita in the theory of regularizing the two-dimensional motion equations.
The above-stated Waldvogel mapping, considering his previous formula, becomes In the scalar form from the last formula, we get "which is exactly the KS transformation in its classical form or up to a permutation of the indices, the Hopf mapping" (statement by Waldvogel [10] ). Note that Waldvogel quaternions x, u, and u * match Vivarelli [2] quaternions x, q, and q * up to a designation of the indices of their components.
In the classical quaternion theory, a three-dimensional vector x is associated with the quaternion x = ix 1 + jx 2 + kx 3 with zero scalar part. In the works by the author of this paper, the quaternion variable u = u 0 + iu 1 + ju 2 + ku 3 , which differs (by meaning) from the Waldvogels quaternion variable, and the quaternion x = ix 1 + jx 2 + kx 3 with zero scalar part are used for regularization. These works used the mapping x = uiu and also the mapping The scalar version of the first of these mappings is exactly the KS transformation, the form of which differs from the above-described Waldvogel transformation. We should note that Waldvogel [9][10] derived the elegant quaternion representation of the Birkhoff spatial mapping [80] which was used in the theory of regularization of the restricted three-body problem equations. The same mapping, known as the Joukowsky mapping, was used in aerodynamics for mapping the cross sections of aerodynamic surfaces to approximately circular shape. In the work of Waldvogel [10] , this mapping was called the Joukowsky-Birkhoff mapping. Waldvogel presented this mapping as an addition to his early works on the regularization theory proposed by Stiefel and Waldvogel [78] and Stiefel et al. [79] .
Saha [11] showed that the KS transformation in the quaternion form can be interpreted as a rotation in three dimensions using the rotation axis and the angle of rotation, who wrote, "The KS transform turns a gravitational two-body problem into a harmonic oscillator, by going to four dimensions. In addition to the mathematical-physics interest, the KS transform has proved very useful in N -body simulations, where it helps handle close encounters. Yet the formalism remains somewhat arcane, with the role of the extra dimension being especially mysterious. This paper shows how the basic transformation can be interpreted as a rotation in three dimensions. For example, if we slew a telescope from zenith to a chosen star in one rotation, we can think of the rotation axis and angle as the KS transform of the star. The non-uniqueness of the rotation axis encodes the extra dimension. This geometrical interpretation becomes evident on writing KS transforms in quaternion form, which also helps derive concise expressions for regularized equations of motion." Zhao [12] presented the KS regularization of a spatial Kepler problem in symplectic and quaternion approaches.
Roa and Pelaez [15] showed two fully regular and universal solutions to the problem of spacecraft relative motion derived from the Sperling-Burdet (SB) and KS regularizations. There are no singularities in the resulting solutions, and their form is not affected by the type of reference orbit (circular, elliptic, parabolic, or hyperbolic). In addition, the solutions to the problem are given in compact tensorial expressions and directly refer to the initial state vector of the leader spacecraft. The SB and KS formulations introduce a fictitious time by means of the Sundman transformation. Because of an alternative independent variable, the solutions are built based on the theory of asynchronous relative motion. This technique simplifies the required derivations. Closed-form expressions of the partial derivatives of orbital motion with respect to the initial state are provided explicitly.
Among the latest works on KS transformations, we should note the works by Roa et al. [14] "Stability and chaos in Kustaanheimo-Stiefel space induced by the Hopf fibration" and Ferrer and Crespo [19] "Alternative angle-based approach to the KS map. An interpretation through symmetry".
Among the latest works on Levi-Civita and KS transformations, we should also note the works by Breiter and Langner [16][17][18] . Breiter and Langner [18] noted that "The KS transformation gained popularity in the matrix-vector formulation of Kustaanheimo and Stiefel (1965), but it is much easier to interpret and generalize in the language of quaternion algebra, very closely related to the original spinor formulation of Kustaanheimo (1964)." Saha [11] described the KS transformation in a quaternion form, proposed its generalized definition, developed the geometric interpretation of the KS variables, and considered the bilinear form and its generalization.
Note, however, that the main goal of the work by Breiter and Langner [18] was "to derive an alternative set of the action-angle variables which is not based upon the notion of an orbital plane (thus avoiding singularities when the orbit degenerates into a straight segment) and to test it on some well-known astronomical problem (Lidov-Kozai problem)." In relation to that, the work stated, "But those who want to benefit from the wealth of canonical formalism require a set of action angle variables of the regularized Kepler problem. The first step in this direction can be found in the monograph by Stiefel and Scheifele (1971), where the symplectic polar coordinates are introduced for each separate degree of freedom. However, this approach does not account for degeneracy of the problem and thus is unfit for the averaging-based perturbation techniques. Moreover, no attempt was made to relate this set to the constraint known as the "bilinear invariant", effectively reducing the system to three degrees of freedom. Both problems have been resolved by Zhao (2015), who proposed the LCF variables [presumably named after Levi-Civita (1906) and Fejoz (2001)]. In his approach, the motion in the KS variables is considered in an osculating "Levi-Civita plane" (Deprit et al. 1994) as a two-degree-of-freedom problem. The third degree of freedom is added by the pair of action-angle variables orienting the plane. The redundant fourth degree is hidden in the definition of the Levi-Civita plane. The transformed Keplerian Hamiltonian depends on a single action variable, the other two actions being closely related to the angular momentum and its projection on the polar axis. Interestingly, the result is identical to the "isoenergetic variables" found by Levi-Civita (1913) without regularization." Breiter and Langner [18] also analyzed the Lidov-Kozai problem in terms of LKS variables, which allow a direct study of stability for all equilibria except the circular equatorial and the polar radial orbits.
Many other works compared the accuracy of numerical solution to the regular equations in KS variables with the solutions to other regular equations proposed in these works.
Pelaez et al. [42] performed a comparison with the equations proposed in this work, including the equations in the Euler parameters, which characterize the orientation of an orbital coordinate frame, and used the distance, its derivative with respect to the new independent variable (and then the values which are the integration constants in the solutions to the unperturbed motion equations for these variables), the Euler parameters, and the moment of orbital momentum.
Bau et al. [43] proposed seven spatial elements as state variables of the new special perturbation method for the two-body problem. These elements included the element inverse to the doubled total energy and two first integrals of the unperturbed motion and the time element. The new elements retain zero eccentricity and inclination and negative total energy. The new motion equations were proposed, written in the orbital coordinate frame, the orientation of which was described with the Euler parameters. Equations for the Euler parameters were written in scalar and matrix forms. The results were compared with the numerical solutions to the proposed equations and other equations, including the equations in the KS variables. From the results of the numeric solutions, we can conclude that for the considered problems, the equations in the KS variables and the equations proposed by the authors provide better accuracy. Moreover, the solutions to the new equations have smaller errors than those to the equations in the KS variables.
Amato et al. [44] showed that special perturbation methods based on regularized formulations can compete and even outperform the semi-analytical methods for the long-term propagations (on the order of decades) of objects orbiting around the Earth. For this kind of applications, the Cowell formulation has never been used because of the required small integration step sizes, which causes strong accumulation of round-off error and long computational time. They developed a FORTRAN code, named THALASSA, which included Cowell's method, EDromo, the KS regularization [52] , and a set of regular elements obtained by Stiefel and Scheifele [37] from the equations in the KS variables. A sophisticated numerical solver, named Livermore solver for ordinary differential equations with automatic root-finding (LSODAR), has been included to integrate the differential equations of motion.
The authors of that work noted that they presented a collection of non-averaged methods based on the integration of existing regularized formulations of the equations of motion through an adaptive solver, and that they showed for the first time that efficient implementations of non-averaged regularized formulations of the equations of motion, and especially of non-singular element methods, are attractive candidates for the long-term study of high-altitude and highly elliptical Earth satellite orbits.
Bau and Roa [45] presented a new method for computing orbits in the perturbed two-body problem, in which the position and velocity vectors of the propagated object in the Cartesian coordinates were replaced by eight orbital elements, i.e., constants of unperturbed motion. Two of them are related to radial motion, and the next four, given by the Euler parameters, define the orientation of the intermediate coordinate frame, the evolution of which follows the orientation of the orbital plane and the direction of reference on it. The total energy and the time element complete the state vector. Numerical tests are included to evaluate the performance of the proposed special perturbation method. On the example of orbital motion of two comets, it is shown that the computation with this method is much more accurate and faster than the classical computation of the orbits in Cartesian coordinates.
6 Quaternion regularization of the equations for the perturbed spatial twobody problem using Levi-Civita variables and Euler parameters In the introduction, we have already noted that Levi-Civita, regarding his attempts to generalize his regularization of equations of two-dimensional two-body problem to the spatial problem, admitted later [48] that the problem in space has long resisted his efforts, as he tried to approach it by similar coordinate changes.
Stiefel and Scheifele [37] also wrote that Levi-Civita tried hard to generalize his method of regularizing the differential equations of two-dimensional motion in the two-body problem to the general spatial two-body problem, but without success.
Aarseth and Zare [58] wrote that because of fundamental difficulties originally clarified by Hopf [59] and Hurwitz [60] , it is impossible to generalize the Levi-Civita transformation to an equivalent set of three-dimensional variables. The Levi-Civita transformation was also covered by Arseth [81] .
However, Chelnokov [31] demonstrated that the Levi-Civita regularization can be used successfully to derive the regular equations of the perturbed spatial two-body problem. We accomplished that by using ideal rectangular Hansen coordinates, regular Levi-Civita variables, the orientation quaternion of ideal coordinate frame (this name for the coordinate frame was apparently introduced by Deprit [82] , in which the differential equations of the perturbed spatial two-body problem were written), the Keplerian energy as an additional variable, and a new independent variable. Let us consider this quaternion regularization of the perturbed spatial two-body problem equations, which was proposed by Chelnokov [31] . 6.1 Two-body problem equations, written in orbital and ideal coordinate frames using the quaternion osculating element The perturbed spatial two-body problem equations, written in orbital and ideal coordinate frames using the quaternion osculating element, are as follows [24][25]31] : 2Λ = Λ•Ω ξ , Ω ξ = Ω 1 i + Ω 2 j = (r/c)p 3 (cos ϕi + sin ϕj).
Equation (88) is the quaternion version of the scalar system of differential equations (86) and (87) in the Euler parameters Λ j , which characterize the orientation of the ideal coordinate frame in the inertial coordinate frame. The quaternion variable Λ has a meaning of quaternion osculating element of an orbit of the second body. When the component p 3 of the perturbing acceleration vector of the center of mass of the second body (body in interest), which is perpendicular to the plane including the radius vector r and the velocity vector v of that body, i.e., the plane of instantaneous orbit of the second body, is zero, the quaternion Λ = const.
Equations (85)- (87) or (85) and (88) are the motion equations of the second body, written in two rotating coordinate frames, i.e., the orbital coordinate frame η and the ideal coordinate frame ξ. Equation (85) is written in the orbital coordinate frame, and Eqs. (86) and (87) (or (88)) are written in the ideal coordinate frame. In these equations, the variables are the distance r from the center of mass of the second body to the center of mass of the first body, the derivativeṙ (the projection v 1 of the velocity vector v of the center of mass of the second body on the direction of radius vector r (on the axis η 1 of the orbital coordinate frame)), the modulus of orbital velocity moment of the second body c = |r × v|, the generalized true anomaly ϕ, and the Euler (Rodrigues-Hamilton) parameters Λ j (j = 0, 3), which characterize the orientation of the ideal coordinate frame ξ in the inertial coordinate frame OXY Z. p k are the projections of p of the perturbing acceleration of the center of mass of the second body on the axes of orbital coordinate frame η.
Ω of the absolute angular velocity of the ideal coordinate frame ξ is parallel to the radius vector r of the center of mass of the second body and is defined by The projections Ω i of Ω on the axes of the ideal coordinate frame ξ are where Ξ i are the projections of r on the axis ξ i .
The Cartesian coordinates x, y, z in the inertial coordinate frame and the projections v k of the velocity vector of the center of mass of the second body on the axes of orbital coordinate frame can be found from the mentioned variables using where λ j are the Euler parameters, which characterize the orientation of the orbital coordinate frame in the inertial coordinate frame and are calculated in advance from the variables Λ j using or using the quaternion formula The projections v x , v y , and v z of the velocity vector of the center of mass of the second body on the axes of the inertial coordinate frame are defined by the quaternion relations, and the relations of the projections p x , p y , and p z of the perturbing acceleration vector p on the axes of the inertial coordinate frame with their projections p k on the axes of the orbital coordinate frame are defined by the quaternion remapping relations, Equations (90)-(92) or (90), (96) are the motion equations of the second body, written in the ideal coordinate frame ξ. In these equations, the variables are the ideal rectangular Hansen coordinates Ξ i , their first derivatives with respect to timeΞ i , and the Euler parameters Λ j (j = 0, 3), which characterize the orientation of the ideal coordinate frame ξ in the inertial coordinate frame OXY Z. p ξk are the projections of p of the perturbing acceleration of the center of mass of the second body on the axes of the ideal coordinate frame ξ.
Note that the scalar equations of the perturbed Keplerian motion in the Hansen coordinates and Euler parameters, written in the ideal coordinate frame, were derived earlier in other forms and by other methods by Stiefel and Waldvogel [78] and Brumberg [53] . Equations (90)-(92), derived independently by the author of this paper, match the equations of Deprit and Brumberg.

Quaternion regular equations of the perturbed spatial two-body problem in
Levi-Civita variables and Euler parameters Let us introduce the rotating coordinate frame Oξ with the origin in the center of attraction (center of mass of the first body) and with the coordinate axes parallel to the axes of the above-introduced ideal coordinate frame ξ. The orientation of the orbital coordinate frame η in the coordinate frame Oξ is characterized by the rotation quaternion Φ, which is as follows: The components Φ j (j = 0, 3) are defined by The above-introduced ideal rectangular Hansen coordinates Ξ i , which are defined by Eq. (89), are the Cartesian coordinates of the center of mass of the second body in the coordinate frame Oξ and are related to r and Φ j by Let us introduce the Levi-Civita variables as follows: which are related to the Hansen coordinates by The projections of velocity vector of the center of mass of the second body on the axes of the coordinate frame Oξ (on the axes of the ideal coordinate frame ξ) are related to the time derivatives of Levi-Civita variables by The Keplerian energy h and the modulus c of the orbital velocity moment vector of the second body are defined by while in new Levi-Civita variables, they are as follows: Equations (116)-(119) must be complemented with the quaternion remapping relations (97)-(99).
Equations (108)-(111) or (116)-(119) are the regular equations of the perturbed spatial twobody problem, derived by the ideal rectangular Hansen coordinates. In the scalar equations (108)-(111), the regular variables are the Levi-Civita variables U 0 and U 3 , which describe the motion of the center of mass of the second body in the ideal coordinate frame Oξ, the Keplerian energy h, time t, and the Euler parameters Λ j (j = 0, 3), which characterize the orientation of the ideal coordinate frame in the inertial coordinate frame. In Eqs. (116)-(119), the regular variables are the quaternion U , which describes the motion of the center of mass of the second body in the ideal coordinate frame Oξ, the Keplerian energy h, time t, and the quaternion Λ (quaternion osculating element), which characterizes the orientation of the ideal coordinate frame in the inertial coordinate frame. p ξk are the projections of p of the perturbing acceleration of the center of mass of the second body on the axes of the ideal coordinate frame. These values are calculated from the projections of p on the axes of inertial coordinate frame using the above-described remapping relations.
The regular equations (108)-(111) or (116)-(119) of the perturbed spatial two-body problem form the system of nonlinear nonstationary differential equations of the tenth order (the KS regular equations have the same order) and have all the advantages of the KS equations.
(i) Unlike the Newtonian equations, they are regular at the attraction center.
(ii) For the unperturbed Keplerian motion, they are linear and are as follows: For the elliptic Keplerian motion, when the Keplerian energy h < 0, these equations are equivalent to the motion equations of the two-dimensional single-frequency harmonic oscillator, the squared frequency of which equals the halved Keplerian energy, taken with a negative sign.
(iii) They allow to develop a unified approach to study all three types of Keplerian motion.
(iv) They are close to the linear equations for the perturbed Keplerian motion.
(v) They allow to present the right-hand sides of differential equations for the motion of celestial and cosmic bodies in the polynomial form, which is convenient for solving them using computers.
In addition, these regular equations have significant differences. (i) For the unperturbed elliptic Keplerian motion of the body in interest, they are equivalent not to the motion equations of four-dimensional single-frequency harmonic oscillator, as in the case of KS, but to the motion equations of two-dimensional single-frequency harmonic oscillator, because the orientation quaternion of an ideal coordinate frame, for which the new motion equations are written, is constant for this case.
(ii) For the perturbed motion of the body in interest, the orientation quaternion of the ideal coordinate frame is a quaternion osculating element (i.e., slowly changing quaternion variable), which is also a useful property of these equations, which allows to efficiently use the methods of nonlinear mechanics.
Note, however, that these equations are not suitable for studying rectilinear orbits, when the modulus of orbital velocity vector of the second body c = |r × v| becomes zero, because the quaternion differential equation of orientation of the ideal coordinate frame degenerates in this case.
Note that Brumberg [53] described the application of Euler parameters for derivation of the perturbed spatial two-body problem motion equations in Hansen coordinates and showed the possibility to further transform the perturbed motion equations, derived by him, using the Levi-Civita parabolic coordinates.
7 Developing the quaternion regularization of the two-body problem equations in the work of the author of this paper and the applications of quaternion regular models 7.1 Developing the quaternion regularization of the two-body problem equations Chelnokov [31] considered other singularities (divisions by zero) resulting from the classical equations in angular variables (in particular, in Euler angles) in celestial mechanics and astrodynamics, and eliminated using the Euler (Rodrigues-Hamilton) parameters and Hamilton quaternions. Fundamental regular in the described sense quaternion models of celestial mechanics and astrodynamics is considered, i.e., the equations of orbital motion written in non-holonomic, orbital, and ideal moving trihedrals, the rotational motion defined using the Euler parameters and rotation quaternions, and the quaternion equations of orientation of instantaneous orbit of a celestial body (spacecraft). Also, the new quaternion regular equations of the perturbed spatial two-body problem (free from singularities of classical equations of this problem, induced by the Newtonian gravitational force) are derived in Levi-Civita variables and Euler parameters, which, along with the known advantages of regular KS equations, have their own additional advantages.
In all works on the problem of regularization of differential equations of the perturbed spatial two-body problem known to the author of this paper, the regularization of equations of motion of the center of mass of the second body (body in interest) relative to the coordinate frame moving translationally in relation to the inertial coordinate frame is considered, i.e., the regularization of equations of absolute motion of the body in interest is considered. In our work [35] , the regularization of equations of relative motion of the body in interest was considered within the context of the perturbed spatial two-body problem. We proposed the regular quaternion differential equations of perturbed motion of the center of mass of the second body (body in interest) relative to the coordinate frame, which rotates in relation to the inertial coordinate frame in accordance to an arbitrarily specified law, and also in relation to the coordinate frame associated with the Earth, which is viewed as the first (central) body. The first integrals and the analytical solution were derived for regular quaternion differential equations of unperturbed motion of a body in interest relative to the Earth using the Stumpff functions. The proposed equations can be used, for example, for studying the motion of an artificial satellite (or a planetary satellite) of the Earth in relation to the coordinate frame associated with the Earth (or another planet). These equations were convenient to study the motion of an Earth artificial satellite, since many forces acting on the satellite, and primarily the Earth's gravity force (the potential of the Earth's gravitational field), depend on relative coordinates of the satellite location (on its geographical coordinates, latitude, and longitude), instead of its absolute coordinates (in the coordinate system moving translationally in relation to the inertial coordinate frame). 7.2 Applications of quaternion regular models in the theory of orbital motion of a solid body (spacecraft and artificial satellite) in the Earth's gravitational field Chelnokov [36] derived the regularized quaternion equations of the perturbed motion of an artificial satellite in the Earth's gravitational field, taking into account its zonal, tesseral, and sectorial harmonics in four-dimensional KS variables and in modified four-dimensional variables, for which the motion equations of a satellite have a simpler and more symmetrical structure compared with those in the KS variables. These equations use energy of motion of a satellite and time as additional variables. The new independent variable is related to time by a differential relation, which includes the distance from a satellite to the Earth's center of mass.
Chelnokov [57] proposed the regularized quaternion models of perturbed orbital motion of a solid body, which are free from native singularities of classical models, for a case when a body is moving in the Newtonian gravitational field and, in the general case, when a body is moving in a central force field, the potential of which is a polynomial of negative powers of the distance to the center of attraction of the fourth order.
Different from the models proposed in Ref. [36], Chelnokov [57] proposed regularized quaternion models of orbital motion of a solid body in the Earth's gravitational field, which was defined considering not only central (Newtonian), but also zonal, tesseral, and sectorial harmonics of the gravitational field potential, which take into account nonsphericity of the Earth. In these models, the negative powers of distance to the center of attraction were lower by several orders in terms that define the effect of zonal, tesseral, and sectorial harmonics of the Earth's gravitational field potential on the orbital motion of a solid body. The main variables are the Euler parameters, the distance from the center of mass of a body to the center of attraction, the total energy of orbital motion of a body, and the squared modulus of orbital velocity moment vector of a body (or the projections of this vector). The models used a new independent variable, related to time by the differential relation which involves the squared distance from the center of mass of a body to the center of attraction.
The convenience and efficiency of using the regularized quaternion models of orbital motion of a solid body, derived by Chelnokov [36,57] , for analytical research of motion of a body, are exemplified by the motion of a body in the Earth's gravitational field, which is defined taking into account the central and zonal harmonics of the gravitational field potential. The first integrals are found for the derived orbital motion equations. The substitutions of variables and transformations are proposed for these equations, which allow to derive, in order to study the motion of a body, the determined systems of differential equations of lower orders, in particular, the system of equations of the third order for distance, for the sine of geocentric latitude and for the squared modulus of orbital velocity moment vector. The equations proposed in these works are convenient for applying the methods of nonlinear mechanics and for high accuracy numeric calculations.

Optimal control of orbital motion of a spacecraft, inertial navigation
Based on the earlier works, Chelnokov [28] analyzed the optimal motion control for a point particle in the Newtonian (central) gravitational field using the maximum principle and the quaternion models of orbital motion, written in various rotating coordinate frames (including the orbital coordinate frame), which include the quaternion equations of angular motion of these coordinate frames in the Euler parameters, and also using the quaternion models in the regular KS variables. In particular, the differential equations were derived for the perturbed elliptic motion of a point particle (taking into account the perturbing and controlling accelerations) in quaternion osculating elements, which correspond to the quaternion regular variables u and u/dτ (dt = rdτ ). The advantages of using the quaternion regular models of orbital motion of a point particle for solving the problems of optimal control of that point particle are listed. The new quaternion first integrals for the differential boundary problems of optimal control of orbital motion of a point particle and the regular quaternion analytical solutions for differential phase and adjoint equations for the passive motion sections are described. The analytical solution is derived for the quaternion equations of orbital motion in KS variables for the case when an external force, which is constant by modulus and direction (in an inertial coordinate frame), is applied to a point particle. It is shown that the efficiency of analytical investigation and numerical solution of the boundary-value problems of optimal control of a point particle can be increased by applying quaternion regular models of orbital motion.
Also note that the book of Chelnokov (in Russian) [29] , in particular, presented the quaternion method for regularizing the differential equations of the perturbed spatial two-body problem and the perturbed central motion of a point particle, described the quaternion regular models for celestial mechanics and astrodynamics, and demonstrated their applications to a number of problems of optimal control of the trajectory motion of a spacecraft.
Chelnokov [30,32] also considered regularization of the models of celestial mechanics and astrodynamics, derived within the context of the perturbed spatial two-body problem and the perturbed central motion of a point particle. The fundamental regular quaternion models in celestial mechanics and astrodynamics, which are free from singularities induced by gravitation and other central forces, are presented. It is shown that the efficiency of analytical investigation and numerical solution of boundary-value problems of optimal control of trajectory motion of spacecrafts center of mass can be increased by applying quaternion models of astrodynamics.
In particular, Chelnokov [32] analyzed the basic problems, which arise when solving the problems of optimal control of spacecraft trajectory motion using the maximum principle (including the Lyapunov instability of solutions for conjugate equations of the maximum principle). The use of quaternion models of astrodynamics is shown to allow the elimination of singularities in the differential phase and conjugate equations and in their partial analytical solutions, the derivation of the new quaternion first integrals, the derivation of general solutions to differential equations for phase and conjugate variables on the sections of spacecraft passive motion in the simplest and most convenient form, which is important for solving the problems of optimal spacecraft impulse transfers, the extension of the possibilities for the analytical investigation of differential equations of boundary-value problems with the purpose of identifying the basic laws of spacecraft optimal control and motion, the improvement of computational stability of the solution of boundary value problems, and the decrease in required amount of computation.
The quaternion regular models of astrodynamics in the KS variables, derived by us, were used to efficiently solve a number of boundary problems of optimal control of orbital motion of a spacecraft [85][86][87] using the maximum principle, and also to derive the new regular quaternion equations and algorithms for inertial navigation in space [88][89] .
In particular, Chelnokov [88][89] proposed the new quaternion regular equations and algorithms for spatial inertial navigation systems with azimuthally stabilized platform and with gyrostabilized platform, which retains its orientation in an inertial reference frame, and also proposed strapdown inertial navigation systems in regular four-dimensional KS variables and in the versions of these variables. These equations take into account the zonal, tesseral, and sectorial harmonics of the Earth's gravitational field and are dynamically analogous to the quaternion regular equations of the perturbed spatial two-body problem in the KS variables, which allows to use the results, obtained in the theory of regular celestial mechanics and astrodynamics, for space inertial navigation.

Perspectives of quaternion regularization development: the restricted three-body problem
Further development of the ideas of quaternion regularization of equations of celestial mechanics and astrodynamics is associated, in our opinion, with the construction of regular quaternion equations of the perturbed restricted spatial three-body problem. Some of them were obtained by Chelnokov [33][34] .
In this paper, we consider the regular quaternion equations of the two-body problem, free from the singularity (division by zero) caused by gravitational forces. Other singularities are caused by using angular variables, in particular, the Euler angles, in classical equations of celestial mechanics and astrodynamics.
In this relation, we should note that the efficiency of solving a number of problems in celestial mechanics and astrodynamics is increased by using the two-body problem equations, written in some rotating coordinate frames. In such equations, the variables are the values which characterize the form, size, and orientation of an instantaneous orbit or orbital plane of the body in interest, or the orientation of a rotating coordinate frame. Traditionally, the Euler angles or direction cosines are used in mechanics and astrodynamics as such variables. Using the Euler angles provides geometrical clarity to the description of orbital motion, but then the equations contain cumbersome trigonometric expressions and additional singularities, where the equations degenerate. For example, the widely used Newton-Euler equations for osculating elements (slowly changing variables) include the differential equations for angular elements, i.e., the longitude of ascending node, the inclination of an orbit, and the angular distance from the pericenter to the node. These equations are highly nonlinear and degenerate when the inclination of an instantaneous orbit of a body in interest is zero or 180 degrees. Using the direction cosines allows to eliminate mentioned singularities of motion equations of a body in interest, but leads to significant increase in the order of the system of motion equations and to the loss of geometrical clarity. The disadvantages of using the Euler angles and direction cosines can be avoided by using Euler parameters as the orientation parameters of the used rotating coordinate frame or the instantaneous orbit of a body in interest. In this case, it is convenient to use the Hamilton rotation quaternion, the components of which are the Euler parameters, to define the orientation of that coordinate frame or orbit of a body in interest. Consequently, the equations of orbital motion include the quaternion differential equation of angular motion of the used rotating coordinate frame or instantaneous orbit or orbital plane of a body in interest, which has a compact, symmetric, and non-degenerating structure.
Till now, we and a number of other researchers have proposed equations of celestial mechanics and astrodynamics that are regular in the indicated sense, obtained using the Euler parameters and Hamiltonian quaternions of rotation on the basis of the differential equations of the perturbed spatial problem of two bodies. We have obtained various new, regular in the above sense, differential equations of the perturbed spatial restricted three-body problem, written in various rotating coordinate frames; in non-holonomic (azimuthally free) coordinate frames as follows: M 0 X nh 0 Y nh 0 Z nh 0 and M 1 X nh 1 Y nh 1 Z nh 1 , the axes M 0 X nh 0 and M 1 X nh 1 of which are directed along radius vectors r 0 and r 1 of a point particle in interest M , which are directed from gravitating points M 0 and M 1 ; in orbital coordinate frames, M 0 X orb 0 Y orb 0 Z orb 0 and M 1 X orb 1 Y orb 1 Z orb 0 ; in orbital and ideal coordinate frames, M 0 X id 0 Y id 0 Z id 0 and M 1 X id 1 Y id 1 Z id 0 ; and in ideal coordinate frames.
Each of the derived regular equation systems consists of two equation groups, i.e., equations defining the motion of a point particle in interest M with negligible mass in the coordinate frame M 0 X 0 Y 0 Z 0 , which has its origin at the gravitating point M 0 , and equations defining the motion of this point in the coordinate frame M 1 X 1 Y 1 Z 1 , which has its origin at the gravitating point M 1 . The axes of these coordinate frames are parallel to the corresponding axes of the inertial coordinate frame.
In these equations, we use the Euler parameters and Hamilton rotation quaternions to describe the orientation of the used rotating coordinate frames. The equations include the corresponding quaternion kinematic equations of rotational motion of these coordinate frames in Euler parameters.
The independent variable in the equations is time t, and the dependent variables are Euler parameters and one of the four groups of the following variables: (i) the distances r 0 and r 1 from the point M to points M 0 and M 1 and the projections c 02 , c 03 and c 12 , c 13 of vectors c 0 and c 1 of velocity moments of the point M in relation to points M 0 and M 1 on the axes of non-holonomic coordinate frames, (ii) the distances r 0 and r 1 and moduli c 0 and c 1 of velocity moment vectors of the point particle M , (iii) the distances r 0 and r 1 , moduli c 0 and c 1 of velocity moment vectors of the point particle M , and the angular variables ϕ 0 and ϕ 1 , which characterize the rotation of the ideal coordinate frames in relation to orbital coordinate frames, (iv) the ideal rectangular Hansen coordinates X i , Y i , Z i = 0 (i = 0, 1) of the point particle M in ideal coordinate frames.
We also obtain the equations of the perturbed spatial restricted three-body problem in complex variables composed from Hansen coordinates and in complex Cayley-Klein parameters. These equations are of the smallest order. They allow the use of a non-positional number system of residual classes, which is self-correcting and which increases the accuracy of numerical solutions of differential equations of orbital motion.
From the above-described (constructed by us) equations of the three-body problem, we obtain various new regular quaternion differential equations of the perturbed spatial restricted three-body problem, written in a nonholonomic, orbital, and ideal coordinate frames. By regular equations, we note that the equations free of singularities are caused by gravitational forces.
Each of the obtained systems of regular equations consists of two groups of equations, i.e., the equations defining the motion of a point particle in interest M with negligible mass in the coordinate frame M 0 X 0 Y 0 Z 0 , which has its origin at the gravitating point particle M 0 , and the equations defining the motion of this point particle in the coordinate frame M 1 X 1 Y 1 Z 1 , which has its origin at the gravitating point particle M 1 . The axes of these coordinate frames are parallel to the corresponding axes of the inertial coordinate frame. The first one of these equation groups (i = 0) is used for studying the motion of the point particle M in the neighborhood of the point particle M 0 when the distances r 0 and r 1 satisfy the inequation m 1 r 2 0 m 0 r 2 1 , and the second one of these equation groups (i = 1) is used for studying the motion of the point particle M in the neighborhood of the point particle are the distances r 0 and r 1 , the moduli c 0 and c 1 of velocity moment vectors of the point particle M , the Keplerian energy h 0 and h 1 of motion of the point particle M in coordinate frames M 0 X 0 Y 0 Z 0 and M 1 X 1 Y 1 Z 1 , the Euler parameters or corresponding rotation quaternions, which characterize the orientation of orbital coordinate frames in the inertial coordinate frame, and also time t. These equation systems use the new independent variables τ 0 and τ 1 defined by dt = r 0 dτ 0 and dt = r 1 dτ 1 . Each system also includes the additional differential equation to correlate time τ 0 and τ 1 .
From the equations written in orbital coordinate frames, we derive the regular equations in which the new independent variables are the angular variables ϕ 0 and ϕ 1 , related to time t by differential relations dt = (r 2 0 /c 0 )dϕ 0 and dt = (r 2 1 /c 1 )dϕ 1 . The equations are also derived with the variables ρ 0 = 1/r 0 and ρ 1 = 1/r 1 , the moduli c 0 and c 1 of velocity moment vectors of the point particle M , the Euler parameters or corresponding rotation quaternions, which characterize the orientation of orbital coordinate frames in the inertial coordinate frame, and also time t.
The variables in the equation systems written in ideal coordinate frames M 0 X id 0 Y id 0 Z id 0 and M 1 X id 1 Y id 1 Z id 0 are the Levi-Civita variables U 00 , U 03 and U 11 , U 13 , the Keplerian energy h 0 and h 1 of motion of the point particle M in coordinate frames M 0 X 0 Y 0 Z 0 and M 1 X 1 Y 1 Z 1 , the Euler parameters or corresponding rotation quaternions, which characterize the orientation of ideal coordinate frames in the inertial coordinate frame, and also time t. The equations use the above-mentioned independent variables τ 0 and τ 1 . Each system includes the additional differential equation to correlate time τ 0 and τ 1 .
We also obtain new regular quaternion differential equations of the perturbed spatial restricted three-body problem, in which some of the variables are the total energy of motion of the point particle M in the coordinate frames M 0 X 0 Y 0 Z 0 and M 1 X 1 Y 1 Z 1 or the variables, which are the first Jacobi integrals for the case of the unperturbed spatial circular restricted three-body problem.

Conclusions
Our work analyzes the quaternion regularization methods and the regular quaternion models of celestial mechanics and astrodynamics, derived on the basis of differential equations of the perturbed spatial two-body problem and the perturbed central motion of a point particle using the KS variables and their modified versions, and also using the Levi-Civita variables and the Euler (Rodrigues-Hamilton) parameters. Some applications of these methods and models are also discussed.
In spite of the negative opinion of Stiefel and Scheifele on the efficiency of using the quaternion matrices for regularizing the two-body problem equations, quaternion matrices and quaternions proved to be convenient and efficient tools for regularizing the perturbed spatial two-body problem equations, for eliminating the singularities (divisions by zero) induced in the two-body problem equations by Newtonian gravitational forces.
Also, in spite of the existing opinion that the Levi-Civita transformation cannot be generalized to the case of three-dimensional space, it turned out that the Levi-Civita variables can be successfully used to derive the regular differential equations of the perturbed spatial two-body problem. We accomplish that by using the ideal rectangular Hansen coordinates, the regular Levi-Civita variables, and the orientation quaternion of the ideal coordinate frame, in which the original differential equations of the perturbed spatial two-body problem are written, and also by using the Keplerian energy as an additional variable and using the new independent variable. The regular equations of the perturbed spatial two-body problem, derived by us, include the equations in the Levi-Civita variables, which have the form of the regular Levi-Civita equations for the two-dimensional problem.
We also show that it is convenient to use the perturbed spatial two-body problem equations, written in the rotating coordinate frame, to efficiently derive regular models of celestial mechanics and astrodynamics. This naturally leads to the introduction of the four-dimensional Euler parameters and Hamilton quaternions in the orbital motion equations to define the orientation of the rotating coordinate frame, and to introduce with their help the four-dimensional KS variables and their modified versions. Therefore, quaternions, which are used widely and efficiently today to define the rotational (angular) motion of a solid body and to derive the analytically and numerically convenient kinematic and dynamic equations of rotational motion of a solid body, and which are also efficiently used in many technical applications (for example, for determining the orientation of various moving objects and for controlling their rotational motion), are naturally included in the dynamics of the orbital (trajectory) motion of a point particle.
The efficiency of applying the regular quaternion models of celestial mechanics and astrodynamics, based on the differential equations of the perturbed spatial two-body problem and the perturbed central motion of a point particle, is demonstrated in the problem of predicting the motion of celestial and cosmic objects, in the problem of studying the motion of the Earth's satellite, considering not only central, but also zonal, tesseral, and sectorial harmonics of the Earth's gravitational field, in the problem of optimal control of spacecraft orbital motion using the maximum principle, and also in the problem of autonomous inertial navigation in space.