Add optional JPL DE440/441 ephemeris support (v0.3.0)
Clean-room DE binary reader (~400 lines C) with Chebyshev/Clenshaw evaluation — no GPL dependency on jpl_eph. Per-backend lazy initialization preserves PARALLEL SAFE. Existing VSOP87/ELP82B functions stay IMMUTABLE; new _de() variants are STABLE with automatic fallback to compiled-in ephemerides on any DE failure. Implementation: - de_reader.c: header parse, record seek, Clenshaw recurrence - eph_provider.c: GUC (pg_orbit.ephemeris_path), lazy init, ICRS-to-ecliptic frame rotation, on_proc_exit cleanup - de_funcs.c: 11 new SQL functions (_de variants + diagnostics) - Constant chain of custody rules 6-8 (frame rotation, same-provider, AU consistency) Extract observe_from_geocentric() to astro_math.h for shared use by planet_funcs.c, moon_funcs.c, and de_funcs.c. 57 → 68 functions, 11 → 12 regression test suites, all passing.
This commit is contained in:
parent
0e09975e08
commit
15fa553c0e
91
CLAUDE.md
91
CLAUDE.md
@ -1,9 +1,9 @@
|
||||
# pg_orbit — Solar System Computation for PostgreSQL
|
||||
|
||||
## What This Is
|
||||
A PostgreSQL extension that moves orbital mechanics inside the database — the way PostGIS did for geography. Native C extension using PGXS, 57 SQL functions, 7 custom types, covering satellites (SGP4/SDP4), planets (VSOP87), Moon (ELP2000-82B), 19 planetary moons (L1.2/TASS17/GUST86/MarsSat), stars, comets, Jupiter radio bursts, and interplanetary Lambert transfers.
|
||||
A PostgreSQL extension that moves orbital mechanics inside the database — the way PostGIS did for geography. Native C extension using PGXS, 68 SQL functions, 7 custom types, covering satellites (SGP4/SDP4), planets (VSOP87 + optional JPL DE441), Moon (ELP2000-82B), 19 planetary moons (L1.2/TASS17/GUST86/MarsSat), stars, comets, Jupiter radio bursts, and interplanetary Lambert transfers.
|
||||
|
||||
**Current version:** 0.2.0 on branch `phase/solar-system-expansion`
|
||||
**Current version:** 0.3.0 on branch `phase/solar-system-expansion`
|
||||
**Repository:** https://git.supported.systems/warehack.ing/pg_orbit
|
||||
**Documentation:** https://pg-orbit.warehack.ing
|
||||
|
||||
@ -11,7 +11,7 @@ A PostgreSQL extension that moves orbital mechanics inside the database — the
|
||||
```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 11 regression test suites
|
||||
make installcheck PG_CONFIG=/usr/bin/pg_config # Run 12 regression test suites
|
||||
```
|
||||
|
||||
Requires: PostgreSQL 17 development headers, GCC, G++ (for sat_code C++), Make.
|
||||
@ -27,16 +27,18 @@ Image: `git.supported.systems/warehack.ing/pg_orbit:pg17`
|
||||
|
||||
## Project Layout
|
||||
```
|
||||
pg_orbit.control # Extension metadata (version 0.2.0)
|
||||
pg_orbit.control # Extension metadata (version 0.3.0)
|
||||
Makefile # PGXS build + Docker targets
|
||||
sql/
|
||||
pg_orbit--0.1.0.sql # v0.1.0: satellite types/functions/operators
|
||||
pg_orbit--0.2.0.sql # v0.2.0: complete extension (all 57 functions)
|
||||
pg_orbit--0.1.0--0.2.0.sql # Migration script (adds solar system)
|
||||
pg_orbit--0.2.0.sql # v0.2.0: solar system (57 functions)
|
||||
pg_orbit--0.3.0.sql # v0.3.0: complete extension (68 functions)
|
||||
pg_orbit--0.1.0--0.2.0.sql # Migration: v0.1.0 → v0.2.0 (adds solar system)
|
||||
pg_orbit--0.2.0--0.3.0.sql # Migration: v0.2.0 → v0.3.0 (adds DE ephemeris)
|
||||
src/
|
||||
pg_orbit.c # PG_MODULE_MAGIC entry point
|
||||
types.h # All struct definitions + constants
|
||||
astro_math.h # Shared astronomical helper functions
|
||||
pg_orbit.c # PG_MODULE_MAGIC + _PG_init() (GUC registration)
|
||||
types.h # All struct definitions + constants + DE body ID mapping
|
||||
astro_math.h # Shared astronomical helpers + observe_from_geocentric()
|
||||
# --- Satellite (v0.1.0) ---
|
||||
tle_type.c # TLE custom type (I/O, binary, 15 accessors)
|
||||
eci_type.c # ECI position type + geodetic/topocentric types
|
||||
@ -62,17 +64,21 @@ src/
|
||||
radio_funcs.c # io_phase_angle(), jupiter_cml(), burst_probability()
|
||||
lambert.c / lambert.h # Lambert transfer solver (Izzo 2015)
|
||||
transfer_funcs.c # lambert_transfer(), lambert_c3()
|
||||
# --- JPL DE Ephemeris (v0.3.0) ---
|
||||
de_reader.h / de_reader.c # Clean-room JPL DE binary reader (Chebyshev/Clenshaw)
|
||||
eph_provider.h / eph_provider.c # Provider dispatch, GUC, lazy init, frame rotation
|
||||
de_funcs.c # All _de() SQL function implementations
|
||||
lib/
|
||||
sat_code/ # Bill Gray's SGP4/SDP4 (MIT, git submodule)
|
||||
test/
|
||||
sql/ # 11 regression test suites
|
||||
sql/ # 12 regression test suites
|
||||
expected/ # Expected output
|
||||
docs/
|
||||
DESIGN.md # Architecture decisions, theory-to-code mappings
|
||||
Dockerfile # Starlight docs site (Astro + Caddy)
|
||||
package.json # Docs site dependencies
|
||||
astro.config.mjs # Starlight configuration
|
||||
src/content/docs/ # 34 MDX documentation pages
|
||||
src/content/docs/ # MDX documentation pages
|
||||
```
|
||||
|
||||
## Type System
|
||||
@ -89,7 +95,7 @@ All types are fixed-size, `STORAGE = plain`, `ALIGNMENT = double`. No TOAST over
|
||||
| `pass_event` | 48 | AOS/MAX/LOS times + max_el + AOS/LOS azimuth |
|
||||
| `heliocentric` | 24 | x, y, z in AU (ecliptic J2000 frame) |
|
||||
|
||||
## Function Domains (57 total)
|
||||
## Function Domains (68 total)
|
||||
|
||||
| Domain | Theory | Key Functions | Count |
|
||||
|--------|--------|---------------|-------|
|
||||
@ -101,9 +107,11 @@ All types are fixed-size, `STORAGE = plain`, `ALIGNMENT = double`. No TOAST over
|
||||
| Comets/asteroids | Two-body Keplerian | `kepler_propagate()`, `comet_observe()` | 2 |
|
||||
| Jupiter radio | Carr et al. (1983) | `jupiter_burst_probability()` | 3 |
|
||||
| Transfers | Lambert (Izzo 2015) | `lambert_transfer()`, `lambert_c3()` | 2 |
|
||||
| DE ephemeris | JPL DE440/441 (optional) | `planet_observe_de()`, `moon_observe_de()` | 11 |
|
||||
| GiST index | Altitude-band approximation | `&&` (overlap), `<->` (distance) | 8 |
|
||||
| Diagnostics | -- | `pg_orbit_ephemeris_info()` | 1 |
|
||||
|
||||
All 57 functions are `PARALLEL SAFE`. No global mutable state.
|
||||
All functions are `PARALLEL SAFE`. VSOP87/ELP82B functions are `IMMUTABLE` (compiled-in coefficients). DE functions are `STABLE` (external file dependency).
|
||||
|
||||
## Body IDs
|
||||
|
||||
@ -132,6 +140,9 @@ All 57 functions are `PARALLEL SAFE`. No global mutable state.
|
||||
3. **TEME frame**: Only 4 of 106 IAU-80 nutation terms (matching SGP4's internal model)
|
||||
4. **Solar system pipeline**: IAU 1976 precession, J2000 obliquity, GMST from Vallado Eq. 3-47
|
||||
5. **Never mix**: WGS-72 propagation + WGS-84 output. No other combination.
|
||||
6. **DE frame rotation**: DE positions (ICRS equatorial) pass through `equatorial_to_ecliptic()` at the provider boundary before entering the observation pipeline
|
||||
7. **Same-provider rule**: Both target and Earth must come from the same provider in any geocentric computation (never mix DE target with VSOP87 Earth)
|
||||
8. **DE AU consistency**: Verify DE header AU matches compiled-in `AU_KM` (149597870.7) at init time
|
||||
|
||||
### WGS-72 Constants (from Hoots & Roehrich STR#3, propagation only)
|
||||
```c
|
||||
@ -155,6 +166,50 @@ All 57 functions are `PARALLEL SAFE`. No global mutable state.
|
||||
#define J2000_JD 2451545.0 /* 2000 Jan 1.5 TT */
|
||||
```
|
||||
|
||||
## JPL DE Ephemeris (Optional)
|
||||
|
||||
v0.3.0 adds optional JPL DE440/441 ephemeris support (~0.1 milliarcsecond accuracy) alongside the existing VSOP87 pipeline (~1 arcsecond). DE functions are separate `_de()` variants — existing VSOP87 functions are completely unchanged.
|
||||
|
||||
### Architecture
|
||||
|
||||
- **Clean-room DE reader** (`de_reader.c`): ~250 lines of C. Parses the JPL binary format, evaluates Chebyshev polynomials via Clenshaw recurrence. No GPL dependency (avoids Bill Gray's `jpl_eph`).
|
||||
- **Per-backend lazy init**: Each PostgreSQL backend opens its own file descriptor on first `_de()` call. Never opens in `_PG_init()` (postmaster context). Safe for `PARALLEL SAFE`.
|
||||
- **VSOP87 fallback**: Every `_de()` function falls back to its VSOP87/ELP82B equivalent when DE is unavailable.
|
||||
- **STABLE volatility**: DE functions are `STABLE` (not `IMMUTABLE`) because output depends on an external file. Existing VSOP87 functions remain `IMMUTABLE`.
|
||||
|
||||
### GUC Configuration
|
||||
|
||||
```sql
|
||||
-- Set the path to a JPL DE binary file (requires superuser)
|
||||
ALTER SYSTEM SET pg_orbit.ephemeris_path = '/var/lib/postgres/de441.bin';
|
||||
SELECT pg_reload_conf();
|
||||
|
||||
-- Check which provider is active
|
||||
SELECT * FROM pg_orbit_ephemeris_info();
|
||||
```
|
||||
|
||||
| GUC | Type | Default | Context |
|
||||
|-----|------|---------|---------|
|
||||
| `pg_orbit.ephemeris_path` | string | `''` (empty = VSOP87 only) | `SIGHUP` (superuser only) |
|
||||
|
||||
### DE Function Variants
|
||||
|
||||
Every `_de()` function mirrors an existing VSOP87 function:
|
||||
|
||||
| DE Function | VSOP87 Equivalent | Volatility |
|
||||
|-------------|-------------------|------------|
|
||||
| `planet_heliocentric_de()` | `planet_heliocentric()` | STABLE |
|
||||
| `planet_observe_de()` | `planet_observe()` | STABLE |
|
||||
| `sun_observe_de()` | `sun_observe()` | STABLE |
|
||||
| `moon_observe_de()` | `moon_observe()` | STABLE |
|
||||
| `lambert_transfer_de()` | `lambert_transfer()` | STABLE |
|
||||
| `lambert_c3_de()` | `lambert_c3()` | STABLE |
|
||||
| `galilean_observe_de()` | `galilean_observe()` | STABLE |
|
||||
| `saturn_moon_observe_de()` | `saturn_moon_observe()` | STABLE |
|
||||
| `uranus_moon_observe_de()` | `uranus_moon_observe()` | STABLE |
|
||||
| `mars_moon_observe_de()` | `mars_moon_observe()` | STABLE |
|
||||
| `pg_orbit_ephemeris_info()` | — | STABLE |
|
||||
|
||||
## sat_code Submodule
|
||||
|
||||
Bill Gray's SGP4/SDP4: https://github.com/Bill-Gray/sat_code (MIT license)
|
||||
@ -169,7 +224,7 @@ Key files: `sgp4.cpp`, `sdp4.cpp`, `deep.cpp`, `common.cpp`, `basics.cpp`, `nora
|
||||
|
||||
## Testing
|
||||
|
||||
11 regression test suites via `make installcheck`:
|
||||
12 regression test suites via `make installcheck`:
|
||||
|
||||
| Suite | What it tests |
|
||||
|-------|--------------|
|
||||
@ -184,6 +239,7 @@ Key files: `sgp4.cpp`, `sdp4.cpp`, `deep.cpp`, `common.cpp`, `basics.cpp`, `nora
|
||||
| planet_observe | VSOP87 planets, sun_observe, moon_observe (ELP2000-82B) |
|
||||
| moon_observe | Galilean/Saturn/Uranus/Mars moons, Io phase, Jupiter CML, burst probability |
|
||||
| lambert_transfer | Lambert solver, lambert_c3, pork chop grid, error handling |
|
||||
| de_ephemeris | DE function fallback to VSOP87, cross-provider consistency, error handling |
|
||||
|
||||
## Error Handling Patterns
|
||||
|
||||
@ -196,9 +252,9 @@ Key files: `sgp4.cpp`, `sdp4.cpp`, `deep.cpp`, `common.cpp`, `basics.cpp`, `nora
|
||||
|
||||
**Live:** https://pg-orbit.warehack.ing
|
||||
|
||||
Starlight docs at `docs/` — 35 MDX pages covering all domains.
|
||||
Starlight docs at `docs/` — 36 MDX pages covering all domains.
|
||||
|
||||
Sections: Getting Started, Guides (8 domain walkthroughs), Workflow Translation (Skyfield/Horizons/GMAT/Radio Jupiter Pro comparisons), Reference (all 57 functions), 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 68 functions incl. DE variants), Architecture (Hamilton's principles, constant custody, observation pipeline), Performance (benchmarks).
|
||||
|
||||
### Local Development
|
||||
```bash
|
||||
@ -249,9 +305,10 @@ cd docs && make prod
|
||||
- `ereport(ERROR, ...)` for user-facing errors, never `elog(ERROR, ...)`
|
||||
- All memory via `palloc`/`pfree` (PostgreSQL memory contexts)
|
||||
- Comments explain "why", not "what"
|
||||
- No global mutable state — all computation from function arguments
|
||||
- No global mutable state — all computation from function arguments (exception: per-backend DE handle, lazily initialized, cleaned up via `on_proc_exit`)
|
||||
- Every function handling SGP4 must check the error return code
|
||||
- All functions marked `PARALLEL SAFE`
|
||||
- DE functions: always fall back to VSOP87/ELP82B on any error
|
||||
|
||||
## Git Conventions
|
||||
- One commit per logical change
|
||||
|
||||
14
Dockerfile
14
Dockerfile
@ -35,7 +35,7 @@ RUN make PG_CONFIG=${PG_CONFIG}
|
||||
# Install to system location (needed for installcheck)
|
||||
RUN make PG_CONFIG=${PG_CONFIG} install
|
||||
|
||||
# Run all 6 regression test suites against a throwaway cluster
|
||||
# Run all 12 regression test suites against a throwaway cluster
|
||||
RUN su postgres -c "/usr/lib/postgresql/${PG_MAJOR}/bin/initdb -D /tmp/pgtest" && \
|
||||
su postgres -c "/usr/lib/postgresql/${PG_MAJOR}/bin/pg_ctl -D /tmp/pgtest -l /tmp/pgtest.log start" && \
|
||||
su postgres -c "/usr/lib/postgresql/${PG_MAJOR}/bin/createuser -s root" && \
|
||||
@ -58,6 +58,18 @@ COPY docker/020_install_pg_orbit.sh /docker-entrypoint-initdb.d/
|
||||
# ── Stage 3: Standalone dev/test image ───────────────────────
|
||||
# Ready-to-run PostgreSQL with pg_orbit pre-installed.
|
||||
# For development, CI, and standalone experiments.
|
||||
#
|
||||
# Optional DE ephemeris at runtime (recommended):
|
||||
# docker run -v /path/to/de440.bin:/var/lib/postgresql/pg_orbit/de440.bin pg_orbit
|
||||
# Then: ALTER SYSTEM SET pg_orbit.ephemeris_path = '/var/lib/postgresql/pg_orbit/de440.bin';
|
||||
#
|
||||
# Or bake into the image (115 MB for DE440, 3.1 GB for DE441):
|
||||
# Place the DE file in the build context, then:
|
||||
# docker build --build-arg DE_FILE=de440.bin -t pg_orbit:de440 .
|
||||
FROM postgres:${PG_MAJOR}-bookworm AS standalone
|
||||
|
||||
COPY --from=artifact / /
|
||||
|
||||
# Create the pg_orbit data directory for DE ephemeris files
|
||||
RUN mkdir -p /var/lib/postgresql/pg_orbit && \
|
||||
chown postgres:postgres /var/lib/postgresql/pg_orbit
|
||||
|
||||
9
Makefile
9
Makefile
@ -1,6 +1,7 @@
|
||||
MODULE_big = pg_orbit
|
||||
EXTENSION = pg_orbit
|
||||
DATA = sql/pg_orbit--0.1.0.sql sql/pg_orbit--0.2.0.sql sql/pg_orbit--0.1.0--0.2.0.sql
|
||||
DATA = sql/pg_orbit--0.1.0.sql sql/pg_orbit--0.2.0.sql sql/pg_orbit--0.1.0--0.2.0.sql \
|
||||
sql/pg_orbit--0.3.0.sql sql/pg_orbit--0.2.0--0.3.0.sql
|
||||
|
||||
# Our extension C sources
|
||||
OBJS = src/pg_orbit.o src/tle_type.o src/eci_type.o src/observer_type.o \
|
||||
@ -10,7 +11,8 @@ OBJS = src/pg_orbit.o src/tle_type.o src/eci_type.o src/observer_type.o \
|
||||
src/precession.o src/sidereal_time.o src/planet_funcs.o \
|
||||
src/tass17.o src/gust86.o src/marssat.o src/l12.o \
|
||||
src/moon_funcs.o src/radio_funcs.o \
|
||||
src/lambert.o src/transfer_funcs.o
|
||||
src/lambert.o src/transfer_funcs.o \
|
||||
src/de_reader.o src/eph_provider.o src/de_funcs.o
|
||||
|
||||
# sat_code C++ sources (compiled with g++, linked with extern "C" symbols)
|
||||
SAT_CODE_DIR = lib/sat_code
|
||||
@ -24,7 +26,8 @@ OBJS += $(SAT_CODE_OBJS)
|
||||
|
||||
# Regression tests
|
||||
REGRESS = tle_parse sgp4_propagate coord_transforms pass_prediction gist_index convenience \
|
||||
star_observe kepler_comet planet_observe moon_observe lambert_transfer
|
||||
star_observe kepler_comet planet_observe moon_observe lambert_transfer \
|
||||
de_ephemeris
|
||||
REGRESS_OPTS = --inputdir=test
|
||||
|
||||
# Need C++ runtime for sat_code
|
||||
|
||||
@ -7,8 +7,8 @@ planet positions, observe 19 planetary moons, predict Jupiter radio bursts, and
|
||||
plan interplanetary trajectories — all from standard SQL. Think PostGIS, but for
|
||||
objects in space.
|
||||
|
||||
57 functions. 7 custom types. All `PARALLEL SAFE`. Zero external dependencies at
|
||||
runtime.
|
||||
68 functions. 7 custom types. All `PARALLEL SAFE`. Optional JPL DE440/441 support
|
||||
for sub-arcsecond accuracy; core functions have zero external dependencies at runtime.
|
||||
|
||||
## Installation
|
||||
|
||||
@ -181,7 +181,7 @@ Measured on PostgreSQL 17, single backend:
|
||||
|
||||
## Testing
|
||||
|
||||
11 regression test suites covering all domains:
|
||||
12 regression test suites covering all domains:
|
||||
|
||||
```bash
|
||||
make installcheck PG_CONFIG=/usr/bin/pg_config
|
||||
@ -189,7 +189,7 @@ make installcheck PG_CONFIG=/usr/bin/pg_config
|
||||
|
||||
Tests: TLE parsing, SGP4/SDP4 propagation, coordinate transforms, pass prediction,
|
||||
GiST indexing, convenience functions, star observation, Keplerian propagation,
|
||||
planet observation, moon observation, and Lambert transfers.
|
||||
planet observation, moon observation, Lambert transfers, and DE ephemeris.
|
||||
|
||||
## Documentation
|
||||
|
||||
|
||||
168
docs/DESIGN.md
168
docs/DESIGN.md
@ -484,11 +484,16 @@ no `new`, no static buffers.
|
||||
|
||||
### No Global Mutable State
|
||||
|
||||
There are no file-scope variables, no static locals that accumulate
|
||||
state, no caches. Every function computes from its arguments alone.
|
||||
This is required for `PARALLEL SAFE` (all pg_orbit functions are
|
||||
declared PARALLEL SAFE) and avoids cross-session contamination in
|
||||
a multi-backend PostgreSQL server.
|
||||
For v0.1.0/v0.2.0 functions, there are no file-scope variables, no
|
||||
static locals that accumulate state, no caches. Every function computes
|
||||
from its arguments alone.
|
||||
|
||||
The v0.3.0 DE ephemeris layer introduces per-backend static state in
|
||||
`eph_provider.c` (file descriptor, coefficient cache, init flags). This
|
||||
is safe because each backend gets its own copy after fork(). The handle
|
||||
is cleaned up via `on_proc_exit()`. All pg_orbit functions remain
|
||||
`PARALLEL SAFE` -- parallel workers each open their own DE handle
|
||||
independently.
|
||||
|
||||
sat_code itself has no global mutable state. The propagator state is
|
||||
entirely in the `params[]` array and the `tle_t` struct, both of which
|
||||
@ -588,3 +593,156 @@ implementation in pg_orbit and sat_code.
|
||||
| Observer to ECEF | Geodesy standard | WGS-84 ellipsoid surface point to Cartesian | `src/coord_funcs.c:observer_to_ecef()` lines 143-156 |
|
||||
| Range rate | Dot product | Projection of relative velocity onto line-of-sight unit vector | `src/coord_funcs.c:eci_to_topocentric()` line 618 |
|
||||
| Near/deep selection | Hoots & Roehrich STR#3 | Period threshold: 225 minutes (n < 2*pi/225 rad/min) | `sat_code/norad.h:select_ephemeris()` |
|
||||
|
||||
|
||||
## 10. JPL DE Ephemeris Architecture (v0.3.0)
|
||||
|
||||
v0.3.0 adds optional JPL DE440/441 ephemeris support (~0.1 milliarcsecond
|
||||
accuracy) without modifying any existing VSOP87/ELP82B code path. This
|
||||
section documents the architectural decisions specific to DE integration.
|
||||
|
||||
### The Fundamental Tension
|
||||
|
||||
pg_orbit's core properties (compiled-in coefficients, no file I/O, no
|
||||
mutable state) are precisely what DE441 challenges. A ~3GB binary file
|
||||
introduces file dependency, per-backend state (file descriptor,
|
||||
coefficient cache), and OS-level file descriptor management across
|
||||
PostgreSQL's multi-process model.
|
||||
|
||||
The architecture resolves this by treating DE as an **additive layer**:
|
||||
new `_de()` function variants alongside existing functions. Existing
|
||||
functions are untouched. VSOP87 is the bulkhead -- always compiled in,
|
||||
always works, always IMMUTABLE.
|
||||
|
||||
### Why Separate Functions (Not a GUC Switch)
|
||||
|
||||
Considered and rejected: a single `planet_observe()` that dispatches to
|
||||
DE or VSOP87 based on a GUC setting.
|
||||
|
||||
The problem is volatility. VSOP87 functions are IMMUTABLE -- their data
|
||||
is compiled into the binary. PostgreSQL can:
|
||||
- Cache results in expression indexes (almanac tables, pork chop grids)
|
||||
- Constant-fold during planning (bake results into prepared statements)
|
||||
|
||||
If `planet_observe()` dispatched to DE based on a GUC, it would need to
|
||||
be STABLE, losing these optimizations for ALL users, even those without
|
||||
a DE file. Separate `_de()` variants let VSOP87 users keep IMMUTABLE
|
||||
and DE users get STABLE -- honest volatility for both.
|
||||
|
||||
### Clean-Room DE Reader Design
|
||||
|
||||
The JPL DE binary format (used by DE405, DE430, DE440, DE441) consists
|
||||
of fixed-size records of `double` values:
|
||||
|
||||
- **Record 1 (header):** Start/end JD, record interval, number of
|
||||
coefficients, AU value, Earth-Moon mass ratio (EMRAT), coefficient
|
||||
layout table (3 values per body: offset, ncoeff, nsub)
|
||||
- **Record 2:** Constant values
|
||||
- **Records 3+:** Chebyshev polynomial coefficients covering a time
|
||||
interval for all 13 body groups
|
||||
|
||||
The reader is implemented in ~250 lines of C (`de_reader.c`), using:
|
||||
- Raw POSIX I/O (`open`/`lseek`/`read`, not `fread`) to avoid libc
|
||||
buffering issues across fork
|
||||
- Clenshaw recurrence for Chebyshev evaluation (~15 lines)
|
||||
- Single-record coefficient cache (most queries hit consecutive times)
|
||||
|
||||
#### Why Not jpl_eph
|
||||
|
||||
Bill Gray's `jpl_eph` (GPL-2+) would be the obvious choice, but:
|
||||
1. GPL-2+ license constrains pg_orbit's licensing flexibility
|
||||
2. Uses global statics (`static int init_err_code`)
|
||||
3. Written in C++ (`jpleph.cpp`); pg_orbit is pure C
|
||||
4. We only need position queries, not velocity or nutation
|
||||
|
||||
The format is well-documented and the algorithm is straightforward.
|
||||
A clean-room implementation in ~250 lines avoids all four issues.
|
||||
|
||||
### Per-Backend Lazy Initialization
|
||||
|
||||
PostgreSQL parallel workers are forked from the **postmaster, not the
|
||||
leader**. They don't inherit the leader's file descriptors or initialized
|
||||
handles.
|
||||
|
||||
Each backend (including parallel workers) gets its own independently-
|
||||
opened handle via lazy init on first `_de()` call:
|
||||
|
||||
```
|
||||
Backend A calls planet_observe_de() -> opens own fd, own cache
|
||||
Worker W1 calls planet_observe_de() -> opens own fd, own cache
|
||||
Worker W2 calls planet_observe_de() -> opens own fd, own cache
|
||||
```
|
||||
|
||||
Static per-backend variables (`de_handle_ptr`, `de_init_attempted`) are
|
||||
safe because after `fork()`, each process gets its own copy-on-write
|
||||
address space.
|
||||
|
||||
**Critical rule:** Never open the DE file in `_PG_init()`. That runs in
|
||||
the postmaster, and the file descriptor would be inherited by all children
|
||||
with undefined stream behavior. Always lazy init.
|
||||
|
||||
Cleanup is via `on_proc_exit(eph_cleanup, 0)`, registered in `_PG_init()`.
|
||||
|
||||
### ICRS-to-Ecliptic Frame Rotation
|
||||
|
||||
DE ephemerides return positions in the ICRS equatorial frame. The
|
||||
pg_orbit observation pipeline expects ecliptic J2000. The conversion
|
||||
happens at the provider boundary in `eph_provider.c`:
|
||||
|
||||
```
|
||||
DE position (ICRS equatorial) -> equatorial_to_ecliptic() -> ecliptic J2000
|
||||
^
|
||||
already in astro_math.h since v0.2.0
|
||||
```
|
||||
|
||||
This rotation is applied to BOTH the target body and Earth positions
|
||||
before they enter `observe_from_geocentric()`. The frame conversion
|
||||
happens exactly once, at the provider boundary.
|
||||
|
||||
### Earth Position from DE
|
||||
|
||||
DE ephemerides don't store Earth directly. They store the Earth-Moon
|
||||
Barycenter (EMB, body 3) and the Moon (body 10). Earth's position is
|
||||
derived:
|
||||
|
||||
```
|
||||
Earth = EMB - Moon / (1 + EMRAT)
|
||||
```
|
||||
|
||||
where EMRAT (~81.300587) is the Earth-Moon mass ratio from the DE
|
||||
header. This calculation happens in `eph_provider.c:de_get_earth_helio_ecliptic()`.
|
||||
|
||||
### Constant Chain of Custody Extension
|
||||
|
||||
Three new rules for the DE pipeline:
|
||||
|
||||
| Rule | What | Why |
|
||||
|------|------|-----|
|
||||
| 6 | DE positions pass through `equatorial_to_ecliptic()` at the provider boundary | DE returns ICRS equatorial; the observation pipeline expects ecliptic J2000 |
|
||||
| 7 | Both target and Earth must come from the same provider | Mixing DE Mars with VSOP87 Earth introduces frame inconsistency |
|
||||
| 8 | DE header AU must match compiled-in `AU_KM` (149597870.7) | AU defines the length scale; a mismatch corrupts all distance calculations |
|
||||
|
||||
### Fallback Strategy
|
||||
|
||||
Every `_de()` function follows this pattern:
|
||||
|
||||
1. Try to get position from DE
|
||||
2. If DE succeeds, use DE result
|
||||
3. If DE fails (no file, JD out of range, I/O error):
|
||||
- If DE was explicitly configured (GUC set), emit `ereport(NOTICE)`
|
||||
- Fall back to VSOP87/ELP82B equivalent
|
||||
- Return the VSOP87 result
|
||||
|
||||
When `pg_orbit.ephemeris_path` is empty (default), DE functions fall
|
||||
back silently -- no NOTICE, no overhead, identical results to the
|
||||
non-DE variants.
|
||||
|
||||
### Lauren Bugs (DE-Specific)
|
||||
|
||||
| Bug | What "Can't Happen" | How It Happens | Mitigation |
|
||||
|-----|---------------------|----------------|------------|
|
||||
| L1 | File won't change mid-query | DBA replaces file; new parallel worker opens new version | Per-backend handle is stable for backend lifetime |
|
||||
| L2 | File will always be there | Docker restart with ephemeral volume | VSOP87 fallback always available |
|
||||
| L3 | File is always valid | Partial download, bit rot, wrong endianness | Canary validation: Earth at J2000.0 ~1 AU from Sun; byte-order detection via AU constant |
|
||||
| L4 | All backends use same version | Primary has DE441, replica has DE440 | AU consistency check; log ephemeris version at NOTICE |
|
||||
| L5 | read() always succeeds | NFS timeout, disk error, fd limit | Every `read()` return checked; propagate as `ereport(ERROR)` |
|
||||
|
||||
@ -59,6 +59,7 @@ export default defineConfig({
|
||||
{ label: "Jupiter Radio Burst Prediction", slug: "guides/jupiter-radio-bursts" },
|
||||
{ label: "Interplanetary Trajectories", slug: "guides/interplanetary-trajectories" },
|
||||
{ label: "Conjunction Screening", slug: "guides/conjunction-screening" },
|
||||
{ label: "JPL DE Ephemeris", slug: "guides/de-ephemeris" },
|
||||
],
|
||||
},
|
||||
{
|
||||
@ -81,6 +82,7 @@ export default defineConfig({
|
||||
{ label: "Functions: Stars & Comets", slug: "reference/functions-stars-comets" },
|
||||
{ label: "Functions: Radio", slug: "reference/functions-radio" },
|
||||
{ label: "Functions: Transfers", slug: "reference/functions-transfers" },
|
||||
{ label: "Functions: DE Ephemeris", slug: "reference/functions-de" },
|
||||
{ label: "Operators & GiST Index", slug: "reference/operators-gist" },
|
||||
{ label: "Body ID Reference", slug: "reference/body-ids" },
|
||||
{ label: "Constants & Accuracy", slug: "reference/constants-accuracy" },
|
||||
|
||||
@ -124,6 +124,39 @@ pg_orbit deliberately does **not** use a higher-precision GMST model (e.g., IAU
|
||||
|
||||
This is the constant chain of custody in action: match the precision of the input, not the precision available in the literature.
|
||||
|
||||
## Rules 6-8: DE Ephemeris Pipeline (v0.3.0)
|
||||
|
||||
The v0.3.0 DE integration adds three rules to the chain of custody.
|
||||
|
||||
### Rule 6: ICRS-to-ecliptic frame rotation at the provider boundary
|
||||
|
||||
DE ephemerides return positions in the ICRS equatorial frame. The observation pipeline expects ecliptic J2000. The conversion uses `equatorial_to_ecliptic()` from `astro_math.h`, applied in `eph_provider.c` before the position enters the shared observation pipeline. This rotation is applied to both the target body and Earth.
|
||||
|
||||
### Rule 7: Same-provider consistency
|
||||
|
||||
Both the target body and Earth must come from the same provider in any geocentric computation. If DE succeeds for the target but fails for Earth (or vice versa), the entire computation falls back to VSOP87. This prevents mixing DE Mars with VSOP87 Earth, which would introduce a frame-dependent systematic error at the arcsecond level.
|
||||
|
||||
```c
|
||||
/* From de_funcs.c: both positions or neither */
|
||||
if (eph_de_planet(body_id, jd, target_xyz) &&
|
||||
eph_de_earth(jd, earth_xyz))
|
||||
{
|
||||
/* Use DE positions for both */
|
||||
}
|
||||
else
|
||||
{
|
||||
/* Fall back to VSOP87 for both */
|
||||
}
|
||||
```
|
||||
|
||||
### Rule 8: AU consistency verification
|
||||
|
||||
The DE header contains an AU value (in km). At init time, `eph_provider.c` verifies this matches pg_orbit's compiled-in `AU_KM` constant (149597870.7 km, IAU 2012). A mismatch would corrupt every distance calculation. If they disagree, the DE file is rejected and fallback to VSOP87 activates with a log message.
|
||||
|
||||
<Aside type="note" title="For maintainers">
|
||||
If you are modifying `eph_provider.c` or `de_funcs.c`, remember that Rule 7 is the critical invariant. Never return a DE position for one body and a VSOP87 position for another within the same geocentric computation. The conditional must gate both positions atomically.
|
||||
</Aside>
|
||||
|
||||
## Verification
|
||||
|
||||
The chain of custody is verified through the Vallado 518 test vectors --- 518 reference propagations across a range of orbit types (LEO, MEO, GEO, deep-space, high-eccentricity). Every test vector must match to machine epsilon before any other development proceeds.
|
||||
|
||||
@ -100,7 +100,7 @@ Hamilton defined ultra-reliable software as software that behaves correctly unde
|
||||
|
||||
### Zero global mutable state
|
||||
|
||||
There are no file-scope variables, no static locals, no caches. Every function computes from its arguments alone. This is not a style preference --- it is required for PostgreSQL's `PARALLEL SAFE` declaration. All 57 pg_orbit functions carry this declaration, meaning the query planner can distribute work across multiple CPU cores without coordination.
|
||||
For v0.1.0/v0.2.0 functions, there are no file-scope variables, no static locals, no caches. Every function computes from its arguments alone. The v0.3.0 DE ephemeris layer introduces per-backend static state (a file descriptor and coefficient cache in `eph_provider.c`), but each backend gets its own copy after `fork()` --- no shared state between processes. All 68 pg_orbit functions carry the `PARALLEL SAFE` declaration, meaning the query planner can distribute work across multiple CPU cores without coordination.
|
||||
|
||||
### Fixed-size types
|
||||
|
||||
@ -179,3 +179,9 @@ A negative time of flight. The Lambert solver might converge on a mathematically
|
||||
When computing the topocentric observation of Earth (body ID 3), the geocentric vector is zero --- the observer is on the body being observed. Division by zero in the range computation. pg_orbit catches this case and returns a clear error rather than NaN or infinity propagating through the rest of the pipeline.
|
||||
|
||||
These are not edge cases in the traditional sense. They are the inputs that a SQL user will inevitably produce when exploring the system with ad hoc queries, and they must produce clear errors rather than silently wrong results.
|
||||
|
||||
### DE ephemeris file dependency (v0.3.0)
|
||||
|
||||
The v0.3.0 DE integration introduces five new failure domains --- external file dependency, per-backend mutable state, OS-level file descriptor management, frame boundary crossings, and volatility contract management --- into a system that previously had none.
|
||||
|
||||
Each is handled by the principle of graceful degradation: DE functions fall back to VSOP87 on any DE-specific failure. The existing VSOP87 pipeline is the bulkhead --- always compiled in, always works, always IMMUTABLE. See the [DE Ephemeris guide](/guides/de-ephemeris/) for details on the fallback strategy.
|
||||
|
||||
@ -87,15 +87,17 @@ pg_tle *result = (pg_tle *) palloc(sizeof(pg_tle));
|
||||
PG_RETURN_POINTER(result);
|
||||
```
|
||||
|
||||
## Zero global mutable state
|
||||
## Minimal global mutable state
|
||||
|
||||
There are no file-scope variables, no static locals that accumulate state, no caches. Every function computes from its arguments alone.
|
||||
The v0.1.0/v0.2.0 functions have zero global mutable state --- no file-scope variables, no static locals, no caches. Every function computes from its arguments alone.
|
||||
|
||||
The v0.3.0 DE ephemeris layer introduces a controlled exception: per-backend static variables in `eph_provider.c` (`de_handle_ptr`, `de_init_attempted`, `de_init_success`). These are safe because each PostgreSQL backend (including parallel workers) gets its own copy after `fork()` --- no shared state between processes. The handle is opened lazily on first `_de()` call and cleaned up via `on_proc_exit()`. See the [DE Ephemeris guide](/guides/de-ephemeris/) for details.
|
||||
|
||||
This guarantee has three consequences:
|
||||
|
||||
### PARALLEL SAFE
|
||||
|
||||
All 57 pg_orbit functions are declared `PARALLEL SAFE` in the SQL definition. This tells PostgreSQL's query planner that the function can be executed in parallel worker processes without coordination. For bulk operations like propagating 12,000 TLEs, the planner can distribute work across multiple CPU cores:
|
||||
All 68 pg_orbit functions are declared `PARALLEL SAFE` in the SQL definition. This tells PostgreSQL's query planner that the function can be executed in parallel worker processes without coordination. For bulk operations like propagating 12,000 TLEs, the planner can distribute work across multiple CPU cores:
|
||||
|
||||
```sql
|
||||
-- PostgreSQL may parallelize this across available cores
|
||||
@ -183,4 +185,6 @@ When `ereport(ERROR)` fires inside a pg_orbit function, PostgreSQL's error recov
|
||||
3. Rolls back the current transaction
|
||||
4. Returns an error message to the client
|
||||
|
||||
Because pg_orbit uses only `palloc` and has no global state, there is nothing to clean up beyond what PostgreSQL's context system handles automatically. No file handles, no sockets, no mutex locks, no C++ destructors. The extension is always in a consistent state after error recovery.
|
||||
Because pg_orbit uses only `palloc` and has no global state for v0.1.0/v0.2.0 functions, there is nothing to clean up beyond what PostgreSQL's context system handles automatically. No sockets, no mutex locks, no C++ destructors.
|
||||
|
||||
The v0.3.0 DE reader holds a file descriptor in per-backend static state. This is cleaned up via `on_proc_exit(eph_cleanup, 0)`, registered during `_PG_init()`. If `ereport(ERROR)` fires during a DE function, the file descriptor persists (it will be reused by the next DE call in the same backend) --- it is not leaked, just kept open for the backend's lifetime. The extension is always in a consistent state after error recovery.
|
||||
|
||||
@ -123,7 +123,7 @@ flowchart TD
|
||||
|
||||
The result is packed into a `pg_topocentric` struct with range computed from the geocentric distance ($\text{range}_\text{km} = d_\text{AU} \times 149597870.7$). Range rate is set to 0.0 for planetary observations (velocity computation is not yet implemented for the VSOP87 pipeline).
|
||||
|
||||
**Code**: `src/planet_funcs.c:observe_from_geocentric()`
|
||||
**Code**: `src/astro_math.h:observe_from_geocentric()` (shared by `planet_funcs.c`, `moon_funcs.c`, and `de_funcs.c`)
|
||||
</Steps>
|
||||
|
||||
## Pipeline variants
|
||||
@ -144,6 +144,16 @@ The seven-stage pipeline applies to planets observed via VSOP87. Other observati
|
||||
3. Transform to heliocentric ecliptic J2000 (parent offset + frame rotation)
|
||||
4. Proceed from stage 3 of the standard pipeline (geocentric ecliptic)
|
||||
</TabItem>
|
||||
<TabItem label="DE ephemeris (v0.3.0)">
|
||||
The `_de()` function variants replace stages 1-2 with JPL DE positions:
|
||||
|
||||
1. DE reader returns target position in ICRS equatorial frame
|
||||
2. DE reader returns Earth position (derived from EMB - Moon/(1+EMRAT))
|
||||
3. Both positions converted from ICRS equatorial to ecliptic J2000 via `equatorial_to_ecliptic()` at the provider boundary
|
||||
4. Proceed from stage 3 of the standard pipeline (geocentric ecliptic)
|
||||
|
||||
Rule 7 of the constant chain of custody requires both target and Earth to come from the same provider. If either DE call fails, both fall back to VSOP87.
|
||||
</TabItem>
|
||||
<TabItem label="Satellites (SGP4)">
|
||||
Satellites use a fundamentally different pipeline. SGP4 outputs TEME (True Equator, Mean Equinox) positions, not heliocentric ecliptic. The satellite pipeline:
|
||||
|
||||
@ -171,7 +181,7 @@ Several simplifications are deliberate.
|
||||
The pipeline uses precession but not nutation. For 1-arcsecond VSOP87 positions, the ~9-arcsecond nutation correction is below the noise floor of the ephemeris. Adding nutation would increase computation cost without improving practical accuracy.
|
||||
</Aside>
|
||||
|
||||
**No aberration correction.** Annual aberration shifts apparent positions by up to 20 arcseconds, but for observation planning (which quadrant of the sky is Jupiter in tonight?) this is irrelevant. Sub-arcsecond work should use SPICE or Skyfield with DE441.
|
||||
**No aberration correction.** Annual aberration shifts apparent positions by up to 20 arcseconds, but for observation planning (which quadrant of the sky is Jupiter in tonight?) this is irrelevant. For sub-arcsecond planet positions, pg_orbit v0.3.0 supports [optional DE440/441 ephemeris files](/guides/de-ephemeris/); for apparent position corrections (aberration, light-time), use SPICE or Skyfield.
|
||||
|
||||
**No light-time iteration.** The positions returned are geometric, not apparent. Light-time corrections of a few minutes for the outer planets shift the apparent position by a fraction of an arcsecond at most --- again, below the VSOP87 accuracy floor.
|
||||
|
||||
@ -184,7 +194,7 @@ To add a new observation target, identify which stages change:
|
||||
| New target | Stages replaced | Stages reused |
|
||||
|-----------|-----------------|---------------|
|
||||
| New moon theory | 1-3 (new theory + parent VSOP87) | 4-7 |
|
||||
| New planet ephemeris | 1 (new theory replaces VSOP87) | 2-7 (Earth still VSOP87) |
|
||||
| DE ephemeris (v0.3.0) | 1-2 (DE positions for both target and Earth) | 3-7 |
|
||||
| Near-Earth asteroid | 1 (Keplerian propagation) | 2-7 |
|
||||
| Distant star | 1-3 (catalog lookup, no distance) | 5-7 (skip obliquity) |
|
||||
|
||||
|
||||
@ -44,7 +44,9 @@ This is the canonical SGP4/SDP4 reference. All subsequent implementations, inclu
|
||||
| Range rate | Dot product | Projection of relative velocity onto line-of-sight unit vector | `src/coord_funcs.c:eci_to_topocentric()` |
|
||||
| Semi-major axis from $n$ | Kepler's third law | $a = (k_e / n)^{2/3}$ in earth radii | `src/tle_type.c:tle_perigee()` |
|
||||
| IAU 1976 precession | Lieske et al. (1977) | Three Euler angles $\zeta_A$, $z_A$, $\theta_A$ for precession from J2000 to date | `src/precession.c:precess_j2000_to_date()` |
|
||||
| Ecliptic to equatorial | IAU | X-axis rotation by obliquity $\varepsilon_0 = 23.4392911\degree$ | `src/planet_funcs.c:ecliptic_to_equatorial()` |
|
||||
| Ecliptic to equatorial | IAU | X-axis rotation by obliquity $\varepsilon_0 = 23.4392911\degree$ | `src/astro_math.h:ecliptic_to_equatorial()` |
|
||||
| Equatorial to ecliptic | IAU | Inverse obliquity rotation (for DE ICRS→ecliptic) | `src/astro_math.h:equatorial_to_ecliptic()` |
|
||||
| Observation pipeline | Standard | Geocentric ecliptic → equatorial → precess → topocentric | `src/astro_math.h:observe_from_geocentric()` |
|
||||
|
||||
### Primary references
|
||||
|
||||
@ -58,6 +60,7 @@ This is the canonical SGP4/SDP4 reference. All subsequent implementations, inclu
|
||||
|--------|--------|--------|----------|---------------|
|
||||
| VSOP87 | Bretagnon & Francou (1988) | Mercury through Neptune | ~1 arcsecond | `src/vsop87.c` |
|
||||
| ELP2000-82B | Chapront-Touze & Chapront (1988) | Moon | ~10 arcseconds | `src/elp82b.c` |
|
||||
| JPL DE440/441 (optional) | JPL (2021) | All bodies + Moon | ~0.1 milliarcsecond | `src/de_reader.c`, `src/eph_provider.c` |
|
||||
|
||||
### VSOP87
|
||||
|
||||
@ -153,7 +156,11 @@ A quick reference for finding the implementation of a specific theory.
|
||||
| `src/coord_funcs.c` | Coordinate transforms | ~650 |
|
||||
| `src/pass_funcs.c` | Pass prediction algorithm | ~550 |
|
||||
| `src/gist_tle.c` | GiST altitude-band index | ~400 |
|
||||
| `src/planet_funcs.c` | Observation pipeline | ~250 |
|
||||
| `src/planet_funcs.c` | VSOP87 observation pipeline | ~200 |
|
||||
| `src/de_reader.c` | Clean-room JPL DE binary reader (Chebyshev/Clenshaw) | ~400 |
|
||||
| `src/eph_provider.c` | DE provider dispatch, GUC, lazy init, frame rotation | ~350 |
|
||||
| `src/de_funcs.c` | All `_de()` SQL function implementations | ~650 |
|
||||
| `src/astro_math.h` | Shared math: obliquity rotation, observation pipeline | ~220 |
|
||||
| `src/radio_funcs.c` | Jupiter radio emission | ~200 |
|
||||
| `sat_code/sgp4.cpp` | SGP4 near-earth propagator | ~300 |
|
||||
| `sat_code/sdp4.cpp` | SDP4 deep-space propagator | ~200 |
|
||||
|
||||
@ -97,20 +97,24 @@ import { Tabs, TabItem, Steps, Aside } from "@astrojs/starlight/components";
|
||||
|
||||
## Running the test suite
|
||||
|
||||
If building from source, the regression tests verify all 57 functions across 11 test suites:
|
||||
If building from source, the regression tests verify all 68 functions across 12 test suites:
|
||||
|
||||
```bash
|
||||
make installcheck PG_CONFIG=/usr/bin/pg_config
|
||||
```
|
||||
|
||||
This runs the tests listed in the `REGRESS` variable: TLE parsing, SGP4 propagation, coordinate transforms, pass prediction, GiST indexing, convenience functions, star observation, Keplerian propagation, planet observation, moon observation, and Lambert transfers.
|
||||
This runs the tests listed in the `REGRESS` variable: TLE parsing, SGP4 propagation, coordinate transforms, pass prediction, GiST indexing, convenience functions, star observation, Keplerian propagation, planet observation, moon observation, Lambert transfers, and DE ephemeris.
|
||||
|
||||
## Upgrading from v0.1.0
|
||||
## Upgrading
|
||||
|
||||
If you have pg_orbit 0.1.0 installed (satellite-only), upgrade to 0.2.0:
|
||||
If you have a previous version installed, upgrade in place:
|
||||
|
||||
```sql
|
||||
-- From v0.1.0 (satellite-only) to v0.2.0 (solar system)
|
||||
ALTER EXTENSION pg_orbit UPDATE TO '0.2.0';
|
||||
|
||||
-- From v0.2.0 to v0.3.0 (DE ephemeris support)
|
||||
ALTER EXTENSION pg_orbit UPDATE TO '0.3.0';
|
||||
```
|
||||
|
||||
This adds all solar system functions (planets, moons, stars, comets, radio, Lambert transfers) while preserving your existing TLE data and satellite functions.
|
||||
Each migration adds new functions while preserving existing data and functions.
|
||||
|
||||
167
docs/src/content/docs/guides/de-ephemeris.mdx
Normal file
167
docs/src/content/docs/guides/de-ephemeris.mdx
Normal file
@ -0,0 +1,167 @@
|
||||
---
|
||||
title: JPL DE Ephemeris
|
||||
sidebar:
|
||||
order: 9
|
||||
---
|
||||
|
||||
import { Steps, Aside, Tabs, TabItem } from "@astrojs/starlight/components";
|
||||
|
||||
pg_orbit v0.3.0 adds optional support for JPL Development Ephemeris files (DE440/DE441), bringing positional accuracy from VSOP87's ~1 arcsecond to DE441's ~0.1 milliarcsecond. The upgrade path is designed so that nothing changes unless you opt in.
|
||||
|
||||
## When you need DE
|
||||
|
||||
VSOP87 is sufficient for most purposes --- visual observation planning, rise/set times, conjunction identification, telescope pointing at optical wavelengths. You need DE when:
|
||||
|
||||
- **GHz dish pointing.** A 10m dish at 10 GHz has a ~0.1 degree beamwidth. VSOP87's 1 arcsecond error is well within the beam, but systematic campaigns benefit from the tighter DE error budget.
|
||||
- **Precision astrometry.** Comparing CCD field positions against computed planet positions at the sub-arcsecond level.
|
||||
- **Occultation timing.** Predicting when a planet or asteroid will occult a star requires milliarcsecond precision in the planet's geocentric position.
|
||||
- **Interplanetary mission design.** Lambert transfer C3 values are sensitive to planet position accuracy. DE441 eliminates VSOP87 as an error source.
|
||||
|
||||
<Aside type="tip" title="Zero-change default">
|
||||
If you don't configure a DE file, all `_de()` functions silently fall back to VSOP87/ELP2000-82B. You get identical results to the non-DE variants, with no performance penalty. The extension ships with VSOP87 compiled in --- DE is strictly opt-in.
|
||||
</Aside>
|
||||
|
||||
## Setup
|
||||
|
||||
<Steps>
|
||||
1. **Obtain a DE file.** Download DE441 (~3.1 GB, covers -13200 to +17191) or DE440 (~115 MB, covers 1550 to 2650) from JPL:
|
||||
|
||||
```bash
|
||||
# DE440 (smaller, covers 1550-2650, sufficient for most uses)
|
||||
curl -O https://ssd.jpl.nasa.gov/ftp/eph/planets/Linux/de440/linux_p1550p2650.440
|
||||
|
||||
# DE441 (full range, large file)
|
||||
curl -O https://ssd.jpl.nasa.gov/ftp/eph/planets/Linux/de441/linux_m13000p17000.441
|
||||
```
|
||||
|
||||
2. **Place the file** where PostgreSQL can read it:
|
||||
|
||||
```bash
|
||||
sudo mkdir -p /var/lib/postgres/pg_orbit
|
||||
sudo cp linux_p1550p2650.440 /var/lib/postgres/pg_orbit/de440.bin
|
||||
sudo chown postgres:postgres /var/lib/postgres/pg_orbit/de440.bin
|
||||
```
|
||||
|
||||
3. **Set the GUC** (requires superuser):
|
||||
|
||||
```sql
|
||||
ALTER SYSTEM SET pg_orbit.ephemeris_path = '/var/lib/postgres/pg_orbit/de440.bin';
|
||||
SELECT pg_reload_conf();
|
||||
```
|
||||
|
||||
4. **Verify:**
|
||||
|
||||
```sql
|
||||
SELECT * FROM pg_orbit_ephemeris_info();
|
||||
```
|
||||
|
||||
You should see `provider = 'JPL_DE'` with the file's JD range and version.
|
||||
</Steps>
|
||||
|
||||
## Using DE functions
|
||||
|
||||
Every VSOP87 function has a `_de()` counterpart with the same signature:
|
||||
|
||||
<Tabs>
|
||||
<TabItem label="VSOP87 (existing)">
|
||||
```sql
|
||||
-- Always VSOP87, always IMMUTABLE
|
||||
SELECT topo_azimuth(planet_observe(4, '40.0N 105.3W 1655m'::observer, now())) AS mars_az;
|
||||
```
|
||||
</TabItem>
|
||||
<TabItem label="DE (new)">
|
||||
```sql
|
||||
-- Uses DE if configured, falls back to VSOP87 otherwise. STABLE.
|
||||
SELECT topo_azimuth(planet_observe_de(4, '40.0N 105.3W 1655m'::observer, now())) AS mars_az;
|
||||
```
|
||||
</TabItem>
|
||||
</Tabs>
|
||||
|
||||
The `_de()` variants share the same parameter types, return types, and body ID conventions. You can swap between them by adding or removing `_de` from the function name.
|
||||
|
||||
### All planets with DE positions
|
||||
|
||||
```sql
|
||||
SELECT body_id,
|
||||
CASE body_id
|
||||
WHEN 1 THEN 'Mercury' WHEN 2 THEN 'Venus'
|
||||
WHEN 4 THEN 'Mars' WHEN 5 THEN 'Jupiter'
|
||||
WHEN 6 THEN 'Saturn' WHEN 7 THEN 'Uranus'
|
||||
WHEN 8 THEN 'Neptune'
|
||||
END AS planet,
|
||||
round(topo_azimuth(planet_observe_de(body_id, '40.0N 105.3W 1655m'::observer, now()))::numeric, 4) AS az,
|
||||
round(topo_elevation(planet_observe_de(body_id, '40.0N 105.3W 1655m'::observer, now()))::numeric, 4) AS el
|
||||
FROM unnest(ARRAY[1,2,4,5,6,7,8]) AS body_id
|
||||
WHERE topo_elevation(planet_observe_de(body_id, '40.0N 105.3W 1655m'::observer, now())) > 0
|
||||
ORDER BY el DESC;
|
||||
```
|
||||
|
||||
### Moon via DE
|
||||
|
||||
```sql
|
||||
SELECT round(topo_azimuth(moon_observe_de('40.0N 105.3W 1655m'::observer, now()))::numeric, 2) AS az,
|
||||
round(topo_elevation(moon_observe_de('40.0N 105.3W 1655m'::observer, now()))::numeric, 2) AS el,
|
||||
round(topo_range(moon_observe_de('40.0N 105.3W 1655m'::observer, now()))::numeric, 0) AS range_km;
|
||||
```
|
||||
|
||||
### Lambert transfers with DE positions
|
||||
|
||||
```sql
|
||||
-- Earth-Mars transfer with DE-quality planet positions
|
||||
SELECT * FROM lambert_transfer_de(3, 4, '2026-05-01 00:00:00+00', '2027-01-15 00:00:00+00');
|
||||
```
|
||||
|
||||
### Cross-validation: DE vs VSOP87
|
||||
|
||||
Compare DE and VSOP87 results to see the sub-arcsecond difference:
|
||||
|
||||
```sql
|
||||
SELECT body_id,
|
||||
round(helio_distance(planet_heliocentric(body_id, '2024-06-21 12:00:00+00'))::numeric, 10) AS vsop87_r,
|
||||
round(helio_distance(planet_heliocentric_de(body_id, '2024-06-21 12:00:00+00'))::numeric, 10) AS de_r,
|
||||
round((helio_distance(planet_heliocentric_de(body_id, '2024-06-21 12:00:00+00'))
|
||||
- helio_distance(planet_heliocentric(body_id, '2024-06-21 12:00:00+00')))::numeric * 149597870.7, 2) AS diff_km
|
||||
FROM generate_series(1, 8) AS body_id;
|
||||
```
|
||||
|
||||
## Accuracy comparison
|
||||
|
||||
| Provider | Theory | Accuracy | Volatility | Data Source |
|
||||
|----------|--------|----------|------------|-------------|
|
||||
| VSOP87 | Fourier series (Bretagnon 1988) | ~1 arcsecond (inner planets), up to 40 arcseconds (Neptune, +/- 4000 yr) | IMMUTABLE | Compiled-in coefficients |
|
||||
| ELP2000-82B | Fourier series (Chapront 1983) | ~2 arcseconds (longitude) | IMMUTABLE | Compiled-in coefficients |
|
||||
| JPL DE440 | Numerical integration | ~0.1 milliarcsecond | STABLE | External binary file (115 MB) |
|
||||
| JPL DE441 | Numerical integration | ~0.1 milliarcsecond | STABLE | External binary file (3.1 GB) |
|
||||
|
||||
<Aside type="note" title="The volatility trade-off">
|
||||
IMMUTABLE functions can be used in expression indexes and are constant-folded during planning. STABLE functions execute once per row per statement. If you have expression indexes on VSOP87 functions (e.g., almanac tables), keep them on the non-DE variants. Use DE functions for queries where accuracy matters more than planning-time optimization.
|
||||
</Aside>
|
||||
|
||||
## Diagnostics
|
||||
|
||||
The `pg_orbit_ephemeris_info()` function reports the current state:
|
||||
|
||||
```sql
|
||||
SELECT * FROM pg_orbit_ephemeris_info();
|
||||
```
|
||||
|
||||
| Column | Type | Description |
|
||||
|--------|------|-------------|
|
||||
| `provider` | text | `'VSOP87'` or `'JPL_DE'` |
|
||||
| `file_path` | text | Path to DE file (empty if VSOP87-only) |
|
||||
| `start_jd` | float8 | First Julian Date in the DE file |
|
||||
| `end_jd` | float8 | Last Julian Date in the DE file |
|
||||
| `version` | int4 | DE version number (440, 441, etc.) |
|
||||
| `au_km` | float8 | AU value from the DE header |
|
||||
|
||||
When no DE file is configured, `provider` returns `'VSOP87'` and the remaining columns are defaults.
|
||||
|
||||
## Fallback behavior
|
||||
|
||||
DE functions never fail silently. The fallback rules:
|
||||
|
||||
1. **No GUC set** (default): Fall back to VSOP87/ELP2000-82B silently. No NOTICE, no overhead.
|
||||
2. **GUC set, file works**: Use DE positions.
|
||||
3. **GUC set, query-time failure** (JD out of range, I/O error): Emit `NOTICE` explaining the fallback, then use VSOP87/ELP2000-82B.
|
||||
|
||||
This means `_de()` functions always return a valid result. They never abort due to DE-specific issues --- at worst, they degrade to VSOP87 accuracy with a NOTICE explaining why.
|
||||
@ -40,7 +40,7 @@ Body IDs follow the VSOP87 convention: 1=Mercury, 2=Venus, 3=Earth, 4=Mars, 5=Ju
|
||||
VSOP87 and ELP2000-82B are analytic theories. They trade the last bits of precision for computational speed and zero external data dependencies.
|
||||
</Aside>
|
||||
|
||||
- **VSOP87 accuracy is about 1 arcsecond.** JPL DE441 (used by Skyfield and SPICE) achieves 0.001 arcsecond. For visual observation planning, 1 arcsecond is more than sufficient. For pointing a dish at GHz frequencies or precision astrometry, use SPICE.
|
||||
- **VSOP87 accuracy is about 1 arcsecond.** JPL DE441 (used by Skyfield and SPICE) achieves 0.001 arcsecond. For visual observation planning, 1 arcsecond is more than sufficient. For precision astrometry or GHz dish pointing, pg_orbit v0.3.0 supports [optional DE440/441 ephemeris files](/guides/de-ephemeris/) with sub-milliarcsecond accuracy.
|
||||
- **ELP2000-82B accuracy is about 10 arcseconds** for the Moon. Good enough for knowing when the Moon is up, what phase it is in, and whether it will interfere with observations. Not sufficient for occultation timing.
|
||||
- **No light-time iteration.** pg_orbit computes geometric positions, not apparent positions. The difference matters at the milliarcsecond level.
|
||||
- **No atmospheric refraction.** Objects near the horizon appear slightly higher than their geometric position. pg_orbit does not apply refraction corrections.
|
||||
|
||||
@ -6,7 +6,7 @@ sidebar:
|
||||
|
||||
import { Aside, Tabs, TabItem } from "@astrojs/starlight/components";
|
||||
|
||||
Physical constants, astronomical constants, and accuracy bounds for every computational theory used in pg_orbit. All constants are compiled from their original sources and embedded at compile time -- no runtime configuration files, no external data dependencies.
|
||||
Physical constants, astronomical constants, and accuracy bounds for every computational theory used in pg_orbit. VSOP87/ELP2000-82B constants are compiled from their original sources and embedded at compile time. Optional JPL DE440/441 ephemeris support (v0.3.0+) reads constants from an external binary file.
|
||||
|
||||
---
|
||||
|
||||
@ -104,6 +104,32 @@ For dates within a few centuries of J2000, VSOP87 is accurate to sub-arcsecond l
|
||||
|
||||
**Valid epoch range:** The series coefficients are fitted over the interval J2000 +/- 4000 years. Outside this range, accuracy degrades unpredictably. For observational planning (present-day queries), accuracy is well within 1 arcsecond.
|
||||
|
||||
### JPL DE440/441 (Optional, v0.3.0+)
|
||||
|
||||
When configured via the `pg_orbit.ephemeris_path` GUC, DE ephemeris positions are available through the `_de()` function variants.
|
||||
|
||||
| Property | DE440 | DE441 |
|
||||
|----------|-------|-------|
|
||||
| Accuracy | ~0.1 milliarcsecond | ~0.1 milliarcsecond |
|
||||
| Valid range | 1550 to 2650 | -13200 to +17191 |
|
||||
| File size | ~115 MB | ~3.1 GB |
|
||||
| Method | Numerical integration | Numerical integration |
|
||||
|
||||
Accuracy relative to VLBI observations:
|
||||
|
||||
| Body | Position Accuracy | Notes |
|
||||
|------|-------------------|-------|
|
||||
| Inner planets (Mercury-Mars) | < 0.3 milliarcsecond | Constrained by radar ranging |
|
||||
| Jupiter, Saturn | < 1 milliarcsecond | Constrained by spacecraft tracking + CCD astrometry |
|
||||
| Uranus, Neptune | < 10 milliarcsecond | Limited by ground-based astrometry |
|
||||
| Moon | ~1 meter | Constrained by Lunar Laser Ranging |
|
||||
|
||||
<Aside type="note">
|
||||
DE positions are in the ICRS equatorial frame. pg_orbit applies the equatorial-to-ecliptic rotation at the provider boundary before feeding positions into the observation pipeline. Both target and Earth positions always come from the same provider (Rule 7 of the constant chain of custody).
|
||||
</Aside>
|
||||
|
||||
**DE440 vs DE441:** DE440 is recommended for present-day and near-future work. It covers 1550-2650 and produces identical results to DE441 within that range, but is ~27x smaller. DE441's extended range (-13200 to +17191) is needed only for paleoclimate studies, historical astronomy, or very long-term solar system dynamics.
|
||||
|
||||
### ELP2000-82B (Lunar Position)
|
||||
|
||||
| Quantity | Accuracy |
|
||||
|
||||
339
docs/src/content/docs/reference/functions-de.mdx
Normal file
339
docs/src/content/docs/reference/functions-de.mdx
Normal file
@ -0,0 +1,339 @@
|
||||
---
|
||||
title: "Functions: DE Ephemeris"
|
||||
sidebar:
|
||||
order: 8
|
||||
---
|
||||
|
||||
import { Aside, Tabs, TabItem } from "@astrojs/starlight/components";
|
||||
|
||||
Optional high-precision function variants that use JPL DE440/441 ephemeris files when configured via the `pg_orbit.ephemeris_path` GUC. Each function mirrors an existing VSOP87/ELP2000-82B counterpart with identical parameters and return types.
|
||||
|
||||
All DE functions are `STABLE STRICT PARALLEL SAFE` (except `pg_orbit_ephemeris_info` which is `STABLE PARALLEL SAFE`, not STRICT). When DE is unavailable, they fall back to their VSOP87/ELP2000-82B equivalents.
|
||||
|
||||
See the [DE Ephemeris guide](/guides/de-ephemeris/) for setup and configuration.
|
||||
|
||||
---
|
||||
|
||||
## planet_heliocentric_de
|
||||
|
||||
Computes the heliocentric ecliptic J2000 position of a planet using DE positions when available, falling back to VSOP87.
|
||||
|
||||
### Signature
|
||||
|
||||
```sql
|
||||
planet_heliocentric_de(body_id int4, t timestamptz) → heliocentric
|
||||
```
|
||||
|
||||
### Parameters
|
||||
|
||||
| Parameter | Type | Description |
|
||||
|-----------|------|-------------|
|
||||
| `body_id` | `int4` | Planet identifier: 0 (Sun), 1-8 (Mercury through Neptune) |
|
||||
| `t` | `timestamptz` | Evaluation time |
|
||||
|
||||
### Returns
|
||||
|
||||
A `heliocentric` position in AU (ecliptic J2000 frame). Identical return type to `planet_heliocentric()`.
|
||||
|
||||
### Example
|
||||
|
||||
```sql
|
||||
-- Compare DE vs VSOP87 heliocentric distances
|
||||
SELECT body_id,
|
||||
round(helio_distance(planet_heliocentric(body_id, '2024-06-21 12:00:00+00'))::numeric, 10) AS vsop87,
|
||||
round(helio_distance(planet_heliocentric_de(body_id, '2024-06-21 12:00:00+00'))::numeric, 10) AS de
|
||||
FROM generate_series(1, 8) AS body_id;
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## planet_observe_de
|
||||
|
||||
Computes the topocentric position of a planet using DE ephemeris for both the target and Earth positions.
|
||||
|
||||
### Signature
|
||||
|
||||
```sql
|
||||
planet_observe_de(body_id int4, obs observer, t timestamptz) → topocentric
|
||||
```
|
||||
|
||||
### Parameters
|
||||
|
||||
| Parameter | Type | Description |
|
||||
|-----------|------|-------------|
|
||||
| `body_id` | `int4` | Planet identifier (1-8, excluding 3/Earth) |
|
||||
| `obs` | `observer` | Observer location on Earth |
|
||||
| `t` | `timestamptz` | Observation time |
|
||||
|
||||
<Aside type="caution">
|
||||
Body IDs 0 (Sun) and 3 (Earth) are not valid. Use `sun_observe_de()` for the Sun. Observing Earth from Earth raises an error.
|
||||
</Aside>
|
||||
|
||||
### Returns
|
||||
|
||||
A `topocentric` with azimuth, elevation, range (km), and range rate (km/s).
|
||||
|
||||
### Example
|
||||
|
||||
```sql
|
||||
-- Mars position from Boulder using DE
|
||||
SELECT round(topo_azimuth(t)::numeric, 4) AS az,
|
||||
round(topo_elevation(t)::numeric, 4) AS el,
|
||||
round(topo_range(t)::numeric, 0) AS range_km
|
||||
FROM planet_observe_de(4, '40.0N 105.3W 1655m'::observer, now()) AS t;
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## sun_observe_de
|
||||
|
||||
Computes the topocentric position of the Sun using DE positions.
|
||||
|
||||
### Signature
|
||||
|
||||
```sql
|
||||
sun_observe_de(obs observer, t timestamptz) → topocentric
|
||||
```
|
||||
|
||||
### Parameters
|
||||
|
||||
| Parameter | Type | Description |
|
||||
|-----------|------|-------------|
|
||||
| `obs` | `observer` | Observer location on Earth |
|
||||
| `t` | `timestamptz` | Observation time |
|
||||
|
||||
### Returns
|
||||
|
||||
A `topocentric` with azimuth, elevation, range (km), and range rate (km/s).
|
||||
|
||||
### Example
|
||||
|
||||
```sql
|
||||
SELECT round(topo_azimuth(t)::numeric, 2) AS az,
|
||||
round(topo_elevation(t)::numeric, 2) AS el
|
||||
FROM sun_observe_de('40.0N 105.3W 1655m'::observer, now()) AS t;
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## moon_observe_de
|
||||
|
||||
Computes the topocentric position of the Moon using DE positions. Falls back to ELP2000-82B when DE is unavailable.
|
||||
|
||||
### Signature
|
||||
|
||||
```sql
|
||||
moon_observe_de(obs observer, t timestamptz) → topocentric
|
||||
```
|
||||
|
||||
### Parameters
|
||||
|
||||
| Parameter | Type | Description |
|
||||
|-----------|------|-------------|
|
||||
| `obs` | `observer` | Observer location on Earth |
|
||||
| `t` | `timestamptz` | Observation time |
|
||||
|
||||
### Returns
|
||||
|
||||
A `topocentric` with azimuth, elevation, range (km), and range rate (km/s).
|
||||
|
||||
### Example
|
||||
|
||||
```sql
|
||||
SELECT round(topo_azimuth(t)::numeric, 2) AS az,
|
||||
round(topo_elevation(t)::numeric, 2) AS el,
|
||||
round(topo_range(t)::numeric, 0) AS range_km
|
||||
FROM moon_observe_de('40.0N 105.3W 1655m'::observer, now()) AS t;
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## lambert_transfer_de
|
||||
|
||||
Computes a Lambert transfer orbit using DE-precision planet positions.
|
||||
|
||||
### Signature
|
||||
|
||||
```sql
|
||||
lambert_transfer_de(
|
||||
dep_body_id int4,
|
||||
arr_body_id int4,
|
||||
dep_time timestamptz,
|
||||
arr_time timestamptz
|
||||
) → RECORD(c3_departure float8, c3_arrival float8, v_inf_dep float8, v_inf_arr float8, tof_days float8)
|
||||
```
|
||||
|
||||
### Parameters
|
||||
|
||||
| Parameter | Type | Description |
|
||||
|-----------|------|-------------|
|
||||
| `dep_body_id` | `int4` | Departure planet (1-8) |
|
||||
| `arr_body_id` | `int4` | Arrival planet (1-8, different from departure) |
|
||||
| `dep_time` | `timestamptz` | Departure epoch |
|
||||
| `arr_time` | `timestamptz` | Arrival epoch (must be after departure) |
|
||||
|
||||
### Returns
|
||||
|
||||
Transfer orbit parameters including departure and arrival C3 (km^2/s^2), v-infinity magnitudes, and time of flight.
|
||||
|
||||
### Example
|
||||
|
||||
```sql
|
||||
-- Earth-Mars transfer window with DE positions
|
||||
SELECT round(c3_departure::numeric, 2) AS c3_dep,
|
||||
round(c3_arrival::numeric, 2) AS c3_arr,
|
||||
round(tof_days::numeric, 1) AS tof
|
||||
FROM lambert_transfer_de(3, 4, '2026-05-01 00:00:00+00', '2027-01-15 00:00:00+00');
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## lambert_c3_de
|
||||
|
||||
Returns only the departure C3 (km^2/s^2) from a Lambert transfer computation using DE positions. A convenience wrapper around `lambert_transfer_de()`.
|
||||
|
||||
### Signature
|
||||
|
||||
```sql
|
||||
lambert_c3_de(dep_body_id int4, arr_body_id int4, dep_time timestamptz, arr_time timestamptz) → float8
|
||||
```
|
||||
|
||||
### Example
|
||||
|
||||
```sql
|
||||
-- Pork chop grid with DE accuracy
|
||||
SELECT dep::date AS departure,
|
||||
arr::date AS arrival,
|
||||
round(lambert_c3_de(3, 4, dep, arr)::numeric, 2) AS c3
|
||||
FROM generate_series('2026-04-01'::timestamptz, '2026-08-01'::timestamptz, '14 days') dep,
|
||||
generate_series('2026-12-01'::timestamptz, '2027-04-01'::timestamptz, '14 days') arr
|
||||
WHERE arr > dep + interval '90 days';
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## galilean_observe_de
|
||||
|
||||
Computes the topocentric position of a Galilean moon of Jupiter. Uses DE for Jupiter's heliocentric position and L1.2 theory for the moon's offset.
|
||||
|
||||
### Signature
|
||||
|
||||
```sql
|
||||
galilean_observe_de(moon_id int4, obs observer, t timestamptz) → topocentric
|
||||
```
|
||||
|
||||
### Parameters
|
||||
|
||||
| Parameter | Type | Description |
|
||||
|-----------|------|-------------|
|
||||
| `moon_id` | `int4` | 0=Io, 1=Europa, 2=Ganymede, 3=Callisto |
|
||||
| `obs` | `observer` | Observer location |
|
||||
| `t` | `timestamptz` | Observation time |
|
||||
|
||||
### Example
|
||||
|
||||
```sql
|
||||
-- All four Galilean moons with DE-precision Jupiter
|
||||
SELECT moon_id,
|
||||
CASE moon_id WHEN 0 THEN 'Io' WHEN 1 THEN 'Europa'
|
||||
WHEN 2 THEN 'Ganymede' WHEN 3 THEN 'Callisto' END AS moon,
|
||||
round(topo_elevation(galilean_observe_de(moon_id, '40.0N 105.3W 1655m'::observer, now()))::numeric, 2) AS el
|
||||
FROM generate_series(0, 3) AS moon_id;
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## saturn_moon_observe_de
|
||||
|
||||
Computes the topocentric position of a Saturn moon. Uses DE for Saturn's heliocentric position and TASS17 theory for the moon's offset.
|
||||
|
||||
### Signature
|
||||
|
||||
```sql
|
||||
saturn_moon_observe_de(moon_id int4, obs observer, t timestamptz) → topocentric
|
||||
```
|
||||
|
||||
### Parameters
|
||||
|
||||
| Parameter | Type | Description |
|
||||
|-----------|------|-------------|
|
||||
| `moon_id` | `int4` | 0=Mimas, 1=Enceladus, 2=Tethys, 3=Dione, 4=Rhea, 5=Titan, 6=Iapetus, 7=Hyperion |
|
||||
| `obs` | `observer` | Observer location |
|
||||
| `t` | `timestamptz` | Observation time |
|
||||
|
||||
---
|
||||
|
||||
## uranus_moon_observe_de
|
||||
|
||||
Computes the topocentric position of a Uranus moon. Uses DE for Uranus's heliocentric position and GUST86 theory for the moon's offset.
|
||||
|
||||
### Signature
|
||||
|
||||
```sql
|
||||
uranus_moon_observe_de(moon_id int4, obs observer, t timestamptz) → topocentric
|
||||
```
|
||||
|
||||
### Parameters
|
||||
|
||||
| Parameter | Type | Description |
|
||||
|-----------|------|-------------|
|
||||
| `moon_id` | `int4` | 0=Miranda, 1=Ariel, 2=Umbriel, 3=Titania, 4=Oberon |
|
||||
| `obs` | `observer` | Observer location |
|
||||
| `t` | `timestamptz` | Observation time |
|
||||
|
||||
---
|
||||
|
||||
## mars_moon_observe_de
|
||||
|
||||
Computes the topocentric position of a Mars moon. Uses DE for Mars's heliocentric position and MarsSat theory for the moon's offset.
|
||||
|
||||
### Signature
|
||||
|
||||
```sql
|
||||
mars_moon_observe_de(moon_id int4, obs observer, t timestamptz) → topocentric
|
||||
```
|
||||
|
||||
### Parameters
|
||||
|
||||
| Parameter | Type | Description |
|
||||
|-----------|------|-------------|
|
||||
| `moon_id` | `int4` | 0=Phobos, 1=Deimos |
|
||||
| `obs` | `observer` | Observer location |
|
||||
| `t` | `timestamptz` | Observation time |
|
||||
|
||||
---
|
||||
|
||||
## pg_orbit_ephemeris_info
|
||||
|
||||
Returns diagnostic information about the current ephemeris provider.
|
||||
|
||||
### Signature
|
||||
|
||||
```sql
|
||||
pg_orbit_ephemeris_info() → RECORD(provider text, file_path text, start_jd float8, end_jd float8, version int4, au_km float8)
|
||||
```
|
||||
|
||||
<Aside type="note">
|
||||
This function is `STABLE PARALLEL SAFE` but **not** `STRICT` --- it takes no arguments and always returns a row.
|
||||
</Aside>
|
||||
|
||||
### Returns
|
||||
|
||||
| Column | Type | Description |
|
||||
|--------|------|-------------|
|
||||
| `provider` | `text` | `'VSOP87'` or `'JPL_DE'` |
|
||||
| `file_path` | `text` | Path to the DE file (empty string if no DE) |
|
||||
| `start_jd` | `float8` | First Julian Date covered by the DE file |
|
||||
| `end_jd` | `float8` | Last Julian Date covered by the DE file |
|
||||
| `version` | `int4` | DE version number (440, 441, etc.) |
|
||||
| `au_km` | `float8` | Astronomical Unit in km from the DE header |
|
||||
|
||||
### Example
|
||||
|
||||
```sql
|
||||
-- Check current provider
|
||||
SELECT (pg_orbit_ephemeris_info()).provider;
|
||||
|
||||
-- Full diagnostic
|
||||
SELECT * FROM pg_orbit_ephemeris_info();
|
||||
```
|
||||
@ -6,7 +6,9 @@ sidebar:
|
||||
|
||||
import { Aside, Tabs, TabItem } from "@astrojs/starlight/components";
|
||||
|
||||
Functions for computing planetary positions, observing the Sun, Moon, and planets from an Earth-based observer. Planetary positions use the VSOP87 theory (Bretagnon & Francou, 1988). Lunar position uses ELP2000-82B (Chapront-Touze & Chapront, 1983).
|
||||
Functions for computing planetary positions, observing the Sun, Moon, and planets from an Earth-based observer. Planetary positions use the VSOP87 theory (Bretagnon & Francou, 1988). Lunar position uses ELP2000-82B (Chapront-Touze & Chapront, 1983). All functions are `IMMUTABLE STRICT PARALLEL SAFE`.
|
||||
|
||||
For higher precision, v0.3.0 adds optional `_de()` variants that use JPL DE440/441 ephemeris files. See [Functions: DE Ephemeris](/reference/functions-de/).
|
||||
|
||||
---
|
||||
|
||||
|
||||
@ -1,4 +1,4 @@
|
||||
comment = 'Orbital mechanics types and functions for PostgreSQL'
|
||||
default_version = '0.2.0'
|
||||
default_version = '0.3.0'
|
||||
module_pathname = '$libdir/pg_orbit'
|
||||
relocatable = true
|
||||
|
||||
86
sql/pg_orbit--0.2.0--0.3.0.sql
Normal file
86
sql/pg_orbit--0.2.0--0.3.0.sql
Normal file
@ -0,0 +1,86 @@
|
||||
-- pg_orbit 0.2.0 -> 0.3.0 migration
|
||||
--
|
||||
-- Adds optional JPL DE440/441 ephemeris functions.
|
||||
-- Existing VSOP87/ELP2000-82B functions are unchanged (still IMMUTABLE).
|
||||
-- New _de() functions are STABLE (depend on external DE binary file).
|
||||
-- When DE is unavailable, _de() functions fall back to VSOP87 silently.
|
||||
|
||||
-- ============================================================
|
||||
-- Phase 5: DE ephemeris functions (optional high-precision)
|
||||
-- ============================================================
|
||||
|
||||
-- Planet observation with DE ephemeris
|
||||
|
||||
CREATE FUNCTION planet_heliocentric_de(int4, timestamptz) RETURNS heliocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION planet_heliocentric_de(int4, timestamptz) IS
|
||||
'Heliocentric ecliptic J2000 position via JPL DE (sub-arcsecond). Falls back to VSOP87 if DE unavailable. Body IDs: 0=Sun, 1-8=Mercury-Neptune.';
|
||||
|
||||
CREATE FUNCTION planet_observe_de(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION planet_observe_de(int4, observer, timestamptz) IS
|
||||
'Observe planet via JPL DE. Falls back to VSOP87. Body IDs: 1-8 (Mercury-Neptune).';
|
||||
|
||||
CREATE FUNCTION sun_observe_de(observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION sun_observe_de(observer, timestamptz) IS
|
||||
'Observe Sun via JPL DE. Falls back to VSOP87.';
|
||||
|
||||
CREATE FUNCTION moon_observe_de(observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION moon_observe_de(observer, timestamptz) IS
|
||||
'Observe Moon via JPL DE. Falls back to ELP2000-82B.';
|
||||
|
||||
|
||||
-- Lambert transfer with DE positions
|
||||
|
||||
CREATE FUNCTION lambert_transfer_de(
|
||||
dep_body_id int4, arr_body_id int4,
|
||||
dep_time timestamptz, arr_time timestamptz,
|
||||
OUT c3_departure float8, OUT c3_arrival float8,
|
||||
OUT v_inf_departure float8, OUT v_inf_arrival float8,
|
||||
OUT tof_days float8, OUT transfer_sma float8
|
||||
) RETURNS RECORD
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lambert_transfer_de(int4, int4, timestamptz, timestamptz) IS
|
||||
'Lambert transfer via JPL DE positions. Falls back to VSOP87. Body IDs 1-8.';
|
||||
|
||||
CREATE FUNCTION lambert_c3_de(int4, int4, timestamptz, timestamptz) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lambert_c3_de(int4, int4, timestamptz, timestamptz) IS
|
||||
'Departure C3 via JPL DE. Falls back to VSOP87. For pork chop plots.';
|
||||
|
||||
|
||||
-- Planetary moon observation with DE parent positions
|
||||
|
||||
CREATE FUNCTION galilean_observe_de(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION galilean_observe_de(int4, observer, timestamptz) IS
|
||||
'Observe Galilean moon with JPL DE parent position. L1.2 moon offsets. Body IDs: 0-3 (Io-Callisto).';
|
||||
|
||||
CREATE FUNCTION saturn_moon_observe_de(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION saturn_moon_observe_de(int4, observer, timestamptz) IS
|
||||
'Observe Saturn moon with JPL DE parent position. TASS 1.7 moon offsets. Body IDs: 0-7 (Mimas-Hyperion).';
|
||||
|
||||
CREATE FUNCTION uranus_moon_observe_de(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION uranus_moon_observe_de(int4, observer, timestamptz) IS
|
||||
'Observe Uranus moon with JPL DE parent position. GUST86 moon offsets. Body IDs: 0-4 (Miranda-Oberon).';
|
||||
|
||||
CREATE FUNCTION mars_moon_observe_de(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION mars_moon_observe_de(int4, observer, timestamptz) IS
|
||||
'Observe Mars moon with JPL DE parent position. MarsSat moon offsets. Body IDs: 0-1 (Phobos-Deimos).';
|
||||
|
||||
|
||||
-- Diagnostic function
|
||||
|
||||
CREATE FUNCTION pg_orbit_ephemeris_info(
|
||||
OUT provider text, OUT file_path text,
|
||||
OUT start_jd float8, OUT end_jd float8,
|
||||
OUT version int4, OUT au_km float8
|
||||
) RETURNS RECORD
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION pg_orbit_ephemeris_info() IS
|
||||
'Returns current ephemeris provider status: VSOP87 or JPL_DE with file path, JD range, version, and AU value.';
|
||||
790
sql/pg_orbit--0.3.0.sql
Normal file
790
sql/pg_orbit--0.3.0.sql
Normal file
@ -0,0 +1,790 @@
|
||||
-- pg_orbit -- Orbital mechanics types and functions for PostgreSQL
|
||||
--
|
||||
-- Types: tle, eci_position, geodetic, topocentric, observer, pass_event
|
||||
-- Provides SGP4/SDP4 propagation, coordinate transforms, pass prediction,
|
||||
-- and GiST indexing on altitude bands for conjunction screening.
|
||||
--
|
||||
-- All propagation uses WGS-72 constants (matching TLE mean element fitting).
|
||||
-- Coordinate output uses WGS-84 (matching modern geodetic standards).
|
||||
|
||||
-- ============================================================
|
||||
-- Shell types (forward declarations)
|
||||
-- ============================================================
|
||||
|
||||
CREATE TYPE tle;
|
||||
CREATE TYPE eci_position;
|
||||
CREATE TYPE geodetic;
|
||||
CREATE TYPE topocentric;
|
||||
CREATE TYPE observer;
|
||||
CREATE TYPE pass_event;
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- TLE type: Two-Line Element set
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION tle_in(cstring) RETURNS tle
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION tle_out(tle) RETURNS cstring
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION tle_recv(internal) RETURNS tle
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION tle_send(tle) RETURNS bytea
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE TYPE tle (
|
||||
INPUT = tle_in,
|
||||
OUTPUT = tle_out,
|
||||
RECEIVE = tle_recv,
|
||||
SEND = tle_send,
|
||||
INTERNALLENGTH = 112,
|
||||
ALIGNMENT = double,
|
||||
STORAGE = plain
|
||||
);
|
||||
|
||||
COMMENT ON TYPE tle IS 'Two-Line Element set — parsed mean orbital elements for SGP4/SDP4 propagation';
|
||||
|
||||
-- TLE accessor functions
|
||||
|
||||
CREATE FUNCTION tle_epoch(tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_epoch(tle) IS 'TLE epoch as Julian date (UTC)';
|
||||
|
||||
CREATE FUNCTION tle_norad_id(tle) RETURNS int4
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_norad_id(tle) IS 'NORAD catalog number';
|
||||
|
||||
CREATE FUNCTION tle_inclination(tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_inclination(tle) IS 'Orbital inclination in degrees';
|
||||
|
||||
CREATE FUNCTION tle_eccentricity(tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_eccentricity(tle) IS 'Orbital eccentricity (dimensionless)';
|
||||
|
||||
CREATE FUNCTION tle_raan(tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_raan(tle) IS 'Right ascension of ascending node in degrees';
|
||||
|
||||
CREATE FUNCTION tle_arg_perigee(tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_arg_perigee(tle) IS 'Argument of perigee in degrees';
|
||||
|
||||
CREATE FUNCTION tle_mean_anomaly(tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_mean_anomaly(tle) IS 'Mean anomaly in degrees';
|
||||
|
||||
CREATE FUNCTION tle_mean_motion(tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_mean_motion(tle) IS 'Mean motion in revolutions per day';
|
||||
|
||||
CREATE FUNCTION tle_bstar(tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_bstar(tle) IS 'B* drag coefficient (1/earth-radii)';
|
||||
|
||||
CREATE FUNCTION tle_period(tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_period(tle) IS 'Orbital period in minutes';
|
||||
|
||||
CREATE FUNCTION tle_age(tle, timestamptz) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_age(tle, timestamptz) IS 'TLE age in days (positive = past epoch, negative = before epoch)';
|
||||
|
||||
CREATE FUNCTION tle_perigee(tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_perigee(tle) IS 'Perigee altitude in km above WGS-72 ellipsoid';
|
||||
|
||||
CREATE FUNCTION tle_apogee(tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_apogee(tle) IS 'Apogee altitude in km above WGS-72 ellipsoid';
|
||||
|
||||
CREATE FUNCTION tle_intl_desig(tle) RETURNS text
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_intl_desig(tle) IS 'International designator (COSPAR ID)';
|
||||
|
||||
CREATE FUNCTION tle_from_lines(text, text) RETURNS tle
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_from_lines(text, text) IS
|
||||
'Construct TLE from separate line1/line2 text columns';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- ECI position type: True Equator Mean Equinox (TEME) frame
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION eci_in(cstring) RETURNS eci_position
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION eci_out(eci_position) RETURNS cstring
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION eci_recv(internal) RETURNS eci_position
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION eci_send(eci_position) RETURNS bytea
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE TYPE eci_position (
|
||||
INPUT = eci_in,
|
||||
OUTPUT = eci_out,
|
||||
RECEIVE = eci_recv,
|
||||
SEND = eci_send,
|
||||
INTERNALLENGTH = 48,
|
||||
ALIGNMENT = double,
|
||||
STORAGE = plain
|
||||
);
|
||||
|
||||
COMMENT ON TYPE eci_position IS 'Earth-Centered Inertial position and velocity in TEME frame (km, km/s)';
|
||||
|
||||
-- ECI accessor functions
|
||||
|
||||
CREATE FUNCTION eci_x(eci_position) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION eci_y(eci_position) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION eci_z(eci_position) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION eci_vx(eci_position) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION eci_vy(eci_position) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION eci_vz(eci_position) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION eci_speed(eci_position) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION eci_speed(eci_position) IS 'Velocity magnitude in km/s';
|
||||
|
||||
CREATE FUNCTION eci_altitude(eci_position) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION eci_altitude(eci_position) IS 'Approximate geocentric altitude in km (radius - WGS72_AE)';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Geodetic type: WGS-84 latitude/longitude/altitude
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION geodetic_in(cstring) RETURNS geodetic
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION geodetic_out(geodetic) RETURNS cstring
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION geodetic_recv(internal) RETURNS geodetic
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION geodetic_send(geodetic) RETURNS bytea
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE TYPE geodetic (
|
||||
INPUT = geodetic_in,
|
||||
OUTPUT = geodetic_out,
|
||||
RECEIVE = geodetic_recv,
|
||||
SEND = geodetic_send,
|
||||
INTERNALLENGTH = 24,
|
||||
ALIGNMENT = double,
|
||||
STORAGE = plain
|
||||
);
|
||||
|
||||
COMMENT ON TYPE geodetic IS 'Geodetic coordinates on WGS-84 ellipsoid (lat/lon in degrees, altitude in km)';
|
||||
|
||||
CREATE FUNCTION geodetic_lat(geodetic) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION geodetic_lon(geodetic) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION geodetic_alt(geodetic) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Topocentric type: observer-relative az/el/range
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION topocentric_in(cstring) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION topocentric_out(topocentric) RETURNS cstring
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION topocentric_recv(internal) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION topocentric_send(topocentric) RETURNS bytea
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE TYPE topocentric (
|
||||
INPUT = topocentric_in,
|
||||
OUTPUT = topocentric_out,
|
||||
RECEIVE = topocentric_recv,
|
||||
SEND = topocentric_send,
|
||||
INTERNALLENGTH = 32,
|
||||
ALIGNMENT = double,
|
||||
STORAGE = plain
|
||||
);
|
||||
|
||||
COMMENT ON TYPE topocentric IS 'Topocentric coordinates relative to observer (azimuth, elevation, range, range rate)';
|
||||
|
||||
CREATE FUNCTION topo_azimuth(topocentric) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION topo_azimuth(topocentric) IS 'Azimuth in degrees (0=N, 90=E, 180=S, 270=W)';
|
||||
|
||||
CREATE FUNCTION topo_elevation(topocentric) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION topo_elevation(topocentric) IS 'Elevation in degrees (0=horizon, 90=zenith)';
|
||||
|
||||
CREATE FUNCTION topo_range(topocentric) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION topo_range(topocentric) IS 'Slant range in km';
|
||||
|
||||
CREATE FUNCTION topo_range_rate(topocentric) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION topo_range_rate(topocentric) IS 'Range rate in km/s (positive = receding)';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Observer type: ground station location
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION observer_in(cstring) RETURNS observer
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION observer_out(observer) RETURNS cstring
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION observer_recv(internal) RETURNS observer
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION observer_send(observer) RETURNS bytea
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE TYPE observer (
|
||||
INPUT = observer_in,
|
||||
OUTPUT = observer_out,
|
||||
RECEIVE = observer_recv,
|
||||
SEND = observer_send,
|
||||
INTERNALLENGTH = 24,
|
||||
ALIGNMENT = double,
|
||||
STORAGE = plain
|
||||
);
|
||||
|
||||
COMMENT ON TYPE observer IS 'Ground observer location (accepts: 40.0N 105.3W 1655m or decimal degrees)';
|
||||
|
||||
CREATE FUNCTION observer_lat(observer) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION observer_lat(observer) IS 'Latitude in degrees (positive = North)';
|
||||
|
||||
CREATE FUNCTION observer_lon(observer) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION observer_lon(observer) IS 'Longitude in degrees (positive = East)';
|
||||
|
||||
CREATE FUNCTION observer_alt(observer) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION observer_alt(observer) IS 'Altitude in meters above WGS-84 ellipsoid';
|
||||
|
||||
CREATE FUNCTION observer_from_geodetic(float8, float8, float8 DEFAULT 0.0) RETURNS observer
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION observer_from_geodetic(float8, float8, float8) IS
|
||||
'Construct observer from lat (deg), lon (deg), altitude (meters). Avoids text formatting round-trips.';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Pass event type: satellite visibility window
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION pass_event_in(cstring) RETURNS pass_event
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION pass_event_out(pass_event) RETURNS cstring
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION pass_event_recv(internal) RETURNS pass_event
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION pass_event_send(pass_event) RETURNS bytea
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE TYPE pass_event (
|
||||
INPUT = pass_event_in,
|
||||
OUTPUT = pass_event_out,
|
||||
RECEIVE = pass_event_recv,
|
||||
SEND = pass_event_send,
|
||||
INTERNALLENGTH = 48,
|
||||
ALIGNMENT = double,
|
||||
STORAGE = plain
|
||||
);
|
||||
|
||||
COMMENT ON TYPE pass_event IS 'Satellite pass event (AOS/MAX/LOS times, max elevation, AOS/LOS azimuths)';
|
||||
|
||||
CREATE FUNCTION pass_aos_time(pass_event) RETURNS timestamptz
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION pass_aos_time(pass_event) IS 'Acquisition of signal time';
|
||||
|
||||
CREATE FUNCTION pass_max_el_time(pass_event) RETURNS timestamptz
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION pass_max_el_time(pass_event) IS 'Maximum elevation time';
|
||||
|
||||
CREATE FUNCTION pass_los_time(pass_event) RETURNS timestamptz
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION pass_los_time(pass_event) IS 'Loss of signal time';
|
||||
|
||||
CREATE FUNCTION pass_max_elevation(pass_event) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION pass_max_elevation(pass_event) IS 'Maximum elevation in degrees';
|
||||
|
||||
CREATE FUNCTION pass_aos_azimuth(pass_event) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION pass_aos_azimuth(pass_event) IS 'AOS azimuth in degrees (0=N)';
|
||||
|
||||
CREATE FUNCTION pass_los_azimuth(pass_event) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION pass_los_azimuth(pass_event) IS 'LOS azimuth in degrees (0=N)';
|
||||
|
||||
CREATE FUNCTION pass_duration(pass_event) RETURNS interval
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION pass_duration(pass_event) IS 'Pass duration (LOS - AOS)';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- SGP4/SDP4 propagation functions
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION sgp4_propagate(tle, timestamptz) RETURNS eci_position
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION sgp4_propagate(tle, timestamptz) IS
|
||||
'Propagate TLE to a point in time using SGP4 (near-earth) or SDP4 (deep-space). Returns TEME ECI position/velocity.';
|
||||
|
||||
CREATE FUNCTION sgp4_propagate_safe(tle, timestamptz) RETURNS eci_position
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION sgp4_propagate_safe(tle, timestamptz) IS
|
||||
'Like sgp4_propagate but returns NULL on error instead of raising an exception. For batch queries with potentially invalid TLEs.';
|
||||
|
||||
CREATE FUNCTION sgp4_propagate_series(tle, timestamptz, timestamptz, interval)
|
||||
RETURNS TABLE(t timestamptz, x float8, y float8, z float8, vx float8, vy float8, vz float8)
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE
|
||||
ROWS 100;
|
||||
COMMENT ON FUNCTION sgp4_propagate_series(tle, timestamptz, timestamptz, interval) IS
|
||||
'Propagate TLE over a time range at regular intervals. Returns time series of TEME positions.';
|
||||
|
||||
CREATE FUNCTION tle_distance(tle, tle, timestamptz) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_distance(tle, tle, timestamptz) IS
|
||||
'Euclidean distance in km between two TLEs at a reference time';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Coordinate transform functions
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION eci_to_geodetic(eci_position, timestamptz) RETURNS geodetic
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION eci_to_geodetic(eci_position, timestamptz) IS
|
||||
'Convert TEME ECI position to WGS-84 geodetic coordinates at given time';
|
||||
|
||||
CREATE FUNCTION eci_to_topocentric(eci_position, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION eci_to_topocentric(eci_position, observer, timestamptz) IS
|
||||
'Convert TEME ECI position to topocentric (az/el/range) relative to observer';
|
||||
|
||||
CREATE FUNCTION subsatellite_point(tle, timestamptz) RETURNS geodetic
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION subsatellite_point(tle, timestamptz) IS
|
||||
'Subsatellite (nadir) point on WGS-84 ellipsoid at given time';
|
||||
|
||||
CREATE FUNCTION ground_track(tle, timestamptz, timestamptz, interval)
|
||||
RETURNS TABLE(t timestamptz, lat float8, lon float8, alt float8)
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE
|
||||
ROWS 100;
|
||||
COMMENT ON FUNCTION ground_track(tle, timestamptz, timestamptz, interval) IS
|
||||
'Ground track as time series of subsatellite points (lat/lon in degrees, alt in km)';
|
||||
|
||||
CREATE FUNCTION observe(tle, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION observe(tle, observer, timestamptz) IS
|
||||
'Propagate TLE and compute observer-relative look angles in one call. Returns topocentric (az/el/range/range_rate).';
|
||||
|
||||
CREATE FUNCTION observe_safe(tle, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION observe_safe(tle, observer, timestamptz) IS
|
||||
'Like observe() but returns NULL on propagation error. For batch queries with potentially invalid/decayed TLEs.';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Pass prediction functions
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION next_pass(tle, observer, timestamptz) RETURNS pass_event
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION next_pass(tle, observer, timestamptz) IS
|
||||
'Find the next satellite pass over observer (searches up to 7 days ahead)';
|
||||
|
||||
CREATE FUNCTION predict_passes(tle, observer, timestamptz, timestamptz, float8 DEFAULT 0.0)
|
||||
RETURNS SETOF pass_event
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE
|
||||
ROWS 10;
|
||||
COMMENT ON FUNCTION predict_passes(tle, observer, timestamptz, timestamptz, float8) IS
|
||||
'Predict all satellite passes over observer in time window. Optional min_elevation in degrees.';
|
||||
|
||||
CREATE FUNCTION pass_visible(tle, observer, timestamptz, timestamptz) RETURNS boolean
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION pass_visible(tle, observer, timestamptz, timestamptz) IS
|
||||
'True if any pass occurs over observer in the time window';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- GiST operator support functions
|
||||
-- ============================================================
|
||||
|
||||
-- Overlap operator: do orbital keys overlap in altitude AND inclination?
|
||||
CREATE FUNCTION tle_overlap(tle, tle) RETURNS boolean
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
-- Altitude distance operator (altitude-only, for KNN ordering)
|
||||
CREATE FUNCTION tle_alt_distance(tle, tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE OPERATOR && (
|
||||
LEFTARG = tle,
|
||||
RIGHTARG = tle,
|
||||
FUNCTION = tle_overlap,
|
||||
COMMUTATOR = &&,
|
||||
RESTRICT = areasel,
|
||||
JOIN = areajoinsel
|
||||
);
|
||||
|
||||
COMMENT ON OPERATOR && (tle, tle) IS 'Orbital key overlap (altitude band AND inclination range) — necessary condition for conjunction';
|
||||
|
||||
CREATE OPERATOR <-> (
|
||||
LEFTARG = tle,
|
||||
RIGHTARG = tle,
|
||||
FUNCTION = tle_alt_distance,
|
||||
COMMUTATOR = <->
|
||||
);
|
||||
|
||||
COMMENT ON OPERATOR <-> (tle, tle) IS 'Minimum altitude-band separation in km (0 if overlapping). Altitude-only — does not account for inclination. Use && for 2-D filtering.';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- GiST operator class for 2-D orbital indexing (altitude + inclination)
|
||||
-- ============================================================
|
||||
|
||||
-- GiST internal support functions
|
||||
CREATE FUNCTION gist_tle_compress(internal) RETURNS internal
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION gist_tle_decompress(internal) RETURNS internal
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION gist_tle_consistent(internal, tle, smallint, oid, internal) RETURNS boolean
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION gist_tle_union(internal, internal) RETURNS internal
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION gist_tle_penalty(internal, internal, internal) RETURNS internal
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION gist_tle_picksplit(internal, internal) RETURNS internal
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION gist_tle_same(internal, internal, internal) RETURNS internal
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION gist_tle_distance(internal, tle, smallint, oid, internal) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE OPERATOR CLASS tle_ops
|
||||
DEFAULT FOR TYPE tle USING gist AS
|
||||
OPERATOR 3 && ,
|
||||
OPERATOR 15 <-> (tle, tle) FOR ORDER BY float_ops,
|
||||
FUNCTION 1 gist_tle_consistent(internal, tle, smallint, oid, internal),
|
||||
FUNCTION 2 gist_tle_union(internal, internal),
|
||||
FUNCTION 3 gist_tle_compress(internal),
|
||||
FUNCTION 4 gist_tle_decompress(internal),
|
||||
FUNCTION 5 gist_tle_penalty(internal, internal, internal),
|
||||
FUNCTION 6 gist_tle_picksplit(internal, internal),
|
||||
FUNCTION 7 gist_tle_same(internal, internal, internal),
|
||||
FUNCTION 8 gist_tle_distance(internal, tle, smallint, oid, internal);
|
||||
-- pg_orbit 0.1.0 -> 0.2.0 migration
|
||||
--
|
||||
-- Phase 1: Stars, comets, and Keplerian propagation.
|
||||
-- Adds heliocentric type, star observation, and two-body propagation.
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Heliocentric type: ecliptic J2000 position in AU
|
||||
-- ============================================================
|
||||
|
||||
CREATE TYPE heliocentric;
|
||||
|
||||
CREATE FUNCTION heliocentric_in(cstring) RETURNS heliocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION heliocentric_out(heliocentric) RETURNS cstring
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION heliocentric_recv(internal) RETURNS heliocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION heliocentric_send(heliocentric) RETURNS bytea
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE TYPE heliocentric (
|
||||
INPUT = heliocentric_in,
|
||||
OUTPUT = heliocentric_out,
|
||||
RECEIVE = heliocentric_recv,
|
||||
SEND = heliocentric_send,
|
||||
INTERNALLENGTH = 24,
|
||||
ALIGNMENT = double,
|
||||
STORAGE = plain
|
||||
);
|
||||
|
||||
COMMENT ON TYPE heliocentric IS 'Heliocentric position in ecliptic J2000 frame (x, y, z in AU)';
|
||||
|
||||
CREATE FUNCTION helio_x(heliocentric) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION helio_x(heliocentric) IS 'X component in AU (ecliptic J2000, vernal equinox direction)';
|
||||
|
||||
CREATE FUNCTION helio_y(heliocentric) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION helio_y(heliocentric) IS 'Y component in AU (ecliptic J2000)';
|
||||
|
||||
CREATE FUNCTION helio_z(heliocentric) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION helio_z(heliocentric) IS 'Z component in AU (ecliptic J2000, north ecliptic pole)';
|
||||
|
||||
CREATE FUNCTION helio_distance(heliocentric) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION helio_distance(heliocentric) IS 'Heliocentric distance in AU';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Star observation functions
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION star_observe(float8, float8, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION star_observe(float8, float8, observer, timestamptz) IS
|
||||
'Observe a star from (ra_hours J2000, dec_degrees J2000, observer, time). Returns topocentric az/el. Range is 0 (infinite distance).';
|
||||
|
||||
CREATE FUNCTION star_observe_safe(float8, float8, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION star_observe_safe(float8, float8, observer, timestamptz) IS
|
||||
'Like star_observe but returns NULL on invalid inputs. For batch queries over star catalogs.';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Keplerian propagation functions
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION kepler_propagate(
|
||||
q_au float8, eccentricity float8,
|
||||
inclination_deg float8, arg_perihelion_deg float8,
|
||||
long_asc_node_deg float8, perihelion_jd float8,
|
||||
t timestamptz
|
||||
) RETURNS heliocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION kepler_propagate(float8, float8, float8, float8, float8, float8, timestamptz) IS
|
||||
'Two-body Keplerian propagation from classical orbital elements. Returns heliocentric ecliptic J2000 position in AU. Handles elliptic, parabolic, and hyperbolic orbits.';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Comet observation
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION comet_observe(
|
||||
q_au float8, eccentricity float8,
|
||||
inclination_deg float8, arg_perihelion_deg float8,
|
||||
long_asc_node_deg float8, perihelion_jd float8,
|
||||
earth_x_au float8, earth_y_au float8, earth_z_au float8,
|
||||
obs observer, t timestamptz
|
||||
) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION comet_observe(float8, float8, float8, float8, float8, float8, float8, float8, float8, observer, timestamptz) IS
|
||||
'Observe a comet/asteroid from orbital elements. Requires Earth heliocentric ecliptic J2000 position (AU). Returns topocentric az/el with geocentric range in km.';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Phase 2: VSOP87 planets, ELP82B Moon, Sun observation
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION planet_heliocentric(int4, timestamptz) RETURNS heliocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION planet_heliocentric(int4, timestamptz) IS
|
||||
'VSOP87 heliocentric ecliptic J2000 position (AU). Body IDs: 0=Sun, 1=Mercury, ..., 8=Neptune.';
|
||||
|
||||
CREATE FUNCTION planet_observe(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION planet_observe(int4, observer, timestamptz) IS
|
||||
'Observe a planet from (body_id 1-8, observer, time). Returns topocentric az/el with geocentric range in km.';
|
||||
|
||||
CREATE FUNCTION sun_observe(observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION sun_observe(observer, timestamptz) IS
|
||||
'Observe the Sun from (observer, time). Returns topocentric az/el with geocentric range in km.';
|
||||
|
||||
CREATE FUNCTION moon_observe(observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION moon_observe(observer, timestamptz) IS
|
||||
'Observe the Moon via ELP2000-82B from (observer, time). Returns topocentric az/el with geocentric range in km.';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Phase 3: Planetary moon observation
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION galilean_observe(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION galilean_observe(int4, observer, timestamptz) IS
|
||||
'Observe a Galilean moon of Jupiter via L1.2 theory. Body IDs: 0=Io, 1=Europa, 2=Ganymede, 3=Callisto.';
|
||||
|
||||
CREATE FUNCTION saturn_moon_observe(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION saturn_moon_observe(int4, observer, timestamptz) IS
|
||||
'Observe a Saturn moon via TASS 1.7. Body IDs: 0=Mimas, 1=Enceladus, 2=Tethys, 3=Dione, 4=Rhea, 5=Titan, 6=Iapetus, 7=Hyperion.';
|
||||
|
||||
CREATE FUNCTION uranus_moon_observe(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION uranus_moon_observe(int4, observer, timestamptz) IS
|
||||
'Observe a Uranus moon via GUST86. Body IDs: 0=Miranda, 1=Ariel, 2=Umbriel, 3=Titania, 4=Oberon.';
|
||||
|
||||
CREATE FUNCTION mars_moon_observe(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION mars_moon_observe(int4, observer, timestamptz) IS
|
||||
'Observe a Mars moon via MarsSat. Body IDs: 0=Phobos, 1=Deimos.';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Phase 3: Jupiter decametric radio burst prediction
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION io_phase_angle(timestamptz) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION io_phase_angle(timestamptz) IS
|
||||
'Io orbital phase angle in degrees [0,360). 0=superior conjunction (behind Jupiter). Standard Radio JOVE convention.';
|
||||
|
||||
CREATE FUNCTION jupiter_cml(observer, timestamptz) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION jupiter_cml(observer, timestamptz) IS
|
||||
'Jupiter Central Meridian Longitude, System III (1965.0), in degrees [0,360). Light-time corrected.';
|
||||
|
||||
CREATE FUNCTION jupiter_burst_probability(float8, float8) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION jupiter_burst_probability(float8, float8) IS
|
||||
'Estimated Jupiter decametric burst probability (0-1) from (io_phase_deg, cml_deg). Based on Carr et al. (1983) source regions A, B, C, D.';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Phase 4: Interplanetary transfer orbits (Lambert solver)
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION lambert_transfer(
|
||||
dep_body_id int4, arr_body_id int4,
|
||||
dep_time timestamptz, arr_time timestamptz,
|
||||
OUT c3_departure float8, OUT c3_arrival float8,
|
||||
OUT v_inf_departure float8, OUT v_inf_arrival float8,
|
||||
OUT tof_days float8, OUT transfer_sma float8
|
||||
) RETURNS RECORD
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lambert_transfer(int4, int4, timestamptz, timestamptz) IS
|
||||
'Solve Lambert transfer between two planets. Returns C3 (km^2/s^2), v_infinity (km/s), TOF (days), SMA (AU). Body IDs 1-8.';
|
||||
|
||||
CREATE FUNCTION lambert_c3(int4, int4, timestamptz, timestamptz) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lambert_c3(int4, int4, timestamptz, timestamptz) IS
|
||||
'Departure C3 (km^2/s^2) for a Lambert transfer. Returns NULL if solver fails. For pork chop plots.';
|
||||
-- pg_orbit 0.2.0 -> 0.3.0 migration
|
||||
--
|
||||
-- Adds optional JPL DE440/441 ephemeris functions.
|
||||
-- Existing VSOP87/ELP2000-82B functions are unchanged (still IMMUTABLE).
|
||||
-- New _de() functions are STABLE (depend on external DE binary file).
|
||||
-- When DE is unavailable, _de() functions fall back to VSOP87 silently.
|
||||
|
||||
-- ============================================================
|
||||
-- Phase 5: DE ephemeris functions (optional high-precision)
|
||||
-- ============================================================
|
||||
|
||||
-- Planet observation with DE ephemeris
|
||||
|
||||
CREATE FUNCTION planet_heliocentric_de(int4, timestamptz) RETURNS heliocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION planet_heliocentric_de(int4, timestamptz) IS
|
||||
'Heliocentric ecliptic J2000 position via JPL DE (sub-arcsecond). Falls back to VSOP87 if DE unavailable. Body IDs: 0=Sun, 1-8=Mercury-Neptune.';
|
||||
|
||||
CREATE FUNCTION planet_observe_de(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION planet_observe_de(int4, observer, timestamptz) IS
|
||||
'Observe planet via JPL DE. Falls back to VSOP87. Body IDs: 1-8 (Mercury-Neptune).';
|
||||
|
||||
CREATE FUNCTION sun_observe_de(observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION sun_observe_de(observer, timestamptz) IS
|
||||
'Observe Sun via JPL DE. Falls back to VSOP87.';
|
||||
|
||||
CREATE FUNCTION moon_observe_de(observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION moon_observe_de(observer, timestamptz) IS
|
||||
'Observe Moon via JPL DE. Falls back to ELP2000-82B.';
|
||||
|
||||
|
||||
-- Lambert transfer with DE positions
|
||||
|
||||
CREATE FUNCTION lambert_transfer_de(
|
||||
dep_body_id int4, arr_body_id int4,
|
||||
dep_time timestamptz, arr_time timestamptz,
|
||||
OUT c3_departure float8, OUT c3_arrival float8,
|
||||
OUT v_inf_departure float8, OUT v_inf_arrival float8,
|
||||
OUT tof_days float8, OUT transfer_sma float8
|
||||
) RETURNS RECORD
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lambert_transfer_de(int4, int4, timestamptz, timestamptz) IS
|
||||
'Lambert transfer via JPL DE positions. Falls back to VSOP87. Body IDs 1-8.';
|
||||
|
||||
CREATE FUNCTION lambert_c3_de(int4, int4, timestamptz, timestamptz) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lambert_c3_de(int4, int4, timestamptz, timestamptz) IS
|
||||
'Departure C3 via JPL DE. Falls back to VSOP87. For pork chop plots.';
|
||||
|
||||
|
||||
-- Planetary moon observation with DE parent positions
|
||||
|
||||
CREATE FUNCTION galilean_observe_de(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION galilean_observe_de(int4, observer, timestamptz) IS
|
||||
'Observe Galilean moon with JPL DE parent position. L1.2 moon offsets. Body IDs: 0-3 (Io-Callisto).';
|
||||
|
||||
CREATE FUNCTION saturn_moon_observe_de(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION saturn_moon_observe_de(int4, observer, timestamptz) IS
|
||||
'Observe Saturn moon with JPL DE parent position. TASS 1.7 moon offsets. Body IDs: 0-7 (Mimas-Hyperion).';
|
||||
|
||||
CREATE FUNCTION uranus_moon_observe_de(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION uranus_moon_observe_de(int4, observer, timestamptz) IS
|
||||
'Observe Uranus moon with JPL DE parent position. GUST86 moon offsets. Body IDs: 0-4 (Miranda-Oberon).';
|
||||
|
||||
CREATE FUNCTION mars_moon_observe_de(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION mars_moon_observe_de(int4, observer, timestamptz) IS
|
||||
'Observe Mars moon with JPL DE parent position. MarsSat moon offsets. Body IDs: 0-1 (Phobos-Deimos).';
|
||||
|
||||
|
||||
-- Diagnostic function
|
||||
|
||||
CREATE FUNCTION pg_orbit_ephemeris_info(
|
||||
OUT provider text, OUT file_path text,
|
||||
OUT start_jd float8, OUT end_jd float8,
|
||||
OUT version int4, OUT au_km float8
|
||||
) RETURNS RECORD
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION pg_orbit_ephemeris_info() IS
|
||||
'Returns current ephemeris provider status: VSOP87 or JPL_DE with file path, JD range, version, and AU value.';
|
||||
@ -167,4 +167,47 @@ cartesian_to_spherical(const double *xyz, double *ra, double *dec, double *dist)
|
||||
*ra += 2.0 * M_PI;
|
||||
}
|
||||
|
||||
/*
|
||||
* Geocentric observation pipeline (shared by all observation functions).
|
||||
*
|
||||
* Takes geocentric ecliptic J2000 position in AU, observer location,
|
||||
* and Julian date. Converts through equatorial, precesses to date,
|
||||
* and computes topocentric az/el.
|
||||
*
|
||||
* This is the canonical path:
|
||||
* ecliptic J2000 -> equatorial J2000 -> precess to date ->
|
||||
* sidereal time -> hour angle -> az/el
|
||||
*/
|
||||
static inline void
|
||||
observe_from_geocentric(const double geo_ecl_au[3], double jd,
|
||||
const pg_observer *obs, pg_topocentric *result)
|
||||
{
|
||||
double geo_equ[3];
|
||||
double ra_j2000, dec_j2000, geo_dist;
|
||||
double ra_date, dec_date;
|
||||
double gmst_val, lst, ha;
|
||||
double az, el;
|
||||
|
||||
/* Ecliptic J2000 -> equatorial J2000 */
|
||||
ecliptic_to_equatorial(geo_ecl_au, geo_equ);
|
||||
|
||||
/* Cartesian -> spherical */
|
||||
cartesian_to_spherical(geo_equ, &ra_j2000, &dec_j2000, &geo_dist);
|
||||
|
||||
/* Precess J2000 -> date */
|
||||
precess_j2000_to_date(jd, ra_j2000, dec_j2000, &ra_date, &dec_date);
|
||||
|
||||
/* Hour angle and az/el */
|
||||
gmst_val = gmst_from_jd(jd);
|
||||
lst = gmst_val + obs->lon;
|
||||
ha = lst - ra_date;
|
||||
|
||||
equatorial_to_horizontal(ha, dec_date, obs->lat, &az, &el);
|
||||
|
||||
result->azimuth = az;
|
||||
result->elevation = el;
|
||||
result->range_km = geo_dist * AU_KM;
|
||||
result->range_rate = 0.0; /* no velocity computation yet */
|
||||
}
|
||||
|
||||
#endif /* PG_ORBIT_ASTRO_MATH_H */
|
||||
|
||||
651
src/de_funcs.c
Normal file
651
src/de_funcs.c
Normal file
@ -0,0 +1,651 @@
|
||||
/*
|
||||
* de_funcs.c -- SQL-facing DE ephemeris function variants
|
||||
*
|
||||
* Each _de() function is a STABLE STRICT PARALLEL SAFE variant of an
|
||||
* existing IMMUTABLE function. On any DE failure, falls back to the
|
||||
* compiled-in VSOP87/ELP2000-82B equivalent with a NOTICE.
|
||||
*
|
||||
* The observation pipeline is identical:
|
||||
* 1. Heliocentric ecliptic J2000 position (DE or fallback)
|
||||
* 2. Geocentric ecliptic (subtract Earth's heliocentric)
|
||||
* 3. observe_from_geocentric() -> topocentric az/el
|
||||
*
|
||||
* Constant chain of custody rule 7:
|
||||
* Both target and Earth ALWAYS come from the same provider.
|
||||
* If DE fails for the target, we don't use DE for Earth either.
|
||||
*/
|
||||
|
||||
#include "postgres.h"
|
||||
#include "fmgr.h"
|
||||
#include "funcapi.h"
|
||||
#include "catalog/pg_type.h"
|
||||
#include "utils/builtins.h"
|
||||
#include "utils/timestamp.h"
|
||||
|
||||
#include "types.h"
|
||||
#include "astro_math.h"
|
||||
#include "eph_provider.h"
|
||||
#include "vsop87.h"
|
||||
#include "elp82b.h"
|
||||
#include "lambert.h"
|
||||
#include "l12.h"
|
||||
#include "tass17.h"
|
||||
#include "gust86.h"
|
||||
#include "marssat.h"
|
||||
|
||||
#include <math.h>
|
||||
#include <string.h>
|
||||
|
||||
/* Forward declarations */
|
||||
PG_FUNCTION_INFO_V1(planet_heliocentric_de);
|
||||
PG_FUNCTION_INFO_V1(planet_observe_de);
|
||||
PG_FUNCTION_INFO_V1(sun_observe_de);
|
||||
PG_FUNCTION_INFO_V1(moon_observe_de);
|
||||
PG_FUNCTION_INFO_V1(lambert_transfer_de);
|
||||
PG_FUNCTION_INFO_V1(lambert_c3_de);
|
||||
PG_FUNCTION_INFO_V1(galilean_observe_de);
|
||||
PG_FUNCTION_INFO_V1(saturn_moon_observe_de);
|
||||
PG_FUNCTION_INFO_V1(uranus_moon_observe_de);
|
||||
PG_FUNCTION_INFO_V1(mars_moon_observe_de);
|
||||
PG_FUNCTION_INFO_V1(pg_orbit_ephemeris_info);
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* planet_heliocentric_de(body_id int, timestamptz) -> heliocentric
|
||||
*
|
||||
* DE variant of planet_heliocentric(). STABLE.
|
||||
* Falls back to VSOP87 if DE is unavailable.
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
planet_heliocentric_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
int64 ts = PG_GETARG_INT64(1);
|
||||
double jd;
|
||||
double xyz[6];
|
||||
pg_heliocentric *result;
|
||||
|
||||
if (body_id == BODY_SUN)
|
||||
{
|
||||
result = (pg_heliocentric *) palloc(sizeof(pg_heliocentric));
|
||||
result->x = 0.0;
|
||||
result->y = 0.0;
|
||||
result->z = 0.0;
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("invalid body_id %d: must be 0 (Sun) or 1-8 (Mercury-Neptune)",
|
||||
body_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
|
||||
/* Try DE first */
|
||||
if (!eph_de_planet(body_id, jd, xyz))
|
||||
{
|
||||
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, xyz);
|
||||
}
|
||||
|
||||
result = (pg_heliocentric *) palloc(sizeof(pg_heliocentric));
|
||||
result->x = xyz[0];
|
||||
result->y = xyz[1];
|
||||
result->z = xyz[2];
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* planet_observe_de(body_id int, observer, timestamptz) -> topocentric
|
||||
*
|
||||
* DE variant of planet_observe(). STABLE.
|
||||
* Rule 7: both planet and Earth from the same provider.
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
planet_observe_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(1);
|
||||
int64 ts = PG_GETARG_INT64(2);
|
||||
double jd;
|
||||
double earth_xyz[6];
|
||||
double planet_xyz[6];
|
||||
double geo_ecl[3];
|
||||
pg_topocentric *result;
|
||||
|
||||
if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("planet_observe_de: body_id %d must be 1-8 (Mercury-Neptune)",
|
||||
body_id)));
|
||||
|
||||
if (body_id == BODY_EARTH)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
|
||||
errmsg("cannot observe Earth from Earth")));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
|
||||
/* Try DE for both planet and Earth (rule 7: same provider) */
|
||||
if (eph_de_planet(body_id, jd, planet_xyz) &&
|
||||
eph_de_earth(jd, earth_xyz))
|
||||
{
|
||||
/* DE succeeded */
|
||||
}
|
||||
else
|
||||
{
|
||||
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, 2, earth_xyz);
|
||||
GetVsop87Coor(jd, vsop_body, planet_xyz);
|
||||
}
|
||||
|
||||
/* Geocentric ecliptic = planet - Earth */
|
||||
geo_ecl[0] = planet_xyz[0] - earth_xyz[0];
|
||||
geo_ecl[1] = planet_xyz[1] - earth_xyz[1];
|
||||
geo_ecl[2] = planet_xyz[2] - earth_xyz[2];
|
||||
|
||||
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
|
||||
observe_from_geocentric(geo_ecl, jd, obs, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* sun_observe_de(observer, timestamptz) -> topocentric
|
||||
*
|
||||
* DE variant of sun_observe(). STABLE.
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
sun_observe_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
|
||||
int64 ts = PG_GETARG_INT64(1);
|
||||
double jd;
|
||||
double earth_xyz[6];
|
||||
double geo_ecl[3];
|
||||
pg_topocentric *result;
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
|
||||
/* Sun geocentric = -Earth_heliocentric */
|
||||
if (eph_de_earth(jd, earth_xyz))
|
||||
{
|
||||
geo_ecl[0] = -earth_xyz[0];
|
||||
geo_ecl[1] = -earth_xyz[1];
|
||||
geo_ecl[2] = -earth_xyz[2];
|
||||
}
|
||||
else
|
||||
{
|
||||
if (eph_de_available())
|
||||
ereport(NOTICE,
|
||||
(errmsg("DE ephemeris unavailable, falling back to VSOP87")));
|
||||
|
||||
GetVsop87Coor(jd, 2, earth_xyz);
|
||||
geo_ecl[0] = -earth_xyz[0];
|
||||
geo_ecl[1] = -earth_xyz[1];
|
||||
geo_ecl[2] = -earth_xyz[2];
|
||||
}
|
||||
|
||||
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
|
||||
observe_from_geocentric(geo_ecl, jd, obs, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* moon_observe_de(observer, timestamptz) -> topocentric
|
||||
*
|
||||
* DE variant of moon_observe(). STABLE.
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
moon_observe_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
|
||||
int64 ts = PG_GETARG_INT64(1);
|
||||
double jd;
|
||||
double moon_ecl[3];
|
||||
pg_topocentric *result;
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
|
||||
/* Moon geocentric ecliptic J2000 */
|
||||
if (!eph_de_moon(jd, moon_ecl))
|
||||
{
|
||||
if (eph_de_available())
|
||||
ereport(NOTICE,
|
||||
(errmsg("DE ephemeris unavailable, falling back to ELP2000-82B")));
|
||||
|
||||
GetElp82bCoor(jd, moon_ecl);
|
||||
}
|
||||
|
||||
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
|
||||
observe_from_geocentric(moon_ecl, jd, obs, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* Lambert transfer functions with DE positions
|
||||
* ================================================================
|
||||
*/
|
||||
|
||||
/*
|
||||
* Compute planet heliocentric velocity via numerical differentiation.
|
||||
* Uses DE positions if available, VSOP87 otherwise.
|
||||
*/
|
||||
static void
|
||||
planet_velocity_de(int body_id, double jd, double vel[3])
|
||||
{
|
||||
double pos_fwd[6], pos_bwd[6];
|
||||
double dt = 0.01; /* days */
|
||||
bool got_fwd, got_bwd;
|
||||
|
||||
got_fwd = eph_de_planet(body_id, jd + dt, pos_fwd);
|
||||
got_bwd = eph_de_planet(body_id, jd - dt, pos_bwd);
|
||||
|
||||
if (!got_fwd || !got_bwd)
|
||||
{
|
||||
/* Fall back to VSOP87 */
|
||||
int vsop_body = body_id - 1;
|
||||
GetVsop87Coor(jd + dt, vsop_body, pos_fwd);
|
||||
GetVsop87Coor(jd - dt, vsop_body, pos_bwd);
|
||||
}
|
||||
|
||||
vel[0] = (pos_fwd[0] - pos_bwd[0]) / (2.0 * dt);
|
||||
vel[1] = (pos_fwd[1] - pos_bwd[1]) / (2.0 * dt);
|
||||
vel[2] = (pos_fwd[2] - pos_bwd[2]) / (2.0 * dt);
|
||||
}
|
||||
|
||||
|
||||
Datum
|
||||
lambert_transfer_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 dep_body = PG_GETARG_INT32(0);
|
||||
int32 arr_body = PG_GETARG_INT32(1);
|
||||
int64 dep_ts = PG_GETARG_INT64(2);
|
||||
int64 arr_ts = PG_GETARG_INT64(3);
|
||||
double dep_jd, arr_jd, tof_days;
|
||||
double r1[6], r2[6];
|
||||
double v_planet_dep[3], v_planet_arr[3];
|
||||
double v_inf_dep[3], v_inf_arr[3];
|
||||
double v_inf_dep_mag, v_inf_arr_mag;
|
||||
double c3_dep, c3_arr;
|
||||
lambert_result lr;
|
||||
TupleDesc tupdesc;
|
||||
Datum values[6];
|
||||
bool nulls[6];
|
||||
HeapTuple tuple;
|
||||
double au_per_day_to_km_per_s;
|
||||
int k;
|
||||
bool have_de;
|
||||
|
||||
if (dep_body < 1 || dep_body > 8)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
|
||||
errmsg("lambert_transfer_de: dep_body_id %d must be 1-8", dep_body)));
|
||||
if (arr_body < 1 || arr_body > 8)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
|
||||
errmsg("lambert_transfer_de: arr_body_id %d must be 1-8", arr_body)));
|
||||
if (dep_body == arr_body)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
|
||||
errmsg("lambert_transfer_de: departure and arrival bodies must differ")));
|
||||
|
||||
dep_jd = timestamptz_to_jd(dep_ts);
|
||||
arr_jd = timestamptz_to_jd(arr_ts);
|
||||
tof_days = arr_jd - dep_jd;
|
||||
|
||||
if (tof_days <= 0.0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
|
||||
errmsg("lambert_transfer_de: arrival must be after departure")));
|
||||
|
||||
/* Try DE for both positions (rule 7: same provider) */
|
||||
have_de = eph_de_planet(dep_body, dep_jd, r1) &&
|
||||
eph_de_planet(arr_body, arr_jd, r2);
|
||||
|
||||
if (!have_de)
|
||||
{
|
||||
int dep_vsop = dep_body - 1;
|
||||
int arr_vsop = arr_body - 1;
|
||||
|
||||
if (eph_de_available())
|
||||
ereport(NOTICE,
|
||||
(errmsg("DE ephemeris unavailable, falling back to VSOP87")));
|
||||
|
||||
GetVsop87Coor(dep_jd, dep_vsop, r1);
|
||||
GetVsop87Coor(arr_jd, arr_vsop, r2);
|
||||
}
|
||||
|
||||
if (!lambert_solve_uv(r1, r2, tof_days, GAUSS_K2, 1, &lr))
|
||||
PG_RETURN_NULL();
|
||||
|
||||
/* Planet velocities */
|
||||
planet_velocity_de(dep_body, dep_jd, v_planet_dep);
|
||||
planet_velocity_de(arr_body, arr_jd, v_planet_arr);
|
||||
|
||||
au_per_day_to_km_per_s = AU_KM / 86400.0;
|
||||
|
||||
for (k = 0; k < 3; k++) {
|
||||
v_inf_dep[k] = (lr.v1[k] - v_planet_dep[k]) * au_per_day_to_km_per_s;
|
||||
v_inf_arr[k] = (lr.v2[k] - v_planet_arr[k]) * au_per_day_to_km_per_s;
|
||||
}
|
||||
|
||||
v_inf_dep_mag = sqrt(v_inf_dep[0]*v_inf_dep[0] +
|
||||
v_inf_dep[1]*v_inf_dep[1] +
|
||||
v_inf_dep[2]*v_inf_dep[2]);
|
||||
v_inf_arr_mag = sqrt(v_inf_arr[0]*v_inf_arr[0] +
|
||||
v_inf_arr[1]*v_inf_arr[1] +
|
||||
v_inf_arr[2]*v_inf_arr[2]);
|
||||
|
||||
c3_dep = v_inf_dep_mag * v_inf_dep_mag;
|
||||
c3_arr = v_inf_arr_mag * v_inf_arr_mag;
|
||||
|
||||
if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
|
||||
errmsg("function returning record called in context "
|
||||
"that cannot accept type record")));
|
||||
|
||||
tupdesc = BlessTupleDesc(tupdesc);
|
||||
|
||||
memset(nulls, 0, sizeof(nulls));
|
||||
values[0] = Float8GetDatum(c3_dep);
|
||||
values[1] = Float8GetDatum(c3_arr);
|
||||
values[2] = Float8GetDatum(v_inf_dep_mag);
|
||||
values[3] = Float8GetDatum(v_inf_arr_mag);
|
||||
values[4] = Float8GetDatum(tof_days);
|
||||
values[5] = Float8GetDatum(lr.sma);
|
||||
|
||||
tuple = heap_form_tuple(tupdesc, values, nulls);
|
||||
PG_RETURN_DATUM(HeapTupleGetDatum(tuple));
|
||||
}
|
||||
|
||||
|
||||
Datum
|
||||
lambert_c3_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 dep_body = PG_GETARG_INT32(0);
|
||||
int32 arr_body = PG_GETARG_INT32(1);
|
||||
int64 dep_ts = PG_GETARG_INT64(2);
|
||||
int64 arr_ts = PG_GETARG_INT64(3);
|
||||
double dep_jd, arr_jd, tof_days;
|
||||
double r1[6], r2[6];
|
||||
double v_planet_dep[3];
|
||||
double v_inf_dep[3];
|
||||
double c3_dep;
|
||||
lambert_result lr;
|
||||
double au_per_day_to_km_per_s;
|
||||
int k;
|
||||
bool have_de;
|
||||
|
||||
if (dep_body < 1 || dep_body > 8 || arr_body < 1 || arr_body > 8)
|
||||
PG_RETURN_NULL();
|
||||
if (dep_body == arr_body)
|
||||
PG_RETURN_NULL();
|
||||
|
||||
dep_jd = timestamptz_to_jd(dep_ts);
|
||||
arr_jd = timestamptz_to_jd(arr_ts);
|
||||
tof_days = arr_jd - dep_jd;
|
||||
|
||||
if (tof_days <= 0.0)
|
||||
PG_RETURN_NULL();
|
||||
|
||||
have_de = eph_de_planet(dep_body, dep_jd, r1) &&
|
||||
eph_de_planet(arr_body, arr_jd, r2);
|
||||
|
||||
if (!have_de)
|
||||
{
|
||||
int dep_vsop = dep_body - 1;
|
||||
int arr_vsop = arr_body - 1;
|
||||
|
||||
GetVsop87Coor(dep_jd, dep_vsop, r1);
|
||||
GetVsop87Coor(arr_jd, arr_vsop, r2);
|
||||
}
|
||||
|
||||
if (!lambert_solve_uv(r1, r2, tof_days, GAUSS_K2, 1, &lr))
|
||||
PG_RETURN_NULL();
|
||||
|
||||
planet_velocity_de(dep_body, dep_jd, v_planet_dep);
|
||||
|
||||
au_per_day_to_km_per_s = AU_KM / 86400.0;
|
||||
|
||||
for (k = 0; k < 3; k++)
|
||||
v_inf_dep[k] = (lr.v1[k] - v_planet_dep[k]) * au_per_day_to_km_per_s;
|
||||
|
||||
c3_dep = v_inf_dep[0]*v_inf_dep[0] +
|
||||
v_inf_dep[1]*v_inf_dep[1] +
|
||||
v_inf_dep[2]*v_inf_dep[2];
|
||||
|
||||
PG_RETURN_FLOAT8(c3_dep);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* Planetary moon observation with DE parent positions
|
||||
*
|
||||
* For each planetary moon, the moon-theory offset (L1.2, TASS17,
|
||||
* GUST86, MarsSat) is computed relative to its parent planet.
|
||||
* The parent's position comes from DE instead of VSOP87, giving
|
||||
* sub-arcsecond accuracy for the parent while keeping the
|
||||
* moon-theory accuracy for the relative offset.
|
||||
* ================================================================
|
||||
*/
|
||||
|
||||
/*
|
||||
* Internal: observe a planetary moon using DE for the parent planet
|
||||
* and Earth positions. Falls back to VSOP87 if DE is unavailable.
|
||||
*
|
||||
* moon_rel[3]: moon position relative to parent (ecliptic J2000, AU)
|
||||
* parent_body_id: pg_orbit body ID of parent (5=Jupiter, 6=Saturn, etc.)
|
||||
*/
|
||||
static void
|
||||
observe_moon_de(const double moon_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); /* VSOP87 body 2 = Earth */
|
||||
}
|
||||
|
||||
/* Moon geocentric = (parent + moon_relative) - Earth */
|
||||
geo_ecl[0] = (parent_xyz[0] + moon_rel[0]) - earth_xyz[0];
|
||||
geo_ecl[1] = (parent_xyz[1] + moon_rel[1]) - earth_xyz[1];
|
||||
geo_ecl[2] = (parent_xyz[2] + moon_rel[2]) - earth_xyz[2];
|
||||
|
||||
observe_from_geocentric(geo_ecl, jd, obs, result);
|
||||
}
|
||||
|
||||
|
||||
Datum
|
||||
galilean_observe_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(1);
|
||||
int64 ts = PG_GETARG_INT64(2);
|
||||
double jd;
|
||||
double moon_xyz[3];
|
||||
pg_topocentric *result;
|
||||
|
||||
if (body_id < L12_IO || body_id > L12_CALLISTO)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
|
||||
errmsg("galilean_observe_de: body_id %d must be 0-3 (Io-Callisto)",
|
||||
body_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetL12Coor(jd, body_id, moon_xyz, NULL);
|
||||
|
||||
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
|
||||
observe_moon_de(moon_xyz, BODY_JUPITER, jd, obs, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
Datum
|
||||
saturn_moon_observe_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(1);
|
||||
int64 ts = PG_GETARG_INT64(2);
|
||||
double jd;
|
||||
double moon_xyz[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_observe_de: body_id %d must be 0-7 (Mimas-Hyperion)",
|
||||
body_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetTass17Coor(jd, body_id, moon_xyz, NULL);
|
||||
|
||||
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
|
||||
observe_moon_de(moon_xyz, BODY_SATURN, jd, obs, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
Datum
|
||||
uranus_moon_observe_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(1);
|
||||
int64 ts = PG_GETARG_INT64(2);
|
||||
double jd;
|
||||
double moon_xyz[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_observe_de: body_id %d must be 0-4 (Miranda-Oberon)",
|
||||
body_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetGust86Coor(jd, body_id, moon_xyz, NULL);
|
||||
|
||||
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
|
||||
observe_moon_de(moon_xyz, BODY_URANUS, jd, obs, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
Datum
|
||||
mars_moon_observe_de(PG_FUNCTION_ARGS)
|
||||
{
|
||||
int32 body_id = PG_GETARG_INT32(0);
|
||||
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(1);
|
||||
int64 ts = PG_GETARG_INT64(2);
|
||||
double jd;
|
||||
double moon_xyz[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_observe_de: body_id %d must be 0-1 (Phobos-Deimos)",
|
||||
body_id)));
|
||||
|
||||
jd = timestamptz_to_jd(ts);
|
||||
GetMarsSatCoor(jd, body_id, moon_xyz, NULL);
|
||||
|
||||
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
|
||||
observe_moon_de(moon_xyz, BODY_MARS, jd, obs, result);
|
||||
|
||||
PG_RETURN_POINTER(result);
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* pg_orbit_ephemeris_info() -> RECORD
|
||||
*
|
||||
* Diagnostic function: returns current ephemeris provider status.
|
||||
* STABLE (not STRICT — no args), PARALLEL SAFE.
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
pg_orbit_ephemeris_info(PG_FUNCTION_ARGS)
|
||||
{
|
||||
TupleDesc tupdesc;
|
||||
Datum values[6];
|
||||
bool nulls[6];
|
||||
HeapTuple tuple;
|
||||
|
||||
if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
|
||||
errmsg("function returning record called in context "
|
||||
"that cannot accept type record")));
|
||||
|
||||
tupdesc = BlessTupleDesc(tupdesc);
|
||||
memset(nulls, 0, sizeof(nulls));
|
||||
|
||||
if (eph_de_available())
|
||||
{
|
||||
const char *path = eph_get_path();
|
||||
|
||||
values[0] = CStringGetTextDatum("JPL_DE");
|
||||
values[1] = path ? CStringGetTextDatum(path) : CStringGetTextDatum("");
|
||||
values[2] = Float8GetDatum(eph_de_start_jd());
|
||||
values[3] = Float8GetDatum(eph_de_end_jd());
|
||||
values[4] = Int32GetDatum(eph_de_version());
|
||||
values[5] = Float8GetDatum(eph_de_au_km());
|
||||
}
|
||||
else
|
||||
{
|
||||
values[0] = CStringGetTextDatum("VSOP87");
|
||||
|
||||
nulls[1] = true; /* no file path */
|
||||
nulls[2] = true; /* no start_jd */
|
||||
nulls[3] = true; /* no end_jd */
|
||||
nulls[4] = true; /* no version */
|
||||
|
||||
values[5] = Float8GetDatum((double)AU_KM);
|
||||
}
|
||||
|
||||
tuple = heap_form_tuple(tupdesc, values, nulls);
|
||||
PG_RETURN_DATUM(HeapTupleGetDatum(tuple));
|
||||
}
|
||||
662
src/de_reader.c
Normal file
662
src/de_reader.c
Normal file
@ -0,0 +1,662 @@
|
||||
/*
|
||||
* de_reader.c -- Clean-room JPL Development Ephemeris reader
|
||||
*
|
||||
* Reads JPL DE binary format (DE430/DE440/DE441).
|
||||
* No GPL dependency: derived from the public format specification.
|
||||
*
|
||||
* The JPL binary ephemeris consists of fixed-size records:
|
||||
* Record 1: Header (titles, constant names, JD range, layout table)
|
||||
* Record 2: Constant values
|
||||
* Records 3+: Chebyshev coefficients for all bodies
|
||||
*
|
||||
* Each data record covers a time interval (32 days for DE441).
|
||||
* Within a record, each body has Chebyshev coefficients for x,y,z
|
||||
* (or longitude/latitude/distance), possibly split into sub-intervals.
|
||||
*
|
||||
* Evaluation uses the Clenshaw recurrence for Chebyshev polynomials.
|
||||
*
|
||||
* Reference:
|
||||
* JPL IOM 312.N-03-009 (Standish 1998)
|
||||
* "Explanatory Supplement to the Astronomical Almanac" Ch. 8
|
||||
*/
|
||||
|
||||
#include "de_reader.h"
|
||||
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <math.h>
|
||||
#include <unistd.h>
|
||||
#include <fcntl.h>
|
||||
|
||||
/* Known AU value for byte-order detection */
|
||||
#define DE_AU_KNOWN 149597870.700
|
||||
|
||||
/* Header offsets (in doubles, 0-based within record 1) */
|
||||
#define HDR_TITLE_WORDS 84 /* 3 title lines x 84 chars = 252 chars = ~32 doubles */
|
||||
#define HDR_CNAME_WORDS 400 /* constant names */
|
||||
#define HDR_SS_OFFSET (HDR_TITLE_WORDS / sizeof(double)) /* not used directly */
|
||||
|
||||
/*
|
||||
* Byte-swap a double (for big-endian DE files on little-endian hosts).
|
||||
*/
|
||||
static void
|
||||
swap_double(double *val)
|
||||
{
|
||||
unsigned char *p = (unsigned char *)val;
|
||||
unsigned char tmp;
|
||||
|
||||
tmp = p[0]; p[0] = p[7]; p[7] = tmp;
|
||||
tmp = p[1]; p[1] = p[6]; p[6] = tmp;
|
||||
tmp = p[2]; p[2] = p[5]; p[5] = tmp;
|
||||
tmp = p[3]; p[3] = p[4]; p[4] = tmp;
|
||||
}
|
||||
|
||||
/*
|
||||
* Swap an array of doubles in-place.
|
||||
*/
|
||||
static void
|
||||
swap_doubles(double *arr, int count)
|
||||
{
|
||||
int i;
|
||||
for (i = 0; i < count; i++)
|
||||
swap_double(&arr[i]);
|
||||
}
|
||||
|
||||
/*
|
||||
* Read exactly n bytes from fd at the given offset.
|
||||
* Returns 0 on success, -1 on failure.
|
||||
*/
|
||||
static int
|
||||
read_at(int fd, void *buf, size_t n, off_t offset)
|
||||
{
|
||||
ssize_t total = 0;
|
||||
ssize_t got;
|
||||
|
||||
if (lseek(fd, offset, SEEK_SET) == (off_t)-1)
|
||||
return -1;
|
||||
|
||||
while ((size_t)total < n)
|
||||
{
|
||||
got = read(fd, (char *)buf + total, n - total);
|
||||
if (got <= 0)
|
||||
return -1;
|
||||
total += got;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Evaluate a Chebyshev polynomial using the Clenshaw recurrence.
|
||||
*
|
||||
* coeffs: array of Chebyshev coefficients (T0, T1, ..., T_{n-1})
|
||||
* n: number of coefficients
|
||||
* x: normalized argument in [-1, +1]
|
||||
*
|
||||
* Returns: sum_{i=0}^{n-1} coeffs[i] * T_i(x)
|
||||
*/
|
||||
static double
|
||||
chebyshev_eval(const double *coeffs, int n, double x)
|
||||
{
|
||||
double bk1 = 0.0, bk2 = 0.0, bk;
|
||||
int i;
|
||||
|
||||
for (i = n - 1; i >= 1; i--)
|
||||
{
|
||||
bk = 2.0 * x * bk1 - bk2 + coeffs[i];
|
||||
bk2 = bk1;
|
||||
bk1 = bk;
|
||||
}
|
||||
return x * bk1 - bk2 + coeffs[0];
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Interpolate a single component (x, y, or z) for a body at a JD.
|
||||
*
|
||||
* h: reader handle (with record already loaded in record_buf)
|
||||
* body: body group index (0-12)
|
||||
* comp: component (0=x, 1=y, 2=z)
|
||||
* jd: Julian date
|
||||
*
|
||||
* Returns the interpolated value.
|
||||
*/
|
||||
static double
|
||||
interp_component(de_handle *h, int body, int comp, double jd)
|
||||
{
|
||||
de_body_layout *lay = &h->layout[body];
|
||||
double rec_start, sub_length, t_sub, x;
|
||||
int sub_idx, coeff_offset;
|
||||
|
||||
/* Record start JD */
|
||||
rec_start = h->record_buf[0];
|
||||
|
||||
/* Sub-interval length in days */
|
||||
sub_length = h->interval_days / lay->nsub;
|
||||
|
||||
/* Which sub-interval? */
|
||||
t_sub = jd - rec_start;
|
||||
sub_idx = (int)(t_sub / sub_length);
|
||||
if (sub_idx >= lay->nsub)
|
||||
sub_idx = lay->nsub - 1;
|
||||
if (sub_idx < 0)
|
||||
sub_idx = 0;
|
||||
|
||||
/* Normalize to [-1, +1] within the sub-interval */
|
||||
x = 2.0 * (t_sub - sub_idx * sub_length) / sub_length - 1.0;
|
||||
|
||||
/*
|
||||
* Coefficient offset in the record buffer.
|
||||
* Layout offset is 1-based (Fortran convention), convert to 0-based.
|
||||
* Each sub-interval has ncoeff coefficients for each of 3 components.
|
||||
*/
|
||||
coeff_offset = (lay->offset - 1)
|
||||
+ sub_idx * lay->ncoeff * 3
|
||||
+ comp * lay->ncoeff;
|
||||
|
||||
return chebyshev_eval(&h->record_buf[coeff_offset], lay->ncoeff, x);
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Load the record containing the given JD into the handle's buffer.
|
||||
* Returns DE_OK or DE_ERR_*.
|
||||
*/
|
||||
static int
|
||||
load_record(de_handle *h, double jd)
|
||||
{
|
||||
int recno;
|
||||
off_t offset;
|
||||
|
||||
if (jd < h->start_jd || jd > h->end_jd)
|
||||
return DE_ERR_RANGE;
|
||||
|
||||
recno = (int)((jd - h->start_jd) / h->interval_days);
|
||||
|
||||
/* Clamp to last record if jd == end_jd */
|
||||
{
|
||||
int max_rec = (int)((h->end_jd - h->start_jd) / h->interval_days) - 1;
|
||||
if (recno > max_rec)
|
||||
recno = max_rec;
|
||||
}
|
||||
|
||||
/* Already cached? */
|
||||
if (recno == h->cached_recno)
|
||||
return DE_OK;
|
||||
|
||||
/* Seek and read the record (skip 2 header records) */
|
||||
offset = h->data_offset + (off_t)recno * h->record_bytes;
|
||||
|
||||
if (read_at(h->fd, h->record_buf, h->record_bytes, offset) != 0)
|
||||
return DE_ERR_READ;
|
||||
|
||||
if (h->swap_bytes)
|
||||
swap_doubles(h->record_buf, h->ncoeff);
|
||||
|
||||
h->cached_recno = recno;
|
||||
return DE_OK;
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Parse the header and validate the ephemeris file.
|
||||
*/
|
||||
static int
|
||||
parse_header(de_handle *h)
|
||||
{
|
||||
/*
|
||||
* The header spans 2 records of ncoeff doubles each.
|
||||
* But we don't know ncoeff yet — it's at a fixed position
|
||||
* in the first record.
|
||||
*
|
||||
* Header layout (byte offsets for first record):
|
||||
* 0..251: 3 title lines (84 chars each)
|
||||
* 252..2651: 400 constant names (6 chars each)
|
||||
* 2652..2675: SS[3] = {start_jd, end_jd, interval_days}
|
||||
* 2676..2679: NCON (int32)
|
||||
* 2680..2687: AU (double)
|
||||
* 2688..2695: EMRAT (double)
|
||||
* 2696..2851: IPT[12][3] = 12 body groups x 3 ints (offset, ncoeff, nsub)
|
||||
* 2852..2855: DE version number (int32)
|
||||
* 2856..2867: IPT[12][3] for 13th body (librations)
|
||||
*
|
||||
* After IPT parsing, we know ncoeff and can size the record buffer.
|
||||
*
|
||||
* All byte offsets assume the first record starts at byte 0.
|
||||
* The record size in bytes = ncoeff * 8.
|
||||
*/
|
||||
|
||||
unsigned char buf[4096]; /* enough for the header area */
|
||||
double ss[3];
|
||||
double au_val;
|
||||
int32_t ncon;
|
||||
int ipt[13][3];
|
||||
int32_t de_ver;
|
||||
int i, j;
|
||||
int max_offset_needed;
|
||||
|
||||
/* Read the first 4096 bytes of the file */
|
||||
if (read_at(h->fd, buf, 4096, 0) != 0)
|
||||
return DE_ERR_READ;
|
||||
|
||||
/* SS: start_jd, end_jd, interval (at byte 2652) */
|
||||
memcpy(ss, buf + 2652, 3 * sizeof(double));
|
||||
|
||||
/* AU (at byte 2680) */
|
||||
memcpy(&au_val, buf + 2680, sizeof(double));
|
||||
|
||||
/*
|
||||
* Byte-order detection: compare AU against known value.
|
||||
* If it doesn't match, try byte-swapping.
|
||||
*/
|
||||
h->swap_bytes = 0;
|
||||
if (fabs(au_val - DE_AU_KNOWN) > 1.0)
|
||||
{
|
||||
swap_double(&au_val);
|
||||
if (fabs(au_val - DE_AU_KNOWN) > 1.0)
|
||||
return DE_ERR_ENDIAN;
|
||||
h->swap_bytes = 1;
|
||||
swap_doubles(ss, 3);
|
||||
}
|
||||
|
||||
h->start_jd = ss[0];
|
||||
h->end_jd = ss[1];
|
||||
h->interval_days = ss[2];
|
||||
h->au_km = au_val;
|
||||
|
||||
/* NCON at byte 2676 */
|
||||
memcpy(&ncon, buf + 2676, sizeof(int32_t));
|
||||
if (h->swap_bytes)
|
||||
{
|
||||
unsigned char *p = (unsigned char *)&ncon;
|
||||
unsigned char tmp;
|
||||
tmp = p[0]; p[0] = p[3]; p[3] = tmp;
|
||||
tmp = p[1]; p[1] = p[2]; p[2] = tmp;
|
||||
}
|
||||
|
||||
/* EMRAT at byte 2688 */
|
||||
memcpy(&h->emrat, buf + 2688, sizeof(double));
|
||||
if (h->swap_bytes)
|
||||
swap_double(&h->emrat);
|
||||
|
||||
/* IPT: 12 body groups at byte 2696, each 3 x int32 */
|
||||
for (i = 0; i < 12; i++)
|
||||
{
|
||||
int32_t vals[3];
|
||||
memcpy(vals, buf + 2696 + i * 12, 12);
|
||||
if (h->swap_bytes)
|
||||
{
|
||||
for (j = 0; j < 3; j++)
|
||||
{
|
||||
unsigned char *p = (unsigned char *)&vals[j];
|
||||
unsigned char tmp;
|
||||
tmp = p[0]; p[0] = p[3]; p[3] = tmp;
|
||||
tmp = p[1]; p[1] = p[2]; p[2] = tmp;
|
||||
}
|
||||
}
|
||||
ipt[i][0] = vals[0];
|
||||
ipt[i][1] = vals[1];
|
||||
ipt[i][2] = vals[2];
|
||||
}
|
||||
|
||||
/* DE version at byte 2840 */
|
||||
memcpy(&de_ver, buf + 2840, sizeof(int32_t));
|
||||
if (h->swap_bytes)
|
||||
{
|
||||
unsigned char *p = (unsigned char *)&de_ver;
|
||||
unsigned char tmp;
|
||||
tmp = p[0]; p[0] = p[3]; p[3] = tmp;
|
||||
tmp = p[1]; p[1] = p[2]; p[2] = tmp;
|
||||
}
|
||||
h->de_version = de_ver;
|
||||
|
||||
/* IPT[12] (librations) at byte 2844 */
|
||||
{
|
||||
int32_t vals[3];
|
||||
memcpy(vals, buf + 2844, 12);
|
||||
if (h->swap_bytes)
|
||||
{
|
||||
for (j = 0; j < 3; j++)
|
||||
{
|
||||
unsigned char *p = (unsigned char *)&vals[j];
|
||||
unsigned char tmp;
|
||||
tmp = p[0]; p[0] = p[3]; p[3] = tmp;
|
||||
tmp = p[1]; p[1] = p[2]; p[2] = tmp;
|
||||
}
|
||||
}
|
||||
ipt[12][0] = vals[0];
|
||||
ipt[12][1] = vals[1];
|
||||
ipt[12][2] = vals[2];
|
||||
}
|
||||
|
||||
/* Store layout and compute ncoeff */
|
||||
max_offset_needed = 0;
|
||||
for (i = 0; i < DE_NUM_BODIES; i++)
|
||||
{
|
||||
h->layout[i].offset = ipt[i][0];
|
||||
h->layout[i].ncoeff = ipt[i][1];
|
||||
h->layout[i].nsub = ipt[i][2];
|
||||
|
||||
if (ipt[i][0] > 0 && ipt[i][1] > 0 && ipt[i][2] > 0)
|
||||
{
|
||||
int end;
|
||||
int ncomp = (i == DE_NUTATION) ? 2 : 3;
|
||||
end = (ipt[i][0] - 1) + ipt[i][1] * ncomp * ipt[i][2];
|
||||
if (end > max_offset_needed)
|
||||
max_offset_needed = end;
|
||||
}
|
||||
}
|
||||
|
||||
/*
|
||||
* ncoeff: the header record size in doubles.
|
||||
* The record must be large enough to hold all coefficient data
|
||||
* plus the 2-double JD range at the start of each data record.
|
||||
*/
|
||||
h->ncoeff = max_offset_needed;
|
||||
if (h->ncoeff < 2)
|
||||
return DE_ERR_HEADER;
|
||||
|
||||
h->record_bytes = h->ncoeff * sizeof(double);
|
||||
|
||||
/* Data starts after 2 header records */
|
||||
h->data_offset = 2L * h->record_bytes;
|
||||
|
||||
/* Validate: start_jd < end_jd, interval > 0 */
|
||||
if (h->start_jd >= h->end_jd || h->interval_days <= 0.0)
|
||||
return DE_ERR_HEADER;
|
||||
|
||||
return DE_OK;
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Canary validation: evaluate Earth at J2000.0 and check against
|
||||
* known DE441 reference position.
|
||||
*
|
||||
* Expected Earth ICRS position at J2000.0 (JD 2451545.0):
|
||||
* x ~ -0.1771 AU, y ~ 0.8873 AU, z ~ 0.3848 AU
|
||||
*
|
||||
* Tolerance: 0.01 AU (very loose — just catches gross file corruption).
|
||||
*/
|
||||
static int
|
||||
canary_check(de_handle *h)
|
||||
{
|
||||
double pos[3];
|
||||
int err;
|
||||
|
||||
err = de_reader_get_pos(h, 2451545.0, DE_EMB, DE_SUN, pos);
|
||||
if (err != DE_OK)
|
||||
return DE_ERR_CANARY;
|
||||
|
||||
/* Rough check: Earth should be ~1 AU from Sun */
|
||||
{
|
||||
double dist = sqrt(pos[0]*pos[0] + pos[1]*pos[1] + pos[2]*pos[2]);
|
||||
if (dist < 0.9 || dist > 1.1)
|
||||
return DE_ERR_CANARY;
|
||||
}
|
||||
|
||||
return DE_OK;
|
||||
}
|
||||
|
||||
|
||||
de_handle *
|
||||
de_reader_open(const char *path, int *errcode)
|
||||
{
|
||||
de_handle *h;
|
||||
int err;
|
||||
|
||||
h = (de_handle *)calloc(1, sizeof(de_handle));
|
||||
if (!h)
|
||||
{
|
||||
*errcode = DE_ERR_OPEN;
|
||||
return NULL;
|
||||
}
|
||||
|
||||
h->fd = -1;
|
||||
h->cached_recno = -1;
|
||||
h->record_buf = NULL;
|
||||
|
||||
/* Open the ephemeris file (read-only, no O_CLOEXEC dependency) */
|
||||
h->fd = open(path, O_RDONLY);
|
||||
if (h->fd < 0)
|
||||
{
|
||||
*errcode = DE_ERR_OPEN;
|
||||
free(h);
|
||||
return NULL;
|
||||
}
|
||||
|
||||
/* Parse header */
|
||||
err = parse_header(h);
|
||||
if (err != DE_OK)
|
||||
{
|
||||
*errcode = err;
|
||||
close(h->fd);
|
||||
free(h);
|
||||
return NULL;
|
||||
}
|
||||
|
||||
/* Allocate record buffer */
|
||||
h->record_buf = (double *)calloc(h->ncoeff, sizeof(double));
|
||||
if (!h->record_buf)
|
||||
{
|
||||
*errcode = DE_ERR_OPEN;
|
||||
close(h->fd);
|
||||
free(h);
|
||||
return NULL;
|
||||
}
|
||||
|
||||
/* Canary check */
|
||||
err = canary_check(h);
|
||||
if (err != DE_OK)
|
||||
{
|
||||
*errcode = err;
|
||||
free(h->record_buf);
|
||||
close(h->fd);
|
||||
free(h);
|
||||
return NULL;
|
||||
}
|
||||
|
||||
*errcode = DE_OK;
|
||||
return h;
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Get raw body position (3 components) from the ephemeris.
|
||||
* No center subtraction. Returns position in AU (ICRS equatorial).
|
||||
*/
|
||||
static int
|
||||
get_raw_pos(de_handle *h, double jd, int body, double pos[3])
|
||||
{
|
||||
int err;
|
||||
|
||||
if (body < 0 || body > DE_SUN)
|
||||
return DE_ERR_BODY;
|
||||
|
||||
/* Nutation and libration are not position queries */
|
||||
if (body == DE_NUTATION || body == DE_LIBRATION)
|
||||
return DE_ERR_BODY;
|
||||
|
||||
err = load_record(h, jd);
|
||||
if (err != DE_OK)
|
||||
return err;
|
||||
|
||||
pos[0] = interp_component(h, body, 0, jd);
|
||||
pos[1] = interp_component(h, body, 1, jd);
|
||||
pos[2] = interp_component(h, body, 2, jd);
|
||||
|
||||
/* Convert from km to AU */
|
||||
pos[0] /= h->au_km;
|
||||
pos[1] /= h->au_km;
|
||||
pos[2] /= h->au_km;
|
||||
|
||||
return DE_OK;
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Derive Earth position from Earth-Moon Barycenter and Moon.
|
||||
*
|
||||
* Earth = EMB - Moon * (1 / (1 + EMRAT))
|
||||
*
|
||||
* where EMRAT = M_earth / M_moon ~ 81.3.
|
||||
* The Moon position in the DE file is geocentric, so:
|
||||
* Earth = EMB - Moon / (1 + EMRAT)
|
||||
*/
|
||||
static int
|
||||
get_earth_pos(de_handle *h, double jd, double pos[3])
|
||||
{
|
||||
double emb[3], moon[3];
|
||||
double factor;
|
||||
int err;
|
||||
|
||||
err = get_raw_pos(h, jd, DE_EMB, emb);
|
||||
if (err != DE_OK)
|
||||
return err;
|
||||
|
||||
err = get_raw_pos(h, jd, DE_MOON, moon);
|
||||
if (err != DE_OK)
|
||||
return err;
|
||||
|
||||
factor = 1.0 / (1.0 + h->emrat);
|
||||
|
||||
pos[0] = emb[0] - moon[0] * factor;
|
||||
pos[1] = emb[1] - moon[1] * factor;
|
||||
pos[2] = emb[2] - moon[2] * factor;
|
||||
|
||||
return DE_OK;
|
||||
}
|
||||
|
||||
|
||||
int
|
||||
de_reader_get_pos(de_handle *h, double jd, int target, int center,
|
||||
double pos[3])
|
||||
{
|
||||
double tpos[3], cpos[3];
|
||||
int err;
|
||||
|
||||
if (!h)
|
||||
return DE_ERR_OPEN;
|
||||
|
||||
/* Get target position */
|
||||
if (target == DE_EMB)
|
||||
{
|
||||
/* EMB is stored directly */
|
||||
err = get_raw_pos(h, jd, DE_EMB, tpos);
|
||||
}
|
||||
else if (target == DE_MOON && center != -1)
|
||||
{
|
||||
/*
|
||||
* Moon in the file is geocentric. For Moon relative to
|
||||
* something other than Earth, we need Earth + geocentric Moon.
|
||||
*/
|
||||
double earth[3], moon_geo[3];
|
||||
|
||||
err = get_earth_pos(h, jd, earth);
|
||||
if (err != DE_OK)
|
||||
return err;
|
||||
|
||||
err = get_raw_pos(h, jd, DE_MOON, moon_geo);
|
||||
if (err != DE_OK)
|
||||
return err;
|
||||
|
||||
/* Moon barycentric = Earth + geocentric_moon */
|
||||
tpos[0] = earth[0] + moon_geo[0];
|
||||
tpos[1] = earth[1] + moon_geo[1];
|
||||
tpos[2] = earth[2] + moon_geo[2];
|
||||
}
|
||||
else if (target >= DE_MERCURY && target <= DE_PLUTO && target != DE_EMB)
|
||||
{
|
||||
/* Planets are stored directly as SSB-relative */
|
||||
err = get_raw_pos(h, jd, target, tpos);
|
||||
}
|
||||
else if (target == DE_SUN)
|
||||
{
|
||||
err = get_raw_pos(h, jd, DE_SUN, tpos);
|
||||
}
|
||||
else
|
||||
{
|
||||
/* Raw mode for other body types */
|
||||
err = get_raw_pos(h, jd, target, tpos);
|
||||
}
|
||||
|
||||
if (err != DE_OK)
|
||||
return err;
|
||||
|
||||
/* No center subtraction requested */
|
||||
if (center < 0)
|
||||
{
|
||||
pos[0] = tpos[0];
|
||||
pos[1] = tpos[1];
|
||||
pos[2] = tpos[2];
|
||||
return DE_OK;
|
||||
}
|
||||
|
||||
/* Get center position */
|
||||
if (center == DE_SUN)
|
||||
{
|
||||
err = get_raw_pos(h, jd, DE_SUN, cpos);
|
||||
}
|
||||
else if (center == DE_EMB)
|
||||
{
|
||||
err = get_raw_pos(h, jd, DE_EMB, cpos);
|
||||
}
|
||||
else if (center == 99)
|
||||
{
|
||||
/* Special code for "Earth" as center */
|
||||
err = get_earth_pos(h, jd, cpos);
|
||||
}
|
||||
else if (center >= DE_MERCURY && center <= DE_PLUTO)
|
||||
{
|
||||
err = get_raw_pos(h, jd, center, cpos);
|
||||
}
|
||||
else
|
||||
{
|
||||
return DE_ERR_BODY;
|
||||
}
|
||||
|
||||
if (err != DE_OK)
|
||||
return err;
|
||||
|
||||
pos[0] = tpos[0] - cpos[0];
|
||||
pos[1] = tpos[1] - cpos[1];
|
||||
pos[2] = tpos[2] - cpos[2];
|
||||
|
||||
return DE_OK;
|
||||
}
|
||||
|
||||
|
||||
void
|
||||
de_reader_close(de_handle *h)
|
||||
{
|
||||
if (!h)
|
||||
return;
|
||||
|
||||
if (h->fd >= 0)
|
||||
close(h->fd);
|
||||
|
||||
if (h->record_buf)
|
||||
free(h->record_buf);
|
||||
|
||||
free(h);
|
||||
}
|
||||
|
||||
|
||||
double
|
||||
de_reader_get_const(de_handle *h, const char *name)
|
||||
{
|
||||
/*
|
||||
* Constant lookup requires reading constant names from record 1
|
||||
* and values from record 2. For now, we expose the key values
|
||||
* that are already parsed from the header.
|
||||
*/
|
||||
if (!h || !name)
|
||||
return 0.0;
|
||||
|
||||
if (strcmp(name, "AU") == 0)
|
||||
return h->au_km;
|
||||
if (strcmp(name, "EMRAT") == 0)
|
||||
return h->emrat;
|
||||
|
||||
return 0.0;
|
||||
}
|
||||
140
src/de_reader.h
Normal file
140
src/de_reader.h
Normal file
@ -0,0 +1,140 @@
|
||||
/*
|
||||
* de_reader.h -- Clean-room JPL Development Ephemeris reader
|
||||
*
|
||||
* Reads JPL DE430/DE440/DE441 binary ephemeris files.
|
||||
* Implements Chebyshev polynomial evaluation via Clenshaw recurrence.
|
||||
*
|
||||
* No GPL dependency: written from the public JPL binary format spec.
|
||||
* No global state: each handle is independently opened/managed.
|
||||
*
|
||||
* Reference:
|
||||
* JPL IOM 312.N-03-009 "The JPL Planetary and Lunar Ephemerides,
|
||||
* DE405/LE405" (Standish 1998) — format description.
|
||||
*/
|
||||
|
||||
#ifndef PG_ORBIT_DE_READER_H
|
||||
#define PG_ORBIT_DE_READER_H
|
||||
|
||||
#include <stdint.h>
|
||||
|
||||
/*
|
||||
* JPL DE body targets.
|
||||
*
|
||||
* The ephemeris stores 13 "body groups" in a specific order.
|
||||
* Each group has its own coefficient layout (offset, ncoeff, nsub).
|
||||
*/
|
||||
#define DE_MERCURY 0
|
||||
#define DE_VENUS 1
|
||||
#define DE_EMB 2 /* Earth-Moon Barycenter */
|
||||
#define DE_MARS 3
|
||||
#define DE_JUPITER 4
|
||||
#define DE_SATURN 5
|
||||
#define DE_URANUS 6
|
||||
#define DE_NEPTUNE 7
|
||||
#define DE_PLUTO 8
|
||||
#define DE_MOON 9 /* Moon (geocentric) */
|
||||
#define DE_SUN 10
|
||||
#define DE_NUTATION 11 /* nutations (dpsi, deps) */
|
||||
#define DE_LIBRATION 12 /* lunar librations */
|
||||
#define DE_NUM_BODIES 13
|
||||
|
||||
/*
|
||||
* DE reader error codes.
|
||||
*/
|
||||
#define DE_OK 0
|
||||
#define DE_ERR_OPEN -1 /* cannot open file */
|
||||
#define DE_ERR_READ -2 /* read() failed or short read */
|
||||
#define DE_ERR_HEADER -3 /* header validation failed */
|
||||
#define DE_ERR_RANGE -4 /* JD out of ephemeris range */
|
||||
#define DE_ERR_BODY -5 /* invalid body target */
|
||||
#define DE_ERR_CANARY -6 /* canary validation failed */
|
||||
#define DE_ERR_ENDIAN -7 /* byte order detection failed */
|
||||
|
||||
/*
|
||||
* Coefficient layout for one body group.
|
||||
* Parsed from the header's IPT array.
|
||||
*/
|
||||
typedef struct de_body_layout
|
||||
{
|
||||
int offset; /* word offset within record (1-based) */
|
||||
int ncoeff; /* number of Chebyshev coefficients per component */
|
||||
int nsub; /* number of sub-intervals per record interval */
|
||||
} de_body_layout;
|
||||
|
||||
/*
|
||||
* DE reader handle.
|
||||
*
|
||||
* One per PostgreSQL backend. Owns its file descriptor and
|
||||
* a pre-allocated coefficient buffer.
|
||||
*/
|
||||
typedef struct de_handle
|
||||
{
|
||||
int fd; /* file descriptor */
|
||||
int swap_bytes; /* 1 if byte-swapping needed */
|
||||
|
||||
/* Header metadata */
|
||||
double start_jd; /* first valid JD */
|
||||
double end_jd; /* last valid JD */
|
||||
double interval_days; /* days per record */
|
||||
double au_km; /* AU in km (from header) */
|
||||
double emrat; /* Earth-Moon mass ratio */
|
||||
int ncoeff; /* number of doubles per record */
|
||||
int de_version; /* e.g. 441 */
|
||||
|
||||
/* Coefficient layout per body group */
|
||||
de_body_layout layout[DE_NUM_BODIES];
|
||||
|
||||
/* Record buffer: re-used across queries */
|
||||
double *record_buf; /* ncoeff doubles */
|
||||
int cached_recno; /* which record is loaded, -1 = none */
|
||||
|
||||
/* Record geometry */
|
||||
int record_bytes; /* ncoeff * sizeof(double) */
|
||||
long data_offset; /* byte offset to first data record */
|
||||
} de_handle;
|
||||
|
||||
|
||||
/*
|
||||
* Open and validate a JPL DE binary file.
|
||||
*
|
||||
* Returns a heap-allocated handle on success, NULL on failure.
|
||||
* Sets *errcode to one of DE_OK / DE_ERR_*.
|
||||
* Caller must eventually call de_reader_close().
|
||||
*
|
||||
* Does NOT use ereport() — caller translates error codes.
|
||||
*/
|
||||
de_handle *de_reader_open(const char *path, int *errcode);
|
||||
|
||||
|
||||
/*
|
||||
* Get the position of a body relative to a center.
|
||||
*
|
||||
* target: DE_MERCURY..DE_SUN (0-10)
|
||||
* center: DE_SUN (10) for heliocentric, DE_EMB (2) for geocentric, etc.
|
||||
* Use -1 for "raw" (no center subtraction)
|
||||
* jd: Julian date (TDB)
|
||||
* pos[3]: output position in AU (ICRS equatorial frame)
|
||||
*
|
||||
* Returns DE_OK on success, DE_ERR_* on failure.
|
||||
*
|
||||
* For Earth position: returns Earth (derived from EMB and Moon)
|
||||
* For Moon position with center=Earth: returns geocentric Moon
|
||||
*/
|
||||
int de_reader_get_pos(de_handle *h, double jd, int target, int center,
|
||||
double pos[3]);
|
||||
|
||||
|
||||
/*
|
||||
* Close the DE reader handle and free all resources.
|
||||
* Safe to call with NULL handle.
|
||||
*/
|
||||
void de_reader_close(de_handle *h);
|
||||
|
||||
|
||||
/*
|
||||
* Look up a named constant from the DE header.
|
||||
* Returns the value, or 0.0 if not found.
|
||||
*/
|
||||
double de_reader_get_const(de_handle *h, const char *name);
|
||||
|
||||
#endif /* PG_ORBIT_DE_READER_H */
|
||||
341
src/eph_provider.c
Normal file
341
src/eph_provider.c
Normal file
@ -0,0 +1,341 @@
|
||||
/*
|
||||
* eph_provider.c -- Ephemeris provider management for pg_orbit
|
||||
*
|
||||
* Manages the lifecycle of the optional JPL DE ephemeris reader
|
||||
* within PostgreSQL's multi-process architecture.
|
||||
*
|
||||
* Critical design:
|
||||
* - Each backend (including parallel workers) opens its own
|
||||
* file descriptor via lazy initialization on first DE call.
|
||||
* - The postmaster never opens the file (_PG_init only registers GUCs).
|
||||
* - on_proc_exit callback ensures cleanup.
|
||||
* - ICRS equatorial -> ecliptic J2000 rotation happens here,
|
||||
* at the provider boundary, before data enters the observation
|
||||
* pipeline.
|
||||
*
|
||||
* Constant chain of custody:
|
||||
* - DE returns positions in ICRS equatorial (AU, km internally)
|
||||
* - We rotate to ecliptic J2000 via equatorial_to_ecliptic()
|
||||
* - Both target and Earth always come from the same provider
|
||||
*/
|
||||
|
||||
#include "postgres.h"
|
||||
#include "fmgr.h"
|
||||
#include "utils/guc.h"
|
||||
#include "storage/ipc.h"
|
||||
|
||||
#include "eph_provider.h"
|
||||
#include "de_reader.h"
|
||||
#include "astro_math.h"
|
||||
|
||||
#include <math.h>
|
||||
#include <string.h>
|
||||
|
||||
/* GUC variable: path to DE ephemeris file */
|
||||
static char *eph_path_guc = NULL;
|
||||
|
||||
/*
|
||||
* Per-backend state (safe: each process gets its own copy after fork).
|
||||
* These are NOT shared memory — they live in each backend's address space.
|
||||
*/
|
||||
static de_handle *de_handle_ptr = NULL;
|
||||
static bool de_init_attempted = false;
|
||||
static bool de_init_success = false;
|
||||
|
||||
|
||||
void
|
||||
eph_register_gucs(void)
|
||||
{
|
||||
DefineCustomStringVariable(
|
||||
"pg_orbit.ephemeris_path",
|
||||
"Path to JPL DE binary ephemeris file (DE430/DE440/DE441).",
|
||||
"When set, enables high-precision _de() function variants. "
|
||||
"Default empty = no DE support (VSOP87 only). "
|
||||
"Relative paths are relative to $PGDATA.",
|
||||
&eph_path_guc,
|
||||
"", /* default: empty = disabled */
|
||||
PGC_SIGHUP, /* changeable by SIGHUP reload */
|
||||
GUC_SUPERUSER_ONLY, /* only superuser can set */
|
||||
NULL, /* check_hook */
|
||||
NULL, /* assign_hook */
|
||||
NULL /* show_hook */
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Lazy initialization: open the DE file on first call.
|
||||
* Never called from _PG_init() (postmaster context).
|
||||
*/
|
||||
static void
|
||||
ensure_de_init(void)
|
||||
{
|
||||
int errcode;
|
||||
|
||||
if (de_init_attempted)
|
||||
return;
|
||||
|
||||
de_init_attempted = true;
|
||||
de_init_success = false;
|
||||
|
||||
/* No path configured? DE is simply unavailable. */
|
||||
if (eph_path_guc == NULL || eph_path_guc[0] == '\0')
|
||||
return;
|
||||
|
||||
de_handle_ptr = de_reader_open(eph_path_guc, &errcode);
|
||||
if (de_handle_ptr == NULL)
|
||||
{
|
||||
ereport(NOTICE,
|
||||
(errmsg("pg_orbit: cannot open DE ephemeris at \"%s\" (error %d), "
|
||||
"DE functions will fall back to VSOP87",
|
||||
eph_path_guc, errcode)));
|
||||
return;
|
||||
}
|
||||
|
||||
/* Verify AU consistency with pg_orbit's compiled-in value */
|
||||
if (fabs(de_handle_ptr->au_km - AU_KM) > 0.01)
|
||||
{
|
||||
ereport(NOTICE,
|
||||
(errmsg("pg_orbit: DE ephemeris AU = %.3f km, expected %.3f km; "
|
||||
"DE functions will fall back to VSOP87",
|
||||
de_handle_ptr->au_km, (double)AU_KM)));
|
||||
de_reader_close(de_handle_ptr);
|
||||
de_handle_ptr = NULL;
|
||||
return;
|
||||
}
|
||||
|
||||
ereport(DEBUG1,
|
||||
(errmsg("pg_orbit: DE%d ephemeris loaded, JD %.1f to %.1f",
|
||||
de_handle_ptr->de_version,
|
||||
de_handle_ptr->start_jd,
|
||||
de_handle_ptr->end_jd)));
|
||||
|
||||
de_init_success = true;
|
||||
}
|
||||
|
||||
|
||||
void
|
||||
eph_cleanup(int code, Datum arg)
|
||||
{
|
||||
if (de_handle_ptr != NULL)
|
||||
{
|
||||
de_reader_close(de_handle_ptr);
|
||||
de_handle_ptr = NULL;
|
||||
}
|
||||
de_init_attempted = false;
|
||||
de_init_success = false;
|
||||
}
|
||||
|
||||
|
||||
bool
|
||||
eph_de_available(void)
|
||||
{
|
||||
ensure_de_init();
|
||||
return de_init_success;
|
||||
}
|
||||
|
||||
|
||||
EphProvider
|
||||
eph_current_provider(void)
|
||||
{
|
||||
if (eph_de_available())
|
||||
return EPH_JPL_DE;
|
||||
return EPH_VSOP87;
|
||||
}
|
||||
|
||||
|
||||
const char *
|
||||
eph_get_path(void)
|
||||
{
|
||||
return eph_path_guc;
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Internal: get heliocentric position of a DE body and convert
|
||||
* from ICRS equatorial to ecliptic J2000.
|
||||
*
|
||||
* Returns true on success.
|
||||
*/
|
||||
static bool
|
||||
de_get_helio_ecliptic(int de_target, double jd, double ecl[3])
|
||||
{
|
||||
double icrs[3];
|
||||
int err;
|
||||
|
||||
if (!de_handle_ptr)
|
||||
return false;
|
||||
|
||||
/* Get position relative to Sun (DE_SUN = 10) */
|
||||
err = de_reader_get_pos(de_handle_ptr, jd, de_target, 10, icrs);
|
||||
if (err != DE_OK)
|
||||
return false;
|
||||
|
||||
/* ICRS equatorial -> ecliptic J2000 */
|
||||
equatorial_to_ecliptic(icrs, ecl);
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Internal: get Earth's heliocentric position from DE.
|
||||
* Earth = EMB - Moon/(1+EMRAT), already handled by de_reader.
|
||||
*/
|
||||
static bool
|
||||
de_get_earth_helio_ecliptic(double jd, double ecl[3])
|
||||
{
|
||||
double icrs[3];
|
||||
int err;
|
||||
|
||||
if (!de_handle_ptr)
|
||||
return false;
|
||||
|
||||
/* Earth relative to Sun: use EMB and subtract Moon component */
|
||||
err = de_reader_get_pos(de_handle_ptr, jd, DE_EMB, 10, icrs);
|
||||
if (err != DE_OK)
|
||||
return false;
|
||||
|
||||
/* Correct EMB to Earth: subtract Moon/(1+EMRAT) */
|
||||
{
|
||||
double moon_icrs[3];
|
||||
double factor;
|
||||
|
||||
err = de_reader_get_pos(de_handle_ptr, jd, 9, -1, moon_icrs);
|
||||
if (err != DE_OK)
|
||||
return false;
|
||||
|
||||
/* Moon in DE is geocentric in km->AU. Factor for EMB->Earth */
|
||||
factor = 1.0 / (1.0 + de_handle_ptr->emrat);
|
||||
icrs[0] -= moon_icrs[0] * factor;
|
||||
icrs[1] -= moon_icrs[1] * factor;
|
||||
icrs[2] -= moon_icrs[2] * factor;
|
||||
}
|
||||
|
||||
equatorial_to_ecliptic(icrs, ecl);
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
bool
|
||||
eph_de_planet(int body_id, double jd, double xyz[6])
|
||||
{
|
||||
int de_target;
|
||||
double ecl[3];
|
||||
|
||||
if (!eph_de_available())
|
||||
return false;
|
||||
|
||||
de_target = pgbody_to_de_target(body_id);
|
||||
if (de_target < 0)
|
||||
return false;
|
||||
|
||||
/* Earth is special: derived from EMB and Moon */
|
||||
if (body_id == BODY_EARTH)
|
||||
{
|
||||
if (!de_get_earth_helio_ecliptic(jd, ecl))
|
||||
return false;
|
||||
}
|
||||
else
|
||||
{
|
||||
if (!de_get_helio_ecliptic(de_target, jd, ecl))
|
||||
return false;
|
||||
}
|
||||
|
||||
xyz[0] = ecl[0];
|
||||
xyz[1] = ecl[1];
|
||||
xyz[2] = ecl[2];
|
||||
xyz[3] = 0.0; /* velocity not yet implemented */
|
||||
xyz[4] = 0.0;
|
||||
xyz[5] = 0.0;
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
bool
|
||||
eph_de_earth(double jd, double xyz[6])
|
||||
{
|
||||
double ecl[3];
|
||||
|
||||
if (!eph_de_available())
|
||||
return false;
|
||||
|
||||
if (!de_get_earth_helio_ecliptic(jd, ecl))
|
||||
return false;
|
||||
|
||||
xyz[0] = ecl[0];
|
||||
xyz[1] = ecl[1];
|
||||
xyz[2] = ecl[2];
|
||||
xyz[3] = 0.0;
|
||||
xyz[4] = 0.0;
|
||||
xyz[5] = 0.0;
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
bool
|
||||
eph_de_moon(double jd, double xyz[3])
|
||||
{
|
||||
double icrs[3];
|
||||
int err;
|
||||
|
||||
if (!eph_de_available())
|
||||
return false;
|
||||
|
||||
/* Moon geocentric: target=Moon(9), center=Earth(99) */
|
||||
err = de_reader_get_pos(de_handle_ptr, jd, 9, 99, icrs);
|
||||
if (err != DE_OK)
|
||||
return false;
|
||||
|
||||
/* ICRS -> ecliptic J2000 */
|
||||
equatorial_to_ecliptic(icrs, xyz);
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
bool
|
||||
eph_de_sun(double jd, double xyz[6])
|
||||
{
|
||||
/* Sun is at the origin of heliocentric coordinates */
|
||||
(void)jd;
|
||||
if (!eph_de_available())
|
||||
return false;
|
||||
|
||||
xyz[0] = xyz[1] = xyz[2] = 0.0;
|
||||
xyz[3] = xyz[4] = xyz[5] = 0.0;
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
double
|
||||
eph_de_start_jd(void)
|
||||
{
|
||||
if (de_handle_ptr)
|
||||
return de_handle_ptr->start_jd;
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
double
|
||||
eph_de_end_jd(void)
|
||||
{
|
||||
if (de_handle_ptr)
|
||||
return de_handle_ptr->end_jd;
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
int
|
||||
eph_de_version(void)
|
||||
{
|
||||
if (de_handle_ptr)
|
||||
return de_handle_ptr->de_version;
|
||||
return 0;
|
||||
}
|
||||
|
||||
double
|
||||
eph_de_au_km(void)
|
||||
{
|
||||
if (de_handle_ptr)
|
||||
return de_handle_ptr->au_km;
|
||||
return 0.0;
|
||||
}
|
||||
110
src/eph_provider.h
Normal file
110
src/eph_provider.h
Normal file
@ -0,0 +1,110 @@
|
||||
/*
|
||||
* eph_provider.h -- Ephemeris provider dispatch for pg_orbit
|
||||
*
|
||||
* Manages the optional JPL DE ephemeris alongside the compiled-in
|
||||
* VSOP87 and ELP2000-82B theories.
|
||||
*
|
||||
* Key design constraints:
|
||||
* - Per-backend lazy initialization (never in _PG_init/postmaster)
|
||||
* - ICRS equatorial -> ecliptic J2000 frame rotation at the boundary
|
||||
* - Automatic VSOP87 fallback on any DE error
|
||||
* - No shared state between PostgreSQL backends
|
||||
*/
|
||||
|
||||
#ifndef PG_ORBIT_EPH_PROVIDER_H
|
||||
#define PG_ORBIT_EPH_PROVIDER_H
|
||||
|
||||
#include "postgres.h"
|
||||
#include "types.h"
|
||||
|
||||
/*
|
||||
* Ephemeris provider identifiers.
|
||||
*/
|
||||
typedef enum
|
||||
{
|
||||
EPH_VSOP87, /* compiled-in VSOP87 / ELP2000-82B (always available) */
|
||||
EPH_JPL_DE /* JPL DE binary file (optional, STABLE) */
|
||||
} EphProvider;
|
||||
|
||||
|
||||
/*
|
||||
* Initialize the GUC variable (pg_orbit.ephemeris_path).
|
||||
* Called from _PG_init(). Does NOT open the DE file.
|
||||
*/
|
||||
void eph_register_gucs(void);
|
||||
|
||||
/*
|
||||
* Cleanup callback for on_proc_exit.
|
||||
* Closes the DE file descriptor if open.
|
||||
*/
|
||||
void eph_cleanup(int code, Datum arg);
|
||||
|
||||
/*
|
||||
* Attempt to initialize the DE reader for this backend.
|
||||
* Returns true if DE is available and ready, false otherwise.
|
||||
* On first call, opens the file and validates. On subsequent calls,
|
||||
* returns cached status.
|
||||
*
|
||||
* Thread-safe per backend (each backend has its own static state).
|
||||
*/
|
||||
bool eph_de_available(void);
|
||||
|
||||
/*
|
||||
* Which provider is currently active?
|
||||
*/
|
||||
EphProvider eph_current_provider(void);
|
||||
|
||||
/*
|
||||
* Get the configured ephemeris file path (or NULL if not set).
|
||||
*/
|
||||
const char *eph_get_path(void);
|
||||
|
||||
|
||||
/*
|
||||
* Get planet heliocentric ecliptic J2000 position via DE.
|
||||
*
|
||||
* body_id: pg_orbit body ID (1=Mercury..8=Neptune)
|
||||
* jd: Julian date (TDB)
|
||||
* xyz[6]: output position (AU) in ecliptic J2000 frame
|
||||
* (xyz[3..5] are zero — velocity not yet implemented)
|
||||
*
|
||||
* Returns true on success. On failure, xyz is unchanged and
|
||||
* caller should fall back to VSOP87.
|
||||
*/
|
||||
bool eph_de_planet(int body_id, double jd, double xyz[6]);
|
||||
|
||||
/*
|
||||
* Get Earth heliocentric ecliptic J2000 position via DE.
|
||||
*
|
||||
* jd: Julian date (TDB)
|
||||
* xyz[6]: output position (AU) in ecliptic J2000 frame
|
||||
*
|
||||
* Returns true on success.
|
||||
*/
|
||||
bool eph_de_earth(double jd, double xyz[6]);
|
||||
|
||||
/*
|
||||
* Get Moon geocentric ecliptic J2000 position via DE.
|
||||
*
|
||||
* jd: Julian date (TDB)
|
||||
* xyz[3]: output position (AU) in ecliptic J2000 frame
|
||||
*
|
||||
* Returns true on success. On failure, caller falls back to ELP2000-82B.
|
||||
*/
|
||||
bool eph_de_moon(double jd, double xyz[3]);
|
||||
|
||||
/*
|
||||
* Get Sun "heliocentric" position (always 0,0,0 but consistent API).
|
||||
*/
|
||||
bool eph_de_sun(double jd, double xyz[6]);
|
||||
|
||||
/*
|
||||
* DE metadata accessors (for pg_orbit_ephemeris_info).
|
||||
* Return 0/0.0 if DE is not loaded.
|
||||
*/
|
||||
double eph_de_start_jd(void);
|
||||
double eph_de_end_jd(void);
|
||||
int eph_de_version(void);
|
||||
double eph_de_au_km(void);
|
||||
|
||||
#endif /* PG_ORBIT_EPH_PROVIDER_H */
|
||||
@ -33,38 +33,10 @@ PG_FUNCTION_INFO_V1(uranus_moon_observe);
|
||||
PG_FUNCTION_INFO_V1(mars_moon_observe);
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* Internal: geocentric observation pipeline
|
||||
*
|
||||
* Duplicated from planet_funcs.c to avoid cross-TU coupling.
|
||||
* (Same logic: ecliptic J2000 -> equatorial -> precess -> az/el)
|
||||
* ================================================================
|
||||
/*
|
||||
* observe_from_geocentric() is now in astro_math.h as a static inline,
|
||||
* shared by planet_funcs.c, moon_funcs.c, and de_funcs.c.
|
||||
*/
|
||||
static void
|
||||
moon_observe_from_geocentric(const double geo_ecl_au[3], double jd,
|
||||
const pg_observer *obs, pg_topocentric *result)
|
||||
{
|
||||
double geo_equ[3];
|
||||
double ra_j2000, dec_j2000, geo_dist;
|
||||
double ra_date, dec_date;
|
||||
double gmst_val, lst, ha;
|
||||
double az, el;
|
||||
|
||||
ecliptic_to_equatorial(geo_ecl_au, geo_equ);
|
||||
cartesian_to_spherical(geo_equ, &ra_j2000, &dec_j2000, &geo_dist);
|
||||
precess_j2000_to_date(jd, ra_j2000, dec_j2000, &ra_date, &dec_date);
|
||||
|
||||
gmst_val = gmst_from_jd(jd);
|
||||
lst = gmst_val + obs->lon;
|
||||
ha = lst - ra_date;
|
||||
|
||||
equatorial_to_horizontal(ha, dec_date, obs->lat, &az, &el);
|
||||
|
||||
result->azimuth = az;
|
||||
result->elevation = el;
|
||||
result->range_km = geo_dist * AU_KM;
|
||||
result->range_rate = 0.0;
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
@ -97,7 +69,7 @@ observe_planetary_moon(const double moon_rel[3], int vsop_parent,
|
||||
geo_ecl[1] = (parent_xyz[1] + moon_rel[1]) - earth_xyz[1];
|
||||
geo_ecl[2] = (parent_xyz[2] + moon_rel[2]) - earth_xyz[2];
|
||||
|
||||
moon_observe_from_geocentric(geo_ecl, jd, obs, result);
|
||||
observe_from_geocentric(geo_ecl, jd, obs, result);
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -4,9 +4,30 @@
|
||||
* PostgreSQL extension for orbital mechanics.
|
||||
* Provides TLE, ECI, geodetic, topocentric, observer, and pass_event types
|
||||
* with SGP4/SDP4 propagation, coordinate transforms, and pass prediction.
|
||||
*
|
||||
* v0.3.0 adds _PG_init for GUC registration (pg_orbit.ephemeris_path)
|
||||
* and on_proc_exit cleanup for the optional DE ephemeris handle.
|
||||
*/
|
||||
|
||||
#include "postgres.h"
|
||||
#include "fmgr.h"
|
||||
#include "storage/ipc.h"
|
||||
|
||||
#include "eph_provider.h"
|
||||
|
||||
PG_MODULE_MAGIC;
|
||||
|
||||
/*
|
||||
* _PG_init -- called when the extension shared library is loaded.
|
||||
*
|
||||
* Runs in the postmaster context. Registers GUC variables and
|
||||
* the process exit cleanup callback. Does NOT open the DE file
|
||||
* (that would create a shared fd inherited by all backends via
|
||||
* fork with undefined stream behavior).
|
||||
*/
|
||||
void
|
||||
_PG_init(void)
|
||||
{
|
||||
eph_register_gucs();
|
||||
on_proc_exit(eph_cleanup, (Datum) 0);
|
||||
}
|
||||
|
||||
@ -31,47 +31,10 @@ PG_FUNCTION_INFO_V1(sun_observe);
|
||||
PG_FUNCTION_INFO_V1(moon_observe);
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* Internal: geocentric observation pipeline
|
||||
*
|
||||
* Takes geocentric ecliptic J2000 position, observer location,
|
||||
* and Julian date. Converts through equatorial, precesses to
|
||||
* date, and computes topocentric az/el.
|
||||
*
|
||||
* geo_ecl_au[3] = geocentric ecliptic J2000 position in AU
|
||||
* ================================================================
|
||||
/*
|
||||
* observe_from_geocentric() is now in astro_math.h as a static inline,
|
||||
* shared by planet_funcs.c, moon_funcs.c, and de_funcs.c.
|
||||
*/
|
||||
static void
|
||||
observe_from_geocentric(const double geo_ecl_au[3], double jd,
|
||||
const pg_observer *obs, pg_topocentric *result)
|
||||
{
|
||||
double geo_equ[3];
|
||||
double ra_j2000, dec_j2000, geo_dist;
|
||||
double ra_date, dec_date;
|
||||
double gmst_val, lst, ha;
|
||||
double az, el;
|
||||
|
||||
/* Ecliptic J2000 -> equatorial J2000 */
|
||||
ecliptic_to_equatorial(geo_ecl_au, geo_equ);
|
||||
|
||||
/* Cartesian -> spherical */
|
||||
cartesian_to_spherical(geo_equ, &ra_j2000, &dec_j2000, &geo_dist);
|
||||
|
||||
/* Precess J2000 -> date */
|
||||
precess_j2000_to_date(jd, ra_j2000, dec_j2000, &ra_date, &dec_date);
|
||||
|
||||
/* Hour angle and az/el */
|
||||
gmst_val = gmst_from_jd(jd);
|
||||
lst = gmst_val + obs->lon;
|
||||
ha = lst - ra_date;
|
||||
|
||||
equatorial_to_horizontal(ha, dec_date, obs->lat, &az, &el);
|
||||
|
||||
result->azimuth = az;
|
||||
result->elevation = el;
|
||||
result->range_km = geo_dist * AU_KM;
|
||||
result->range_rate = 0.0; /* no velocity computation yet */
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
|
||||
28
src/types.h
28
src/types.h
@ -222,4 +222,32 @@ typedef struct pg_heliocentric
|
||||
#define BODY_NEPTUNE 8
|
||||
#define BODY_MOON 10
|
||||
|
||||
/*
|
||||
* JPL DE body target IDs (for de_reader.h)
|
||||
* Used by eph_provider.c to translate pg_orbit body IDs to DE targets.
|
||||
*
|
||||
* pg_orbit: 0=Sun, 1=Mercury, ..., 8=Neptune, 10=Moon
|
||||
* DE file: 0=Mercury, ..., 7=Neptune, 8=Pluto, 9=Moon, 10=Sun, 2=EMB
|
||||
*/
|
||||
#define DE_CENTER_EARTH 99 /* virtual: Earth derived from EMB-Moon */
|
||||
|
||||
static inline int
|
||||
pgbody_to_de_target(int body_id)
|
||||
{
|
||||
switch (body_id)
|
||||
{
|
||||
case BODY_SUN: return 10; /* DE_SUN */
|
||||
case BODY_MERCURY: return 0; /* DE_MERCURY */
|
||||
case BODY_VENUS: return 1; /* DE_VENUS */
|
||||
case BODY_EARTH: return 2; /* DE_EMB (Earth-Moon Barycenter) */
|
||||
case BODY_MARS: return 3; /* DE_MARS */
|
||||
case BODY_JUPITER: return 4; /* DE_JUPITER */
|
||||
case BODY_SATURN: return 5; /* DE_SATURN */
|
||||
case BODY_URANUS: return 6; /* DE_URANUS */
|
||||
case BODY_NEPTUNE: return 7; /* DE_NEPTUNE */
|
||||
case BODY_MOON: return 9; /* DE_MOON */
|
||||
default: return -1; /* invalid */
|
||||
}
|
||||
}
|
||||
|
||||
#endif /* PG_ORBIT_TYPES_H */
|
||||
|
||||
186
test/expected/de_ephemeris.out
Normal file
186
test/expected/de_ephemeris.out
Normal file
@ -0,0 +1,186 @@
|
||||
-- de_ephemeris regression tests
|
||||
--
|
||||
-- Tests the _de() function variants and VSOP87 fallback behavior.
|
||||
-- Since DE ephemeris files are not available in CI, these tests
|
||||
-- verify that fallback behavior is correct and produces identical
|
||||
-- results to the VSOP87 variants.
|
||||
\set boulder '''40.015N 105.270W 1655m'''::observer
|
||||
-- ============================================================
|
||||
-- Test 1: pg_orbit_ephemeris_info() returns VSOP87 when no DE file
|
||||
-- The provider should be 'VSOP87' since no ephemeris_path is set.
|
||||
-- ============================================================
|
||||
SELECT 'eph_info_vsop87' AS test,
|
||||
(pg_orbit_ephemeris_info()).provider AS provider;
|
||||
test | provider
|
||||
-----------------+----------
|
||||
eph_info_vsop87 | VSOP87
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Test 2: planet_heliocentric_de falls back to VSOP87
|
||||
-- Should produce identical results to planet_heliocentric()
|
||||
-- when DE is unavailable.
|
||||
-- ============================================================
|
||||
SELECT 'helio_fallback' AS test,
|
||||
round(helio_x(planet_heliocentric(3, '2024-06-21 12:00:00+00'))::numeric, 8) =
|
||||
round(helio_x(planet_heliocentric_de(3, '2024-06-21 12:00:00+00'))::numeric, 8) AS x_match,
|
||||
round(helio_y(planet_heliocentric(3, '2024-06-21 12:00:00+00'))::numeric, 8) =
|
||||
round(helio_y(planet_heliocentric_de(3, '2024-06-21 12:00:00+00'))::numeric, 8) AS y_match,
|
||||
round(helio_z(planet_heliocentric(3, '2024-06-21 12:00:00+00'))::numeric, 8) =
|
||||
round(helio_z(planet_heliocentric_de(3, '2024-06-21 12:00:00+00'))::numeric, 8) AS z_match;
|
||||
test | x_match | y_match | z_match
|
||||
----------------+---------+---------+---------
|
||||
helio_fallback | t | t | t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Test 3: planet_heliocentric_de Sun at origin
|
||||
-- ============================================================
|
||||
SELECT 'sun_origin_de' AS test,
|
||||
round(helio_x(planet_heliocentric_de(0, '2024-06-21 12:00:00+00'))::numeric, 10) AS x,
|
||||
round(helio_y(planet_heliocentric_de(0, '2024-06-21 12:00:00+00'))::numeric, 10) AS y,
|
||||
round(helio_z(planet_heliocentric_de(0, '2024-06-21 12:00:00+00'))::numeric, 10) AS z;
|
||||
test | x | y | z
|
||||
---------------+--------------+--------------+--------------
|
||||
sun_origin_de | 0.0000000000 | 0.0000000000 | 0.0000000000
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Test 4: planet_observe_de falls back to VSOP87
|
||||
-- Elevation and azimuth should match planet_observe().
|
||||
-- ============================================================
|
||||
SELECT 'observe_fallback' AS test,
|
||||
round(topo_azimuth(planet_observe(5, :boulder, '2024-03-15 03:00:00+00'))::numeric, 4) =
|
||||
round(topo_azimuth(planet_observe_de(5, :boulder, '2024-03-15 03:00:00+00'))::numeric, 4) AS az_match,
|
||||
round(topo_elevation(planet_observe(5, :boulder, '2024-03-15 03:00:00+00'))::numeric, 4) =
|
||||
round(topo_elevation(planet_observe_de(5, :boulder, '2024-03-15 03:00:00+00'))::numeric, 4) AS el_match;
|
||||
test | az_match | el_match
|
||||
------------------+----------+----------
|
||||
observe_fallback | t | t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Test 5: sun_observe_de falls back to VSOP87
|
||||
-- ============================================================
|
||||
SELECT 'sun_fallback' AS test,
|
||||
round(topo_azimuth(sun_observe(:boulder, '2024-06-21 18:00:00+00'))::numeric, 4) =
|
||||
round(topo_azimuth(sun_observe_de(:boulder, '2024-06-21 18:00:00+00'))::numeric, 4) AS az_match,
|
||||
round(topo_elevation(sun_observe(:boulder, '2024-06-21 18:00:00+00'))::numeric, 4) =
|
||||
round(topo_elevation(sun_observe_de(:boulder, '2024-06-21 18:00:00+00'))::numeric, 4) AS el_match;
|
||||
test | az_match | el_match
|
||||
--------------+----------+----------
|
||||
sun_fallback | t | t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Test 6: moon_observe_de falls back to ELP2000-82B
|
||||
-- ============================================================
|
||||
SELECT 'moon_fallback' AS test,
|
||||
round(topo_azimuth(moon_observe(:boulder, '2024-06-21 18:00:00+00'))::numeric, 4) =
|
||||
round(topo_azimuth(moon_observe_de(:boulder, '2024-06-21 18:00:00+00'))::numeric, 4) AS az_match,
|
||||
round(topo_range(moon_observe(:boulder, '2024-06-21 18:00:00+00'))::numeric, 0) =
|
||||
round(topo_range(moon_observe_de(:boulder, '2024-06-21 18:00:00+00'))::numeric, 0) AS range_match;
|
||||
test | az_match | range_match
|
||||
---------------+----------+-------------
|
||||
moon_fallback | t | t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Test 7: lambert_c3_de falls back to VSOP87
|
||||
-- Earth-Mars 2024 window should match lambert_c3().
|
||||
-- ============================================================
|
||||
SELECT 'lambert_fallback' AS test,
|
||||
round(lambert_c3(3, 4, '2024-05-01 00:00:00+00', '2025-02-01 00:00:00+00')::numeric, 2) =
|
||||
round(lambert_c3_de(3, 4, '2024-05-01 00:00:00+00', '2025-02-01 00:00:00+00')::numeric, 2) AS c3_match;
|
||||
test | c3_match
|
||||
------------------+----------
|
||||
lambert_fallback | t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Test 8: lambert_transfer_de falls back to VSOP87
|
||||
-- Should produce identical departure C3.
|
||||
-- ============================================================
|
||||
SELECT 'transfer_fallback' AS test,
|
||||
round((lambert_transfer(3, 4, '2024-05-01 00:00:00+00', '2025-02-01 00:00:00+00')).c3_departure::numeric, 2) =
|
||||
round((lambert_transfer_de(3, 4, '2024-05-01 00:00:00+00', '2025-02-01 00:00:00+00')).c3_departure::numeric, 2) AS c3_match;
|
||||
test | c3_match
|
||||
-------------------+----------
|
||||
transfer_fallback | t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Test 9: galilean_observe_de falls back to VSOP87
|
||||
-- ============================================================
|
||||
SELECT 'galilean_fallback' AS test,
|
||||
round(topo_elevation(galilean_observe(0, :boulder, '2024-03-15 03:00:00+00'))::numeric, 4) =
|
||||
round(topo_elevation(galilean_observe_de(0, :boulder, '2024-03-15 03:00:00+00'))::numeric, 4) AS el_match;
|
||||
test | el_match
|
||||
-------------------+----------
|
||||
galilean_fallback | t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Test 10: saturn_moon_observe_de falls back to VSOP87
|
||||
-- ============================================================
|
||||
SELECT 'saturn_moon_fallback' AS test,
|
||||
round(topo_elevation(saturn_moon_observe(5, :boulder, '2024-06-15 04:00:00+00'))::numeric, 4) =
|
||||
round(topo_elevation(saturn_moon_observe_de(5, :boulder, '2024-06-15 04:00:00+00'))::numeric, 4) AS el_match;
|
||||
test | el_match
|
||||
----------------------+----------
|
||||
saturn_moon_fallback | t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Test 11: uranus_moon_observe_de falls back to VSOP87
|
||||
-- ============================================================
|
||||
SELECT 'uranus_moon_fallback' AS test,
|
||||
round(topo_elevation(uranus_moon_observe(3, :boulder, '2024-06-15 04:00:00+00'))::numeric, 4) =
|
||||
round(topo_elevation(uranus_moon_observe_de(3, :boulder, '2024-06-15 04:00:00+00'))::numeric, 4) AS el_match;
|
||||
test | el_match
|
||||
----------------------+----------
|
||||
uranus_moon_fallback | t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Test 12: mars_moon_observe_de falls back to VSOP87
|
||||
-- ============================================================
|
||||
SELECT 'mars_moon_fallback' AS test,
|
||||
round(topo_elevation(mars_moon_observe(0, :boulder, '2024-06-15 04:00:00+00'))::numeric, 4) =
|
||||
round(topo_elevation(mars_moon_observe_de(0, :boulder, '2024-06-15 04:00:00+00'))::numeric, 4) AS el_match;
|
||||
test | el_match
|
||||
--------------------+----------
|
||||
mars_moon_fallback | t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Test 13: All DE planet functions work (fallback mode)
|
||||
-- Mercury through Neptune, all should return valid positions
|
||||
-- matching their VSOP87 counterparts.
|
||||
-- ============================================================
|
||||
SELECT 'all_planets_de' AS test,
|
||||
body_id,
|
||||
round(helio_distance(planet_heliocentric_de(body_id, '2024-06-21 12:00:00+00'))::numeric, 2) AS dist_au
|
||||
FROM generate_series(1, 8) AS body_id;
|
||||
test | body_id | dist_au
|
||||
----------------+---------+---------
|
||||
all_planets_de | 1 | 0.33
|
||||
all_planets_de | 2 | 0.72
|
||||
all_planets_de | 3 | 1.02
|
||||
all_planets_de | 4 | 1.40
|
||||
all_planets_de | 5 | 5.02
|
||||
all_planets_de | 6 | 9.69
|
||||
all_planets_de | 7 | 19.59
|
||||
all_planets_de | 8 | 29.90
|
||||
(8 rows)
|
||||
|
||||
-- ============================================================
|
||||
-- Test 14: Error handling - invalid body_id
|
||||
-- ============================================================
|
||||
SELECT 'invalid_body_de' AS test, planet_heliocentric_de(9, now());
|
||||
ERROR: invalid body_id 9: must be 0 (Sun) or 1-8 (Mercury-Neptune)
|
||||
-- ============================================================
|
||||
-- Test 15: Error handling - cannot observe Earth from Earth
|
||||
-- ============================================================
|
||||
SELECT 'earth_error_de' AS test, planet_observe_de(3, :boulder, now());
|
||||
ERROR: cannot observe Earth from Earth
|
||||
128
test/sql/de_ephemeris.sql
Normal file
128
test/sql/de_ephemeris.sql
Normal file
@ -0,0 +1,128 @@
|
||||
-- de_ephemeris regression tests
|
||||
--
|
||||
-- Tests the _de() function variants and VSOP87 fallback behavior.
|
||||
-- Since DE ephemeris files are not available in CI, these tests
|
||||
-- verify that fallback behavior is correct and produces identical
|
||||
-- results to the VSOP87 variants.
|
||||
|
||||
\set boulder '''40.015N 105.270W 1655m'''::observer
|
||||
|
||||
-- ============================================================
|
||||
-- Test 1: pg_orbit_ephemeris_info() returns VSOP87 when no DE file
|
||||
-- The provider should be 'VSOP87' since no ephemeris_path is set.
|
||||
-- ============================================================
|
||||
SELECT 'eph_info_vsop87' AS test,
|
||||
(pg_orbit_ephemeris_info()).provider AS provider;
|
||||
|
||||
-- ============================================================
|
||||
-- Test 2: planet_heliocentric_de falls back to VSOP87
|
||||
-- Should produce identical results to planet_heliocentric()
|
||||
-- when DE is unavailable.
|
||||
-- ============================================================
|
||||
SELECT 'helio_fallback' AS test,
|
||||
round(helio_x(planet_heliocentric(3, '2024-06-21 12:00:00+00'))::numeric, 8) =
|
||||
round(helio_x(planet_heliocentric_de(3, '2024-06-21 12:00:00+00'))::numeric, 8) AS x_match,
|
||||
round(helio_y(planet_heliocentric(3, '2024-06-21 12:00:00+00'))::numeric, 8) =
|
||||
round(helio_y(planet_heliocentric_de(3, '2024-06-21 12:00:00+00'))::numeric, 8) AS y_match,
|
||||
round(helio_z(planet_heliocentric(3, '2024-06-21 12:00:00+00'))::numeric, 8) =
|
||||
round(helio_z(planet_heliocentric_de(3, '2024-06-21 12:00:00+00'))::numeric, 8) AS z_match;
|
||||
|
||||
-- ============================================================
|
||||
-- Test 3: planet_heliocentric_de Sun at origin
|
||||
-- ============================================================
|
||||
SELECT 'sun_origin_de' AS test,
|
||||
round(helio_x(planet_heliocentric_de(0, '2024-06-21 12:00:00+00'))::numeric, 10) AS x,
|
||||
round(helio_y(planet_heliocentric_de(0, '2024-06-21 12:00:00+00'))::numeric, 10) AS y,
|
||||
round(helio_z(planet_heliocentric_de(0, '2024-06-21 12:00:00+00'))::numeric, 10) AS z;
|
||||
|
||||
-- ============================================================
|
||||
-- Test 4: planet_observe_de falls back to VSOP87
|
||||
-- Elevation and azimuth should match planet_observe().
|
||||
-- ============================================================
|
||||
SELECT 'observe_fallback' AS test,
|
||||
round(topo_azimuth(planet_observe(5, :boulder, '2024-03-15 03:00:00+00'))::numeric, 4) =
|
||||
round(topo_azimuth(planet_observe_de(5, :boulder, '2024-03-15 03:00:00+00'))::numeric, 4) AS az_match,
|
||||
round(topo_elevation(planet_observe(5, :boulder, '2024-03-15 03:00:00+00'))::numeric, 4) =
|
||||
round(topo_elevation(planet_observe_de(5, :boulder, '2024-03-15 03:00:00+00'))::numeric, 4) AS el_match;
|
||||
|
||||
-- ============================================================
|
||||
-- Test 5: sun_observe_de falls back to VSOP87
|
||||
-- ============================================================
|
||||
SELECT 'sun_fallback' AS test,
|
||||
round(topo_azimuth(sun_observe(:boulder, '2024-06-21 18:00:00+00'))::numeric, 4) =
|
||||
round(topo_azimuth(sun_observe_de(:boulder, '2024-06-21 18:00:00+00'))::numeric, 4) AS az_match,
|
||||
round(topo_elevation(sun_observe(:boulder, '2024-06-21 18:00:00+00'))::numeric, 4) =
|
||||
round(topo_elevation(sun_observe_de(:boulder, '2024-06-21 18:00:00+00'))::numeric, 4) AS el_match;
|
||||
|
||||
-- ============================================================
|
||||
-- Test 6: moon_observe_de falls back to ELP2000-82B
|
||||
-- ============================================================
|
||||
SELECT 'moon_fallback' AS test,
|
||||
round(topo_azimuth(moon_observe(:boulder, '2024-06-21 18:00:00+00'))::numeric, 4) =
|
||||
round(topo_azimuth(moon_observe_de(:boulder, '2024-06-21 18:00:00+00'))::numeric, 4) AS az_match,
|
||||
round(topo_range(moon_observe(:boulder, '2024-06-21 18:00:00+00'))::numeric, 0) =
|
||||
round(topo_range(moon_observe_de(:boulder, '2024-06-21 18:00:00+00'))::numeric, 0) AS range_match;
|
||||
|
||||
-- ============================================================
|
||||
-- Test 7: lambert_c3_de falls back to VSOP87
|
||||
-- Earth-Mars 2024 window should match lambert_c3().
|
||||
-- ============================================================
|
||||
SELECT 'lambert_fallback' AS test,
|
||||
round(lambert_c3(3, 4, '2024-05-01 00:00:00+00', '2025-02-01 00:00:00+00')::numeric, 2) =
|
||||
round(lambert_c3_de(3, 4, '2024-05-01 00:00:00+00', '2025-02-01 00:00:00+00')::numeric, 2) AS c3_match;
|
||||
|
||||
-- ============================================================
|
||||
-- Test 8: lambert_transfer_de falls back to VSOP87
|
||||
-- Should produce identical departure C3.
|
||||
-- ============================================================
|
||||
SELECT 'transfer_fallback' AS test,
|
||||
round((lambert_transfer(3, 4, '2024-05-01 00:00:00+00', '2025-02-01 00:00:00+00')).c3_departure::numeric, 2) =
|
||||
round((lambert_transfer_de(3, 4, '2024-05-01 00:00:00+00', '2025-02-01 00:00:00+00')).c3_departure::numeric, 2) AS c3_match;
|
||||
|
||||
-- ============================================================
|
||||
-- Test 9: galilean_observe_de falls back to VSOP87
|
||||
-- ============================================================
|
||||
SELECT 'galilean_fallback' AS test,
|
||||
round(topo_elevation(galilean_observe(0, :boulder, '2024-03-15 03:00:00+00'))::numeric, 4) =
|
||||
round(topo_elevation(galilean_observe_de(0, :boulder, '2024-03-15 03:00:00+00'))::numeric, 4) AS el_match;
|
||||
|
||||
-- ============================================================
|
||||
-- Test 10: saturn_moon_observe_de falls back to VSOP87
|
||||
-- ============================================================
|
||||
SELECT 'saturn_moon_fallback' AS test,
|
||||
round(topo_elevation(saturn_moon_observe(5, :boulder, '2024-06-15 04:00:00+00'))::numeric, 4) =
|
||||
round(topo_elevation(saturn_moon_observe_de(5, :boulder, '2024-06-15 04:00:00+00'))::numeric, 4) AS el_match;
|
||||
|
||||
-- ============================================================
|
||||
-- Test 11: uranus_moon_observe_de falls back to VSOP87
|
||||
-- ============================================================
|
||||
SELECT 'uranus_moon_fallback' AS test,
|
||||
round(topo_elevation(uranus_moon_observe(3, :boulder, '2024-06-15 04:00:00+00'))::numeric, 4) =
|
||||
round(topo_elevation(uranus_moon_observe_de(3, :boulder, '2024-06-15 04:00:00+00'))::numeric, 4) AS el_match;
|
||||
|
||||
-- ============================================================
|
||||
-- Test 12: mars_moon_observe_de falls back to VSOP87
|
||||
-- ============================================================
|
||||
SELECT 'mars_moon_fallback' AS test,
|
||||
round(topo_elevation(mars_moon_observe(0, :boulder, '2024-06-15 04:00:00+00'))::numeric, 4) =
|
||||
round(topo_elevation(mars_moon_observe_de(0, :boulder, '2024-06-15 04:00:00+00'))::numeric, 4) AS el_match;
|
||||
|
||||
-- ============================================================
|
||||
-- Test 13: All DE planet functions work (fallback mode)
|
||||
-- Mercury through Neptune, all should return valid positions
|
||||
-- matching their VSOP87 counterparts.
|
||||
-- ============================================================
|
||||
SELECT 'all_planets_de' AS test,
|
||||
body_id,
|
||||
round(helio_distance(planet_heliocentric_de(body_id, '2024-06-21 12:00:00+00'))::numeric, 2) AS dist_au
|
||||
FROM generate_series(1, 8) AS body_id;
|
||||
|
||||
-- ============================================================
|
||||
-- Test 14: Error handling - invalid body_id
|
||||
-- ============================================================
|
||||
SELECT 'invalid_body_de' AS test, planet_heliocentric_de(9, now());
|
||||
|
||||
-- ============================================================
|
||||
-- Test 15: Error handling - cannot observe Earth from Earth
|
||||
-- ============================================================
|
||||
SELECT 'earth_error_de' AS test, planet_observe_de(3, :boulder, now());
|
||||
Loading…
x
Reference in New Issue
Block a user