Functions: Orbit Determination

Functions for fitting TLE orbital elements from observations via differential correction (DC). The solver uses equinoctial elements internally to avoid singularities at zero eccentricity and inclination, with LAPACK SVD (dgelss_) for the least-squares solve.

Three observation types are supported: ECI position/velocity, topocentric (az/el/range with optional range rate), and angles-only (RA/Dec). Each has single-observer and multi-observer variants where applicable.

All OD functions share the same 8-column output record.

---

Shared Output Record

Every OD function returns a RECORD with these fields:

Field	Type	Description
fitted_tle	tle	The fitted TLE. NULL if the solver did not converge.
iterations	int4	Number of DC iterations performed.
rms_final	float8	RMS residual after final iteration. Units depend on the observation type (km for ECI/topocentric, radians for angles-only).
rms_initial	float8	RMS residual before the first iteration. Compare with rms_final to assess improvement.
status	text	Convergence status: 'converged', 'max_iterations', or an error description.
condition_number	float8	Condition number of the normal equations matrix. Values above ~1e10 suggest poorly-constrained geometry.
covariance	float8[]	Formal covariance matrix (H^T H)^{-1} from the final Jacobian, stored as a flat array in row-major order. Length is nstate * nstate.
nstate	int4	Number of estimated parameters (6 for equinoctial elements, 7 if fit_bstar is true).

Covariance interpretation

The covariance matrix is formal — it reflects the linear approximation at the solution and does not account for systematic errors in the observation model. It is optimistic. Use it for relative comparisons between solutions, not as an absolute uncertainty bound.

---

tle_from_eci

Fit a TLE from ECI position/velocity observations. This is the simplest OD mode — the observations are already in the SGP4 propagation frame (TEME), so no observer geometry is involved.

Signature

tle_from_eci(
    positions eci_position[],
    times timestamptz[],
    seed tle DEFAULT NULL,
    fit_bstar boolean DEFAULT false,
    max_iter int4 DEFAULT 15,
    weights float8[] DEFAULT NULL,
    OUT fitted_tle tle,
    OUT iterations int4,
    OUT rms_final float8,
    OUT rms_initial float8,
    OUT status text,
    OUT condition_number float8,
    OUT covariance float8[],
    OUT nstate int4
) RETURNS RECORD

Parameters

Parameter	Type	Default	Description
positions	eci_position[]	—	Array of ECI position/velocity observations. Requires >= 6 observations.
times	timestamptz[]	—	Observation timestamps. Must be same length as positions.
seed	tle	NULL	Initial TLE estimate. When NULL, the solver uses the Gibbs or Herrick-Gibbs method to bootstrap an initial orbit from three position vectors.
fit_bstar	boolean	false	When true, estimates B* drag coefficient as a 7th parameter. Requires a well-distributed observation arc.
max_iter	int4	15	Maximum differential correction iterations.
weights	float8[]	NULL	Per-observation weights. NULL means uniform weighting. Length must equal length of positions. Higher weight = more influence on the solution.

Seed-free fitting

When no seed TLE is provided, the solver uses the Gibbs method (or Herrick-Gibbs for closely-spaced observations) to derive an initial velocity estimate from three position vectors. This makes tle_from_eci() fully seed-free for most use cases.

Example

-- Round-trip test: propagate a TLE, then fit it back
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,
       round(rms_initial::numeric, 3) AS rms_init_km,
       status,
       round(condition_number::numeric, 1) AS cond
FROM obs, tle_from_eci(obs.positions, obs.times);

---

tle_from_topocentric (single observer)

Fit a TLE from topocentric (azimuth/elevation/range) observations collected by a single ground station.

Signature

tle_from_topocentric(
    observations topocentric[],
    times timestamptz[],
    obs observer,
    seed tle DEFAULT NULL,
    fit_bstar boolean DEFAULT false,
    max_iter int4 DEFAULT 15,
    fit_range_rate boolean DEFAULT false,
    weights float8[] DEFAULT NULL,
    OUT fitted_tle tle,
    OUT iterations int4,
    OUT rms_final float8,
    OUT rms_initial float8,
    OUT status text,
    OUT condition_number float8,
    OUT covariance float8[],
    OUT nstate int4
) RETURNS RECORD

Parameters

Parameter	Type	Default	Description
observations	topocentric[]	—	Array of topocentric observations (az, el, range, range_rate). Requires >= 6 observations.
times	timestamptz[]	—	Observation timestamps.
obs	observer	—	Ground station location.
seed	tle	NULL	Initial TLE estimate. Unlike ECI fitting, topocentric fitting typically needs a seed for convergence.
fit_bstar	boolean	false	Estimate B* drag coefficient.
max_iter	int4	15	Maximum iterations.
fit_range_rate	boolean	false	When true, includes range rate as a 4th residual component per observation. Use when you have Doppler or radar range-rate data.
weights	float8[]	NULL	Per-observation weights. NULL = uniform.

Range rate scaling

Range rate residuals are internally scaled by OD_RR_SCALE = 10.0, mapping 1 km/s of range-rate error to 10 km equivalent in the residual vector. This prevents the typically small range-rate values from being swamped by position residuals.

Example

-- Observe a satellite, then fit back from 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 obs,
     tle_from_topocentric(
         obs.observations, obs.times,
         '40.0N 105.3W 1655m'::observer,
         seed := iss.tle
     );

---

tle_from_topocentric (multi-observer)

Fit a TLE from topocentric observations collected by multiple ground stations. The observer_ids array maps each observation to its originating station.

Signature

tle_from_topocentric(
    observations topocentric[],
    times timestamptz[],
    observers observer[],
    observer_ids int4[],
    seed tle DEFAULT NULL,
    fit_bstar boolean DEFAULT false,
    max_iter int4 DEFAULT 15,
    fit_range_rate boolean DEFAULT false,
    weights float8[] DEFAULT NULL,
    OUT fitted_tle tle,
    OUT iterations int4,
    OUT rms_final float8,
    OUT rms_initial float8,
    OUT status text,
    OUT condition_number float8,
    OUT covariance float8[],
    OUT nstate int4
) RETURNS RECORD

Parameters

Parameter	Type	Default	Description
observations	topocentric[]	—	All observations from all stations, interleaved in time order.
times	timestamptz[]	—	Observation timestamps.
observers	observer[]	—	Array of ground station locations (1-indexed).
observer_ids	int4[]	—	Per-observation index into observers[]. observer_ids[i] identifies which station produced observations[i].
seed	tle	NULL	Initial TLE estimate.
fit_bstar	boolean	false	Estimate B* drag coefficient.
max_iter	int4	15	Maximum iterations.
fit_range_rate	boolean	false	Include range rate in residuals.
weights	float8[]	NULL	Per-observation weights. Useful when stations have different measurement accuracies.

Example

-- Two ground stations observe the same satellite
SELECT iterations, round(rms_final::numeric, 4) AS rms_km, status
FROM tle_from_topocentric(
    observations := ARRAY[obs1_t1, obs1_t2, obs2_t1, obs2_t2]::topocentric[],
    times := ARRAY['2024-01-01 12:00+00', '2024-01-01 12:05+00',
                    '2024-01-01 12:02+00', '2024-01-01 12:07+00']::timestamptz[],
    observers := ARRAY['40.0N 105.3W 1655m', '34.1N 118.3W 100m']::observer[],
    observer_ids := ARRAY[1, 1, 2, 2],
    seed := seed_tle
);

Observer ID indexing

observer_ids uses 1-based indexing into the observers array. Station 1 is observers[1], station 2 is observers[2], etc. This matches PostgreSQL’s native array indexing convention.

---

tle_from_angles (single observer)

Fit a TLE from angles-only (RA/Dec) observations collected by a single ground station. When no seed TLE is provided, the Gauss method derives an initial orbit from three observations — making this function fully seed-free.

Signature

tle_from_angles(
    ra_hours float8[],
    dec_degrees float8[],
    times timestamptz[],
    obs observer,
    seed tle DEFAULT NULL,
    fit_bstar boolean DEFAULT false,
    max_iter int4 DEFAULT 15,
    weights float8[] DEFAULT NULL,
    OUT fitted_tle tle,
    OUT iterations int4,
    OUT rms_final float8,
    OUT rms_initial float8,
    OUT status text,
    OUT condition_number float8,
    OUT covariance float8[],
    OUT nstate int4
) RETURNS RECORD

Parameters

Parameter	Type	Default	Description
ra_hours	float8[]	—	Right ascension values in hours, range [0, 24). Matches the star_observe() convention.
dec_degrees	float8[]	—	Declination values in degrees, range [-90, 90].
times	timestamptz[]	—	Observation timestamps. Requires >= 3 observations.
obs	observer	—	Ground station location.
seed	tle	NULL	Initial TLE estimate. When NULL, the Gauss method bootstraps an initial orbit from 3 observations.
fit_bstar	boolean	false	Estimate B* drag coefficient.
max_iter	int4	15	Maximum iterations.
weights	float8[]	NULL	Per-observation weights.

RMS units

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

RA/Dec conventions

Right ascension is in hours [0, 24), matching the star_observe() output convention. Declination is in degrees [-90, 90]. Internally, both are converted to radians for the residual computation.

Example

-- Angles-only OD: RA/Dec observations of a satellite
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
);

---

tle_from_angles (multi-observer)

Fit a TLE from angles-only (RA/Dec) observations collected by multiple ground stations. Uses the same Gauss IOD bootstrap as the single-observer variant when no seed is provided.

Signature

tle_from_angles(
    ra_hours float8[],
    dec_degrees float8[],
    times timestamptz[],
    observers observer[],
    observer_ids int4[],
    seed tle DEFAULT NULL,
    fit_bstar boolean DEFAULT false,
    max_iter int4 DEFAULT 15,
    weights float8[] DEFAULT NULL,
    OUT fitted_tle tle,
    OUT iterations int4,
    OUT rms_final float8,
    OUT rms_initial float8,
    OUT status text,
    OUT condition_number float8,
    OUT covariance float8[],
    OUT nstate int4
) RETURNS RECORD

Parameters

Parameter	Type	Default	Description
ra_hours	float8[]	—	Right ascension values in hours [0, 24).
dec_degrees	float8[]	—	Declination values in degrees [-90, 90].
times	timestamptz[]	—	Observation timestamps from all stations, interleaved in time order.
observers	observer[]	—	Array of ground station locations (1-indexed).
observer_ids	int4[]	—	Per-observation index into observers[].
seed	tle	NULL	Initial TLE estimate. NULL = Gauss IOD bootstrap.
fit_bstar	boolean	false	Estimate B* drag coefficient.
max_iter	int4	15	Maximum iterations.
weights	float8[]	NULL	Per-observation weights. Useful when stations have different apertures or sky conditions.

Example

-- Two optical stations observe the same satellite
SELECT iterations, round(rms_final::numeric, 6) AS rms_rad, status
FROM tle_from_angles(
    ra_hours := ARRAY[12.3, 12.5, 12.7, 12.4, 12.6, 12.8],
    dec_degrees := ARRAY[45.0, 44.5, 44.0, 44.8, 44.3, 43.7],
    times := ARRAY[
        '2024-01-01 12:00+00', '2024-01-01 12:02+00',
        '2024-01-01 12:04+00', '2024-01-01 12:01+00',
        '2024-01-01 12:03+00', '2024-01-01 12:05+00'
    ]::timestamptz[],
    observers := ARRAY['40.0N 105.3W 1655m', '34.1N 118.3W 100m']::observer[],
    observer_ids := ARRAY[1, 1, 1, 2, 2, 2]
);

---

tle_fit_residuals

Compute per-observation position residuals between a fitted TLE and the original ECI observations. Returns one row per observation with the XYZ and total position error in km. Use this to identify outlier observations or assess spatial error distribution.

Signature

tle_fit_residuals(
    fitted tle,
    positions eci_position[],
    times timestamptz[]
) RETURNS TABLE (
    t timestamptz,
    dx_km float8,
    dy_km float8,
    dz_km float8,
    pos_err_km float8
)

Parameters

Parameter	Type	Description
fitted	tle	The fitted TLE (typically the fitted_tle output from tle_from_eci()).
positions	eci_position[]	The original ECI observations used for fitting.
times	timestamptz[]	The original observation timestamps.

Returns

One row per observation:

Column	Type	Unit	Description
t	timestamptz	—	Observation time
dx_km	float8	km	X-axis residual (observed - computed)
dy_km	float8	km	Y-axis residual
dz_km	float8	km	Z-axis residual
pos_err_km	float8	km	Total position error: sqrt(dx^2 + dy^2 + dz^2)

Example

-- After fitting, inspect per-observation residuals
WITH fit AS (
    SELECT fitted_tle, positions, times
    FROM obs, tle_from_eci(obs.positions, obs.times)
)
SELECT t,
       round(dx_km::numeric, 4) AS dx,
       round(dy_km::numeric, 4) AS dy,
       round(dz_km::numeric, 4) AS dz,
       round(pos_err_km::numeric, 4) AS total_err
FROM fit, tle_fit_residuals(fit.fitted_tle, fit.positions, fit.times)
ORDER BY pos_err_km DESC;

Outlier detection

Observations with residuals significantly larger than the RMS are potential outliers. Consider removing them and re-fitting, or using the weights parameter to down-weight them.