I want to document this simple code before adding any of several pieces of complexity such as:
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.)