For my project i have made a ROOT-program that calculates orbits around a black holes,
or other heavy objects where relativistic effects are important.
The program shows the trajectory in a TCanvas after the calculation is done, and some numbers describing the orbit is printed in the ROOT-console.
The numbers are:
click here if you want to inspect the code just using the browser.
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.
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.
The proper time of the orbit is the period of one orbit as measured by an observer sitting on the orbiting object.
That time is first and foremost something you measure with a watch, and not something universal for all observers is a simple, but important idea in relativity.
The coordinate time of the orbit is the period of one orbit as measured by an observer sitting outside the gravitational field, far away from the black hole.
The general theory of relativity tells us that the proper time and the coordinate time are different. The observer on the orbiting object and the observer outside the gravity-field each measure a different period of the orbit!
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.
To run the program, bring the program into the pwd in ROOT, and type ".x orbits.cpp".
You will then be prompted to input the parameters for the calculation, which are the following:
After you have input the parameters, the calculation starts. 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.
You can't always trust the result of the calculation. The example with the orbit of Mercury demonstrates this.
Further down on the page there is more discussion about when you can trust the results, and why.
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.
Neil Degrasse Tyson may have told you that you shouldn't go too close to a black hole, which is also the case for this program!
If the program gives you nonsense results, like a NaN time-dilation factor and a weird plot, you have probably 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:
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.
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 precession of the apoapsis of Mercury's orbit was one of the first hints that Newtonian gravity was not the final theory of gravity. 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 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, this new theory predicted that relativistic effects also contribute to the precession. When you do the math again using general relativity, all of the precession of Mercurys orbit is accounted for. Nice!
We know the parameters for Mercurys orbit, and the precession is known. Then we can try to calculate the precession with the program to see how it performs.
We get the following plot:
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 where the object needs two full revolutions beetween two apoapses.
This means that the apoapsis and the periapsis happen at the same angle!
I got this plot with the following parameters:
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 that the calculation of dr/dphi makes sense 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 a choice that has to be made as to what kinds of parameters the user should put into 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 per mass 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 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 some other method. Another accuracy problem is finding a more clever way of switching between the object moving outward and inward, and recording the angle of the turning point, so that the error in the precession can be reduced.