8.14 Where All the Work is Done

We now arrive at the core function of nbody_sh1.C, where all the hard work is being done. In addition, this function is both the longest and the most complicated among the ten functions in the file. The main reason for the complexity is that we are trying to accomplish four things in one function, as the name indicates. While calculating accelerations and jerks are logically related, the calculation of the potential energy and the collision time is more a matter of convenience with little natural or logical relation to the calculation of the first two. The main reason for bundling these four operations is efficiency. Here is the code:

Notice the distribution of const declarations here, which is just the opposite from what we saw in predict_step() and correct_step(). In the latter two accelerations and jerk were const while positions and velocities were updated. Here the roles are reversed. In addition, there are two variables that are called by reference, epot and coll_time, which enable the information about potential energy and collision time to flow back to the calling function evolve_step() and from there back to evolve(), where they are used, as we have seen above.

After preparing the proper initial values for the four variables of interest, we enter the $\{i,j\}$ loop running over all particle pairs. As we have seen in the previous two chapters, we first compute a number of auxiliary quantities before we are ready to calculate first the contribution of a pair of particles to the potential energy and then their mutual contributions to each others acceleration and jerk.

At the end of the loop, we compute the two different collision time step estimates, in the same way we discovered at the end of the previous chapter. The first estimate follows the approximate of unperturbed linear motion, extrapolating current separation and rate of change of separation in order to guess when the particles will change their relative configuration substantially. The second estimate neglects the current rate of change of the pairwise separation, estimating instead the free-fall time of the two particles, in case they would start off at rest. In practice, the smaller of the two estimates provides a reasonably safe estimate for the time scale on which significant changes in configuration can occur.