Euler’s Forgotten Equation of the Center

Sylvio R. Bistafa

University of São Paulo, São Paulo, Brazil.

DOI: 10.4236/ahs.2021.101005                

Abstract

In a 1778 publication in Latin, titled Nova Methodvs Motvm Planetarvm Determinandi (New method to determine the motion of planets), Euler derives an equation of the center, which, apparently, has been forgotten. In the present work, the developments that led to Euler’s equation of the center are revisited, and applied to three planets of the Solar System. These are then compared with results obtained from an equation of center that has been proposed, showing good agreement for planets with not so high eccentricities. Nonetheless, Euler’s derivation was not influential, and since then, the resulting equation of the center has been neglected by scholars and by specialized publications alike.

Keywords

Euler’s Works on Astronomy, Equation of the Center, History of Orbital Calculations, Investigations on Planets’ Orbits, Astronomical Calculations, Keplerian Orbital Mechanics

1. Introduction

Since antiquity, the problem of predicting the motions of the heavenly bodies has been simplified by reducing it to one of a single body in orbit about another. In calculating the position of the body around its orbit, it is often convenient to begin by assuming circular motion. This first approximation is then simply a constant angular rate multiplied by an amount of time. There are various methods of proceeding to correct the approximate circular position to that produced by elliptical motion, many of them complexes, and many involving solutions of Kepler's equation. In contrast, the equation of the center is one of the easiest methods to apply. In cases of small eccentricity, the position given by the equation of the center can be nearly as accurate as any other method of solving the problem. Many orbits of interest, such as those of bodies in the Solar System or of artificial Earth satellites, have these nearly-circular orbits. As eccentricity becomes greater, and orbits more elliptical, the equation's accuracy declines, failing completely at the highest values, hence it is not used for such orbits. The equation in its modern form can be truncated at any arbitrary level of accuracy, and when limited to just the most important terms, it can produce an easily calculated approximation of the true position when full accuracy is not important. The ancient Greeks, in particular Hipparchus, knew the equation of the center as prostaphaeresis, although their understanding of the geometry of the planets’ motion was not the same. It was specified and used by Kepler, as that variable quantity determined by calculation which must be added or subtracted from the mean motion to obtain the true motion. Equation of the center. (n.d.). In Wikipedia. Retrieved November 25, 2020, from https://en.wikipedia.org/wiki/Equation_of_the_center.

An orbiting body’s mean longitude is calculated $l=\Omega +\omega +M$, where Ω is the longitude of the ascending node, ω is the argument of the pericenter and Mis the mean anomaly, the body's angular distance from the pericenter as if it moved with constant speed rather than with the variable speed of an elliptical orbit (Figure 1). Its true longitude is calculated similarly, $L=\Omega +\omega +\nu$, where ν is the true anomaly. The equation that gives the difference between the true anomaly ν and the mean anomaly M is called the equation of the center, and is typically expressed a function of mean anomaly, M, and orbital eccentricity, e. Mean longitude. (n.d.). In Wikipedia. Retrieved November 25, 2020, from https://en.wikipedia.org/wiki/Mean_longitude.

Narrien (1833) traces the origin of the equation of the center to Hipparchus (190 BC-120 BC), who made extensive investigations of orbits of the Sun and the Moon, and to Ptolemy (100 - 170), who extended these works in the determination of the orbits of Venus and Mercury. Historically, graphical methods and various sets of tables exist giving the true anomaly for various eccentricities Roy (2005). A common practice is to develop the equation of the center as a series in mean anomaly whose coefficients are given in terms of powers of the eccentricity of the orbit. Fourier series, logarithmic series Smart (1956) and Bessel functions Brower and Clemence (1961) have been used to this end. Nowadays,

Figure 1. Orbital angular parameters.

modern computer power has replaced the need for such tables and procedures.

2. Euler’s Derivation of the Equation of the Center

Euler’s equation of the center was developed in E519—Nova methodus motum planetarum determinandi (New method to determine the motion of planets) Euler (1778) as follows. The Planet is supposed to move on the plane of the figure (Figure 2), where S is the center of the Sun, with the axis S♈ pointing to the first point of Aries. After the time τ, expressed in days, has elapsed, the Planet appears in P, from where the perpendicular PQ to the axis S♈ is drawn, thus defining the two coordinates $SQ=x$ and $QP=y$ ; and with the distance from the Planet to the Sun $SP=v$, which can be written as $v^2=x^2+y^2$.

From these, the principle of motion then gives these two equations:

$\frac{1}{\Delta }\frac{\text{dd}x}{\text{d}{\tau }^2}=-\frac{x}{v^3}\text{and}\frac{1}{\Delta }\frac{\text{dd}y}{\text{d}{\tau }^2}=-\frac{y}{v^3},$

where Δ is a constant quantity which will be soon defined from the motion of the Earth.

Let be drawn the line SM in the mean location in which the Planet, at the same time, is about to occupy; certainly that location can be easily obtained from tables for the mean motion. Let us then call its longitude, given by the angle ♈ $\hat{S}M=\zeta$, which is proportional to the time, with the segment SM and the axis being observed at this time, from which the location of the Planet is referred to by the orthogonal coordinates, $Sq=X$, $qP=Y$, and then $v^2=X^2+Y^2$ ; and from these, the previous coordinates are defined as

$x=X\cos \zeta -Y\sin \zeta \text{and}y=X\sin \zeta -Y\cos \zeta .$.

By differentiating twice the above two expressions, for $\text{d}\zeta$ constant, results in the following

$\text{dd}x\cos \zeta +\text{dd}y\sin \zeta =\text{dd}X-2\text{d}\zeta \text{d}Y-X\text{d}{\zeta }^3,$

$\text{dd}y\cos \zeta +\text{dd}x\sin \zeta =\text{dd}Y+2\text{d}\zeta \text{d}X-Y\text{d}{\zeta }^3.$

From the fundamental equations we have that

$\frac{\text{dd}x\cos \zeta +\text{dd}y\sin \zeta }{\Delta \text{d}{\tau }^2}=\frac{-x\cos \zeta -y\sin \zeta }{v^3}=-\frac{-X}{v^3},$

$\frac{\text{dd}y\cos \zeta +\text{dd}x\sin \zeta }{\Delta \text{d}{\tau }^2}=\frac{-y\cos \zeta +x\sin \zeta }{v^3}=\frac{-Y}{v^3},$

Figure 2. Sun S and Planet P locations, with the axis S♈ pointing to the first point of Aries.

which, in face of the above results, can be expressed as functions of X and Y as

$\frac{\text{dd}X\cos \zeta -2\text{d}\zeta \text{d}Y-X\text{d}{\zeta }^2}{\Delta \text{d}{\tau }^2}=\frac{-X}{v^3},$

$\frac{\text{dd}Y\cos \zeta -2\text{d}\zeta \text{d}X-Y\text{d}{\zeta }^2}{\Delta \text{d}{\tau }^2}=\frac{-Y}{v^3},$

and these two equations allow the determination of the new coordinates X and Y.

For the determination of the quantity Δ, let us assume that the Planet P is the Earth, and that the mean distance from the Earth to the Sun = 1. And considering that the mean motion of the Earth is circular, then $X=1$ and $Y=0$, therefore, $v=1$, and then

$-\frac{\text{d}{\zeta }^2}{\Delta \text{d}{\tau }^2}=-1\text{and}0=0.$

Therefore,

$\text{d}{\zeta }^2=\Delta \text{d}{\tau }^2,\text{and}\text{then}\Delta =\frac{\text{d}{\zeta }^2}{\text{d}{\tau }^2},$

where ζ indicates the mean longitude of the Earth, which, if for the time τ is expressed in days, we call it t, then, $\Delta =\frac{\text{d}{t}^2}{\text{d}{\tau }^2}$, such that in our formulas, in place of $\Delta \text{d}{\tau }^2$, we may write $\text{d}{t}^2$. Since now, the time τ corresponds to the mean motion of the Earth, then for the time of one day $t=59.8\text{'},19\text{'}\text{'}$, and once this value is taken into account, then, for any Planet, we have the following equations

$\frac{\text{dd}X\cos \zeta -2\text{d}\zeta \text{d}Y-X\text{d}{\zeta }^2}{\text{d}{t}^2}=\frac{-X}{v^3},$

$\frac{\text{dd}Y\cos \zeta -2\text{d}\zeta \text{d}X-Y\text{d}{\zeta }^2}{\text{d}{t}^2}=\frac{-Y}{v^3}.$

Since ζ indicates the medium longitude of any Planet, it certainly holds a relation with the angle t, and if we put $\zeta =nt$, such that $\zeta =n\text{d}t$, our equations assume the following forms

$\frac{\text{dd}X\cos \zeta }{\text{d}{t}^2}-\frac{2n\text{d}Y}{\text{d}t}-n^{2}X=\frac{-X}{v^3},$

$\frac{\text{dd}Y\cos \zeta }{\text{d}{t}^2}-\frac{2n\text{d}X}{\text{d}t}-n^{2}Y=\frac{-Y}{v^3}.$

Since $\zeta =nt$, and knowing that for the time of one year $t={360}^{\circ }$, then $\zeta =n\cdot {360}^{\circ }$, which is the mean motion of the Planet for the time of one year, then for λ years, the mean motion of the Planet will be $\lambda n\cdot {360}^{\circ }$, whence, the Planet whole revolution is completed when $\lambda n\cdot {360}^{\circ }=360$, thus $\lambda =\frac{1}{n}$, such that the periodic time of the Planet will be $=\frac{1}{n}$ years.

Let us now also consider that the mean distance from our planet to the Sun =a, and let us concede a mean motion to this Planet, with no eccentricity; then $Y=0$ and $X=v=a$, and in this case our formulas give $n^2=\frac{1}{a^3}$ and $0=0$ ; whence, since the periodic time of the Planet is known to be $=\frac{1}{n}$ years, now it will be $=a\sqrt{a}$ years. In this way, it exposes another Keplerian rule, which establishes that the ratios of the periodic times of the Planets are in the three halves power of the mean distances. Since indeed we have taken the mean distance of the Earth to the Sun = 1, because its periodic time is one year, by this rule, the periodic time of the Planet, whose mean distance to the Sun = a ought to be $a\sqrt{a}$ years.

Let be introduced now, the Planet’s orbit eccentricity, considering that it s not very pronounced, such that the real location of the Planet does not differ too much from the mean location; then the abscissa $Sq=X$ does not differ too much from the mean distance a, which we establish such that $X=a\left(1+x\right)$ and $qP=Y=ay$. It is observed that both x and y are new quantities and should not be confused with the coordinates x and y previously defined. And, recognizing that now $v=a\sqrt{{\left(1+x\right)}^2+y^2}$, results in the following two equations in terms of the new variables x and y

$\frac{\text{dd}x}{\text{d}{t}^2}-\frac{2n\text{d}y}{\text{d}t}-n^2\left(1+x\right)=\frac{-\left(1+x\right)}{a^3{\left[{\left(1+x\right)}^2+y^2\right]}^{3/2}},$

$\frac{\text{dd}y}{\text{d}{t}^2}-\frac{2n\text{d}x}{\text{d}t}-n^{2}y=\frac{-y}{a^3{\left[{\left(1+x\right)}^2+y^2\right]}^{3/2}},$

and recalling that $n^2=\frac{1}{a^3}$, then our equations are reduced to the following forms

$\frac{\text{dd}x}{\text{d}{t}^2}-\frac{2n\text{d}y}{\text{d}t}-n^2\left(1+x\right)=\frac{-n^2\left(1+x\right)}{{\left[{\left(1+x\right)}^2+y^2\right]}^{3/2}},$

$\frac{\text{dd}y}{\text{d}{t}^2}-\frac{2n\text{d}x}{\text{d}t}-n^{2}y=\frac{-n^{2}y}{{\left[{\left(1+x\right)}^2+y^2\right]}^{3/2}}.$

Recognizing that $y^2$ is very small when compared to ${\left(1+x\right)}^2$, the irrational expression in the second member can be expanded into a series, resulting in the following two equations

$\begin{array}{l}\frac{\text{dd}x}{\text{d}{t}^2}-\frac{2n\text{d}y}{\text{d}t}-n^2\left(1+x\right)\\ =\frac{-n^2\left(1+x\right)}{{\left(1+x\right)}^2}+\frac{3\cdot n^2{y}^2}{2\cdot {\left(1+x\right)}^4}-\frac{3\cdot 5\cdot n^2{y}^4}{2\cdot 4{\left(1+x\right)}^6}+\frac{3\cdot 5\cdot n^2{y}^6}{2\cdot 4\cdot 6{\left(1+x\right)}^8}etc,\end{array}$

$\begin{array}{l}\frac{\text{dd}y}{\text{d}{t}^2}-\frac{2n\text{d}x}{\text{d}t}-n^{2}y\\ =\frac{-n^{2}y}{{\left(1+x\right)}^3}+\frac{3\cdot n^2{y}^3}{2\cdot {\left(1+x\right)}^5}-\frac{3\cdot 5\cdot n^2{y}^5}{2\cdot 4{\left(1+x\right)}^7}+\frac{3\cdot 5\cdot 7\cdot n^2{y}^7}{2\cdot 4\cdot 6{\left(1+x\right)}^9}etc.\end{array}$

By expanding the irrational terms in the second members of these equations into series, and introducing back the mean longitude of the Planet ζ, such that $nt=\zeta$, $n\text{d}t=\text{d}\zeta$, and $n^2\text{d}{t}^2=\text{d}{\zeta }^2$, gives the following simpler equations

$\begin{array}{c}\frac{\text{dd}x}{\text{d}{\zeta }^2}-\frac{2\text{d}y}{\text{d}\zeta }=3x-3x^2+4x^3-5x^4+6x^5-7x^{6}etc.\\ +\frac{3}{2}{y}^2-6x{y}^2+15x^2{y}^2-30x^3{y}^2+\frac{105}{2}{x}^4{y}^{2}etc.\\ -\frac{15}{8}{y}^4+\frac{45}{4}x{y}^4-\frac{315}{8}{x}^2{y}^{4}etc.\\ +\frac{35}{16}{y}^{6}etc.\end{array}$

$\begin{array}{c}\frac{\text{dd}y}{\text{d}{\zeta }^2}-\frac{2\text{d}x}{\text{d}\zeta }=3xy-6x^{2y}+10x^{3}y-15x^{4}y+21x^{5}yetc.\\ +\frac{3}{2}{y}^3-\frac{15}{2}x{y}^3+\frac{45}{2}{x}^2{y}^3-\frac{105}{2}{x}^3{y}^{3}etc.\\ -\frac{15}{8}{y}^5+\frac{105}{8}x{y}^{5}etc.\end{array}$

It should be observed that if the orbit of the Planet lacks eccentricity, then $x=0$ and $y=0$. As long as these quantities are not zero, to what extent the eccentricity is present. Let us then consider that the eccentricity = e, and since we expect that it is rather small, it can be easily realized that both x and y can be expressed by the following convergent series

$x=eP+e^{2}Q+e^{3}R+e^{4}S+e^{5}Tetc.$

$y=ep+e^{2}q+e^{3}r+e^{4}s+e^{5}t+e^{6}uetc.$

Next, the various terms that appear in the equations are written in terms of the new variables $e,P,Q,R,S,T,p,q,r,s,t,u$, and substituted into these equations. Then, by considering that the unknown quantities $P,Q,R,S$ and $p,q,r,s$ should be independent of the eccentricity e, allows six pairs of independent differential equations to be constructed with powers of e varying from one to six. The first four pairs of differential equations up to power four in e are the following.

First Order in e (affected by e).

1) $\frac{\text{dd}P}{\text{d}{\zeta }^2}-\frac{2\text{d}p}{\text{d}\zeta }=3P,$

2) $\frac{\text{dd}p}{\text{d}{\zeta }^2}+\frac{2\text{d}P}{\text{d}\zeta }=0.$

Second Order in e (affected by e²).

1) $\frac{\text{dd}Q}{\text{d}{\zeta }^2}-\frac{2\text{d}q}{\text{d}\zeta }=3Q-3P^2+\frac{3}{2}{p}^2,$

2) $\frac{\text{dd}q}{\text{d}{\zeta }^2}+\frac{2\text{d}Q}{\text{d}\zeta }=3Pp.$

Third Order in e (affected by e³).

1) $\frac{\text{dd}R}{\text{d}{\zeta }^2}-\frac{2\text{d}r}{\text{d}\zeta }=3R-6PQ+4P^3+3pq-6p^{2}P,$

2) $\frac{\text{dd}r}{\text{d}{\zeta }^2}+\frac{2\text{d}R}{\text{d}\zeta }=3pQ+3qP-6P^2+\frac{3}{2}{p}^3.$

Fourth Order in e (affected by e⁴).

1) $\begin{array}{l}\frac{\text{dd}S}{\text{d}{\zeta }^2}-\frac{2\text{d}s}{\text{d}\zeta }=3S-6PR-3Q^2+3pr+\frac{3}{2}{q}^2-6p^{2}Q+12P^2\\ -12pqP-5P^4+15p^2{P}^2-\frac{15}{8}{p}^4,\end{array}$

2) $\frac{\text{dd}s}{\text{d}{\zeta }^2}+\frac{2\text{d}S}{\text{d}\zeta }=3pR+3qQ+3pP-12pPQ-6q{P}^2+\frac{9}{2}{p}^{2}q+10p{P}^3-\frac{15}{2}{p}^{3}P.$

The solutions of these equations are given by sine and cosine functions of a new variable θ, which is identified as the mean anomaly. The solutions are the following

$P=\cos \theta \text{and}p=-2\sin \theta ,$

$Q=-\frac{1}{2}+\frac{1}{2}\cos 2\theta \text{and}q=\frac{1}{4}\sin 2\theta ,$

$R=-\frac{3}{8}\cos 3\theta \text{and}r=+\frac{9}{8}\sin \theta -\frac{7}{24}\sin 3\theta ,$

$S=\frac{407}{8}-\frac{209}{24}\cos 2\theta +\frac{219}{965}\cos 4\theta \text{and}s=\frac{355}{48}\sin 2\theta +\frac{49}{480}\sin 4\theta .$

By summing up these solutions the complete solution for x and y are the following

$\begin{array}{c}x=e\cos \theta +e^2\left(\frac{1}{2}-\frac{1}{2}\cos 2\theta \right)-\frac{3}{8}{e}^3\cos 3\theta \\ +e^4\left(\frac{407}{8}-\frac{209}{24}\cos 2\theta +\frac{219}{965}\cos 4\theta \right),\end{array}$

$\begin{array}{c}y=2e\sin \theta +\frac{1}{4}{e}^2\sin e^{2}2\theta +e^3\left(\frac{9}{8}\sin \theta -\frac{7}{24}\sin 3\theta \right)\\ +e^4\left(\frac{355}{48}\sin 2\theta +\frac{49}{480}\sin 4\theta \right).\end{array}$

Having been found the values for x and y, the angle $M\hat{S}P$ becomes known, which it is called equation of the center ω, and then the true longitude of the Planet will be given by

♈ $\hat{S}P=\zeta +\omega ,$

where ζ is the mean longitude of the Planet, and $\omega ={\tan }^{-1}\frac{y}{1+x}$.

And the distance of the Planet to the Sun will be given by $SP=a\left(1+x\right)\sec \omega$, where a is the mean distance from the Planet to the Sun.

3. Application

In Keplerian motion, the coordinates of the body retrace the same values with

Figure 3. Comparative results of Euler’s equation of the center and an equation of the center obtained via Fourier expansion for three planets with different eccentricities.

each orbit, which is the definition of a periodic function. Such functions can be expressed as periodic series of any continuously increasing angular variable and the variable of most interest is the mean anomaly, M. Because it increases uniformly with time, expressing any other variable as a series in mean anomaly is essentially the same as expressing it in terms of time. Because the eccentricity, e, of the orbit is small in value, the coefficients of the series can be developed in terms of powers of e. Such is the case of an equation of the center obtained via Fourier expansion that has been proposed Smart (1956):

$\nu -M=\left(2e-\frac{1}{4}{e}^3\right)\sin M+\frac{5}{4}{e}^2\sin 2M+\frac{13}{12}{e}^3\sin 3M+O\left(e^4\right).$

This equation was compared to Euler’s equation of the center by applying them to three planets, namely: Venus ( $e=0.006773$ —low eccentricity), Mars ( $e=0.093405$ —medium eccentricity), and Mercury ( $e=0.205635$ —high eccentricity). The plot in Figure 3 shows that they are in good agreement for planets with not so high eccentricities.

4. Conclusion

Although Euler’s equation of the center is also given as series expansions, its derivation is based on a completely different approach, than the ones that we have mentioned here. This approach consists on applying the equation of motion (Newton’s 2^nd law) for a Planet attracted to the Sun by central forces inversely proportional to the square of the separation distance. The solutions of the resulting differential equations were then given by sine and cosine functions of mean anomaly, whose coefficients are powers of the eccentricity. Euler’s equation of the center gives results that compare favorably with other methods that have been proposed, at least for small eccentricities. Nonetheless, Euler’s derivation was not influential, and since then, the resulting equation of the center has been neglected by scholars and by specialized publications alike.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Brower, D., & Clemence, G. M. (1961). Methods of Celestial Mechanics. New York: Academic Press.
[2]	Euler, L. (1778). Nova methodus motum planetarum determinandi. Acta Academiae Scientarum Imperialis Petropolitinae, 1781, 277-302.
[3]	Narrien, J. (1833). An Historical Account of the Origin and Progress of Astronomy. London: Baldwin and Cradock.
[4]	Roy, A. E. (2005). Orbital Motion (4th ed.). Bristol: IOP Publishing Ltd. https://doi.org/10.1201/9781420056884
[5]	Smart, W. M. (1956). Textbook on Spherical Astronomy. London: Cambridge University Press.