2017ApJS..230...19B

We construct an advanced model for interacting multiple stellar systems in which we compute all trajectories with a numerical N-body integrator, namely the Bulirsch-Stoer from the SWIFT package. We can then derive various observables: astrometric positions, radial velocities, minima timings (TTVs), eclipse durations, interferometric visibilities, closure phases, synthetic spectra, spectral energy distribution, and even complete light curves. We use a modified version of the Wilson-Devinney code for the latter, in which the instantaneous true phase and inclination of the eclipsing binary are governed by the N-body integration. If all of these types of observations are at one's disposal, a joint χ² metric and an optimization algorithm (a simplex or simulated annealing) allow one to search for a global minimum and construct very robust models of stellar systems. At the same time, our N-body model is free from artifacts that may arise if mutual gravitational interactions among all components are not self-consistently accounted for. Finally, we present a number of examples showing dynamical effects that can be studied with our code and we discuss how systematic errors may affect the results (and how to prevent this from happening).