3.2 Writing a Code

After some further discussion among the three, here is the code that Alice typed in. Actually, the first version contained a few errors, not surprisingly, but to speed up a bit, we will just show here the final product, by printing the file forward_euler1.C:

Let us look at each line in turn. The first #include statement is needed in C++ in order to gain access to the I/O library, to make sense of output statements like ``cout <<'' at the end of the code. The second #include statement gives access to the math library that contains functions such as the square root sqrt(). The third line is a default C++ way of indicating that we use the standard name space. Later, when our software environment will have grown significantly, it will make sense to introduce our own name spaces, but for now the default choice suffices.

There is only one function in this file, the main() function which has to be present in every C++ program. It is defined as type int, which means that it returns an integer value (by default defined to be 0 upon successful completion, and nonzero otherwise). The three arrays r[3], v[3], a[3] store the values for the three Cartesian components of the relative position, velocity and acceleration of one particle with respect to the other. These are declared as double, the C++ way to indicate 8-byte floating point variables.

In this simple program there is no need for any input: all initial values are chosen already in the program. The time step dt is fixed to be $0.01$. The relative position is chosen along the $x$ axis at a distance of unity from the origin where the other particle resides in the relative coordinate system that we use. The velocity is chosen to be $0.5$ in the direction of the positive $y$ axis. This is less than the value $v_y = 1$ needed to sustain a circular orbit, as you can check in any book on celestial mechanics or stellar dynamics. The choice of a velocity less than the circular velocity means that we have captured the relative motion at apocenter, which is the point in the orbit at which the two particles are furthest away from each other. The word is derived from the Greek $\alpha \pi o$ (apo) meaning `far (away) from'. The subsequent motion will bring the two particles closer together, until the pericenter is reached on the negative $x$ axis, which is the point in the orbit at which the two particles are closest. This word is derived from the Greek $\pi \epsilon \rho \iota$ (peri) meaning `around' or `near'. After passing through pericenter the particles will move away from each other again.

The for loop takes 1000 time steps that are counted by the variable ns, covering a time interval of 10 units. The body of the loop shows how we first compute the square of the scalar distance between the two particles by summing the squares of the Cartesian components. We then take the square root, and multiply that with the square of the distance to get the third power of the distance, which we have seen in the denominator of Eq. 2.8 on p. , which we repeat here for convenience:

$$\frac{d^2}{dt^2}{\bf r}= - \frac{{\bf r}}{r^3}$$ (3.3)

In our code we find this equation back in component form as:

        a[k] = - r[k] / (r2 * sqrt(r2))

Finally, the integration step is completed by executing Eqs. 3.1, 3.2 on p.  above, to update the relative position and velocity. At the end of each time step, we print the three position and three velocity components, all on a single line.

To compile the code, we must invoke a C++ compiler. In a Linux environment, the natural compiler is the GNU C++ compiler that is called by the g++ command:

g++ -o forward_euler1 forward_euler1.C

which will produce an executable file forward_euler1, if no errors are detected during compilation. Let us follow the discussion of our three friends, from the moment that their code first compiled.