Co-orbitals are top-of-mind right now! As-is they can be simulated without issue (the system is defined as a DAG where each body defines the parent bodies it is gravitationally influenced by, e.g. Sun+Jupiter+Saturn for the trojans) but I haven't yet figured out how to get the proper starting position+velocity for them in the current epoch. I foresee a deep dive down the ephemeris rabbit hole in my future...
Initial conditions are modeled by Keplerian elements around the main parent body. These are 6 scalars that describe the size, shape, and orientation of the orbital ellipse as well as the body's position along it at a given time (epoch). I scraped these values from various places including research papers, JPL databases, Wikipedia pages, and university web pages.
From there everything is mapped into a consistent epoch (now) and the initial position/velocity are calculated using a Keplerian->Cartesian transformation (the math here is a more or less straightforward coordinate transformation). The physical simulation is then run by updating these 2 Cartesian vectors, applying gravitational acceleration over ∆t.
This works pretty well, all things considered, but there's a lot of room for improvement to account for solar wind, relativistic effects, planets not being perfect spheres, etc. The Keplerian elements are also just an approximation of the true orbit, and this approximation can show up at high zoom levels (it's why planets don't always line up perfectly with their ellipses when you zoom in). I'm also still figuring out how best to get the initial position/velocity for objects that aren't in simple elliptical orbits (co-orbitals like the Trojans, objects on escape trajectories like Voyager probes). There's a lot for me to learn, maybe later I will write a blog post!
Initial positions and velocities are interesting. How did you get about that info?
I would love to read about how the orbital trajectories are calculated, and how is it done for co orbitals. Do you plan to write a blog on this?