Simulating Star Systems with REBOUND

Simulating Star Systems with REBOUND

For my current project, a novel set in an imaginary star system ('somewhere on a planet far far away'), I need two planets that are near enough to each other once a year (give or take) so that the one planet is clearly visible in the sky of the other (and the other way around of course). One of the planets is inhabited. Let's call that one Kuayot. The other one, Kuauat, is utterly dead. Every year, seen from Kuayot, the planet Kuauat would grow from a bright star to become a visible disk about half the size of our moon. The people would rejoice and festivities would ensue.

This is all nice and dandy, but I do strive for a bit of realism in my writing. I don't just want to make up stuff that's impossible — I write sci-fi not fantasy. So I took to the internet and tried to find out whether such a star system can exist?

My search led to a slew of interesting articles, but an actual and clear-cut answer was not to be found. Clearly, I needed to dig deeper. And what better way to find out whether such a configuration is possible than by observing an actual system with the desirable properties. Sadly, my budget for exo-planet research is a bit limited. All I have is a cheap telescope from ALDI. That's wasn't going to cut it. The next best thing was well within my reach though. So I went online again and found an excellent bit of software that allowed me to model and simulate a star system.

REBOUND, the open source N-body integrator

REBOUND is an N-body integrator. Or, in other words (and I've ripped this right off their site): it's a software library that can integrate the motion of particles under the influence of gravity. It can do much more than simulate planets orbiting a star and it's open source as well. The library has bindings for the programming languages C and Python.

The library itself just operates on numbers — coordinates, orbital parameters, particle masses. Some examples to plot orbits from Python using matplotlib (another nifty little library that can produce plots from data) are provided, but I wanted something more interactive and real-time.

I don't usually create graphical user interfaces, so I needed something simple. My first attempt used graphics.py, which turned out to be a bit too simple. It did allow me to draw orbits, but it's more suited for static images (much like matplotlib actually) and not so much for live action. As the simulation progressed, drawing to the screen became slower and slower (and there's a good technical reason for this).

My next bet was pygame, a library for building games in Python. While it is a bit more complex, it is much more suited for my purposes. With its blitting engine and layered graphics, it let me easily build a live view of the planets orbiting a sun while leaving a trail to draw their paths. Pretty soon, I had a animation of two planets orbiting a sun and I was ready to experiment. To finish things off I added a screenshot function (leveraging another great Python library, opencv2) and expanded that into a video recording function as well.

Playing with planets

All seven planets after running the simulation for a simulated time of about 164 years

All seven planets and their orbits (click for video on youtube)

Once I had the visualisation, it was time to define the system. I started out with a star much like our own and the two earth-like planets — one a bit smaller and the other a bit larger than actual earth. Their orbits are slightly eccentric and their semi-major axes are rotated relative to each other to have their orbits extend out in different directions. To make the system a bit more interesting, I sprinkled two smaller planets closer to the star, two gas giants (or ice giants, I'm not entirely sure yet) in wider orbits and finally a Neptune-like massive rocky planet even farther out.

It took a bit of tweaking to get things exactly how I wanted them, but that was exactly the point. The simulation allowed me to see the effects of my tweaks play out in the simulation. Refining the parameters got me close enough to my desired system that I was confident that the imagined setup is possible.

Just the inner four planets (Kuayot is the pink one, Kuauat the red one)

Inner planets only, Kuayot is pink, Kuauat red (click for video on youtube)

The whole process is fascinating. Seeing the planets orbit their star, their relative speeds, is hypnotic. The difference between the frantic movement of the inner planets compared to the glacial tempo of the outer-most one puts life on earth in an interesting perspective. By the time that far-away cold bit of rock completes a single orbit, more than 160 years have passed on Kuayot.

See it for yourself on youtube, where I posted three videos of the latest experiment:

Other parameters

Of course, this is a huge simplification of an actual star system. For example, I've not bothered to incline the orbital planes at all. All of the planets are restricted to the same two-dimensional plane. I haven't put in any moons, even though there probably would be.

Also, the precision of my computer is limited and finite. In computing, representing irrational numbers comes with limitations. To properly represent an irrational number, an infinite amount of storage is needed. Naturally, computers don't have infinite amounts of storage. As a result of the way we store approximations of irrational numbers in computer systems, there is a trade-off. Either your number is big and not very precise, or your number is small and very precise. You can't have both.

One problem I could not solve in my simulation due to this trade-off was that one of my two planets inevitably falls behind the other. I just can't set the orbital parameters with enough precision that the planets have the exact same orbital period. I would need to tweak either the distance to the sun or the mass wich more precision than is afforded for these rather large numbers.

I could probably compensate by scaling things one way or the other, or applying other more mathematical tricks. I'm happy with the result though and don't need it to be absolutely right.

One final thing to note is that, because the planets are point masses, this simulation does not take into account the tidal forces acting on the planets when they are close enough to be visible in the sky. I did some quick calculations, estimating the required distance of a roughly earth-sized planet to another earth-sized planet to appear as a disc about half the size of our moon, then ran some quick numbers to get a feel for the magnitude of the tidal forces. It's enough to pull on the oceans, but far from too much. The planets won't tear each other apart.

Try it yourself

If you are confident with Python and want to try this for yourself, you can find the code available on my gitlab instance. It's a bit rough around the edges though. Remember, this was basically just a quick experiment to validate one of many details that go into building the world for my upcoming novel. Enjoy!

follow