This code reads no input, even from the command line. It is meant for those who can change it in ways that I do not anticipate. I describe here some of its logic and switches which may be useful. Words below in bold occur in the source code and are used here to point into the code.

Routine qt is how the code quits in an organized way. The first argument is printed to identify the call site.

The compile time parameter D is the number of dimensions the code processes. N is the number of masses. It has been tested for N=2 and N=3 while D=2. Some testing has been done for D=3. N must be adjusted to match the definition of z and initialization code checks for this match.

The defined type R should probably be defined as either double or long double. All of the dynamic variables are of this type. At least on the x86 architecture for Linux and Mac OS X, long double provides 64 significant bits. It is necessary in this case to enable the definition of sqrtl which is otherwise suppressed.

The structure Dv (struct{R x, xd;}) holds a dynamic variable and its first time derivative. pt (struct {Dv dv[D]; R m;}) carries all the information for one mass. ss (struct {pt pos[N]; R t, dt;}) is a description of the entire system state.

Routine spin has all of the physics and performs an entire run of the physics. 

void spin(pt z[N], R tlim, char draw, int phop)
spin may be called twice to determine sensitivity to initial conditions. Parameter z points to a block that is the initial system state and and spin mutates it to track the running state of the dynamic system. The tlim argument tells spin how long to run. spin’s argument draw is 0 or 1 and controls whether the spin creates files with data to be rendered as .png files by the Java program. The argument phop is the frequency, in cycles, that the numerical system state is sent to the window. Here are short descriptions of several invocations of spin.

Routine ri is called at the beginning and the end of a run. It computes and prints the known invariants:

Routines e, and E compute the potential energy between two masses, and the kinetic energy of one respectively. They are used only for reports. The routine pr reports the system state and is typically called periodically.

step is called once each time step to advance t and the dynamic variables in z. Step reads Δt and updates z[n].dv[d], the dynamic variables. It also adjusts Δt down if any mass pair are too close and adjusts it up if all pairs are far enough apart. f is at the core of the Runge Kutta method which proposes a way to solve the differential equation x' = f(t, x). f(t, x) computes the forces for change given the state x. The third argument to f, dtc, merely tells f whether to pay the cost of worrying about Δt control. kt (struct {Dv dv[N][D];}), defined within step, is a convenient package which describes an argument and the returned value of routine f(x, t). lc is an auxiliary routine to produce the linear combination of states, required by the method.

ML is a block from which the entire execution is controlled. A px (struct{short x, y;}) remembers where a particular mass was depicted so as not to record it there again in the plot output file. The array v has a px per mass. v is the memory for routine cx which reports changes.

The following variables are involved in the same pattern to provide periodic functions thruout the calculation:

The variables are each initialized in block ML. pc is a scaling function to pixel coordinates for the problems at hand.

Loop is the compute loop proper. ph emits 16 bit fixed binary coordinates into the plot file. The left bits go first for I plot on the Mac. You will probably need to reverse the two statements in ph if you run Java on a little endian machine.

If gurg is not zero we compute a state for each integral value of t. Variable Δt makes this non trivial for the natural logic skips over most integral t values. The expression cpnt - t computes a negative value for Δt with which we back up to the integral value. The running system state in z is copied into the temporary tmp during this detour.

The main routine has master switch MS which selects one of two modes:

For the sensitivity analysis we report for each integral time value, the distance in phase space between points on the two compared trajectories.