-- 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;

-- ============================================================
-- 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;

-- ============================================================
-- 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;

-- ============================================================
-- 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;

-- ============================================================
-- 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;

-- ============================================================
-- 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 $$;

-- ============================================================
-- 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);

-- ============================================================
-- 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;

-- ============================================================
-- 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;

-- ============================================================
-- 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;

-- ============================================================
-- 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 $$;

-- ============================================================
-- 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;

-- ============================================================
-- 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;

-- ============================================================
-- 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;

-- ============================================================
-- 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;

-- ============================================================
-- 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;

-- ============================================================
-- 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;
