Orbits of One Body around Another

Joachim Köppen Kiel/Strasbourg/Illkirch April 2001

Contents

About the Applet and what it does

The Applet simulates the movement of a body around another under the influence of an attractive force which each body exerts onto the other. This force could be the normal gravity which keeps a star orbiting another, a planet around the sun, a moon around a planet or a satellite around a planet. So you can see what types of orbits are possible in the Universe, and how they depend on the velocity the planet had at the start of the simulation. This force could also be an elastic string connecting a ball to a fixed central post, and you'll be able see what kind of path the ball takes when you release it from its initial position with some velocity. The force can also model the reflection of a ball by a circular wall surrounding the centre. And finally, this force can be described by a law that does not really exist in nature, but is fun to see what would happen to our planets and stars if the World would have a gravity law different from what it really is.

The smaller body's movements are shown by a red path which it leaves behind. You can change the law for the force, the velocity at the beginning of the simulation - the body will always start at the same position. The simulation is done by computing how much the body moves during a short intervall in time. If you make this time step smaller, the calculations will be more accurate, but also slower! You have the choice between several numerical methods of performing the calculations, and between a number of different representations of the situation. Once you've made your choices, you may start the simulation, halt it at any time, and then continue with it.

When you start up the applet, I recommend that you hit the Clear button, and you see something like the following screen.

On the left hand side there is the large white screen showing the central body as the small circle with a cross in the centre and the smaller, moving body as a small triangle. I got this view by letting the simulation pass a few very small time steps, so you'll see essentially the initial situation with the moving body in its start position, and just beginning to move. At the upper right hand corner the current time is shown. The time unit is not seconds or years, but chosen so that one circular orbit is done in one time units. As you'll see below, it is easy to scale results to any problem from atoms to galaxies.

The next control field tells us that the initial velocity had been chosen to be 0.7, that is 70 percent of the velocity necessary for a circular orbit. Then comes the time step, also in units of one circular orbit. With the one shown here, it takes about 100 steps to complete one circle.

The next control is the numerical method; Here we have chosen the simple and quick Leap-Frog. And we had chosen the simple plain view of the neighbourhood around the central body.

Now let's get something done: let's do a circular orbit! Change the initial velocity to 1.0 and hit the Start button. This is what you get when Stopping after about 27 orbits

Since this looks quite nice, such as a circle should look like, let's continue by hitting the Carry on button. After many an orbit, we get this:

There is something weird happening: the circle does not remain a circle. In fact, if you were to measure carefully, the orbit is not a real circle, but kind of an ellipse which wobbles about. What is happening?

Let us investigate this problem a bit more closely by doing something more extreme: launch the satellite by a velocity smaller than circular velocity, say 70 percent. So we should expect a nice ellipse. See what we get:

It is sort of an ellipse, but one whose large axis is slowly rotated abot the central body. Such a movement is called precession. Isn't the precession of Mercury's orbit about the Sun one of the proofs for Einstein's Theory of Relativity? Before we can claim that our applet has given another explanation and shattered this proof, let's check better whether this effect could have other reasons. Increase the accuracy by taking a smaller time step

and we see now a nice ellipse which keeps its orientation for a long while. This shows that the precession was simply due to a poor accuracy of the calculations! (But see here that all is not yet well)

Now you are ready to use the applet for your own experiments. But better mind the accuracy, if you get results that look weird or are unexpected from the textbooks!

The controls:
Force exponent
is the power with which the force increases with distance from the central body. Normal gravitation is -2, and +1 is the harmonic oscillator, i.e. as if the satellite was tied with an elastic spring to the centre
Initial Y-velocity
is the speed with which the satellite is put into orbit (positive velocities point upwards), measured in units of the circular velocity. The x-component of the velocity is zero
Time step
is the time intervall between two successive calculations of the positions. In the units used here, one complete circular orbit takes time units. Influences the accuracy
Leap Frog ...
selects the numerical method
Plot ...
selects what is plotted: normal x-y-views of the satellite as seen from above at three different scales, a Poincaré plot of x and vx when the satellite crosses the x-plane towards positive y, and the Potential as a function of distance from the centre
Start
starts a new model, with the shown parameter values
Stop
halts the simulation
Carry On
continues with the current model. If the time step width was changed, its new value is taken
Clear
wipes the screen (also during a simulation)
normalPlot / quickPlot
this button switches between the normal way of plotting and the quicker method where the plot will be wiped out when it had been overlaid by another window
time
when one stops a simulation, the elapsed time - in our units - is displayed in the upper right hand corner of the plot

If one follows the orbit for a long time, the screen gets filled with overlapping tracks, and eventually one cannot see anything at all. There is a neat way to study the long-term behaviour by making a short-hand description of the perodic movement of the satellite: imagine that we set up a screen perpendicular to the orbital plane, passing through the central body and our starting point. Every time the satellite crosses this screen, we make a record of the position and the velocity. In the applet we show the x-position and the x-component of the velocity. The small square indicates the initial situation. Since circular and elliptical orbits (as in normal gravity field) always return to the initial position, their Poincaré-plot consists of but a single point, in the little square. Hence, this type of plot shows in a very compact way, whether orbits are closed or not. So the information from a great number of orbits can be compressed in a single geometric figure; a nice plot may take a while, of course.

One can prove that for the movement of a single body around another one, the Poincaré-plots must always consist of simple curves and points that lie on simple curves. On the other hand, chaotic orbits would show up by filling an area or may be the whole plot with dots.

Let us look first at the surface-of-section for the elliptical orbit we computed above , using a time step of 0.1 units. This is how it looks after some time, still being built up:

One notices that the points are organized in some groups quite regularly arranged, and as time continues, filling a closed band:

From what was said earlier, the fact that the structure is not a curve but an area, would indicate that we deal with a chaotic orbit?! But let's check the accuracy, and repeat the experiment with a smaller time step . Now all points lie on a single curve which starts from the initial situation (the little square)

and after long time finally make a complete closed curve;

What is the explanation of all this? First of all, the surface-of-section is also affected by the numerical accuracy. This is why in the first experiment we got something that looked only like chaos, but never was! Secondly, the form of a closed curve in the surface-of-section is characteristic of a precession of the orbit, with the curve showing how the intersection of the elliptic orbit with the screen passes through all distances from the centre between apocentre (largest distance) and pericentre (smallest distance). The curve found in the more accurate simulation is the same as the middle of the band from the first one. This means that both simulations lead to the same precession, independent of the time step. This also means that our first impression from the direct picture of the elliptical orbit that the smaller time step of 0.001 units had cured the precession, was simply wrong! The precession only slowed down, but remained the same.

What is the reason for that? Since we know from Kepler's laws anmd it can be proven from first principles that the orbit of a single planet around a perfectly spherical sun cannot precess, it can only be something to do with the details of the calculations......

The computation consists of calculating how much the position of the body changes during a small time step, due to the force that is exerted on it from the central body. This is done in this sequence

• from the present position compute the x- and y-component of the force (and accerelation) which depends on the distance r = sqrt(x^2 + y^2) with the force exponent n

a_x = f_x(r) = -x r^(n-1)

a_y = f_y(r) = -y r^(n-1)

The minus-sign of the acceleration shows that the force is attractive, pulling the bodies together.

• assuming that the force does not change during this small time step dt, we can compute how the components of the velocity have changed:

v_x(t+dt) = v_x(t) + a_x dt

v_y(t+dt) = v_y(t) + a_y dt

• then, by assuming that the "new" velocity is a good approximation for the velocity during this intervall, we compute how the coordinates of the body have changed:

x(t+dt) = x(t) + v_x(t+dt) dt

y(t+dt) = y(t) + v_y(t+dt) dt

To do our simulation, we simply apply this sequence many times over and so follow the movement of the body.

The method descibed above is what is called the Symplectic Euler method. It is a modification of the simple Euler method, because we use the new velocities to recompute the positions. This trick is very useful, as one may show that the angular momentum is conserved during every single time step with machine precision. Since the law of conservation of angular momentum is one of the fundamental laws of nature and strictly obeyed by any moving thing, application of this method makes the simulation a much better approximation of nature than if there remained some numerical errors. For anyone who likes to work this out, you can show with the above formulae that

x(t+dt)v_y(t+dt)-y(t+dt) v_x(t+dt) - x(t)v_y(t)-y(t) v_x(t) = 0

and that with the straight Euler's method (x(t+dt) = x(t) + v_x(t) dt) does not do so neatly!

These methods are available:
Explicit/implicit or Symplectic Euler
is explained above
Leap Frog
is a second-order method. Positions and velocities are computed on times staggered by a half step. With the velocity one half step earlier, the new position is computed. With the accelerations at that position one computes the velocity one half-step later ...
Trapez predictor-corrector
this is a simple method of the perdictor-corrector variety. First, one takes the information from the last and present point to estimate what the next position would be. Then, this guess is recalculated using the information of the new point as well. In this way, one takes into account that during one step the forces have of course changed
Runge-Kutta 2
The Runge-Kutta method of second order makes one intermediate step and therefore also makes a better estimate of the changing forces during a time step
Runge-Kutta 4
The classical fourth-order RK method takes 3 intermediate points and is more accurate
Symplectic 2
is similar to 2nd order RK, but by a slightly different formulation it conserves angular momentum. It is more accurate than the Leap Frog method, although it needs for the accelerations only one calculation per time step!
Symplectic 4
is fourth-order Runge-Kutta, but with conservation of angular momentum
Milne 4 PC
one of the classical fourth-order predictor corrector schemes (see Maths books for details, e.g. Abramowitz and Stegun)
Milne 6 PC
sixth-order predictor-corrector scheme with Milne's formulae
as above, but of fourth-order

Now let's have a look how some of these methods work with our elliptic orbit and its problems. The classical Runge-Kutta 4 method with a fairly large time step is rather useless:

The orbit now shrinks into a circle which eventually will collapse into the central body. The reason is simply that with this method there is a tendency to loose energy. And if the time step is too large, this becomes really catastrophic! But taking a smaller time step, the RK4 method works rather nicely. It is better than the Leap Frog at the same time step.

But how it behaves in the long run, I leave to you to explore with the Surface-of-Section plots. All I'll say here that they are different from the ones we got with the Leap Frog.

Also, I leave it to your exploratory spirit to have a look at the other methods. Finally, I want to show this very boring plot:

This was done with the fourth order symplectic method. Even after a long time the orbit still remains an elliptical one, and it precesses only very very very slowly.

How the attracting force between the two bodies depends on the distance between them, is described by the force law. We consider a simple law where the force depends on the n th power of the distance. This force exponent we can change from its normal value of -2 for normal gravity. The Plot ... choice lets you look at the potential curves for any of these laws.

Changing the exponent is not entirely a game, as the movement of a star under the influence of a large number of other stars in a galaxy can well be described by an average force law with an exponent different from -2; and different from a simple power law as well. So what are the characteristics of the motion of such a star within a galaxy, is one of the important questions.

Here, we'll just play a bit with Nature's law of gravity to find out what would be if the World was different! What kind of satellites we could launch if the gravity exponent were -2.5, so the force depends a bit more strongly on distance:

If you try a circular orbit, you'll a circular one. But what had been an elliptical orbit, now becomes a spiral! The satellite spirals down towards the centre, and up again, and reaches another apocentre on the other side. What interesting possibilities for Earth observation satellites!

The surface-of section shows that the orbit evidently precesses. We also see from the scattered points in the upper right corner that the accuracy of the calculation is insufficient.

If we use a smaller time step, we see that the surface-of-section starts out from three points corresponding to the three intersections of the spiralling orbit. See what kind of orbit you get with a better method and a smaller time step. You will notice that the taming of the precession takes an effort, and one suspects that it is purely numerical, so that the orbit with its interesting shape does not move about.

Becoming more daring, try an exponent of -2.75! Here the spiraling takes the satellite even closer to the planet, and on a more complicated path. The Poincaré plot shows again the signature of a precessing orbit. The precessing has now become so strong that it causes the orbit to look more like a spiral. Note that the speed near pericentre becomes too large for our plot:

Here is a picture done at very small time steps to show the details of the spiral path: from apocentre to apocentre, we have circle the planet three times - in normal gravity we do that only once.

Even more daring is to use an exponent of -3: You will find that you have a hard time keeping a satellite on a circular orbit. Making the time step smaller and smaller will still lengthen the time before it either plunges into the planet or escapes from it. Here the satellite crashes into the planet, and because of numerical difficulties which occur when it comes too close to the centre itself, the satellite gets a (numerical) kick and dashes off:

Mathematical analysis shows for such a force law circular orbits can exist, but they are instable. The smallest perturbance, like throwing something out of the satellite, will either send it towards the planet or away. Elliptical orbits don't exist any more, they all crash into the planet, like this one here, no matter how small you make the time step. In such a Universe, space travel would be simply catastrophic.

There are some general features of orbits around a central body: First, there is always a characteristic speed that inserts the satellite into a circular orbit. If the initial speed is larger than this circular velocity, one gets orbits that take us further away from the planet, and we might even escape completely from its gravitational pull.

If the insertion velocity is smaller than the circular velocity, we get in normal gravity eliptical orbits that take us closer to the planet, in general we get rosette-type orbits. The motion of the satellite is governed by two effects: it moves around the planet, but it also moves from apocentre to pericentre and back to apocentre. The time period for one circumnavigation and the one between two apocentres do not need to be the same. It is only in the inverse square-law of normal gravity (force exponent = -2) and with the force exponent +1 which corresponds to the satellite tied to the centre by an elastic string, that the two period are equal. Only in these two cases, one gets elliptic orbits and no orbit will ever precess. With our normal gravity, we are lucky to live in a special Universe!

For all other force exponents, the two periods are not equal, and one gets rosette orbits. They can either be closed, when the ratio of the period is a rational number, like 2 or 3/4 or 551/8. Then the satellite returns to its initial position and initial velocity after no-many passes around the planet. The orbit is closed. If this ratio is an irrational number, which cannot be expressed by a simple ratio of integers (examples are 1.4142...., 0.3010299..., 3.1415...), the orbit will not be closed. The satellite will never come back to exactly the same position, though come as close as you want - if you wait long enough. See what happens to our elliptical orbit near the normal gravity:

With an exponent of -1.9 the orbit is a bit wider than an ellipse, so one precesses opposite to the sense of circling the planet, but ...

... with -2.1 the orbit is a bit more narrow, and thus the precession goes in the same sense.

Here is an example for a closed orbit which takes just a few passes around the planet to come back. The Poincaré plot shows that there are only four intersection points. There is still some drift of these points, so either the numerics is not yet accurate enough, or the parameters - force exponent or initial velocity - are not yet the exact values for this rosette orbit

The simulations are done in a normalized way, we use as a unit of length the distance of the initial position of the satellite, and as a unit of velocity the speed for the circular orbit with that radius. All times are given in units of the period of the circular orbit.

The physical problem of the movement of a single satellite around a central body permits to scale these results to any true situation: If we compute the orbit of a planet or asteroid around the sun, with an initial distance of say 1 astronomical unit (orbit of the Earth) we compute from the gravitational constant and the mass of the sun that the circular speed is about 30 km/s and the period for one orbit is 1 year. Our elliptical orbit model with initial speed of 70 percent thus means a true velocity of 21 km/s, and when you'll measure the time for one orbit you'll get about 0.54 time units, which means that the true period for the elliptical orbit will be 197 days.

For a satellite around the Earth we only need to keep in mind that usually the altitude above ground is given, so we need to add the Earth's radius of 6371 km to get the distance of the satellite from the Earth centre. The circular velocity which is about 8 km/s for a satellite at zero altitude, decreases like the inverse square root of the distance from the Earth centre.

In general, all you need to know is the radius, the velocity, and the period for a circular orbit of the initial distance. This you have to compute from the physics involved (gravity, electrostatic forces, etc) and the property of the central body (mass, electric charge, ...).

| Top of the Page | Controls | Applet | Applet Index |