Orbit Determination

Orbit determination (OD) is the inverse of propagation: instead of computing where a satellite will be given its orbital elements, you compute the orbital elements given where the satellite was observed. pg_orrery’s OD solver fits SGP4/SDP4 mean elements (as a TLE) from observations using differential correction with an SVD least-squares solve.

Three observation types are supported:

• ECI — position/velocity vectors in the TEME frame
• Topocentric — azimuth, elevation, range (and optionally range rate) from a ground station
• Angles-only — right ascension and declination (RA/Dec) from optical observations

How differential correction works

The solver starts with an initial orbit estimate (the “seed”) and iteratively refines it to minimize the residuals between observed and computed positions. Each iteration:

1. Propagates the current TLE to each observation time
2. Computes residuals (observed minus computed)
3. Builds the Jacobian matrix (partial derivatives of residuals with respect to orbital elements)
4. Solves the normal equations via SVD to get a correction vector
5. Applies the correction to the equinoctial elements and converts back to a TLE

The solver uses equinoctial elements rather than classical Keplerian elements to avoid singularities at zero eccentricity and zero inclination — both common for real satellites.

ECI fitting

The simplest case: you have ECI position/velocity vectors (perhaps from another propagator, a simulation, or a precise ephemeris) and want to fit a TLE.

-- Generate synthetic observations by propagating a known TLE
WITH iss AS (
    SELECT '1 25544U 98067A   24001.50000000  .00016717  00000-0  10270-3 0  9025
2 25544  51.6400 208.9163 0006703  30.1694  61.7520 15.50100486 00001'::tle AS tle
),
obs AS (
    SELECT array_agg(sgp4_propagate(iss.tle, t) ORDER BY t) AS positions,
           array_agg(t ORDER BY t) AS times
    FROM iss,
         generate_series(
             '2024-01-01 12:00:00+00'::timestamptz,
             '2024-01-01 13:30:00+00'::timestamptz,
             interval '5 minutes') t
)
SELECT iterations,
       round(rms_final::numeric, 6) AS rms_km,
       status,
       round(condition_number::numeric, 1) AS cond
FROM obs, tle_from_eci(obs.positions, obs.times);

With clean synthetic data (propagated from the same TLE model), this converges in a few iterations with sub-meter RMS. With real observations containing measurement noise, expect RMS in the 0.1-10 km range depending on data quality.

No seed needed

tle_from_eci() can bootstrap its own seed using the Gibbs method (three position vectors to derive an initial velocity). You only need to provide a seed TLE when the Gibbs geometry is degenerate (nearly colinear positions) or when you want to start from a known approximate orbit.

Topocentric fitting

When observations come from a ground station (radar tracks, optical telescopes with range data), use tle_from_topocentric(). The solver accounts for the observer’s position on the rotating Earth when computing the observation geometry.

-- Observe a satellite, then fit from the topocentric data
WITH iss AS (
    SELECT '1 25544U 98067A   24001.50000000  .00016717  00000-0  10270-3 0  9025
2 25544  51.6400 208.9163 0006703  30.1694  61.7520 15.50100486 00001'::tle AS tle
),
obs AS (
    SELECT array_agg(
               observe(iss.tle, '40.0N 105.3W 1655m'::observer, t) ORDER BY t
           ) AS observations,
           array_agg(t ORDER BY t) AS times
    FROM iss,
         generate_series(
             '2024-01-01 12:00:00+00'::timestamptz,
             '2024-01-01 12:10:00+00'::timestamptz,
             interval '30 seconds') t
)
SELECT iterations,
       round(rms_final::numeric, 4) AS rms_km,
       status
FROM iss, obs,
     tle_from_topocentric(
         obs.observations, obs.times,
         '40.0N 105.3W 1655m'::observer,
         seed := iss.tle
     );

Seed usually required

Unlike ECI fitting, topocentric fitting typically needs a seed TLE. The azimuth/elevation/range → orbit mapping is more nonlinear than the ECI case, and the solver may not converge from a poor initial guess.

Range rate fitting

When your observations include Doppler or radar range-rate data, enable fit_range_rate to use it as an additional constraint:

SELECT iterations, round(rms_final::numeric, 4) AS rms_km, status
FROM tle_from_topocentric(
    observations := topo_obs,
    times := obs_times,
    obs := '40.0N 105.3W 1655m'::observer,
    seed := seed_tle,
    fit_range_rate := true
);

Range rate adds a 4th residual component per observation (alongside azimuth, elevation, and range). This is particularly valuable for short observation arcs where position-only data doesn’t constrain the velocity well. The range-rate residuals are internally scaled by a factor of 10 (1 km/s maps to 10 km equivalent) to balance their magnitude against position residuals.

Multi-observer fitting

When observations come from multiple ground stations, use the multi-observer variant. Each observation is tagged with its originating station via the observer_ids array:

-- Define two ground stations
-- Station 1: Boulder, CO
-- Station 2: Los Angeles, CA

SELECT iterations, round(rms_final::numeric, 4) AS rms_km, status
FROM tle_from_topocentric(
    observations := all_topo_obs,
    times := all_times,
    observers := ARRAY[
        '40.0N 105.3W 1655m',
        '34.1N 118.3W 100m'
    ]::observer[],
    observer_ids := station_ids,    -- e.g., ARRAY[1,1,1,2,2,2]
    seed := seed_tle
);

The observer_ids array uses 1-based indexing into observers[]. Observations from different stations can be interleaved in any order — the solver uses observer_ids[i] to look up the correct station geometry for each observation.

Multi-observer data improves orbit determination geometry, especially for short arcs. Two well-separated stations observing the same pass can constrain an orbit that a single station cannot.

Weighted observations

When observations have different accuracies (e.g., radar range vs. optical angles, or a high-quality station vs. a noisy one), use the weights parameter:

SELECT iterations, round(rms_final::numeric, 4) AS rms_km, status
FROM tle_from_eci(
    positions := obs_positions,
    times := obs_times,
    weights := ARRAY[1.0, 1.0, 0.5, 1.0, 0.2, 1.0]
);

Higher weights give observations more influence on the solution. A weight of 0.5 means “half as trustworthy as weight 1.0.” This is useful for:

• Down-weighting observations near the horizon (higher atmospheric refraction)
• Down-weighting observations from less accurate sensors
• Implementing iterative reweighting after identifying outliers via tle_fit_residuals()

Angles-only (Gauss method)

Optical observations often provide only RA/Dec — no range information. The tle_from_angles() function handles this case, using the Gauss method for initial orbit determination when no seed TLE is available:

SELECT iterations,
       round(rms_final::numeric, 6) AS rms_rad,
       status,
       round(condition_number::numeric, 1) AS cond
FROM tle_from_angles(
    ra_hours := ARRAY[12.345, 12.567, 12.789, 13.012, 13.234, 13.456],
    dec_degrees := ARRAY[45.1, 44.8, 44.3, 43.6, 42.8, 41.9],
    times := ARRAY[
        '2024-01-01 12:00+00', '2024-01-01 12:01+00',
        '2024-01-01 12:02+00', '2024-01-01 12:03+00',
        '2024-01-01 12:04+00', '2024-01-01 12:05+00'
    ]::timestamptz[],
    obs := '40.0N 105.3W 1655m'::observer
);

RMS units

For angles-only fitting, rms_final and rms_initial are in radians, not km. A typical good fit produces RMS around 0.001 radians (~3.4 arcminutes).

The Gauss method requires at least 3 observations with sufficient angular separation and time spread. The method works by solving for the geocentric range at the middle observation time, then deriving a full state vector. This bootstrap orbit is then refined by the standard DC solver.

Gauss limitations

The Gauss method can fail or produce spurious results when:

• Observations span too short an arc (insufficient angular motion)
• The satellite is near zenith for all observations (poor parallax geometry)
• The time gaps between observations are too large (the linear approximation breaks down)

When Gauss fails, provide a seed TLE from a catalog or from a separate IOD method.

RA is in hours [0, 24) and declination in degrees [-90, 90], matching the output convention of star_observe(). This lets you use the same coordinate system for both stellar calibration and satellite observations.

Interpreting results

Every OD function returns the same 8-column record. Here’s how to use each field:

Field	What to check
fitted_tle	NULL means the solver failed completely. Check status for the reason.
iterations	Should be well below max_iter. If it equals max_iter, the solver didn’t converge — increase max_iter or improve the seed.
rms_final	The fit quality. For ECI: sub-km is good. For topocentric: 1-10 km is typical. For angles-only: < 0.01 rad is good.
rms_initial	Compare with rms_final. A large ratio (rms_initial / rms_final) means the solver made significant improvement. If they’re similar, the seed was already good or the solver made little progress.
status	'converged' is success. 'max_iterations' means more iterations needed. Other strings describe specific failures.
condition_number	Measures how well-conditioned the geometry is. Values below ~1e6 are well-conditioned. Above ~1e10 suggests the observations don’t constrain all orbital elements.
covariance	Formal uncertainty. Extract diagonal elements for per-parameter variance. Off-diagonal elements show correlations between parameters.
nstate	6 (equinoctial elements) or 7 (with B*). The covariance matrix is nstate x nstate.

Tips

Seed selection. For ECI fitting, the Gibbs IOD bootstrap usually works. For topocentric and angles-only fitting, start with a TLE from a catalog (Space-Track, CelesTrak) for the target object. The seed doesn’t need to be precise — it just needs to be in the right ballpark (correct orbit regime, approximate inclination).

Minimum observations. The solver needs at least 6 observations for a 6-parameter fit (or 7 for B*). In practice, more is better — 10-20 well-distributed observations across at least one orbit give the solver enough leverage to constrain all elements.

Observation spacing. Spread observations across the orbit. Clustered observations from a short arc constrain some elements well (position) but others poorly (period, eccentricity). A full orbit of data is ideal.

Convergence troubleshooting. If status shows 'max_iterations':

1. Try a better seed TLE (closer to the true orbit)
2. Increase max_iter to 30 or 50
3. Check if observations contain outliers using tle_fit_residuals()
4. Verify that the observation timestamps are correct (off-by-one-hour timezone errors are common)

B fitting.* Only enable fit_bstar when you have observations spanning multiple orbits. B* captures atmospheric drag, which manifests as a slow period decay — a short arc doesn’t have enough signal to constrain it. With insufficient data, fitting B* degrades the solution by introducing an under-determined 7th parameter.