Quasi-Analytical Solution of Kepler’s Equation as an Explicit Function of Time

1. Introduction

Kepler’s laws provide an empirical mathematical model to describe the dynamics of the motion of the planets around the Sun (later proven by Newton’s analytical results [1,2]. Therefore, they are of fundamental importance in astrophysics, celestial mechanics, and earth science [3]. In this way and under the Copernican hypothesis, where the Earth rotates on itself and around the Sun, phenomena such as the apparent motion of the stars and the flattening of the Earth, among others, could be explained. Nevertheless, Kepler’s fundamental problem is the motion of a planet over time as it traverses through its orbit. As we shall see, it is the so-called Kepler’s equation (KE) that provides the explicit time dependence of the fundamental quantities of the problem. Since this equation was first established by Kepler in 1609, until the present day, several solutions of different types have been proposed [2,4,5,6]. Due to its inherent nature, it has not yet been possible to obtain an exact explicit solution, and a large number of calculations must be made to improve the approximate solutions. It is hardly possible to refer here to all the relevant works. Some of the procedures give satisfactory approximations, but they are very complex and often require many terms in order to obtain a solution with a given degree of accuracy and may involve high computational cost [3,7]. These solutions normally consist of calculating infinite series, or high-order iterative approaches [3,7,8,9,10,12,14] and by integral solutions [13]. Recently, some methods which use artificial intelligence were proposed [11].

In this work, we present a quasi-analytical solution that gives the explicit equations of motion as a function of time. This solution is quasi-analytical because even though it gives explicit formulas, it depends on six numerical coefficients, which in turn depend on the eccentricity. As mentioned above, this result provides an explicit solution to KE.

1.1. The Orbit Equation

The classical problem begins then, with the description of the motion of the planets according to the following equations [15]:

$$t\left(\theta \right)=\frac{{\alpha }^2\mu }{l}{\int }_0^{\theta }\frac{1}{{(1+\epsilon cos\phi )}^2}d\phi ,$$

(1)

where

$$r\left(\theta \right)=\frac{\alpha }{1+\epsilon cos\theta },$$

(2)

has been used. Here, $(r,\theta )$ are the coordinates of the planet in the polar coordinate system, with the origin at the Sun, t is the time, $\epsilon$ is the eccentricity, and $\alpha$, $\mu$, and l are the focal parameter, the reduced mass, and the angular momentum, respectively [2]. The angle $\theta$ is called the true anomaly, defined as $\theta =0$ at the pericenter [14].

The inverse of the solution of Eq. (1), $\theta \left(t\right)$ is closely related to KE:

$$E-\epsilon sinE=M,$$

(3)

where E is the eccentric anomaly and M is the mean anomaly, expressed as:

$$M=\frac{2\pi }{T}t,t\in [0,T],$$

(4)

T being the orbital period. As was pointed out before, usually, E needs to be estimated by series expansions [16] or iterative methods [6]. Thus, the accuracy of the position and/or velocity of the celestial object moving in the Keplerian orbit, which may be obtained from the solution of KE, depends on the approximation method used. Figure 1 illustrates the variables involved in the problem.

Equations (1) and (3) are both transcendental. The solution of Eq. (1) gives the solution of the Eq. (3), and vice-versa, which has led to a vast literature concerning alternative approaches [3]. A detailed description of the different works presented in the literature on this problem can be found in [3,14].

The remainder of this paper is structured as follows. In Section 2 we present in detail, step by step, the quasi-analytical solution of Eq. (1). In Section 3 the results are presented and compared with the real orbit of each planet of the solar system. In Section 4 we present the solution to KE based on the solution of Eq. (1). Finally, in Section 5 we give a brief discussion and conclusions about the results obtained in this work.

2. Methodology

In the present section, the quasi-analytical method is established for solving Eq. (1). For clarity, we divide the method into sections that we will call steps. All calculations and graphics were performed with MATLAB 2020b. Now we briefly describe each step.

Step 0. In this step we integrate Eq. (1) and introduce a normalized time $\tau$ which maps the real time t from the interval $[0,T]$ to the interval $[0,\pi ]$, which leads us to Eq. (7) for $\tau \left(\theta \right)$.

Step 1. In this step we separate $\tau \left(\theta \right)$ into two parts, ${\tau }_a\left(\theta \right)$ and ${\tau }_b\left(\theta \right)$. Based on the observation that ${\tau }_a\left(\theta \right)$ resembles the behavior of $\tau \left(\theta \right)$ (especially for small values of $\epsilon$), and, taking into account its analytical form, we find the inverse of ${\tau }_a\left(\theta \right)$, ${\theta }_0\left(\tau \right)$ as our zeroth approximation for $\theta \left(\tau \right)$.

Step 2. Assuming the exact solution has the form of Eq. (14), we introduce the true argument of the arctangent function, $\phi \left(\tau \right)$, which gives the exact solution. Then, we introduce the function ${\psi }_0\left(\tau \right)$ (Eq. (17)) as the quotient of the true argument $\phi \left(\tau \right)$ and the argument of the zeroth approximation, ${\phi }_0\left(\tau \right)$. Finally, we define the first approximation ${\theta }_1\left(\tau \right)$ by adjusting the argument ${\phi }_0\left(\tau \right)$ of the arctangent function which appears in the zeroth approximation ${\theta }_0\left(\tau \right)$. The main idea is to make a series of adjustments in such a way that the quotient mentioned above approaches 1 as closely as possible.

Step 3. Following the logic of Step 2, we introduce the function ${\psi }_1\left(\tau \right)$ (Eq. (22)) and the argument ${\phi }_2\left(\tau \right)$ (Eq. (25)) in order to obtain the second approximation ${\theta }_2\left(\tau \right)$ (Eq. (26)). First, we propose two kinds of approximations for ${\psi }_1\left(\tau \right)$: linear and harmonic. This gives the approximations ${\theta }_{2.1}\left(\tau \right)$ and ${\theta }_{2.2}\left(\tau \right)$ (Eqs. (27) and (28)), respectively. This approximations give the corresponding precision of 3,400 km and 470 km, respectively, for the position of the Earth. This formulas are completely analytical, and may be used in some applications which do not require high precision.

Next, we proceed to improve the function ${\psi }_1\left(\tau \right)$ as an analytical expression (Eq. (29)), involving an arctangent function, with an argument which, in turn, is an expression that depends on six numerical coefficients (Eq. (30)). In order to calculate these coefficients, we propose two methods, namely, Method A and Method B.

Method A consists in calculating the coefficients for fixed values of the eccentricity, while Method B allows us to express these coefficients as an analytic function with other numeric coefficients calculated for ranges of the eccentricity.

2.1. Step 0

Our first step is integrating Eq. (1) to obtain $t\left(\theta \right)$. The solution is (see Appendix A):

$$t\left(\theta \right)=\frac{{\alpha }^2\mu }{l}\left[\frac{2}{{(1-{\epsilon }^2)}^{\frac{3}{2}}}arctan\left(\sqrt{\frac{1-\epsilon }{1+\epsilon }}tan\frac{\theta }{2}\right)-\frac{2\epsilon }{1-{\epsilon }^2}\frac{tan\frac{\theta }{2}}{1+\epsilon +(1-\epsilon ){tan}^2\frac{\theta }{2}}\right].$$

(5)

This solution is valid on the interval $\theta \in [-\pi ,\pi ]$. In order to make the solution valid on the interval $\theta \in (-\infty ,\infty )$, we need to apply the following correction:

$$t_1\left(\theta \right)=2t\left(\pi \right)⌊\frac{(1+(\theta -\pi \left)\right)}{2\pi }⌋+t\left(\theta \right).$$

(5a)

In order to avoid complications related with the periodicity of the functions of $\theta$ involved, and the monotony of time, in what follows, we focus our work only on the half-period of the orbit, namely, $\theta \in [0,\pi ]$. The second half of the orbit may be obtained from the first part using the symmetry of the orbit.

We now propose the following change of variable:

$$\tau =\frac{l}{{\alpha }^2\mu }\frac{{(1-{\epsilon }^2)}^{\frac{3}{2}}}{2}t.$$

(6)

So, Eq. (5), in terms of $\tau$ and $\theta$ has the form:

$$\begin{array}{c}\hfill \tau \left(\theta \right)=\left\{\begin{array}{c}arctan\left(\sqrt{\frac{1-\epsilon }{1+\epsilon }}tan\frac{\theta }{2}\right)-\epsilon \sqrt{1-{\epsilon }^2}\frac{tan\frac{\theta }{2}}{1+\epsilon +(1-\epsilon ){tan}^2\frac{\theta }{2}},\theta \in [0,\pi ),\hfill \\ \frac{\pi }{2},\theta =\pi .\hfill \end{array}\right.\end{array}$$

(7)

Let us seek for the solution to Eq. (7) in the form $\theta \left(\tau \right)$ on the interval $\theta \in [0,\pi ]$, and, as deduced from Eq. (7), in $\tau \in [0,\frac{\pi }{2}]$.

Note the relationship between $\tau$ and the mean anomaly M:

$$\tau =\frac{M}{2}.$$

(8)

2.2. Step 1

We rewrite Eq. (7) as follows:

$$\begin{array}{c}\hfill \tau \left(\theta \right)=\left\{\begin{array}{c}{\tau }_a\left(\theta \right)-{\tau }_b\left(\theta \right),\theta \in [0,\pi ),\hfill \\ \frac{\pi }{2},\theta =\pi ,\hfill \end{array}\right.\end{array}$$

(9)

where:

$${\tau }_a\left(\theta \right)=arctan\left(\sqrt{\frac{1-\epsilon }{1+\epsilon }}tan\frac{\theta }{2}\right),$$

(10)

$${\tau }_b\left(\theta \right)=\epsilon \sqrt{1-{\epsilon }^2}\frac{tan\frac{\theta }{2}}{1+\epsilon +(1-\epsilon ){tan}^2\frac{\theta }{2}}.$$

(11)

Figure 2 shows the behavior of $\tau \left(\theta \right)$, ${\tau }_a\left(\theta \right)$ and ${\tau }_b\left(\theta \right)$ with $\epsilon$ as parameter.

As we can see, Figure 2 shows how ${\tau }_a\left(\theta \right)$ is similar to $\tau \left(\theta \right)$, therefore, as our zeroth approximation ${\theta }_0\left(\tau \right)$, we will solve the following equation:

$$\tau \left(\theta \right)={\tau }_a\left(\theta \right)=arctan\left(\sqrt{\frac{1-\epsilon }{1+\epsilon }}tan\frac{\theta }{2}\right).$$

(12)

To solve the Eq. (12) we apply the tangent function to both parts. Solving for $\theta$, we obtain:

$$\begin{array}{c}\hfill {\theta }_0\left(\tau \right)=\left\{\begin{array}{c}2arctan\left(\sqrt{\frac{1+\epsilon }{1-\epsilon }}tan\tau \right)\tau \in [0,\frac{\pi }{2}),\hfill \\ \pi ,\tau =\frac{\pi }{2}.\hfill \end{array}\right.\end{array}$$

(13)

Figure 3 shows plots of $\theta \left(\tau \right)$ and ${\theta }_0\left(\tau \right)$ with $\epsilon$ as parameter. $\theta \left(\tau \right)$ was obtained using Eq. (7) and interchanging axes.

2.3. Step 2

We introduce the functions $\phi \left(\tau \right)$ and ${\phi }_0\left(\tau \right)$ as follows:

$$\theta \left(\tau \right)=2arctan\left[\phi \right(\tau \left)\right],$$

(14)

$${\phi }_0\left(\tau \right)=\sqrt{\frac{1+\epsilon }{1-\epsilon }}tan\tau ,$$

(15)

so that:

$${\theta }_0\left(\tau \right)=2arctan\left[{\phi }_0\left(\tau \right)\right].$$

(16)

Now, we introduce another function ${\psi }_0\left(\tau \right)$, defined on the open interval $\tau \in (0,\frac{\pi }{2})$:

$${\psi }_0\left(\tau \right)=\frac{\phi \left(\tau \right)}{{\phi }_0\left(\tau \right)}.$$

(17)

Appendix B proves the following limits (the former is used immediately to introduce the function ${\phi }_1\left(\tau \right)$, and the latter is used in step 3):

$$\underset{\tau \to 0}{lim}{\psi }_0\left(\tau \right)=\frac{1}{1-\epsilon },$$

(18)

and

$$\underset{\tau \to \frac{\pi }{2}}{lim}{\psi }_0\left(\tau \right)=1+\epsilon .$$

(19)

Now, we introduce, using Eq. (15) and Eq. (18), the function ${\phi }_1\left(\tau \right)$ as follows:

$${\phi }_1\left(\tau \right)=\underset{\tau \to 0}{lim}{\psi }_0\left(\tau \right){\phi }_0\left(\tau \right)=\frac{1}{1-\epsilon }{\phi }_0\left(\tau \right)=\frac{\sqrt{1+\epsilon }}{{(1-\epsilon )}^{\frac{3}{2}}}tan\tau .$$

(20)

We define the approximation one, ${\theta }_1\left(\tau \right)$ as follows:

$$\begin{array}{c}\hfill {\theta }_1\left(\tau \right)=\left\{\begin{array}{c}2arctan\left[{\phi }_1\left(\tau \right)\right]=2arctan\left(\frac{\sqrt{1+\epsilon }}{{(1-\epsilon )}^{\frac{3}{2}}}tan\tau \right),\tau \in (0,\frac{\pi }{2}),\hfill \\ 0,\tau =0,\hfill \\ \pi ,\tau =\frac{\pi }{2}.\hfill \end{array}\right.\end{array}$$

(21)

Figure 4 shows plots of $\theta \left(\tau \right)$ and ${\theta }_1\left(\tau \right)$ with $\epsilon$ as parameter.

In Figure 4 we observe an evident improvement of the ${\theta }_1\left(\tau \right)$ approximation compared to ${\theta }_0\left(\tau \right)$ shown in Figure 3.

2.4. Step 3

We introduce a new function, ${\psi }_1\left(\tau \right)$, defined on the open interval $\tau \in (0,\frac{\pi }{2})$:

$${\psi }_1\left(\tau \right)=\frac{\phi \left(\tau \right)}{{\phi }_1\left(\tau \right)}.$$

(22)

As followed from Eqs. (18), (19) and (20):

$$\underset{\tau \to 0}{lim}{\psi }_1\left(\tau \right)=1,$$

(23)

$$\underset{\tau \to \frac{\pi }{2}}{lim}{\psi }_1\left(\tau \right)=1-{\epsilon }^2.$$

(24)

Figure 5 shows the plot of ${\psi }_1\left(\tau \right)$ with $\epsilon$ as parameter.

We introduce the function ${\phi }_2\left(\tau \right)$ as follows:

$${\phi }_2\left(\tau \right)={\psi }_1\left(\tau \right){\phi }_1\left(\tau \right)={\psi }_1\left(\tau \right)\frac{\sqrt{1+\epsilon }}{{(1-\epsilon )}^{\frac{3}{2}}}tan\tau .$$

(25)

The definition of ${\phi }_2\left(\tau \right)$ as the product of ${\psi }_1\left(\tau \right)$ and ${\phi }_1\left(\tau \right)$ leads us to an apparent logical paradox, since, as it follows from Eq. (22), ${\phi }_2\left(\tau \right)=\phi \left(\tau \right)$. The solution to this paradox is the fact that ${\psi }_1\left(\tau \right)$ is an unknown function, which we are trying to approximate with the best precision possible. The existence of such a function is obvious, for it can be theoretically constructed point-wise. However, in this work we are looking for a quasi-analytic expression for ${\psi }_1\left(\tau \right)$, since the exact analytic expression for it, in principle, may not even exist.

Now, we define the approximation two, ${\theta }_2\left(\tau \right)$ as follows:

$$\begin{array}{c}\hfill {\theta }_2\left(\tau \right)=\left\{\begin{array}{c}2arctan\left[{\phi }_2\left(\tau \right)\right]=2arctan\left({\psi }_1\left(\tau \right)\frac{\sqrt{1+\epsilon }}{{(1-\epsilon )}^{\frac{3}{2}}}tan\tau \right),\tau \in (0,\frac{\pi }{2}),\hfill \\ 0,\tau =0,\hfill \\ \pi ,\tau =\frac{\pi }{2}.\hfill \end{array}\right.\end{array}$$

(26)

It should be noted that the approximation ${\theta }_1\left(\tau \right)$ (Eq. (21)) can be considered as the approximation ${\theta }_2\left(\tau \right)$ with ${\psi }_1\left(\tau \right)=1$. For small values of $\epsilon$, such as Earth’s ($\epsilon =0.0167$) ${\psi }_1\left(\tau \right)$ has a small range (for the Earth, it will be of [0.9997, 1]) so that the approximation ${\psi }_1\left(\tau \right)=1$ already gives us an acceptable result for certain applications. Thus, for the Earth, the maximum error of ${\theta }_1\left(\tau \right)$ is of $1.8\times {10}^{-4}$ rad, which, translated to the error in the position, taking into account the average distance from the Earth to the Sun as $1.5\times {10}^8$ km, is equivalent to 27,000 km, which is a little more than two diameters of the Earth. Now, approximating ${\psi }_1\left(\tau \right)$ in a linear form, and, as a cosine with period $\pi$ and amplitude adjusted to the range [$1-{\epsilon }^2$, 1] the following corresponding approximations are obtained:

$$\begin{array}{c}\hfill {\theta }_{2.1}\left(\tau \right)=\left\{\begin{array}{c}2arctan\left[\left(1-\frac{2{\epsilon }^2}{\pi }\tau \right)\frac{\sqrt{1+\epsilon }}{{(1-\epsilon )}^{\frac{3}{2}}}tan\tau \right],\tau \in [0,\frac{\pi }{2}),\hfill \\ \pi ,\tau =\frac{\pi }{2},\hfill \end{array}\right.\end{array}$$

(27)

$$\begin{array}{c}\hfill {\theta }_{2.2}\left(\tau \right)=\left\{\begin{array}{c}2arctan\left[\left(1+\frac{{\epsilon }^2}{2}(cos2\tau -1)\right)\frac{\sqrt{1+\epsilon }}{{(1-\epsilon )}^{\frac{3}{2}}}tan\tau \right],\tau \in [0,\frac{\pi }{2}),\hfill \\ \pi ,\tau =\frac{\pi }{2}.\hfill \end{array}\right.\end{array}$$

(28)

The approximation ${\theta }_{2.1}\left(\tau \right)$ gives a maximum error of $2.24\times {10}^{-5}$ rad, which translates to 3400 km in the position (slightly more than half the Earth’s radius), and the approximation ${\theta }_{2.2}\left(\tau \right)$ gives a maximum error of $3.11\times {10}^{-6}$ rad, which corresponds to approximately 470 km in the position.

In Figure 5, we notice that the shape of the curves resembles an arctangent function, with a mapping of the argument from ($-\infty$, ∞) to ($\frac{\pi }{2}$,0), with the range mapping from ($-\frac{\pi }{2}$, $\frac{\pi }{2}$) to ($1-{\epsilon }^2,1$), and with the behavior of the argument being asymmetric and nonlinear. By the aforementioned, we search for the approximation of the function ${\psi }_1\left(\tau \right)$ in the following form:

$${\psi }_1\left(\tau \right)=1+\frac{{\epsilon }^2}{2}\left(\frac{2}{\pi }arctan\left[\xi (\epsilon ,\tau )\right]-1\right),$$

(29)

where:

$$\xi (\epsilon ,\tau )=a_1{\tau }^{-2}+a_2{\tau }^{-1}+a_3\tau +b_1{\left(\tau -\frac{\pi }{2}\right)}^{-2}+b_2{\left(\tau -\frac{\pi }{2}\right)}^{-1}+b_3\left(\tau -\frac{\pi }{2}\right),$$

(30)

with $a_i$ y $b_i$ ($i=1,2,3$) being functions of $\epsilon$.

In order to evaluate the coefficients $a_i$ and $b_i$ (the six numerical coefficients mentioned in the introduction), two methods were developed, Method A and Method B. In Method A, the coefficients are calculated for each value of $\epsilon$, while in Method B, they are expressed as analytic functions which, in turn, depend on other numeric coefficients, that are calculated for ranges of $\epsilon$. As will be seen further, Method A is more accurate, but less general, while Method B is slightly less accurate (especially after a certain value of $\epsilon$), but more general and easier to use in applications.

2.4.1. Method A

The coefficients $a_i\left(\epsilon \right)$ and $b_i\left(\epsilon \right)$ were calculated numerically for each eccentricity value, corresponding to every planet of the solar system, by means of RMSE (Root Mean Square Error) minimization, using MATLAB 2020b.

Figure 6 shows the behavior of $a_i\left(\epsilon \right)$ and $b_i\left(\epsilon \right)$ in the range of $\epsilon$∈ (0, 0.5) (this range was selected because some absolute values of $a_i\left(\epsilon \right)$ and $b_i\left(\epsilon \right)$ grow drastically after $\epsilon =0.5$, which does not allow us to appreciate the behavior before $\epsilon =0.5$).

Table 1 shows the values of $a_i\left(\epsilon \right)$ and $b_i\left(\epsilon \right)$, as well as the error of $\theta \left(\tau \right)$ for the Earth and Pluto, in radians, where ME is the Maximum absolute Error, MAE is the Mean Absolute Error and RMSE is the Root Mean Square Error.

Since Table 1 serves just for illustrative purposes, the values of $a_i\left(\epsilon \right)$ and $b_i\left(\epsilon \right)$ were rounded to three decimal places. A detailed table including more accurate values for all planets will be shown in the section of Results. For now let us mention that the absolute maximum error in the position of the Earth is 915 meters and approximately 39,000 kilometers (close to the value of the circumference of the Earth, or 16 Pluto diameters) for the position of Pluto (the average distance of Pluto to the Sun is approximately 5.9 billion kilometers).

2.4.2. Method B

Analyzing Figure 6, it can be seen that the shapes of the curves have a behavior of a cubic polynomial, so we search for the coefficients $a_i\left(\epsilon \right)$ and $b_i\left(\epsilon \right)$ in the following form:

$$\begin{array}{c}\hfill a_i\left(\epsilon \right)=a_{ij}{\epsilon }^j,\\ \hfill b_i\left(\epsilon \right)=b_{ij}{\epsilon }^j,\\ \hfill i=1,2,3,\\ \hfill j=0,1,2,3,\end{array}$$

(31)

where repeated indices imply summation.

We introduce the approximation ${\theta }_3\left(\tau \right)$ with the same form of ${\theta }_2\left(\tau \right)$ given by Eqs. (26), (29), (30), where now Eq. (31) is used to approximate the $a_i$ and $b_i$ in Eq. (30).

The attempt to calculate the coefficients $a_{ij}$ and $b_{ij}$ in the complete range of $\epsilon$∈ (0,1) was not successful. Thus, it was decided to divide the range into five intervals, $\epsilon$∈ (0,0.1], $\epsilon$∈ (0.1,0.25], $\epsilon$∈ (0.25,0.5], $\epsilon$∈ (0.5,0.7] and $\epsilon$∈ (0.7,1). In the next section, we will present the values of $a_{ij}$ and $b_{ij}$ by ranges, and compare the precision of the approximation ${\theta }_2\left(\tau \right)$, which uses values of ${\psi }_1\left(\tau \right)$ calculated using Method A, with that of ${\theta }_3\left(\tau \right)$.

3. Results

We now compare the results of our quasi-analytical solution with the real orbits of the planets of the solar system. First, we resume the results obtained above in the following equations. By the change of variable:

$$\tau =\frac{1}{{\alpha }^2\mu }\frac{{(1-{ϵ}^2)}^{\frac{3}{2}}}{2}t,$$

(32)

the equations of motion in the polar coordinate system with the origin at the Sun, and in the interval of time $\tau$∈$[0,\frac{\pi }{2}]$$(t\in [0,\frac{T}{2}])$, $\theta$∈$[0,\pi ]$ are the following:

$$\begin{array}{c}\hfill \theta \left(\tau \right)=\left\{\begin{array}{c}2arctan\left[\left(1+\frac{{\epsilon }^2}{2}(\frac{2}{\pi }arctan\left[\xi (\epsilon ,\tau )\right]-1)\right)\frac{\sqrt{1+\epsilon }}{{(1-\epsilon )}^{\frac{3}{2}}}tan\tau \right],\tau \in (0,\frac{\pi }{2}),\hfill \\ 0,\tau =0,\hfill \\ \pi ,\tau =\frac{\pi }{2},\hfill \end{array}\right.\end{array}$$

(33)

$$r\left(\tau \right)=\frac{\alpha }{1+\epsilon cos\theta \left(\tau \right)},$$

(34)

where $\xi (\epsilon ,\tau )$ is given by the following expression:

$$\xi (\epsilon ,\tau )=a_1{\tau }^{-2}+a_2{\tau }^{-1}+a_3\tau +b_1{\left(\tau -\frac{\pi }{2}\right)}^{-2}+b_2{\left(\tau -\frac{\pi }{2}\right)}^{-1}+b_3\left(\tau -\frac{\pi }{2}\right),$$

(35)

being $a_i$ and $b_i$$(i=1,2,3)$ functions of $\epsilon$. In what follows we present the two methods developed for the calculation of the coefficients $a_i$ and $b_i$.

3.1. Method A

As mentioned above, in this method, the values of $a_i$ and $b_i$ were calculated for each eccentricity by means of RMSE minimization in MATLAB 2020b. Table 2 shows the values of $a_i$ and $b_i$, ME, MAE, and RMSE errors in $\theta \left(\tau \right)$, in radians, $D_s$ being the average distance of the planet to the Sun in kilometers and $E_p$ is the maximum error of the position in kilometers with respect to the real orbit.

The results of Table 2 show that for small values of $\epsilon$, the accuracy is high. However, for values of $\epsilon \sim 0.25$ the error increases significantly.

3.2. Method B

As was pointed out before, in this method we calculate the coefficients of a third-degree polynomial (Eq. (31)) by minimizing the RMSE, for the ranges of $\epsilon$. So, in this case, it is not necessary to calculate the $a_i$ and $b_i$ for each value of the eccentricity. Table 3 shows the values of the coefficients $a_{ij}$ and $b_{ij}$ for the different ranges of $\epsilon$. Obviously, using this method results in errors greater than those of Method A.

Figure 7 shows the comparison of the errors between methods A and B. As can be seen, the error in Method B is slightly higher (especially for ME) up to values of $\epsilon =0.8$. For values of $\epsilon >0.8$ the errors of Method B increase significantly with respect to those of Method A.

3.3. Equations of Motion for the Planets

Finally, let us express the equations of motion of both position and velocity in the Cartesian coordinate system, with the origin at the Sun.

The equations for the position are:

$$\begin{array}{c}\hfill x\left(\tau \right)=r\left(\tau \right)cos\theta \left(\tau \right),\\ \hfill y\left(\tau \right)=r\left(\tau \right)sin\theta \left(\tau \right),\end{array}$$

(36)

where $\theta \left(\tau \right)$ and $r\left(\tau \right)$ are given by Eqs. (33) and (34), respectively.

The equations for the velocity are:

$$\begin{array}{c}\hfill v_x=\dot{r}\left(\tau \right)cos\theta \left(\tau \right)-r\dot{\theta }\left(\tau \right)sin\theta \left(\tau \right),\\ \hfill v_y=\dot{r}\left(\tau \right)sin\theta \left(\tau \right)+r\dot{\theta }\left(\tau \right)cos\theta \left(\tau \right),\end{array}$$

(37)

where:

$$\begin{array}{c}\hfill \dot{\theta }\left(\tau \right)=2\frac{\sqrt{1+\epsilon }}{{(1-\epsilon )}^{\frac{3}{2}}}\frac{\dot{\psi }\left(\tau \right)tan\tau +\psi \left(\tau \right){sec}^2\tau }{1+\frac{1+\epsilon }{{(1-\epsilon )}^3}{\psi }^2\left(\tau \right){tan}^2\tau },\\ \hfill \dot{r}\left(\tau \right)=\frac{\alpha \epsilon \dot{\theta }\left(\tau \right)sin\theta \left(\tau \right)}{{(1+\epsilon cos\theta \left(\tau \right))}^2},\end{array}$$

(38)

and

$$\begin{array}{c}\hfill \psi \left(\tau \right)=\left\{\begin{array}{c}1+\frac{{\epsilon }^2}{2}(\frac{2}{\pi }arctan\left[\xi (\epsilon ,\tau )\right]-1),\tau \in (0,\frac{\pi }{2}),\hfill \\ 1,\tau =0,\hfill \\ 1-,{\epsilon }^2,\tau =\frac{\pi }{2},\hfill \end{array}\right.\end{array}$$

(39)

$$\begin{array}{c}\hfill \dot{\psi }\left(\tau \right)=\left\{\begin{array}{c}\frac{{\epsilon }^2}{\pi }\frac{\dot{\xi }\left(\tau \right)}{1+{\xi }^2\left(\tau \right)},\tau \in (0,\frac{\pi }{2}),\hfill \\ 0,\tau =0,\hfill \\ 0-,\tau =\frac{\pi }{2}.\hfill \end{array}\right.\end{array}$$

(40)

The function $\xi \left(\tau \right)$ defined in the interval $\tau$∈$(0,\frac{\pi }{2})$ is given by Eq. (35) and its derivative is the following:

$$\dot{\xi }\left(\tau \right)=-2a_1{\tau }^{-3}-a_2{\tau }^{-2}+a_3-2b_1{\left(\tau -\frac{\pi }{2}\right)}^{-3}-b_2{\left(\tau -\frac{\pi }{2}\right)}^{-2}+b_3.$$

(41)

Figure 8 shows the comparison plot between the real orbit, calculated numerically using Eq. (7), and the orbit obtained by numerical integration of Eq. (37) (where Method A was used). Additionally, the velocity vectors are shown at an arbitrary scale at some points of the trajectory.

As can be seen in Figure 8, the integrated orbit closely resembles the behavior of the real orbit. This shows that our Method A is sufficiently accurate even for eccentricities of the order up to $\epsilon =0.5$.

4. Solution to Kepler’s Equation

As can be seen from Figure 1b, it is trivial to obtain the following expression, which gives the formula for E (the eccentric anomaly), appearing in Eq. (3) as a function of the polar coordinates $(\rho ,\theta )$ with the center in the focus F. In what follows, we set the major semi-axis equal to one ($a=1$):

$$E=arccos(F+\rho cos\theta ),$$

(42)

where $F=\left|OF\right|$ (Figure 1b). As is well known from the theory of conic sections, $F=\epsilon a=\epsilon$, since we set $a=1$. Thus, Eq. (42) acquires the following form:

$$E=arccos(\epsilon +\rho cos\theta ).$$

(43)

Due to the normalization $a=1$, the expression for $\rho$ is as follows:

$$\rho =\frac{r}{a}.$$

(44)

The functions $\theta \left(\tau \right)$ and $r\left(\tau \right)$ are given by our final quasi-analytical solution (Eqs. (33) and (34)).

Now, in order to obtain the final explicit expression for the solution of KE, $E\left(M\right)$, we need to express $\theta$ and $\rho$ as functions of M (the mean anomaly given by Eq. (4)). So far, $\theta$ and r are functions of $\tau$ (where $\tau$ is given by Eq. (6)). Thus, by expressing $\tau$ as a function of M, $\theta$ and r automatically become functions of M. The expression for $\tau \left(M\right)$ is the following:

$$\tau \left(M\right)=\frac{M}{2}.$$

(45)

In order to obtain $\rho$ as a function of M, we use Kepler’s third law to express a as follows:

$$a={\left(\frac{l}{2\pi \mu \sqrt{1-{\epsilon }^2}}T\right)}^{\frac{1}{2}}.$$

(46)

The final solution to KE is given by Eq. (43), using the quasi-analytical solution for the motion given by Eqs. (33) and (34) and taking into account the relations (44), (45) and (46).

5. Remarks and Conclusions

In this work, we obtained a quasi-analytical solution for the motion of the celestial bodies as an explicit function of time in four steps. We called this solution quasi-analytical, due to dependency on some numerical coefficients, which, in turn, themselves depend on the eccentricity. We proposed two methods for evaluating these coefficients, Method A and Method B. The former gives a higher degree of accuracy, but involves calculations for each value of the eccentricity, while the ladder is less accurate (especially for values of $\epsilon >0.8$), but more practical, since it works for ranges of the eccentricity. All of our calculations involve low computational cost, and provide a high degree of accuracy for the motion of celestial bodies. Using our results, we obtained an explicit solution to KE without the need of calculating infinite series, or using numeric iterative methods.