My bonding period working for GSoC’18 poliastro project

The community bonding period is ending soon and I would like to tell what I’ve managed to achieve and contribute so far.

My major at the university is physics and my research is devoted to supercomputer simulations of strongly-coupled systems such as quark-gluon plasma (the one that people study at colliders) and graphene (the well-known one-atom-thick carbon material that has fantastic conductivity). I am a PhD student now and will get my degree soon. During my research at the lab, I became quite skillful in programming and calculus, that’s why I’ve chosen poliastro to be my organisation during this summer. I decided that this project will utilize all my skills. In addition, in school I was very interested in astronomy and I even wanted to become an astronomer, so this project will somewhat my dreams come true :)

The poliastro project, part of the OpenAstronomy community, aims at writing purely pythonic high-performance code for orbital movement simulations. In short, the idea is to predict the position and velocity of any object at any time, given its initial position and velocity. But the project is much more!

I could now try to create a well-connected and linked text about my contributions, but I will better describe all the pull requests one-by-one, and a reader will understand the structure of poliastro :)

Fixing an issue with non-smoothness in orbit plotting

In poliastro, people plot orbits to have a look at them, that’s simple. The idea is that the more points choose on the orbit and connect them with lines, the more smooth the final orbit will look. One could say that 150 points is enough to plot a smooth ellipse, however, not always. When the perigee distance of the orbit is much smaller than the apogee distance, the velocity at the perigee is very high.

Previously, the poliastro sampled from the true anomaly, that is, the same angles intervals between consequent points plotted on the ellipse. That caused the following issue. One can see that the orbit is very non-smooth at the perigee:

Blue orbit at the perigee is non-smooth, people don’t like non-smoothness.

The issue was resolved in this pull request. One just need to uniformly sample not true anomaly (the angles, you remember), but eccentric anomaly. The eccentric anomaly, by definition, reflects non-uniformity of velocity on orbit and varies faster as the body passes perigee. So, automatically, the algorithm would sample more points at the perigee and less at the apogee.

So, happy ending: without changing the number of sampled points (that would decrease performance of course) we made orbits smooth again.

Adding new analytical propagator

The humanity was interested in shapes of orbits long time ago. At first, people thought that these were circles. Then, when people discovered that the orbit is not a cycle, they started to think that this is an epicycle: the circle moving on the other circle. Soon they realized that even this is on enough and added one more epicycle: cycle moving on the cycle on the cycle. You guess what happened next? They added more epicycles.

Trajectories of planets described using epicycles, as thought before Kepler.

Their number reached 21 when people understood that perhaps something is wrong, and there comes Kepler! He discovered that orbits are not terrible epicycles, but neat ellipses, parabolas and hyperbolas, each of them has a beautiful describing formula.

Kepler orbits of three types: ellipse (red and orange), parabola (green), hyperbola (blue).

So, the shape of the orbits (unless there are no other forces but the gravitation of the attractor) is known analytically. The poliastro project already had two orbit propagators. One of them, Cowell propagator, used iterative solution of the second Newton law. The other, Kepler propagator, utilized the analytical shape of orbits.

The issue was that some highly-eccentric orbits (that look more like a line than something else) propagate incorrectly: the final values of velocity and position do not match what we really want. This issue happened with the Kepler propagator. But at that point I thought that Kepler propagator is also approximate (like Cowell) and so I decided to implement my own propagator called Mean_Motion, that resolved an issue.

You may ask, what was the problem with Kepler? Say, I can understand why iterative approximate propagator fails for difficult orbits, but what’s the problem with the analytical one? The Kepler convergence problem was resolved later in the other pull request :)

What is actually the mean_motion propagator? The second Kepler’s law states that

A line segment joining a planet and the Sun sweeps out equal areas during equal intervals of time.

An illustration to the second Kepler’s law. The sectorial velocity is equal during orbit propagation, so radius-vector swipes equal areas during equal time intervals.

So, there is a quantity (named mean motion) that grows linearly with time. Given the initial coordinates and velocity, we recalculate them to mean motion, shift mean motion to its final value (just adding propagation time times some factor) and recalculate to velocity and coordinates back. This is the mean motion propagator.

Surprisingly, happy end again: the pull request added the new propagator and the difficult orbits started to propagate correctly again. Using the mean motion propagator, in the other pull request I added a possibility to propagate orbit to certain angle (not propagate for some time, but to certain value of true anomaly).

Perturbations

One part of my project goals was to add pre-defined perturbation models as stated in this issue. The idea is that, of course, in space there is not only gravitational force of the attractor present. There is also gravity from some other bodies (which is smaller), solar drag force, atmospheric drag and so on. The idea of a pre-defined perturbation is that we add some extra force to Cowell propagator and integrate both gravitational force and perturbation.

During the community bonding period, I added two pre-defined perturbations. First of them is the J2-perturbation here. You see, the Earth is not a sphere, but some very complicated body. This fact can be accounted by theso-called J2-factor. The J2-perturbations, for example, lead to nodal precession of a satellite, which makes an orbit non-elliptic.

The other perturbation I managed to add was the atmospheric drag force here. Simple it is, the atmosphere drags the body if it is too close to Earth, which is true for some low-orbit satellites.

Orbital decay of the first Chinese prototype space station due to the atmospheric drag

Fixing Kepler propagator

Finally, if one remembers, Kepler propagator failed to analytically propagate some highly-eccentric orbits. After some research I’ve managed to locate the issue: at some point one during calculation one needs to get the eccentric anomaly from mean anomaly, that is, solve an equation M = -F + e sinh F. The equation is transcendental and can not be solved analytically, so the Newton method is applied.

But the Newton method, though very strong, relies heavily on the choice of the initial approximation. And, if the orbit was highly eccentric, the initial approximation F = M did not work anymore: the second term was very large. I changed the initial approximation to F = arcsinh (M / e), and the Newton method started to converge in all cases.

To sum up, the community bonding period was quite a hard work, but there is a lot of new things to do: other perturbations, returning DOPRI835 integrator to the project, benchmarking… See the next update!