warehack.ing/pg_orrery
2
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

pg_orbit 0.2.0: Full solar system computation at the SQL layer

Browse Source 
Phase 1 — Stars, comets, Keplerian propagation:
- star_observe() / star_observe_safe(): fixed star alt/az via IAU 1976
  precession, equatorial-to-horizontal transform
- kepler_propagate(): two-body Keplerian orbit propagation for
  elliptic, parabolic, and hyperbolic orbits
- comet_observe(): observe comets/asteroids from orbital elements
- heliocentric type: ecliptic J2000 position (x, y, z in AU)

Phase 2 — VSOP87 planets, ELP82B Moon, Sun:
- planet_heliocentric(): VSOP87 heliocentric ecliptic J2000 positions
  for Mercury through Neptune (Bretagnon & Francou, MIT)
- planet_observe(): full observation pipeline for any planet
- sun_observe(): Sun position from negated Earth VSOP87
- moon_observe(): ELP2000-82B lunar position (Chapront-Touzé, MIT)
- Clean-room precession (IAU 2006) and sidereal time (IERS 2010)
- elliptic_to_rectangular utility (Stellarium, MIT)

All Stellarium extractions are MIT-licensed, thread-safe (static
caching removed for PARALLEL SAFE), zero external data files.

All 9 regression tests pass (90ms total).
This commit is contained in:
Ryan Malloy 2026-02-16 01:36:27 -07:00
parent 07bc4e47c6
commit 0544a78276
25 changed files with 303415 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
Makefile
View File

@ -1,10 +1,13 @@
MODULE_big = pg_orbit MODULE_big = pg_orbit
EXTENSION = pg_orbit EXTENSION = pg_orbit
DATA = sql/pg_orbit--0.1.0.sql DATA = sql/pg_orbit--0.1.0.sql sql/pg_orbit--0.2.0.sql sql/pg_orbit--0.1.0--0.2.0.sql
# Our extension C sources # Our extension C sources
OBJS = src/pg_orbit.o src/tle_type.o src/eci_type.o src/observer_type.o \ OBJS = src/pg_orbit.o src/tle_type.o src/eci_type.o src/observer_type.o \
src/sgp4_funcs.o src/coord_funcs.o src/pass_funcs.o src/gist_tle.o src/sgp4_funcs.o src/coord_funcs.o src/pass_funcs.o src/gist_tle.o \
src/star_funcs.o src/kepler_funcs.o \
src/vsop87.o src/elp82b.o src/elliptic_to_rectangular.o \
src/precession.o src/sidereal_time.o src/planet_funcs.o
# sat_code C++ sources (compiled with g++, linked with extern "C" symbols) # sat_code C++ sources (compiled with g++, linked with extern "C" symbols)
SAT_CODE_DIR = lib/sat_code SAT_CODE_DIR = lib/sat_code
@ -17,7 +20,8 @@ SAT_CODE_OBJS = $(SAT_CODE_SRCS:.cpp=.o)
OBJS += $(SAT_CODE_OBJS) OBJS += $(SAT_CODE_OBJS)
# Regression tests # Regression tests
REGRESS = tle_parse sgp4_propagate coord_transforms pass_prediction gist_index convenience REGRESS = tle_parse sgp4_propagate coord_transforms pass_prediction gist_index convenience \
star_observe kepler_comet planet_observe
REGRESS_OPTS = --inputdir=test REGRESS_OPTS = --inputdir=test
# Need C++ runtime for sat_code # Need C++ runtime for sat_code

2
pg_orbit.control
View File

@ -1,4 +1,4 @@
comment = 'Orbital mechanics types and functions for PostgreSQL' comment = 'Orbital mechanics types and functions for PostgreSQL'
default_version = '0.1.0' default_version = '0.2.0'
module_pathname = '$libdir/pg_orbit' module_pathname = '$libdir/pg_orbit'
relocatable = true relocatable = true

122
sql/pg_orbit--0.1.0--0.2.0.sql Normal file
View File

@ -0,0 +1,122 @@
-- pg_orbit 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.';

638
sql/pg_orbit--0.2.0.sql Normal file
View File

@ -0,0 +1,638 @@
-- pg_orbit -- 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_orbit 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.';

170
src/astro_math.h Normal file
View File

@ -0,0 +1,170 @@
/*
* astro_math.h -- Shared astronomical math for pg_orbit
*
* Static inline functions used by star_funcs.c, kepler_funcs.c,
* and future planet/moon observation code.
*
* Using static inline preserves the project convention of no
* cross-translation-unit symbol coupling.
*/
#ifndef PG_ORBIT_ASTRO_MATH_H
#define PG_ORBIT_ASTRO_MATH_H
#include <math.h>
#include "types.h"
#define DEG_TO_RAD (M_PI / 180.0)
#define RAD_TO_DEG (180.0 / M_PI)
#define ARCSEC_TO_RAD (M_PI / (180.0 * 3600.0))
/*
* Greenwich Mean Sidereal Time from Julian date.
* Vallado, "Fundamentals of Astrodynamics", Eq. 3-47.
* Returns angle in radians.
*/
static inline double
gmst_from_jd(double jd)
{
double t_ut1 = (jd - 2451545.0) / 36525.0;
double gmst = 67310.54841
+ (876600.0 * 3600.0 + 8640184.812866) * t_ut1
+ 0.093104 * t_ut1 * t_ut1
- 6.2e-6 * t_ut1 * t_ut1 * t_ut1;
gmst = fmod(gmst * M_PI / 43200.0, 2.0 * M_PI);
if (gmst < 0.0)
gmst += 2.0 * M_PI;
return gmst;
}
/*
* IAU 1976 precession: rotate J2000 equatorial to mean equatorial of date.
*
* Source: Lieske (1979), A&A 73, 282.
* Three angles zeta_A, z_A, theta_A define the precession rotation.
* Valid within ~1 arcsecond for several centuries around J2000.
*/
static inline void
precess_j2000_to_date(double jd, double ra_j2000, double dec_j2000,
double *ra_date, double *dec_date)
{
double T, T2, T3;
double zeta_A, z_A, theta_A;
double cos_dec, sin_dec, cos_ra_zeta, sin_ra_zeta;
double cos_theta, sin_theta;
double A, B, C;
T = (jd - J2000_JD) / 36525.0;
T2 = T * T;
T3 = T2 * T;
/* Precession angles in arcseconds (Lieske 1979) */
zeta_A = (2306.2181 + 1.39656 * T - 0.000139 * T2) * T
+ (0.30188 - 0.000344 * T) * T2 + 0.017998 * T3;
z_A = (2306.2181 + 1.39656 * T - 0.000139 * T2) * T
+ (1.09468 + 0.000066 * T) * T2 + 0.018203 * T3;
theta_A = (2004.3109 - 0.85330 * T - 0.000217 * T2) * T
- (0.42665 + 0.000217 * T) * T2 - 0.041833 * T3;
/* Arcseconds to radians */
zeta_A *= ARCSEC_TO_RAD;
z_A *= ARCSEC_TO_RAD;
theta_A *= ARCSEC_TO_RAD;
/* Direct formula: R3(-z_A) R2(theta_A) R3(-zeta_A) applied to (ra, dec) */
cos_dec = cos(dec_j2000);
sin_dec = sin(dec_j2000);
cos_ra_zeta = cos(ra_j2000 + zeta_A);
sin_ra_zeta = sin(ra_j2000 + zeta_A);
cos_theta = cos(theta_A);
sin_theta = sin(theta_A);
A = cos_dec * sin_ra_zeta;
B = cos_theta * cos_dec * cos_ra_zeta - sin_theta * sin_dec;
C = sin_theta * cos_dec * cos_ra_zeta + cos_theta * sin_dec;
*dec_date = asin(C);
*ra_date = atan2(A, B) + z_A;
if (*ra_date < 0.0)
*ra_date += 2.0 * M_PI;
if (*ra_date >= 2.0 * M_PI)
*ra_date -= 2.0 * M_PI;
}
/*
* Equatorial (hour angle, declination) to horizontal (azimuth, elevation).
* All angles in radians.
* Azimuth convention: 0=N, pi/2=E, pi=S, 3*pi/2=W.
*/
static inline void
equatorial_to_horizontal(double ha, double dec, double lat,
double *az, double *el)
{
double sin_lat = sin(lat);
double cos_lat = cos(lat);
double sin_dec = sin(dec);
double cos_dec = cos(dec);
double cos_ha = cos(ha);
double y, x;
*el = asin(sin_lat * sin_dec + cos_lat * cos_dec * cos_ha);
y = -cos_dec * sin(ha);
x = cos_lat * sin_dec - sin_lat * cos_dec * cos_ha;
*az = atan2(y, x);
if (*az < 0.0)
*az += 2.0 * M_PI;
}
/*
* Ecliptic J2000 to equatorial J2000.
* Simple rotation around X-axis by -obliquity.
*/
static inline void
ecliptic_to_equatorial(const double *ecl, double *equ)
{
double cos_eps = cos(OBLIQUITY_J2000);
double sin_eps = sin(OBLIQUITY_J2000);
equ[0] = ecl[0];
equ[1] = ecl[1] * cos_eps - ecl[2] * sin_eps;
equ[2] = ecl[1] * sin_eps + ecl[2] * cos_eps;
}
/*
* Equatorial J2000 to ecliptic J2000.
* Rotation around X-axis by +obliquity.
*/
static inline void
equatorial_to_ecliptic(const double *equ, double *ecl)
{
double cos_eps = cos(OBLIQUITY_J2000);
double sin_eps = sin(OBLIQUITY_J2000);
ecl[0] = equ[0];
ecl[1] = equ[1] * cos_eps + equ[2] * sin_eps;
ecl[2] = -equ[1] * sin_eps + equ[2] * cos_eps;
}
/*
* Cartesian to spherical: (x, y, z) -> (ra, dec, dist).
* ra in [0, 2*pi), dec in [-pi/2, pi/2], dist in same units as input.
*/
static inline void
cartesian_to_spherical(const double *xyz, double *ra, double *dec, double *dist)
{
*dist = sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1] + xyz[2] * xyz[2]);
*dec = asin(xyz[2] / *dist);
*ra = atan2(xyz[1], xyz[0]);
if (*ra < 0.0)
*ra += 2.0 * M_PI;
}
#endif /* PG_ORBIT_ASTRO_MATH_H */

149
src/elliptic_to_rectangular.c Normal file
View File

@ -0,0 +1,149 @@
/************************************************************************
The code in this file is heavily inspired by the TASS17 and GUST86 theories
found on
ftp://ftp.imcce.fr/pub/ephem/satel
I (Johannes Gajdosik) have just taken the Fortran code and data
obtained from above and rearranged it into this piece of software.
I can neither allow nor forbid the above theories.
The copyright notice below covers just my work,
that is the compilation of the data obtained from above
into the software supplied in this file.
Copyright (c) 2005 Johannes Gajdosik
Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the "Software"),
to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense,
and/or sell copies of the Software, and to permit persons to whom the
Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included
in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
****************************************************************/
#include "elliptic_to_rectangular.h"
#include <math.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
/*
Given the orbital elements at some time t0 calculate the
rectangular coordinates at time (t0+dt).
mu = G*(m1+m2) .. gravitational constant of the two body problem
a .. semi major axis
n = mean motion = 2*M_PI/(orbit period)
elem[0] .. irrelevant (either a (called by EllipticToRectangularA()) or n (called by EllipticToRectangularN())
elem[1] .. L
elem[2] .. K=e*cos(Omega+omega)
elem[3] .. H=e*sin(Omega+omega)
elem[4] .. Q=sin(i/2)*cos(Omega)
elem[5] .. P=sin(i/2)*sin(Omega)
Omega = longitude of ascending node
omega = argument of pericenter
L = mean longitude = Omega + omega + M
M = mean anomaly
i = inclination
e = eccentricity
Units are suspected to be: Julian days, AU, rad
*/
static void
EllipticToRectangular(const double a,const double n,
const double elem[6],const double dt,double xyz[]) {
const double L = fmod(elem[1]+n*dt,2.0*M_PI);
/* solve Keplers equation
x = L - elem[2]*sin(x) + elem[3]*cos(x)
not by trivially iterating
x_0 = L
x_{j+1} = L - elem[2]*sin(x_j) + elem[3]*cos(x_j)
but instead by Newton approximation:
0 = f(x) = x - L - elem[2]*sin(x) + elem[3]*cos(x)
f'(x) = 1 - elem[2]*cos(x) - elem[3]*sin(x)
x_0 = L or whatever, perhaps first step of trivial iteration
x_{j+1} = x_j - f(x_j)/f'(x_j)
*/
double Le = L - elem[2]*sin(L) + elem[3]*cos(L);
for (;;) {
const double cLe = cos(Le);
const double sLe = sin(Le);
/* for eccentricity < 1 we have denominator > 0 */
const double dLe = (L - Le + elem[2]*sLe - elem[3]*cLe)
/ (1.0 - elem[2]*cLe - elem[3]*sLe);
Le += dLe;
if (fabs(dLe) <= 1e-14) break; /* L1: <1e-12 */
}
{
const double cLe = cos(Le);
const double sLe = sin(Le);
const double dlf = -elem[2]*sLe + elem[3]*cLe;
const double phi = sqrt(1.0 - elem[2]*elem[2] - elem[3]*elem[3]);
const double psi = 1.0 / (1.0 + phi);
const double x1 = a * (cLe - elem[2] - psi*dlf*elem[3]);
const double y1 = a * (sLe - elem[3] + psi*dlf*elem[2]);
const double elem_4q = elem[4] * elem[4]; // Q²
const double elem_5q = elem[5] * elem[5]; // P²
const double dwho = 2.0 * sqrt(1.0 - elem_4q - elem_5q);
const double rtp = 1.0 - elem_5q - elem_5q;
const double rtq = 1.0 - elem_4q - elem_4q;
const double rdg = 2.0 * elem[5] * elem[4];
xyz[0] = x1 * rtp + y1 * rdg;
xyz[1] = x1 * rdg + y1 * rtq;
xyz[2] = (-x1 * elem[5] + y1 * elem[4]) * dwho;
// GZ 2017-11: Re-enable these lines, they provide velocity!
const double rsam1 = -elem[2]*cLe - elem[3]*sLe;
const double h = a*n / (1.0 + rsam1);
const double vx1 = h * (-sLe - psi*rsam1*elem[3]);
const double vy1 = h * ( cLe + psi*rsam1*elem[2]);
xyz[3] = vx1 * rtp + vy1 * rdg;
xyz[4] = vx1 * rdg + vy1 * rtq;
xyz[5] = (-vx1 * elem[5] + vy1 * elem[4]) * dwho;
}
}
void EllipticToRectangularN(double mu,const double elem[6],double dt,
double xyz[]) {
const double n = elem[0];
#if defined __USE_MISC || defined __USE_XOPEN_EXTENDED || defined __USE_ISOC99
/* linux math.h declares cbrt: */
const double a = cbrt(mu/(n*n));
#else
const double a = exp(log(mu/(n*n))/3.0);
#endif
EllipticToRectangular(a,n,elem,dt,xyz);
}
void EllipticToRectangularA(double mu,const double elem[6],double dt,
double xyz[]) {
const double a = elem[0];
const double n = sqrt(mu/(a*a*a)); // mean motion
EllipticToRectangular(a,n,elem,dt,xyz);
}

84
src/elliptic_to_rectangular.h Normal file
View File

@ -0,0 +1,84 @@
/************************************************************************
The code in this file is heavily inspired by the TASS17 and GUST86 theories
found on
ftp://ftp.imcce.fr/pub/ephem/satel
I (Johannes Gajdosik) have just taken the Fortran code and data
obtained from above and rearranged it into this piece of software.
I can neither allow nor forbid the above theories.
The copyright notice below covers just my work,
that is the compilation of the data obtained from above
into the software supplied in this file.
Copyright (c) 2005 Johannes Gajdosik
Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the "Software"),
to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense,
and/or sell copies of the Software, and to permit persons to whom the
Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included
in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
****************************************************************/
#ifndef ELLIPTIC_TO_RECTANGULAR_H
#define ELLIPTIC_TO_RECTANGULAR_H
#ifdef __cplusplus
extern "C" {
#endif
/*
Given the orbital elements at some time t0 calculate the
rectangular coordinates at time (t0+dt).
mu = G*(m1+m2) .. gravitational constant of the two body problem
a .. semi major axis
n = mean motion = 2*M_PI/(orbit period)
elem[0] .. either a (EllipticToRectangularA()) or n (EllipticToRectangularN())
elem[1] .. L
elem[2] .. K=e*cos(Omega+omega)
elem[3] .. H=e*sin(Omega+omega)
elem[4] .. Q=sin(i/2)*cos(Omega)
elem[5] .. P=sin(i/2)*sin(Omega)
Omega = longitude of ascending node
omega = argument of pericenter
L = mean longitude = Omega + omega + M
M = mean anomaly
i = inclination
e = eccentricity
Units are suspected to be: Julian days, AU, rad
Results:
xyz[0,1,2]=Position [AU]
xyz[3,4,5]=Velocity [AU/d]
*/
void EllipticToRectangularN(double mu,const double elem[6],double dt,
double xyz[]);
void EllipticToRectangularA(double mu,const double elem[6],double dt,
double xyz[]);
#ifdef __cplusplus
}
#endif
#endif

162761
src/elp82b.c Normal file
View File

File diff suppressed because it is too large Load Diff

92
src/elp82b.h Normal file
View File

@ -0,0 +1,92 @@
/************************************************************************
LUNAR SOLUTION ELP2000-82B
by Chapront-Touze M., Chapront J.
ftp://ftp.imcce.fr/pub/ephem/moon/elp82b
I (Johannes Gajdosik) have just taken the Fortran code and data
obtained from above and used it to create this piece of software.
I can neither allow nor forbid the usage of ELP2000-82B.
The copyright notice below covers not the works of
Chapront-Touze M. and Chapront J., but just my work,
that is the compilation and rearrangement of
the Fortran code and data obtained from above
into the software supplied in this file.
Copyright (c) 2005 Johannes Gajdosik
Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the "Software"),
to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense,
and/or sell copies of the Software, and to permit persons to whom the
Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included
in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
My implementation of ELP2000-82B has the following modifications compared to
the original Fortran code:
1) fundamentally rearrange the series into optimized instructions
for fast calculation of the results
2) units are: julian day, AU
This file is derived from the Stellarium project
(https://github.com/Stellarium/stellarium). Modified for pg_orbit:
removed static mutable state for thread safety (PostgreSQL PARALLEL SAFE).
****************************************************************/
#ifndef ELP82B_H
#define ELP82B_H
#ifdef __cplusplus
extern "C" {
#endif
void GetElp82bCoor(double jd,double xyz[3]);
/* Return the rectangular coordinates of the earths moon
on the given julian date jd expressed in dynamical time (TAI+32.184s).
The origin of the xyz-coordinates is the center of the earth.
The reference frame is "dynamical equinox and ecliptic J2000",
which is the reference frame in VSOP87 and VSOP87A.
According to vsop87.doc VSOP87 coordinates can be transformed to FK5 by
X cos(psi) -sin(psi) 0 1 0 0 X
Y = sin(psi) cos(psi) 0 * 0 cos(eps) -sin(eps) * Y
Z FK5 0 0 1 0 sin(eps) cos(eps) Z VSOP87
with psi = -0.0000275 degrees = -0.099 arcsec and
eps = 23.4392803055556 degrees = 23d26m21.4091sec.
http://ssd.jpl.nasa.gov/horizons_doc.html#frames says:
"J2000" selects an Earth Mean-Equator and dynamical Equinox of
Epoch J2000.0 inertial reference system, where the Epoch of J2000.0
is the Julian date 2451545.0. "Mean" indicates nutation effects are
ignored in the frame definition. The system is aligned with the
IAU-sponsored J2000 frame of the Radio Source Catalog of the
International Earth Rotational Service (ICRF).
The ICRF is thought to differ from FK5 by at most 0.01 arcsec.
From this I conclude that in the context of stellarium
ICRF, J2000 and FK5 are the same, while the transformation
ICRF <-> VSOP87 must be done with the matrix given above.
*/
#ifdef __cplusplus
}
#endif
#endif

427
src/kepler_funcs.c Normal file
View File

@ -0,0 +1,427 @@
/*
* kepler_funcs.c -- Keplerian propagation and heliocentric type
*
* Two-body propagation from classical orbital elements for comets
* and asteroids. Handles elliptic (e<1), near-parabolic (e~1),
* and hyperbolic (e>1) orbits.
*
* Also provides the heliocentric type (ecliptic J2000 position in AU)
* and comet_observe() which computes topocentric coordinates given
* the comet's orbital elements and Earth's heliocentric position.
*/
#include "postgres.h"
#include "fmgr.h"
#include "funcapi.h"
#include "utils/timestamp.h"
#include "libpq/pqformat.h"
#include "types.h"
#include "astro_math.h"
#include <math.h>
/* Heliocentric type I/O */
PG_FUNCTION_INFO_V1(heliocentric_in);
PG_FUNCTION_INFO_V1(heliocentric_out);
PG_FUNCTION_INFO_V1(heliocentric_recv);
PG_FUNCTION_INFO_V1(heliocentric_send);
PG_FUNCTION_INFO_V1(helio_x);
PG_FUNCTION_INFO_V1(helio_y);
PG_FUNCTION_INFO_V1(helio_z);
PG_FUNCTION_INFO_V1(helio_distance);
/* Propagation functions */
PG_FUNCTION_INFO_V1(kepler_propagate);
PG_FUNCTION_INFO_V1(comet_observe);
/* ================================================================
* Kepler equation solvers
* ================================================================
*/
/*
* Elliptic Kepler equation: M = E - e*sin(E).
* Newton-Raphson, converges in 3-5 iterations for e < 0.99.
*/
static double
solve_kepler_elliptic(double M, double e)
{
double E = M;
int i;
/* Better initial guess for high eccentricity */
if (e > 0.8)
E = M_PI;
for (i = 0; i < 30; i++)
{
double dE = (E - e * sin(E) - M) / (1.0 - e * cos(E));
E -= dE;
if (fabs(dE) < 1.0e-15)
break;
}
return E;
}
/*
* Hyperbolic Kepler equation: M = e*sinh(H) - H.
* Newton-Raphson.
*/
static double
solve_kepler_hyperbolic(double M, double e)
{
/* Initial guess: asinh(M/e) for large M, M for small */
double H = (fabs(M) > 1.0) ? copysign(log(fabs(M) / e + 1.0), M) : M;
int i;
for (i = 0; i < 30; i++)
{
double dH = (e * sinh(H) - H - M) / (e * cosh(H) - 1.0);
H -= dH;
if (fabs(dH) < 1.0e-15)
break;
}
return H;
}
/*
* Near-parabolic: Barker's equation W^3 + 3W - 3M = 0.
* Returns true anomaly directly.
*/
static double
solve_kepler_parabolic(double M)
{
double W, f, fp;
int i;
/* Cardano's formula for W^3 + 3W - 3M = 0 gives decent initial guess */
double disc = sqrt(1.0 + M * M);
W = cbrt(3.0 * M + 3.0 * disc) - cbrt(-3.0 * M + 3.0 * disc);
for (i = 0; i < 30; i++)
{
f = W * W * W + 3.0 * W - 3.0 * M;
fp = 3.0 * W * W + 3.0;
W -= f / fp;
if (fabs(f / fp) < 1.0e-15)
break;
}
return 2.0 * atan(W);
}
/* ================================================================
* Internal: Keplerian position from orbital elements
*
* Input: q (AU), e, inc (rad), omega (rad), Omega (rad),
* T_peri (JD), jd (observation JD)
* Output: pos[3] in AU, ecliptic J2000 frame
* ================================================================
*/
static void
kepler_position(double q, double e, double inc, double omega, double Omega,
double T_peri, double jd, double pos[3])
{
double dt = jd - T_peri;
double v = 0.0, r = 0.0;
double x_orb, y_orb;
double cos_om, sin_om, cos_Om, sin_Om, cos_i, sin_i;
double Px, Py, Pz, Qx, Qy, Qz;
if (e < 0.99)
{
double a = q / (1.0 - e);
double n = GAUSS_K / (a * sqrt(a));
double M = fmod(n * dt, 2.0 * M_PI);
double E;
if (M < 0.0)
M += 2.0 * M_PI;
E = solve_kepler_elliptic(M, e);
v = 2.0 * atan2(sqrt(1.0 + e) * sin(E / 2.0),
sqrt(1.0 - e) * cos(E / 2.0));
r = a * (1.0 - e * cos(E));
}
else if (e > 1.01)
{
double a = q / (e - 1.0);
double n = GAUSS_K / (a * sqrt(a));
double M = n * dt;
double H = solve_kepler_hyperbolic(M, e);
v = 2.0 * atan2(sqrt(e + 1.0) * tanh(H / 2.0),
sqrt(e - 1.0));
r = a * (e * cosh(H) - 1.0);
}
else
{
double n = GAUSS_K * sqrt(1.0 / (2.0 * q * q * q));
double M = n * dt;
v = solve_kepler_parabolic(M);
r = q * (1.0 + tan(v / 2.0) * tan(v / 2.0));
}
x_orb = r * cos(v);
y_orb = r * sin(v);
/* Perifocal-to-ecliptic rotation (P, Q vectors) */
cos_om = cos(omega);
sin_om = sin(omega);
cos_Om = cos(Omega);
sin_Om = sin(Omega);
cos_i = cos(inc);
sin_i = sin(inc);
Px = cos_Om * cos_om - sin_Om * sin_om * cos_i;
Py = sin_Om * cos_om + cos_Om * sin_om * cos_i;
Pz = sin_om * sin_i;
Qx = -cos_Om * sin_om - sin_Om * cos_om * cos_i;
Qy = -sin_Om * sin_om + cos_Om * cos_om * cos_i;
Qz = cos_om * sin_i;
pos[0] = Px * x_orb + Qx * y_orb;
pos[1] = Py * x_orb + Qy * y_orb;
pos[2] = Pz * x_orb + Qz * y_orb;
}
/* ================================================================
* Heliocentric type I/O
*
* Text format: (x_au, y_au, z_au)
* ================================================================
*/
Datum
heliocentric_in(PG_FUNCTION_ARGS)
{
char *str = PG_GETARG_CSTRING(0);
pg_heliocentric *result;
double x, y, z;
int nfields;
result = (pg_heliocentric *) palloc(sizeof(pg_heliocentric));
nfields = sscanf(str, " ( %lf , %lf , %lf )", &x, &y, &z);
if (nfields != 3)
ereport(ERROR,
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
errmsg("invalid input syntax for type heliocentric: \"%s\"", str),
errhint("Expected (x_au,y_au,z_au).")));
result->x = x;
result->y = y;
result->z = z;
PG_RETURN_POINTER(result);
}
Datum
heliocentric_out(PG_FUNCTION_ARGS)
{
pg_heliocentric *h = (pg_heliocentric *) PG_GETARG_POINTER(0);
PG_RETURN_CSTRING(psprintf("(%.10f,%.10f,%.10f)", h->x, h->y, h->z));
}
Datum
heliocentric_recv(PG_FUNCTION_ARGS)
{
StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
pg_heliocentric *result;
result = (pg_heliocentric *) palloc(sizeof(pg_heliocentric));
result->x = pq_getmsgfloat8(buf);
result->y = pq_getmsgfloat8(buf);
result->z = pq_getmsgfloat8(buf);
PG_RETURN_POINTER(result);
}
Datum
heliocentric_send(PG_FUNCTION_ARGS)
{
pg_heliocentric *h = (pg_heliocentric *) PG_GETARG_POINTER(0);
StringInfoData buf;
pq_begintypsend(&buf);
pq_sendfloat8(&buf, h->x);
pq_sendfloat8(&buf, h->y);
pq_sendfloat8(&buf, h->z);
PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
}
/* --- heliocentric accessors --- */
Datum
helio_x(PG_FUNCTION_ARGS)
{
pg_heliocentric *h = (pg_heliocentric *) PG_GETARG_POINTER(0);
PG_RETURN_FLOAT8(h->x);
}
Datum
helio_y(PG_FUNCTION_ARGS)
{
pg_heliocentric *h = (pg_heliocentric *) PG_GETARG_POINTER(0);
PG_RETURN_FLOAT8(h->y);
}
Datum
helio_z(PG_FUNCTION_ARGS)
{
pg_heliocentric *h = (pg_heliocentric *) PG_GETARG_POINTER(0);
PG_RETURN_FLOAT8(h->z);
}
Datum
helio_distance(PG_FUNCTION_ARGS)
{
pg_heliocentric *h = (pg_heliocentric *) PG_GETARG_POINTER(0);
PG_RETURN_FLOAT8(sqrt(h->x * h->x + h->y * h->y + h->z * h->z));
}
/* ================================================================
* kepler_propagate(q, e, i_deg, omega_deg, Omega_deg, T_peri_jd,
* timestamptz) -> heliocentric
*
* Propagate a body from classical orbital elements via two-body
* Keplerian dynamics. Returns heliocentric ecliptic J2000 position.
* ================================================================
*/
Datum
kepler_propagate(PG_FUNCTION_ARGS)
{
double q_au = PG_GETARG_FLOAT8(0);
double ecc = PG_GETARG_FLOAT8(1);
double inc_deg = PG_GETARG_FLOAT8(2);
double omega_deg = PG_GETARG_FLOAT8(3);
double Omega_deg = PG_GETARG_FLOAT8(4);
double T_peri_jd = PG_GETARG_FLOAT8(5);
int64 ts = PG_GETARG_INT64(6);
double jd;
double pos[3];
pg_heliocentric *result;
if (q_au <= 0.0)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("perihelion distance must be positive: %.6f", q_au)));
if (ecc < 0.0)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("eccentricity must be non-negative: %.6f", ecc)));
jd = timestamptz_to_jd(ts);
kepler_position(q_au, ecc,
inc_deg * DEG_TO_RAD,
omega_deg * DEG_TO_RAD,
Omega_deg * DEG_TO_RAD,
T_peri_jd, jd, pos);
result = (pg_heliocentric *) palloc(sizeof(pg_heliocentric));
result->x = pos[0];
result->y = pos[1];
result->z = pos[2];
PG_RETURN_POINTER(result);
}
/* ================================================================
* comet_observe(q, e, i, omega, Omega, T_peri,
* earth_x_au, earth_y_au, earth_z_au,
* observer, timestamptz) -> topocentric
*
* Full pipeline: Keplerian propagation -> geocentric position ->
* equatorial J2000 -> precession to date -> topocentric az/el.
*
* Earth's heliocentric ecliptic J2000 position must be provided.
* In Phase 1, Skyfield or similar supplies this.
* In Phase 2, planet_heliocentric(BODY_EARTH, ts) provides it.
* ================================================================
*/
Datum
comet_observe(PG_FUNCTION_ARGS)
{
double q_au = PG_GETARG_FLOAT8(0);
double ecc = PG_GETARG_FLOAT8(1);
double inc_deg = PG_GETARG_FLOAT8(2);
double omega_deg = PG_GETARG_FLOAT8(3);
double Omega_deg = PG_GETARG_FLOAT8(4);
double T_peri_jd = PG_GETARG_FLOAT8(5);
double earth_x = PG_GETARG_FLOAT8(6);
double earth_y = PG_GETARG_FLOAT8(7);
double earth_z = PG_GETARG_FLOAT8(8);
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(9);
int64 ts = PG_GETARG_INT64(10);
double jd;
double comet_helio[3];
double geo_ecl[3];
double geo_equ[3];
double ra_j2000, dec_j2000, geo_dist;
double ra_date, dec_date;
double gmst_val, lst, ha;
double az, el;
pg_topocentric *result;
if (q_au <= 0.0)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("perihelion distance must be positive: %.6f", q_au)));
if (ecc < 0.0)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("eccentricity must be non-negative: %.6f", ecc)));
jd = timestamptz_to_jd(ts);
/* Comet heliocentric ecliptic J2000 position */
kepler_position(q_au, ecc,
inc_deg * DEG_TO_RAD,
omega_deg * DEG_TO_RAD,
Omega_deg * DEG_TO_RAD,
T_peri_jd, jd, comet_helio);
/* Geocentric ecliptic position = comet_helio - earth_helio */
geo_ecl[0] = comet_helio[0] - earth_x;
geo_ecl[1] = comet_helio[1] - earth_y;
geo_ecl[2] = comet_helio[2] - earth_z;
/* Ecliptic J2000 -> equatorial J2000 */
ecliptic_to_equatorial(geo_ecl, geo_equ);
/* Cartesian -> spherical (RA, Dec, distance) */
cartesian_to_spherical(geo_equ, &ra_j2000, &dec_j2000, &geo_dist);
/* Precess J2000 RA/Dec to date */
precess_j2000_to_date(jd, ra_j2000, dec_j2000, &ra_date, &dec_date);
/* Hour angle and az/el */
gmst_val = gmst_from_jd(jd);
lst = gmst_val + obs->lon;
ha = lst - ra_date;
equatorial_to_horizontal(ha, dec_date, obs->lat, &az, &el);
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
result->azimuth = az;
result->elevation = el;
result->range_km = geo_dist * AU_KM;
result->range_rate = 0.0; /* no velocity computation in Phase 1 */
PG_RETURN_POINTER(result);
}

236
src/planet_funcs.c Normal file
View File

@ -0,0 +1,236 @@
/*
* planet_funcs.c -- VSOP87 planet and ELP82B Moon observation
*
* SQL functions for planetary/Sun/Moon position and observation.
* Wraps VSOP87 (Bretagnon & Francou) for 8 planets + Sun, and
* ELP2000-82B (Chapront-Touze) for the Moon.
*
* The observation pipeline:
* 1. Heliocentric ecliptic J2000 position (VSOP87 or ELP82B)
* 2. Geocentric ecliptic (subtract Earth's heliocentric)
* 3. Ecliptic -> equatorial J2000 (obliquity rotation)
* 4. Spherical coordinates (RA, Dec, distance)
* 5. Precession J2000 -> date (IAU 1976, Lieske)
* 6. Sidereal time -> hour angle -> az/el
* 7. Return topocentric result
*/
#include "postgres.h"
#include "fmgr.h"
#include "funcapi.h"
#include "utils/timestamp.h"
#include "types.h"
#include "astro_math.h"
#include "vsop87.h"
#include "elp82b.h"
#include <math.h>
PG_FUNCTION_INFO_V1(planet_heliocentric);
PG_FUNCTION_INFO_V1(planet_observe);
PG_FUNCTION_INFO_V1(sun_observe);
PG_FUNCTION_INFO_V1(moon_observe);
/* ================================================================
* Internal: geocentric observation pipeline
*
* Takes geocentric ecliptic J2000 position, observer location,
* and Julian date. Converts through equatorial, precesses to
* date, and computes topocentric az/el.
*
* geo_ecl_au[3] = geocentric ecliptic J2000 position in AU
* ================================================================
*/
static void
observe_from_geocentric(const double geo_ecl_au[3], double jd,
const pg_observer *obs, pg_topocentric *result)
{
double geo_equ[3];
double ra_j2000, dec_j2000, geo_dist;
double ra_date, dec_date;
double gmst_val, lst, ha;
double az, el;
/* Ecliptic J2000 -> equatorial J2000 */
ecliptic_to_equatorial(geo_ecl_au, geo_equ);
/* Cartesian -> spherical */
cartesian_to_spherical(geo_equ, &ra_j2000, &dec_j2000, &geo_dist);
/* Precess J2000 -> date */
precess_j2000_to_date(jd, ra_j2000, dec_j2000, &ra_date, &dec_date);
/* Hour angle and az/el */
gmst_val = gmst_from_jd(jd);
lst = gmst_val + obs->lon;
ha = lst - ra_date;
equatorial_to_horizontal(ha, dec_date, obs->lat, &az, &el);
result->azimuth = az;
result->elevation = el;
result->range_km = geo_dist * AU_KM;
result->range_rate = 0.0; /* no velocity computation yet */
}
/* ================================================================
* planet_heliocentric(body_id int, timestamptz) -> heliocentric
*
* Returns the heliocentric ecliptic J2000 position of a planet
* (or the Sun, as the barycenter proxy at (0,0,0)).
*
* Body IDs: 0=Sun, 1=Mercury, ..., 8=Neptune
* ================================================================
*/
Datum
planet_heliocentric(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
int64 ts = PG_GETARG_INT64(1);
double jd;
double xyz[6];
int vsop_body;
pg_heliocentric *result;
if (body_id == BODY_SUN)
{
/* Sun is at the origin of heliocentric coordinates */
result = (pg_heliocentric *) palloc(sizeof(pg_heliocentric));
result->x = 0.0;
result->y = 0.0;
result->z = 0.0;
PG_RETURN_POINTER(result);
}
if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("invalid body_id %d: must be 0 (Sun) or 1-8 (Mercury-Neptune)",
body_id)));
jd = timestamptz_to_jd(ts);
vsop_body = body_id - 1; /* pg_orbit 1-based -> VSOP87 0-based */
GetVsop87Coor(jd, vsop_body, xyz);
result = (pg_heliocentric *) palloc(sizeof(pg_heliocentric));
result->x = xyz[0];
result->y = xyz[1];
result->z = xyz[2];
PG_RETURN_POINTER(result);
}
/* ================================================================
* planet_observe(body_id int, observer, timestamptz) -> topocentric
*
* Full pipeline: VSOP87 position -> geocentric -> equatorial ->
* precession -> topocentric az/el.
*
* Body IDs: 1=Mercury, ..., 8=Neptune (not Sun or Moon)
* ================================================================
*/
Datum
planet_observe(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(1);
int64 ts = PG_GETARG_INT64(2);
double jd;
double earth_xyz[6];
double planet_xyz[6];
double geo_ecl[3];
int vsop_body;
pg_topocentric *result;
if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("planet_observe: body_id %d must be 1-8 (Mercury-Neptune)",
body_id)));
if (body_id == BODY_EARTH)
ereport(ERROR,
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
errmsg("cannot observe Earth from Earth")));
jd = timestamptz_to_jd(ts);
/* Earth's heliocentric position */
GetVsop87Coor(jd, 2, earth_xyz); /* VSOP87 body 2 = Earth */
/* Target planet heliocentric position */
vsop_body = body_id - 1;
GetVsop87Coor(jd, vsop_body, planet_xyz);
/* Geocentric ecliptic = planet - Earth */
geo_ecl[0] = planet_xyz[0] - earth_xyz[0];
geo_ecl[1] = planet_xyz[1] - earth_xyz[1];
geo_ecl[2] = planet_xyz[2] - earth_xyz[2];
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
observe_from_geocentric(geo_ecl, jd, obs, result);
PG_RETURN_POINTER(result);
}
/* ================================================================
* sun_observe(observer, timestamptz) -> topocentric
*
* The Sun's geocentric position is the negation of Earth's
* heliocentric VSOP87 position.
* ================================================================
*/
Datum
sun_observe(PG_FUNCTION_ARGS)
{
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
int64 ts = PG_GETARG_INT64(1);
double jd;
double earth_xyz[6];
double geo_ecl[3];
pg_topocentric *result;
jd = timestamptz_to_jd(ts);
/* Earth heliocentric -> Sun geocentric by negation */
GetVsop87Coor(jd, 2, earth_xyz);
geo_ecl[0] = -earth_xyz[0];
geo_ecl[1] = -earth_xyz[1];
geo_ecl[2] = -earth_xyz[2];
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
observe_from_geocentric(geo_ecl, jd, obs, result);
PG_RETURN_POINTER(result);
}
/* ================================================================
* moon_observe(observer, timestamptz) -> topocentric
*
* ELP2000-82B returns geocentric ecliptic J2000 Moon position.
* No Earth subtraction needed.
* ================================================================
*/
Datum
moon_observe(PG_FUNCTION_ARGS)
{
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
int64 ts = PG_GETARG_INT64(1);
double jd;
double moon_ecl[3];
pg_topocentric *result;
jd = timestamptz_to_jd(ts);
/* Moon geocentric ecliptic J2000 in AU */
GetElp82bCoor(jd, moon_ecl);
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
observe_from_geocentric(moon_ecl, jd, obs, result);
PG_RETURN_POINTER(result);
}

205
src/precession.c Normal file
View File

@ -0,0 +1,205 @@
/*
* precession.c -- IAU 2006 precession and simplified nutation
*
* Clean-room implementation from published standards:
*
* Precession polynomials:
* Capitaine, Wallace & Chapront (2003), "Expressions for IAU 2000
* precession quantities", A&A 412, 567-586.
* Hilton et al. (2006), "Report of the IAU Division I Working Group
* on Precession and the Ecliptic", Celestial Mechanics 94, 351-367.
* Coefficients are the IAU 2006 values adopted by resolution B1.
*
* Nutation fundamental arguments:
* Simon et al. (1994), "Numerical expressions for precession formulae
* and mean elements for the Moon and the planets", A&A 282, 663-683.
* IERS Conventions (2010), IERS Technical Note 36, Ch. 5.
*
* The 4-term nutation uses the dominant terms tabulated in the
* Explanatory Supplement to the Astronomical Almanac (3rd ed.),
* Seidelmann (ed.), Table 6.1. These are the same terms used
* by Meeus, "Astronomical Algorithms" (2nd ed.), Ch. 22.
*/
#include "precession.h"
#include <math.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
/* Julian centuries from J2000.0 */
#define T_FROM_JDE(jde) (((jde) - 2451545.0) / 36525.0)
/*
* get_precession_angles_vondrak
*
* IAU 2006 precession: polynomial expressions for the four
* fundamental precession quantities.
*
* T = Julian centuries of TDB from J2000.0
*
* Polynomial coefficients from Capitaine et al. (2003) Table 1
* as adopted by IAU 2006 (Hilton et al. 2006).
*
* psi_A: luni-solar precession
* omega_A: mean obliquity of date (inclination of equator on
* the fixed ecliptic of J2000)
* chi_A: planetary precession on the ecliptic
* epsilon_A: obliquity of the ecliptic (mean obliquity of date)
*/
void
get_precession_angles_vondrak(double jde,
double *epsilon_A,
double *chi_A,
double *omega_A,
double *psi_A)
{
double T, T2, T3, T4, T5;
T = T_FROM_JDE(jde);
T2 = T * T;
T3 = T2 * T;
T4 = T3 * T;
T5 = T4 * T;
/*
* psi_A: luni-solar precession in longitude
* Units: arcseconds
* Capitaine et al. (2003), Eq. (37); Hilton et al. (2006), Table 1
*/
*psi_A = 5038.481507 * T
- 1.0790069 * T2
- 0.00114045 * T3
+ 0.000132851 * T4
- 0.0000000951 * T5;
/*
* omega_A: inclination of the equator on the fixed ecliptic
* of J2000. This equals the mean obliquity at J2000 plus
* a small secular drift.
* Capitaine et al. (2003), Eq. (37)
*/
*omega_A = 84381.406
- 0.025754 * T
+ 0.0512623 * T2
- 0.00772503 * T3
- 0.000000467 * T4
+ 0.0000003337 * T5;
/*
* chi_A: planetary precession along the ecliptic
* Capitaine et al. (2003), Eq. (37)
*/
*chi_A = 10.556403 * T
- 2.3814292 * T2
- 0.00121197 * T3
+ 0.000170663 * T4
- 0.0000000560 * T5;
/*
* epsilon_A: mean obliquity of the ecliptic
* IAU 2006 value: starts at the IAU 2000 value 84381.406"
* and includes the secular variation.
* Capitaine et al. (2003), Eq. (37); Hilton et al. (2006)
*/
*epsilon_A = 84381.406
- 46.836769 * T
- 0.0001831 * T2
+ 0.00200340 * T3
- 0.000000576 * T4
- 0.0000000434 * T5;
}
/*
* get_nutation_angles_iau2000b
*
* Simplified nutation using the 4 dominant lunisolar terms.
*
* Fundamental arguments (IERS Conventions 2010, Eq. 5.43):
* Omega = longitude of Moon's ascending node
* L_sun = mean longitude of the Sun (mean anomaly not needed
* for these dominant terms; we use 2*L forms)
* L_moon = mean longitude of the Moon
*
* The 4-term approximation:
* delta_psi = -17.2 sin(Omega)
* - 1.32 sin(2*L_sun)
* - 0.23 sin(2*L_moon)
* + 0.21 sin(2*Omega)
*
* delta_epsilon = 9.2 cos(Omega)
* + 0.57 cos(2*L_sun)
* + 0.10 cos(2*L_moon)
* - 0.09 cos(2*Omega)
*
* Amplitudes in arcseconds. Accurate to ~0.5" for delta_psi
* and ~0.1" for delta_epsilon.
*
* Reference: Explanatory Supplement to the Astronomical Almanac
* (3rd ed.), Table 6.1; Meeus (1998), Ch. 22 Eq. 22.1.
*/
void
get_nutation_angles_iau2000b(double jde,
double *delta_psi,
double *delta_epsilon)
{
double T;
double omega, L_sun, L_moon;
double sin_omega, cos_omega;
double sin_2Lsun, cos_2Lsun;
double sin_2Lmoon, cos_2Lmoon;
double sin_2omega, cos_2omega;
double two_Lsun, two_Lmoon, two_omega;
T = T_FROM_JDE(jde);
/*
* Moon's ascending node longitude.
* Simon et al. (1994); IERS Conventions (2010), Eq. 5.43.
*/
omega = 125.04452 - 1934.136261 * T;
/*
* Mean longitude of the Sun (geometric, referred to mean equinox of date).
* Meeus (1998), Eq. 25.2, truncated.
*/
L_sun = 280.4665 + 36000.7698 * T;
/*
* Mean longitude of the Moon.
* Meeus (1998), Eq. 47.1, truncated.
*/
L_moon = 218.3165 + 481267.8813 * T;
/* Convert to radians */
omega *= M_PI / 180.0;
L_sun *= M_PI / 180.0;
L_moon *= M_PI / 180.0;
two_Lsun = 2.0 * L_sun;
two_Lmoon = 2.0 * L_moon;
two_omega = 2.0 * omega;
sin_omega = sin(omega);
cos_omega = cos(omega);
sin_2Lsun = sin(two_Lsun);
cos_2Lsun = cos(two_Lsun);
sin_2Lmoon = sin(two_Lmoon);
cos_2Lmoon = cos(two_Lmoon);
sin_2omega = sin(two_omega);
cos_2omega = cos(two_omega);
/* Nutation in longitude (arcseconds) */
*delta_psi = -17.20 * sin_omega
- 1.32 * sin_2Lsun
- 0.23 * sin_2Lmoon
+ 0.21 * sin_2omega;
/* Nutation in obliquity (arcseconds) */
*delta_epsilon = 9.20 * cos_omega
+ 0.57 * cos_2Lsun
+ 0.10 * cos_2Lmoon
- 0.09 * cos_2omega;
}

55
src/precession.h Normal file
View File

@ -0,0 +1,55 @@
/*
* precession.h -- IAU 2006 precession and simplified nutation
*
* Clean-room implementation from published standards:
* Precession: Capitaine, Wallace & Chapront (2003), A&A 412, 567
* (IAU 2006 polynomial expressions)
* Nutation: Simplified 4-term model from fundamental arguments
* documented in IERS Technical Note 32, Ch. 5.
*
* No PostgreSQL dependencies -- pure math, suitable for standalone use.
*/
#ifndef PG_ORBIT_PRECESSION_H
#define PG_ORBIT_PRECESSION_H
/*
* get_precession_angles_vondrak -- IAU 2006 precession angles
*
* Computes the four fundamental precession quantities at JDE:
* psi_A -- luni-solar precession (arcsec)
* omega_A -- inclination of equator on J2000 ecliptic (arcsec)
* chi_A -- planetary precession (arcsec)
* epsilon_A -- obliquity of the ecliptic (arcsec)
*
* These are the IAU 2006 polynomial expressions valid to ~1 mas
* for several centuries around J2000. The full Vondrak (2011)
* periodic terms are not included here; they extend validity
* to +/-200,000 years but add sub-arcsecond corrections that
* are negligible for our use case.
*
* Output is in arcseconds.
*/
void get_precession_angles_vondrak(double jde,
double *epsilon_A,
double *chi_A,
double *omega_A,
double *psi_A);
/*
* get_nutation_angles_iau2000b -- simplified nutation
*
* Computes nutation in longitude (delta_psi) and nutation in
* obliquity (delta_epsilon) using the dominant 4 lunisolar terms.
*
* This matches the level of nutation correction used in SGP4
* coordinate transforms and is suitable for "what's up?" queries
* where ~1 arcsecond accuracy is sufficient.
*
* Output is in arcseconds.
*/
void get_nutation_angles_iau2000b(double jde,
double *delta_psi,
double *delta_epsilon);
#endif /* PG_ORBIT_PRECESSION_H */

140
src/sidereal_time.c Normal file
View File

@ -0,0 +1,140 @@
/*
* sidereal_time.c -- Greenwich sidereal time (mean and apparent)
*
* Clean-room implementation from published standards:
*
* Mean sidereal time:
* Capitaine, Guinot & McCarthy (2000), A&A 355, 398-405.
* IERS Conventions (2010), IERS Technical Note 36, Ch. 5,
* Eq. 5.29 (GMST as polynomial in UT1).
*
* The polynomial form is:
* GMST(UT1) = theta_0 + theta_1*T + theta_2*T^2 + theta_3*T^3
* where T is Julian centuries of UT1 from J2000.0, and the
* result is in seconds of time. Coefficients from Capitaine
* et al. (2000), Eq. (42), which are the values adopted by
* the IERS for consistency with IAU 2000/2006 precession.
*
* Apparent sidereal time:
* GAST = GMST + equation_of_equinoxes
* EqEq = delta_psi * cos(epsilon_A)
* IERS Conventions (2010), Eq. 5.30.
*/
#include "sidereal_time.h"
#include "precession.h"
#include <math.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
/* Arcseconds to radians */
#define ARCSEC_TO_RAD (M_PI / (180.0 * 3600.0))
/*
* get_mean_sidereal_time
*
* GMST in radians from the classical polynomial expression.
*
* The polynomial (Capitaine et al. 2000, Eq. 42):
* GMST = 67310.54841
* + (876600h + 8640184.812866) * T
* + 0.093104 * T^2
* - 6.2e-6 * T^3
*
* where T = (JD_UT1 - 2451545.0) / 36525.0
* and the result is in seconds of sidereal time.
*
* The constant 876600h = 876600 * 3600 = 3155760000 seconds
* accounts for complete rotations.
*
* To convert seconds of time to radians:
* 1 day = 86400 seconds of time = 2*pi radians
* 1 second of time = 2*pi/86400 = pi/43200 radians
*
* JDE is accepted for interface consistency but unused here;
* GMST is fundamentally a UT1 quantity.
*/
double
get_mean_sidereal_time(double JD, double JDE)
{
double T, T2, T3;
double gmst_sec;
double gmst_rad;
(void)JDE; /* GMST depends on UT1, not TDB */
T = (JD - 2451545.0) / 36525.0;
T2 = T * T;
T3 = T2 * T;
/*
* Polynomial in seconds of sidereal time.
* The large linear coefficient includes whole rotations
* accumulated since J2000.0.
*/
gmst_sec = 67310.54841
+ (876600.0 * 3600.0 + 8640184.812866) * T
+ 0.093104 * T2
- 6.2e-6 * T3;
/* Seconds of time -> radians, then normalize to [0, 2*pi) */
gmst_rad = fmod(gmst_sec * (M_PI / 43200.0), 2.0 * M_PI);
if (gmst_rad < 0.0)
gmst_rad += 2.0 * M_PI;
return gmst_rad;
}
/*
* get_apparent_sidereal_time
*
* GAST = GMST + equation of the equinoxes.
*
* The equation of the equinoxes (IERS 2010, Eq. 5.30):
* EqEq = delta_psi * cos(epsilon_A)
*
* where delta_psi is nutation in longitude and epsilon_A is
* the mean obliquity of the ecliptic, both in arcseconds.
*
* For the full IAU 2000A model, there is an additional
* "complementary terms" correction to EqEq (Eq. 5.31),
* but it is at the 0.1 mas level and negligible for our
* accuracy requirements.
*/
double
get_apparent_sidereal_time(double JD, double JDE)
{
double gmst;
double epsilon_A, chi_A, omega_A, psi_A;
double delta_psi, delta_epsilon;
double eq_eq;
double gast;
gmst = get_mean_sidereal_time(JD, JDE);
/* Get mean obliquity from precession (in arcseconds) */
get_precession_angles_vondrak(JDE, &epsilon_A, &chi_A, &omega_A, &psi_A);
/* Get nutation in longitude (in arcseconds) */
get_nutation_angles_iau2000b(JDE, &delta_psi, &delta_epsilon);
/*
* Equation of the equinoxes:
* delta_psi is in arcseconds, epsilon_A is in arcseconds.
* Convert delta_psi to radians, multiply by cos(epsilon_A in radians),
* result is in radians.
*/
eq_eq = (delta_psi * ARCSEC_TO_RAD) * cos(epsilon_A * ARCSEC_TO_RAD);
gast = gmst + eq_eq;
/* Normalize to [0, 2*pi) */
gast = fmod(gast, 2.0 * M_PI);
if (gast < 0.0)
gast += 2.0 * M_PI;
return gast;
}

47
src/sidereal_time.h Normal file
View File

@ -0,0 +1,47 @@
/*
* sidereal_time.h -- Greenwich sidereal time (mean and apparent)
*
* Clean-room implementation from published standards:
* IERS Conventions (2010), IERS Technical Note 36, Ch. 5.
* Capitaine, Guinot & McCarthy (2000), "Definition of the
* Celestial Ephemeris Origin and of UT1 in the International
* Celestial Reference Frame", A&A 355, 398-405.
*
* No PostgreSQL dependencies -- pure math, suitable for standalone use.
*/
#ifndef PG_ORBIT_SIDEREAL_TIME_H
#define PG_ORBIT_SIDEREAL_TIME_H
/*
* get_mean_sidereal_time -- Greenwich Mean Sidereal Time
*
* Computes GMST using the classical polynomial in UT1,
* consistent with IAU 2006 precession.
*
* JD = Julian Date (UT1 timescale)
* JDE = Julian Ephemeris Date (TDB/TT timescale)
*
* For pg_orbit Phase 2, JD and JDE are treated as identical
* (delta-T correction not yet applied; the difference is <1s
* for recent epochs).
*
* Returns GMST in radians, normalized to [0, 2*pi).
*/
double get_mean_sidereal_time(double JD, double JDE);
/*
* get_apparent_sidereal_time -- Greenwich Apparent Sidereal Time
*
* GAST = GMST + equation of the equinoxes.
*
* The equation of the equinoxes accounts for nutation:
* EqEq = delta_psi * cos(epsilon_A)
* where delta_psi is nutation in longitude and epsilon_A is
* the mean obliquity, both from precession.h.
*
* Returns GAST in radians, normalized to [0, 2*pi).
*/
double get_apparent_sidereal_time(double JD, double JDE);
#endif /* PG_ORBIT_SIDEREAL_TIME_H */

122
src/star_funcs.c Normal file
View File

@ -0,0 +1,122 @@
/*
* star_funcs.c -- Star and fixed-position object observation
*
* Takes J2000 catalog coordinates (RA in hours, Dec in degrees),
* applies IAU 1976 precession to date of observation, computes
* local hour angle, and converts to topocentric azimuth/elevation.
*
* Range and range_rate are zero -- stars are effectively at infinity.
* For objects with known proper motion, apply it to (RA, Dec) before
* calling star_observe.
*/
#include "postgres.h"
#include "fmgr.h"
#include "utils/timestamp.h"
#include "types.h"
#include "astro_math.h"
PG_FUNCTION_INFO_V1(star_observe);
PG_FUNCTION_INFO_V1(star_observe_safe);
/*
* star_observe(ra_hours, dec_degrees, observer, timestamptz) -> topocentric
*
* Compute az/el of a fixed celestial object from an observer at a time.
* Uses IAU 1976 precession (~1 arcsecond accuracy for centuries near J2000).
*/
Datum
star_observe(PG_FUNCTION_ARGS)
{
double ra_hours = PG_GETARG_FLOAT8(0);
double dec_deg = PG_GETARG_FLOAT8(1);
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2);
int64 ts = PG_GETARG_INT64(3);
double jd;
double ra_j2000, dec_j2000;
double ra_date, dec_date;
double gmst, lst, ha;
double az, el;
pg_topocentric *result;
if (ra_hours < 0.0 || ra_hours >= 24.0)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("right ascension out of range: %.6f", ra_hours),
errhint("RA must be in [0, 24) hours.")));
if (dec_deg < -90.0 || dec_deg > 90.0)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("declination out of range: %.6f", dec_deg),
errhint("Declination must be between -90 and +90 degrees.")));
jd = timestamptz_to_jd(ts);
ra_j2000 = ra_hours * (M_PI / 12.0);
dec_j2000 = dec_deg * DEG_TO_RAD;
precess_j2000_to_date(jd, ra_j2000, dec_j2000, &ra_date, &dec_date);
gmst = gmst_from_jd(jd);
lst = gmst + obs->lon;
ha = lst - ra_date;
equatorial_to_horizontal(ha, dec_date, obs->lat, &az, &el);
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
result->azimuth = az;
result->elevation = el;
result->range_km = 0.0;
result->range_rate = 0.0;
PG_RETURN_POINTER(result);
}
/*
* star_observe_safe -- returns NULL if inputs are out of range.
* For batch queries over star catalogs.
*/
Datum
star_observe_safe(PG_FUNCTION_ARGS)
{
double ra_hours = PG_GETARG_FLOAT8(0);
double dec_deg = PG_GETARG_FLOAT8(1);
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2);
int64 ts = PG_GETARG_INT64(3);
double jd;
double ra_j2000, dec_j2000;
double ra_date, dec_date;
double gmst, lst, ha;
double az, el;
pg_topocentric *result;
if (ra_hours < 0.0 || ra_hours >= 24.0)
PG_RETURN_NULL();
if (dec_deg < -90.0 || dec_deg > 90.0)
PG_RETURN_NULL();
jd = timestamptz_to_jd(ts);
ra_j2000 = ra_hours * (M_PI / 12.0);
dec_j2000 = dec_deg * DEG_TO_RAD;
precess_j2000_to_date(jd, ra_j2000, dec_j2000, &ra_date, &dec_date);
gmst = gmst_from_jd(jd);
lst = gmst + obs->lon;
ha = lst - ra_date;
equatorial_to_horizontal(ha, dec_date, obs->lat, &az, &el);
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
result->azimuth = az;
result->elevation = el;
result->range_km = 0.0;
result->range_rate = 0.0;
PG_RETURN_POINTER(result);
}

35
src/types.h
View File

@ -187,4 +187,39 @@ jd_to_minutes_since_epoch(double jd, double tle_epoch_jd)
return (jd - tle_epoch_jd) * 1440.0; return (jd - tle_epoch_jd) * 1440.0;
} }
/*
* Heliocentric position -- ecliptic J2000 frame
*
* Position in AU relative to the Sun.
* Used for planets (VSOP87), comets (Keplerian), and Lambert transfers.
*/
typedef struct pg_heliocentric
{
double x, y, z; /* AU, ecliptic J2000 */
} pg_heliocentric;
/*
* Astronomical constants
*/
#define AU_KM 149597870.7 /* km per AU (IAU 2012) */
#define GAUSS_K 0.01720209895 /* gravitational constant, AU^(3/2)/day */
#define GAUSS_K2 (GAUSS_K * GAUSS_K)
#define OBLIQUITY_J2000 0.40909280422232897 /* 23.4392911 deg in radians */
/*
* Solar system body IDs (VSOP87 convention, extended)
*/
#define BODY_SUN 0
#define BODY_MERCURY 1
#define BODY_VENUS 2
#define BODY_EARTH 3
#define BODY_MARS 4
#define BODY_JUPITER 5
#define BODY_SATURN 6
#define BODY_URANUS 7
#define BODY_NEPTUNE 8
#define BODY_MOON 10
#endif /* PG_ORBIT_TYPES_H */ #endif /* PG_ORBIT_TYPES_H */

137337
src/vsop87.c Normal file
View File

File diff suppressed because it is too large Load Diff

83
src/vsop87.h Normal file
View File

@ -0,0 +1,83 @@
/************************************************************************
Derived from Stellarium's VSOP87 implementation.
Modified for pg_orbit: removed static caching for thread safety
(PostgreSQL PARALLEL SAFE).
The PLANETARY SOLUTION VSOP87 by Bretagnon P. and Francou G. can be found at
ftp://ftp.imcce.fr/pub/ephem/planets/vsop87
I (Johannes Gajdosik) have just taken the data obtained from above
(VSOP87.mer,...,VSOP87.nep) and rearranged it into this piece of software.
I can neither allow nor forbid the usage of VSOP87.
The copyright notice below covers not the work of Bretagnon P. and Francou G.
but just my work, that is the compilation of the VSOP87 data
into the software supplied in this file.
Copyright (c) 2006 Johannes Gajdosik
Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the "Software"),
to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense,
and/or sell copies of the Software, and to permit persons to whom the
Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included
in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
This is my implementation of the VSOP87 planetary solution.
I tried to optimize for speed by rearranging the terms so that
for a given touple of (a[0],a[1],...,a[11]) the values
(cos,sin)(a[0]*lambda[0](T)+...a[11]*lambda[11](T))
have only to be calculated once.
Furthermore I used the addition formulas
(cos,sin)(x+y) = ...
so that for given T the functions cos and sin have only to be called 12 times.
Thread-safe: all functions are reentrant with no static mutable state.
****************************************************************/
#ifndef VSOP87_H
#define VSOP87_H
#ifdef __cplusplus
extern "C" {
#endif
void GetVsop87Coor(double jde,int body,double *xyz);
/* Return the rectangular coordinates of the given planet
and the given julian date jd expressed in dynamical time (TAI+32.184s).
The origin of the xyz-coordinates is the center of the sun.
The reference frame is "dynamical equinox and ecliptic J2000",
which is the reference frame in VSOP87 and VSOP87A.
body: 0=Mercury, 1=Venus, 2=Earth, 3=Mars,
4=Jupiter, 5=Saturn, 6=Uranus, 7=Neptune
xyz must point to an array of at least 6 doubles:
xyz[0..2] = position (AU)
xyz[3..5] = velocity (AU/day)
*/
void GetVsop87OsculatingCoor(const double jde0,const double jde, const int body,double *xyz);
/* The osculating orbit of epoch jde0, evaluated at jde, is returned.
*/
#ifdef __cplusplus
}
#endif
#endif

135
test/expected/kepler_comet.out Normal file
View File

@ -0,0 +1,135 @@
-- kepler_comet regression tests
--
-- Tests Keplerian propagation and comet observation.
-- Verifies Kepler equation solver for elliptic, parabolic,
-- and hyperbolic orbits.
\set boulder '''40.015N 105.270W 1655m'''::observer
-- ============================================================
-- Test 1: Circular orbit (e=0) at 1 AU
-- A body in a circular orbit at 1 AU with i=0, omega=0, Omega=0,
-- perihelion at J2000.0 should complete one orbit in ~365.25 days.
-- At T=0 (perihelion), position should be (1, 0, 0).
-- ============================================================
SELECT 'circular_t0' AS test,
round(helio_x(kepler_propagate(1.0, 0.0, 0.0, 0.0, 0.0,
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS x,
round(helio_y(kepler_propagate(1.0, 0.0, 0.0, 0.0, 0.0,
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS y,
round(helio_z(kepler_propagate(1.0, 0.0, 0.0, 0.0, 0.0,
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS z;
test | x | y | z
-------------+----------+----------+----------
circular_t0 | 1.000000 | 0.000000 | 0.000000
(1 row)
-- ============================================================
-- Test 2: Circular orbit at 1 AU, quarter orbit (~91.3 days)
-- Should be at approximately (0, 1, 0) for i=0 prograde orbit.
-- ============================================================
SELECT 'circular_quarter' AS test,
round(helio_x(kepler_propagate(1.0, 0.0, 0.0, 0.0, 0.0,
2451545.0, '2000-04-01 18:00:00+00'))::numeric, 1) AS x_approx,
round(helio_y(kepler_propagate(1.0, 0.0, 0.0, 0.0, 0.0,
2451545.0, '2000-04-01 18:00:00+00'))::numeric, 1) AS y_approx;
test | x_approx | y_approx
------------------+----------+----------
circular_quarter | 0.0 | 1.0
(1 row)
-- ============================================================
-- Test 3: Eccentric elliptic orbit (e=0.5)
-- q=0.5 AU, e=0.5 -> a=1.0 AU, same period as Earth.
-- At perihelion (t=0), position should be (0.5, 0, 0).
-- ============================================================
SELECT 'eccentric_t0' AS test,
round(helio_x(kepler_propagate(0.5, 0.5, 0.0, 0.0, 0.0,
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS x,
round(helio_y(kepler_propagate(0.5, 0.5, 0.0, 0.0, 0.0,
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS y;
test | x | y
--------------+----------+----------
eccentric_t0 | 0.500000 | 0.000000
(1 row)
-- ============================================================
-- Test 4: Inclined orbit (i=90 deg)
-- With i=90, the orbit is polar. At perihelion the position
-- is still in the orbital plane but rotated by inclination.
-- q=1.0 AU, e=0, i=90, omega=0, Omega=0 at t=0:
-- position should be (1, 0, 0) regardless of inclination
-- (perihelion direction is in the node line).
-- ============================================================
SELECT 'inclined_t0' AS test,
round(helio_x(kepler_propagate(1.0, 0.0, 90.0, 0.0, 0.0,
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS x,
round(helio_z(kepler_propagate(1.0, 0.0, 90.0, 0.0, 0.0,
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS z;
test | x | z
-------------+----------+----------
inclined_t0 | 1.000000 | 0.000000
(1 row)
-- ============================================================
-- Test 5: Heliocentric distance conservation for circular orbit
-- Distance should remain ~1.0 AU at all times.
-- ============================================================
SELECT 'distance_conservation' AS test,
round(helio_distance(kepler_propagate(1.0, 0.0, 0.0, 0.0, 0.0,
2451545.0, '2000-07-01 12:00:00+00'))::numeric, 6) AS dist_au;
test | dist_au
-----------------------+----------
distance_conservation | 1.000000
(1 row)
-- ============================================================
-- Test 6: Hyperbolic orbit (e=1.5)
-- q=1.0 AU, e=1.5, i=0, omega=0, Omega=0
-- At perihelion, position = (1, 0, 0). Should propagate without error.
-- ============================================================
SELECT 'hyperbolic_t0' AS test,
round(helio_x(kepler_propagate(1.0, 1.5, 0.0, 0.0, 0.0,
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS x,
round(helio_distance(kepler_propagate(1.0, 1.5, 0.0, 0.0, 0.0,
2451545.0, '2000-07-01 12:00:00+00'))::numeric, 2) AS dist_later;
test | x | dist_later
---------------+----------+------------
hyperbolic_t0 | 1.000000 | 3.54
(1 row)
-- ============================================================
-- Test 7: Near-parabolic orbit (e=1.0)
-- Comet on parabolic trajectory, q=1.0 AU.
-- At perihelion, position = (1, 0, 0).
-- ============================================================
SELECT 'parabolic_t0' AS test,
round(helio_x(kepler_propagate(1.0, 1.0, 0.0, 0.0, 0.0,
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS x;
test | x
--------------+----------
parabolic_t0 | 1.000000
(1 row)
-- ============================================================
-- Test 8: Comet observation with known Earth position
-- Place a "comet" at (2, 0, 0) AU heliocentric ecliptic.
-- Earth at (1, 0, 0) AU. Geocentric comet is at (1, 0, 0) AU,
-- which in equatorial J2000 is RA=0h, Dec=0 (on ecliptic plane
-- after rotation, dec ~ -obliquity... actually RA=0, Dec=0
-- because ecliptic-to-equatorial of (1,0,0) = (1,0,0)).
-- Geocentric distance should be 1 AU = 149597870.7 km.
-- ============================================================
SELECT 'comet_obs_range' AS test,
round(topo_range(comet_observe(
2.0, 0.0, 0.0, 0.0, 0.0, 2451545.0,
1.0, 0.0, 0.0,
:boulder, '2000-01-01 12:00:00+00'))::numeric, 0) AS range_km;
test | range_km
-----------------+-----------
comet_obs_range | 149597871
(1 row)
-- ============================================================
-- Test 9: Input validation
-- ============================================================
SELECT 'negative_q' AS test, kepler_propagate(-1.0, 0.5, 0.0, 0.0, 0.0, 2451545.0, now());
ERROR: perihelion distance must be positive: -1.000000

166
test/expected/planet_observe.out Normal file
View File

@ -0,0 +1,166 @@
-- planet_observe regression tests
--
-- Tests VSOP87 planetary positions, Sun and Moon observation,
-- and the planet_heliocentric function.
-- Reference values cross-checked against JPL Horizons.
\set boulder '''40.015N 105.270W 1655m'''::observer
-- ============================================================
-- Test 1: Sun at heliocentric origin
-- planet_heliocentric(0, t) should always return (0, 0, 0)
-- because the Sun IS the origin of heliocentric coordinates.
-- ============================================================
SELECT 'sun_origin' AS test,
round(helio_x(planet_heliocentric(0, '2024-06-21 12:00:00+00'))::numeric, 10) AS x,
round(helio_y(planet_heliocentric(0, '2024-06-21 12:00:00+00'))::numeric, 10) AS y,
round(helio_z(planet_heliocentric(0, '2024-06-21 12:00:00+00'))::numeric, 10) AS z;
test | x | y | z
------------+--------------+--------------+--------------
sun_origin | 0.0000000000 | 0.0000000000 | 0.0000000000
(1 row)
-- ============================================================
-- Test 2: Earth heliocentric distance ~ 1 AU
-- Earth (body_id=3) should be ~0.983-1.017 AU from the Sun
-- throughout the year (eccentricity ~0.017).
-- ============================================================
SELECT 'earth_distance' AS test,
round(helio_distance(planet_heliocentric(3, '2024-01-03 12:00:00+00'))::numeric, 3) AS perihelion,
round(helio_distance(planet_heliocentric(3, '2024-07-05 12:00:00+00'))::numeric, 3) AS aphelion;
test | perihelion | aphelion
----------------+------------+----------
earth_distance | 0.983 | 1.017
(1 row)
-- ============================================================
-- Test 3: Mars heliocentric distance ~ 1.38-1.67 AU
-- Mars (body_id=4) average distance is ~1.524 AU.
-- ============================================================
SELECT 'mars_distance' AS test,
round(helio_distance(planet_heliocentric(4, '2024-06-21 12:00:00+00'))::numeric, 2) AS dist_au;
test | dist_au
---------------+---------
mars_distance | 1.40
(1 row)
-- ============================================================
-- Test 4: Jupiter heliocentric distance ~ 4.95-5.46 AU
-- Jupiter (body_id=5) average distance is ~5.203 AU.
-- ============================================================
SELECT 'jupiter_distance' AS test,
round(helio_distance(planet_heliocentric(5, '2024-06-21 12:00:00+00'))::numeric, 1) AS dist_au;
test | dist_au
------------------+---------
jupiter_distance | 5.0
(1 row)
-- ============================================================
-- Test 5: planet_observe for Jupiter from Boulder
-- Jupiter should be observable (check we get a valid result,
-- not an error). Elevation will vary by time.
-- ============================================================
SELECT 'jupiter_observe' AS test,
round(topo_azimuth(planet_observe(5, :boulder,
'2024-03-15 03:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(planet_observe(5, :boulder,
'2024-03-15 03:00:00+00'))::numeric, 0) AS el_deg;
test | az_deg | el_deg
-----------------+--------+--------
jupiter_observe | 270 | 24
(1 row)
-- ============================================================
-- Test 6: planet_observe for Venus from Boulder
-- Venus (body_id=2) should return valid az/el.
-- ============================================================
SELECT 'venus_observe' AS test,
round(topo_azimuth(planet_observe(2, :boulder,
'2024-06-15 02:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(planet_observe(2, :boulder,
'2024-06-15 02:00:00+00'))::numeric, 0) AS el_deg;
test | az_deg | el_deg
---------------+--------+--------
venus_observe | 295 | 7
(1 row)
-- ============================================================
-- Test 7: sun_observe from Boulder at local noon (~18:00 UTC)
-- At summer solstice noon in Boulder (MDT = UTC-6),
-- Sun should be high in the south (az ~180, el ~73 deg).
-- ============================================================
SELECT 'sun_noon' AS test,
round(topo_azimuth(sun_observe(:boulder,
'2024-06-21 18:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(sun_observe(:boulder,
'2024-06-21 18:00:00+00'))::numeric, 0) AS el_deg;
test | az_deg | el_deg
----------+--------+--------
sun_noon | 137 | 69
(1 row)
-- ============================================================
-- Test 8: sun_observe range should be ~1 AU = 149,597,871 km
-- ============================================================
SELECT 'sun_range' AS test,
round(topo_range(sun_observe(:boulder,
'2024-06-21 18:00:00+00'))::numeric, -4) AS range_km;
test | range_km
-----------+-----------
sun_range | 152030000
(1 row)
-- ============================================================
-- Test 9: moon_observe from Boulder
-- Should return valid az/el with range ~ 356,000-407,000 km.
-- ============================================================
SELECT 'moon_observe' AS test,
round(topo_azimuth(moon_observe(:boulder,
'2024-06-21 18:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(moon_observe(:boulder,
'2024-06-21 18:00:00+00'))::numeric, 0) AS el_deg,
round(topo_range(moon_observe(:boulder,
'2024-06-21 18:00:00+00'))::numeric, -3) AS range_km;
test | az_deg | el_deg | range_km
--------------+--------+--------+----------
moon_observe | 317 | -75 | 381000
(1 row)
-- ============================================================
-- Test 10: Moon range sanity (should be 356k-407k km)
-- ============================================================
SELECT 'moon_range_check' AS test,
topo_range(moon_observe(:boulder, '2024-06-21 18:00:00+00')) BETWEEN 350000 AND 410000 AS in_range;
test | in_range
------------------+----------
moon_range_check | t
(1 row)
-- ============================================================
-- Test 11: Error handling - cannot observe Earth from Earth
-- ============================================================
SELECT 'earth_error' AS test, planet_observe(3, :boulder, now());
ERROR: cannot observe Earth from Earth
-- ============================================================
-- Test 12: Error handling - invalid body_id
-- ============================================================
SELECT 'invalid_body' AS test, planet_heliocentric(9, now());
ERROR: invalid body_id 9: must be 0 (Sun) or 1-8 (Mercury-Neptune)
-- ============================================================
-- Test 13: All planets return valid heliocentric positions
-- Mercury through Neptune, none should error.
-- ============================================================
SELECT 'all_planets' AS test,
body_id,
round(helio_distance(planet_heliocentric(body_id, '2024-06-21 12:00:00+00'))::numeric, 2) AS dist_au
FROM generate_series(1, 8) AS body_id;
test | body_id | dist_au
-------------+---------+---------
all_planets | 1 | 0.33
all_planets | 2 | 0.72
all_planets | 3 | 1.02
all_planets | 4 | 1.40
all_planets | 5 | 5.02
all_planets | 6 | 9.69
all_planets | 7 | 19.59
all_planets | 8 | 29.90
(8 rows)

105
test/expected/star_observe.out Normal file
View File

@ -0,0 +1,105 @@
-- star_observe regression tests
--
-- Tests IAU 1976 precession and equatorial-to-horizontal transform.
-- Reference time: J2000.0 (2000-01-01 12:00:00 UTC) simplifies
-- precession checks (should be identity at epoch).
-- Setup: common observer location (Boulder, CO)
\set boulder '''40.015N 105.270W 1655m'''::observer
-- ============================================================
-- Test 1: Star at J2000.0 epoch (precession should be ~identity)
-- Polaris: RA 2h 31m 49.09s = 2.530303h, Dec +89.2641 deg
-- At J2000.0 from Boulder, should be at high elevation (~49 deg)
-- because observer lat ~40N and Polaris dec ~89.26 deg
-- ============================================================
SELECT 'polaris_j2000' AS test,
round(topo_azimuth(star_observe(2.530303, 89.2641, :boulder,
'2000-01-01 12:00:00+00'))::numeric, 1) AS az_deg,
round(topo_elevation(star_observe(2.530303, 89.2641, :boulder,
'2000-01-01 12:00:00+00'))::numeric, 1) AS el_deg;
test | az_deg | el_deg
---------------+--------+--------
polaris_j2000 | 359.4 | 39.5
(1 row)
-- ============================================================
-- Test 2: Star on the celestial equator at RA=LST should transit
-- at elevation = (90 - observer_lat) from the south.
-- For Boulder (lat ~40): transit elevation ~50 deg, az ~180
-- We pick RA such that it's on the meridian at our test time.
-- ============================================================
-- Test 3: Polaris from different times (precession moves it)
-- At 2025-06-15, precession has shifted by ~0.35 deg over 25 years.
-- Elevation should change by a fraction of a degree.
SELECT 'polaris_2025' AS test,
round(topo_elevation(star_observe(2.530303, 89.2641, :boulder,
'2025-06-15 04:00:00+00'))::numeric, 1) AS el_deg;
test | el_deg
--------------+--------
polaris_2025 | 39.4
(1 row)
-- ============================================================
-- Test 4: Sirius (brightest star)
-- RA 6h 45m 08.92s = 6.752478h, Dec -16.7161 deg
-- From Boulder at 2024-01-15 03:00 UTC (evening, winter)
-- Should be visible (positive elevation)
-- ============================================================
SELECT 'sirius_winter' AS test,
round(topo_azimuth(star_observe(6.752478, -16.7161, :boulder,
'2024-01-15 03:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(star_observe(6.752478, -16.7161, :boulder,
'2024-01-15 03:00:00+00'))::numeric, 0) AS el_deg;
test | az_deg | el_deg
---------------+--------+--------
sirius_winter | 132 | 18
(1 row)
-- ============================================================
-- Test 5: Vega
-- RA 18h 36m 56.34s = 18.615650h, Dec +38.7837 deg
-- ============================================================
SELECT 'vega' AS test,
round(topo_elevation(star_observe(18.615650, 38.7837, :boulder,
'2024-07-15 04:00:00+00'))::numeric, 0) AS el_deg;
test | el_deg
------+--------
vega | 66
(1 row)
-- ============================================================
-- Test 6: star_observe_safe with invalid RA returns NULL
-- ============================================================
SELECT 'safe_null' AS test,
star_observe_safe(25.0, 45.0, :boulder, '2024-01-01 00:00:00+00') IS NULL AS is_null;
test | is_null
-----------+---------
safe_null | t
(1 row)
-- ============================================================
-- Test 7: Heliocentric type I/O round-trip
-- ============================================================
SELECT 'helio_io' AS test,
'(1.0000000000,0.0000000000,0.0000000000)'::heliocentric::text AS val;
test | val
----------+------------------------------------------
helio_io | (1.0000000000,0.0000000000,0.0000000000)
(1 row)
SELECT 'helio_accessors' AS test,
round(helio_x(h)::numeric, 4) AS x,
round(helio_y(h)::numeric, 4) AS y,
round(helio_z(h)::numeric, 4) AS z,
round(helio_distance(h)::numeric, 4) AS dist
FROM (SELECT '(0.3000000000,0.4000000000,0.0000000000)'::heliocentric AS h) sub;
test | x | y | z | dist
-----------------+--------+--------+--------+--------
helio_accessors | 0.3000 | 0.4000 | 0.0000 | 0.5000
(1 row)
-- ============================================================
-- Test 8: star_observe with error input
-- ============================================================
SELECT 'error_ra' AS test, star_observe(30.0, 0.0, :boulder, now());
ERROR: right ascension out of range: 30.000000
HINT: RA must be in [0, 24) hours.

104
test/sql/kepler_comet.sql Normal file
View File

@ -0,0 +1,104 @@
-- kepler_comet regression tests
--
-- Tests Keplerian propagation and comet observation.
-- Verifies Kepler equation solver for elliptic, parabolic,
-- and hyperbolic orbits.
\set boulder '''40.015N 105.270W 1655m'''::observer
-- ============================================================
-- Test 1: Circular orbit (e=0) at 1 AU
-- A body in a circular orbit at 1 AU with i=0, omega=0, Omega=0,
-- perihelion at J2000.0 should complete one orbit in ~365.25 days.
-- At T=0 (perihelion), position should be (1, 0, 0).
-- ============================================================
SELECT 'circular_t0' AS test,
round(helio_x(kepler_propagate(1.0, 0.0, 0.0, 0.0, 0.0,
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS x,
round(helio_y(kepler_propagate(1.0, 0.0, 0.0, 0.0, 0.0,
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS y,
round(helio_z(kepler_propagate(1.0, 0.0, 0.0, 0.0, 0.0,
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS z;
-- ============================================================
-- Test 2: Circular orbit at 1 AU, quarter orbit (~91.3 days)
-- Should be at approximately (0, 1, 0) for i=0 prograde orbit.
-- ============================================================
SELECT 'circular_quarter' AS test,
round(helio_x(kepler_propagate(1.0, 0.0, 0.0, 0.0, 0.0,
2451545.0, '2000-04-01 18:00:00+00'))::numeric, 1) AS x_approx,
round(helio_y(kepler_propagate(1.0, 0.0, 0.0, 0.0, 0.0,
2451545.0, '2000-04-01 18:00:00+00'))::numeric, 1) AS y_approx;
-- ============================================================
-- Test 3: Eccentric elliptic orbit (e=0.5)
-- q=0.5 AU, e=0.5 -> a=1.0 AU, same period as Earth.
-- At perihelion (t=0), position should be (0.5, 0, 0).
-- ============================================================
SELECT 'eccentric_t0' AS test,
round(helio_x(kepler_propagate(0.5, 0.5, 0.0, 0.0, 0.0,
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS x,
round(helio_y(kepler_propagate(0.5, 0.5, 0.0, 0.0, 0.0,
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS y;
-- ============================================================
-- Test 4: Inclined orbit (i=90 deg)
-- With i=90, the orbit is polar. At perihelion the position
-- is still in the orbital plane but rotated by inclination.
-- q=1.0 AU, e=0, i=90, omega=0, Omega=0 at t=0:
-- position should be (1, 0, 0) regardless of inclination
-- (perihelion direction is in the node line).
-- ============================================================
SELECT 'inclined_t0' AS test,
round(helio_x(kepler_propagate(1.0, 0.0, 90.0, 0.0, 0.0,
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS x,
round(helio_z(kepler_propagate(1.0, 0.0, 90.0, 0.0, 0.0,
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS z;
-- ============================================================
-- Test 5: Heliocentric distance conservation for circular orbit
-- Distance should remain ~1.0 AU at all times.
-- ============================================================
SELECT 'distance_conservation' AS test,
round(helio_distance(kepler_propagate(1.0, 0.0, 0.0, 0.0, 0.0,
2451545.0, '2000-07-01 12:00:00+00'))::numeric, 6) AS dist_au;
-- ============================================================
-- Test 6: Hyperbolic orbit (e=1.5)
-- q=1.0 AU, e=1.5, i=0, omega=0, Omega=0
-- At perihelion, position = (1, 0, 0). Should propagate without error.
-- ============================================================
SELECT 'hyperbolic_t0' AS test,
round(helio_x(kepler_propagate(1.0, 1.5, 0.0, 0.0, 0.0,
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS x,
round(helio_distance(kepler_propagate(1.0, 1.5, 0.0, 0.0, 0.0,
2451545.0, '2000-07-01 12:00:00+00'))::numeric, 2) AS dist_later;
-- ============================================================
-- Test 7: Near-parabolic orbit (e=1.0)
-- Comet on parabolic trajectory, q=1.0 AU.
-- At perihelion, position = (1, 0, 0).
-- ============================================================
SELECT 'parabolic_t0' AS test,
round(helio_x(kepler_propagate(1.0, 1.0, 0.0, 0.0, 0.0,
2451545.0, '2000-01-01 12:00:00+00'))::numeric, 6) AS x;
-- ============================================================
-- Test 8: Comet observation with known Earth position
-- Place a "comet" at (2, 0, 0) AU heliocentric ecliptic.
-- Earth at (1, 0, 0) AU. Geocentric comet is at (1, 0, 0) AU,
-- which in equatorial J2000 is RA=0h, Dec=0 (on ecliptic plane
-- after rotation, dec ~ -obliquity... actually RA=0, Dec=0
-- because ecliptic-to-equatorial of (1,0,0) = (1,0,0)).
-- Geocentric distance should be 1 AU = 149597870.7 km.
-- ============================================================
SELECT 'comet_obs_range' AS test,
round(topo_range(comet_observe(
2.0, 0.0, 0.0, 0.0, 0.0, 2451545.0,
1.0, 0.0, 0.0,
:boulder, '2000-01-01 12:00:00+00'))::numeric, 0) AS range_km;
-- ============================================================
-- Test 9: Input validation
-- ============================================================
SELECT 'negative_q' AS test, kepler_propagate(-1.0, 0.5, 0.0, 0.0, 0.0, 2451545.0, now());

116
test/sql/planet_observe.sql Normal file
View File

@ -0,0 +1,116 @@
-- planet_observe regression tests
--
-- Tests VSOP87 planetary positions, Sun and Moon observation,
-- and the planet_heliocentric function.
-- Reference values cross-checked against JPL Horizons.
\set boulder '''40.015N 105.270W 1655m'''::observer
-- ============================================================
-- Test 1: Sun at heliocentric origin
-- planet_heliocentric(0, t) should always return (0, 0, 0)
-- because the Sun IS the origin of heliocentric coordinates.
-- ============================================================
SELECT 'sun_origin' AS test,
round(helio_x(planet_heliocentric(0, '2024-06-21 12:00:00+00'))::numeric, 10) AS x,
round(helio_y(planet_heliocentric(0, '2024-06-21 12:00:00+00'))::numeric, 10) AS y,
round(helio_z(planet_heliocentric(0, '2024-06-21 12:00:00+00'))::numeric, 10) AS z;
-- ============================================================
-- Test 2: Earth heliocentric distance ~ 1 AU
-- Earth (body_id=3) should be ~0.983-1.017 AU from the Sun
-- throughout the year (eccentricity ~0.017).
-- ============================================================
SELECT 'earth_distance' AS test,
round(helio_distance(planet_heliocentric(3, '2024-01-03 12:00:00+00'))::numeric, 3) AS perihelion,
round(helio_distance(planet_heliocentric(3, '2024-07-05 12:00:00+00'))::numeric, 3) AS aphelion;
-- ============================================================
-- Test 3: Mars heliocentric distance ~ 1.38-1.67 AU
-- Mars (body_id=4) average distance is ~1.524 AU.
-- ============================================================
SELECT 'mars_distance' AS test,
round(helio_distance(planet_heliocentric(4, '2024-06-21 12:00:00+00'))::numeric, 2) AS dist_au;
-- ============================================================
-- Test 4: Jupiter heliocentric distance ~ 4.95-5.46 AU
-- Jupiter (body_id=5) average distance is ~5.203 AU.
-- ============================================================
SELECT 'jupiter_distance' AS test,
round(helio_distance(planet_heliocentric(5, '2024-06-21 12:00:00+00'))::numeric, 1) AS dist_au;
-- ============================================================
-- Test 5: planet_observe for Jupiter from Boulder
-- Jupiter should be observable (check we get a valid result,
-- not an error). Elevation will vary by time.
-- ============================================================
SELECT 'jupiter_observe' AS test,
round(topo_azimuth(planet_observe(5, :boulder,
'2024-03-15 03:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(planet_observe(5, :boulder,
'2024-03-15 03:00:00+00'))::numeric, 0) AS el_deg;
-- ============================================================
-- Test 6: planet_observe for Venus from Boulder
-- Venus (body_id=2) should return valid az/el.
-- ============================================================
SELECT 'venus_observe' AS test,
round(topo_azimuth(planet_observe(2, :boulder,
'2024-06-15 02:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(planet_observe(2, :boulder,
'2024-06-15 02:00:00+00'))::numeric, 0) AS el_deg;
-- ============================================================
-- Test 7: sun_observe from Boulder at local noon (~18:00 UTC)
-- At summer solstice noon in Boulder (MDT = UTC-6),
-- Sun should be high in the south (az ~180, el ~73 deg).
-- ============================================================
SELECT 'sun_noon' AS test,
round(topo_azimuth(sun_observe(:boulder,
'2024-06-21 18:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(sun_observe(:boulder,
'2024-06-21 18:00:00+00'))::numeric, 0) AS el_deg;
-- ============================================================
-- Test 8: sun_observe range should be ~1 AU = 149,597,871 km
-- ============================================================
SELECT 'sun_range' AS test,
round(topo_range(sun_observe(:boulder,
'2024-06-21 18:00:00+00'))::numeric, -4) AS range_km;
-- ============================================================
-- Test 9: moon_observe from Boulder
-- Should return valid az/el with range ~ 356,000-407,000 km.
-- ============================================================
SELECT 'moon_observe' AS test,
round(topo_azimuth(moon_observe(:boulder,
'2024-06-21 18:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(moon_observe(:boulder,
'2024-06-21 18:00:00+00'))::numeric, 0) AS el_deg,
round(topo_range(moon_observe(:boulder,
'2024-06-21 18:00:00+00'))::numeric, -3) AS range_km;
-- ============================================================
-- Test 10: Moon range sanity (should be 356k-407k km)
-- ============================================================
SELECT 'moon_range_check' AS test,
topo_range(moon_observe(:boulder, '2024-06-21 18:00:00+00')) BETWEEN 350000 AND 410000 AS in_range;
-- ============================================================
-- Test 11: Error handling - cannot observe Earth from Earth
-- ============================================================
SELECT 'earth_error' AS test, planet_observe(3, :boulder, now());
-- ============================================================
-- Test 12: Error handling - invalid body_id
-- ============================================================
SELECT 'invalid_body' AS test, planet_heliocentric(9, now());
-- ============================================================
-- Test 13: All planets return valid heliocentric positions
-- Mercury through Neptune, none should error.
-- ============================================================
SELECT 'all_planets' AS test,
body_id,
round(helio_distance(planet_heliocentric(body_id, '2024-06-21 12:00:00+00'))::numeric, 2) AS dist_au
FROM generate_series(1, 8) AS body_id;

78
test/sql/star_observe.sql Normal file
View File

@ -0,0 +1,78 @@
-- star_observe regression tests
--
-- Tests IAU 1976 precession and equatorial-to-horizontal transform.
-- Reference time: J2000.0 (2000-01-01 12:00:00 UTC) simplifies
-- precession checks (should be identity at epoch).
-- Setup: common observer location (Boulder, CO)
\set boulder '''40.015N 105.270W 1655m'''::observer
-- ============================================================
-- Test 1: Star at J2000.0 epoch (precession should be ~identity)
-- Polaris: RA 2h 31m 49.09s = 2.530303h, Dec +89.2641 deg
-- At J2000.0 from Boulder, should be at high elevation (~49 deg)
-- because observer lat ~40N and Polaris dec ~89.26 deg
-- ============================================================
SELECT 'polaris_j2000' AS test,
round(topo_azimuth(star_observe(2.530303, 89.2641, :boulder,
'2000-01-01 12:00:00+00'))::numeric, 1) AS az_deg,
round(topo_elevation(star_observe(2.530303, 89.2641, :boulder,
'2000-01-01 12:00:00+00'))::numeric, 1) AS el_deg;
-- ============================================================
-- Test 2: Star on the celestial equator at RA=LST should transit
-- at elevation = (90 - observer_lat) from the south.
-- For Boulder (lat ~40): transit elevation ~50 deg, az ~180
-- We pick RA such that it's on the meridian at our test time.
-- ============================================================
-- Test 3: Polaris from different times (precession moves it)
-- At 2025-06-15, precession has shifted by ~0.35 deg over 25 years.
-- Elevation should change by a fraction of a degree.
SELECT 'polaris_2025' AS test,
round(topo_elevation(star_observe(2.530303, 89.2641, :boulder,
'2025-06-15 04:00:00+00'))::numeric, 1) AS el_deg;
-- ============================================================
-- Test 4: Sirius (brightest star)
-- RA 6h 45m 08.92s = 6.752478h, Dec -16.7161 deg
-- From Boulder at 2024-01-15 03:00 UTC (evening, winter)
-- Should be visible (positive elevation)
-- ============================================================
SELECT 'sirius_winter' AS test,
round(topo_azimuth(star_observe(6.752478, -16.7161, :boulder,
'2024-01-15 03:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(star_observe(6.752478, -16.7161, :boulder,
'2024-01-15 03:00:00+00'))::numeric, 0) AS el_deg;
-- ============================================================
-- Test 5: Vega
-- RA 18h 36m 56.34s = 18.615650h, Dec +38.7837 deg
-- ============================================================
SELECT 'vega' AS test,
round(topo_elevation(star_observe(18.615650, 38.7837, :boulder,
'2024-07-15 04:00:00+00'))::numeric, 0) AS el_deg;
-- ============================================================
-- Test 6: star_observe_safe with invalid RA returns NULL
-- ============================================================
SELECT 'safe_null' AS test,
star_observe_safe(25.0, 45.0, :boulder, '2024-01-01 00:00:00+00') IS NULL AS is_null;
-- ============================================================
-- Test 7: Heliocentric type I/O round-trip
-- ============================================================
SELECT 'helio_io' AS test,
'(1.0000000000,0.0000000000,0.0000000000)'::heliocentric::text AS val;
SELECT 'helio_accessors' AS test,
round(helio_x(h)::numeric, 4) AS x,
round(helio_y(h)::numeric, 4) AS y,
round(helio_z(h)::numeric, 4) AS z,
round(helio_distance(h)::numeric, 4) AS dist
FROM (SELECT '(0.3000000000,0.4000000000,0.0000000000)'::heliocentric AS h) sub;
-- ============================================================
-- Test 8: star_observe with error input
-- ============================================================
SELECT 'error_ra' AS test, star_observe(30.0, 0.0, :boulder, now());