8.10 Orbit Integration

We now come to the function that manages the orbit evolution, driving the Hermite integrator and scheduling the various output operations:

Starting again with the argument list, we see that the mass array, as always, is defined as const, since we do not model a mechanism for mass loss for stars, nor do we (yet) allow collisions between stars, which could be followed by mergers that would produce a merger remnant with a mass equal to the sum of the masses of the two stars. The only place where we do not define the mass array as const is in the function get_snapshot, where the mass values are read in from the standard input stream. Note that the time t is passed by reference. In our current program, this is not necessary, since the value of t is not used in main(), where execution is halted immediately upon completion of the call to evolve(). However, in future extensions we may well add further commands in main(), and in that case it would be useful to have the value of the current time available.

As we have seen before, before we can enter the integration loop we have to start with an initial call to the function computing the accelerations and jerks. This function, get_acc_jerk_pot_coll() does what its name suggest: besides calculating accelerations and jerks, it also reports the value of the total potential energy of the system as well as the value of the time scale on which a `collision' between particles can occur, i.e. a significant change of order unity in the local configuration of at least two particles. The latter information, stored in the variable coll_time, will be needed in the main integration loop in order to determine the size of the first time step. Accelerations and jerks are needed for the first part of the first integration time step, and the potential energy is used in the initial call to write_diagnostics(), following the first call to get_acc_jerk_pot_coll().

In addition, if the user has specified the init_out flag to be true, the input values of the $N$-body system are echoed as they are on the output stream; the default behavior is to wait with output until some integration steps have been taken. This is a sensible default, since in many cases we are only interested in one final output snapshot, which can then served as the input for a later invocation of the integrator. If we invoke our program with the same value for the snapshot output interval as the duration of the run, we guarantee that only one final output will be made. An example usage of this type is:

|gravity> nbody_sh1 -d 0.01 -e 2 -o 40 -t 40 < data.in > data.out

Before entering the main integration loop, we schedule the next times for diagnostics and snapshot output, as well as the final halting time. The loop itself is an infinite loop, governed by the tautological while (true), which is obviously always the case. The standard C/C++ trick to define an infinite loop uses an empty for loop, in the form for(;;), but that expression is less transparent, whereas while (true) leaves no doubt as to it being an infinite loop. The only way to jump out of this infinite loop is at the end of the loop: when time progresses past the halting time t_end, the break statement causes control flow to continue past the loop.

The first time we enter the loop, the second while argument will be evaluated as true, unless one of the three values dt_dia, dt_out or dt_tot would be zero or negative, which would be nonsensical values. Ideally, we should check somewhere that all command line option arguments fall within reasonable ranges. Since in the present code we have already introduced so many new features, we will not include such a defensive programming style at this point. However, later on we will insist on checking all values which reach a program through an interface, such as presented by command line options. For now, we will live with the danger of a non-positive value for either dt_dia or dt_out, which combined with a positive value for dt_tot would lead to an infinite number of output attempts, without the time ever advancing.

With natural choices of parameters, the majority of loop cycles will not lead to any output. In those cases a new time step size is determined, and the function evolve_step() is called, which as the name implies will advance the system by one integration step, and in addition update the time by an amount dt. Sooner or later it will be time for output or for ending the run. In either case, the second while statement will evaluate as false, no integration time step will be taken and therefore the time will not be advanced either. Instead, the required output will be done and/or the integration will be finished altogether. If the run is not yet finished, the next cycle in the infinite loop will lead to another integration step, and so on.

Note the freeing up of memory for acceleration and jerk arrays, at the end of evolve(). As in the case of the memory allocation in main(), this is not strictly necessary, since the program is about to finish, but again it is certainly good form to include these statements here.