LambertFit: Python Angles-Only Orbit Determination via Lambert Solver
I often think of HelioLinC as a zeroth-order orbit determination. You guess heliocentric ranges for tracklets and then find tracklets that cluster when you propagate to a common epoch. Eventually you filter these candidate links with real orbit determination. I've used (and will continue to use) find_orb for orbit determination, but I've always thought it'd be interesting to try to leverage a Lambert solver for OD; technically this would probably be considered IOD (initial orbit determination). I thought it might be possible, in this era of fast Lambert solvers (Izzo, 2014), to get pretty good orbits quickly, so I wrote LambertFit to calculate those orbits in Python. The 'quickly' part isn't panning out yet, but I think the solver is interesting enough to merit publishing the code and writeup.
Conceptualization
LambertFit takes equatorial angular observations of solar system objects in right ascension and declination and finds an orbit that minimizes the RMS residuals. Conceptually, LambertFit is probably most closely related to other IOD techniques like Gooding IOD or Double-r integration IOD. As the name implies, LambertFit leverages the solution to Lambert's problem to fit the orbit. The outline of the orbit determination process is as follows:
- Lambert solve an initial orbit between two observation endpoints (e.g. the first and the last observation) by guessing a number of constant, scalar heliocentric ranges at those two endpoints.
- Take the Lambert solution for the singular heliocentric range guess that yields the best (i.e. smallest) residuals as your starting orbit.
- Update each endpoint range estimate, one at a time, by stepping in the direction that reduces the residuals (as measured by the RMS error) until the residuals stop getting smaller.
That's it. Besides being sometimes slow to converge, this technique works pretty well.
Results
The Github repo for LambertFit has some example code that can get you started, so I'll just jump straight into the results here. The plots below show the LambertFit solution for the first 5 numbered small bodies (Ceres, Pallas, Juno, Vesta, Astraea) for varying observation durations and counts. There are also links to the same plots for the first 99 numbered objects below each observation duration/count group.
The observation data that serves as the input RA/DEC data for LamberFit comes from JPL Horizons. The blue orbit is the earth. The true orbit of the body we're trying to fit an orbit to is in green. The left orange orbit is the initial guess orbit that LambertFit starts with (generated internally by LambertFit) and then refines. The LambertFit solution orbit is in orange on the right. The observations from JPL are equally spaced in time between the green diamond (the first) and the green circle (the last). The diamond and circle on the blue orbit represent the earth's position at the starting and ending observation times. The diamond and circle on the orange orbit represent LambertFit's estimate of the position of the solar system object at the starting and ending observation times. The reported RMS errors are in arc seconds.
6 observations over 3 days
30 observations over 14 days
30 observations over 90 days
Observations
Initial guess orbits can have relatively small residuals
I haven't worked on orbit determination algorithms in a while, and I've never implemented an angles-only OD solver before this. With precision orbit determination your observations include measured ranges which the angles-only problem (as the name implies) lacks. Thus for the angles-only technique you have to estimate the range in the solution process. To estimate the range with LambertFit I test a bunch of range guesses for the endpoint observations with a Lambert solver and then calculate the residuals of the orbits it generates. If you look at the 14 day solution for Ceres above, for instance, you can see that the initial guess orbit generated by this process has a ~4.6" RMS error. Not all of the objects have initial guess orbit RMS values that low, but it's interesting that course range guessing can sometimes yield reasonably good fits without any refinement. Of course the number and duration of the observations play a part in this too as the 90 day solves have much higher initial guess RMS values.
Phase angle and solution quality
Another interesting thing to note, when you look at all 99 of the 14 day observation duration solves, objects with a phase angle at the endpoints less than about 60° tend to consistently solve well - both visually and measured by RMS error. This is convenient because solar system objects are usually and best observed at lower phase angles. Observations with phase angles greater than 90° generally aren't even possible for earth-based observers as it would require daytime viewing. I show a bunch of these greater than 90° phase angle objects (and some of them solve well), but they're probably not going to be encountered very much in practice. You can see this with Astraea in the plots above (object 5). It's at roughly a 90 phase angle at the observation times for the 3 and 14 day observation durations and it doesn't have as good of a solution as the other objects. The interesting question that I don't have an answer to is why objects with phase angles around 60-90 degrees don't always solve quite as well. This is probably common knowledge for astronomers, but as I said, this is my first angles-only OD implementation and I don't know the reason yet.
Conclusion
LambertFit turned into an accidental two month detour from my solar system object detection work, but I think I've learned a few things about the angles-only OD process that will fundamentally help with my HelioLinC work and object detection work in general. Will LambertFit replace find_orb in my OD toolbox; I don't think so. But I'm hoping that LambertFit (or the initial orbit guessing algorithm it uses) will help me filter non-physical candidate links more quickly (and in Python!).
I also hope to revisit the LambertFit code in the future to speed up the dumb gradient descent technique I'm using. Publishing this initial release at least lets me walk away with something concrete that I can come back to later rather than abandoning the code in a messy pile because it's not fast enough. So now back to HelioLinC and object detection - hopefully with reports of successful LambertFit filtering added to the mix.
Published: 1/5/2023