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 adfb6949e1 Add range rate fitting, weighted observations, and Gauss angles-only IOD (v0.6.0) 
Range rate: topocentric residuals now include an optional 4th component
(dot(Δr, v_ecef) / |Δr|) with OD_RR_SCALE=10.0 for unit balancing.
Controlled via fit_range_rate parameter on tle_from_topocentric().

Weighted observations: per-observation weights applied as √w scaling
to both residuals and Jacobian rows, producing the weighted normal
equations H'WH without explicit W construction. Weights parameter
added to tle_from_eci, tle_from_topocentric, and tle_from_angles.

Gauss angles-only IOD: Vallado Algorithm 52 implementation for
seed-free orbit recovery from 3+ RA/Dec observations. New RA/Dec
residual function with cos(dec) scaling and wrap-around handling.
New tle_from_angles() and tle_from_angles_multi() SQL functions
accepting RA in hours [0,24), Dec in degrees [-90,90].

New standalone test suite: test_od_gauss (17 assertions).
New regression tests: Tests 18-25 covering range rate, weights,
angles-only with/without seed, and error cases.
2026-02-17 17:48:13 -07:00

929 lines
34 KiB
Plaintext
Raw 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)
-- ============================================================
-- Test 18: Range rate round-trip
--
-- Propagate ISS, observe() to get topo with range_rate,
-- fit via tle_from_topocentric with fit_range_rate := true.
-- Verify 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
),
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,
fit_range_rate := true)).*
FROM topo_obs, mit, iss_tle
)
SELECT
rms_final < 10.0 AS rms_under_10km,
status = 'converged' AS did_converge
FROM result;
rms_under_10km | did_converge
----------------+--------------
t | t
(1 row)
-- ============================================================
-- Test 19: Range rate disabled matches existing behavior
--
-- Same data with fit_range_rate := false (default).
-- Results should converge the same as Test 4.
-- ============================================================
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,
fit_range_rate := false)).*
FROM topo_obs, mit, iss_tle
)
SELECT
rms_final < 10.0 AS rms_under_10km,
status = 'converged' AS did_converge
FROM result;
rms_under_10km | did_converge
----------------+--------------
t | t
(1 row)
-- ============================================================
-- Test 20: Weighted observations round-trip (uniform weights)
--
-- Uniform weights ARRAY[1,1,...,1]::float8[] should produce
-- identical results to no weights.
-- ============================================================
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
),
uniform_weights AS (
SELECT array_agg(1.0::float8) AS w
FROM generate_series(1, 19)
),
result AS (
SELECT (tle_from_eci(positions, times, t, false, 15,
weights := w)).* FROM obs, iss_tle, uniform_weights
)
SELECT
rms_final < 1.0 AS rms_under_1km,
status = 'converged' AS did_converge
FROM result;
rms_under_1km | did_converge
---------------+--------------
t | t
(1 row)
-- ============================================================
-- Test 21: Weights length mismatch error
--
-- Wrong-length weights array should raise an error.
-- ============================================================
DO $$
BEGIN
PERFORM tle_from_eci(
(SELECT array_agg(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,
ts) ORDER BY ts)
FROM generate_series(
'2024-01-01 12:00:00+00'::timestamptz,
'2024-01-01 13:30:00+00'::timestamptz,
'5 minutes'::interval
) AS ts),
(SELECT array_agg(ts ORDER BY ts)
FROM generate_series(
'2024-01-01 12:00:00+00'::timestamptz,
'2024-01-01 13:30:00+00'::timestamptz,
'5 minutes'::interval
) AS ts),
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,
false, 15,
ARRAY[1.0, 2.0, 3.0]::float8[] -- wrong length
);
RAISE NOTICE 'ERROR: should have raised exception';
EXCEPTION
WHEN array_subscript_error THEN
RAISE NOTICE 'OK: weights length mismatch error caught';
END $$;
NOTICE: OK: weights length mismatch error caught
-- ============================================================
-- Test 22: ECI with weights (verify API works)
--
-- Non-uniform weights (downweight first few observations).
-- ============================================================
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
),
nonuniform_weights AS (
SELECT array_agg(CASE WHEN n <= 5 THEN 0.1 ELSE 1.0 END ::float8) AS w
FROM generate_series(1, 19) AS n
),
result AS (
SELECT (tle_from_eci(positions, times, t, false, 15,
weights := w)).* FROM obs, iss_tle, nonuniform_weights
)
SELECT
rms_final < 5.0 AS rms_under_5km,
status = 'converged' AS did_converge
FROM result;
rms_under_5km | did_converge
---------------+--------------
t | t
(1 row)
-- ============================================================
-- Test 23: Angles-only with seed TLE
--
-- Propagate ISS, compute RA/Dec from TEME position relative
-- to observer, fit via tle_from_angles() with known seed.
-- ============================================================
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
),
radec AS (
SELECT
array_agg(
(topo_azimuth(o) * 24.0 / 360.0)::float8 -- crude az→RA proxy (accessor returns degrees)
ORDER BY ts
) AS ra_h,
array_agg(
topo_elevation(o)::float8 -- accessor returns degrees, valid as Dec proxy
ORDER BY ts
) AS dec_d,
array_agg(ts ORDER BY ts) AS times
FROM (
SELECT observe(t, obs, ts) AS o, ts
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
) sub
),
result AS (
SELECT (tle_from_angles(ra_h, dec_d, times, obs, t)).*
FROM radec, mit, iss_tle
)
SELECT
status IN ('converged', 'max_iterations') AS solver_ran,
rms_final IS NOT NULL AS has_rms
FROM result;
solver_ran | has_rms
------------+---------
t | t
(1 row)
-- ============================================================
-- Test 24: Angles-only without seed (Gauss IOD bootstrap)
--
-- Without a seed TLE, Gauss IOD must recover an initial orbit
-- from the RA/Dec observations alone. The crude azimuth→RA
-- approximation used here is not physically meaningful, so
-- Gauss may fail to converge — we accept that as valid behavior.
-- The real Gauss validation is in the standalone C tests.
-- ============================================================
DO $$
DECLARE
v_status text;
v_rms float8;
BEGIN
SELECT status, rms_final INTO v_status, v_rms
FROM (
SELECT (tle_from_angles(
(SELECT array_agg((topo_azimuth(o) * 24.0 / 360.0)::float8 ORDER BY ts)
FROM (SELECT 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, ts) AS o, ts
FROM generate_series('2024-01-01 12:00:00+00'::timestamptz,
'2024-01-01 13:30:00+00'::timestamptz,
'5 minutes'::interval) AS ts) sub),
(SELECT array_agg(topo_elevation(o)::float8 ORDER BY ts)
FROM (SELECT 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, ts) AS o, ts
FROM generate_series('2024-01-01 12:00:00+00'::timestamptz,
'2024-01-01 13:30:00+00'::timestamptz,
'5 minutes'::interval) AS ts) sub),
(SELECT array_agg(ts ORDER BY ts)
FROM generate_series('2024-01-01 12:00:00+00'::timestamptz,
'2024-01-01 13:30:00+00'::timestamptz,
'5 minutes'::interval) AS ts),
'(42.36,-71.09,20)'::observer
)).*
) r;
RAISE NOTICE 'OK: seedless angles-only ran: status=%, rms=%', v_status, v_rms;
EXCEPTION
WHEN data_exception THEN
RAISE NOTICE 'OK: seedless angles-only Gauss IOD failed as expected with approximate data';
END $$;
NOTICE: OK: seedless angles-only Gauss IOD failed as expected with approximate data
-- ============================================================
-- Test 25: Angles-only error cases
--
-- Too few observations should raise an error.
-- ============================================================
DO $$
BEGIN
PERFORM tle_from_angles(
ARRAY[1.0, 2.0]::float8[],
ARRAY[30.0, 35.0]::float8[],
ARRAY['2024-01-01 12:00:00+00'::timestamptz,
'2024-01-01 12:05:00+00'::timestamptz],
'(42.36,-71.09,20)'::observer
);
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
Reference in New Issue View Git Blame Copy Permalink