-- 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.
CREATE EXTENSION IF NOT EXISTS pg_orrery;
ALTER EXTENSION pg_orrery UPDATE TO '0.5.0';
-- ============================================================
-- 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