The end of May in GSoC’18 poliastro project
Writing our own ODE integrator
The poliastro’s main feature is the ability to propagate bodies on their orbits during time. However, if there are some forces other than the usual Newtonian gravity, one needs to integrate the equations of motion numerically. The equations that we integrate are called the ODE (ordinary differential equations). The ODE describes how some function evolves with time.
The poliastro, of course, already was able to do this, however, it used the old SciPy API for ODE solving. There were two main problems with this old API: first of all, the core functions were in fortran and, most importantly, the solver returned only the solution at the final timevalue. This was not perfect, cause poliastro not only wants to propagate orbits, but also to plot them. In order to plot an orbit, one needs to know the body position at around 100 points. And, if one wanted to plot an orbit, one had to first evolve from time 0 to time 1, then from time 0 to time 2, from 0 to 3 et cetera. That was of course very inefficient, it would be much better to just propagate once and remember the intermediate body positions.
Recently the SciPy project has released the new ODE solver API. It also provided the dense_output feature that was essentially what was lacking in the old API. But, unfortunately, that was not simply changing the API: in the new API the solver that poliastro loves very much, was not automatically supported.
The ODE x’(t) = f(x, t) can be solved quite easily with the so-called Newton method.
How to get the x(t) function value at some time t + dt? Simply say x(t + dt) = x(t) + dt f(x(t), t). Simply coordinate equals previous coordinate plus velocity times time, that’s essentially it. This method is pretty simple, but it is not very precise: the error of the solution depends very much on dt, and in order to get the decent error, one needs to make dt extremely small, and the integration will take almost infinite time.
Luckily, around 1900, Runge and Kutta developed new integration methods. The dt step could be done much bigger, but the solution was still very precise. Much later, in 1980, Dormand and Prince developed the RK-like method of 8th order. That is, the relative error of x is proportional to dt to the power of 8, which is really small. The problem is that this method was really hard to implement.
With the large help of Juan (my poliastro mentor!) we managed to get the readable version of DOP835 method written in 90-fortran. I’ve spent quite much time making it purely pythonic, but finally it worked. And later we begun the support of dense_output, which made orbit plotting with Cowell much faster.
Last, the best part of it all, is that with the help of Juan it may be possible to add this integrator to the SciPy project after the summer, and this make the world a little better place.
Third body perturbation
Finally, I’ve started working on implementation of the 3rd-body perturbation. The idea is that, say, Moon can distract the usual satellite orbit and for some purposes this should be accounted. What is interesting about this case is that this perturbation is time-dependent, that is, the vector of the perturbation force depends on the position of the Moon. To extract it, we use some astropy features. But about this thing, I will write in the next post, when we are done.