diff --git a/Dockerfile b/Dockerfile index 2798449..3b395a8 100644 --- a/Dockerfile +++ b/Dockerfile @@ -46,6 +46,7 @@ RUN su postgres -c "/usr/lib/postgresql/${PG_MAJOR}/bin/initdb -D /tmp/pgtest" & # Standalone unit tests (no PostgreSQL dependency) RUN make test-de-reader RUN make test-od-math +RUN make test-od-iod # Capture artifacts under /pg_orrery prefix for the next stage RUN make PG_CONFIG=${PG_CONFIG} DESTDIR=/pg_orrery install diff --git a/Makefile b/Makefile index 936bce8..8104e3f 100644 --- a/Makefile +++ b/Makefile @@ -3,7 +3,7 @@ EXTENSION = pg_orrery DATA = sql/pg_orrery--0.1.0.sql sql/pg_orrery--0.2.0.sql sql/pg_orrery--0.1.0--0.2.0.sql \ sql/pg_orrery--0.3.0.sql sql/pg_orrery--0.2.0--0.3.0.sql \ sql/pg_orrery--0.4.0.sql sql/pg_orrery--0.3.0--0.4.0.sql \ - sql/pg_orrery--0.4.0--0.5.0.sql + sql/pg_orrery--0.5.0.sql sql/pg_orrery--0.4.0--0.5.0.sql # Our extension C sources OBJS = src/pg_orrery.o src/tle_type.o src/eci_type.o src/observer_type.o \ diff --git a/TODO b/TODO index ada691e..6393dd3 100644 --- a/TODO +++ b/TODO @@ -1,10 +1,7 @@ - - IOD bootstrap: Gauss/Gibbs/double-r methods to generate - initial TLE from angles-only (eliminates seed requirement) - - Covariance output: Return (A^TWA)^{-1} matrix for - uncertainty estimates / conjunction screening - - Multi-observer: Accept observations from multiple ground - stations - - Adaptive step limiting: Dynamically tune correction limits - per Vallado's observation about different strategies failing - on different satellite subsets + - Gauss method for angles-only initial orbit determination + (eliminates seed requirement for sensors without ranging) + - Weighted observations (per-obs covariance weighting for + heterogeneous sensor fusion) + - Range rate fitting in topocentric mode (currently reserved + via vel_ecef in residual computation) diff --git a/pg_orrery.control b/pg_orrery.control index 89b48bc..677a08c 100644 --- a/pg_orrery.control +++ b/pg_orrery.control @@ -1,4 +1,4 @@ comment = 'A database orrery — celestial mechanics types and functions for PostgreSQL' -default_version = '0.4.0' +default_version = '0.5.0' module_pathname = '$libdir/pg_orrery' relocatable = true diff --git a/sql/pg_orrery--0.4.0--0.5.0.sql b/sql/pg_orrery--0.4.0--0.5.0.sql index 388fdcb..bc5764a 100644 --- a/sql/pg_orrery--0.4.0--0.5.0.sql +++ b/sql/pg_orrery--0.4.0--0.5.0.sql @@ -2,9 +2,49 @@ -- -- Adds multi-observer support, IOD bootstrap (seed-free fitting), -- and covariance output for uncertainty estimation. +-- +-- Covariance changes the return type of tle_from_eci and +-- tle_from_topocentric (5 → 8 OUT params), which requires +-- DROP + re-CREATE. -- ============================================================ --- Multi-observer topocentric fitting +-- Drop old 5-column OD functions +-- ============================================================ + +DROP FUNCTION IF EXISTS tle_from_eci(eci_position[], timestamptz[], tle, boolean, int4); +DROP FUNCTION IF EXISTS tle_from_topocentric(topocentric[], timestamptz[], observer, tle, boolean, int4); + +-- ============================================================ +-- Re-create with 8-column output (adds covariance) +-- ============================================================ + +CREATE FUNCTION tle_from_eci( + positions eci_position[], times timestamptz[], + seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false, + max_iter int4 DEFAULT 15, + 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 + AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE; +COMMENT ON FUNCTION tle_from_eci(eci_position[], timestamptz[], tle, boolean, int4) IS + 'Fit a TLE from ECI position/velocity observations via differential correction. Returns fitted TLE, iteration count, RMS residuals, convergence status, condition number, and formal covariance matrix.'; + +CREATE FUNCTION tle_from_topocentric( + observations topocentric[], times timestamptz[], + obs observer, + seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false, + max_iter int4 DEFAULT 15, + 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 + AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE; +COMMENT ON FUNCTION tle_from_topocentric(topocentric[], timestamptz[], observer, tle, boolean, int4) IS + 'Fit a TLE from topocentric (az/el/range) observations via differential correction. Returns fitted TLE, RMS residuals, convergence status, condition number, and formal covariance matrix.'; + +-- ============================================================ +-- Multi-observer topocentric fitting (new overload) -- ============================================================ CREATE FUNCTION tle_from_topocentric( @@ -13,9 +53,10 @@ CREATE FUNCTION tle_from_topocentric( seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false, max_iter int4 DEFAULT 15, OUT fitted_tle tle, OUT iterations int4, - OUT rms_final float8, OUT rms_initial float8, OUT status text + OUT rms_final float8, OUT rms_initial float8, OUT status text, + OUT condition_number float8, OUT covariance float8[], OUT nstate int4 ) RETURNS RECORD AS 'MODULE_PATHNAME', 'tle_from_topocentric_multi' LANGUAGE C STABLE PARALLEL SAFE; COMMENT ON FUNCTION tle_from_topocentric(topocentric[], timestamptz[], observer[], int4[], tle, boolean, int4) IS - 'Fit a TLE from topocentric observations collected by multiple ground stations. observer_ids[i] indexes into observers[]. Requires seed TLE and >= 6 observations.'; + 'Fit a TLE from topocentric observations collected by multiple ground stations. observer_ids[i] indexes into observers[]. Returns convergence status, condition number, and formal covariance matrix.'; diff --git a/sql/pg_orrery--0.5.0.sql b/sql/pg_orrery--0.5.0.sql new file mode 100644 index 0000000..9d56206 --- /dev/null +++ b/sql/pg_orrery--0.5.0.sql @@ -0,0 +1,862 @@ +-- pg_orrery -- Orbital mechanics types and functions for PostgreSQL +-- +-- Types: tle, eci_position, geodetic, topocentric, observer, pass_event +-- Provides SGP4/SDP4 propagation, coordinate transforms, pass prediction, +-- and GiST indexing on altitude bands for conjunction screening. +-- +-- All propagation uses WGS-72 constants (matching TLE mean element fitting). +-- Coordinate output uses WGS-84 (matching modern geodetic standards). + +-- ============================================================ +-- Shell types (forward declarations) +-- ============================================================ + +CREATE TYPE tle; +CREATE TYPE eci_position; +CREATE TYPE geodetic; +CREATE TYPE topocentric; +CREATE TYPE observer; +CREATE TYPE pass_event; + + +-- ============================================================ +-- TLE type: Two-Line Element set +-- ============================================================ + +CREATE FUNCTION tle_in(cstring) RETURNS tle + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION tle_out(tle) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION tle_recv(internal) RETURNS tle + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION tle_send(tle) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE tle ( + INPUT = tle_in, + OUTPUT = tle_out, + RECEIVE = tle_recv, + SEND = tle_send, + INTERNALLENGTH = 112, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE tle IS 'Two-Line Element set — parsed mean orbital elements for SGP4/SDP4 propagation'; + +-- TLE accessor functions + +CREATE FUNCTION tle_epoch(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_epoch(tle) IS 'TLE epoch as Julian date (UTC)'; + +CREATE FUNCTION tle_norad_id(tle) RETURNS int4 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_norad_id(tle) IS 'NORAD catalog number'; + +CREATE FUNCTION tle_inclination(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_inclination(tle) IS 'Orbital inclination in degrees'; + +CREATE FUNCTION tle_eccentricity(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_eccentricity(tle) IS 'Orbital eccentricity (dimensionless)'; + +CREATE FUNCTION tle_raan(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_raan(tle) IS 'Right ascension of ascending node in degrees'; + +CREATE FUNCTION tle_arg_perigee(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_arg_perigee(tle) IS 'Argument of perigee in degrees'; + +CREATE FUNCTION tle_mean_anomaly(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_mean_anomaly(tle) IS 'Mean anomaly in degrees'; + +CREATE FUNCTION tle_mean_motion(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_mean_motion(tle) IS 'Mean motion in revolutions per day'; + +CREATE FUNCTION tle_bstar(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_bstar(tle) IS 'B* drag coefficient (1/earth-radii)'; + +CREATE FUNCTION tle_period(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_period(tle) IS 'Orbital period in minutes'; + +CREATE FUNCTION tle_age(tle, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_age(tle, timestamptz) IS 'TLE age in days (positive = past epoch, negative = before epoch)'; + +CREATE FUNCTION tle_perigee(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_perigee(tle) IS 'Perigee altitude in km above WGS-72 ellipsoid'; + +CREATE FUNCTION tle_apogee(tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_apogee(tle) IS 'Apogee altitude in km above WGS-72 ellipsoid'; + +CREATE FUNCTION tle_intl_desig(tle) RETURNS text + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_intl_desig(tle) IS 'International designator (COSPAR ID)'; + +CREATE FUNCTION tle_from_lines(text, text) RETURNS tle + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_from_lines(text, text) IS + 'Construct TLE from separate line1/line2 text columns'; + + +-- ============================================================ +-- ECI position type: True Equator Mean Equinox (TEME) frame +-- ============================================================ + +CREATE FUNCTION eci_in(cstring) RETURNS eci_position + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_out(eci_position) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_recv(internal) RETURNS eci_position + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_send(eci_position) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE eci_position ( + INPUT = eci_in, + OUTPUT = eci_out, + RECEIVE = eci_recv, + SEND = eci_send, + INTERNALLENGTH = 48, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE eci_position IS 'Earth-Centered Inertial position and velocity in TEME frame (km, km/s)'; + +-- ECI accessor functions + +CREATE FUNCTION eci_x(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_y(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_z(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_vx(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_vy(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_vz(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION eci_speed(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eci_speed(eci_position) IS 'Velocity magnitude in km/s'; + +CREATE FUNCTION eci_altitude(eci_position) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eci_altitude(eci_position) IS 'Approximate geocentric altitude in km (radius - WGS72_AE)'; + + +-- ============================================================ +-- Geodetic type: WGS-84 latitude/longitude/altitude +-- ============================================================ + +CREATE FUNCTION geodetic_in(cstring) RETURNS geodetic + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION geodetic_out(geodetic) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION geodetic_recv(internal) RETURNS geodetic + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION geodetic_send(geodetic) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE geodetic ( + INPUT = geodetic_in, + OUTPUT = geodetic_out, + RECEIVE = geodetic_recv, + SEND = geodetic_send, + INTERNALLENGTH = 24, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE geodetic IS 'Geodetic coordinates on WGS-84 ellipsoid (lat/lon in degrees, altitude in km)'; + +CREATE FUNCTION geodetic_lat(geodetic) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION geodetic_lon(geodetic) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION geodetic_alt(geodetic) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + + +-- ============================================================ +-- Topocentric type: observer-relative az/el/range +-- ============================================================ + +CREATE FUNCTION topocentric_in(cstring) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION topocentric_out(topocentric) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION topocentric_recv(internal) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION topocentric_send(topocentric) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE topocentric ( + INPUT = topocentric_in, + OUTPUT = topocentric_out, + RECEIVE = topocentric_recv, + SEND = topocentric_send, + INTERNALLENGTH = 32, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE topocentric IS 'Topocentric coordinates relative to observer (azimuth, elevation, range, range rate)'; + +CREATE FUNCTION topo_azimuth(topocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION topo_azimuth(topocentric) IS 'Azimuth in degrees (0=N, 90=E, 180=S, 270=W)'; + +CREATE FUNCTION topo_elevation(topocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION topo_elevation(topocentric) IS 'Elevation in degrees (0=horizon, 90=zenith)'; + +CREATE FUNCTION topo_range(topocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION topo_range(topocentric) IS 'Slant range in km'; + +CREATE FUNCTION topo_range_rate(topocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION topo_range_rate(topocentric) IS 'Range rate in km/s (positive = receding)'; + + +-- ============================================================ +-- Observer type: ground station location +-- ============================================================ + +CREATE FUNCTION observer_in(cstring) RETURNS observer + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION observer_out(observer) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION observer_recv(internal) RETURNS observer + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION observer_send(observer) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE observer ( + INPUT = observer_in, + OUTPUT = observer_out, + RECEIVE = observer_recv, + SEND = observer_send, + INTERNALLENGTH = 24, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE observer IS 'Ground observer location (accepts: 40.0N 105.3W 1655m or decimal degrees)'; + +CREATE FUNCTION observer_lat(observer) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION observer_lat(observer) IS 'Latitude in degrees (positive = North)'; + +CREATE FUNCTION observer_lon(observer) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION observer_lon(observer) IS 'Longitude in degrees (positive = East)'; + +CREATE FUNCTION observer_alt(observer) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION observer_alt(observer) IS 'Altitude in meters above WGS-84 ellipsoid'; + +CREATE FUNCTION observer_from_geodetic(float8, float8, float8 DEFAULT 0.0) RETURNS observer + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION observer_from_geodetic(float8, float8, float8) IS + 'Construct observer from lat (deg), lon (deg), altitude (meters). Avoids text formatting round-trips.'; + + +-- ============================================================ +-- Pass event type: satellite visibility window +-- ============================================================ + +CREATE FUNCTION pass_event_in(cstring) RETURNS pass_event + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION pass_event_out(pass_event) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION pass_event_recv(internal) RETURNS pass_event + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION pass_event_send(pass_event) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE pass_event ( + INPUT = pass_event_in, + OUTPUT = pass_event_out, + RECEIVE = pass_event_recv, + SEND = pass_event_send, + INTERNALLENGTH = 48, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE pass_event IS 'Satellite pass event (AOS/MAX/LOS times, max elevation, AOS/LOS azimuths)'; + +CREATE FUNCTION pass_aos_time(pass_event) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_aos_time(pass_event) IS 'Acquisition of signal time'; + +CREATE FUNCTION pass_max_el_time(pass_event) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_max_el_time(pass_event) IS 'Maximum elevation time'; + +CREATE FUNCTION pass_los_time(pass_event) RETURNS timestamptz + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_los_time(pass_event) IS 'Loss of signal time'; + +CREATE FUNCTION pass_max_elevation(pass_event) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_max_elevation(pass_event) IS 'Maximum elevation in degrees'; + +CREATE FUNCTION pass_aos_azimuth(pass_event) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_aos_azimuth(pass_event) IS 'AOS azimuth in degrees (0=N)'; + +CREATE FUNCTION pass_los_azimuth(pass_event) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_los_azimuth(pass_event) IS 'LOS azimuth in degrees (0=N)'; + +CREATE FUNCTION pass_duration(pass_event) RETURNS interval + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_duration(pass_event) IS 'Pass duration (LOS - AOS)'; + + +-- ============================================================ +-- SGP4/SDP4 propagation functions +-- ============================================================ + +CREATE FUNCTION sgp4_propagate(tle, timestamptz) RETURNS eci_position + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sgp4_propagate(tle, timestamptz) IS + 'Propagate TLE to a point in time using SGP4 (near-earth) or SDP4 (deep-space). Returns TEME ECI position/velocity.'; + +CREATE FUNCTION sgp4_propagate_safe(tle, timestamptz) RETURNS eci_position + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE PARALLEL SAFE; +COMMENT ON FUNCTION sgp4_propagate_safe(tle, timestamptz) IS + 'Like sgp4_propagate but returns NULL on error instead of raising an exception. For batch queries with potentially invalid TLEs.'; + +CREATE FUNCTION sgp4_propagate_series(tle, timestamptz, timestamptz, interval) + RETURNS TABLE(t timestamptz, x float8, y float8, z float8, vx float8, vy float8, vz float8) + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE + ROWS 100; +COMMENT ON FUNCTION sgp4_propagate_series(tle, timestamptz, timestamptz, interval) IS + 'Propagate TLE over a time range at regular intervals. Returns time series of TEME positions.'; + +CREATE FUNCTION tle_distance(tle, tle, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_distance(tle, tle, timestamptz) IS + 'Euclidean distance in km between two TLEs at a reference time'; + + +-- ============================================================ +-- Coordinate transform functions +-- ============================================================ + +CREATE FUNCTION eci_to_geodetic(eci_position, timestamptz) RETURNS geodetic + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eci_to_geodetic(eci_position, timestamptz) IS + 'Convert TEME ECI position to WGS-84 geodetic coordinates at given time'; + +CREATE FUNCTION eci_to_topocentric(eci_position, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION eci_to_topocentric(eci_position, observer, timestamptz) IS + 'Convert TEME ECI position to topocentric (az/el/range) relative to observer'; + +CREATE FUNCTION subsatellite_point(tle, timestamptz) RETURNS geodetic + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION subsatellite_point(tle, timestamptz) IS + 'Subsatellite (nadir) point on WGS-84 ellipsoid at given time'; + +CREATE FUNCTION ground_track(tle, timestamptz, timestamptz, interval) + RETURNS TABLE(t timestamptz, lat float8, lon float8, alt float8) + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE + ROWS 100; +COMMENT ON FUNCTION ground_track(tle, timestamptz, timestamptz, interval) IS + 'Ground track as time series of subsatellite points (lat/lon in degrees, alt in km)'; + +CREATE FUNCTION observe(tle, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION observe(tle, observer, timestamptz) IS + 'Propagate TLE and compute observer-relative look angles in one call. Returns topocentric (az/el/range/range_rate).'; + +CREATE FUNCTION observe_safe(tle, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE PARALLEL SAFE; +COMMENT ON FUNCTION observe_safe(tle, observer, timestamptz) IS + 'Like observe() but returns NULL on propagation error. For batch queries with potentially invalid/decayed TLEs.'; + + +-- ============================================================ +-- Pass prediction functions +-- ============================================================ + +CREATE FUNCTION next_pass(tle, observer, timestamptz) RETURNS pass_event + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION next_pass(tle, observer, timestamptz) IS + 'Find the next satellite pass over observer (searches up to 7 days ahead)'; + +CREATE FUNCTION predict_passes(tle, observer, timestamptz, timestamptz, float8 DEFAULT 0.0) + RETURNS SETOF pass_event + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE + ROWS 10; +COMMENT ON FUNCTION predict_passes(tle, observer, timestamptz, timestamptz, float8) IS + 'Predict all satellite passes over observer in time window. Optional min_elevation in degrees.'; + +CREATE FUNCTION pass_visible(tle, observer, timestamptz, timestamptz) RETURNS boolean + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION pass_visible(tle, observer, timestamptz, timestamptz) IS + 'True if any pass occurs over observer in the time window'; + + +-- ============================================================ +-- GiST operator support functions +-- ============================================================ + +-- Overlap operator: do orbital keys overlap in altitude AND inclination? +CREATE FUNCTION tle_overlap(tle, tle) RETURNS boolean + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +-- Altitude distance operator (altitude-only, for KNN ordering) +CREATE FUNCTION tle_alt_distance(tle, tle) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OPERATOR && ( + LEFTARG = tle, + RIGHTARG = tle, + FUNCTION = tle_overlap, + COMMUTATOR = &&, + RESTRICT = areasel, + JOIN = areajoinsel +); + +COMMENT ON OPERATOR && (tle, tle) IS 'Orbital key overlap (altitude band AND inclination range) — necessary condition for conjunction'; + +CREATE OPERATOR <-> ( + LEFTARG = tle, + RIGHTARG = tle, + FUNCTION = tle_alt_distance, + COMMUTATOR = <-> +); + +COMMENT ON OPERATOR <-> (tle, tle) IS 'Minimum altitude-band separation in km (0 if overlapping). Altitude-only — does not account for inclination. Use && for 2-D filtering.'; + + +-- ============================================================ +-- GiST operator class for 2-D orbital indexing (altitude + inclination) +-- ============================================================ + +-- GiST internal support functions +CREATE FUNCTION gist_tle_compress(internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_tle_decompress(internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_tle_consistent(internal, tle, smallint, oid, internal) RETURNS boolean + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_tle_union(internal, internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_tle_penalty(internal, internal, internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_tle_picksplit(internal, internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_tle_same(internal, internal, internal) RETURNS internal + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION gist_tle_distance(internal, tle, smallint, oid, internal) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE OPERATOR CLASS tle_ops + DEFAULT FOR TYPE tle USING gist AS + OPERATOR 3 && , + OPERATOR 15 <-> (tle, tle) FOR ORDER BY float_ops, + FUNCTION 1 gist_tle_consistent(internal, tle, smallint, oid, internal), + FUNCTION 2 gist_tle_union(internal, internal), + FUNCTION 3 gist_tle_compress(internal), + FUNCTION 4 gist_tle_decompress(internal), + FUNCTION 5 gist_tle_penalty(internal, internal, internal), + FUNCTION 6 gist_tle_picksplit(internal, internal), + FUNCTION 7 gist_tle_same(internal, internal, internal), + FUNCTION 8 gist_tle_distance(internal, tle, smallint, oid, internal); +-- pg_orrery 0.1.0 -> 0.2.0 migration +-- +-- Phase 1: Stars, comets, and Keplerian propagation. +-- Adds heliocentric type, star observation, and two-body propagation. + + +-- ============================================================ +-- Heliocentric type: ecliptic J2000 position in AU +-- ============================================================ + +CREATE TYPE heliocentric; + +CREATE FUNCTION heliocentric_in(cstring) RETURNS heliocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION heliocentric_out(heliocentric) RETURNS cstring + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION heliocentric_recv(internal) RETURNS heliocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE FUNCTION heliocentric_send(heliocentric) RETURNS bytea + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; + +CREATE TYPE heliocentric ( + INPUT = heliocentric_in, + OUTPUT = heliocentric_out, + RECEIVE = heliocentric_recv, + SEND = heliocentric_send, + INTERNALLENGTH = 24, + ALIGNMENT = double, + STORAGE = plain +); + +COMMENT ON TYPE heliocentric IS 'Heliocentric position in ecliptic J2000 frame (x, y, z in AU)'; + +CREATE FUNCTION helio_x(heliocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION helio_x(heliocentric) IS 'X component in AU (ecliptic J2000, vernal equinox direction)'; + +CREATE FUNCTION helio_y(heliocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION helio_y(heliocentric) IS 'Y component in AU (ecliptic J2000)'; + +CREATE FUNCTION helio_z(heliocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION helio_z(heliocentric) IS 'Z component in AU (ecliptic J2000, north ecliptic pole)'; + +CREATE FUNCTION helio_distance(heliocentric) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION helio_distance(heliocentric) IS 'Heliocentric distance in AU'; + + +-- ============================================================ +-- Star observation functions +-- ============================================================ + +CREATE FUNCTION star_observe(float8, float8, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION star_observe(float8, float8, observer, timestamptz) IS + 'Observe a star from (ra_hours J2000, dec_degrees J2000, observer, time). Returns topocentric az/el. Range is 0 (infinite distance).'; + +CREATE FUNCTION star_observe_safe(float8, float8, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE PARALLEL SAFE; +COMMENT ON FUNCTION star_observe_safe(float8, float8, observer, timestamptz) IS + 'Like star_observe but returns NULL on invalid inputs. For batch queries over star catalogs.'; + + +-- ============================================================ +-- Keplerian propagation functions +-- ============================================================ + +CREATE FUNCTION kepler_propagate( + q_au float8, eccentricity float8, + inclination_deg float8, arg_perihelion_deg float8, + long_asc_node_deg float8, perihelion_jd float8, + t timestamptz +) RETURNS heliocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION kepler_propagate(float8, float8, float8, float8, float8, float8, timestamptz) IS + 'Two-body Keplerian propagation from classical orbital elements. Returns heliocentric ecliptic J2000 position in AU. Handles elliptic, parabolic, and hyperbolic orbits.'; + + +-- ============================================================ +-- Comet observation +-- ============================================================ + +CREATE FUNCTION comet_observe( + q_au float8, eccentricity float8, + inclination_deg float8, arg_perihelion_deg float8, + long_asc_node_deg float8, perihelion_jd float8, + earth_x_au float8, earth_y_au float8, earth_z_au float8, + obs observer, t timestamptz +) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION comet_observe(float8, float8, float8, float8, float8, float8, float8, float8, float8, observer, timestamptz) IS + 'Observe a comet/asteroid from orbital elements. Requires Earth heliocentric ecliptic J2000 position (AU). Returns topocentric az/el with geocentric range in km.'; + + +-- ============================================================ +-- Phase 2: VSOP87 planets, ELP82B Moon, Sun observation +-- ============================================================ + +CREATE FUNCTION planet_heliocentric(int4, timestamptz) RETURNS heliocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_heliocentric(int4, timestamptz) IS + 'VSOP87 heliocentric ecliptic J2000 position (AU). Body IDs: 0=Sun, 1=Mercury, ..., 8=Neptune.'; + +CREATE FUNCTION planet_observe(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_observe(int4, observer, timestamptz) IS + 'Observe a planet from (body_id 1-8, observer, time). Returns topocentric az/el with geocentric range in km.'; + +CREATE FUNCTION sun_observe(observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_observe(observer, timestamptz) IS + 'Observe the Sun from (observer, time). Returns topocentric az/el with geocentric range in km.'; + +CREATE FUNCTION moon_observe(observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_observe(observer, timestamptz) IS + 'Observe the Moon via ELP2000-82B from (observer, time). Returns topocentric az/el with geocentric range in km.'; + + +-- ============================================================ +-- Phase 3: Planetary moon observation +-- ============================================================ + +CREATE FUNCTION galilean_observe(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION galilean_observe(int4, observer, timestamptz) IS + 'Observe a Galilean moon of Jupiter via L1.2 theory. Body IDs: 0=Io, 1=Europa, 2=Ganymede, 3=Callisto.'; + +CREATE FUNCTION saturn_moon_observe(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION saturn_moon_observe(int4, observer, timestamptz) IS + 'Observe a Saturn moon via TASS 1.7. Body IDs: 0=Mimas, 1=Enceladus, 2=Tethys, 3=Dione, 4=Rhea, 5=Titan, 6=Iapetus, 7=Hyperion.'; + +CREATE FUNCTION uranus_moon_observe(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION uranus_moon_observe(int4, observer, timestamptz) IS + 'Observe a Uranus moon via GUST86. Body IDs: 0=Miranda, 1=Ariel, 2=Umbriel, 3=Titania, 4=Oberon.'; + +CREATE FUNCTION mars_moon_observe(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION mars_moon_observe(int4, observer, timestamptz) IS + 'Observe a Mars moon via MarsSat. Body IDs: 0=Phobos, 1=Deimos.'; + + +-- ============================================================ +-- Phase 3: Jupiter decametric radio burst prediction +-- ============================================================ + +CREATE FUNCTION io_phase_angle(timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION io_phase_angle(timestamptz) IS + 'Io orbital phase angle in degrees [0,360). 0=superior conjunction (behind Jupiter). Standard Radio JOVE convention.'; + +CREATE FUNCTION jupiter_cml(observer, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION jupiter_cml(observer, timestamptz) IS + 'Jupiter Central Meridian Longitude, System III (1965.0), in degrees [0,360). Light-time corrected.'; + +CREATE FUNCTION jupiter_burst_probability(float8, float8) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION jupiter_burst_probability(float8, float8) IS + 'Estimated Jupiter decametric burst probability (0-1) from (io_phase_deg, cml_deg). Based on Carr et al. (1983) source regions A, B, C, D.'; + + +-- ============================================================ +-- Phase 4: Interplanetary transfer orbits (Lambert solver) +-- ============================================================ + +CREATE FUNCTION lambert_transfer( + dep_body_id int4, arr_body_id int4, + dep_time timestamptz, arr_time timestamptz, + OUT c3_departure float8, OUT c3_arrival float8, + OUT v_inf_departure float8, OUT v_inf_arrival float8, + OUT tof_days float8, OUT transfer_sma float8 +) RETURNS RECORD + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lambert_transfer(int4, int4, timestamptz, timestamptz) IS + 'Solve Lambert transfer between two planets. Returns C3 (km^2/s^2), v_infinity (km/s), TOF (days), SMA (AU). Body IDs 1-8.'; + +CREATE FUNCTION lambert_c3(int4, int4, timestamptz, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lambert_c3(int4, int4, timestamptz, timestamptz) IS + 'Departure C3 (km^2/s^2) for a Lambert transfer. Returns NULL if solver fails. For pork chop plots.'; +-- pg_orrery 0.2.0 -> 0.3.0 migration +-- +-- Adds optional JPL DE440/441 ephemeris functions. +-- Existing VSOP87/ELP2000-82B functions are unchanged (still IMMUTABLE). +-- New _de() functions are STABLE (depend on external DE binary file). +-- When DE is unavailable, _de() functions fall back to VSOP87 silently. + +-- ============================================================ +-- Phase 5: DE ephemeris functions (optional high-precision) +-- ============================================================ + +-- Planet observation with DE ephemeris + +CREATE FUNCTION planet_heliocentric_de(int4, timestamptz) RETURNS heliocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_heliocentric_de(int4, timestamptz) IS + 'Heliocentric ecliptic J2000 position via JPL DE (sub-arcsecond). Falls back to VSOP87 if DE unavailable. Body IDs: 0=Sun, 1-8=Mercury-Neptune.'; + +CREATE FUNCTION planet_observe_de(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION planet_observe_de(int4, observer, timestamptz) IS + 'Observe planet via JPL DE. Falls back to VSOP87. Body IDs: 1-8 (Mercury-Neptune).'; + +CREATE FUNCTION sun_observe_de(observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION sun_observe_de(observer, timestamptz) IS + 'Observe Sun via JPL DE. Falls back to VSOP87.'; + +CREATE FUNCTION moon_observe_de(observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION moon_observe_de(observer, timestamptz) IS + 'Observe Moon via JPL DE. Falls back to ELP2000-82B.'; + + +-- Lambert transfer with DE positions + +CREATE FUNCTION lambert_transfer_de( + dep_body_id int4, arr_body_id int4, + dep_time timestamptz, arr_time timestamptz, + OUT c3_departure float8, OUT c3_arrival float8, + OUT v_inf_departure float8, OUT v_inf_arrival float8, + OUT tof_days float8, OUT transfer_sma float8 +) RETURNS RECORD + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lambert_transfer_de(int4, int4, timestamptz, timestamptz) IS + 'Lambert transfer via JPL DE positions. Falls back to VSOP87. Body IDs 1-8.'; + +CREATE FUNCTION lambert_c3_de(int4, int4, timestamptz, timestamptz) RETURNS float8 + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION lambert_c3_de(int4, int4, timestamptz, timestamptz) IS + 'Departure C3 via JPL DE. Falls back to VSOP87. For pork chop plots.'; + + +-- Planetary moon observation with DE parent positions + +CREATE FUNCTION galilean_observe_de(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION galilean_observe_de(int4, observer, timestamptz) IS + 'Observe Galilean moon with JPL DE parent position. L1.2 moon offsets. Body IDs: 0-3 (Io-Callisto).'; + +CREATE FUNCTION saturn_moon_observe_de(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION saturn_moon_observe_de(int4, observer, timestamptz) IS + 'Observe Saturn moon with JPL DE parent position. TASS 1.7 moon offsets. Body IDs: 0-7 (Mimas-Hyperion).'; + +CREATE FUNCTION uranus_moon_observe_de(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION uranus_moon_observe_de(int4, observer, timestamptz) IS + 'Observe Uranus moon with JPL DE parent position. GUST86 moon offsets. Body IDs: 0-4 (Miranda-Oberon).'; + +CREATE FUNCTION mars_moon_observe_de(int4, observer, timestamptz) RETURNS topocentric + AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION mars_moon_observe_de(int4, observer, timestamptz) IS + 'Observe Mars moon with JPL DE parent position. MarsSat moon offsets. Body IDs: 0-1 (Phobos-Deimos).'; + + +-- Diagnostic function + +CREATE FUNCTION pg_orrery_ephemeris_info( + OUT provider text, OUT file_path text, + OUT start_jd float8, OUT end_jd float8, + OUT version int4, OUT au_km float8 +) RETURNS RECORD + AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE; +COMMENT ON FUNCTION pg_orrery_ephemeris_info() IS + 'Returns current ephemeris provider status: VSOP87 or JPL_DE with file path, JD range, version, and AU value.'; +-- pg_orrery 0.3.0 -> 0.4.0 migration +-- +-- Adds observation-to-TLE fitting via batch weighted least-squares +-- differential correction (Vallado & Crawford 2008, AIAA 2008-6770). +-- Uses equinoctial elements internally for singularity-free optimization. +-- LAPACK dgelss_() for SVD solve. + +-- ============================================================ +-- Phase 6: Orbit determination (TLE fitting from observations) +-- ============================================================ + +-- Fit TLE from ECI position/velocity ephemeris + +CREATE FUNCTION tle_from_eci( + positions eci_position[], times timestamptz[], + seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false, + max_iter int4 DEFAULT 15, + 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 + AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE; +COMMENT ON FUNCTION tle_from_eci(eci_position[], timestamptz[], tle, boolean, int4) IS + 'Fit a TLE from ECI position/velocity observations via differential correction. Returns fitted TLE, iteration count, RMS residuals, convergence status, condition number, and formal covariance matrix.'; + +-- Fit TLE from topocentric observations (az/el/range) — single observer + +CREATE FUNCTION tle_from_topocentric( + observations topocentric[], times timestamptz[], + obs observer, + seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false, + max_iter int4 DEFAULT 15, + 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 + AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE; +COMMENT ON FUNCTION tle_from_topocentric(topocentric[], timestamptz[], observer, tle, boolean, int4) IS + 'Fit a TLE from topocentric (az/el/range) observations via differential correction. Returns fitted TLE, RMS residuals, convergence status, condition number, and formal covariance matrix.'; + +-- Fit TLE from topocentric observations — multiple observers + +CREATE FUNCTION 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, + 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 + AS 'MODULE_PATHNAME', 'tle_from_topocentric_multi' + LANGUAGE C STABLE PARALLEL SAFE; +COMMENT ON FUNCTION tle_from_topocentric(topocentric[], timestamptz[], observer[], int4[], tle, boolean, int4) IS + 'Fit a TLE from topocentric observations collected by multiple ground stations. observer_ids[i] indexes into observers[]. Returns convergence status, condition number, and formal covariance matrix.'; + +-- Per-observation residuals diagnostic + +CREATE FUNCTION 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 +) + AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; +COMMENT ON FUNCTION tle_fit_residuals(tle, eci_position[], timestamptz[]) IS + 'Compute per-observation position residuals (km) between a TLE and ECI observations. Useful for fit quality assessment.'; diff --git a/src/od_funcs.c b/src/od_funcs.c index 2403ed2..378e7f0 100644 --- a/src/od_funcs.c +++ b/src/od_funcs.c @@ -24,6 +24,9 @@ #include #include +/* For construct_array() — covariance output */ +#include "utils/lsyscache.h" + PG_FUNCTION_INFO_V1(tle_from_eci); PG_FUNCTION_INFO_V1(tle_from_topocentric); PG_FUNCTION_INFO_V1(tle_from_topocentric_multi); @@ -113,8 +116,8 @@ tle_from_eci(PG_FUNCTION_ARGS) int i, rc; TupleDesc tupdesc; - Datum values[5]; - bool nulls[5]; + Datum values[8]; + bool nulls[8]; HeapTuple tuple; /* Deconstruct arrays. @@ -197,6 +200,28 @@ tle_from_eci(PG_FUNCTION_ARGS) values[2] = Float8GetDatum(result.rms_final); values[3] = Float8GetDatum(result.rms_initial); values[4] = CStringGetTextDatum(result.status); + values[5] = Float8GetDatum(result.condition_number); + + /* Covariance: float8[] or NULL */ + if (result.cov_size > 0) + { + Datum *cov_datums = (Datum *) palloc(sizeof(Datum) * result.cov_size); + int ci; + for (ci = 0; ci < result.cov_size; ci++) + cov_datums[ci] = Float8GetDatum(result.covariance[ci]); + values[6] = PointerGetDatum( + construct_array(cov_datums, result.cov_size, + FLOAT8OID, sizeof(float8), FLOAT8PASSBYVAL, + TYPALIGN_DOUBLE)); + } + else + { + nulls[6] = true; + } + + values[7] = Int32GetDatum(result.cov_size > 0 + ? (result.cov_size == 28 ? 7 : 6) + : 0); tuple = heap_form_tuple(tupdesc, values, nulls); PG_RETURN_DATUM(HeapTupleGetDatum(tuple)); @@ -229,8 +254,8 @@ tle_from_topocentric(PG_FUNCTION_ARGS) int i, rc; TupleDesc tupdesc; - Datum values[5]; - bool nulls[5]; + Datum values[8]; + bool nulls[8]; HeapTuple tuple; /* Build observer */ @@ -321,6 +346,27 @@ tle_from_topocentric(PG_FUNCTION_ARGS) values[2] = Float8GetDatum(result.rms_final); values[3] = Float8GetDatum(result.rms_initial); values[4] = CStringGetTextDatum(result.status); + values[5] = Float8GetDatum(result.condition_number); + + if (result.cov_size > 0) + { + Datum *cov_datums = (Datum *) palloc(sizeof(Datum) * result.cov_size); + int ci; + for (ci = 0; ci < result.cov_size; ci++) + cov_datums[ci] = Float8GetDatum(result.covariance[ci]); + values[6] = PointerGetDatum( + construct_array(cov_datums, result.cov_size, + FLOAT8OID, sizeof(float8), FLOAT8PASSBYVAL, + TYPALIGN_DOUBLE)); + } + else + { + nulls[6] = true; + } + + values[7] = Int32GetDatum(result.cov_size > 0 + ? (result.cov_size == 28 ? 7 : 6) + : 0); tuple = heap_form_tuple(tupdesc, values, nulls); PG_RETURN_DATUM(HeapTupleGetDatum(tuple)); @@ -358,8 +404,8 @@ tle_from_topocentric_multi(PG_FUNCTION_ARGS) int i, rc; TupleDesc tupdesc; - Datum values[5]; - bool nulls[5]; + Datum values[8]; + bool nulls[8]; HeapTuple tuple; /* Deconstruct all arrays */ @@ -473,6 +519,27 @@ tle_from_topocentric_multi(PG_FUNCTION_ARGS) values[2] = Float8GetDatum(result.rms_final); values[3] = Float8GetDatum(result.rms_initial); values[4] = CStringGetTextDatum(result.status); + values[5] = Float8GetDatum(result.condition_number); + + if (result.cov_size > 0) + { + Datum *cov_datums = (Datum *) palloc(sizeof(Datum) * result.cov_size); + int ci; + for (ci = 0; ci < result.cov_size; ci++) + cov_datums[ci] = Float8GetDatum(result.covariance[ci]); + values[6] = PointerGetDatum( + construct_array(cov_datums, result.cov_size, + FLOAT8OID, sizeof(float8), FLOAT8PASSBYVAL, + TYPALIGN_DOUBLE)); + } + else + { + nulls[6] = true; + } + + values[7] = Int32GetDatum(result.cov_size > 0 + ? (result.cov_size == 28 ? 7 : 6) + : 0); tuple = heap_form_tuple(tupdesc, values, nulls); PG_RETURN_DATUM(HeapTupleGetDatum(tuple)); diff --git a/src/od_solver.c b/src/od_solver.c index 49b6381..f26bc96 100644 --- a/src/od_solver.c +++ b/src/od_solver.c @@ -57,6 +57,14 @@ extern void dgelss_(int *m, int *n, int *nrhs, double *s, double *rcond, int *rank, double *work, int *lwork, int *info); +/* + * dpotrf_: Cholesky factorization of a symmetric positive-definite matrix. + * dpotri_: Invert using the Cholesky factor from dpotrf_. + * Used together to compute (H^T H)^{-1} for covariance output. + */ +extern void dpotrf_(char *uplo, int *n, double *a, int *lda, int *info); +extern void dpotri_(char *uplo, int *n, double *a, int *lda, int *info); + /* ── Internal propagation wrapper ─────────────────────── */ @@ -565,6 +573,7 @@ od_fit_tle(const od_observation_t *obs, int n_obs, int iter; int max_iter; int converged = 0; + int j_cov, k_cov; /* covariance loop indices */ /* Validate inputs */ nstate = config->fit_bstar ? OD_NSTATE_7 : OD_NSTATE_6; @@ -652,22 +661,14 @@ od_fit_tle(const od_observation_t *obs, int n_obs, result->rms_initial = rms_cur; rms_prev = rms_cur; - /* Already converged (perfect seed)? Skip iteration. */ + /* Already converged (perfect seed)? Skip DC loop but still + * compute covariance — users need uncertainty estimates even + * when the initial guess is exact. */ if (rms_cur < OD_RMS_ABS_TOL) { converged = 1; - memcpy(&result->fitted_tle, ¤t_tle, sizeof(tle_t)); - result->iterations = 0; - result->rms_final = rms_cur; - result->converged = 1; - snprintf(result->status, sizeof(result->status), "converged"); - - od_free(residuals); - od_free(jacobian); - od_free(dx); - od_free(sv); - od_free(work); - return 0; + iter = 0; + goto compute_covariance; } /* ── DC iteration loop ─────────────────────────────── */ @@ -848,6 +849,146 @@ od_fit_tle(const od_observation_t *obs, int n_obs, rms_prev = rms_cur; } +compute_covariance: + /* ── Covariance computation ─────────────────────────── */ + + /* + * Compute formal covariance matrix (H^T H)^{-1} at the converged + * solution. This requires rebuilding the Jacobian at the final state + * (the DC loop overwrites it each iteration via LAPACK). + * + * The covariance gives parameter uncertainty estimates — diagonal + * elements are variances, off-diagonals are correlations. Essential + * for conjunction screening probability computations. + */ + result->cov_size = 0; + result->condition_number = condition_number; + + if (converged || iter >= max_iter) + { + /* Rebuild Jacobian at final state (same finite-diff as DC loop) */ + for (j_cov = 0; j_cov < nstate; j_cov++) + { + od_equinoctial_t eq_pert; + tle_t tle_pert; + double *resid_pert; + double h, elem_val; + + memcpy(&eq_pert, &eq, sizeof(od_equinoctial_t)); + + switch (j_cov) + { + case 0: elem_val = eq.n; break; + case 1: elem_val = eq.af; break; + case 2: elem_val = eq.ag; break; + case 3: elem_val = eq.chi; break; + case 4: elem_val = eq.psi; break; + case 5: elem_val = eq.L0; break; + case 6: elem_val = eq.bstar; break; + default: elem_val = 0.0; break; + } + + h = fabs(elem_val) * 1e-6; + if (h < 1e-10) + h = 1e-10; + + switch (j_cov) + { + case 0: eq_pert.n += h; break; + case 1: eq_pert.af += h; break; + case 2: eq_pert.ag += h; break; + case 3: eq_pert.chi += h; break; + case 4: eq_pert.psi += h; break; + case 5: eq_pert.L0 += h; break; + case 6: eq_pert.bstar += h; break; + } + + equinoctial_to_tle(&eq_pert, &seed_tle, &tle_pert); + + resid_pert = (double *)od_alloc(sizeof(double) * nrows); + + if (config->obs_type == OD_OBS_ECI) + compute_residuals_eci(&tle_pert, obs, n_obs, resid_pert); + else + compute_residuals_topo(&tle_pert, obs, n_obs, + config->observers, resid_pert); + + for (k_cov = 0; k_cov < nrows; k_cov++) + jacobian[j_cov * nrows + k_cov] = + (residuals[k_cov] - resid_pert[k_cov]) / h; + + od_free(resid_pert); + } + + /* Compute H^T H before SVD destroys the Jacobian */ + { + double *hth; + int ci, cj, ri, info_cov; + char uplo = 'U'; + + hth = (double *)od_alloc(sizeof(double) * nstate * nstate); + memset(hth, 0, sizeof(double) * nstate * nstate); + + for (ci = 0; ci < nstate; ci++) + { + for (cj = ci; cj < nstate; cj++) + { + double dot = 0.0; + for (ri = 0; ri < nrows; ri++) + dot += jacobian[ci * nrows + ri] * + jacobian[cj * nrows + ri]; + hth[ci + cj * nstate] = dot; + hth[cj + ci * nstate] = dot; + } + } + + /* SVD of Jacobian for condition number (destroys jacobian) */ + { + int m = nrows, n = nstate, nrhs = 1; + int lda = nrows, ldb = nrows, rank, info_svd; + double rcond = -1.0; + + memcpy(dx, residuals, sizeof(double) * nrows); + + od_free(work); + { + int lw = -1; + double wq; + dgelss_(&m, &n, &nrhs, jacobian, &lda, dx, &ldb, + sv, &rcond, &rank, &wq, &lw, &info_svd); + lw = (int)wq + 1; + work = (double *)od_alloc(sizeof(double) * lw); + dgelss_(&m, &n, &nrhs, jacobian, &lda, dx, &ldb, + sv, &rcond, &rank, work, &lw, &info_svd); + } + + if (info_svd == 0 && sv[nstate - 1] > 1e-30) + condition_number = sv[0] / sv[nstate - 1]; + } + + result->condition_number = condition_number; + + /* Cholesky factorize H^T H, then invert for covariance */ + dpotrf_(&uplo, &nstate, hth, &nstate, &info_cov); + + if (info_cov == 0) + { + dpotri_(&uplo, &nstate, hth, &nstate, &info_cov); + + if (info_cov == 0) + { + int idx = 0; + for (ci = 0; ci < nstate; ci++) + for (cj = ci; cj < nstate; cj++) + result->covariance[idx++] = hth[ci + cj * nstate]; + result->cov_size = idx; /* 21 or 28 */ + } + } + + od_free(hth); + } + } + /* ── Build result ──────────────────────────────────── */ memcpy(&result->fitted_tle, ¤t_tle, sizeof(tle_t)); @@ -860,8 +1001,6 @@ od_fit_tle(const od_observation_t *obs, int n_obs, else if (iter >= max_iter && result->status[0] == '\0') snprintf(result->status, sizeof(result->status), "max_iterations"); - (void)condition_number; /* used by covariance output (v0.5.0) */ - /* Cleanup */ od_free(residuals); od_free(jacobian); diff --git a/src/od_solver.h b/src/od_solver.h index c56b03c..3ab9d96 100644 --- a/src/od_solver.h +++ b/src/od_solver.h @@ -79,6 +79,9 @@ typedef struct double rms_final; /* RMS residual after fitting (km) */ int converged; /* 1 if converged, 0 if hit max_iter */ char status[64]; /* human-readable status */ + double covariance[28]; /* upper triangle of (H^T H)^{-1}, row-major */ + int cov_size; /* 21 (6-state) or 28 (7-state), 0 if failed */ + double condition_number; /* s_max / s_min from final SVD */ } od_result_t; /* diff --git a/test/expected/od_fit.out b/test/expected/od_fit.out index 9d0ac61..39e8454 100644 --- a/test/expected/od_fit.out +++ b/test/expected/od_fit.out @@ -2,8 +2,9 @@ -- -- 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; -ALTER EXTENSION pg_orrery UPDATE TO '0.5.0'; +SET client_min_messages TO notice; -- ============================================================ -- Test 1: ECI round-trip (ISS-like LEO orbit) -- @@ -482,3 +483,143 @@ FROM result; 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) + diff --git a/test/sql/od_fit.sql b/test/sql/od_fit.sql index d6f5261..0b6a9e5 100644 --- a/test/sql/od_fit.sql +++ b/test/sql/od_fit.sql @@ -3,8 +3,9 @@ -- 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; -ALTER EXTENSION pg_orrery UPDATE TO '0.5.0'; +SET client_min_messages TO notice; -- ============================================================ -- Test 1: ECI round-trip (ISS-like LEO orbit) @@ -452,3 +453,131 @@ SELECT 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;