As you probably know, I am working (slowly) on an astronomically accurate calendar generator. All the orbital calculations involved are quite challenging, and I am discovering a lot. It is a lot of fun (except for the all the times I need to do some trigonometric magic to make some formula work). Anyway, during this process, I am reshaping and producing many many formulas. I am sure that in six months I will forget all the motivations behind them. For this reason, I want to try to save some of them here. In this way, I will have a good place where to look back at my notes and, moreover, I can be useful to other people trying to do some low-accuracy orbital calculations. I want to start from the beginning: orbital period and orbital trajectory.
Position of the planet on the orbit
First of all, we need to have a clear idea of what are the orbital elements. If we consider a true three-object system (e.g., sun-earth-moon) things get very, very messy. A complete orbital diagram for such system involves an impossible number of angles, axises and intersections between imaginary lines. For this reason, for simplicity, we want to start from a very basic two-objects system: a point object orbiting around another (more massive) point object. In this situation, there are no rotation axises, there are no orbital inclination. There is only one object orbiting the other.
![A diagram showing the basic measures for an object P on an elliptic orbit around another object in F [source].](https://www.davideaversa.it/blog/computing-planetary-orbits-two-celestial-objects/Eccentric_and_True_Anomaly.png)
We also assume that the mass of the star is much greater than the planet mass. In this way we can consider the star “fixed” and put there our system origin. As we know, an orbiting system like this obeys to the Kepler’s laws of planetary motion:
- The orbit of a planet is an ellipse with the star is at one of the two foci (F in the above diagram).
- A line segment joining a planet and the star sweeps out equal areas during equal intervals of time.
- The square of the orbital period of a planet is proportional to the cube of the semi-major axis of its orbit.
We will look at Law 2 and 3 later, for now the first law is enough. Orbits are elliptical. This seems easy but keeping track of non-uniform movement along an ellipse is more complicated than you think. Assuming that you know how an ellipse is defined (eccentricity, major axis, minor axis, and so on) we are interested in one measure: the true anomaly. This is the angle between the planet P and the star in F. In the diagram, this is the angle f. As you can see, if we know this angle, and we know the orbit dimensions, we know where the planet is on the orbit. Because orbit dimensions are usually given (there are many sets of equivalent dimensions, for the moment, we don’t care), the only thing we need to determine the planet position over time is the true anomaly over time.
Computing the True Anomaly
Computing the true anomaly over time is not a trivial task. Before we can compute the true anomaly, we need some intermediate steps. The first stop is computing the mean anomaly. Imagine a circle with radius $$ a$ (the semi-axis major of the ellipse, in other words, half the “width” of an ellipse). Imagine now a body orbiting this circular orbit with constant speed (because circle) with the same period (\( P \) ) of the real planet. The angle of this imaginary body at time \( t \) is the mean anomaly. Computing the mean anomaly over time is very easy, after all, it is moving with constant speed. Therefore:
$$ M(t) = \frac{2\pi}{P}$$Now, we introduce a new angular measure, the eccentric anomaly. This is the angle $$ E$ shown in the above diagram. It represents the angle between the circular orbit center C and the kind of projection of the real planet position on the elliptic orbit on the imaginary circular orbit we defined before. To compute this value we use the Kepler’s Equation:
$$ M(t) = E(t) - \epsilon \sin(E(t))$$Note that \( \epsilon \) is just the ellipse eccentricity. This is a wonderful relation between the mean anomaly and the eccentric anomaly. However, give this formula, it is impossible to express \( E(t) \) in closed form. For this reason, there are a thousand of methods for approximating this value. The one I prefer is the recursive one. In short, we rewrite the formula as
$$ E(t) = M(t) + \epsilon \sin(E(t))$$You can note that we have expressed a kind of recursive function. Thus, we can replace \( E(t)\) on the right with the same formula again
$$ E(t) = M(t) + \epsilon \sin(M(t) + \epsilon \sin(E(T))$$And we can continue the process over and over, every time obtaining a closer approximation. I usually do this process 3 times, and, with some trigonometric magic trick, I get the following formula
$$ E = M + \epsilon \sin(M) + \epsilon^2 \sin(M)\cos(M) + \frac{1}{2} \epsilon^3 \sin(M) \left( 3 \cos^2(M) - 1) \right)$$Finally, we can rewrite this formula in a way that it is easier for the computer (it is faster, it performs less sin and cos, and reduce the floating-point errors):
$$ E = M + \epsilon \sin(M) \left(1 - \epsilon \left( \frac{\epsilon}{2} + \cos(M) \left(1 + \frac{3 \epsilon}{2} \cos(M) \right) \right) \right) $$Cool. Now it is time for the last step of our process. We have seen that the eccentric anomaly is just the “angle between the center of the circular orbit and the projection of the real planet on it” (this is a very loose definition, but I think it is the most intuitive approximation). Therefore, to find the real angle we just to “project the planet back” to the elliptical orbit. Fortunately, this is much less painful problem. The relation between the two angles is given by the following formula (where $$ \theta$ is the true anomaly).
$$ (1 - \epsilon) \tan^2\left(\frac{\theta}{2}\right) = (1 + \epsilon) \tan^2\left(\frac{E}{2}\right)$$This can be rewritten as
$$ \theta = 2 \, \mathop{\mathrm{arg}}\left(\sqrt{1-e} \, \cos\frac{E}{2} , \sqrt{1+e}\sin\frac{E}{2}\right)$$A formula that finally concludes our journey.
Conclusions
This is everything for now. I’ve put many concepts on the table and, before we move on, we need to be perfectly at ease with the tree angular orbit measurements, the anomalies. In the next articles we will see how to use them to compute some basic orbital events (such as equinoxes) using these formulae.
Bibliography
I have some book that was extremely useful in this task:
- Astronomical Algorithms by Jean Meeus. This book is a 20+ year old book for writing astronomical algorithms on a computer. For this reason the code examples are in BASIC. Yep. But the math and the core of the problems are all there. This book is still a gold mine if you are interested in this kind of works.
- Introduction to Geodetic Astronomy by T.B. Thompson. There is no link for this, but you can use your web search ability to find a PDF version, somewhere. This is an even older book (December 1981) but, at least, there is no obsolete code. The book is all about the celestial math governing the Earth and the Sun (and other orbiting spheres).
- If you want something more recent, there is Astronomy on the Personal Computer by Montenbruck and Pfleger, but I didn’t read it. I know many people suggested me this book and the examples involve a more recent C++, but it costs a lot for a book that I need for a side project.
