Orbital motion of three or more bodies is generally unsolvable by purely arithmetic methods, and therefore requires numerical analysis. That is to say, there's no easy way around it: if you want to predict orbital motion, you'll need an orbital simulation that starts at a known state and continues forward step by step, accumulating small errors along the way.
While this is unfortunate in that it makes modeling orbital motion harder, it also makes such modeling easier to implement: at such a large scale (and with the help of some assumptions), the laws of motion are fairly simple, boiling down to a set of Ordinary Differential Equations (ODEs). In particular, each body's position relative to all other bodies determines its acceleration, its acceleration changes its velocity, and its velocity changes its position, closing the loop.
Solving ODEs using numerical methods works, but introduces some error, since they use discrete timesteps instead of following the smooth curve that perfectly describes this motion. Early versions of my simulation (implemented in Java - download the JAR or source code) simply accepted this error, but the current version goes a step further, working to minimize error through some common numerical methods.
At its most basic, the simulation includes an implementation of a Runge-Kutta solver, which uses sub-sampling within each timestep to get a better idea of how the functions are changing. This is extended into a Runge-Kutta-Fehlberg solver by adding an extra subsampling step to determine roughly how much error is present in the calculation. If this error is deemed to be unacceptably high, the timestep is reduced before performing the calculation again.
While reducing the error is beneficial, it also slows down computation - not only are more timesteps needed, but in practice there also needs to be a maximum timestep to prevent the whole problem falling out of the scope of the solver. This tradeoff must always be considered when implementing such a solver, and can never be 100% perfect.
This project is broken down into two JavaScript scripts: the simulation script, which runs in the background as a web worker, and the display script, which runs in the page and handles all input and output. Each time the user hits play/pause, adds/removes a body, or changes the settings, the display script sends messages to the simulation script asking it to update itself. In return, when the simulation script is running, it keeps a timer going to maintain a constant framerate, and once every frame it sends information on the bodies to the display script.
Of course, sometimes the user asks for too much, and the simulation script can't finish within the allotted time. As a result, the frame is "tardy", resulting in a slow-down of the simulation. This doesn't affect the accuracy of the simulation, but does reduce both the effective framerate and the speed of the simulation.
When the display script receives new information from the simulation script, it updates the positions and sizes of all bodies displayed on the screen. Additionally, if the user changes the display settings, such as the view angle or padding, the display script updates the display without needing new input from the simulation script.
The display itself shows every body in the simulation's 3D space, projected directly onto a 2D plane (the screen). There is no perspective, but there is occlusion, so objects "nearer" to the screen will eclipse those "farther away". Objects displayed as circles are all to scale, but those displayed as small squares appear larger than they really are, so as to make them visible. Additionally, if one mouses over or taps on a body's entry in the display list, it will slowly blink to help highlight its position.
To provide a sense of scale and speed, a grid is placed in the background of the display, and the side length of its squares is shown in the top left. It scales with the display, but attempts to maintain a certain range of pixel sizes, making its squares half or twice as long, as needed. The invariant grid intersection is always at (0,0,0) in the simulation's coordinate system.
The plane of the screen can be rotated along three different axes using the sliders in the control panel, to provide whatever viewing angle is desired. In the back-end, these are simple 3x3 matrix rotational transforms performed on the positions of the bodies when calculating their display positions.
Each time the page is loaded, three large bodies are always added to the simulation in the same places, while four smaller bodies are added with random characteristics. This helps to keep a fairly consistent display of the simulation, while also allowing relatively small differences in each example cause large disturbances, displaying that orbital mechanics are relatively chaotic.
The user can remove any of these bodies at any time by clicking the associated "Remove" button in the display list. This removes the body from the simulation instantly, and puts its values into the fields for adding a body. Bodies can be added to the simulation by filling out these fields and clicking "Add Element", which adds them instantly.
Collisions are modeled as perfectly inelastic in the simulation: the two bodies merge to form one, and their momentum, mass, and volume are conserved. This is achieved by checking to see whether two bodys' Axis-Aligned Bounding Boxes over the course of a timestep overlap - if so, a more rigorous check of their closest approach distance is calculated, and if it is less than the sum of their radii, they are considered to have collided.
Firstly, this simulation relies purely on Newtonian mechanics, and as such deviates from reality in relativistic conditions such as high speeds and large gravity wells. Programming a relativistic simulation using Lorentz Transforms and Einstein Field Equations seems like an interesting project, but would definitely take a lot of time and effort.
Secondly, this simulation doesn't take into account mass distributions. While perfectly spherical bodies are reasonable on some scales, this assumption breaks down quickly in fairly common real-world scenarios: both low Earth orbits and low Lunar orbits suffer from significant perturbation dues to mass concentrations in the orbited bodies. Integrating that type of data into this simulation is definitely possible, but would require some slightly more complex math.
Finally, collisions are all-or-nothing; the simulation doesn't allow for aerobraking or glancing blows. Much of the work required to implement these features already exists, but some new work with regard to structural strength and collision energy would be needed. In the extreme, having a body within the Roche limit break up and spread out would be interesting, if rather computationally expensive.
And of course, if you have any ideas about how this simulation could be improved, feel free to let me know!