warehack.ing/pg_orrery
2
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
pg_orrery/test/expected/od_fit.out
Ryan Malloy bca8b3e7eb Add covariance output and condition number to OD solver (v0.5.0) 
Computes formal covariance (H^T·H)^{-1} via LAPACK dpotrf_/dpotri_
after DC convergence. Returns upper-triangle array (21 elements for
6-state, 28 for 7-state with B*), condition number from SVD, and
nstate count. Covariance is computed even for perfect-seed fits.

Bumps extension to v0.5.0 with full install SQL and migration path.
2026-02-17 16:15:44 -07:00

626 lines
23 KiB
Plaintext
Raw Permalink Blame History

-- od_fit.sql -- Regression tests for TLE fitting (orbit determination)
--
-- Tests tle_from_eci(), tle_from_topocentric(), and tle_fit_residuals().
-- Uses round-trip methodology: propagate known TLE → fit from obs → compare.
SET client_min_messages TO warning;
CREATE EXTENSION IF NOT EXISTS pg_orrery;
SET client_min_messages TO notice;
-- ============================================================
-- Test 1: ECI round-trip (ISS-like LEO orbit)
--
-- Generate 20 observations over 90 minutes from ISS TLE,
-- then fit a TLE from those observations.
-- Expected: RMS < 1 km, convergence.
-- ============================================================
WITH iss_tle AS (
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
),
obs AS (
SELECT
array_agg(sgp4_propagate(t, ts) ORDER BY ts) AS positions,
array_agg(ts ORDER BY ts) AS times
FROM iss_tle,
generate_series(
'2024-01-01 12:00:00+00'::timestamptz,
'2024-01-01 13:30:00+00'::timestamptz,
'5 minutes'::interval
) AS ts
),
result AS (
SELECT (tle_from_eci(positions, times, t)).* FROM obs, iss_tle
)
SELECT
iterations > 0 AS has_iterations,
rms_final < 1.0 AS rms_under_1km,
status = 'converged' AS did_converge
FROM result;
has_iterations | rms_under_1km | did_converge
----------------+---------------+--------------
f | t | t
(1 row)
-- ============================================================
-- Test 2: ECI round-trip (GPS-like MEO orbit)
--
-- Higher altitude, lower mean motion.
-- ============================================================
WITH gps_tle AS (
SELECT E'1 28874U 05038A 24001.50000000 .00000003 00000-0 00000+0 0 9999\n2 28874 55.4000 200.0000 0057000 250.0000 110.0000 2.00567486 00001'::tle AS t
),
obs AS (
SELECT
array_agg(sgp4_propagate(t, ts) ORDER BY ts) AS positions,
array_agg(ts ORDER BY ts) AS times
FROM gps_tle,
generate_series(
'2024-01-01 00:00:00+00'::timestamptz,
'2024-01-01 12:00:00+00'::timestamptz,
'30 minutes'::interval
) AS ts
),
result AS (
SELECT (tle_from_eci(positions, times, t)).* FROM obs, gps_tle
)
SELECT
iterations > 0 AS has_iterations,
rms_final < 1.0 AS rms_under_1km,
status = 'converged' AS did_converge
FROM result;
has_iterations | rms_under_1km | did_converge
----------------+---------------+--------------
t | t | t
(1 row)
-- ============================================================
-- Test 3: 7-state with B* fitting
--
-- Fit ISS with B* as a free parameter.
-- Verify B* is recovered within 50% of original.
-- ============================================================
WITH iss_tle AS (
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
),
obs AS (
SELECT
array_agg(sgp4_propagate(t, ts) ORDER BY ts) AS positions,
array_agg(ts ORDER BY ts) AS times
FROM iss_tle,
generate_series(
'2024-01-01 12:00:00+00'::timestamptz,
'2024-01-01 14:00:00+00'::timestamptz,
'5 minutes'::interval
) AS ts
),
result AS (
SELECT (tle_from_eci(positions, times, t, true, 20)).* FROM obs, iss_tle
)
SELECT
status = 'converged' AS did_converge,
rms_final < 1.0 AS rms_under_1km
FROM result;
did_converge | rms_under_1km
--------------+---------------
t | t
(1 row)
-- ============================================================
-- Test 4: Topocentric round-trip (ISS via observe())
--
-- Observe ISS from MIT ground station, then fit TLE from
-- the topocentric data. Wider tolerance (5 km) due to
-- the topo→ECI information loss.
-- ============================================================
WITH iss_tle AS (
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
),
mit AS (
SELECT '(42.36,-71.09,20)'::observer AS obs
),
topo_obs AS (
SELECT
array_agg(observe(t, obs, ts) ORDER BY ts) AS observations,
array_agg(ts ORDER BY ts) AS times
FROM iss_tle, mit,
generate_series(
'2024-01-01 12:00:00+00'::timestamptz,
'2024-01-01 13:30:00+00'::timestamptz,
'5 minutes'::interval
) AS ts
),
result AS (
SELECT (tle_from_topocentric(observations, times, obs, t, false, 20)).*
FROM topo_obs, mit, iss_tle
)
SELECT
iterations > 0 AS has_iterations,
rms_final < 10.0 AS rms_under_10km,
status = 'converged' AS did_converge
FROM result;
has_iterations | rms_under_10km | did_converge
----------------+----------------+--------------
f | t | t
(1 row)
-- ============================================================
-- Test 5: No seed (auto initial guess from first observation)
--
-- For ECI mode, the solver can derive an initial guess from
-- the first observation without needing a seed TLE.
-- ============================================================
WITH iss_tle AS (
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
),
obs AS (
SELECT
array_agg(sgp4_propagate(t, ts) ORDER BY ts) AS positions,
array_agg(ts ORDER BY ts) AS times
FROM iss_tle,
generate_series(
'2024-01-01 12:00:00+00'::timestamptz,
'2024-01-01 13:30:00+00'::timestamptz,
'5 minutes'::interval
) AS ts
),
result AS (
SELECT (tle_from_eci(positions, times)).* FROM obs
)
SELECT
iterations > 0 AS has_iterations,
rms_final < 5.0 AS rms_under_5km
FROM result;
has_iterations | rms_under_5km
----------------+---------------
t | t
(1 row)
-- ============================================================
-- Test 6: Error - insufficient observations
--
-- Less than 6 observations should raise an error.
-- ============================================================
DO $$
BEGIN
PERFORM tle_from_eci(
ARRAY[
sgp4_propagate(
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
'2024-01-01 12:00:00+00'::timestamptz
)
],
ARRAY['2024-01-01 12:00:00+00'::timestamptz]
);
RAISE NOTICE 'ERROR: should have raised exception';
EXCEPTION
WHEN data_exception THEN
RAISE NOTICE 'OK: insufficient observations error caught';
WHEN invalid_parameter_value THEN
RAISE NOTICE 'OK: insufficient observations error caught';
END $$;
NOTICE: OK: insufficient observations error caught
-- ============================================================
-- Test 7: tle_fit_residuals() diagnostic output
--
-- Verify the residual function returns per-observation data.
-- ============================================================
WITH iss_tle AS (
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
),
obs AS (
SELECT
array_agg(sgp4_propagate(t, ts) ORDER BY ts) AS positions,
array_agg(ts ORDER BY ts) AS times
FROM iss_tle,
generate_series(
'2024-01-01 12:00:00+00'::timestamptz,
'2024-01-01 12:30:00+00'::timestamptz,
'5 minutes'::interval
) AS ts
)
SELECT
count(*) AS n_residuals,
count(*) = 7 AS correct_count,
max(pos_err_km) < 0.001 AS residuals_near_zero
FROM iss_tle, obs, tle_fit_residuals(t, positions, times);
n_residuals | correct_count | residuals_near_zero
-------------+---------------+---------------------
7 | t | t
(1 row)
-- ============================================================
-- Test 8: Nearly circular orbit (equinoctial singularity-free)
--
-- Very low eccentricity tests that equinoctial elements handle
-- the e→0 case without blowing up.
-- ============================================================
WITH circ_tle AS (
SELECT E'1 00001U 24001A 24001.50000000 .00000000 00000-0 00000+0 0 9999\n2 00001 0.0100 000.0000 0000100 000.0000 000.0000 15.00000000 00001'::tle AS t
),
obs AS (
SELECT
array_agg(sgp4_propagate(t, ts) ORDER BY ts) AS positions,
array_agg(ts ORDER BY ts) AS times
FROM circ_tle,
generate_series(
'2024-01-01 12:00:00+00'::timestamptz,
'2024-01-01 13:30:00+00'::timestamptz,
'5 minutes'::interval
) AS ts
),
result AS (
SELECT (tle_from_eci(positions, times, t)).* FROM obs, circ_tle
)
SELECT
rms_final < 2.0 AS rms_under_2km,
status = 'converged' AS did_converge
FROM result;
rms_under_2km | did_converge
---------------+--------------
t | t
(1 row)
-- ============================================================
-- Test 9: Multi-observer topocentric (MIT + Boulder observe ISS)
--
-- Two ground stations observe ISS, fit via multi-observer API.
-- Uses the overloaded tle_from_topocentric(topo[], ts[], observer[], int4[], ...).
-- ============================================================
WITH iss_tle AS (
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
),
stations AS (
SELECT
'(42.36,-71.09,20)'::observer AS mit,
'(40.01,-105.27,1655)'::observer AS boulder
),
topo_obs AS (
SELECT
array_agg(topo ORDER BY rn) AS observations,
array_agg(ts ORDER BY rn) AS times,
array_agg(obs_id ORDER BY rn) AS observer_ids
FROM (
-- MIT observations (observer_id = 0)
SELECT observe(t, mit, ts) AS topo, ts, 0 AS obs_id,
row_number() OVER (ORDER BY ts) AS rn
FROM iss_tle, stations,
generate_series(
'2024-01-01 12:00:00+00'::timestamptz,
'2024-01-01 13:00:00+00'::timestamptz,
'5 minutes'::interval
) AS ts
UNION ALL
-- Boulder observations (observer_id = 1)
SELECT observe(t, boulder, ts) AS topo, ts, 1 AS obs_id,
row_number() OVER (ORDER BY ts) + 100 AS rn
FROM iss_tle, stations,
generate_series(
'2024-01-01 12:00:00+00'::timestamptz,
'2024-01-01 13:00:00+00'::timestamptz,
'5 minutes'::interval
) AS ts
) sub
),
result AS (
SELECT (tle_from_topocentric(
observations, times,
ARRAY[mit, boulder], observer_ids,
t, false, 20
)).*
FROM topo_obs, stations, iss_tle
)
SELECT
status = 'converged' AS did_converge,
rms_final < 10.0 AS rms_under_10km
FROM result;
did_converge | rms_under_10km
--------------+----------------
t | t
(1 row)
-- ============================================================
-- Test 10: Single observer via multi-observer API
--
-- Verify the multi-observer path produces equivalent results
-- when all observations come from one station.
-- ============================================================
WITH iss_tle AS (
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
),
mit AS (
SELECT '(42.36,-71.09,20)'::observer AS obs
),
topo_obs AS (
SELECT
array_agg(observe(t, obs, ts) ORDER BY ts) AS observations,
array_agg(ts ORDER BY ts) AS times,
array_agg(0) AS observer_ids
FROM iss_tle, mit,
generate_series(
'2024-01-01 12:00:00+00'::timestamptz,
'2024-01-01 13:30:00+00'::timestamptz,
'5 minutes'::interval
) AS ts
),
result AS (
SELECT (tle_from_topocentric(
observations, times,
ARRAY[obs], observer_ids,
t, false, 20
)).*
FROM topo_obs, mit, iss_tle
)
SELECT
status = 'converged' AS did_converge,
rms_final < 10.0 AS rms_under_10km
FROM result;
did_converge | rms_under_10km
--------------+----------------
t | t
(1 row)
-- ============================================================
-- Test 11: Error - observer_id out of range
-- ============================================================
DO $$
BEGIN
PERFORM tle_from_topocentric(
ARRAY[
observe(
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
'(42.36,-71.09,20)'::observer,
'2024-01-01 12:00:00+00'::timestamptz
),
observe(
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
'(42.36,-71.09,20)'::observer,
'2024-01-01 12:05:00+00'::timestamptz
),
observe(
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
'(42.36,-71.09,20)'::observer,
'2024-01-01 12:10:00+00'::timestamptz
),
observe(
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
'(42.36,-71.09,20)'::observer,
'2024-01-01 12:15:00+00'::timestamptz
),
observe(
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
'(42.36,-71.09,20)'::observer,
'2024-01-01 12:20:00+00'::timestamptz
),
observe(
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
'(42.36,-71.09,20)'::observer,
'2024-01-01 12:25:00+00'::timestamptz
)
],
ARRAY[
'2024-01-01 12:00:00+00'::timestamptz,
'2024-01-01 12:05:00+00'::timestamptz,
'2024-01-01 12:10:00+00'::timestamptz,
'2024-01-01 12:15:00+00'::timestamptz,
'2024-01-01 12:20:00+00'::timestamptz,
'2024-01-01 12:25:00+00'::timestamptz
],
ARRAY['(42.36,-71.09,20)'::observer],
ARRAY[5, 0, 0, 0, 0, 0], -- observer_id 5 out of range
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle
);
RAISE NOTICE 'ERROR: should have raised exception';
EXCEPTION
WHEN invalid_parameter_value THEN
RAISE NOTICE 'OK: observer_id out of range error caught';
END $$;
NOTICE: OK: observer_id out of range error caught
-- ============================================================
-- Test 12: Topocentric without seed TLE (IOD bootstrap)
--
-- The Gibbs method provides an initial orbit from 3 well-spaced
-- topocentric observations, eliminating the seed requirement.
-- Wider tolerance (~20 km) because the IOD initial guess is rough.
-- ============================================================
WITH iss_tle AS (
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
),
mit AS (
SELECT '(42.36,-71.09,20)'::observer AS obs
),
topo_obs AS (
SELECT
array_agg(observe(t, obs, ts) ORDER BY ts) AS observations,
array_agg(ts ORDER BY ts) AS times
FROM iss_tle, mit,
generate_series(
'2024-01-01 12:00:00+00'::timestamptz,
'2024-01-01 13:30:00+00'::timestamptz,
'5 minutes'::interval
) AS ts
),
result AS (
SELECT (tle_from_topocentric(observations, times, obs)).*
FROM topo_obs, mit
)
SELECT
rms_final < 20.0 AS rms_under_20km,
status IN ('converged', 'max_iterations') AS solver_ran
FROM result;
rms_under_20km | solver_ran
----------------+------------
t | t
(1 row)
-- ============================================================
-- Test 13: ECI without seed, Gibbs-improved initial guess
--
-- With Gibbs, the ECI seedless path uses 3-position geometry
-- instead of a single r,v pair. Should converge.
-- ============================================================
WITH iss_tle AS (
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
),
obs AS (
SELECT
array_agg(sgp4_propagate(t, ts) ORDER BY ts) AS positions,
array_agg(ts ORDER BY ts) AS times
FROM iss_tle,
generate_series(
'2024-01-01 12:00:00+00'::timestamptz,
'2024-01-01 14:00:00+00'::timestamptz,
'5 minutes'::interval
) AS ts
),
result AS (
SELECT (tle_from_eci(positions, times)).* FROM obs
)
SELECT
iterations > 0 AS has_iterations,
rms_final < 5.0 AS rms_under_5km,
status = 'converged' AS did_converge
FROM result;
has_iterations | rms_under_5km | did_converge
----------------+---------------+--------------
t | t | t
(1 row)
-- ============================================================
-- Test 14: Covariance is returned when converged
--
-- A converged 6-state fit should produce a non-NULL covariance
-- array of length 21 (upper triangle of 6x6).
-- ============================================================
WITH iss_tle AS (
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
),
obs AS (
SELECT
array_agg(sgp4_propagate(t, ts) ORDER BY ts) AS positions,
array_agg(ts ORDER BY ts) AS times
FROM iss_tle,
generate_series(
'2024-01-01 12:00:00+00'::timestamptz,
'2024-01-01 13:30:00+00'::timestamptz,
'5 minutes'::interval
) AS ts
),
result AS (
SELECT (tle_from_eci(positions, times, t)).* FROM obs, iss_tle
)
SELECT
covariance IS NOT NULL AS has_covariance,
array_length(covariance, 1) = 21 AS cov_length_21,
nstate = 6 AS is_6state,
condition_number > 0 AS cond_positive
FROM result;
has_covariance | cov_length_21 | is_6state | cond_positive
----------------+---------------+-----------+---------------
t | t | t | t
(1 row)
-- ============================================================
-- Test 15: Covariance diagonal elements positive (variances > 0)
--
-- For a well-conditioned fit, all diagonal elements of the
-- covariance matrix must be positive (they are variances).
-- Upper triangle row-major: diag indices are 0,6,11,15,18,20.
-- ============================================================
WITH iss_tle AS (
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
),
obs AS (
SELECT
array_agg(sgp4_propagate(t, ts) ORDER BY ts) AS positions,
array_agg(ts ORDER BY ts) AS times
FROM iss_tle,
generate_series(
'2024-01-01 12:00:00+00'::timestamptz,
'2024-01-01 13:30:00+00'::timestamptz,
'5 minutes'::interval
) AS ts
),
result AS (
SELECT (tle_from_eci(positions, times, t)).* FROM obs, iss_tle
)
SELECT
covariance[1] > 0 AS var_n_positive,
covariance[7] > 0 AS var_af_positive,
covariance[12] > 0 AS var_ag_positive,
covariance[16] > 0 AS var_chi_positive,
covariance[19] > 0 AS var_psi_positive,
covariance[21] > 0 AS var_L0_positive
FROM result;
var_n_positive | var_af_positive | var_ag_positive | var_chi_positive | var_psi_positive | var_l0_positive
----------------+-----------------+-----------------+------------------+------------------+-----------------
t | t | t | t | t | t
(1 row)
-- ============================================================
-- Test 16: 7-state covariance (with B* fitting)
--
-- When fit_bstar = true, covariance should have 28 elements
-- (upper triangle of 7x7) and nstate = 7.
-- ============================================================
WITH iss_tle AS (
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
),
obs AS (
SELECT
array_agg(sgp4_propagate(t, ts) ORDER BY ts) AS positions,
array_agg(ts ORDER BY ts) AS times
FROM iss_tle,
generate_series(
'2024-01-01 12:00:00+00'::timestamptz,
'2024-01-01 13:30:00+00'::timestamptz,
'5 minutes'::interval
) AS ts
),
result AS (
SELECT (tle_from_eci(positions, times, t, fit_bstar := true)).* FROM obs, iss_tle
)
SELECT
covariance IS NOT NULL AS has_covariance,
array_length(covariance, 1) = 28 AS cov_length_28,
nstate = 7 AS is_7state
FROM result;
has_covariance | cov_length_28 | is_7state
----------------+---------------+-----------
t | t | t
(1 row)
-- ============================================================
-- Test 17: Condition number positive
--
-- The condition number (s_max / s_min from SVD) should be
-- positive for any non-degenerate fit. Topocentric mode.
-- ============================================================
WITH iss_tle AS (
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
),
mit AS (
SELECT '(42.36,-71.09,20)'::observer AS obs
),
topo_obs AS (
SELECT
array_agg(observe(t, obs, ts) ORDER BY ts) AS observations,
array_agg(ts ORDER BY ts) AS times
FROM iss_tle, mit,
generate_series(
'2024-01-01 12:00:00+00'::timestamptz,
'2024-01-01 13:30:00+00'::timestamptz,
'5 minutes'::interval
) AS ts
),
result AS (
SELECT (tle_from_topocentric(observations, times, obs, t)).* FROM topo_obs, mit, iss_tle
)
SELECT
condition_number > 0 AS cond_positive,
covariance IS NOT NULL AS has_covariance,
nstate = 6 AS is_6state
FROM result;
cond_positive | has_covariance | is_6state
---------------+----------------+-----------
t | t | t
(1 row)
Reference in New Issue View Git Blame Copy Permalink