Generalized physics simulator for asteroids and comets

[cosinekitty]

It would be nice to adapt the existing gravity integrator used for calculating Pluto's orbit to allow calculating the trajectory of asteroids, comets, coasting spacecraft, or other small bodies traveling ballistically through the Solar System. This would be focused on interplanetary movement, not satellites of any planet.

The caller will provide the state vector (position and velocity) of a small body. The StateVector objects also include a time value, which will be needed to calculate the state of the Sun and planets.

I need to think more about the details. Here are some ideas to ponder:

• The caller could pass in a time-step value Δt, and the output would be a state vector at the updated time t+Δt. The caller would be responsible for guessing a reasonable Δt value that does not lead to too much inaccuracy. Or Astronomy Engine could internally break up excessive time steps into smaller time steps as needed.
• The caller could ask Astronomy Engine to choose a suitable Δt, in which case the returned object would be an interpolation model that is valid over the Δt it chooses. (More complicated, but much more efficient for some applications.)
• Maybe there could be one function for calculating the state at t+Δt, and a separate function for interpolating state vectors between pairs of adjacent state vectors calculated by the first function.
• All calculations should be done using J2000 orientation with the SSB (Solar System Barycenter) as the origin. It would be nice to allow coordinate transforms by passing in some enum value, but this would complicate the interface and documentation.
• Should I calculate only the Sun and major outer planets (Jupiter..Neptune)? This is fine for Pluto, but probably not for the asteroid belt or comets that reach the inner solar system. My hunch is to keep the interface simpler by assuming I need to calculate all planets. However, I would calculate only the Earth's position, not the Moon's, but would include the Moon's mass with the Earth's as an approximate EMB. The idea is that asteroids and comets do not come close enough to the Earth/Moon system to justify the extreme cost of calculating the Moon's state vector.
• For unit testing, I would use JPL Horizons data for a few asteroids and/or comets. It would be nice to find a variety of orbital eccentricities. A demo program could provide ephemerides including RA, DEC and visual magnitude.

[pakastin]

Wow, cool! 🤩👌

[cosinekitty]

Here are some more thoughts on this. I'm mostly documenting for my own use. Also, I am thinking in terms of the C implementation. I usually start with C because it requires the most consideration of low-level concerns. Then the other languages tend to follow naturally.

Introduction

Unlike the existing C functions, this one requires memory allocation and deallocation, because (a) the structures are getting larger, and (b) I want to make an opaque data type to decouple outside callers with Astronomy Engine internals. I'm not trying to keep secrets; this is an open source project after all. The purpose is to leave the door open for changing my mind later on how to implement the internals, without breaking dependent code.

Implementation Plan

The following functions will be added:

• GravSimInit(&sim, state) initializes the gravity simulation stepper object sim. The state is an initial state vector for a small body that contains the body's position vector, velocity vector, and the time those vectors are valid.
• GravSimStep(&sim, dt) numerically integrates the gravitational forces of the Sun and planets on the small body over a small time step dt. The function returns an updated state vector of the body at the updated time. Astronomy Engine will not attempt to figure out how large or small a time step should be. It is up to the developer to experiment and figure out an appropriate time step to compromise between accuracy and efficiency. This will depend mostly on the known nature of the body's orbit. The value of dt may be positive for simulating into the future, or negative for simulating into the past. This is a consequence of the time symmetry inherent in Newtonian mechanics.
• GravSimFree(&sim) frees all memory for a simulation stepper object that has been initialized by GravSimInit. (This will only be needed for the C version. The other supported languages provide automatic memory management.)

This is the bare minimum of a useful gravity simulator. More features could be added later. The prime candidate I can think of so far is:

• GravSimExtrap(&sim, time) would provide a fast extrapolation of the small body's state vector at time near the current time step, i.e. smaller than the chosen dt value. This could be used for smooth animation or plotting of a trajectory between the discrete time steps.

Testing and Validation

To verify the simulator is working correctly, I will use JPL Horizons to generate test data containing state vectors for major asteroids and/or comets. I would like a mix of asteroids with low-eccentricity orbits (like Ceres, Pallas, Vesta, and Juno) and a couple of comets having more eccentric orbits.

Limitations

This will not be the most general possible gravity simulator. It will be designed primarily for asteroids and comets that never pass very close to any planet. Therefore, it will not work for calculating planetary satellite orbits.

By default I will model gravitational effects of the Earth/Moon system by ignoring the Moon as a separate object and adding its mass to the Earth's mass.

In fact, there might be an option passed into GravSimInit that specifies whether the inner planets (Mercury, Venus, Earth, Mars) should be included at all. This would optimize execution speed for outer solar system bodies, with negligible effects on accuracy of the results.

In the future I could add another option to calculate the Earth and Moon separately. (Calculating the Moon's state vector is very time consuming.) So I will design the interface to allow for that option to be added later, without breaking the existing API. This would allow calculation of artificial satellite orbits or trajectories of Earth-grazing objects.

[cosinekitty]

A couple more notes about the gravity simulator design:

• All positions and velocities will be with respect to the Solar System Barycenter (SSB), because it is an inertial frame of reference, which makes the Newtonian mechanics simplest.
• At each step of the calculation, because I have to internally calculate the position of the Sun and planets, it makes sense to provide an additional method that can extract the barycentric position of the Sun or a planet from the simulation data structure. The data structure is opaque, but calling code should be able to save significant overhead if it also needs to know where the Sun or planet(s) are.
• As an example of why this would be helpful, getting the Sun's barycentric state can be used to convert the simulated small body's barycentric state into heliocentric coordinates. Similarly, availability of the Earth's barycentric state can convert the body's state into geocentric coordinates.

[cosinekitty]

Progress update

In the new gravsim branch, I have a prototype that simulates the movement of the asteroids Ceres, Pallas, Vesta, and Juno. It does so over a 20-year time span from 2020 to 2040, and compares the calculations with test data from JPL Horizons. The results are encouraging. There is less than 0.6 arcminutes of position error over that 20-year span.

Accuracy considerations

The 0.6 arcminutes of error persists no matter how small I make the time increment, which suggests that the error comes not from the finite time increment but from something else. My prime suspect is the fact that my planet position calculations are already approximated due to truncated VSOP87 series. If this is true, the most serious source of error would be Jupiter, which has more mass than all the other planets combined. Later today I will do an experiment to see if using more VSOP87 terms for Jupiter causes a decrease in the error for the asteroids. If so, it may be worth trading a modest increase in code size and execution time for better accuracy of the gravity simulator.

Another possible source of error is that even when I include the inner 4 planets (Mercury, Venus, Earth, Mars) in the gravity simulation, I do not use them to adjust wobble of the Sun (that is, the SSB). I have been assuming these planets are too small and too close to the Sun to matter. But I may also experiment with using them to refine the SSB position to see if that decreases error in asteroid positions.

Design change: multi-body simulation

Because I have to calculate the positions of the Sun and planets with respect to the Solar System Barycenter on each time step, I realized anyone who wanted to calculate more than one small body would be burdened with redundant calculations. Therefore, I changed the design to allow calculating any number of small bodies simultaneously with each simulation step. Once I add the additional method to retrieve the state vector for the Sun and planets from the simulator object, this will enable an efficient simulation of the entire Solar System, with as many additional bodies as desired.

Inconvenience of barycentric inputs/outputs?

To use the gravity simulator, the calling code must pass in the initial state of each body. This means you have to provide the position vector and velocity vector for each body. In the current design, these vectors are all with respect to the Solar System Barycenter (SSB). It occurs to me that this may be inconvenient for developers. If you have heliocentric coordinates, in order to convert them to barycentric coordinates, you would have to call HelioVector and pass it BODY_SSB to find the barycentric position of the Sun. Then you would have to add the Sun's barycentric vectors to the body's heliocentric vectors to obtain the body's barycentric vectors.

This is not so bad, but it is redundant, because the GravSimInit function already has to calculate the SSB.

So now I'm thinking it may be helpful to provide another parameter to the gravsim functions that specifies what coordinate system to use: barycentric, heliocentric, or geocentric. Then the inputs would be provided in the selected system, and the outputs would be reported in that same system. Not only will this be more efficient, it will simplify using the code when state vectors are known in a non-barycentric reference frame.

[cosinekitty]

Completed in b4f485eeccd0bf61c367a93ade78370710f71d0c