Project: Orbits around Black Holes

Project description

For my project i have made a ROOT-program that calculates trajectories for objects orbiting a black hole, and other heavy objects where the General theory of relativity is necessary for calculation.
The program shows the trajectory in a TCanvas after the calculation is done, and information about the orbit is printed in the ROOT-console.
The information shown is:

These expressions are explained further down on this page.
On this page you can find instructions on how to use the program, and some parameters for known objects and orbits.


click here if you want to inspect the code using the browser.

Explanation of expressions


The apoapsis of the orbit is the radius of the point on the orbit that is furthest from the black hole.


The periapsis of the orbit is the radius of the point on the orbit that is closest to the black hole.

Precession of the orbit

When an orbit precesses, the angle where the orbit hits its highest point changes with each cycle. The precession is quantified by the difference in angle between two successive highpoints.

Proper time

The proper time of the orbit is the time of one orbit as measured by an observer sitting on the orbiting object.
The general theory of relativity tells us that the proper time and the coordinate time are different!

Coordinate time

The coordinate time of the orbit is the time of one orbit as measured by an observer sitting outside the gravitationals field, far away from the black hole.
The general theory of relativity tells us that the proper time and the coordinate time are different!

Schwarzschild Radius and Event Horizon

The schwarzschild radius of an object is given by

The event horizon is the sphere of radius r s around the center of the black hole. Objects that cross the event horizon will never be able to come out again.

Using the program

To run the program, bring the program into the pwd in ROOT, and type ".x orbits.cpp".
Doing this will bring up a message telling you that you can go to this page for more information about the program.
You will then be prompted to input the parameters for the calculation, which are:

All these parameters are to be entered in kilometers.

After you have input the parameters, the calculation starts, and when it is done, the program will open a TCanvas with a plot of the orbit, and some information about the orbit will be printed in the ROOT-console. When looking at the results, you should not trust the calculation of the precession if it is very small. See the part about the orbit of Mercury further down on the page for an example of this.
The plot will also show a circle representing the event horizon of the black hole, which may or may not be visible, depending on the scale of orbit.

You may have heard from Neil Degrasse Tyson that you shouldn't go too close to a black hole, and the same is the case for this program! If the program gives you nonsense results, like a NaN time-dilation factor and a weird plot, when you try to go very close to the event horizon, you probably have tried to plot an orbit that is not physically possible.
For example, if you try to plot an orbit that goes between 2 times and 4 times the Schwarzschild radius, you'll get nonsense.

Some options, like setting the number of steps in the calculation of one rotation, or the number of rotations that are to be calculated must be set by going into the script file. These options are set at the beginning of the code.
The default values are:

If the orbits look strange, for example if the apoapsis of the orbit looks "pointy", you should try to increase the number of steps in the calculation.

Some limitations of the program:

The program is not meant for circular orbits! If you try to plot a completely or very nearly circular orbit, the part of the code keeping track of when the object is moving inward or outward may go nuts. Also, the problem of solving for circular orbits are probably best handled by solving algebraic equations anyway.

When you try to calculate the precession of a very non-relativistic orbit, you cant trust the answer. If the precession is too small, the actual value will drown in the uncertainties from the numerical methods. An example of this is given in the example of the orbit of mercury below.

Parameters some known objects and orbits:

Schwarzschild radii:


Examples of calculations

I will show three examples of calculations done by the program. The first and last examples shown are orbits where relativistic effects are clearly present, especially the precession of the apoapsis.
A precession of the apoapsis means that the angle at which the orbit reaches it's highest point changes with each rotation.
I got the following plot with these parameters:

Here you can very clearly see that the apoapsis of the orbit precesses with each revolution, which is a consequence of general relativity!
The program calculated that the orbit precesses about 155 degrees each revolution.
Another calculation done by the program tells us that the orbital time as measured by an outside observer far away from the gravitational field is about 1.06 times longer than the orbital time as measured by a clock sitting on the object, another relativistic effect!

The orbit of Mercury

The precession of the apoapsis of Mercury's orbit was one of the first hints that Newtonian gravity was wrong. A Newtonian gravitational interaction between two perfectly spherical bodies, with no other forces at work, produces perfectly elliptic orbits, where there is no precession. To calculate the actual orbits of the planets in our solar system, you have to take additional influences into account, like the gravitational interactions between the planets, and the fact that the sun is not perfectly spherical. These other influences creates precessing orbits. When you do the math, you find that not all of the precession of Mercurys orbit can be explained by these non-relativistic effects. When Einstein came along, with his general theory of relativity, the theory predicted that relativistic effects also contribute to the precession, and when you do the math again, using general relativity, all of the precession of Mercurys orbit is accounted for!

We know the parameters for Mercurys orbit, and the precession is known, so we can put it into the program as a test of its performance. We get the following picture:

And a calculated precession of -0.0288 degrees per cycle, using 100'000 steps in the calculation of one rotation.
When we convert this into arcseconds per century, we get 43'200 arcseconds/century. When we compare this to the calculated precession of 43 arcseconds/century [Source], we see that this isn't even close. The uncertainties from the numerical calculation are too big to give a usefull answer for too small precessions.
We can also check the calculation of the orbital time, which is 7.604x10^6 seconds = 0.2411 Years, which fits much better with the actual value of 0.2408 Years[Source].

Finally, here is another cool plot of a closed orbit. The cool thing about this one is that the object needs two full revolutions beetween two apoapses, instead of just one, which is the case for all orbits in Newtonian mechanics. This means that the apoapsis and the periapsis happen at the same angle!

I got this plot with the following parameters:

Computer methods

To solve the equations of motion numerically, i used the second order Taylor-method, using the second order Taylor-approximation for ri+1 around ri.
The theory gives us explicit equations for these two as a function of r, the energy per mc^2, and the angular momentum per mass, aaalmost. The theory only gives us the absolute value for the first derivative. Then the program must keep track of when the sign of dr/dphi is negative and positive, and when the sign changes. That is: when is the object moving inward and outward, and when does it hit a turning point?
I put in two different mechanisms for changing the sign of the derivative when this is needed. One is to check if ri + 1 would lie outside the area where dr/dphi is properly defined with the current sign of the first derivative (that is, will the part of the expression for |dr/dphi| inside the square root still be positive, so the calculation of dr/dphi makes sence in the next calculation). If this is the case, the program switches the sign before it calulates ri+1.
The other mechanism checks if dr/dphi changes sign along the path between ri and ri+1, using the second order taylor-approximations ri. If this is the case, the sign changes before the calculation of ri+2.

When making the user interface, there is some choice in what kind of parameters the user puts in to the program. The mass of the gravitational source is an obvious choice(given by the Schwarzschild Radius), but we also need two other parameters characterizing the orbit.
At first i tried to use the energy and angular momentum per mass, but in the end i decided to use the radii of the apoapsis and periapsis. These seamed much more intuitive and "fun" to work with, as the energy and angular momentum of the orbiting object become more abstract quantities when we are working in the regime of general relativity. This also makes the program much easier to use for people who haven't done any general relativity.


In the examples i discussed the orbit of Mercury, and tried to calculate the precession with the program. The result was not very good, beeing several orders of magnitude bigger than the actual number.
The reason that we didn't get a good result is that the actual value drowned in the uncertainties from the numerical calculation. If we calculate the precession of an orbit with a much larger precession, we should be able to trust the result much more (that is, trust it at all!).
The other parts of the calculation for on-relativistic orbits should be fine, as long as the orbit isn't too elliptical. For the calculations of the very relativistic orbits, i am very happy with the result, again, as long as the orbit isn't too elliptical

Finally, some ideas for extension of the program.
One possibility is making the user able to plot two or more orbits in the same plot, so one can compare them easily.
An obvious extension is improving the accuracy of the calculation, maybe by using a higher order Taylor-method, or other kind of method entirely. Another accuracy problem is finding a more clever way of switching between the object moving outward and inward, and marking the angle of the turning point, that decreases the uncertainty coming from this.