This is code (ISO C) for the Pythagorean three body problem. It is first order and variable Δt but with care taken to preserve 2nd order properties that this sort of differential equation enjoys. The program writes a file for plotting. I change the input parameters by modifying the source. There is little excess hair. I have used gcc extensions to C that are unavailable in some other compilers.

I want to document this simple code before adding any of several pieces of complexity such as:

The struct “pt” holds the parameters and state of one mass. The array called “z” has one pt per mass. The compile parameter N must be the count of entries in the array z. Initialization logic will complain if this is not so. Within a pt are: x & y, the location of the mass; xd & yd, the velocity times Δt; m, the mass; and mid which is Δt2/m which serves to speed the inner loop. The factor of Δt in the velocity is introduced during initiation and mid is also computed then. The first 5 fields are compiled in initial values. There is no code to modify m.

The initial value for Δt is nominal for the automatic Δt adjustment soon makes it what it wants. The two lines below the definition of array z are for alternate orbit problems with known solutions. These can be used to test modifications to the logic. The circle test case is for two equal masses orbiting each other. The parabola is also for two equal masses with the left mass describing the parabola x = y2.

For any of several reasons the program stops politely by calling routine qt. The argument is the reason code and is designed for someone who has and can read the source. qt reports the current configuration and then quits.

Routines e reports the potential energy between two masses and e reports the kinetic energy for one mass. These routines are used only when printing the system state. Failure to conserve total energy is helpful for indicating numerical problems.

pr prints the state. The argument is merely a comment.

Routine conf calculates and applies the force between two masses. It is called each step for each distinct pair of masses. fx & fy are the components of the vector force between the masses. They contribute to the velocities of the two masses. b is the square of the distance these two masses move relative to each other in this time step. It is computed to control Δt. Δt is adjusted so that there are at least n steps per encounter. The comparisons “if(s2 < n2*b)” test to see if Δt should be adjusted. The call adt(x) multiples Δt and its dependents, by x while fussing over centering. Δt is halved when one of the mass pairs is too close and it is doubled when all are sufficiently separated.

The last five lines of code merely control when to print, when to plot and when to quit. The variable “ol” is a limit on the size of the output plot file.

This version has been adapted for the Pentium’s 80 bit floating format which provides 64 significant bits. “ex” is a type that is normally “double” but is here “long double” for the Pentium’s extended precision. (Also long double arguments to printf are cast to double and sqrtl is substituted for sqrt for the Pentium.)