An experimental Python-based angles-only initial orbit determination tool for solar system objects by way of a Lambert solver.

LambertFit takes equatorial angular observations of solar system objects in right ascension and declination and finds an orbit that minimizes the RMS residuals. As the name implies, LambertFit leverages the solution to Lambert's problem to fit the orbit. The conceptual outline of the orbit determination process is as follows:

1. 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.
2. Take the Lambert solution for the singular heliocentric range guess that yields the best (i.e. smallest) residuals as your starting orbit.
3. 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.

There's a little more LambertFit detail on my website including plots for solves for the first 99 numbered objects for various observation durations and counts.

pip install git+https://github.com/bengebre/lambertfit

I've tested the pip install method above with recent Anaconda installs, but if you're having issues, this is the environment LambertFit was developed and tested in.

• python=3.8.8
• numpy=1.20.1
• astroquery=0.4.3
• astropy=4.2.1
• poliastro=[0.16.0,0.17.0]

import numpy as np
import lambertfit as lf
from astroquery.jplhorizons import Horizons
from astropy.time import Time

#define observer location and body to observe
loc = 'G96' #MPC observatory location code (eg. 'G96' or '500@0')
body_id = '2000002' #small body designation (e.g. '2000002' == Pallas)

#how many (simulated) observations do we want and at what time
tstart = 2459000.5 #julian start date
dt = 14 #number of days of observations
nobs = 30 #number of observations

#observation times
times = np.linspace(tstart,tstart+dt,nobs) #julian date observation times
tdbs = Time(times,format='jd').tdb.value #barycenter dynamical times for the observations

#get simulated observation data from Horizons
ephs = Horizons(id=body_id,location=loc,epochs=times,id_type='designation').ephemerides()
radecs = np.array([[eph['RA'],eph['DEC']] for eph in ephs]) #RA and DEC in degrees

#sun->observer position vector at observation times
r_so = lf.sun_observer(loc,tdbs)

#LambertFit OD solution for RA/DEC observations
rvm_sb, rvn_sb, fit_rms = lf.fit(radecs,times,r_so)

Note: rvm_sb and rvn_sb are heliocentric position and velocity vectors of the LambertFit solution with units of AU and km/s at the endpoints (the first and last observation by default). fit_rms is the final RMS error of the orbit fit in arc seconds.

LambertFit can be quite slow to converge. By default it stops solving when the step size reaches 100km. To prevent LamberFit from running for a half an hour while it chips away at that 5th decimal place of RMS, you can use the following optional arguments on the fit() function:

• tol: The minimum improvement (in arc seconds) at which the solver stops solving. Default: 1e-5
• max_iter: The maximum number of iterations the solver will attempt. Default: 100,000
• min_rms: When the solver achieves this RMS threshold it will stop. Default: 0.001

To specify different endpoints to fit around than the first and last observations:

• pts: The indicies of the observation endpoints that will be fit around. Default: [0,-1]

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.

• LambertFit uses two body propagation around the solar system barycenter. In my testing, solving for observation durations greater than 90 days or so starts to be limited by the fidelity of the orbit propagator.
• LambertFit requires that you choose two observations (which I call endpoints) to solve around. By default these are the first and last observations. More than any other observations, the quality of the orbit LambertFit finds is dependent on how good these two endpoint observations are. The orbit is in a sense 'pinned' at these two endpoint observations and fit to the remaining observations.
• LambertFit is SLOW. I'm essentially using the dumbest implementation of gradient descent possible right now. I'm optimistic that there is room for improvement here.

• LambertFit uses Poliastro for both orbit propagation and for Lambert solving (Izzo, 2014). It's an excellent astrodynamics package.
• While my approach is not the same as the authors, I did find the paper Initial orbit determination methods for track-to-track association (Pastor et. al, 2021) very helpful for thinking about how to use a Lambert solver for initial orbit determination.