This Algol68 version 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.

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

The constant d is the number of dimensions the code processes and n is the number of masses. It has been tested for 2 and 3. n must be adjusted to match the definition of z and compiler logic checks for this agreement.

The defined type FL should probably be defined as REAL, LONG REAL or LONG LONG REAL. All of the physics variables are of this type. In the a68g compiler LONG LONG REAL is of variable precision. It is necessary in this case to add as many “long”s before “sqrt”, “sin” and “cos” as there are LONGs in FL

The structure MODE DV = STRUCT(FL x, xd) holds a dynamic variable and its first time derivative. MODE PT = STRUCT([d]DV dv, FL m) carries all the information for one mass. MODE SS = STRUCT([n]PT pos, FL 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. It may be called twice to determine sensitivity to initial conditions. Its first argument, gurg, optionally locates an array of instances of ss into which it will deposit periodic system snapshots for subsequent comparisons. del is a small modification to the initial value of some dynamic variable. Modifying the definition of z is how to change which. This logic has not been tested in the Algol 68 version. spin’s argument draw is 0 or 1 and controls whether the spin creates files with data to be rendered on the screen by the Java program.

z tracks the running state of the dynamic system. Its definition provides the initial values. Three alternate definitions are provided for z with known answers. This can be used to debug modifications to the code. The circle definition is for two orbiting masses and parabolas yields two masses that escape to infinity with 0 total energy.

Routine ri computes and reports the known invariants:

Routines e, and ke 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.