Orbit Determination Program (ODP) vs ORBIT14 vs ESOC: What Public Papers Say, Plus a Python Batch Fit

ORBIT14 is the University of Pisa's orbit determination code, not ESA's, and ESA's operations centre (ESOC) runs a separate system. This post compares NASA/JPL's orbit determination program (ODP), Orbit14 and ESOC's system strictly from their published papers, then runs the math they share in a Python script and turns the output into a heuristic for choosing between normal equations, QR and a sequential filter. I had no access to any of the three programs and benchmarked none of them; every number about the demo is from a run on this machine.

ODP vs ORBIT14 vs ESOC: what each program is

Orbit14 comes from the Celestial Mechanics Group at the University of Pisa (with the spin-off SpaceDyS). Its reference paper, Lari et al., "Orbit determination methods for interplanetary missions: development and use of the Orbit14 software", says in its Introduction (Sec. 1) that the code is written in Fortran90 and was built for radio-science data from BepiColombo and, from 2011, Juno. The same section says "its source code is fully accessible to our research group, so that it can be directly modified and updated", in contrast to "other OD programs, where generally users have a software license". The paragraph just before names JPL's MONTE and CNES's GINS as examples of sophisticated OD software, but it does not say which programs the license remark refers to. The paper doesn't benchmark against either, and it doesn't mention ODP.

ESOC's system is described in Budnik et al., ISSFD 2004. The paper doesn't name the program, so neither do I.

JPL's ODP: the public source I read in full is Moyer, JPL TR 32-1527 (1971) (I read the PDF, not just the catalogue page). Its summary calls the Double-Precision Orbit Determination Program (DPODP) "a third-generation program", developed from 1964 to 1968 and first used operationally on Mariner VI and VII. It is an early report, so the DPODP column below describes the program as of 1971, not today's ODP. The abstract of the 1989 NASA Tech Brief says the programs DPTRAJ and ODP were used "extensively in support of NASA's Voyager project" and lists a UNIVAC version and a VAX/VMS version. I read only the catalogue abstract of the 1974 intro report.

Two facts from these papers matter more than the table:

1. ESOC built on ODP's published math. The ESOC paper says its starting point was the program used for Giotto in the 1980s, and that "the mathematical formulation of the Orbit Determination Program (ODP) from NASA/JPL served as the main reference document" (Sec. 1).
2. They cross-checked it against ODP (Sec. 3.1, "Cross-Verification Tests with JPL"). For a 90-day heliocentric cruise, end-of-arc position differences were "ca. 1m in all 5 tests". For 20 days around a Rosetta Mars swing-by they were "ca. 44 m", which the authors attribute to sensitivity to the pericentre state (differences there: about 0.1 m). For a 2-day Mars orbit they were "less than 2.2 mm". With each agency's normal operational set-up on Stardust data, Table 2 (Sec. 3.1) lists 55.8 mm for two-way range and 0.003 mm/s for two-way Doppler. These are ESOC-versus-JPL agreement figures on chosen scenarios, not an accuracy ranking, and they say nothing about Orbit14.

ESOC's GODOT page says "While still in active development, GODOT is also already free to use for external entities in ESA Member States" and that Solar Orbiter's most recent Venus flyby orbit determination "was used to successfully validate GODOT". The page is undated, so check its current terms.

Batch least squares vs Kalman filter for orbit determination

All of these programs do the same loop: propagate a state with a force model, predict measurements, compare with real ones, adjust parameters to shrink the weighted residuals.

• Batch least squares. Linearize around a guess, solve the normal equations for a correction, iterate. The Pisa paper calls it a differential correction.
• SRIF / square-root form. The ESOC paper (Sec. 2.1) says the SRIF "is mathematically exactly equivalent to weighted least squares but it is numerically superior"; Moyer says the same of the square-root form (quoted in the table). I test the QR version of this claim below.
• Kalman filtering. Sequential updating. With no process noise and an epoch state, it should land on the batch answer. Tested below too.
• Force models. Central gravity plus perturbations. Whatever you leave out shows up as residual signal.

A runnable batch orbit fit from range and angles

Save the script as od_demo.py. It simulates a 7000 km, 51.6° circular orbit, a fixed point on Earth's surface as the station (no Earth rotation), and 51 observations over 50 minutes. Measurements are range with 10 m noise and right ascension/declination with 20 µrad noise. The fit starts 5 km and 5 m/s off, using SciPy's DOP853 integrator and finite-difference Jacobians. Cases A, B and B2 use noise seed 42; the arc-length table uses 50 other seeds (1000-1049). Run it with python -I od_demo.py.

Environment and tolerance. I ran it on Python 3.12.10 on Windows (NumPy 2.5.3, SciPy 1.18.1, printed in the output). The whole script took 75 s there, nearly all of it the 50-seed table. Different NumPy/SciPy versions or platforms can shift trailing digits, so small differences in the costs, metres and condition numbers are normal; I did not test other versions, so I can't say how large the shift gets.

Limits of this demo. A fixed, non-rotating station and a 50-minute arc make the problem far better conditioned than real tracking, where Earth rotation and multi-day arcs change the geometry. I did not test those.

import sys
import numpy as np
import scipy
from scipy.integrate import solve_ivp

MU = 398600.4418; RE = 6378.137; J2 = 1.08262668e-3     # km^3/s^2, km, -
SEED = 42                                              # noise seed for cases A, B, J2-fit
SIG = np.array([0.010, 2e-5, 2e-5])                    # range 10 m, RA/Dec 20 urad

def acc(r, j2):
    rn = np.linalg.norm(r); a = -MU*r/rn**3
    if j2:
        z2 = (r[2]/rn)**2; k = 1.5*J2*MU*RE**2/rn**5
        a += k*np.array([r[0]*(5*z2-1), r[1]*(5*z2-1), r[2]*(5*z2-3)])
    return a

def prop(x0, ts, j2=False):
    f = lambda t, y: np.r_[y[3:], acc(y[:3], j2)]
    s = solve_ivp(f, (0, ts[-1]), x0, t_eval=ts, rtol=1e-11, atol=1e-12, method="DOP853")
    return s.y.T

STN = np.array([RE*np.cos(np.radians(40)), 0, RE*np.sin(np.radians(40))])  # fixed inertial point

def meas(X):                                           # range, right ascension, declination
    d = X[:, :3]-STN; rho = np.linalg.norm(d, axis=1)
    return np.c_[rho, np.arctan2(d[:, 1], d[:, 0]), np.arcsin(d[:, 2]/rho)]

r0 = 7000.; v0 = np.sqrt(MU/r0); inc = np.radians(51.6)
x_true = np.array([r0, 0, 0, 0, v0*np.cos(inc), v0*np.sin(inc)])
ts = np.arange(0, 3000.1, 60.)                         # 51 observations over 50 min
wrap = lambda a: (a+np.pi) % (2*np.pi)-np.pi

def run(rng, j2_truth, j2_model=False, nobs=len(ts), verbose=True):
    t = ts[:nobs]
    Z = meas(prop(x_true, t, j2_truth)) + rng.normal(size=(nobs, 3))*SIG
    x = x_true + np.array([5., -5., 3., 0.005, -0.004, 0.003])   # start off the truth
    W = 1/SIG; cost_old = np.inf
    for it in range(8):
        h = meas(prop(x, t, j2_model))
        res = Z-h; res[:, 1:] = wrap(res[:, 1:])
        J = np.zeros((3*nobs, 6))
        for k in range(6):                             # finite-difference partials
            dx = np.zeros(6); dx[k] = 1e-3 if k < 3 else 1e-6
            dd = meas(prop(x+dx, t, j2_model))-h; dd[:, 1:] = wrap(dd[:, 1:])
            J[:, k] = ((dd/dx[k])*W).ravel()
        r = (res*W).ravel(); cost = np.sum(r**2)
        dxs = np.linalg.solve(J.T@J, J.T@r)
        if verbose: print(f"  iter {it+1}: cost={cost:.4f}  |dx|={np.linalg.norm(dxs):.3e}")
        if abs(cost_old-cost) < 1e-6*cost: break       # stop on cost change, not |dx|
        cost_old = cost
        x = x + dxs
    return x, np.linalg.inv(J.T@J), J, r

def report(x, P, r):
    e = x-x_true; sd = np.sqrt(np.diag(P))
    print(" |pos err| m = %.2f ; formal 1-sigma pos (m) = %s ; chi2/dof = %.3f"
          % (np.linalg.norm(e[:3])*1e3, np.round(sd[:3]*1e3, 2), np.sum(r**2)/(len(r)-6)))

if __name__ == "__main__":
    import time; t0 = time.time()
    print("Python %s, NumPy %s, SciPy %s" % (sys.version.split()[0], np.__version__, scipy.__version__))
    print("Case A: truth = two-body, model = two-body (seed %d)" % SEED)
    x, P, J, r = run(np.random.default_rng(SEED), False); report(x, P, r)
    dxn = np.linalg.solve(J.T@J, J.T@r)
    Q, R = np.linalg.qr(J); dxq = np.linalg.solve(R, Q.T@r)
    print(" cond(J)=%.3e cond(J^T J)=%.3e ; QR vs normal-eq step diff = %.3e"
          % (np.linalg.cond(J), np.linalg.cond(J.T@J), np.linalg.norm(dxq-dxn)))
    Pk = np.diag([100.**2]*3+[0.1**2]*3); xk = np.zeros(6)   # sequential filter, no process noise
    for i in range(len(ts)):
        Hi = J[3*i:3*i+3]; zi = r[3*i:3*i+3]
        K = [email protected]@np.linalg.inv(Hi@[email protected]+np.eye(3))
        xk = xk+K@(zi-Hi@xk); Pk = (np.eye(6)-K@Hi)@Pk
    print(" sequential filter vs batch correction diff: %.3e" % np.linalg.norm(xk-dxn))
    print("\nCase B: truth includes J2, model = two-body")
    x, P, J, r = run(np.random.default_rng(SEED), True, False, verbose=False); report(x, P, r)
    print("\nCase B2: truth includes J2, model includes J2 (same seed)")
    x, P, J, r = run(np.random.default_rng(SEED), True, True, verbose=False); report(x, P, r)
    print("\nCase C: arc length, mean over 50 seeds (two-body truth)")
    for n in (5, 10, 20, 51):
        act, frm, cnd = [], [], []
        for s in range(50):
            x, P, J, r = run(np.random.default_rng(1000+s), False, nobs=n, verbose=False)
            act.append(np.linalg.norm((x-x_true)[:3])*1e3); frm.append(np.sqrt(np.trace(P[:3, :3]))*1e3)
            cnd.append(np.linalg.cond(J))
        act, frm = np.array(act), np.array(frm)
        print(" nobs=%2d arc=%2.0f min  actual=%6.2f +/- %.2f m  formal=%6.2f m  actual/formal=%.2f  mean cond(J)=%.2e"
              % (n, ts[n-1]/60, act.mean(), act.std(ddof=1)/np.sqrt(50), frm.mean(), act.mean()/frm.mean(), np.mean(cnd)))
    print("runtime: %.0f s" % (time.time()-t0))

Output (from od_output.txt, the run described above):

Python 3.12.10, NumPy 2.5.3, SciPy 1.18.1
Case A: truth = two-body, model = two-body (seed 42)
  iter 1: cost=31668679.1523  |dx|=7.729e+00
  iter 2: cost=467.8266  |dx|=7.267e-02
  iter 3: cost=112.3031  |dx|=5.393e-06
  iter 4: cost=112.3031  |dx|=9.095e-06
 |pos err| m = 10.29 ; formal 1-sigma pos (m) = [10.49 13.38  9.84] ; chi2/dof = 0.764
 cond(J)=1.079e+04 cond(J^T J)=1.164e+08 ; QR vs normal-eq step diff = 1.788e-17
 sequential filter vs batch correction diff: 4.296e-11

Case B: truth includes J2, model = two-body
 |pos err| m = 2786.25 ; formal 1-sigma pos (m) = [10.5  13.39  9.84] ; chi2/dof = 595.794

Case B2: truth includes J2, model includes J2 (same seed)
 |pos err| m = 10.28 ; formal 1-sigma pos (m) = [10.49 13.39  9.84] ; chi2/dof = 0.764

Case C: arc length, mean over 50 seeds (two-body truth)
 nobs= 5 arc= 4 min  actual= 67.95 +/- 5.81 m  formal= 72.86 m  actual/formal=0.93  mean cond(J)=2.13e+03
 nobs=10 arc= 9 min  actual= 43.62 +/- 3.80 m  formal= 50.59 m  actual/formal=0.86  mean cond(J)=4.40e+03
 nobs=20 arc=19 min  actual= 30.69 +/- 2.70 m  formal= 35.12 m  actual/formal=0.87  mean cond(J)=7.35e+03
 nobs=51 arc=50 min  actual= 20.18 +/- 1.66 m  formal= 19.65 m  actual/formal=1.03  mean cond(J)=1.08e+04
runtime: 75 s

Case A: two iterations took the cost from 3.2e7 to 468, and the third to 112 (chi²/dof 0.764). The true position error of 10.3 m sits inside the formal 1σ values of roughly 10-13 m per axis, which is what a healthy fit looks like.

Gotcha 1: a step-size stop criterion never fires

Stopping when |dx| drops below a tight tolerance such as 1e-9 doesn't work here. The printed |dx| falls to 5.4e-6 at iteration 3 and is 9.1e-6 at iteration 4, so it plateaus instead of going to zero, while the cost is unchanged at 112.3031. I didn't test why; finite-difference Jacobian error is the likely suspect, but that is a guess. The script stops on the change in cost instead. For production codes, the Pisa paper says it integrates the partial derivatives of the dynamics and the ESOC paper says it integrates the variational equations alongside the state.

Gotcha 2: formal uncertainty lies when the force model is wrong

Case B puts J2 in the truth but fits with two-body: position error 2786 m, chi²/dof 596, while the formal 1σ still says about 10 m. Case B2 adds the J2 term to the fitted model, same seed, and chi²/dof returns to 0.764 with a 10.28 m error. Check chi²/dof before believing a covariance. The ESOC paper states a related requirement from the other side: modelled observations should be "approximately one order of magnitude better than the accuracy of the actual measurements" (Sec. 2.3).

Gotcha 3: actual versus formal error needs averaging

One noise draw is not a comparison, so the table averages 50 seeds and prints the standard error of the mean actual error. Actual/formal is 0.93 at 5 observations, 0.86 at 10, 0.87 at 20 and 1.03 at 51. So in my test the formal sigma over-estimated the real error by 7-14% on the shorter arcs and was about right on the longest. The mean actual errors carry standard errors of 1.7-5.8 m, so the ratios are rough. This is with an exactly right model (Gotcha 2); I did not test anything else.

Which solver: a heuristic from one case

• Normal equations vs QR/SRIF. QR and normal equations gave the same step (difference 1.8e-17), but forming JᵀJ squared the condition number, from cond(J)=1.1e4 to cond(JᵀJ)=1.2e8. Double precision carries about 16 digits, so by simple arithmetic that costs roughly 8 of them here. My heuristic: compute cond(J) once; stay with normal equations while it is small, and switch to QR when it is large or you can't tell. It rests on this one setup, so don't treat 1e4 as a threshold. The arc-length table shows cond(J) moving with the setup: mean 2.1e3 at 5 observations up to 1.1e4 at 51. I did not run Earth rotation or multi-day arcs, so I can't say how much higher it goes there. I tested the QR identity only, not a full SRIF with consider parameters.
• Batch vs sequential. With no process noise, my sequential filter landed within 4.3e-11 of the batch correction (a priori sigma 100 km, 100 m/s). With a deterministic model and an epoch state the answers match, so the choice rests on whether you need estimates as data arrive. I did not test process noise or maneuvers.

Is ODP or ORBIT14 available to download?

• ODP: the sources I read don't say how to obtain it or on what terms. Ask JPL.
• Orbit14: per its paper, the source is accessible to the Pisa group; whether outsiders can get it isn't stated. Contact the authors.
• ESOC: the GODOT page says it's free for external entities in ESA Member States, with the undated-page caveat above.
• Learning the math or a small problem: write it yourself as above. I haven't evaluated any open libraries.

None of the sources I read compares the three programs on the same data, so any accuracy ranking between them would be invented, and I found no pricing for any of them.

What to do next

1. Replace the finite-difference Jacobian with integrated variational equations, as both papers describe.
2. Add process noise to the filter loop and see where it stops matching the batch answer.
3. Read the ESOC paper's Sec. 3.1 if you need a template for validating your own OD code against a trusted reference.