7.3 A More Modular Approach: Functions

Now that we have a general-purpose $N$-body integrator, it is high time to rewrite our code in more modular fashion. The first step is to offload much of the work done in main() into a few functions. Obvious candidates are the loop where we calculate the accelerations and jerks of all particles, and the loop where we determine the total energy of the system. Both of these occur twice in hermite3.C, so by bundling each of these loops in a separate function we can shorten the length of the program. More importantly, using functions makes the program easier to understand, and therefore also easier to extend.

Here is the first part of our new version of the Hermite code, hermite4.C, which contains the first function that calculates the accelerations and jerks:

Note that we can pass the arrays for masses, positions, etc., without having to specify the total length of these arrays. This is a good thing, since at compile time we do not yet know the $N$ value(s) with which the program will be run. In a two-dimensional array, however, we do have to specify the length of the second dimension, in order to allow the compiler to make the correct pointer calculations to map the memory. Fortunately, we will almost always work in three dimensions, so that number can be specified in, e.g. r[][3].

The type of the function acc_and_jerk is void, which means that the function does not provide a return value. All the work done is stored in the arrays that contain the accelerations and jerks. Since in C++ array variables are effectively pointers, any change made locally in the function will directly affect the arrays with which the function is called further on in the program where we will see:

    acc_and_jerk(m, r, v, a, jk, n);

The first five arguments will be directly visible in the body of the function, and all changes made there will affect the five variables used to call the function. The sixth and last variable n is passed by value. This means that changes in n made in acc_and_jerk() will not affect the original value of n in the function main() that calls acc_and_jerk(). In our particular case, there is no need to change n. If there would be, we could have passed n by reference, rather than by value. We will see an instance of passing by reference in the next section.

Here is the second function which calculates the total energy of the system:

The return type is double, which means that calling the function gives an immediate handle on the returned value, the total energy. When we call this function with

    double e_in = energy(m, r, v);

the result of the calculation is directly assigned to the variable e_in.

Why do we add so many arguments to our functions, six and four, respectively? An alternative would have been to declare the five main data structures, from masses to jerks, as global variables at the top of our program, as follows.

The use of global variables is considered bad form, for good reasons. The whole point of our exercise of splitting our program into functions is to isolate functionality. This will make it easier to understand and debug the program, and to modify or extend it later. For now the program is small enough that it is easy to keep track of the global variables. When we add more and more features, chances are that we lose track of exactly which global variables are floating around in the program. Also, it will be easy to get name conflicts, since individual functions are likely to use similar names. Finally, there is a very practical reason to stay with local variables only: some time soon we may well want to manipulate more than one $N$-body system, to compare and analyse them. At that time it will be impossible to use a single global set of variables to store the data. In anticipation of such complications, it is more prudent to stick to the use of local variables right from the beginning.

Here is the remainder of the program, containing the main() function of hermite4.C, which is now shortened to half the length it was in hermite3.C:

Exercise 7.1 (Pythagorean problem: constant time steps)
The Pythagorean problem was introduced in the late sixties, as a severe test for the $N$-body codes and the computer hardware available at that time, The name stems from the way the initial conditions are set up for the three particles involved: the particles are initially positioned in a right triangle, with masses and positions chosen such that the center of mass lies in the origin of the coordinate system. The velocities are all chosen to be zero. The code mk_pyth.C below generates the initial conditions.

Experiment with the hermite4.C code, to see how small the step size has to be, in order to get reasonable results. For example, starting with dt=0.1, while seemingly a small number, is clearly too large, given the large errors reported, below. How small is small enough, for dt, if we want to follow the system for 10 time units?

|gravity> mk_pyth | hermite4
mk_pyth | hermite4 0.1 10
Initial total energy E_in = -12.8167
3
0.1
3 0.998076 2.99118 0 -0.0192721 -0.088297 0 
4 -1.98742 -0.998076 0 0.126071 0.0192722 0 
5 0.991091 -0.996247 0 -0.0892936 0.0375604 0 
3
0.2
3 0.995662 2.98013 0 -0.0290448 -0.132818 0 
4 -1.97162 -0.995662 0 0.190175 0.029046 0 
5 0.979897 -0.991547 0 -0.134713 0.0564539 0 
3
 . . . .
3
10
3 -1.74295 -7.58878 0 -0.307547 -1.1837 0 
4 -353.379 265.193 0 -42.5655 31.6683 0 
5 283.749 -207.601 0 34.2369 -24.6244 0 
Final total energy E_out = 10077.9
absolute energy error: E_out - E_in = 10090.7
relative energy error: (E_out - E_in) / E_in = -787.31
|gravity>

Here is the code listing.