Numerical integration of stellar orbits

An efficient integrator that uses Gauss-Radau spacings

This describes our integrator RADAU, which has been used by several groups in the U.S.A., in Italy, and in the U.S.S.R. over the past 10 years in the numerical integration of orbits and other problems involving numerical solution of systems of ordinary differential equations. First- and second-order equations are solved directly, including the general second-order case. A self-starting integrator, RADAU proceeds by sequences within which the substeps are taken at Gauss-Radau spacings. This allows rather high orders of accuracy with relatively few function evaluations. After the first sequence the information from previous sequences is used to improve the accuracy. The integrator itself chooses the next sequence size. When a 64-bit double word is available in double precision, a 15th-order version is often appropriate, and the FORTRAN code for this case is included here. RADAU is at least comparable with the best of other integrators in speed and accuracy, and it is often superior, particularly at high accuracies.