5.3 Setting up Three Stars on a Circle

Although our code is written so as to function for an arbitrary number of particles, we initialize only three particles at the beginning of leapfrog2.C. This can be easily changed to a larger number of particles in different configurations. A few chapters later we will see how to read in an arbitrary $N$-body configuration.

It was Bob who was curious to see what happened when we place three stars equally spaced along a circle. After he had coded up the two-body leapfrog, which showed him how stable a two-body orbit of two particles was, he asked Alice whether a three body configuration would be equally stable. She said it wouldn't, and they all agreed that it would be fun to see for themselves. For simplicity, they gave all particles the same mass m = 1.

Anticipating that the instability might take some time to develop, they decided to add one other input option, besides the time step dt, namely the total duration of the run t_end, as we can see in the code. Directly below that request, we see how they spaced the three stars around the unit circle. The first star was chosen to reside at distance unity along the positive $x$ axis, implying a position vector ${\bf r}= \{1,0,0\}$. With the condition that the motion takes place in the $\{x,y\}$ plane, the position of the other two stars is then fixed, as shown.

Following the position assignments, the accelerations are calculated. As we saw before, we have to spin the wheels of the gravitational acceleration module once, before we can safely enter the main integration loop. As far as the code is concerned, there is no surprise here, given that we already figured out the overall structure in the two-body case. The only thing new, apart from the use of two-dimensional arrays, is the double loop over all particles in the system. The outer loop once visits each particle $i$, while the inner loop visits only those particles $j$ for which $j>i$. At the end of the loop the accelerations of particles $i$ and $j$ on each other are calculated. The reason to avoid $j=i$ is that the particles show no self-attraction. The reason to skip $j<i$ is a matter of efficiency: once we have calculated the intermediate quantities needed to calculate a[j][k], the acceleration of particle $j$ caused by the gravitational attraction of particle $i$, we may as well calculate a[j][k] on the spot, which is the acceleration of particle $i$ caused by the gravitational attraction of particle $j$. In our equal-mass case the two are in fact equal in magnitude and opposite in sign, but in the more general case where the two masses are unequal, we would introduce a mass array, and write:

                a[i][k] += m[j] * rji[k] / r3;
                a[j][k] -= m[i] * rji[k] / r3;

The next group of lines employs a shortcut in order to determine the initial velocities of the three particles. Since an $N$-body system is governed by second order differential equations, we have to supply positions and velocities for all particles as initial conditions before we can start to solve the equations of motion, which give us the accelerations. In order to guarantee that our three stars start off on a circular motion, we use the fact that the centrifugal acceleration in a uniform circular orbit is given by $v^2/r$. Since the gravitational centripetal acceleration has to balance the centrifugal acceleration, we can solve for $v$, the magnitude of the velocity, once we know $a$, the magnitude of the gravitational acceleration, a quantity that is equal by symmetry for all three particles, as is $v$ and $r$. We find $v = \sqrt{ar}$.

In our case, the magnitude of the position of each star was chosen to be unity. Therefore, the magnitude of the velocities are $v = \sqrt{a}$. The direction of the velocity of the first particle is chosen to lie along the positive $y$ axis. This choice determines the other two velocities v[1][] and v[2][] by successive rotations over 120 degrees, or $2 \pi / 3$ (by the way, while it is nice to let the computer do all the calculating, it is not difficult to derive with pen and paper that in our case $a = 1\sqrt{3}$ and thus $v = 3^{-1/4}$).

The rest of the code is again a straightforward generalization of the two-body case. Like we did before, we compute the total energy before and after the full orbit integration. Note that we could have folded the calculation of the potential energy in with the calculation of the initial accelerations. Doing so would make the program somewhat shorter, but it would make the logical structure less easy to understand. There would be no gain in overall efficiency either, since almost all the compute happens in the for loop which gets traversed thousands or more times. In comparison, what happens only once before and after the orbit integration can be as inefficient as we want and still not affect the total running time significantly.