Merge feature/lagrange-pg: v0.20.0 Lagrange point SQL functions
37 new SQL objects (188 → 225 total): Sun-planet, Earth-Moon, planetary moon L1-L5, Hill radius, DE variants, regression tests. All 31 tests pass.
This commit is contained in:
commit
a5123295d7
1
.gitignore
vendored
1
.gitignore
vendored
@ -24,6 +24,7 @@ test/test_de_reader
|
||||
test/test_od_math
|
||||
test/test_od_iod
|
||||
test/test_od_gauss
|
||||
test/test_lagrange
|
||||
|
||||
# Bench — downloaded TLE catalogs (large, ephemeral)
|
||||
# Already-tracked files (active.tle, spacetrack_full*.tle) are unaffected.
|
||||
|
||||
24
CLAUDE.md
24
CLAUDE.md
@ -1,9 +1,9 @@
|
||||
# pg_orrery — A Database Orrery for PostgreSQL
|
||||
|
||||
## What This Is
|
||||
A database orrery — celestial mechanics types and functions for PostgreSQL. Native C extension using PGXS, 188 SQL objects (172 user-visible functions + 16 GiST support), 9 custom types, covering satellites (SGP4/SDP4), planets (VSOP87 + optional JPL DE441), Moon (ELP2000-82B), 19 planetary moons (L1.2/TASS17/GUST86/MarsSat), stars (with proper motion and annual parallax), comets, asteroids (MPC catalog), Jupiter radio bursts, interplanetary Lambert transfers, equatorial RA/Dec coordinates with GiST-indexed angular separation, atmospheric refraction, annual stellar aberration, light-time correction, rise/set prediction (geometric + refracted + event windows + sun almanac) with status diagnostics, IAU constellation identification with full name lookup (Roman 1987), twilight dawn/dusk (civil/nautical/astronomical), lunar phase (angle, illumination, name, age), planet apparent magnitude with Saturn ring correction (Mallama & Hilton 2018), solar elongation, planet phase fraction, conjunction detection, satellite eclipse prediction (conical shadow with penumbral fraction), observing night quality assessment, lunar libration (optical + physical, Meeus Ch. 53 + p. 373), and angular separation rate.
|
||||
A database orrery — celestial mechanics types and functions for PostgreSQL. Native C extension using PGXS, 225 SQL objects (209 user-visible functions + 16 GiST support), 9 custom types, covering satellites (SGP4/SDP4), planets (VSOP87 + optional JPL DE441), Moon (ELP2000-82B), 19 planetary moons (L1.2/TASS17/GUST86/MarsSat), stars (with proper motion and annual parallax), comets, asteroids (MPC catalog), Jupiter radio bursts, interplanetary Lambert transfers, equatorial RA/Dec coordinates with GiST-indexed angular separation, atmospheric refraction, annual stellar aberration, light-time correction, rise/set prediction (geometric + refracted + event windows + sun almanac) with status diagnostics, IAU constellation identification with full name lookup (Roman 1987), twilight dawn/dusk (civil/nautical/astronomical), lunar phase (angle, illumination, name, age), planet apparent magnitude with Saturn ring correction (Mallama & Hilton 2018), solar elongation, planet phase fraction, conjunction detection, satellite eclipse prediction (conical shadow with penumbral fraction), observing night quality assessment, lunar libration (optical + physical, Meeus Ch. 53 + p. 373), angular separation rate, and Lagrange point equilibrium positions (CR3BP L1-L5 for Sun-planet, Earth-Moon, and 19 planetary moon systems).
|
||||
|
||||
**Current version:** 0.19.0
|
||||
**Current version:** 0.20.0
|
||||
**Repository:** https://git.supported.systems/warehack.ing/pg_orrery
|
||||
**Documentation:** https://pg-orrery.warehack.ing
|
||||
|
||||
@ -11,7 +11,7 @@ A database orrery — celestial mechanics types and functions for PostgreSQL. Na
|
||||
```bash
|
||||
make PG_CONFIG=/usr/bin/pg_config # Compile with PGXS
|
||||
sudo make install PG_CONFIG=/usr/bin/pg_config # Install extension
|
||||
make installcheck PG_CONFIG=/usr/bin/pg_config # Run 30 regression test suites
|
||||
make installcheck PG_CONFIG=/usr/bin/pg_config # Run 31 regression test suites
|
||||
```
|
||||
|
||||
Requires: PostgreSQL 17 development headers, GCC, Make.
|
||||
@ -27,7 +27,7 @@ Image: `git.supported.systems/warehack.ing/pg_orrery:pg17`
|
||||
|
||||
## Project Layout
|
||||
```
|
||||
pg_orrery.control # Extension metadata (version 0.19.0)
|
||||
pg_orrery.control # Extension metadata (version 0.20.0)
|
||||
Makefile # PGXS build + Docker targets
|
||||
sql/
|
||||
pg_orrery--0.1.0.sql # v0.1.0: satellite types/functions/operators
|
||||
@ -49,6 +49,7 @@ sql/
|
||||
pg_orrery--0.17.0.sql # v0.17.0: elongation, phase, eclipse, night quality, libration (174 objects)
|
||||
pg_orrery--0.18.0.sql # v0.18.0: ring tilt, penumbral eclipse, rise/set windows, angular rate (184 objects)
|
||||
pg_orrery--0.19.0.sql # v0.19.0: sun almanac, conjunctions, penumbral fraction, physical libration (188 objects)
|
||||
pg_orrery--0.20.0.sql # v0.20.0: Lagrange point equilibrium positions (225 objects)
|
||||
pg_orrery--0.1.0--0.2.0.sql # Migration: v0.1.0 → v0.2.0 (adds solar system)
|
||||
pg_orrery--0.2.0--0.3.0.sql # Migration: v0.2.0 → v0.3.0 (adds DE ephemeris)
|
||||
pg_orrery--0.3.0--0.4.0.sql # Migration: v0.3.0 → v0.4.0
|
||||
@ -67,6 +68,7 @@ sql/
|
||||
pg_orrery--0.16.0--0.17.0.sql # Migration: v0.16.0 → v0.17.0 (elongation, phase, eclipse, night quality, libration)
|
||||
pg_orrery--0.17.0--0.18.0.sql # Migration: v0.17.0 → v0.18.0 (ring tilt, penumbral eclipse, rise/set windows, angular rate)
|
||||
pg_orrery--0.18.0--0.19.0.sql # Migration: v0.18.0 → v0.19.0 (sun almanac, conjunctions, penumbral fraction, physical libration)
|
||||
pg_orrery--0.19.0--0.20.0.sql # Migration: v0.19.0 → v0.20.0 (Lagrange points)
|
||||
src/
|
||||
pg_orrery.c # PG_MODULE_MAGIC + _PG_init() (GUC registration)
|
||||
types.h # All struct definitions + constants + DE body ID mapping
|
||||
@ -100,6 +102,8 @@ src/
|
||||
magnitude_funcs.c # planet_magnitude() (with Saturn ring correction), solar_elongation(), planet_phase(), saturn_ring_tilt()
|
||||
eclipse_funcs.c # satellite eclipse prediction (conical shadow with penumbral fraction, Vallado §5.3)
|
||||
libration.h / libration_funcs.c # lunar libration (optical Meeus Ch. 53 + physical p. 373)
|
||||
lagrange.h # CR3BP solver (header-only): quintic solver, co-rotating frame, Hill radius
|
||||
lagrange_funcs.c # Lagrange point SQL functions (Sun-planet, Earth-Moon, planetary moons)
|
||||
l12.c / l12.h # L1.2 Galilean moon theory (Lieske 1998)
|
||||
tass17.c / tass17.h # TASS 1.7 Saturn moon theory (Vienne & Duriez 1995)
|
||||
gust86.c / gust86.h # GUST86 Uranus moon theory (Laskar & Jacobson 1987)
|
||||
@ -124,7 +128,7 @@ src/
|
||||
PROVENANCE.md # Vendoring decision, modifications, verification
|
||||
LICENSE # MIT license (Bill Gray / Project Pluto)
|
||||
test/
|
||||
sql/ # 30 regression test suites
|
||||
sql/ # 31 regression test suites
|
||||
expected/ # Expected output
|
||||
data/vallado_518.json # 518 Vallado test vectors (AIAA 2006-6753-Rev1)
|
||||
docs/
|
||||
@ -151,7 +155,7 @@ All types are fixed-size, `STORAGE = plain`, `ALIGNMENT = double`. No TOAST over
|
||||
| `orbital_elements` | 72 | Classical Keplerian elements for comets/asteroids (epoch, q, e, inc, omega, Omega, tp, H, G) |
|
||||
| `equatorial` | 24 | Apparent RA (hours), Dec (degrees), distance (km) — of date |
|
||||
|
||||
## Function Domains (188 SQL objects)
|
||||
## Function Domains (225 SQL objects)
|
||||
|
||||
| Domain | Theory | Key Functions | Count |
|
||||
|--------|--------|---------------|-------|
|
||||
@ -179,6 +183,7 @@ All types are fixed-size, `STORAGE = plain`, `ALIGNMENT = double`. No TOAST over
|
||||
| Observing quality | Composite (twilight+Moon) | `observing_night_quality()` | 1 |
|
||||
| Lunar libration | Meeus (1998) Ch. 53 + p. 373 | `moon_libration_longitude()`, `moon_libration()`, `moon_subsolar_longitude()`, `moon_physical_libration()` | 6 |
|
||||
| Constellation | Roman (1987) CDS VI/42 | `constellation()`, `constellation_full_name()` | 3 |
|
||||
| Lagrange points | CR3BP quintic + VSOP87 | `lagrange_heliocentric()`, `lagrange_observe()`, `hill_radius()` | 37 |
|
||||
| Diagnostics | -- | `pg_orrery_ephemeris_info()` | 1 |
|
||||
|
||||
All functions are `PARALLEL SAFE`. VSOP87/ELP82B functions are `IMMUTABLE` (compiled-in coefficients). DE functions are `STABLE` (external file dependency).
|
||||
@ -312,7 +317,7 @@ All numerical logic is byte-identical to upstream. Verified against 518 Vallado
|
||||
|
||||
## Testing
|
||||
|
||||
30 regression test suites via `make installcheck`:
|
||||
31 regression test suites via `make installcheck`:
|
||||
|
||||
| Suite | What it tests |
|
||||
|-------|--------------|
|
||||
@ -346,10 +351,11 @@ All numerical logic is byte-identical to upstream. Verified against 518 Vallado
|
||||
| v017_features | Solar elongation ranges/errors, planet phase ranges, satellite eclipse, observing night quality, lunar libration ranges, subsolar longitude |
|
||||
| v018_features | Saturn ring tilt range/variation, penumbral eclipse (shadow state, penumbra precedes umbra), rise/set event windows (Sun/Moon/planet, refracted vs geometric), angular separation rate (generic + planet convenience) |
|
||||
| v019_features | Sun almanac events (count/order/types/polar/refraction/window guard), conjunction detection (Jupiter-Saturn 2020, Moon-Venus, same-body error, threshold filter), penumbral fraction (range/bounds/eclipse consistency), physical libration (small corrections, time variation, total libration range) |
|
||||
| v020_features | Lagrange L1-L5 heliocentric/observe/equatorial, Hill radius, zone radius, mass ratio, DE fallback, all planet + moon families, input validation |
|
||||
|
||||
### PG Version Matrix
|
||||
|
||||
Test all 30 regression suites + DE reader unit test across PostgreSQL 14-18 using Docker:
|
||||
Test all 31 regression suites + DE reader unit test across PostgreSQL 14-18 using Docker:
|
||||
|
||||
```bash
|
||||
make test-matrix # Full matrix (PG 14-18)
|
||||
@ -375,7 +381,7 @@ Logs saved to `test/matrix-logs/pg${ver}.log`. The script reuses the Dockerfile
|
||||
|
||||
Starlight docs at `docs/` — 44+ MDX pages covering all domains.
|
||||
|
||||
Sections: Getting Started, Guides (9 domain walkthroughs incl. DE ephemeris), Workflow Translation (Skyfield/Horizons/GMAT/Radio Jupiter Pro comparisons), Reference (all 188 SQL objects incl. DE variants, equatorial GiST, refraction, rise/set, sun almanac, conjunction detection, constellation, twilight, lunar phase, planet magnitude, Saturn ring tilt, solar elongation, planet phase, satellite eclipse with penumbral fraction, observing quality, lunar libration with physical corrections, angular separation rate), Architecture (Hamilton's principles, constant custody, observation pipeline), Performance (benchmarks).
|
||||
Sections: Getting Started, Guides (9 domain walkthroughs incl. DE ephemeris), Workflow Translation (Skyfield/Horizons/GMAT/Radio Jupiter Pro comparisons), Reference (all 225 SQL objects incl. DE variants, equatorial GiST, refraction, rise/set, sun almanac, conjunction detection, constellation, Lagrange points, twilight, lunar phase, planet magnitude, Saturn ring tilt, solar elongation, planet phase, satellite eclipse with penumbral fraction, observing quality, lunar libration with physical corrections, angular separation rate), Architecture (Hamilton's principles, constant custody, observation pipeline), Performance (benchmarks).
|
||||
|
||||
### Local Development
|
||||
```bash
|
||||
|
||||
9
Makefile
9
Makefile
@ -17,7 +17,8 @@ DATA = sql/pg_orrery--0.1.0.sql sql/pg_orrery--0.2.0.sql sql/pg_orrery--0.1.0--0
|
||||
sql/pg_orrery--0.16.0.sql sql/pg_orrery--0.15.0--0.16.0.sql \
|
||||
sql/pg_orrery--0.17.0.sql sql/pg_orrery--0.16.0--0.17.0.sql \
|
||||
sql/pg_orrery--0.18.0.sql sql/pg_orrery--0.17.0--0.18.0.sql \
|
||||
sql/pg_orrery--0.19.0.sql sql/pg_orrery--0.18.0--0.19.0.sql
|
||||
sql/pg_orrery--0.19.0.sql sql/pg_orrery--0.18.0--0.19.0.sql \
|
||||
sql/pg_orrery--0.20.0.sql sql/pg_orrery--0.19.0--0.20.0.sql
|
||||
|
||||
# Our extension C sources
|
||||
OBJS = src/pg_orrery.o src/tle_type.o src/eci_type.o src/observer_type.o \
|
||||
@ -38,7 +39,8 @@ OBJS = src/pg_orrery.o src/tle_type.o src/eci_type.o src/observer_type.o \
|
||||
src/rise_set_funcs.o \
|
||||
src/constellation_data.o src/constellation_funcs.o \
|
||||
src/lunar_phase_funcs.o src/magnitude_funcs.o \
|
||||
src/eclipse_funcs.o src/libration_funcs.o
|
||||
src/eclipse_funcs.o src/libration_funcs.o \
|
||||
src/lagrange_funcs.o
|
||||
|
||||
# Vendored SGP4/SDP4 sources (pure C, from Bill Gray's sat_code, MIT license)
|
||||
SGP4_DIR = src/sgp4
|
||||
@ -62,7 +64,8 @@ REGRESS = tle_parse sgp4_propagate coord_transforms pass_prediction gist_index c
|
||||
v016_features \
|
||||
v017_features \
|
||||
v018_features \
|
||||
v019_features
|
||||
v019_features \
|
||||
v020_features
|
||||
REGRESS_OPTS = --inputdir=test
|
||||
|
||||
# Pure C — no C++ runtime needed. LAPACK for OD solver (dgelss_).
|
||||
|
||||
@ -1,4 +1,4 @@
|
||||
comment = 'A database orrery — celestial mechanics types and functions for PostgreSQL'
|
||||
default_version = '0.19.0'
|
||||
default_version = '0.20.0'
|
||||
module_pathname = '$libdir/pg_orrery'
|
||||
relocatable = true
|
||||
|
||||
244
sql/pg_orrery--0.19.0--0.20.0.sql
Normal file
244
sql/pg_orrery--0.19.0--0.20.0.sql
Normal file
@ -0,0 +1,244 @@
|
||||
-- pg_orrery 0.19.0 -> 0.20.0: Lagrange point support
|
||||
-- CR3BP equilibrium positions for Sun-planet, Earth-Moon, and planetary moon systems.
|
||||
|
||||
-- ============================================================
|
||||
-- Sun-planet Lagrange functions (5)
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION lagrange_heliocentric(int4, int4, timestamptz) RETURNS heliocentric
|
||||
AS 'MODULE_PATHNAME', 'lagrange_heliocentric'
|
||||
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lagrange_heliocentric(int4, int4, timestamptz) IS
|
||||
'Heliocentric ecliptic J2000 position of a Sun-planet Lagrange point. body_id: 1-8 (Mercury-Neptune), point_id: 1-5 (L1-L5).';
|
||||
|
||||
CREATE FUNCTION lagrange_observe(int4, int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME', 'lagrange_observe'
|
||||
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lagrange_observe(int4, int4, observer, timestamptz) IS
|
||||
'Observe a Sun-planet Lagrange point from a ground station. body_id: 1-8, point_id: 1-5.';
|
||||
|
||||
CREATE FUNCTION lagrange_equatorial(int4, int4, timestamptz) RETURNS equatorial
|
||||
AS 'MODULE_PATHNAME', 'lagrange_equatorial'
|
||||
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lagrange_equatorial(int4, int4, timestamptz) IS
|
||||
'Geocentric RA/Dec of a Sun-planet Lagrange point. body_id: 1-8, point_id: 1-5.';
|
||||
|
||||
CREATE FUNCTION lagrange_distance(int4, int4, heliocentric, timestamptz) RETURNS float8
|
||||
AS 'MODULE_PATHNAME', 'lagrange_distance'
|
||||
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lagrange_distance(int4, int4, heliocentric, timestamptz) IS
|
||||
'Distance (AU) from a heliocentric position to a Sun-planet Lagrange point.';
|
||||
|
||||
CREATE FUNCTION lagrange_distance_oe(int4, int4, orbital_elements, timestamptz) RETURNS float8
|
||||
AS 'MODULE_PATHNAME', 'lagrange_distance_oe'
|
||||
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lagrange_distance_oe(int4, int4, orbital_elements, timestamptz) IS
|
||||
'Distance (AU) from an asteroid/comet (orbital_elements) to a Sun-planet Lagrange point.';
|
||||
|
||||
-- ============================================================
|
||||
-- Earth-Moon Lagrange functions (2)
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION lunar_lagrange_observe(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME', 'lunar_lagrange_observe'
|
||||
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lunar_lagrange_observe(int4, observer, timestamptz) IS
|
||||
'Observe an Earth-Moon Lagrange point. point_id: 1-5 (L1-L5).';
|
||||
|
||||
CREATE FUNCTION lunar_lagrange_equatorial(int4, timestamptz) RETURNS equatorial
|
||||
AS 'MODULE_PATHNAME', 'lunar_lagrange_equatorial'
|
||||
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lunar_lagrange_equatorial(int4, timestamptz) IS
|
||||
'Geocentric RA/Dec of an Earth-Moon Lagrange point. point_id: 1-5 (L1-L5).';
|
||||
|
||||
-- ============================================================
|
||||
-- Planetary moon Lagrange functions (8)
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION galilean_lagrange_observe(int4, int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME', 'galilean_lagrange_observe'
|
||||
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION galilean_lagrange_observe(int4, int4, observer, timestamptz) IS
|
||||
'Observe a Jupiter-Galilean moon Lagrange point. body_id: 0-3 (Io-Callisto), point_id: 1-5.';
|
||||
|
||||
CREATE FUNCTION galilean_lagrange_equatorial(int4, int4, timestamptz) RETURNS equatorial
|
||||
AS 'MODULE_PATHNAME', 'galilean_lagrange_equatorial'
|
||||
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION galilean_lagrange_equatorial(int4, int4, timestamptz) IS
|
||||
'Geocentric RA/Dec of a Jupiter-Galilean moon Lagrange point. body_id: 0-3, point_id: 1-5.';
|
||||
|
||||
CREATE FUNCTION saturn_moon_lagrange_observe(int4, int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME', 'saturn_moon_lagrange_observe'
|
||||
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION saturn_moon_lagrange_observe(int4, int4, observer, timestamptz) IS
|
||||
'Observe a Saturn moon Lagrange point. body_id: 0-7 (Mimas-Hyperion), point_id: 1-5.';
|
||||
|
||||
CREATE FUNCTION saturn_moon_lagrange_equatorial(int4, int4, timestamptz) RETURNS equatorial
|
||||
AS 'MODULE_PATHNAME', 'saturn_moon_lagrange_equatorial'
|
||||
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION saturn_moon_lagrange_equatorial(int4, int4, timestamptz) IS
|
||||
'Geocentric RA/Dec of a Saturn moon Lagrange point. body_id: 0-7, point_id: 1-5.';
|
||||
|
||||
CREATE FUNCTION uranus_moon_lagrange_observe(int4, int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME', 'uranus_moon_lagrange_observe'
|
||||
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION uranus_moon_lagrange_observe(int4, int4, observer, timestamptz) IS
|
||||
'Observe a Uranus moon Lagrange point. body_id: 0-4 (Miranda-Oberon), point_id: 1-5.';
|
||||
|
||||
CREATE FUNCTION uranus_moon_lagrange_equatorial(int4, int4, timestamptz) RETURNS equatorial
|
||||
AS 'MODULE_PATHNAME', 'uranus_moon_lagrange_equatorial'
|
||||
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION uranus_moon_lagrange_equatorial(int4, int4, timestamptz) IS
|
||||
'Geocentric RA/Dec of a Uranus moon Lagrange point. body_id: 0-4, point_id: 1-5.';
|
||||
|
||||
CREATE FUNCTION mars_moon_lagrange_observe(int4, int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME', 'mars_moon_lagrange_observe'
|
||||
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION mars_moon_lagrange_observe(int4, int4, observer, timestamptz) IS
|
||||
'Observe a Mars moon Lagrange point. body_id: 0-1 (Phobos-Deimos), point_id: 1-5.';
|
||||
|
||||
CREATE FUNCTION mars_moon_lagrange_equatorial(int4, int4, timestamptz) RETURNS equatorial
|
||||
AS 'MODULE_PATHNAME', 'mars_moon_lagrange_equatorial'
|
||||
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION mars_moon_lagrange_equatorial(int4, int4, timestamptz) IS
|
||||
'Geocentric RA/Dec of a Mars moon Lagrange point. body_id: 0-1, point_id: 1-5.';
|
||||
|
||||
-- ============================================================
|
||||
-- Hill radius / zone / convenience (5)
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION hill_radius(int4, timestamptz) RETURNS float8
|
||||
AS 'MODULE_PATHNAME', 'hill_radius_func'
|
||||
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION hill_radius(int4, timestamptz) IS
|
||||
'Hill sphere radius (AU) for a Sun-planet system. body_id: 1-8.';
|
||||
|
||||
CREATE FUNCTION hill_radius_lunar(timestamptz) RETURNS float8
|
||||
AS 'MODULE_PATHNAME', 'hill_radius_lunar'
|
||||
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION hill_radius_lunar(timestamptz) IS
|
||||
'Hill sphere radius (AU) for the Earth-Moon system.';
|
||||
|
||||
CREATE FUNCTION lagrange_zone_radius(int4, int4, timestamptz) RETURNS float8
|
||||
AS 'MODULE_PATHNAME', 'lagrange_zone_radius_func'
|
||||
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lagrange_zone_radius(int4, int4, timestamptz) IS
|
||||
'Approximate libration zone radius (AU) for a Sun-planet Lagrange point.';
|
||||
|
||||
CREATE FUNCTION lagrange_mass_ratio(int4) RETURNS float8
|
||||
AS 'MODULE_PATHNAME', 'lagrange_mass_ratio'
|
||||
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lagrange_mass_ratio(int4) IS
|
||||
'CR3BP mass parameter mu = M_planet / (M_sun + M_planet) for debugging. body_id: 1-8.';
|
||||
|
||||
CREATE FUNCTION lagrange_point_name(int4) RETURNS text
|
||||
AS 'MODULE_PATHNAME', 'lagrange_point_name'
|
||||
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lagrange_point_name(int4) IS
|
||||
'Human-readable name for a Lagrange point ID (1->''L1'', ..., 5->''L5'').';
|
||||
|
||||
-- ============================================================
|
||||
-- DE variant functions (17) -- STABLE
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION lagrange_heliocentric_de(int4, int4, timestamptz) RETURNS heliocentric
|
||||
AS 'MODULE_PATHNAME', 'lagrange_heliocentric_de'
|
||||
LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lagrange_heliocentric_de(int4, int4, timestamptz) IS
|
||||
'DE variant of lagrange_heliocentric(). Falls back to VSOP87 if DE unavailable.';
|
||||
|
||||
CREATE FUNCTION lagrange_observe_de(int4, int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME', 'lagrange_observe_de'
|
||||
LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lagrange_observe_de(int4, int4, observer, timestamptz) IS
|
||||
'DE variant of lagrange_observe(). Falls back to VSOP87 if DE unavailable.';
|
||||
|
||||
CREATE FUNCTION lagrange_equatorial_de(int4, int4, timestamptz) RETURNS equatorial
|
||||
AS 'MODULE_PATHNAME', 'lagrange_equatorial_de'
|
||||
LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lagrange_equatorial_de(int4, int4, timestamptz) IS
|
||||
'DE variant of lagrange_equatorial(). Falls back to VSOP87 if DE unavailable.';
|
||||
|
||||
CREATE FUNCTION lagrange_distance_de(int4, int4, heliocentric, timestamptz) RETURNS float8
|
||||
AS 'MODULE_PATHNAME', 'lagrange_distance_de'
|
||||
LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lagrange_distance_de(int4, int4, heliocentric, timestamptz) IS
|
||||
'DE variant of lagrange_distance(). Falls back to VSOP87 if DE unavailable.';
|
||||
|
||||
CREATE FUNCTION lagrange_distance_oe_de(int4, int4, orbital_elements, timestamptz) RETURNS float8
|
||||
AS 'MODULE_PATHNAME', 'lagrange_distance_oe_de'
|
||||
LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lagrange_distance_oe_de(int4, int4, orbital_elements, timestamptz) IS
|
||||
'DE variant of lagrange_distance_oe(). Falls back to VSOP87 if DE unavailable.';
|
||||
|
||||
CREATE FUNCTION lunar_lagrange_observe_de(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME', 'lunar_lagrange_observe_de'
|
||||
LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lunar_lagrange_observe_de(int4, observer, timestamptz) IS
|
||||
'DE variant of lunar_lagrange_observe(). Falls back to ELP2000-82B if DE unavailable.';
|
||||
|
||||
CREATE FUNCTION lunar_lagrange_equatorial_de(int4, timestamptz) RETURNS equatorial
|
||||
AS 'MODULE_PATHNAME', 'lunar_lagrange_equatorial_de'
|
||||
LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lunar_lagrange_equatorial_de(int4, timestamptz) IS
|
||||
'DE variant of lunar_lagrange_equatorial(). Falls back to ELP2000-82B if DE unavailable.';
|
||||
|
||||
CREATE FUNCTION galilean_lagrange_observe_de(int4, int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME', 'galilean_lagrange_observe_de'
|
||||
LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION galilean_lagrange_observe_de(int4, int4, observer, timestamptz) IS
|
||||
'DE variant of galilean_lagrange_observe(). Falls back to VSOP87 if DE unavailable.';
|
||||
|
||||
CREATE FUNCTION galilean_lagrange_equatorial_de(int4, int4, timestamptz) RETURNS equatorial
|
||||
AS 'MODULE_PATHNAME', 'galilean_lagrange_equatorial_de'
|
||||
LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION galilean_lagrange_equatorial_de(int4, int4, timestamptz) IS
|
||||
'DE variant of galilean_lagrange_equatorial(). Falls back to VSOP87 if DE unavailable.';
|
||||
|
||||
CREATE FUNCTION saturn_moon_lagrange_observe_de(int4, int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME', 'saturn_moon_lagrange_observe_de'
|
||||
LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION saturn_moon_lagrange_observe_de(int4, int4, observer, timestamptz) IS
|
||||
'DE variant of saturn_moon_lagrange_observe(). Falls back to VSOP87 if DE unavailable.';
|
||||
|
||||
CREATE FUNCTION saturn_moon_lagrange_equatorial_de(int4, int4, timestamptz) RETURNS equatorial
|
||||
AS 'MODULE_PATHNAME', 'saturn_moon_lagrange_equatorial_de'
|
||||
LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION saturn_moon_lagrange_equatorial_de(int4, int4, timestamptz) IS
|
||||
'DE variant of saturn_moon_lagrange_equatorial(). Falls back to VSOP87 if DE unavailable.';
|
||||
|
||||
CREATE FUNCTION uranus_moon_lagrange_observe_de(int4, int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME', 'uranus_moon_lagrange_observe_de'
|
||||
LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION uranus_moon_lagrange_observe_de(int4, int4, observer, timestamptz) IS
|
||||
'DE variant of uranus_moon_lagrange_observe(). Falls back to VSOP87 if DE unavailable.';
|
||||
|
||||
CREATE FUNCTION uranus_moon_lagrange_equatorial_de(int4, int4, timestamptz) RETURNS equatorial
|
||||
AS 'MODULE_PATHNAME', 'uranus_moon_lagrange_equatorial_de'
|
||||
LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION uranus_moon_lagrange_equatorial_de(int4, int4, timestamptz) IS
|
||||
'DE variant of uranus_moon_lagrange_equatorial(). Falls back to VSOP87 if DE unavailable.';
|
||||
|
||||
CREATE FUNCTION mars_moon_lagrange_observe_de(int4, int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME', 'mars_moon_lagrange_observe_de'
|
||||
LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION mars_moon_lagrange_observe_de(int4, int4, observer, timestamptz) IS
|
||||
'DE variant of mars_moon_lagrange_observe(). Falls back to VSOP87 if DE unavailable.';
|
||||
|
||||
CREATE FUNCTION mars_moon_lagrange_equatorial_de(int4, int4, timestamptz) RETURNS equatorial
|
||||
AS 'MODULE_PATHNAME', 'mars_moon_lagrange_equatorial_de'
|
||||
LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION mars_moon_lagrange_equatorial_de(int4, int4, timestamptz) IS
|
||||
'DE variant of mars_moon_lagrange_equatorial(). Falls back to VSOP87 if DE unavailable.';
|
||||
|
||||
CREATE FUNCTION hill_radius_de(int4, timestamptz) RETURNS float8
|
||||
AS 'MODULE_PATHNAME', 'hill_radius_de'
|
||||
LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION hill_radius_de(int4, timestamptz) IS
|
||||
'DE variant of hill_radius(). Falls back to VSOP87 if DE unavailable.';
|
||||
|
||||
CREATE FUNCTION lagrange_zone_radius_de(int4, int4, timestamptz) RETURNS float8
|
||||
AS 'MODULE_PATHNAME', 'lagrange_zone_radius_de'
|
||||
LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lagrange_zone_radius_de(int4, int4, timestamptz) IS
|
||||
'DE variant of lagrange_zone_radius(). Falls back to VSOP87 if DE unavailable.';
|
||||
2202
sql/pg_orrery--0.20.0.sql
Normal file
2202
sql/pg_orrery--0.20.0.sql
Normal file
File diff suppressed because it is too large
Load Diff
968
src/de_funcs.c
968
src/de_funcs.c
@ -33,6 +33,7 @@
|
||||
#include "tass17.h"
|
||||
#include "gust86.h"
|
||||
#include "marssat.h"
|
||||
#include "lagrange.h"
|
||||
|
||||
#include <math.h>
|
||||
#include <string.h>
|
||||
@ -62,6 +63,25 @@ PG_FUNCTION_INFO_V1(uranus_moon_equatorial_de);
|
||||
PG_FUNCTION_INFO_V1(mars_moon_equatorial_de);
|
||||
PG_FUNCTION_INFO_V1(pg_orrery_ephemeris_info);
|
||||
|
||||
/* Lagrange DE variants */
|
||||
PG_FUNCTION_INFO_V1(lagrange_heliocentric_de);
|
||||
PG_FUNCTION_INFO_V1(lagrange_observe_de);
|
||||
PG_FUNCTION_INFO_V1(lagrange_equatorial_de);
|
||||
PG_FUNCTION_INFO_V1(lagrange_distance_de);
|
||||
PG_FUNCTION_INFO_V1(lagrange_distance_oe_de);
|
||||
PG_FUNCTION_INFO_V1(lunar_lagrange_observe_de);
|
||||
PG_FUNCTION_INFO_V1(lunar_lagrange_equatorial_de);
|
||||
PG_FUNCTION_INFO_V1(galilean_lagrange_observe_de);
|
||||
PG_FUNCTION_INFO_V1(galilean_lagrange_equatorial_de);
|
||||
PG_FUNCTION_INFO_V1(saturn_moon_lagrange_observe_de);
|
||||
PG_FUNCTION_INFO_V1(saturn_moon_lagrange_equatorial_de);
|
||||
PG_FUNCTION_INFO_V1(uranus_moon_lagrange_observe_de);
|
||||
PG_FUNCTION_INFO_V1(uranus_moon_lagrange_equatorial_de);
|
||||
PG_FUNCTION_INFO_V1(mars_moon_lagrange_observe_de);
|
||||
PG_FUNCTION_INFO_V1(mars_moon_lagrange_equatorial_de);
|
||||
PG_FUNCTION_INFO_V1(hill_radius_de);
|
||||
PG_FUNCTION_INFO_V1(lagrange_zone_radius_de);
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* planet_heliocentric_de(body_id int, timestamptz) -> heliocentric
|
||||
@ -1306,3 +1326,951 @@ pg_orrery_ephemeris_info(PG_FUNCTION_ARGS)
|
||||
tuple = heap_form_tuple(tupdesc, values, nulls);
|
||||
PG_RETURN_DATUM(HeapTupleGetDatum(tuple));
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* Lagrange point functions with DE ephemeris
|
||||
*
|
||||
* DE variants of the Lagrange functions in lagrange_funcs.c.
|
||||
* Each uses DE for planet/Earth positions where possible,
|
||||
* falling back to VSOP87/ELP2000-82B on any DE failure.
|
||||
* Rule 7 always holds: both target and Earth from the same provider.
|
||||
* ================================================================
|
||||
*/
|
||||
|
||||
|
||||
/*
|
||||
* Validate body_id and point_id for Sun-planet Lagrange functions.
|
||||
* Duplicated from lagrange_funcs.c (no cross-TU symbols).
|
||||
*/
|
||||
static void
|
||||
validate_lagrange_args_de(int body_id, int point_id, const char *func_name)
|
||||
{
|
||||
if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("%s: body_id %d must be 1-8 (Mercury-Neptune)",
|
||||
func_name, body_id)));
|
||||
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("%s: point_id %d must be 1-5 (L1-L5)",
|
||||
func_name, point_id)));
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Compute Sun-planet Lagrange point using DE for planet position+velocity.
|
||||
* Falls back to VSOP87 if DE unavailable (rule 7: consistent provider).
|
||||
*
|
||||
* Returns 0 on success, -1 on solver failure.
|
||||
* Sets *used_de to indicate which provider was used.
|
||||
*/
|
||||
static int
|
||||
compute_sun_planet_lagrange_de(int body_id, int point_id, double jd,
|
||||
double result[3], bool *used_de)
|
||||
{
|
||||
double planet_xyz[6];
|
||||
double sun[3] = {0.0, 0.0, 0.0};
|
||||
double planet_vel[3];
|
||||
double ratio, mu;
|
||||
|
||||
/* Try DE for planet position */
|
||||
if (eph_de_planet(body_id, jd, planet_xyz))
|
||||
{
|
||||
/* Velocity via central difference (DE provides position only) */
|
||||
double pos_fwd[6], pos_bwd[6];
|
||||
double dt = 0.01; /* days */
|
||||
|
||||
bool got_fwd = eph_de_planet(body_id, jd + dt, pos_fwd);
|
||||
bool got_bwd = eph_de_planet(body_id, jd - dt, pos_bwd);
|
||||
|
||||
if (got_fwd && got_bwd)
|
||||
{
|
||||
planet_vel[0] = (pos_fwd[0] - pos_bwd[0]) / (2.0 * dt);
|
||||
planet_vel[1] = (pos_fwd[1] - pos_bwd[1]) / (2.0 * dt);
|
||||
planet_vel[2] = (pos_fwd[2] - pos_bwd[2]) / (2.0 * dt);
|
||||
*used_de = true;
|
||||
}
|
||||
else
|
||||
{
|
||||
/* DE boundary straddled — fall through to VSOP87 */
|
||||
goto vsop87_fallback;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
vsop87_fallback:
|
||||
{
|
||||
int vsop_body = body_id - 1;
|
||||
|
||||
if (eph_de_available())
|
||||
ereport(NOTICE,
|
||||
(errmsg("DE ephemeris unavailable for this query, falling back to VSOP87")));
|
||||
|
||||
GetVsop87Coor(jd, vsop_body, planet_xyz);
|
||||
/* VSOP87 provides analytic velocity in xyz[3..5] */
|
||||
planet_vel[0] = planet_xyz[3];
|
||||
planet_vel[1] = planet_xyz[4];
|
||||
planet_vel[2] = planet_xyz[5];
|
||||
*used_de = false;
|
||||
}
|
||||
}
|
||||
|
||||
ratio = sun_planet_ratio(body_id);
|
||||
mu = mu_from_ratio(ratio);
|
||||
|
||||
return lagrange_position(sun, planet_xyz, planet_vel, mu, point_id, result);
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Compute Earth-Moon Lagrange point using DE for Moon position.
|
||||
* Falls back to ELP2000-82B if DE unavailable.
|
||||
*
|
||||
* Moon velocity always via central difference (neither DE nor ELP
|
||||
* provide direct velocity). Result is geocentric ecliptic J2000.
|
||||
*/
|
||||
static int
|
||||
compute_lunar_lagrange_de(int point_id, double jd, double result_geo[3],
|
||||
bool *used_de)
|
||||
{
|
||||
double moon_pos[3], moon_fwd[3], moon_bwd[3];
|
||||
double moon_vel[3];
|
||||
double earth[3] = {0.0, 0.0, 0.0};
|
||||
double mu;
|
||||
double dt = 0.001; /* days */
|
||||
|
||||
if (eph_de_moon(jd, moon_pos))
|
||||
{
|
||||
bool got_fwd = eph_de_moon(jd + dt, moon_fwd);
|
||||
bool got_bwd = eph_de_moon(jd - dt, moon_bwd);
|
||||
|
||||
if (got_fwd && got_bwd)
|
||||
{
|
||||
*used_de = true;
|
||||
}
|
||||
else
|
||||
{
|
||||
/* DE boundary straddled */
|
||||
goto elp_fallback;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
elp_fallback:
|
||||
if (eph_de_available())
|
||||
ereport(NOTICE,
|
||||
(errmsg("DE ephemeris unavailable, falling back to ELP2000-82B")));
|
||||
|
||||
GetElp82bCoor(jd, moon_pos);
|
||||
GetElp82bCoor(jd + dt, moon_fwd);
|
||||
GetElp82bCoor(jd - dt, moon_bwd);
|
||||
*used_de = false;
|
||||
}
|
||||
|
||||
moon_vel[0] = (moon_fwd[0] - moon_bwd[0]) / (2.0 * dt);
|
||||
moon_vel[1] = (moon_fwd[1] - moon_bwd[1]) / (2.0 * dt);
|
||||
moon_vel[2] = (moon_fwd[2] - moon_bwd[2]) / (2.0 * dt);
|
||||
|
||||
mu = mu_from_ratio(EARTH_MOON_EMRAT);
|
||||
|
||||
return lagrange_position(earth, moon_pos, moon_vel, mu, point_id,
|
||||
result_geo);
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Observe a planetary moon Lagrange point using DE for parent planet
|
||||
* and Earth positions. Falls back to VSOP87 if DE unavailable.
|
||||
*
|
||||
* lp_rel[3]: L-point offset relative to parent (ecliptic J2000, AU)
|
||||
* parent_body_id: pg_orrery body ID (5=Jupiter, 6=Saturn, etc.)
|
||||
*/
|
||||
static void
|
||||
observe_pmoon_lagrange_de(const double lp_rel[3], int parent_body_id,
|
||||
double jd, const pg_observer *obs,
|
||||
pg_topocentric *result)
|
||||
{
|
||||
double parent_xyz[6];
|
||||
double earth_xyz[6];
|
||||
double geo_ecl[3];
|
||||
bool have_de;
|
||||
|
||||
/* Rule 7: both parent and Earth from same provider */
|
||||
have_de = eph_de_planet(parent_body_id, jd, parent_xyz) &&
|
||||
eph_de_earth(jd, earth_xyz);
|
||||
|
||||
if (!have_de)
|
||||
{
|
||||
int vsop_parent = parent_body_id - 1;
|
||||
|
||||
if (eph_de_available())
|
||||
ereport(NOTICE,
|
||||
(errmsg("DE ephemeris unavailable, falling back to VSOP87")));
|
||||
|
||||
GetVsop87Coor(jd, vsop_parent, parent_xyz);
|
||||
GetVsop87Coor(jd, 2, earth_xyz);
|
||||
}
|
||||
|
||||
geo_ecl[0] = (parent_xyz[0] + lp_rel[0]) - earth_xyz[0];
|
||||
geo_ecl[1] = (parent_xyz[1] + lp_rel[1]) - earth_xyz[1];
|
||||
geo_ecl[2] = (parent_xyz[2] + lp_rel[2]) - earth_xyz[2];
|
||||
|
||||
observe_from_geocentric(geo_ecl, jd, obs, result);
|
||||
}
|
||||
|
||||
|
||||
static void
|
||||
equatorial_pmoon_lagrange_de(const double lp_rel[3], int parent_body_id,
|
||||
double jd, pg_equatorial *result)
|
||||
{
|
||||
double parent_xyz[6];
|
||||
double earth_xyz[6];
|
||||
double geo_ecl[3];
|
||||
bool have_de;
|
||||
|
||||
have_de = eph_de_planet(parent_body_id, jd, parent_xyz) &&
|
||||
eph_de_earth(jd, earth_xyz);
|
||||
|
||||
if (!have_de)
|
||||
{
|
||||
int vsop_parent = parent_body_id - 1;
|
||||
|
||||
if (eph_de_available())
|
||||
ereport(NOTICE,
|
||||
(errmsg("DE ephemeris unavailable, falling back to VSOP87")));
|
||||
|
||||
GetVsop87Coor(jd, vsop_parent, parent_xyz);
|
||||
GetVsop87Coor(jd, 2, earth_xyz);
|
||||
}
|
||||
|
||||
geo_ecl[0] = (parent_xyz[0] + lp_rel[0]) - earth_xyz[0];
|
||||
geo_ecl[1] = (parent_xyz[1] + lp_rel[1]) - earth_xyz[1];
|
||||
geo_ecl[2] = (parent_xyz[2] + lp_rel[2]) - earth_xyz[2];
|
||||
|
||||
geocentric_to_equatorial(geo_ecl, jd, result);
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Compute planetary moon Lagrange offset (no cross-TU call).
|
||||
* Duplicated from lagrange_funcs.c.
|
||||
*/
|
||||
static int
|
||||
compute_planetary_moon_lagrange_de(const double moon_rel[3],
|
||||
const double moon_vel[3],
|
||||
char family, int moon_id,
|
||||
int point_id,
|
||||
double lp_rel[3])
|
||||
{
|
||||
double planet_origin[3] = {0.0, 0.0, 0.0};
|
||||
double ratio, mu;
|
||||
|
||||
ratio = planet_moon_ratio(family, moon_id);
|
||||
if (ratio < 0.0)
|
||||
return -1;
|
||||
|
||||
mu = mu_from_ratio(ratio);
|
||||
|
||||
return lagrange_position(planet_origin, moon_rel, moon_vel, mu,
|
||||
point_id, lp_rel);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* 1. lagrange_heliocentric_de(body_id, point_id, timestamptz)
|
||||
* -> heliocentric
|
||||
*
|
||||
* DE variant of lagrange_heliocentric(). STABLE.
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
lagrange_heliocentric_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
int64 ts = PG_GETARG_INT64(2);
|
||||
double jd;
|
||||
double lp[3];
|
||||
bool used_de;
|
||||
pg_heliocentric *result;
|
||||
|
||||
validate_lagrange_args_de(body_id, point_id, "lagrange_heliocentric_de");
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
|
||||
if (compute_sun_planet_lagrange_de(body_id, point_id, jd, lp, &used_de) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("lagrange_heliocentric_de: solver failed for body %d L%d",
|
||||
body_id, point_id)));
|
||||
|
||||
result = (pg_heliocentric *) palloc(sizeof(pg_heliocentric));
|
||||
result->x = lp[0];
|
||||
result->y = lp[1];
|
||||
result->z = lp[2];
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* 2. lagrange_observe_de(body_id, point_id, observer, timestamptz)
|
||||
* -> topocentric
|
||||
*
|
||||
* DE variant of lagrange_observe(). STABLE.
|
||||
* Rule 7: L-point + Earth from same provider.
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
lagrange_observe_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2);
|
||||
int64 ts = PG_GETARG_INT64(3);
|
||||
double jd;
|
||||
double lp[3], earth_xyz[6], geo_ecl[3];
|
||||
bool used_de;
|
||||
pg_topocentric *result;
|
||||
|
||||
validate_lagrange_args_de(body_id, point_id, "lagrange_observe_de");
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
|
||||
if (compute_sun_planet_lagrange_de(body_id, point_id, jd, lp, &used_de) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("lagrange_observe_de: solver failed")));
|
||||
|
||||
/* Earth from same provider as the L-point (rule 7) */
|
||||
if (used_de && eph_de_earth(jd, earth_xyz))
|
||||
{
|
||||
/* DE succeeded for both */
|
||||
}
|
||||
else
|
||||
{
|
||||
GetVsop87Coor(jd, 2, earth_xyz);
|
||||
}
|
||||
|
||||
geo_ecl[0] = lp[0] - earth_xyz[0];
|
||||
geo_ecl[1] = lp[1] - earth_xyz[1];
|
||||
geo_ecl[2] = lp[2] - earth_xyz[2];
|
||||
|
||||
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
|
||||
observe_from_geocentric(geo_ecl, jd, obs, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* 3. lagrange_equatorial_de(body_id, point_id, timestamptz)
|
||||
* -> equatorial
|
||||
*
|
||||
* DE variant of lagrange_equatorial(). STABLE.
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
lagrange_equatorial_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
int64 ts = PG_GETARG_INT64(2);
|
||||
double jd;
|
||||
double lp[3], earth_xyz[6], geo_ecl[3];
|
||||
bool used_de;
|
||||
pg_equatorial *result;
|
||||
|
||||
validate_lagrange_args_de(body_id, point_id, "lagrange_equatorial_de");
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
|
||||
if (compute_sun_planet_lagrange_de(body_id, point_id, jd, lp, &used_de) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("lagrange_equatorial_de: solver failed")));
|
||||
|
||||
if (used_de && eph_de_earth(jd, earth_xyz))
|
||||
{
|
||||
/* DE succeeded */
|
||||
}
|
||||
else
|
||||
{
|
||||
GetVsop87Coor(jd, 2, earth_xyz);
|
||||
}
|
||||
|
||||
geo_ecl[0] = lp[0] - earth_xyz[0];
|
||||
geo_ecl[1] = lp[1] - earth_xyz[1];
|
||||
geo_ecl[2] = lp[2] - earth_xyz[2];
|
||||
|
||||
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
|
||||
geocentric_to_equatorial(geo_ecl, jd, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* 4. lagrange_distance_de(body_id, point_id, heliocentric, timestamptz)
|
||||
* -> float8
|
||||
*
|
||||
* Distance from a heliocentric position to a Sun-planet L-point (AU).
|
||||
* Uses DE for planet position. STABLE.
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
lagrange_distance_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
pg_heliocentric *pos = (pg_heliocentric *) PG_GETARG_POINTER(2);
|
||||
int64 ts = PG_GETARG_INT64(3);
|
||||
double jd;
|
||||
double lp[3];
|
||||
double dx, dy, dz, dist;
|
||||
bool used_de;
|
||||
|
||||
validate_lagrange_args_de(body_id, point_id, "lagrange_distance_de");
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
|
||||
if (compute_sun_planet_lagrange_de(body_id, point_id, jd, lp, &used_de) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("lagrange_distance_de: solver failed")));
|
||||
|
||||
dx = pos->x - lp[0];
|
||||
dy = pos->y - lp[1];
|
||||
dz = pos->z - lp[2];
|
||||
dist = sqrt(dx*dx + dy*dy + dz*dz);
|
||||
|
||||
PG_RETURN_FLOAT8(dist);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* 5. lagrange_distance_oe_de(body_id, point_id, orbital_elements, timestamptz)
|
||||
* -> float8
|
||||
*
|
||||
* Distance from an asteroid/comet to a Sun-planet L-point (AU).
|
||||
* Uses DE for L-point computation. STABLE.
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
lagrange_distance_oe_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
pg_orbital_elements *oe = (pg_orbital_elements *) PG_GETARG_POINTER(2);
|
||||
int64 ts = PG_GETARG_INT64(3);
|
||||
double jd;
|
||||
double lp[3], body_pos[3];
|
||||
double dx, dy, dz, dist;
|
||||
bool used_de;
|
||||
|
||||
validate_lagrange_args_de(body_id, point_id, "lagrange_distance_oe_de");
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
|
||||
if (compute_sun_planet_lagrange_de(body_id, point_id, jd, lp, &used_de) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("lagrange_distance_oe_de: solver failed")));
|
||||
|
||||
kepler_position(oe->q, oe->e, oe->inc, oe->arg_peri, oe->raan,
|
||||
oe->tp, jd, body_pos);
|
||||
|
||||
dx = body_pos[0] - lp[0];
|
||||
dy = body_pos[1] - lp[1];
|
||||
dz = body_pos[2] - lp[2];
|
||||
dist = sqrt(dx*dx + dy*dy + dz*dz);
|
||||
|
||||
PG_RETURN_FLOAT8(dist);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* 6. lunar_lagrange_observe_de(point_id, observer, timestamptz)
|
||||
* -> topocentric
|
||||
*
|
||||
* Earth-Moon L-point using DE Moon position. STABLE.
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
lunar_lagrange_observe_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 point_id = PG_GETARG_INT32(0);
|
||||
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(1);
|
||||
int64 ts = PG_GETARG_INT64(2);
|
||||
double jd;
|
||||
double lp_geo[3];
|
||||
bool used_de;
|
||||
pg_topocentric *result;
|
||||
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("lunar_lagrange_observe_de: point_id %d must be 1-5",
|
||||
point_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
|
||||
if (compute_lunar_lagrange_de(point_id, jd, lp_geo, &used_de) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("lunar_lagrange_observe_de: solver failed")));
|
||||
|
||||
/* lp_geo is already geocentric ecliptic J2000 (AU) */
|
||||
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
|
||||
observe_from_geocentric(lp_geo, jd, obs, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* 7. lunar_lagrange_equatorial_de(point_id, timestamptz) -> equatorial
|
||||
*
|
||||
* Earth-Moon L-point using DE Moon position. STABLE.
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
lunar_lagrange_equatorial_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 point_id = PG_GETARG_INT32(0);
|
||||
int64 ts = PG_GETARG_INT64(1);
|
||||
double jd;
|
||||
double lp_geo[3];
|
||||
bool used_de;
|
||||
pg_equatorial *result;
|
||||
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("lunar_lagrange_equatorial_de: point_id %d must be 1-5",
|
||||
point_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
|
||||
if (compute_lunar_lagrange_de(point_id, jd, lp_geo, &used_de) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("lunar_lagrange_equatorial_de: solver failed")));
|
||||
|
||||
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
|
||||
geocentric_to_equatorial(lp_geo, jd, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* Planetary moon Lagrange points with DE parent positions
|
||||
*
|
||||
* Moon-theory offset (L1.2, TASS17, GUST86, MarsSat) computed
|
||||
* relative to parent, then Lagrange solver applied. Parent planet
|
||||
* and Earth positions from DE (with VSOP87 fallback).
|
||||
* ================================================================
|
||||
*/
|
||||
|
||||
|
||||
/* ── 8. Galilean Lagrange observe DE ─────────────────────── */
|
||||
|
||||
Datum
|
||||
galilean_lagrange_observe_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2);
|
||||
int64 ts = PG_GETARG_INT64(3);
|
||||
double jd;
|
||||
double moon_xyz[3], moon_vel[3], lp_rel[3];
|
||||
pg_topocentric *result;
|
||||
|
||||
if (body_id < L12_IO || body_id > L12_CALLISTO)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("galilean_lagrange_observe_de: body_id %d must be 0-3",
|
||||
body_id)));
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("galilean_lagrange_observe_de: point_id %d must be 1-5",
|
||||
point_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetL12Coor(jd, body_id, moon_xyz, moon_vel);
|
||||
|
||||
if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 'g', body_id,
|
||||
point_id, lp_rel) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("galilean_lagrange_observe_de: solver failed")));
|
||||
|
||||
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
|
||||
observe_pmoon_lagrange_de(lp_rel, BODY_JUPITER, jd, obs, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ── 9. Galilean Lagrange equatorial DE ──────────────────── */
|
||||
|
||||
Datum
|
||||
galilean_lagrange_equatorial_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
int64 ts = PG_GETARG_INT64(2);
|
||||
double jd;
|
||||
double moon_xyz[3], moon_vel[3], lp_rel[3];
|
||||
pg_equatorial *result;
|
||||
|
||||
if (body_id < L12_IO || body_id > L12_CALLISTO)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("galilean_lagrange_equatorial_de: body_id %d must be 0-3",
|
||||
body_id)));
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("galilean_lagrange_equatorial_de: point_id %d must be 1-5",
|
||||
point_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetL12Coor(jd, body_id, moon_xyz, moon_vel);
|
||||
|
||||
if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 'g', body_id,
|
||||
point_id, lp_rel) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("galilean_lagrange_equatorial_de: solver failed")));
|
||||
|
||||
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
|
||||
equatorial_pmoon_lagrange_de(lp_rel, BODY_JUPITER, jd, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ── 10. Saturn moon Lagrange observe DE ─────────────────── */
|
||||
|
||||
Datum
|
||||
saturn_moon_lagrange_observe_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2);
|
||||
int64 ts = PG_GETARG_INT64(3);
|
||||
double jd;
|
||||
double moon_xyz[3], moon_vel[3], lp_rel[3];
|
||||
pg_topocentric *result;
|
||||
|
||||
if (body_id < TASS17_MIMAS || body_id > TASS17_HYPERION)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("saturn_moon_lagrange_observe_de: body_id %d must be 0-7",
|
||||
body_id)));
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("saturn_moon_lagrange_observe_de: point_id %d must be 1-5",
|
||||
point_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetTass17Coor(jd, body_id, moon_xyz, moon_vel);
|
||||
|
||||
if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 's', body_id,
|
||||
point_id, lp_rel) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("saturn_moon_lagrange_observe_de: solver failed")));
|
||||
|
||||
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
|
||||
observe_pmoon_lagrange_de(lp_rel, BODY_SATURN, jd, obs, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ── 11. Saturn moon Lagrange equatorial DE ──────────────── */
|
||||
|
||||
Datum
|
||||
saturn_moon_lagrange_equatorial_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
int64 ts = PG_GETARG_INT64(2);
|
||||
double jd;
|
||||
double moon_xyz[3], moon_vel[3], lp_rel[3];
|
||||
pg_equatorial *result;
|
||||
|
||||
if (body_id < TASS17_MIMAS || body_id > TASS17_HYPERION)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("saturn_moon_lagrange_equatorial_de: body_id %d must be 0-7",
|
||||
body_id)));
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("saturn_moon_lagrange_equatorial_de: point_id %d must be 1-5",
|
||||
point_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetTass17Coor(jd, body_id, moon_xyz, moon_vel);
|
||||
|
||||
if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 's', body_id,
|
||||
point_id, lp_rel) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("saturn_moon_lagrange_equatorial_de: solver failed")));
|
||||
|
||||
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
|
||||
equatorial_pmoon_lagrange_de(lp_rel, BODY_SATURN, jd, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ── 12. Uranus moon Lagrange observe DE ─────────────────── */
|
||||
|
||||
Datum
|
||||
uranus_moon_lagrange_observe_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2);
|
||||
int64 ts = PG_GETARG_INT64(3);
|
||||
double jd;
|
||||
double moon_xyz[3], moon_vel[3], lp_rel[3];
|
||||
pg_topocentric *result;
|
||||
|
||||
if (body_id < GUST86_MIRANDA || body_id > GUST86_OBERON)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("uranus_moon_lagrange_observe_de: body_id %d must be 0-4",
|
||||
body_id)));
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("uranus_moon_lagrange_observe_de: point_id %d must be 1-5",
|
||||
point_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetGust86Coor(jd, body_id, moon_xyz, moon_vel);
|
||||
|
||||
if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 'u', body_id,
|
||||
point_id, lp_rel) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("uranus_moon_lagrange_observe_de: solver failed")));
|
||||
|
||||
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
|
||||
observe_pmoon_lagrange_de(lp_rel, BODY_URANUS, jd, obs, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ── 13. Uranus moon Lagrange equatorial DE ──────────────── */
|
||||
|
||||
Datum
|
||||
uranus_moon_lagrange_equatorial_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
int64 ts = PG_GETARG_INT64(2);
|
||||
double jd;
|
||||
double moon_xyz[3], moon_vel[3], lp_rel[3];
|
||||
pg_equatorial *result;
|
||||
|
||||
if (body_id < GUST86_MIRANDA || body_id > GUST86_OBERON)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("uranus_moon_lagrange_equatorial_de: body_id %d must be 0-4",
|
||||
body_id)));
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("uranus_moon_lagrange_equatorial_de: point_id %d must be 1-5",
|
||||
point_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetGust86Coor(jd, body_id, moon_xyz, moon_vel);
|
||||
|
||||
if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 'u', body_id,
|
||||
point_id, lp_rel) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("uranus_moon_lagrange_equatorial_de: solver failed")));
|
||||
|
||||
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
|
||||
equatorial_pmoon_lagrange_de(lp_rel, BODY_URANUS, jd, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ── 14. Mars moon Lagrange observe DE ───────────────────── */
|
||||
|
||||
Datum
|
||||
mars_moon_lagrange_observe_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2);
|
||||
int64 ts = PG_GETARG_INT64(3);
|
||||
double jd;
|
||||
double moon_xyz[3], moon_vel[3], lp_rel[3];
|
||||
pg_topocentric *result;
|
||||
|
||||
if (body_id < MARS_SAT_PHOBOS || body_id > MARS_SAT_DEIMOS)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("mars_moon_lagrange_observe_de: body_id %d must be 0-1",
|
||||
body_id)));
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("mars_moon_lagrange_observe_de: point_id %d must be 1-5",
|
||||
point_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetMarsSatCoor(jd, body_id, moon_xyz, moon_vel);
|
||||
|
||||
if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 'm', body_id,
|
||||
point_id, lp_rel) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("mars_moon_lagrange_observe_de: solver failed")));
|
||||
|
||||
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
|
||||
observe_pmoon_lagrange_de(lp_rel, BODY_MARS, jd, obs, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ── 15. Mars moon Lagrange equatorial DE ────────────────── */
|
||||
|
||||
Datum
|
||||
mars_moon_lagrange_equatorial_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
int64 ts = PG_GETARG_INT64(2);
|
||||
double jd;
|
||||
double moon_xyz[3], moon_vel[3], lp_rel[3];
|
||||
pg_equatorial *result;
|
||||
|
||||
if (body_id < MARS_SAT_PHOBOS || body_id > MARS_SAT_DEIMOS)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("mars_moon_lagrange_equatorial_de: body_id %d must be 0-1",
|
||||
body_id)));
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("mars_moon_lagrange_equatorial_de: point_id %d must be 1-5",
|
||||
point_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetMarsSatCoor(jd, body_id, moon_xyz, moon_vel);
|
||||
|
||||
if (compute_planetary_moon_lagrange_de(moon_xyz, moon_vel, 'm', body_id,
|
||||
point_id, lp_rel) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("mars_moon_lagrange_equatorial_de: solver failed")));
|
||||
|
||||
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
|
||||
equatorial_pmoon_lagrange_de(lp_rel, BODY_MARS, jd, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* 16. hill_radius_de(body_id, timestamptz) -> float8 (AU)
|
||||
*
|
||||
* Hill radius using DE for planet heliocentric distance. STABLE.
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
hill_radius_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int64 ts = PG_GETARG_INT64(1);
|
||||
double jd;
|
||||
double planet_xyz[6], sep, ratio, mu;
|
||||
|
||||
if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("hill_radius_de: body_id %d must be 1-8", body_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
|
||||
if (!eph_de_planet(body_id, jd, planet_xyz))
|
||||
{
|
||||
if (eph_de_available())
|
||||
ereport(NOTICE,
|
||||
(errmsg("DE ephemeris unavailable, falling back to VSOP87")));
|
||||
|
||||
GetVsop87Coor(jd, body_id - 1, planet_xyz);
|
||||
}
|
||||
|
||||
sep = sqrt(planet_xyz[0]*planet_xyz[0] +
|
||||
planet_xyz[1]*planet_xyz[1] +
|
||||
planet_xyz[2]*planet_xyz[2]);
|
||||
|
||||
ratio = sun_planet_ratio(body_id);
|
||||
mu = mu_from_ratio(ratio);
|
||||
|
||||
PG_RETURN_FLOAT8(lagrange_hill_radius(sep, mu));
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* 17. lagrange_zone_radius_de(body_id, point_id, timestamptz)
|
||||
* -> float8 (AU)
|
||||
*
|
||||
* Libration zone radius using DE for planet distance. STABLE.
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
lagrange_zone_radius_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
int64 ts = PG_GETARG_INT64(2);
|
||||
double jd;
|
||||
double planet_xyz[6], sep, ratio, mu, zr;
|
||||
|
||||
validate_lagrange_args_de(body_id, point_id, "lagrange_zone_radius_de");
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
|
||||
if (!eph_de_planet(body_id, jd, planet_xyz))
|
||||
{
|
||||
if (eph_de_available())
|
||||
ereport(NOTICE,
|
||||
(errmsg("DE ephemeris unavailable, falling back to VSOP87")));
|
||||
|
||||
GetVsop87Coor(jd, body_id - 1, planet_xyz);
|
||||
}
|
||||
|
||||
sep = sqrt(planet_xyz[0]*planet_xyz[0] +
|
||||
planet_xyz[1]*planet_xyz[1] +
|
||||
planet_xyz[2]*planet_xyz[2]);
|
||||
|
||||
ratio = sun_planet_ratio(body_id);
|
||||
mu = mu_from_ratio(ratio);
|
||||
|
||||
zr = lagrange_zone_radius(sep, mu, point_id);
|
||||
|
||||
PG_RETURN_FLOAT8(zr);
|
||||
}
|
||||
|
||||
916
src/lagrange_funcs.c
Normal file
916
src/lagrange_funcs.c
Normal file
@ -0,0 +1,916 @@
|
||||
/*
|
||||
* lagrange_funcs.c -- Lagrange point observation functions
|
||||
*
|
||||
* SQL functions for computing positions and observations of the five
|
||||
* Lagrange equilibrium points in the circular restricted three-body
|
||||
* problem (CR3BP). Covers Sun-planet, Earth-Moon, and planetary moon
|
||||
* systems.
|
||||
*
|
||||
* The pipeline mirrors planet_funcs.c / moon_funcs.c:
|
||||
* 1. Primary + secondary heliocentric positions (VSOP87)
|
||||
* 2. Secondary velocity via analytic (VSOP87) or central difference
|
||||
* 3. lagrange_position() -> L-point heliocentric ecliptic J2000
|
||||
* 4. Geocentric ecliptic (subtract Earth)
|
||||
* 5. observe_from_geocentric() or geocentric_to_equatorial()
|
||||
*
|
||||
* All functions are IMMUTABLE STRICT PARALLEL SAFE (compiled-in
|
||||
* mass ratios, VSOP87/ELP82B coefficients).
|
||||
*/
|
||||
|
||||
#include "postgres.h"
|
||||
#include "fmgr.h"
|
||||
#include "funcapi.h"
|
||||
#include "utils/timestamp.h"
|
||||
#include "utils/builtins.h"
|
||||
#include "types.h"
|
||||
#include "astro_math.h"
|
||||
#include "lagrange.h"
|
||||
#include "vsop87.h"
|
||||
#include "kepler.h"
|
||||
#include "elp82b.h"
|
||||
#include "l12.h"
|
||||
#include "tass17.h"
|
||||
#include "gust86.h"
|
||||
#include "marssat.h"
|
||||
#include <math.h>
|
||||
|
||||
/* ── Forward declarations ──────────────────────────────── */
|
||||
|
||||
/* Sun-planet */
|
||||
PG_FUNCTION_INFO_V1(lagrange_heliocentric);
|
||||
PG_FUNCTION_INFO_V1(lagrange_observe);
|
||||
PG_FUNCTION_INFO_V1(lagrange_equatorial);
|
||||
PG_FUNCTION_INFO_V1(lagrange_distance);
|
||||
PG_FUNCTION_INFO_V1(lagrange_distance_oe);
|
||||
|
||||
/* Earth-Moon */
|
||||
PG_FUNCTION_INFO_V1(lunar_lagrange_observe);
|
||||
PG_FUNCTION_INFO_V1(lunar_lagrange_equatorial);
|
||||
|
||||
/* Planetary moons */
|
||||
PG_FUNCTION_INFO_V1(galilean_lagrange_observe);
|
||||
PG_FUNCTION_INFO_V1(galilean_lagrange_equatorial);
|
||||
PG_FUNCTION_INFO_V1(saturn_moon_lagrange_observe);
|
||||
PG_FUNCTION_INFO_V1(saturn_moon_lagrange_equatorial);
|
||||
PG_FUNCTION_INFO_V1(uranus_moon_lagrange_observe);
|
||||
PG_FUNCTION_INFO_V1(uranus_moon_lagrange_equatorial);
|
||||
PG_FUNCTION_INFO_V1(mars_moon_lagrange_observe);
|
||||
PG_FUNCTION_INFO_V1(mars_moon_lagrange_equatorial);
|
||||
|
||||
/* Hill/zone/convenience */
|
||||
PG_FUNCTION_INFO_V1(hill_radius_func);
|
||||
PG_FUNCTION_INFO_V1(hill_radius_lunar);
|
||||
PG_FUNCTION_INFO_V1(lagrange_zone_radius_func);
|
||||
PG_FUNCTION_INFO_V1(lagrange_mass_ratio);
|
||||
PG_FUNCTION_INFO_V1(lagrange_point_name);
|
||||
|
||||
|
||||
/* ── Internal helpers ──────────────────────────────────── */
|
||||
|
||||
/*
|
||||
* Validate body_id and point_id common to Sun-planet functions.
|
||||
*/
|
||||
static void
|
||||
validate_lagrange_args(int body_id, int point_id, const char *func_name)
|
||||
{
|
||||
if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("%s: body_id %d must be 1-8 (Mercury-Neptune)",
|
||||
func_name, body_id)));
|
||||
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("%s: point_id %d must be 1-5 (L1-L5)",
|
||||
func_name, point_id)));
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Compute Sun-planet Lagrange point in heliocentric ecliptic J2000 (AU).
|
||||
*
|
||||
* Uses VSOP87 for planet position and velocity. The Sun is at (0,0,0).
|
||||
*/
|
||||
static int
|
||||
compute_sun_planet_lagrange(int body_id, int point_id, double jd,
|
||||
double result[3])
|
||||
{
|
||||
double planet_xyz[6];
|
||||
double sun[3] = {0.0, 0.0, 0.0};
|
||||
double ratio, mu;
|
||||
int vsop_body = body_id - 1;
|
||||
|
||||
GetVsop87Coor(jd, vsop_body, planet_xyz);
|
||||
|
||||
ratio = sun_planet_ratio(body_id);
|
||||
mu = mu_from_ratio(ratio);
|
||||
|
||||
/* planet_xyz[3..5] is velocity in AU/day from VSOP87 */
|
||||
return lagrange_position(sun, planet_xyz, &planet_xyz[3], mu, point_id,
|
||||
result);
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Compute Earth-Moon Lagrange point in geocentric ecliptic J2000 (AU).
|
||||
*
|
||||
* ELP2000-82B provides position only. Velocity via central difference
|
||||
* (dt=0.001 day ~ 86 seconds). Error ~40 km at lunar distance — negligible
|
||||
* vs libration zone scale (~60,000 km for L4/L5).
|
||||
*
|
||||
* Returns position relative to Earth (geocentric), not heliocentric.
|
||||
*/
|
||||
static int
|
||||
compute_lunar_lagrange(int point_id, double jd, double result_geo[3])
|
||||
{
|
||||
double moon_pos[3], moon_fwd[3], moon_bwd[3];
|
||||
double moon_vel[3];
|
||||
double earth[3] = {0.0, 0.0, 0.0}; /* geocentric frame: Earth at origin */
|
||||
double mu;
|
||||
double dt = 0.001; /* days */
|
||||
int rc;
|
||||
|
||||
/* Moon geocentric position */
|
||||
GetElp82bCoor(jd, moon_pos);
|
||||
|
||||
/* Velocity via central difference */
|
||||
GetElp82bCoor(jd + dt, moon_fwd);
|
||||
GetElp82bCoor(jd - dt, moon_bwd);
|
||||
moon_vel[0] = (moon_fwd[0] - moon_bwd[0]) / (2.0 * dt);
|
||||
moon_vel[1] = (moon_fwd[1] - moon_bwd[1]) / (2.0 * dt);
|
||||
moon_vel[2] = (moon_fwd[2] - moon_bwd[2]) / (2.0 * dt);
|
||||
|
||||
mu = mu_from_ratio(EARTH_MOON_EMRAT);
|
||||
|
||||
rc = lagrange_position(earth, moon_pos, moon_vel, mu, point_id,
|
||||
result_geo);
|
||||
return rc;
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Compute planet-moon Lagrange point heliocentric ecliptic J2000 (AU).
|
||||
*
|
||||
* The moon theory provides both position and velocity relative to the
|
||||
* parent planet. We work in a frame centered on the planet.
|
||||
*/
|
||||
static int
|
||||
compute_planetary_moon_lagrange(const double moon_rel[3],
|
||||
const double moon_vel[3],
|
||||
char family, int moon_id,
|
||||
int point_id,
|
||||
double lp_rel[3])
|
||||
{
|
||||
double planet_origin[3] = {0.0, 0.0, 0.0};
|
||||
double ratio, mu;
|
||||
|
||||
ratio = planet_moon_ratio(family, moon_id);
|
||||
if (ratio < 0.0)
|
||||
return -1;
|
||||
|
||||
mu = mu_from_ratio(ratio);
|
||||
|
||||
return lagrange_position(planet_origin, moon_rel, moon_vel, mu,
|
||||
point_id, lp_rel);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* lagrange_heliocentric(body_id int4, point_id int4, timestamptz)
|
||||
* -> heliocentric
|
||||
*
|
||||
* Sun-planet Lagrange point in heliocentric ecliptic J2000.
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
lagrange_heliocentric(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
int64 ts = PG_GETARG_INT64(2);
|
||||
double jd;
|
||||
double lp[3];
|
||||
pg_heliocentric *result;
|
||||
|
||||
validate_lagrange_args(body_id, point_id, "lagrange_heliocentric");
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
|
||||
if (compute_sun_planet_lagrange(body_id, point_id, jd, lp) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("lagrange_heliocentric: solver failed for body %d L%d",
|
||||
body_id, point_id)));
|
||||
|
||||
result = (pg_heliocentric *) palloc(sizeof(pg_heliocentric));
|
||||
result->x = lp[0];
|
||||
result->y = lp[1];
|
||||
result->z = lp[2];
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* lagrange_observe(body_id, point_id, observer, timestamptz)
|
||||
* -> topocentric
|
||||
*
|
||||
* Full pipeline: L-point heliocentric -> geocentric -> az/el.
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
lagrange_observe(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2);
|
||||
int64 ts = PG_GETARG_INT64(3);
|
||||
double jd;
|
||||
double lp[3], earth_xyz[6], geo_ecl[3];
|
||||
pg_topocentric *result;
|
||||
|
||||
validate_lagrange_args(body_id, point_id, "lagrange_observe");
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
|
||||
if (compute_sun_planet_lagrange(body_id, point_id, jd, lp) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("lagrange_observe: solver failed")));
|
||||
|
||||
/* Geocentric */
|
||||
GetVsop87Coor(jd, 2, earth_xyz);
|
||||
geo_ecl[0] = lp[0] - earth_xyz[0];
|
||||
geo_ecl[1] = lp[1] - earth_xyz[1];
|
||||
geo_ecl[2] = lp[2] - earth_xyz[2];
|
||||
|
||||
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
|
||||
observe_from_geocentric(geo_ecl, jd, obs, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* lagrange_equatorial(body_id, point_id, timestamptz) -> equatorial
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
lagrange_equatorial(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
int64 ts = PG_GETARG_INT64(2);
|
||||
double jd;
|
||||
double lp[3], earth_xyz[6], geo_ecl[3];
|
||||
pg_equatorial *result;
|
||||
|
||||
validate_lagrange_args(body_id, point_id, "lagrange_equatorial");
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
|
||||
if (compute_sun_planet_lagrange(body_id, point_id, jd, lp) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("lagrange_equatorial: solver failed")));
|
||||
|
||||
GetVsop87Coor(jd, 2, earth_xyz);
|
||||
geo_ecl[0] = lp[0] - earth_xyz[0];
|
||||
geo_ecl[1] = lp[1] - earth_xyz[1];
|
||||
geo_ecl[2] = lp[2] - earth_xyz[2];
|
||||
|
||||
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
|
||||
geocentric_to_equatorial(geo_ecl, jd, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* lagrange_distance(body_id, point_id, heliocentric, timestamptz)
|
||||
* -> float8
|
||||
*
|
||||
* Distance from a heliocentric position to a Sun-planet L-point (AU).
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
lagrange_distance(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
pg_heliocentric *pos = (pg_heliocentric *) PG_GETARG_POINTER(2);
|
||||
int64 ts = PG_GETARG_INT64(3);
|
||||
double jd;
|
||||
double lp[3];
|
||||
double dx, dy, dz, dist;
|
||||
|
||||
validate_lagrange_args(body_id, point_id, "lagrange_distance");
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
|
||||
if (compute_sun_planet_lagrange(body_id, point_id, jd, lp) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("lagrange_distance: solver failed")));
|
||||
|
||||
dx = pos->x - lp[0];
|
||||
dy = pos->y - lp[1];
|
||||
dz = pos->z - lp[2];
|
||||
dist = sqrt(dx*dx + dy*dy + dz*dz);
|
||||
|
||||
PG_RETURN_FLOAT8(dist);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* lagrange_distance_oe(body_id, point_id, orbital_elements, timestamptz)
|
||||
* -> float8
|
||||
*
|
||||
* Distance from an asteroid/comet to a Sun-planet L-point (AU).
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
lagrange_distance_oe(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
pg_orbital_elements *oe = (pg_orbital_elements *) PG_GETARG_POINTER(2);
|
||||
int64 ts = PG_GETARG_INT64(3);
|
||||
double jd;
|
||||
double lp[3], body_pos[3];
|
||||
double dx, dy, dz, dist;
|
||||
|
||||
validate_lagrange_args(body_id, point_id, "lagrange_distance_oe");
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
|
||||
if (compute_sun_planet_lagrange(body_id, point_id, jd, lp) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("lagrange_distance_oe: solver failed")));
|
||||
|
||||
kepler_position(oe->q, oe->e, oe->inc, oe->arg_peri, oe->raan,
|
||||
oe->tp, jd, body_pos);
|
||||
|
||||
dx = body_pos[0] - lp[0];
|
||||
dy = body_pos[1] - lp[1];
|
||||
dz = body_pos[2] - lp[2];
|
||||
dist = sqrt(dx*dx + dy*dy + dz*dz);
|
||||
|
||||
PG_RETURN_FLOAT8(dist);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* lunar_lagrange_observe(point_id, observer, timestamptz) -> topocentric
|
||||
*
|
||||
* Earth-Moon Lagrange point observation.
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
lunar_lagrange_observe(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 point_id = PG_GETARG_INT32(0);
|
||||
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(1);
|
||||
int64 ts = PG_GETARG_INT64(2);
|
||||
double jd;
|
||||
double lp_geo[3];
|
||||
pg_topocentric *result;
|
||||
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("lunar_lagrange_observe: point_id %d must be 1-5",
|
||||
point_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
|
||||
if (compute_lunar_lagrange(point_id, jd, lp_geo) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("lunar_lagrange_observe: solver failed")));
|
||||
|
||||
/* lp_geo is already geocentric ecliptic J2000 (AU) */
|
||||
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
|
||||
observe_from_geocentric(lp_geo, jd, obs, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* lunar_lagrange_equatorial(point_id, timestamptz) -> equatorial
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
lunar_lagrange_equatorial(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 point_id = PG_GETARG_INT32(0);
|
||||
int64 ts = PG_GETARG_INT64(1);
|
||||
double jd;
|
||||
double lp_geo[3];
|
||||
pg_equatorial *result;
|
||||
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("lunar_lagrange_equatorial: point_id %d must be 1-5",
|
||||
point_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
|
||||
if (compute_lunar_lagrange(point_id, jd, lp_geo) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("lunar_lagrange_equatorial: solver failed")));
|
||||
|
||||
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
|
||||
geocentric_to_equatorial(lp_geo, jd, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* Planetary moon Lagrange points
|
||||
*
|
||||
* Each family duplicates the observe/equatorial static helpers from
|
||||
* moon_funcs.c (per no-cross-TU-symbol convention), but routes through
|
||||
* the Lagrange solver instead of returning the moon position directly.
|
||||
* ================================================================
|
||||
*/
|
||||
|
||||
/*
|
||||
* Internal: observe a planetary moon Lagrange point.
|
||||
* L-point offset is relative to planet, same as moon offset.
|
||||
*/
|
||||
static void
|
||||
observe_pmoon_lagrange(const double lp_rel[3], int vsop_parent,
|
||||
double jd, const pg_observer *obs,
|
||||
pg_topocentric *result)
|
||||
{
|
||||
double parent_xyz[6];
|
||||
double earth_xyz[6];
|
||||
double geo_ecl[3];
|
||||
|
||||
GetVsop87Coor(jd, vsop_parent, parent_xyz);
|
||||
GetVsop87Coor(jd, 2, earth_xyz);
|
||||
|
||||
geo_ecl[0] = (parent_xyz[0] + lp_rel[0]) - earth_xyz[0];
|
||||
geo_ecl[1] = (parent_xyz[1] + lp_rel[1]) - earth_xyz[1];
|
||||
geo_ecl[2] = (parent_xyz[2] + lp_rel[2]) - earth_xyz[2];
|
||||
|
||||
observe_from_geocentric(geo_ecl, jd, obs, result);
|
||||
}
|
||||
|
||||
static void
|
||||
equatorial_pmoon_lagrange(const double lp_rel[3], int vsop_parent,
|
||||
double jd, pg_equatorial *result)
|
||||
{
|
||||
double parent_xyz[6];
|
||||
double earth_xyz[6];
|
||||
double geo_ecl[3];
|
||||
|
||||
GetVsop87Coor(jd, vsop_parent, parent_xyz);
|
||||
GetVsop87Coor(jd, 2, earth_xyz);
|
||||
|
||||
geo_ecl[0] = (parent_xyz[0] + lp_rel[0]) - earth_xyz[0];
|
||||
geo_ecl[1] = (parent_xyz[1] + lp_rel[1]) - earth_xyz[1];
|
||||
geo_ecl[2] = (parent_xyz[2] + lp_rel[2]) - earth_xyz[2];
|
||||
|
||||
geocentric_to_equatorial(geo_ecl, jd, result);
|
||||
}
|
||||
|
||||
|
||||
/* ── Galilean moons ────────────────────────────────────── */
|
||||
|
||||
Datum
|
||||
galilean_lagrange_observe(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2);
|
||||
int64 ts = PG_GETARG_INT64(3);
|
||||
double jd;
|
||||
double moon_xyz[3], moon_vel[3], lp_rel[3];
|
||||
pg_topocentric *result;
|
||||
|
||||
if (body_id < L12_IO || body_id > L12_CALLISTO)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("galilean_lagrange_observe: body_id %d must be 0-3",
|
||||
body_id)));
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("galilean_lagrange_observe: point_id %d must be 1-5",
|
||||
point_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetL12Coor(jd, body_id, moon_xyz, moon_vel);
|
||||
|
||||
if (compute_planetary_moon_lagrange(moon_xyz, moon_vel, 'g', body_id,
|
||||
point_id, lp_rel) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("galilean_lagrange_observe: solver failed")));
|
||||
|
||||
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
|
||||
observe_pmoon_lagrange(lp_rel, 4, jd, obs, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
Datum
|
||||
galilean_lagrange_equatorial(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
int64 ts = PG_GETARG_INT64(2);
|
||||
double jd;
|
||||
double moon_xyz[3], moon_vel[3], lp_rel[3];
|
||||
pg_equatorial *result;
|
||||
|
||||
if (body_id < L12_IO || body_id > L12_CALLISTO)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("galilean_lagrange_equatorial: body_id %d must be 0-3",
|
||||
body_id)));
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("galilean_lagrange_equatorial: point_id %d must be 1-5",
|
||||
point_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetL12Coor(jd, body_id, moon_xyz, moon_vel);
|
||||
|
||||
if (compute_planetary_moon_lagrange(moon_xyz, moon_vel, 'g', body_id,
|
||||
point_id, lp_rel) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("galilean_lagrange_equatorial: solver failed")));
|
||||
|
||||
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
|
||||
equatorial_pmoon_lagrange(lp_rel, 4, jd, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ── Saturn moons ──────────────────────────────────────── */
|
||||
|
||||
Datum
|
||||
saturn_moon_lagrange_observe(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2);
|
||||
int64 ts = PG_GETARG_INT64(3);
|
||||
double jd;
|
||||
double moon_xyz[3], moon_vel[3], lp_rel[3];
|
||||
pg_topocentric *result;
|
||||
|
||||
if (body_id < TASS17_MIMAS || body_id > TASS17_HYPERION)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("saturn_moon_lagrange_observe: body_id %d must be 0-7",
|
||||
body_id)));
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("saturn_moon_lagrange_observe: point_id %d must be 1-5",
|
||||
point_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetTass17Coor(jd, body_id, moon_xyz, moon_vel);
|
||||
|
||||
if (compute_planetary_moon_lagrange(moon_xyz, moon_vel, 's', body_id,
|
||||
point_id, lp_rel) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("saturn_moon_lagrange_observe: solver failed")));
|
||||
|
||||
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
|
||||
observe_pmoon_lagrange(lp_rel, 5, jd, obs, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
Datum
|
||||
saturn_moon_lagrange_equatorial(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
int64 ts = PG_GETARG_INT64(2);
|
||||
double jd;
|
||||
double moon_xyz[3], moon_vel[3], lp_rel[3];
|
||||
pg_equatorial *result;
|
||||
|
||||
if (body_id < TASS17_MIMAS || body_id > TASS17_HYPERION)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("saturn_moon_lagrange_equatorial: body_id %d must be 0-7",
|
||||
body_id)));
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("saturn_moon_lagrange_equatorial: point_id %d must be 1-5",
|
||||
point_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetTass17Coor(jd, body_id, moon_xyz, moon_vel);
|
||||
|
||||
if (compute_planetary_moon_lagrange(moon_xyz, moon_vel, 's', body_id,
|
||||
point_id, lp_rel) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("saturn_moon_lagrange_equatorial: solver failed")));
|
||||
|
||||
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
|
||||
equatorial_pmoon_lagrange(lp_rel, 5, jd, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ── Uranus moons ──────────────────────────────────────── */
|
||||
|
||||
Datum
|
||||
uranus_moon_lagrange_observe(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2);
|
||||
int64 ts = PG_GETARG_INT64(3);
|
||||
double jd;
|
||||
double moon_xyz[3], moon_vel[3], lp_rel[3];
|
||||
pg_topocentric *result;
|
||||
|
||||
if (body_id < GUST86_MIRANDA || body_id > GUST86_OBERON)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("uranus_moon_lagrange_observe: body_id %d must be 0-4",
|
||||
body_id)));
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("uranus_moon_lagrange_observe: point_id %d must be 1-5",
|
||||
point_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetGust86Coor(jd, body_id, moon_xyz, moon_vel);
|
||||
|
||||
if (compute_planetary_moon_lagrange(moon_xyz, moon_vel, 'u', body_id,
|
||||
point_id, lp_rel) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("uranus_moon_lagrange_observe: solver failed")));
|
||||
|
||||
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
|
||||
observe_pmoon_lagrange(lp_rel, 6, jd, obs, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
Datum
|
||||
uranus_moon_lagrange_equatorial(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
int64 ts = PG_GETARG_INT64(2);
|
||||
double jd;
|
||||
double moon_xyz[3], moon_vel[3], lp_rel[3];
|
||||
pg_equatorial *result;
|
||||
|
||||
if (body_id < GUST86_MIRANDA || body_id > GUST86_OBERON)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("uranus_moon_lagrange_equatorial: body_id %d must be 0-4",
|
||||
body_id)));
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("uranus_moon_lagrange_equatorial: point_id %d must be 1-5",
|
||||
point_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetGust86Coor(jd, body_id, moon_xyz, moon_vel);
|
||||
|
||||
if (compute_planetary_moon_lagrange(moon_xyz, moon_vel, 'u', body_id,
|
||||
point_id, lp_rel) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("uranus_moon_lagrange_equatorial: solver failed")));
|
||||
|
||||
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
|
||||
equatorial_pmoon_lagrange(lp_rel, 6, jd, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ── Mars moons ────────────────────────────────────────── */
|
||||
|
||||
Datum
|
||||
mars_moon_lagrange_observe(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(2);
|
||||
int64 ts = PG_GETARG_INT64(3);
|
||||
double jd;
|
||||
double moon_xyz[3], moon_vel[3], lp_rel[3];
|
||||
pg_topocentric *result;
|
||||
|
||||
if (body_id < MARS_SAT_PHOBOS || body_id > MARS_SAT_DEIMOS)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("mars_moon_lagrange_observe: body_id %d must be 0-1",
|
||||
body_id)));
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("mars_moon_lagrange_observe: point_id %d must be 1-5",
|
||||
point_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetMarsSatCoor(jd, body_id, moon_xyz, moon_vel);
|
||||
|
||||
if (compute_planetary_moon_lagrange(moon_xyz, moon_vel, 'm', body_id,
|
||||
point_id, lp_rel) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("mars_moon_lagrange_observe: solver failed")));
|
||||
|
||||
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
|
||||
observe_pmoon_lagrange(lp_rel, 3, jd, obs, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
Datum
|
||||
mars_moon_lagrange_equatorial(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
int64 ts = PG_GETARG_INT64(2);
|
||||
double jd;
|
||||
double moon_xyz[3], moon_vel[3], lp_rel[3];
|
||||
pg_equatorial *result;
|
||||
|
||||
if (body_id < MARS_SAT_PHOBOS || body_id > MARS_SAT_DEIMOS)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("mars_moon_lagrange_equatorial: body_id %d must be 0-1",
|
||||
body_id)));
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("mars_moon_lagrange_equatorial: point_id %d must be 1-5",
|
||||
point_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetMarsSatCoor(jd, body_id, moon_xyz, moon_vel);
|
||||
|
||||
if (compute_planetary_moon_lagrange(moon_xyz, moon_vel, 'm', body_id,
|
||||
point_id, lp_rel) != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INTERNAL_ERROR),
|
||||
errmsg("mars_moon_lagrange_equatorial: solver failed")));
|
||||
|
||||
result = (pg_equatorial *) palloc(sizeof(pg_equatorial));
|
||||
equatorial_pmoon_lagrange(lp_rel, 3, jd, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* Hill radius / zone / convenience functions
|
||||
* ================================================================
|
||||
*/
|
||||
|
||||
/* hill_radius(body_id int4, timestamptz) -> float8 (AU) */
|
||||
Datum
|
||||
hill_radius_func(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int64 ts = PG_GETARG_INT64(1);
|
||||
double jd;
|
||||
double planet_xyz[6], sep, ratio, mu;
|
||||
|
||||
if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("hill_radius: body_id %d must be 1-8", body_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetVsop87Coor(jd, body_id - 1, planet_xyz);
|
||||
|
||||
sep = sqrt(planet_xyz[0]*planet_xyz[0] +
|
||||
planet_xyz[1]*planet_xyz[1] +
|
||||
planet_xyz[2]*planet_xyz[2]);
|
||||
|
||||
ratio = sun_planet_ratio(body_id);
|
||||
mu = mu_from_ratio(ratio);
|
||||
|
||||
PG_RETURN_FLOAT8(lagrange_hill_radius(sep, mu));
|
||||
}
|
||||
|
||||
|
||||
/* hill_radius_lunar(timestamptz) -> float8 (AU) */
|
||||
Datum
|
||||
hill_radius_lunar(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int64 ts = PG_GETARG_INT64(0);
|
||||
double jd;
|
||||
double moon_pos[3], sep, mu;
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetElp82bCoor(jd, moon_pos);
|
||||
|
||||
sep = sqrt(moon_pos[0]*moon_pos[0] +
|
||||
moon_pos[1]*moon_pos[1] +
|
||||
moon_pos[2]*moon_pos[2]);
|
||||
|
||||
mu = mu_from_ratio(EARTH_MOON_EMRAT);
|
||||
|
||||
PG_RETURN_FLOAT8(lagrange_hill_radius(sep, mu));
|
||||
}
|
||||
|
||||
|
||||
/* lagrange_zone_radius(body_id, point_id, timestamptz) -> float8 (AU) */
|
||||
Datum
|
||||
lagrange_zone_radius_func(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int32 point_id = PG_GETARG_INT32(1);
|
||||
int64 ts = PG_GETARG_INT64(2);
|
||||
double jd;
|
||||
double planet_xyz[6], sep, ratio, mu, zr;
|
||||
|
||||
validate_lagrange_args(body_id, point_id, "lagrange_zone_radius");
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetVsop87Coor(jd, body_id - 1, planet_xyz);
|
||||
|
||||
sep = sqrt(planet_xyz[0]*planet_xyz[0] +
|
||||
planet_xyz[1]*planet_xyz[1] +
|
||||
planet_xyz[2]*planet_xyz[2]);
|
||||
|
||||
ratio = sun_planet_ratio(body_id);
|
||||
mu = mu_from_ratio(ratio);
|
||||
|
||||
zr = lagrange_zone_radius(sep, mu, point_id);
|
||||
|
||||
PG_RETURN_FLOAT8(zr);
|
||||
}
|
||||
|
||||
|
||||
/* lagrange_mass_ratio(body_id int4) -> float8 */
|
||||
Datum
|
||||
lagrange_mass_ratio(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
double ratio;
|
||||
|
||||
if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("lagrange_mass_ratio: body_id %d must be 1-8",
|
||||
body_id)));
|
||||
|
||||
ratio = sun_planet_ratio(body_id);
|
||||
|
||||
PG_RETURN_FLOAT8(mu_from_ratio(ratio));
|
||||
}
|
||||
|
||||
|
||||
/* lagrange_point_name(point_id int4) -> text */
|
||||
Datum
|
||||
lagrange_point_name(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 point_id = PG_GETARG_INT32(0);
|
||||
|
||||
if (point_id < LAGRANGE_L1 || point_id > LAGRANGE_L5)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("lagrange_point_name: point_id %d must be 1-5",
|
||||
point_id)));
|
||||
|
||||
switch (point_id)
|
||||
{
|
||||
case LAGRANGE_L1: PG_RETURN_TEXT_P(cstring_to_text("L1"));
|
||||
case LAGRANGE_L2: PG_RETURN_TEXT_P(cstring_to_text("L2"));
|
||||
case LAGRANGE_L3: PG_RETURN_TEXT_P(cstring_to_text("L3"));
|
||||
case LAGRANGE_L4: PG_RETURN_TEXT_P(cstring_to_text("L4"));
|
||||
case LAGRANGE_L5: PG_RETURN_TEXT_P(cstring_to_text("L5"));
|
||||
default: PG_RETURN_NULL();
|
||||
}
|
||||
}
|
||||
323
test/expected/v020_features.out
Normal file
323
test/expected/v020_features.out
Normal file
@ -0,0 +1,323 @@
|
||||
-- v020_features: Lagrange point support
|
||||
-- Tests Sun-planet, Earth-Moon, planetary moon Lagrange points,
|
||||
-- Hill radius, zone radius, DE fallback, and input validation.
|
||||
-- Reference observer: Greenwich, UK
|
||||
\set obs '''(51.4769,-0.0005,0)'''
|
||||
-- Reference time: J2000 epoch (2000-01-01 12:00:00 UTC)
|
||||
\set t '''2000-01-01 12:00:00+00'''
|
||||
-- ============================================================
|
||||
-- Sun-Earth L1/L2: should be ~0.01 AU from Earth (~1.5 million km)
|
||||
-- SOHO is at L1, JWST at L2.
|
||||
-- ============================================================
|
||||
-- L1 heliocentric: should be close to Earth's heliocentric (~1 AU from Sun)
|
||||
SELECT
|
||||
round(helio_distance(lagrange_heliocentric(3, 1, :t ::timestamptz))::numeric, 2) AS sun_dist_au;
|
||||
sun_dist_au
|
||||
-------------
|
||||
0.97
|
||||
(1 row)
|
||||
|
||||
-- L2 heliocentric: also ~1 AU from Sun, slightly further than L1
|
||||
SELECT
|
||||
round(helio_distance(lagrange_heliocentric(3, 2, :t ::timestamptz))::numeric, 2) AS sun_dist_au;
|
||||
sun_dist_au
|
||||
-------------
|
||||
0.99
|
||||
(1 row)
|
||||
|
||||
-- L1 between Sun and Earth (closer to Sun than L2)
|
||||
SELECT
|
||||
helio_distance(lagrange_heliocentric(3, 1, :t ::timestamptz))
|
||||
<
|
||||
helio_distance(lagrange_heliocentric(3, 2, :t ::timestamptz))
|
||||
AS l1_closer_than_l2;
|
||||
l1_closer_than_l2
|
||||
-------------------
|
||||
t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Sun-Jupiter L4/L5: ~60 degrees from Jupiter, ~5.2 AU from Sun
|
||||
-- These are the Trojan asteroid zones.
|
||||
-- ============================================================
|
||||
-- L4/L5 should be ~5.2 AU from Sun
|
||||
SELECT
|
||||
round(helio_distance(lagrange_heliocentric(5, 4, :t ::timestamptz))::numeric, 1) AS l4_sun_dist;
|
||||
l4_sun_dist
|
||||
-------------
|
||||
5.0
|
||||
(1 row)
|
||||
|
||||
SELECT
|
||||
round(helio_distance(lagrange_heliocentric(5, 5, :t ::timestamptz))::numeric, 1) AS l5_sun_dist;
|
||||
l5_sun_dist
|
||||
-------------
|
||||
5.0
|
||||
(1 row)
|
||||
|
||||
-- L4 and L5 equidistant from Sun (within 0.001 AU)
|
||||
SELECT
|
||||
abs(
|
||||
helio_distance(lagrange_heliocentric(5, 4, :t ::timestamptz))
|
||||
-
|
||||
helio_distance(lagrange_heliocentric(5, 5, :t ::timestamptz))
|
||||
) < 0.001 AS l4_l5_equidistant;
|
||||
l4_l5_equidistant
|
||||
-------------------
|
||||
t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Earth-Moon L1: ~326,000 km from Earth
|
||||
-- ============================================================
|
||||
-- lunar_lagrange_equatorial returns distance in km
|
||||
SELECT
|
||||
round(eq_distance(lunar_lagrange_equatorial(1, :t ::timestamptz))::numeric, -3)
|
||||
BETWEEN 300000 AND 360000 AS em_l1_in_range;
|
||||
em_l1_in_range
|
||||
----------------
|
||||
t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- lagrange_observe returns valid az/el
|
||||
-- ============================================================
|
||||
SELECT
|
||||
topo_elevation(lagrange_observe(3, 2, :obs ::observer, :t ::timestamptz))
|
||||
BETWEEN -90 AND 90 AS valid_elevation;
|
||||
valid_elevation
|
||||
-----------------
|
||||
t
|
||||
(1 row)
|
||||
|
||||
-- lagrange_equatorial returns valid RA/Dec
|
||||
SELECT
|
||||
eq_ra(lagrange_equatorial(3, 1, :t ::timestamptz)) BETWEEN 0 AND 24 AS valid_ra,
|
||||
eq_dec(lagrange_equatorial(3, 1, :t ::timestamptz)) BETWEEN -90 AND 90 AS valid_dec;
|
||||
valid_ra | valid_dec
|
||||
----------+-----------
|
||||
t | t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- lagrange_distance self-test: L-point distance to itself ≈ 0
|
||||
-- ============================================================
|
||||
SELECT
|
||||
round(lagrange_distance(
|
||||
5, 4,
|
||||
lagrange_heliocentric(5, 4, :t ::timestamptz),
|
||||
:t ::timestamptz
|
||||
)::numeric, 10) AS self_distance;
|
||||
self_distance
|
||||
---------------
|
||||
0.0000000000
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Hill radius
|
||||
-- ============================================================
|
||||
-- Jupiter Hill radius ~0.35 AU
|
||||
SELECT
|
||||
round(hill_radius(5, :t ::timestamptz)::numeric, 2)
|
||||
BETWEEN 0.30 AND 0.40 AS jupiter_hill_ok;
|
||||
jupiter_hill_ok
|
||||
-----------------
|
||||
t
|
||||
(1 row)
|
||||
|
||||
-- Earth Hill radius ~0.01 AU
|
||||
SELECT
|
||||
round(hill_radius(3, :t ::timestamptz)::numeric, 3)
|
||||
BETWEEN 0.008 AND 0.012 AS earth_hill_ok;
|
||||
earth_hill_ok
|
||||
---------------
|
||||
t
|
||||
(1 row)
|
||||
|
||||
-- Lunar Hill radius (much smaller, AU)
|
||||
SELECT
|
||||
hill_radius_lunar(:t ::timestamptz) > 0 AS lunar_hill_positive;
|
||||
lunar_hill_positive
|
||||
---------------------
|
||||
t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Zone radius
|
||||
-- ============================================================
|
||||
SELECT
|
||||
lagrange_zone_radius(5, 4, :t ::timestamptz) > 0 AS jup_l4_zone_positive;
|
||||
jup_l4_zone_positive
|
||||
----------------------
|
||||
t
|
||||
(1 row)
|
||||
|
||||
SELECT
|
||||
lagrange_zone_radius(5, 1, :t ::timestamptz) > 0 AS jup_l1_zone_positive;
|
||||
jup_l1_zone_positive
|
||||
----------------------
|
||||
t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Convenience functions
|
||||
-- ============================================================
|
||||
-- lagrange_mass_ratio returns small positive number
|
||||
SELECT
|
||||
lagrange_mass_ratio(5) > 0 AND lagrange_mass_ratio(5) < 0.01 AS jupiter_mu_ok;
|
||||
jupiter_mu_ok
|
||||
---------------
|
||||
t
|
||||
(1 row)
|
||||
|
||||
SELECT
|
||||
lagrange_mass_ratio(3) > 0 AND lagrange_mass_ratio(3) < 0.001 AS earth_mu_ok;
|
||||
earth_mu_ok
|
||||
-------------
|
||||
t
|
||||
(1 row)
|
||||
|
||||
-- lagrange_point_name
|
||||
SELECT lagrange_point_name(1) AS l1_name;
|
||||
l1_name
|
||||
---------
|
||||
L1
|
||||
(1 row)
|
||||
|
||||
SELECT lagrange_point_name(5) AS l5_name;
|
||||
l5_name
|
||||
---------
|
||||
L5
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- All planets produce valid results
|
||||
-- ============================================================
|
||||
SELECT body_id,
|
||||
round(helio_distance(lagrange_heliocentric(body_id, 1, :t ::timestamptz))::numeric, 2) AS sun_dist_au
|
||||
FROM generate_series(1, 8) AS body_id
|
||||
ORDER BY body_id;
|
||||
body_id | sun_dist_au
|
||||
---------+-------------
|
||||
1 | 0.46
|
||||
2 | 0.71
|
||||
3 | 0.97
|
||||
4 | 1.38
|
||||
5 | 4.63
|
||||
6 | 8.77
|
||||
7 | 19.44
|
||||
8 | 29.35
|
||||
(8 rows)
|
||||
|
||||
-- ============================================================
|
||||
-- Planetary moon Lagrange points
|
||||
-- ============================================================
|
||||
-- Galilean: Io L4 (body=0, point=4)
|
||||
SELECT
|
||||
eq_ra(galilean_lagrange_equatorial(0, 4, :t ::timestamptz)) BETWEEN 0 AND 24
|
||||
AS io_l4_valid_ra;
|
||||
io_l4_valid_ra
|
||||
----------------
|
||||
t
|
||||
(1 row)
|
||||
|
||||
-- Saturn: Titan L1 (body=5, point=1)
|
||||
SELECT
|
||||
eq_ra(saturn_moon_lagrange_equatorial(5, 1, :t ::timestamptz)) BETWEEN 0 AND 24
|
||||
AS titan_l1_valid_ra;
|
||||
titan_l1_valid_ra
|
||||
-------------------
|
||||
t
|
||||
(1 row)
|
||||
|
||||
-- Uranus: Titania L2 (body=3, point=2)
|
||||
SELECT
|
||||
eq_ra(uranus_moon_lagrange_equatorial(3, 2, :t ::timestamptz)) BETWEEN 0 AND 24
|
||||
AS titania_l2_valid_ra;
|
||||
titania_l2_valid_ra
|
||||
---------------------
|
||||
t
|
||||
(1 row)
|
||||
|
||||
-- Mars: Phobos L5 (body=0, point=5)
|
||||
SELECT
|
||||
eq_ra(mars_moon_lagrange_equatorial(0, 5, :t ::timestamptz)) BETWEEN 0 AND 24
|
||||
AS phobos_l5_valid_ra;
|
||||
phobos_l5_valid_ra
|
||||
--------------------
|
||||
t
|
||||
(1 row)
|
||||
|
||||
-- Galilean observe returns valid topocentric
|
||||
SELECT
|
||||
topo_elevation(galilean_lagrange_observe(2, 3, :obs ::observer, :t ::timestamptz))
|
||||
BETWEEN -90 AND 90 AS ganymede_l3_valid_el;
|
||||
ganymede_l3_valid_el
|
||||
----------------------
|
||||
t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- DE fallback (no DE loaded, should produce same results as VSOP87)
|
||||
-- ============================================================
|
||||
SELECT
|
||||
round(helio_distance(lagrange_heliocentric_de(3, 1, :t ::timestamptz))::numeric, 4) AS de_l1_dist;
|
||||
de_l1_dist
|
||||
------------
|
||||
0.9735
|
||||
(1 row)
|
||||
|
||||
SELECT
|
||||
round(helio_distance(lagrange_heliocentric(3, 1, :t ::timestamptz))::numeric, 4) AS vsop_l1_dist;
|
||||
vsop_l1_dist
|
||||
--------------
|
||||
0.9735
|
||||
(1 row)
|
||||
|
||||
-- DE hill_radius fallback
|
||||
SELECT
|
||||
round(hill_radius_de(5, :t ::timestamptz)::numeric, 4) =
|
||||
round(hill_radius(5, :t ::timestamptz)::numeric, 4)
|
||||
AS hill_de_matches_vsop;
|
||||
hill_de_matches_vsop
|
||||
----------------------
|
||||
t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Input validation
|
||||
-- ============================================================
|
||||
-- Bad body_id
|
||||
SELECT lagrange_heliocentric(0, 1, :t ::timestamptz); -- Sun not valid
|
||||
ERROR: lagrange_heliocentric: body_id 0 must be 1-8 (Mercury-Neptune)
|
||||
SELECT lagrange_heliocentric(9, 1, :t ::timestamptz); -- body 9 invalid
|
||||
ERROR: lagrange_heliocentric: body_id 9 must be 1-8 (Mercury-Neptune)
|
||||
-- Bad point_id
|
||||
SELECT lagrange_heliocentric(3, 0, :t ::timestamptz); -- point 0 invalid
|
||||
ERROR: lagrange_heliocentric: point_id 0 must be 1-5 (L1-L5)
|
||||
SELECT lagrange_heliocentric(3, 6, :t ::timestamptz); -- point 6 invalid
|
||||
ERROR: lagrange_heliocentric: point_id 6 must be 1-5 (L1-L5)
|
||||
-- Bad lunar point_id
|
||||
SELECT lunar_lagrange_equatorial(0, :t ::timestamptz); -- point 0 invalid
|
||||
ERROR: lunar_lagrange_equatorial: point_id 0 must be 1-5
|
||||
SELECT lunar_lagrange_equatorial(6, :t ::timestamptz); -- point 6 invalid
|
||||
ERROR: lunar_lagrange_equatorial: point_id 6 must be 1-5
|
||||
-- Bad planetary moon body_id
|
||||
SELECT galilean_lagrange_equatorial(4, 1, :t ::timestamptz); -- Galilean 4 invalid
|
||||
ERROR: galilean_lagrange_equatorial: body_id 4 must be 0-3
|
||||
SELECT saturn_moon_lagrange_equatorial(8, 1, :t ::timestamptz); -- Saturn 8 invalid
|
||||
ERROR: saturn_moon_lagrange_equatorial: body_id 8 must be 0-7
|
||||
SELECT uranus_moon_lagrange_equatorial(5, 1, :t ::timestamptz); -- Uranus 5 invalid
|
||||
ERROR: uranus_moon_lagrange_equatorial: body_id 5 must be 0-4
|
||||
SELECT mars_moon_lagrange_equatorial(2, 1, :t ::timestamptz); -- Mars 2 invalid
|
||||
ERROR: mars_moon_lagrange_equatorial: body_id 2 must be 0-1
|
||||
-- lagrange_mass_ratio bad body
|
||||
SELECT lagrange_mass_ratio(0);
|
||||
ERROR: lagrange_mass_ratio: body_id 0 must be 1-8
|
||||
SELECT lagrange_mass_ratio(9);
|
||||
ERROR: lagrange_mass_ratio: body_id 9 must be 1-8
|
||||
-- lagrange_point_name bad id
|
||||
SELECT lagrange_point_name(0);
|
||||
ERROR: lagrange_point_name: point_id 0 must be 1-5
|
||||
SELECT lagrange_point_name(6);
|
||||
ERROR: lagrange_point_name: point_id 6 must be 1-5
|
||||
209
test/sql/v020_features.sql
Normal file
209
test/sql/v020_features.sql
Normal file
@ -0,0 +1,209 @@
|
||||
-- v020_features: Lagrange point support
|
||||
-- Tests Sun-planet, Earth-Moon, planetary moon Lagrange points,
|
||||
-- Hill radius, zone radius, DE fallback, and input validation.
|
||||
|
||||
-- Reference observer: Greenwich, UK
|
||||
\set obs '''(51.4769,-0.0005,0)'''
|
||||
|
||||
-- Reference time: J2000 epoch (2000-01-01 12:00:00 UTC)
|
||||
\set t '''2000-01-01 12:00:00+00'''
|
||||
|
||||
-- ============================================================
|
||||
-- Sun-Earth L1/L2: should be ~0.01 AU from Earth (~1.5 million km)
|
||||
-- SOHO is at L1, JWST at L2.
|
||||
-- ============================================================
|
||||
|
||||
-- L1 heliocentric: should be close to Earth's heliocentric (~1 AU from Sun)
|
||||
SELECT
|
||||
round(helio_distance(lagrange_heliocentric(3, 1, :t ::timestamptz))::numeric, 2) AS sun_dist_au;
|
||||
|
||||
-- L2 heliocentric: also ~1 AU from Sun, slightly further than L1
|
||||
SELECT
|
||||
round(helio_distance(lagrange_heliocentric(3, 2, :t ::timestamptz))::numeric, 2) AS sun_dist_au;
|
||||
|
||||
-- L1 between Sun and Earth (closer to Sun than L2)
|
||||
SELECT
|
||||
helio_distance(lagrange_heliocentric(3, 1, :t ::timestamptz))
|
||||
<
|
||||
helio_distance(lagrange_heliocentric(3, 2, :t ::timestamptz))
|
||||
AS l1_closer_than_l2;
|
||||
|
||||
-- ============================================================
|
||||
-- Sun-Jupiter L4/L5: ~60 degrees from Jupiter, ~5.2 AU from Sun
|
||||
-- These are the Trojan asteroid zones.
|
||||
-- ============================================================
|
||||
|
||||
-- L4/L5 should be ~5.2 AU from Sun
|
||||
SELECT
|
||||
round(helio_distance(lagrange_heliocentric(5, 4, :t ::timestamptz))::numeric, 1) AS l4_sun_dist;
|
||||
|
||||
SELECT
|
||||
round(helio_distance(lagrange_heliocentric(5, 5, :t ::timestamptz))::numeric, 1) AS l5_sun_dist;
|
||||
|
||||
-- L4 and L5 equidistant from Sun (within 0.001 AU)
|
||||
SELECT
|
||||
abs(
|
||||
helio_distance(lagrange_heliocentric(5, 4, :t ::timestamptz))
|
||||
-
|
||||
helio_distance(lagrange_heliocentric(5, 5, :t ::timestamptz))
|
||||
) < 0.001 AS l4_l5_equidistant;
|
||||
|
||||
-- ============================================================
|
||||
-- Earth-Moon L1: ~326,000 km from Earth
|
||||
-- ============================================================
|
||||
|
||||
-- lunar_lagrange_equatorial returns distance in km
|
||||
SELECT
|
||||
round(eq_distance(lunar_lagrange_equatorial(1, :t ::timestamptz))::numeric, -3)
|
||||
BETWEEN 300000 AND 360000 AS em_l1_in_range;
|
||||
|
||||
-- ============================================================
|
||||
-- lagrange_observe returns valid az/el
|
||||
-- ============================================================
|
||||
|
||||
SELECT
|
||||
topo_elevation(lagrange_observe(3, 2, :obs ::observer, :t ::timestamptz))
|
||||
BETWEEN -90 AND 90 AS valid_elevation;
|
||||
|
||||
-- lagrange_equatorial returns valid RA/Dec
|
||||
SELECT
|
||||
eq_ra(lagrange_equatorial(3, 1, :t ::timestamptz)) BETWEEN 0 AND 24 AS valid_ra,
|
||||
eq_dec(lagrange_equatorial(3, 1, :t ::timestamptz)) BETWEEN -90 AND 90 AS valid_dec;
|
||||
|
||||
-- ============================================================
|
||||
-- lagrange_distance self-test: L-point distance to itself ≈ 0
|
||||
-- ============================================================
|
||||
|
||||
SELECT
|
||||
round(lagrange_distance(
|
||||
5, 4,
|
||||
lagrange_heliocentric(5, 4, :t ::timestamptz),
|
||||
:t ::timestamptz
|
||||
)::numeric, 10) AS self_distance;
|
||||
|
||||
-- ============================================================
|
||||
-- Hill radius
|
||||
-- ============================================================
|
||||
|
||||
-- Jupiter Hill radius ~0.35 AU
|
||||
SELECT
|
||||
round(hill_radius(5, :t ::timestamptz)::numeric, 2)
|
||||
BETWEEN 0.30 AND 0.40 AS jupiter_hill_ok;
|
||||
|
||||
-- Earth Hill radius ~0.01 AU
|
||||
SELECT
|
||||
round(hill_radius(3, :t ::timestamptz)::numeric, 3)
|
||||
BETWEEN 0.008 AND 0.012 AS earth_hill_ok;
|
||||
|
||||
-- Lunar Hill radius (much smaller, AU)
|
||||
SELECT
|
||||
hill_radius_lunar(:t ::timestamptz) > 0 AS lunar_hill_positive;
|
||||
|
||||
-- ============================================================
|
||||
-- Zone radius
|
||||
-- ============================================================
|
||||
|
||||
SELECT
|
||||
lagrange_zone_radius(5, 4, :t ::timestamptz) > 0 AS jup_l4_zone_positive;
|
||||
|
||||
SELECT
|
||||
lagrange_zone_radius(5, 1, :t ::timestamptz) > 0 AS jup_l1_zone_positive;
|
||||
|
||||
-- ============================================================
|
||||
-- Convenience functions
|
||||
-- ============================================================
|
||||
|
||||
-- lagrange_mass_ratio returns small positive number
|
||||
SELECT
|
||||
lagrange_mass_ratio(5) > 0 AND lagrange_mass_ratio(5) < 0.01 AS jupiter_mu_ok;
|
||||
|
||||
SELECT
|
||||
lagrange_mass_ratio(3) > 0 AND lagrange_mass_ratio(3) < 0.001 AS earth_mu_ok;
|
||||
|
||||
-- lagrange_point_name
|
||||
SELECT lagrange_point_name(1) AS l1_name;
|
||||
SELECT lagrange_point_name(5) AS l5_name;
|
||||
|
||||
-- ============================================================
|
||||
-- All planets produce valid results
|
||||
-- ============================================================
|
||||
|
||||
SELECT body_id,
|
||||
round(helio_distance(lagrange_heliocentric(body_id, 1, :t ::timestamptz))::numeric, 2) AS sun_dist_au
|
||||
FROM generate_series(1, 8) AS body_id
|
||||
ORDER BY body_id;
|
||||
|
||||
-- ============================================================
|
||||
-- Planetary moon Lagrange points
|
||||
-- ============================================================
|
||||
|
||||
-- Galilean: Io L4 (body=0, point=4)
|
||||
SELECT
|
||||
eq_ra(galilean_lagrange_equatorial(0, 4, :t ::timestamptz)) BETWEEN 0 AND 24
|
||||
AS io_l4_valid_ra;
|
||||
|
||||
-- Saturn: Titan L1 (body=5, point=1)
|
||||
SELECT
|
||||
eq_ra(saturn_moon_lagrange_equatorial(5, 1, :t ::timestamptz)) BETWEEN 0 AND 24
|
||||
AS titan_l1_valid_ra;
|
||||
|
||||
-- Uranus: Titania L2 (body=3, point=2)
|
||||
SELECT
|
||||
eq_ra(uranus_moon_lagrange_equatorial(3, 2, :t ::timestamptz)) BETWEEN 0 AND 24
|
||||
AS titania_l2_valid_ra;
|
||||
|
||||
-- Mars: Phobos L5 (body=0, point=5)
|
||||
SELECT
|
||||
eq_ra(mars_moon_lagrange_equatorial(0, 5, :t ::timestamptz)) BETWEEN 0 AND 24
|
||||
AS phobos_l5_valid_ra;
|
||||
|
||||
-- Galilean observe returns valid topocentric
|
||||
SELECT
|
||||
topo_elevation(galilean_lagrange_observe(2, 3, :obs ::observer, :t ::timestamptz))
|
||||
BETWEEN -90 AND 90 AS ganymede_l3_valid_el;
|
||||
|
||||
-- ============================================================
|
||||
-- DE fallback (no DE loaded, should produce same results as VSOP87)
|
||||
-- ============================================================
|
||||
|
||||
SELECT
|
||||
round(helio_distance(lagrange_heliocentric_de(3, 1, :t ::timestamptz))::numeric, 4) AS de_l1_dist;
|
||||
|
||||
SELECT
|
||||
round(helio_distance(lagrange_heliocentric(3, 1, :t ::timestamptz))::numeric, 4) AS vsop_l1_dist;
|
||||
|
||||
-- DE hill_radius fallback
|
||||
SELECT
|
||||
round(hill_radius_de(5, :t ::timestamptz)::numeric, 4) =
|
||||
round(hill_radius(5, :t ::timestamptz)::numeric, 4)
|
||||
AS hill_de_matches_vsop;
|
||||
|
||||
-- ============================================================
|
||||
-- Input validation
|
||||
-- ============================================================
|
||||
|
||||
-- Bad body_id
|
||||
SELECT lagrange_heliocentric(0, 1, :t ::timestamptz); -- Sun not valid
|
||||
SELECT lagrange_heliocentric(9, 1, :t ::timestamptz); -- body 9 invalid
|
||||
|
||||
-- Bad point_id
|
||||
SELECT lagrange_heliocentric(3, 0, :t ::timestamptz); -- point 0 invalid
|
||||
SELECT lagrange_heliocentric(3, 6, :t ::timestamptz); -- point 6 invalid
|
||||
|
||||
-- Bad lunar point_id
|
||||
SELECT lunar_lagrange_equatorial(0, :t ::timestamptz); -- point 0 invalid
|
||||
SELECT lunar_lagrange_equatorial(6, :t ::timestamptz); -- point 6 invalid
|
||||
|
||||
-- Bad planetary moon body_id
|
||||
SELECT galilean_lagrange_equatorial(4, 1, :t ::timestamptz); -- Galilean 4 invalid
|
||||
SELECT saturn_moon_lagrange_equatorial(8, 1, :t ::timestamptz); -- Saturn 8 invalid
|
||||
SELECT uranus_moon_lagrange_equatorial(5, 1, :t ::timestamptz); -- Uranus 5 invalid
|
||||
SELECT mars_moon_lagrange_equatorial(2, 1, :t ::timestamptz); -- Mars 2 invalid
|
||||
|
||||
-- lagrange_mass_ratio bad body
|
||||
SELECT lagrange_mass_ratio(0);
|
||||
SELECT lagrange_mass_ratio(9);
|
||||
|
||||
-- lagrange_point_name bad id
|
||||
SELECT lagrange_point_name(0);
|
||||
SELECT lagrange_point_name(6);
|
||||
Loading…
x
Reference in New Issue
Block a user