diff --git a/CLAUDE.md b/CLAUDE.md
index c0868c8..977bd03 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -1,9 +1,9 @@
# pg_orrery — A Database Orrery for PostgreSQL
## What This Is
-A database orrery — celestial mechanics types and functions for PostgreSQL. Native C extension using PGXS, 151 SQL objects (135 user-visible functions + 16 GiST support), 9 custom types, covering satellites (SGP4/SDP4), planets (VSOP87 + optional JPL DE441), Moon (ELP2000-82B), 19 planetary moons (L1.2/TASS17/GUST86/MarsSat), stars (with proper motion and annual parallax), comets, asteroids (MPC catalog), Jupiter radio bursts, interplanetary Lambert transfers, equatorial RA/Dec coordinates with GiST-indexed angular separation, atmospheric refraction, annual stellar aberration, light-time correction, rise/set prediction (geometric + refracted) with status diagnostics, and IAU constellation identification with full name lookup (Roman 1987).
+A database orrery — celestial mechanics types and functions for PostgreSQL. Native C extension using PGXS, 162 SQL objects (146 user-visible functions + 16 GiST support), 9 custom types, covering satellites (SGP4/SDP4), planets (VSOP87 + optional JPL DE441), Moon (ELP2000-82B), 19 planetary moons (L1.2/TASS17/GUST86/MarsSat), stars (with proper motion and annual parallax), comets, asteroids (MPC catalog), Jupiter radio bursts, interplanetary Lambert transfers, equatorial RA/Dec coordinates with GiST-indexed angular separation, atmospheric refraction, annual stellar aberration, light-time correction, rise/set prediction (geometric + refracted) with status diagnostics, IAU constellation identification with full name lookup (Roman 1987), twilight dawn/dusk (civil/nautical/astronomical), lunar phase (angle, illumination, name, age), and planet apparent magnitude (Mallama & Hilton 2018).
-**Current version:** 0.15.0
+**Current version:** 0.16.0
**Repository:** https://git.supported.systems/warehack.ing/pg_orrery
**Documentation:** https://pg-orrery.warehack.ing
@@ -11,7 +11,7 @@ A database orrery — celestial mechanics types and functions for PostgreSQL. Na
```bash
make PG_CONFIG=/usr/bin/pg_config # Compile with PGXS
sudo make install PG_CONFIG=/usr/bin/pg_config # Install extension
-make installcheck PG_CONFIG=/usr/bin/pg_config # Run 26 regression test suites
+make installcheck PG_CONFIG=/usr/bin/pg_config # Run 27 regression test suites
```
Requires: PostgreSQL 17 development headers, GCC, Make.
@@ -27,7 +27,7 @@ Image: `git.supported.systems/warehack.ing/pg_orrery:pg17`
## Project Layout
```
-pg_orrery.control # Extension metadata (version 0.15.0)
+pg_orrery.control # Extension metadata (version 0.16.0)
Makefile # PGXS build + Docker targets
sql/
pg_orrery--0.1.0.sql # v0.1.0: satellite types/functions/operators
@@ -45,6 +45,7 @@ sql/
pg_orrery--0.13.0.sql # v0.13.0: nutation, make_equatorial, rise/set (141 objects)
pg_orrery--0.14.0.sql # v0.14.0: refracted rise/set, constellation ID (147 objects)
pg_orrery--0.15.0.sql # v0.15.0: constellation full name, rise/set status (151 objects)
+ pg_orrery--0.16.0.sql # v0.16.0: twilight, lunar phase, planet magnitude (162 objects)
pg_orrery--0.1.0--0.2.0.sql # Migration: v0.1.0 → v0.2.0 (adds solar system)
pg_orrery--0.2.0--0.3.0.sql # Migration: v0.2.0 → v0.3.0 (adds DE ephemeris)
pg_orrery--0.3.0--0.4.0.sql # Migration: v0.3.0 → v0.4.0
@@ -59,6 +60,7 @@ sql/
pg_orrery--0.12.0--0.13.0.sql # Migration: v0.12.0 → v0.13.0 (nutation, make_equatorial, rise/set)
pg_orrery--0.13.0--0.14.0.sql # Migration: v0.13.0 → v0.14.0 (refracted rise/set, constellation ID)
pg_orrery--0.14.0--0.15.0.sql # Migration: v0.14.0 → v0.15.0 (constellation full name, rise/set status)
+ pg_orrery--0.15.0--0.16.0.sql # Migration: v0.15.0 → v0.16.0 (twilight, lunar phase, planet magnitude)
src/
pg_orrery.c # PG_MODULE_MAGIC + _PG_init() (GUC registration)
types.h # All struct definitions + constants + DE body ID mapping
@@ -85,9 +87,11 @@ src/
orbital_elements_type.c # orbital_elements type, MPC parser, small_body_observe/equatorial/apparent()
equatorial_funcs.c # equatorial type I/O, accessors, satellite/planet/sun/moon RA/Dec
refraction_funcs.c # atmospheric_refraction(), _ext(), topo_elevation_apparent()
- rise_set_funcs.c # planet/sun/moon rise/set (geometric + refracted)
+ rise_set_funcs.c # planet/sun/moon rise/set (geometric + refracted) + twilight dawn/dusk
constellation_data.h / .c # Roman (1987) IAU boundary table (CDS VI/42, 357 segments)
constellation_funcs.c # constellation() from equatorial or RA/Dec
+ lunar_phase_funcs.c # moon_phase_angle(), moon_illumination(), moon_phase_name(), moon_age()
+ magnitude_funcs.c # planet_magnitude() (Mallama & Hilton 2018)
l12.c / l12.h # L1.2 Galilean moon theory (Lieske 1998)
tass17.c / tass17.h # TASS 1.7 Saturn moon theory (Vienne & Duriez 1995)
gust86.c / gust86.h # GUST86 Uranus moon theory (Laskar & Jacobson 1987)
@@ -112,7 +116,7 @@ src/
PROVENANCE.md # Vendoring decision, modifications, verification
LICENSE # MIT license (Bill Gray / Project Pluto)
test/
- sql/ # 26 regression test suites
+ sql/ # 27 regression test suites
expected/ # Expected output
data/vallado_518.json # 518 Vallado test vectors (AIAA 2006-6753-Rev1)
docs/
@@ -139,7 +143,7 @@ All types are fixed-size, `STORAGE = plain`, `ALIGNMENT = double`. No TOAST over
| `orbital_elements` | 72 | Classical Keplerian elements for comets/asteroids (epoch, q, e, inc, omega, Omega, tp, H, G) |
| `equatorial` | 24 | Apparent RA (hours), Dec (degrees), distance (km) — of date |
-## Function Domains (151 SQL objects)
+## Function Domains (162 SQL objects)
| Domain | Theory | Key Functions | Count |
|--------|--------|---------------|-------|
@@ -157,6 +161,9 @@ All types are fixed-size, `STORAGE = plain`, `ALIGNMENT = double`. No TOAST over
| GiST index (TLE) | Altitude-band approximation | `&&` (overlap), `<->` (distance) | 8 |
| GiST index (equatorial) | Spherical bounding box | `<->` (KNN ordering) | 8 |
| Rise/set | Bisection (60s scan) | `planet_next_rise()`, `sun_next_rise_refracted()`, `*_rise_set_status()` | 15 |
+| Twilight | Sun depression angles | `sun_civil_dawn()`, `sun_nautical_dusk()`, `sun_astronomical_dawn()` | 6 |
+| Lunar phase | VSOP87 + ELP2000-82B geometry | `moon_phase_angle()`, `moon_illumination()`, `moon_phase_name()`, `moon_age()` | 4 |
+| Planet magnitude | Mallama & Hilton (2018) | `planet_magnitude()` | 1 |
| Constellation | Roman (1987) CDS VI/42 | `constellation()`, `constellation_full_name()` | 3 |
| Diagnostics | -- | `pg_orrery_ephemeris_info()` | 1 |
@@ -291,7 +298,7 @@ All numerical logic is byte-identical to upstream. Verified against 518 Vallado
## Testing
-26 regression test suites via `make installcheck`:
+27 regression test suites via `make installcheck`:
| Suite | What it tests |
|-------|--------------|
@@ -321,10 +328,11 @@ All numerical logic is byte-identical to upstream. Verified against 518 Vallado
| rise_set | Planet/Sun/Moon rise/set (geometric + refracted), circumpolar, polar night |
| constellation | Roman (1987) boundary lookup, known stars, solar system objects, edge cases |
| v015_features | constellation_full_name lookup, rise_set_status diagnostics (circumpolar/never_rises) |
+| v016_features | Twilight ordering/offset/polar, lunar phase at known events, planet magnitude ranges/errors |
### PG Version Matrix
-Test all 26 regression suites + DE reader unit test across PostgreSQL 14-18 using Docker:
+Test all 27 regression suites + DE reader unit test across PostgreSQL 14-18 using Docker:
```bash
make test-matrix # Full matrix (PG 14-18)
@@ -350,7 +358,7 @@ Logs saved to `test/matrix-logs/pg${ver}.log`. The script reuses the Dockerfile
Starlight docs at `docs/` — 44+ MDX pages covering all domains.
-Sections: Getting Started, Guides (9 domain walkthroughs incl. DE ephemeris), Workflow Translation (Skyfield/Horizons/GMAT/Radio Jupiter Pro comparisons), Reference (all 151 SQL objects incl. DE variants, equatorial GiST, refraction, rise/set, constellation), 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 162 SQL objects incl. DE variants, equatorial GiST, refraction, rise/set, constellation, twilight, lunar phase, planet magnitude), Architecture (Hamilton's principles, constant custody, observation pipeline), Performance (benchmarks).
### Local Development
```bash
diff --git a/Makefile b/Makefile
index 955a366..6d7e773 100644
--- a/Makefile
+++ b/Makefile
@@ -13,7 +13,8 @@ DATA = sql/pg_orrery--0.1.0.sql sql/pg_orrery--0.2.0.sql sql/pg_orrery--0.1.0--0
sql/pg_orrery--0.12.0.sql sql/pg_orrery--0.11.0--0.12.0.sql \
sql/pg_orrery--0.13.0.sql sql/pg_orrery--0.12.0--0.13.0.sql \
sql/pg_orrery--0.14.0.sql sql/pg_orrery--0.13.0--0.14.0.sql \
- sql/pg_orrery--0.15.0.sql sql/pg_orrery--0.14.0--0.15.0.sql
+ sql/pg_orrery--0.15.0.sql sql/pg_orrery--0.14.0--0.15.0.sql \
+ sql/pg_orrery--0.16.0.sql sql/pg_orrery--0.15.0--0.16.0.sql
# Our extension C sources
OBJS = src/pg_orrery.o src/tle_type.o src/eci_type.o src/observer_type.o \
@@ -32,7 +33,8 @@ OBJS = src/pg_orrery.o src/tle_type.o src/eci_type.o src/observer_type.o \
src/refraction_funcs.o \
src/gist_equatorial.o \
src/rise_set_funcs.o \
- src/constellation_data.o src/constellation_funcs.o
+ src/constellation_data.o src/constellation_funcs.o \
+ src/lunar_phase_funcs.o src/magnitude_funcs.o
# Vendored SGP4/SDP4 sources (pure C, from Bill Gray's sat_code, MIT license)
SGP4_DIR = src/sgp4
@@ -52,7 +54,8 @@ REGRESS = tle_parse sgp4_propagate coord_transforms pass_prediction gist_index c
gist_equatorial v012_features \
v013_features rise_set \
constellation \
- v015_features
+ v015_features \
+ v016_features
REGRESS_OPTS = --inputdir=test
# Pure C — no C++ runtime needed. LAPACK for OD solver (dgelss_).
diff --git a/docs/TODO-v0.10.0.md b/docs/TODO-v0.10.0.md
deleted file mode 100644
index 8633541..0000000
--- a/docs/TODO-v0.10.0.md
+++ /dev/null
@@ -1,138 +0,0 @@
-# pg_orrery — Post-v0.10.0 Roadmap
-
-## Current State
-
-**Version:** v0.10.0 (tagged 2026-02-22)
-**Branch:** `phase/spgist-orbital-trie` merged to `main`
-**Functions:** 114 SQL functions, 9 custom types, 19 test suites, 1 new operator
-**Docs:** https://pg-orrery.warehack.ing
-
-### What v0.10.0 shipped
-- Annual stellar aberration in all `_apparent()` functions (~20 arcsec)
-- 6 new `_apparent_de()` variants with VSOP87 fallback
-- `eq_angular_distance()` + `eq_within_cone()` + `<->` operator on equatorial
-- Stellar annual parallax in `star_observe_pm()` / `star_equatorial_pm()`
-
-### Astrolock integration status
-Thread: `docs/agent-threads/v090-astrolock-upgrade/`
-- v0.9.0 fully deployed to both local and production servers
-- v0.10.0 upgrade path communicated (message 003)
-- Pending their upgrade — aberration improvement is automatic
-
-## Remaining Housekeeping
-
-- [x] Merge `phase/spgist-orbital-trie` to `main`
-- [ ] Clean up `bench/` — gitignore the untracked TLE catalog files (alpha5, celestrak, satnogs, spacetrack, supgp, tle_api, merged, mega)
-- [ ] Update "From Skyfield" workflow page for v0.9.0/v0.10.0 RA/Dec + aberration parity
-- [ ] Add timing numbers for equatorial, refraction, aberration functions to benchmarks page
-- [ ] Update CLAUDE.md function count: 106 -> 114, test suites: 18 -> 19
-- [ ] Update docs llms.txt and llms-full.txt for v0.10.0 functions
-
-## Feature Candidates — Next Version
-
-### Tier 1 — High value, low effort
-
-#### A. `make_orbital_elements()` constructor
-**Requested by:** astrolock-api (message 002, question 1)
-
-SQL constructor from 9 floats. Lets users compose orbital_elements from individual table columns without `format()`/cast workaround.
-
-```sql
-make_orbital_elements(epoch_jd, q_au, e, inc_rad, omega_rad, node_rad, tp_jd, h_mag, g_slope)
- -> orbital_elements
-```
-
-Complexity: ~30 lines in `orbital_elements_type.c`. One new function.
-
-#### B. `galilean_equatorial()` and moon family equatorial functions
-**Requested by:** astrolock-api (message 002, question 2)
-
-Geocentric RA/Dec for planetary moons. Follows `planet_equatorial()` pattern — convert geocentric ecliptic position to equatorial J2000, precess to date.
-
-New functions (~4):
-- `galilean_equatorial(int4, timestamptz) -> equatorial`
-- `saturn_moon_equatorial(int4, timestamptz) -> equatorial`
-- `uranus_moon_equatorial(int4, timestamptz) -> equatorial`
-- `mars_moon_equatorial(int4, timestamptz) -> equatorial`
-
-Plus DE variants (~4 more).
-
-Complexity: ~100 lines. Follows established pattern.
-
-#### C. GiST/SP-GiST index on equatorial type
-The `<->` operator and `eq_within_cone()` exist but have no index support. For cone-search queries over large catalogs, an index would enable:
-
-```sql
--- Indexed: "what's within 10 deg of Jupiter?"
-SELECT * FROM star_catalog
-WHERE position <-> planet_equatorial(5, NOW()) < 10.0;
-```
-
-Approach: GiST with bounding-box approximation in RA/Dec space, or SP-GiST with HEALPix-style recursive decomposition.
-
-Complexity: Medium (~300-500 lines). The SP-GiST infrastructure from TLE index is reusable.
-
-### Tier 2 — Medium value, medium risk
-
-#### D. Nutation correction (~9 arcsec)
-IAU 1980 nutation (106 terms) or simplified IAU 2000B.
-
-Currently: TEME uses 4 of 106 terms. Equatorial output uses IAU 1976 precession only (no nutation).
-
-Value: ~9 arcsec correction in equatorial coordinates. Matters for sub-arcminute accuracy — telescope GoTo mounts and catalog cross-matching.
-
-Scope: New `nutation.c` + modify `precess_j2000_to_date()` to include nutation matrix.
-Risk: Touches the precession pipeline used by every equatorial function.
-
-#### E. Delta T (TDB - UTC)
-The "affects everything" change. Currently all time is treated as UTC with no TT/TDB distinction.
-
-Requires IERS lookup table or polynomial approximation (Espenak & Meeus 2006).
-
-Scope: Touch `sidereal_time.h`, propagation pipelines, all observation functions.
-Risk: High — affects every time conversion. Needs careful regression testing.
-Value: Improves accuracy for historical epochs (pre-2000) and future predictions (post-2030).
-
-Already noted as deferred at `sidereal_time.h:22-26`.
-
-#### F. Rise/set prediction for solar system objects
-Like `predict_passes()` but for planets, Sun, and Moon. Binary search for horizon crossings.
-
-Use cases: sunrise/sunset, moonrise/moonset, planet visibility windows.
-Complexity: Medium. The pass prediction binary search machinery exists but needs adaptation for much slower angular rates.
-
-### Tier 3 — Future / deferred
-
-- **Perturbed asteroid propagation** — secular perturbation terms for orbital_elements (currently two-body Keplerian)
-- **Eclipse prediction** — Moon shadow cone intersection with observer
-- **Satellite sunlit visibility** — extend `pass_visible()` with Earth shadow geometry
-- **Constellation identification** — equatorial position to IAU constellation boundary lookup
-- **Coordinate frame transforms** — ICRS/FK5/galactic/ecliptic conversion functions
-
-## Suggested Next Phase
-
-```
-Housekeeping (bench cleanup, docs, CLAUDE.md)
- |
- v
-Feature A: make_orbital_elements() — 30 lines, unblocks Craft comets
-Feature B: moon family equatorial — 100 lines, unblocks Craft Galilean moons
- |
- v
-Feature C: equatorial GiST index — enables indexed cone search
-Feature D: nutation — closes largest remaining accuracy gap (~9 arcsec)
- |
- v
-Feature E: Delta T — high risk, needs its own careful phase
-Feature F: rise/set — new domain, independent
-```
-
-A + B could ship as v0.11.0. C + D as v0.12.0.
-
-## Verification
-
-- All 19 existing regression suites must continue to pass
-- New test suites for each feature
-- PG 14-18 matrix (`make test-matrix`)
-- Cross-check against JPL Horizons for nutation accuracy
-- Astrolock integration smoke test after each db upgrade
diff --git a/docs/astro.config.mjs b/docs/astro.config.mjs
index 50104af..918ba54 100644
--- a/docs/astro.config.mjs
+++ b/docs/astro.config.mjs
@@ -72,6 +72,8 @@ export default defineConfig({
{ label: "Orbit Determination", slug: "guides/orbit-determination" },
{ label: "Satellite Pass Prediction", slug: "guides/pass-prediction" },
{ label: "Building TLE Catalogs", slug: "guides/catalog-management" },
+ { label: "Rise/Set Prediction", slug: "guides/rise-set-prediction" },
+ { label: "Constellation Identification", slug: "guides/constellation-identification" },
],
},
{
diff --git a/docs/src/content/docs/guides/constellation-identification.mdx b/docs/src/content/docs/guides/constellation-identification.mdx
new file mode 100644
index 0000000..01c980a
--- /dev/null
+++ b/docs/src/content/docs/guides/constellation-identification.mdx
@@ -0,0 +1,206 @@
+---
+title: Constellation Identification
+sidebar:
+ order: 13
+---
+
+import { Steps, Aside, Tabs, TabItem } from "@astrojs/starlight/components";
+
+pg_orrery identifies which of the 88 IAU constellations contains a given sky position. You pass an equatorial coordinate (or raw RA and Dec values) and get back a three-letter IAU abbreviation. A companion function expands abbreviations to full names. The boundary lookup uses the Roman (1987) definitive boundary table, with internal precession from J2000 to the B1875.0 epoch in which the boundaries were originally defined.
+
+## How you do it today
+
+Determining which constellation an object is in usually involves:
+
+- **Stellarium**: Click on an object, read the constellation from the info panel. One object at a time, not scriptable.
+- **Astropy + regions**: Load the IAU constellation boundary polygons, precess coordinates to B1875.0, and run a point-in-polygon test. Correct, but requires assembling coordinate transforms and boundary data yourself.
+- **SIMBAD/CDS**: Query by object name and the constellation is in the metadata. Only works for cataloged objects, not arbitrary coordinates.
+- **Manual lookup**: Find the RA/Dec on a star chart and visually identify the constellation. Error-prone near boundaries.
+
+The constraint is always the same: the boundary data and the precession step live outside your database. If your star catalog, observation log, or scheduling system is in PostgreSQL, you export coordinates, look up constellations externally, and import the labels.
+
+## What changes with pg_orrery
+
+Three functions handle constellation identification:
+
+| Function | Returns | What it does |
+|---|---|---|
+| `constellation(equatorial)` | `text` (3-letter IAU code) | Identifies constellation from an equatorial coordinate |
+| `constellation(ra_hours, dec_deg)` | `text` (3-letter IAU code) | Convenience overload for raw J2000 RA (hours) and Dec (degrees) |
+| `constellation_full_name(abbrev)` | `text` (full name or NULL) | Expands a 3-letter abbreviation to the full IAU name |
+
+The first overload accepts the `equatorial` type returned by any of pg_orrery's equatorial functions -- `planet_equatorial()`, `sun_equatorial()`, `star_equatorial_pm()`, and so on. This makes the constellation lookup composable with the rest of the observation pipeline.
+
+The second overload takes raw RA in hours [0, 24) and Dec in degrees [-90, 90]. Use it for catalog data stored as individual columns.
+
+Both `constellation()` overloads precess the input J2000 coordinates to B1875.0 internally before testing against the ~357 boundary segments from the Roman (1987) table (CDS catalog VI/42). This precession is necessary because the IAU constellation boundaries were defined at epoch B1875.0.
+
+All three functions are `IMMUTABLE` -- the Roman boundary data is compiled into the extension, so there are no external dependencies. This means they are safe for use in indexes, generated columns, and materialized views.
+
+## What pg_orrery does not replace
+
+
+
+- **Not a substitute for detailed boundary analysis.** Objects within a few arcseconds of a constellation boundary may land on different sides depending on the precession model used. For critical applications (e.g., variable star designations), consult the CDS boundary data directly.
+- **No constellation figures or asterisms.** The function returns the IAU region that contains the coordinate, not any information about the traditional figure or its stars.
+- **88 constellations only.** The function returns the standard IAU three-letter abbreviation. Historical constellations (Argo Navis, Quadrans Muralis) are not represented.
+
+## Try it
+
+### Which constellation is Jupiter in?
+
+The simplest case -- combine `planet_equatorial()` with `constellation()`:
+
+```sql
+SELECT constellation(planet_equatorial(5, '2025-06-15 04:00:00+00')) AS jupiter_constellation;
+```
+
+This returns a three-letter IAU abbreviation like `Tau` or `Gem`, depending on where Jupiter is at the specified time.
+
+### Full composability chain
+
+Chain the equatorial, constellation, and full-name functions together to produce a readable result:
+
+```sql
+SELECT constellation_full_name(
+ constellation(
+ planet_equatorial(5, '2025-06-15 04:00:00+00')
+ )
+ ) AS jupiter_in;
+```
+
+This returns the full name -- something like `Taurus` or `Gemini`. The three functions nest cleanly because each takes the output of the previous one.
+
+### All planets and their constellations
+
+Sweep all seven visible planets at a given time:
+
+```sql
+SELECT 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,
+ constellation(planet_equatorial(body_id, '2025-06-15 04:00:00+00')) AS abbrev,
+ constellation_full_name(
+ constellation(planet_equatorial(body_id, '2025-06-15 04:00:00+00'))
+ ) AS constellation
+FROM generate_series(1, 8) AS body_id
+WHERE body_id != 3 -- cannot observe Earth from Earth
+ORDER BY body_id;
+```
+
+Seven rows, one query. Each planet gets its constellation abbreviation and full name.
+
+### Star catalog enrichment
+
+Add a constellation column to a catalog table using the raw-coordinate overload:
+
+```sql
+-- Existing star catalog
+CREATE TABLE star_catalog (
+ hip_id integer PRIMARY KEY,
+ name text,
+ ra_hours float8 NOT NULL,
+ dec_deg float8 NOT NULL,
+ vmag float8
+);
+
+INSERT INTO star_catalog VALUES
+ (11767, 'Polaris', 2.5303, 89.264, 1.98),
+ (32349, 'Sirius', 6.7525, -16.716, -1.46),
+ (27989, 'Betelgeuse', 5.9195, 7.407, 0.42),
+ (91262, 'Vega', 18.6156, 38.784, 0.03);
+
+-- Query with constellation identification
+SELECT name,
+ vmag,
+ constellation(ra_hours, dec_deg) AS iau_abbrev,
+ constellation_full_name(constellation(ra_hours, dec_deg)) AS constellation
+FROM star_catalog
+ORDER BY vmag;
+```
+
+Sirius returns `CMa` (Canis Major), Betelgeuse returns `Ori` (Orion), Vega returns `Lyr` (Lyra), and Polaris returns `UMi` (Ursa Minor). Because `constellation()` is `IMMUTABLE`, you could also store the result in a generated column:
+
+```sql
+ALTER TABLE star_catalog
+ ADD COLUMN iau_constellation text
+ GENERATED ALWAYS AS (constellation(ra_hours, dec_deg)) STORED;
+```
+
+### Sun through the zodiac
+
+The Sun's constellation changes roughly once a month as it moves along the ecliptic. Sample at monthly intervals to see the progression:
+
+```sql
+SELECT t::date AS date,
+ constellation(sun_equatorial(t)) AS abbrev,
+ constellation_full_name(constellation(sun_equatorial(t))) AS constellation
+FROM generate_series(
+ '2025-01-15'::timestamptz,
+ '2025-12-15'::timestamptz,
+ interval '1 month'
+) AS t;
+```
+
+The Sun passes through 13 constellations over the course of a year -- the 12 traditional zodiac constellations plus Ophiuchus, which the ecliptic crosses between Scorpius and Sagittarius. The IAU boundaries do not match the astrological 30-degree divisions, so the Sun spends significantly different amounts of time in each constellation.
+
+### What constellation is Polaris in?
+
+Using the `(ra_hours, dec_deg)` overload directly with known coordinates:
+
+```sql
+SELECT constellation(2.5303, 89.264) AS polaris_abbrev,
+ constellation_full_name(constellation(2.5303, 89.264)) AS polaris_constellation;
+```
+
+This returns `UMi` and `Ursa Minor`. The raw-coordinate overload is useful when you have RA and Dec values from an external source and do not need to construct an `equatorial` type first.
+
+### Full name display for UI
+
+A common pattern for user-facing output: show the abbreviation alongside the full name.
+
+```sql
+WITH stars(name, ra_h, dec_d) AS (VALUES
+ ('Polaris', 2.5303, 89.264),
+ ('Sirius', 6.7525, -16.716),
+ ('Betelgeuse', 5.9195, 7.407),
+ ('Vega', 18.6156, 38.784)
+)
+SELECT name,
+ constellation(ra_h, dec_d)
+ || ' (' || constellation_full_name(constellation(ra_h, dec_d)) || ')'
+ AS constellation_display
+FROM stars;
+```
+
+This produces strings like `CMa (Canis Major)` and `Ori (Orion)`. The `constellation_full_name()` function returns NULL for unrecognized abbreviations, so if you are working with external data, wrap the concatenation in a `COALESCE` or check for NULL:
+
+```sql
+SELECT COALESCE(
+ constellation_full_name('XYZ'),
+ 'Unknown'
+) AS result;
+-- Returns: Unknown
+```
+
+### Moon and Sun constellations over a week
+
+Track where the Moon and Sun are relative to the constellations:
+
+```sql
+SELECT t::date AS date,
+ constellation(sun_equatorial(t)) AS sun_in,
+ constellation(moon_equatorial(t)) AS moon_in
+FROM generate_series(
+ '2025-03-01'::timestamptz,
+ '2025-03-07'::timestamptz,
+ interval '1 day'
+) AS t;
+```
+
+The Moon moves roughly 13 degrees per day, crossing a constellation boundary every two to three days. The Sun barely moves -- it stays in the same constellation for the entire week.
diff --git a/docs/src/content/docs/guides/rise-set-prediction.mdx b/docs/src/content/docs/guides/rise-set-prediction.mdx
new file mode 100644
index 0000000..3ab20a8
--- /dev/null
+++ b/docs/src/content/docs/guides/rise-set-prediction.mdx
@@ -0,0 +1,309 @@
+---
+title: Rise/Set Prediction
+sidebar:
+ order: 12
+---
+
+import { Steps, Aside, Tabs, TabItem } from "@astrojs/starlight/components";
+
+pg_orrery computes rise and set times for the Sun, Moon, and all eight planets. You pass an observer and a starting timestamp, and get back a `timestamptz` for the next crossing of the horizon. Refracted variants account for atmospheric bending, matching what you actually see. Status diagnostics explain why a body might not cross the horizon at all.
+
+## How you do it today
+
+Finding when things rise and set involves a few approaches:
+
+- **Stellarium**: Shows rise/set times in the object info panel. One object at a time, not scriptable or queryable from your database.
+- **Skyfield**: Computes rise/set events by searching for horizon crossings. Well-designed API, but finding events for many bodies across many days means writing nested loops.
+- **Astropy + astroplan**: The `Observer` class computes rise/set/transit times. Handles refraction and horizon altitude. Per-object Python calls; batch queries over a table of targets need iteration.
+- **JPL Horizons**: Outputs rise/transit/set tables as part of its ephemeris service. One request per body, rate-limited, results live outside your database.
+
+The common pattern: compute rise/set times externally, then import the results into your scheduling table or observation log. If you want to answer "which planets are visible tonight?" in a single SQL query, you assemble the answer from pieces.
+
+## What changes with pg_orrery
+
+Rise and set times are SQL function calls. The functions search forward from a given timestamp using bisection (0.1-second precision) adapted from the satellite pass prediction algorithm. Three tiers cover different needs:
+
+| Tier | Functions | Threshold | Version |
+|---|---|---|---|
+| Geometric | `sun_next_rise`, `sun_next_set`, `moon_next_rise`, `moon_next_set`, `planet_next_rise`, `planet_next_set` | 0.0 deg | v0.13.0 |
+| Refracted | `sun_next_rise_refracted`, `sun_next_set_refracted`, `moon_next_rise_refracted`, `moon_next_set_refracted`, `planet_next_rise_refracted`, `planet_next_set_refracted` | varies | v0.14.0 |
+| Diagnostic | `sun_rise_set_status`, `moon_rise_set_status`, `planet_rise_set_status` | -- | v0.15.0 |
+
+All functions are `STABLE STRICT PARALLEL SAFE`. The search window is 7 days. If the body does not cross the threshold within that window, the function returns NULL.
+
+### Refraction thresholds
+
+Refracted functions use the same thresholds as the USNO and most almanacs:
+
+| Body | Threshold | Why |
+|---|---|---|
+| Sun | -0.833 deg | 34 arcmin atmospheric refraction + 16 arcmin solar semidiameter |
+| Moon | -0.833 deg | 34 arcmin refraction + ~15.5 arcmin mean lunar semidiameter |
+| Planets | -0.569 deg | 34 arcmin refraction only (point sources -- even Jupiter at opposition subtends just 0.4 arcmin) |
+
+The difference between geometric and refracted sunrise is typically 2-5 minutes. At extreme latitudes near the solstices, the difference can be much larger because the Sun's path intersects the horizon at a shallow angle.
+
+## What pg_orrery does not replace
+
+
+
+- **No topographic horizon.** All thresholds assume a flat, unobstructed horizon at sea level. Mountains, buildings, and terrain features are not considered.
+- **No civil/nautical/astronomical twilight.** The functions compute when the Sun's center crosses the specified threshold, not when it reaches -6, -12, or -18 degrees. You can approximate these by sampling `topo_elevation(sun_observe(...))` at the desired threshold.
+- **No lunar limb correction.** Moon rise/set uses a fixed mean semidiameter (15.5 arcmin). The actual semidiameter varies by about 1.5 arcmin between perigee and apogee, introducing up to ~15 seconds of timing error.
+- **No UT1-UTC correction.** Times are in UTC. The difference from UT1 (which governs the actual rotation of the Earth) is at most 0.9 seconds.
+
+## Try it
+
+### Basic sunrise and sunset
+
+The simplest rise/set query. Eagle, Idaho on the 2024 summer solstice:
+
+```sql
+SELECT sun_next_rise('43.7N 116.4W 800m'::observer,
+ '2024-06-21 12:00:00+00') AS sunrise,
+ sun_next_set('43.7N 116.4W 800m'::observer,
+ '2024-06-21 12:00:00+00') AS sunset;
+```
+
+The observer format is `lat lon altitude` -- latitude as degrees with N/S suffix, longitude with E/W suffix, altitude in meters. The start time should be before the expected event. Starting at noon UTC (6am MDT) catches the next sunrise and sunset for a Mountain Time observer.
+
+### Geometric vs refracted
+
+Atmospheric refraction lifts the Sun's apparent position near the horizon. Refracted sunrise is earlier; refracted sunset is later:
+
+
+
+ ```sql
+ SELECT sun_next_rise('43.7N 116.4W 800m'::observer,
+ '2024-06-21 12:00:00+00') AS geometric_rise,
+ sun_next_rise_refracted('43.7N 116.4W 800m'::observer,
+ '2024-06-21 12:00:00+00') AS refracted_rise,
+ sun_next_rise('43.7N 116.4W 800m'::observer,
+ '2024-06-21 12:00:00+00')
+ - sun_next_rise_refracted('43.7N 116.4W 800m'::observer,
+ '2024-06-21 12:00:00+00') AS difference;
+ ```
+
+ The difference is typically 2-5 minutes for mid-latitude locations. This is why newspaper sunrise times differ from the geometric horizon crossing.
+
+
+ ```sql
+ SELECT moon_next_rise('43.7N 116.4W 800m'::observer,
+ '2024-06-21 12:00:00+00') AS geometric_rise,
+ moon_next_rise_refracted('43.7N 116.4W 800m'::observer,
+ '2024-06-21 12:00:00+00') AS refracted_rise,
+ moon_next_rise('43.7N 116.4W 800m'::observer,
+ '2024-06-21 12:00:00+00')
+ - moon_next_rise_refracted('43.7N 116.4W 800m'::observer,
+ '2024-06-21 12:00:00+00') AS difference;
+ ```
+
+ The Moon uses the same -0.833 deg threshold as the Sun because its mean semidiameter is nearly identical.
+
+
+ ```sql
+ SELECT planet_next_rise(5, '43.7N 116.4W 800m'::observer,
+ '2024-06-21 12:00:00+00') AS geometric_rise,
+ planet_next_rise_refracted(5, '43.7N 116.4W 800m'::observer,
+ '2024-06-21 12:00:00+00') AS refracted_rise,
+ planet_next_rise(5, '43.7N 116.4W 800m'::observer,
+ '2024-06-21 12:00:00+00')
+ - planet_next_rise_refracted(5, '43.7N 116.4W 800m'::observer,
+ '2024-06-21 12:00:00+00') AS difference;
+ ```
+
+ Planets use a shallower threshold (-0.569 deg) because they are point sources -- no semidiameter to add.
+
+
+
+### What is visible tonight?
+
+Sweep all planets and check which ones rise before midnight local time:
+
+```sql
+WITH obs AS (SELECT '43.7N 116.4W 800m'::observer AS o),
+ t0 AS (SELECT '2024-06-21 12:00:00+00'::timestamptz AS t)
+SELECT 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,
+ planet_next_rise_refracted(body_id, obs.o, t0.t) AS rises,
+ planet_next_set_refracted(body_id, obs.o, t0.t) AS sets
+FROM generate_series(1, 8) AS body_id,
+ obs, t0
+WHERE body_id != 3 -- cannot observe Earth from Earth
+ AND planet_next_rise_refracted(body_id, obs.o, t0.t) IS NOT NULL
+ORDER BY planet_next_rise_refracted(body_id, obs.o, t0.t);
+```
+
+Body IDs follow the VSOP87 convention: 1=Mercury through 8=Neptune. Body 3 (Earth) and body 0 (Sun) are invalid for `planet_next_rise` and will raise an error.
+
+### Extreme latitude: midnight sun
+
+At 70 degrees north during the June solstice, the Sun never sets. Both rise/set functions express this by returning NULL:
+
+```sql
+SELECT sun_next_rise('70.0N 25.0E 10m'::observer,
+ '2024-06-21 12:00:00+00') AS sunrise,
+ sun_next_set('70.0N 25.0E 10m'::observer,
+ '2024-06-21 12:00:00+00') AS sunset,
+ sun_rise_set_status('70.0N 25.0E 10m'::observer,
+ '2024-06-21 12:00:00+00') AS status;
+```
+
+The `sunrise` column will have a value (the Sun is up and will "rise" again after the brief polar dip near midnight), but `sunset` returns NULL. The status function returns `'circumpolar'`, explaining why.
+
+
+
+### Polar night
+
+The inverse case. At 70 degrees north during the December solstice, the Sun never rises:
+
+```sql
+SELECT sun_next_rise('70.0N 25.0E 10m'::observer,
+ '2024-12-21 12:00:00+00') AS sunrise,
+ sun_next_set('70.0N 25.0E 10m'::observer,
+ '2024-12-21 12:00:00+00') AS sunset,
+ sun_rise_set_status('70.0N 25.0E 10m'::observer,
+ '2024-12-21 12:00:00+00') AS status;
+```
+
+Here `sunrise` is NULL, `sunset` may also be NULL (Sun is already below the horizon), and status returns `'never_rises'`.
+
+### Observation window planning
+
+How long can you observe Jupiter tonight? Compute rise-to-set duration:
+
+```sql
+WITH obs AS (SELECT '43.7N 116.4W 800m'::observer AS o),
+ t0 AS (SELECT '2024-06-21 12:00:00+00'::timestamptz AS t),
+ events AS (
+ SELECT planet_next_rise_refracted(5, obs.o, t0.t) AS jupiter_rise,
+ planet_next_set_refracted(5, obs.o, t0.t) AS jupiter_set
+ FROM obs, t0
+ )
+SELECT jupiter_rise,
+ jupiter_set,
+ jupiter_set - jupiter_rise AS observation_window
+FROM events
+WHERE jupiter_rise IS NOT NULL
+ AND jupiter_set IS NOT NULL
+ AND jupiter_set > jupiter_rise;
+```
+
+The `WHERE jupiter_set > jupiter_rise` clause handles the case where the body is already above the horizon -- the next set comes before the next rise. In that case, you would swap the logic: observe from now until the next set.
+
+### Moonrise for the next week
+
+Generate a daily moonrise table using `generate_series`:
+
+```sql
+SELECT day::date AS date,
+ moon_next_rise_refracted('43.7N 116.4W 800m'::observer, day) AS moonrise
+FROM generate_series(
+ '2024-06-21 12:00:00+00'::timestamptz,
+ '2024-06-28 12:00:00+00'::timestamptz,
+ interval '1 day'
+) AS day;
+```
+
+The Moon's orbital period is about 24 hours and 50 minutes, so moonrise drifts later by roughly 50 minutes each day. Some days may show NULL if the Moon does not rise during the search window starting from noon UTC.
+
+### All rise/set events for one night
+
+Combine Sun, Moon, and all visible planets into a single timeline:
+
+```sql
+WITH obs AS (SELECT '43.7N 116.4W 800m'::observer AS o),
+ t0 AS (SELECT '2024-06-21 12:00:00+00'::timestamptz AS t),
+ events AS (
+ SELECT 'Sun' AS body, 'rise' AS event,
+ sun_next_rise_refracted(obs.o, t0.t) AS time
+ FROM obs, t0
+ UNION ALL
+ SELECT 'Sun', 'set',
+ sun_next_set_refracted(obs.o, t0.t)
+ FROM obs, t0
+ UNION ALL
+ SELECT 'Moon', 'rise',
+ moon_next_rise_refracted(obs.o, t0.t)
+ FROM obs, t0
+ UNION ALL
+ SELECT 'Moon', 'set',
+ moon_next_set_refracted(obs.o, t0.t)
+ FROM obs, t0
+ UNION ALL
+ SELECT 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,
+ 'rise',
+ planet_next_rise_refracted(body_id, obs.o, t0.t)
+ FROM generate_series(1, 8) AS body_id, obs, t0
+ WHERE body_id != 3
+ UNION ALL
+ SELECT 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,
+ 'set',
+ planet_next_set_refracted(body_id, obs.o, t0.t)
+ FROM generate_series(1, 8) AS body_id, obs, t0
+ WHERE body_id != 3
+ )
+SELECT body, event, time
+FROM events
+WHERE time IS NOT NULL
+ORDER BY time;
+```
+
+This produces a chronological timeline of every rise and set event for the Sun, Moon, and all seven visible planets from Eagle, Idaho. NULL events (circumpolar or never-rising bodies) are filtered out.
+
+### Diagnosing NULL results across all bodies
+
+When planning observations at extreme latitudes, check every body's status at once:
+
+```sql
+WITH obs AS (SELECT '70.0N 25.0E 10m'::observer AS o),
+ t0 AS (SELECT '2024-06-21 12:00:00+00'::timestamptz AS t)
+SELECT 'Sun' AS body,
+ sun_rise_set_status(obs.o, t0.t) AS status
+FROM obs, t0
+
+UNION ALL
+
+SELECT 'Moon',
+ moon_rise_set_status(obs.o, t0.t)
+FROM obs, t0
+
+UNION ALL
+
+SELECT 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,
+ planet_rise_set_status(body_id, obs.o, t0.t)
+FROM generate_series(1, 8) AS body_id, obs, t0
+WHERE body_id != 3
+ORDER BY body;
+```
+
+At 70 degrees north on the summer solstice, the Sun is circumpolar. Some planets may also be circumpolar or never-rising depending on their current declination. The status functions classify each body with a single 24-hour scan (48 samples at 30-minute spacing), so this query is lightweight.
diff --git a/pg_orrery.control b/pg_orrery.control
index fa3fe02..9b863f8 100644
--- a/pg_orrery.control
+++ b/pg_orrery.control
@@ -1,4 +1,4 @@
comment = 'A database orrery — celestial mechanics types and functions for PostgreSQL'
-default_version = '0.15.0'
+default_version = '0.16.0'
module_pathname = '$libdir/pg_orrery'
relocatable = true
diff --git a/sql/pg_orrery--0.15.0--0.16.0.sql b/sql/pg_orrery--0.15.0--0.16.0.sql
new file mode 100644
index 0000000..5224cd6
--- /dev/null
+++ b/sql/pg_orrery--0.15.0--0.16.0.sql
@@ -0,0 +1,79 @@
+-- pg_orrery 0.15.0 -> 0.16.0: twilight, lunar phase, planet magnitude
+
+-- ============================================================
+-- Twilight functions (6)
+-- ============================================================
+
+CREATE FUNCTION sun_civil_dawn(observer, timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME', 'sun_civil_dawn'
+ LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_civil_dawn(observer, timestamptz) IS
+ 'Next civil dawn (Sun crosses -6 deg rising). Outdoor activities without artificial light.';
+
+CREATE FUNCTION sun_civil_dusk(observer, timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME', 'sun_civil_dusk'
+ LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_civil_dusk(observer, timestamptz) IS
+ 'Next civil dusk (Sun crosses -6 deg setting). Artificial light needed.';
+
+CREATE FUNCTION sun_nautical_dawn(observer, timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME', 'sun_nautical_dawn'
+ LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_nautical_dawn(observer, timestamptz) IS
+ 'Next nautical dawn (Sun crosses -12 deg rising). Horizon visible at sea.';
+
+CREATE FUNCTION sun_nautical_dusk(observer, timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME', 'sun_nautical_dusk'
+ LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_nautical_dusk(observer, timestamptz) IS
+ 'Next nautical dusk (Sun crosses -12 deg setting). Horizon no longer visible at sea.';
+
+CREATE FUNCTION sun_astronomical_dawn(observer, timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME', 'sun_astronomical_dawn'
+ LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_astronomical_dawn(observer, timestamptz) IS
+ 'Next astronomical dawn (Sun crosses -18 deg rising). Sky was fully dark.';
+
+CREATE FUNCTION sun_astronomical_dusk(observer, timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME', 'sun_astronomical_dusk'
+ LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_astronomical_dusk(observer, timestamptz) IS
+ 'Next astronomical dusk (Sun crosses -18 deg setting). Sky becomes fully dark.';
+
+-- ============================================================
+-- Lunar phase functions (4)
+-- ============================================================
+
+CREATE FUNCTION moon_phase_angle(timestamptz) RETURNS float8
+ AS 'MODULE_PATHNAME', 'moon_phase_angle'
+ LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION moon_phase_angle(timestamptz) IS
+ 'Sun-Earth-Moon phase angle in degrees [0,360). 0=new, 90=first quarter, 180=full, 270=last quarter.';
+
+CREATE FUNCTION moon_illumination(timestamptz) RETURNS float8
+ AS 'MODULE_PATHNAME', 'moon_illumination'
+ LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION moon_illumination(timestamptz) IS
+ 'Illuminated fraction of the Moon disk [0.0, 1.0].';
+
+CREATE FUNCTION moon_phase_name(timestamptz) RETURNS text
+ AS 'MODULE_PATHNAME', 'moon_phase_name'
+ LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION moon_phase_name(timestamptz) IS
+ 'Moon phase name: new_moon, waxing_crescent, first_quarter, waxing_gibbous, full_moon, waning_gibbous, last_quarter, waning_crescent.';
+
+CREATE FUNCTION moon_age(timestamptz) RETURNS float8
+ AS 'MODULE_PATHNAME', 'moon_age'
+ LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION moon_age(timestamptz) IS
+ 'Days since last new moon [0, ~29.53), approximated from phase angle.';
+
+-- ============================================================
+-- Planet magnitude (1)
+-- ============================================================
+
+CREATE FUNCTION planet_magnitude(int4, timestamptz) RETURNS float8
+ AS 'MODULE_PATHNAME', 'planet_magnitude'
+ LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION planet_magnitude(int4, timestamptz) IS
+ 'Apparent visual magnitude of a planet (Mallama & Hilton 2018). Body IDs 1-8. Saturn ring tilt not modeled.';
diff --git a/sql/pg_orrery--0.16.0.sql b/sql/pg_orrery--0.16.0.sql
new file mode 100644
index 0000000..23d1b70
--- /dev/null
+++ b/sql/pg_orrery--0.16.0.sql
@@ -0,0 +1,1674 @@
+-- pg_orrery -- Orbital mechanics types and functions for PostgreSQL
+--
+-- Types: tle, eci_position, geodetic, topocentric, observer, pass_event
+-- Provides SGP4/SDP4 propagation, coordinate transforms, pass prediction,
+-- and GiST indexing on altitude bands for conjunction screening.
+--
+-- All propagation uses WGS-72 constants (matching TLE mean element fitting).
+-- Coordinate output uses WGS-84 (matching modern geodetic standards).
+
+-- ============================================================
+-- Shell types (forward declarations)
+-- ============================================================
+
+CREATE TYPE tle;
+CREATE TYPE eci_position;
+CREATE TYPE geodetic;
+CREATE TYPE topocentric;
+CREATE TYPE observer;
+CREATE TYPE pass_event;
+
+
+-- ============================================================
+-- TLE type: Two-Line Element set
+-- ============================================================
+
+CREATE FUNCTION tle_in(cstring) RETURNS tle
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION tle_out(tle) RETURNS cstring
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION tle_recv(internal) RETURNS tle
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION tle_send(tle) RETURNS bytea
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE TYPE tle (
+ INPUT = tle_in,
+ OUTPUT = tle_out,
+ RECEIVE = tle_recv,
+ SEND = tle_send,
+ INTERNALLENGTH = 112,
+ ALIGNMENT = double,
+ STORAGE = plain
+);
+
+COMMENT ON TYPE tle IS 'Two-Line Element set — parsed mean orbital elements for SGP4/SDP4 propagation';
+
+-- TLE accessor functions
+
+CREATE FUNCTION tle_epoch(tle) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION tle_epoch(tle) IS 'TLE epoch as Julian date (UTC)';
+
+CREATE FUNCTION tle_norad_id(tle) RETURNS int4
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION tle_norad_id(tle) IS 'NORAD catalog number';
+
+CREATE FUNCTION tle_inclination(tle) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION tle_inclination(tle) IS 'Orbital inclination in degrees';
+
+CREATE FUNCTION tle_eccentricity(tle) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION tle_eccentricity(tle) IS 'Orbital eccentricity (dimensionless)';
+
+CREATE FUNCTION tle_raan(tle) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION tle_raan(tle) IS 'Right ascension of ascending node in degrees';
+
+CREATE FUNCTION tle_arg_perigee(tle) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION tle_arg_perigee(tle) IS 'Argument of perigee in degrees';
+
+CREATE FUNCTION tle_mean_anomaly(tle) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION tle_mean_anomaly(tle) IS 'Mean anomaly in degrees';
+
+CREATE FUNCTION tle_mean_motion(tle) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION tle_mean_motion(tle) IS 'Mean motion in revolutions per day';
+
+CREATE FUNCTION tle_bstar(tle) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION tle_bstar(tle) IS 'B* drag coefficient (1/earth-radii)';
+
+CREATE FUNCTION tle_period(tle) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION tle_period(tle) IS 'Orbital period in minutes';
+
+CREATE FUNCTION tle_age(tle, timestamptz) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION tle_age(tle, timestamptz) IS 'TLE age in days (positive = past epoch, negative = before epoch)';
+
+CREATE FUNCTION tle_perigee(tle) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION tle_perigee(tle) IS 'Perigee altitude in km above WGS-72 ellipsoid';
+
+CREATE FUNCTION tle_apogee(tle) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION tle_apogee(tle) IS 'Apogee altitude in km above WGS-72 ellipsoid';
+
+CREATE FUNCTION tle_intl_desig(tle) RETURNS text
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION tle_intl_desig(tle) IS 'International designator (COSPAR ID)';
+
+CREATE FUNCTION tle_from_lines(text, text) RETURNS tle
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION tle_from_lines(text, text) IS
+ 'Construct TLE from separate line1/line2 text columns';
+
+
+-- ============================================================
+-- ECI position type: True Equator Mean Equinox (TEME) frame
+-- ============================================================
+
+CREATE FUNCTION eci_in(cstring) RETURNS eci_position
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION eci_out(eci_position) RETURNS cstring
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION eci_recv(internal) RETURNS eci_position
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION eci_send(eci_position) RETURNS bytea
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE TYPE eci_position (
+ INPUT = eci_in,
+ OUTPUT = eci_out,
+ RECEIVE = eci_recv,
+ SEND = eci_send,
+ INTERNALLENGTH = 48,
+ ALIGNMENT = double,
+ STORAGE = plain
+);
+
+COMMENT ON TYPE eci_position IS 'Earth-Centered Inertial position and velocity in TEME frame (km, km/s)';
+
+-- ECI accessor functions
+
+CREATE FUNCTION eci_x(eci_position) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION eci_y(eci_position) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION eci_z(eci_position) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION eci_vx(eci_position) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION eci_vy(eci_position) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION eci_vz(eci_position) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION eci_speed(eci_position) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION eci_speed(eci_position) IS 'Velocity magnitude in km/s';
+
+CREATE FUNCTION eci_altitude(eci_position) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION eci_altitude(eci_position) IS 'Approximate geocentric altitude in km (radius - WGS72_AE)';
+
+
+-- ============================================================
+-- Geodetic type: WGS-84 latitude/longitude/altitude
+-- ============================================================
+
+CREATE FUNCTION geodetic_in(cstring) RETURNS geodetic
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION geodetic_out(geodetic) RETURNS cstring
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION geodetic_recv(internal) RETURNS geodetic
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION geodetic_send(geodetic) RETURNS bytea
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE TYPE geodetic (
+ INPUT = geodetic_in,
+ OUTPUT = geodetic_out,
+ RECEIVE = geodetic_recv,
+ SEND = geodetic_send,
+ INTERNALLENGTH = 24,
+ ALIGNMENT = double,
+ STORAGE = plain
+);
+
+COMMENT ON TYPE geodetic IS 'Geodetic coordinates on WGS-84 ellipsoid (lat/lon in degrees, altitude in km)';
+
+CREATE FUNCTION geodetic_lat(geodetic) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION geodetic_lon(geodetic) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION geodetic_alt(geodetic) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+
+-- ============================================================
+-- Topocentric type: observer-relative az/el/range
+-- ============================================================
+
+CREATE FUNCTION topocentric_in(cstring) RETURNS topocentric
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION topocentric_out(topocentric) RETURNS cstring
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION topocentric_recv(internal) RETURNS topocentric
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION topocentric_send(topocentric) RETURNS bytea
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE TYPE topocentric (
+ INPUT = topocentric_in,
+ OUTPUT = topocentric_out,
+ RECEIVE = topocentric_recv,
+ SEND = topocentric_send,
+ INTERNALLENGTH = 32,
+ ALIGNMENT = double,
+ STORAGE = plain
+);
+
+COMMENT ON TYPE topocentric IS 'Topocentric coordinates relative to observer (azimuth, elevation, range, range rate)';
+
+CREATE FUNCTION topo_azimuth(topocentric) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION topo_azimuth(topocentric) IS 'Azimuth in degrees (0=N, 90=E, 180=S, 270=W)';
+
+CREATE FUNCTION topo_elevation(topocentric) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION topo_elevation(topocentric) IS 'Elevation in degrees (0=horizon, 90=zenith)';
+
+CREATE FUNCTION topo_range(topocentric) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION topo_range(topocentric) IS 'Slant range in km';
+
+CREATE FUNCTION topo_range_rate(topocentric) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION topo_range_rate(topocentric) IS 'Range rate in km/s (positive = receding)';
+
+
+-- ============================================================
+-- Observer type: ground station location
+-- ============================================================
+
+CREATE FUNCTION observer_in(cstring) RETURNS observer
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION observer_out(observer) RETURNS cstring
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION observer_recv(internal) RETURNS observer
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION observer_send(observer) RETURNS bytea
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE TYPE observer (
+ INPUT = observer_in,
+ OUTPUT = observer_out,
+ RECEIVE = observer_recv,
+ SEND = observer_send,
+ INTERNALLENGTH = 24,
+ ALIGNMENT = double,
+ STORAGE = plain
+);
+
+COMMENT ON TYPE observer IS 'Ground observer location (accepts: 40.0N 105.3W 1655m or decimal degrees)';
+
+CREATE FUNCTION observer_lat(observer) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION observer_lat(observer) IS 'Latitude in degrees (positive = North)';
+
+CREATE FUNCTION observer_lon(observer) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION observer_lon(observer) IS 'Longitude in degrees (positive = East)';
+
+CREATE FUNCTION observer_alt(observer) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION observer_alt(observer) IS 'Altitude in meters above WGS-84 ellipsoid';
+
+CREATE FUNCTION observer_from_geodetic(float8, float8, float8 DEFAULT 0.0) RETURNS observer
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION observer_from_geodetic(float8, float8, float8) IS
+ 'Construct observer from lat (deg), lon (deg), altitude (meters). Avoids text formatting round-trips.';
+
+
+-- ============================================================
+-- Pass event type: satellite visibility window
+-- ============================================================
+
+CREATE FUNCTION pass_event_in(cstring) RETURNS pass_event
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION pass_event_out(pass_event) RETURNS cstring
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION pass_event_recv(internal) RETURNS pass_event
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION pass_event_send(pass_event) RETURNS bytea
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE TYPE pass_event (
+ INPUT = pass_event_in,
+ OUTPUT = pass_event_out,
+ RECEIVE = pass_event_recv,
+ SEND = pass_event_send,
+ INTERNALLENGTH = 48,
+ ALIGNMENT = double,
+ STORAGE = plain
+);
+
+COMMENT ON TYPE pass_event IS 'Satellite pass event (AOS/MAX/LOS times, max elevation, AOS/LOS azimuths)';
+
+CREATE FUNCTION pass_aos_time(pass_event) RETURNS timestamptz
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION pass_aos_time(pass_event) IS 'Acquisition of signal time';
+
+CREATE FUNCTION pass_max_el_time(pass_event) RETURNS timestamptz
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION pass_max_el_time(pass_event) IS 'Maximum elevation time';
+
+CREATE FUNCTION pass_los_time(pass_event) RETURNS timestamptz
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION pass_los_time(pass_event) IS 'Loss of signal time';
+
+CREATE FUNCTION pass_max_elevation(pass_event) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION pass_max_elevation(pass_event) IS 'Maximum elevation in degrees';
+
+CREATE FUNCTION pass_aos_azimuth(pass_event) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION pass_aos_azimuth(pass_event) IS 'AOS azimuth in degrees (0=N)';
+
+CREATE FUNCTION pass_los_azimuth(pass_event) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION pass_los_azimuth(pass_event) IS 'LOS azimuth in degrees (0=N)';
+
+CREATE FUNCTION pass_duration(pass_event) RETURNS interval
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION pass_duration(pass_event) IS 'Pass duration (LOS - AOS)';
+
+
+-- ============================================================
+-- SGP4/SDP4 propagation functions
+-- ============================================================
+
+CREATE FUNCTION sgp4_propagate(tle, timestamptz) RETURNS eci_position
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sgp4_propagate(tle, timestamptz) IS
+ 'Propagate TLE to a point in time using SGP4 (near-earth) or SDP4 (deep-space). Returns TEME ECI position/velocity.';
+
+CREATE FUNCTION sgp4_propagate_safe(tle, timestamptz) RETURNS eci_position
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE PARALLEL SAFE;
+COMMENT ON FUNCTION sgp4_propagate_safe(tle, timestamptz) IS
+ 'Like sgp4_propagate but returns NULL on error instead of raising an exception. For batch queries with potentially invalid TLEs.';
+
+CREATE FUNCTION sgp4_propagate_series(tle, timestamptz, timestamptz, interval)
+ RETURNS TABLE(t timestamptz, x float8, y float8, z float8, vx float8, vy float8, vz float8)
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE
+ ROWS 100;
+COMMENT ON FUNCTION sgp4_propagate_series(tle, timestamptz, timestamptz, interval) IS
+ 'Propagate TLE over a time range at regular intervals. Returns time series of TEME positions.';
+
+CREATE FUNCTION tle_distance(tle, tle, timestamptz) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION tle_distance(tle, tle, timestamptz) IS
+ 'Euclidean distance in km between two TLEs at a reference time';
+
+
+-- ============================================================
+-- Coordinate transform functions
+-- ============================================================
+
+CREATE FUNCTION eci_to_geodetic(eci_position, timestamptz) RETURNS geodetic
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION eci_to_geodetic(eci_position, timestamptz) IS
+ 'Convert TEME ECI position to WGS-84 geodetic coordinates at given time';
+
+CREATE FUNCTION eci_to_topocentric(eci_position, observer, timestamptz) RETURNS topocentric
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION eci_to_topocentric(eci_position, observer, timestamptz) IS
+ 'Convert TEME ECI position to topocentric (az/el/range) relative to observer';
+
+CREATE FUNCTION subsatellite_point(tle, timestamptz) RETURNS geodetic
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION subsatellite_point(tle, timestamptz) IS
+ 'Subsatellite (nadir) point on WGS-84 ellipsoid at given time';
+
+CREATE FUNCTION ground_track(tle, timestamptz, timestamptz, interval)
+ RETURNS TABLE(t timestamptz, lat float8, lon float8, alt float8)
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE
+ ROWS 100;
+COMMENT ON FUNCTION ground_track(tle, timestamptz, timestamptz, interval) IS
+ 'Ground track as time series of subsatellite points (lat/lon in degrees, alt in km)';
+
+CREATE FUNCTION observe(tle, observer, timestamptz) RETURNS topocentric
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION observe(tle, observer, timestamptz) IS
+ 'Propagate TLE and compute observer-relative look angles in one call. Returns topocentric (az/el/range/range_rate).';
+
+CREATE FUNCTION observe_safe(tle, observer, timestamptz) RETURNS topocentric
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE PARALLEL SAFE;
+COMMENT ON FUNCTION observe_safe(tle, observer, timestamptz) IS
+ 'Like observe() but returns NULL on propagation error. For batch queries with potentially invalid/decayed TLEs.';
+
+
+-- ============================================================
+-- Pass prediction functions
+-- ============================================================
+
+CREATE FUNCTION next_pass(tle, observer, timestamptz) RETURNS pass_event
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION next_pass(tle, observer, timestamptz) IS
+ 'Find the next satellite pass over observer (searches up to 7 days ahead)';
+
+CREATE FUNCTION predict_passes(tle, observer, timestamptz, timestamptz, float8 DEFAULT 0.0)
+ RETURNS SETOF pass_event
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE
+ ROWS 10;
+COMMENT ON FUNCTION predict_passes(tle, observer, timestamptz, timestamptz, float8) IS
+ 'Predict all satellite passes over observer in time window. Optional min_elevation in degrees.';
+
+CREATE FUNCTION pass_visible(tle, observer, timestamptz, timestamptz) RETURNS boolean
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION pass_visible(tle, observer, timestamptz, timestamptz) IS
+ 'True if any pass occurs over observer in the time window';
+
+
+-- ============================================================
+-- GiST operator support functions
+-- ============================================================
+
+-- Overlap operator: do orbital keys overlap in altitude AND inclination?
+CREATE FUNCTION tle_overlap(tle, tle) RETURNS boolean
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+-- Altitude distance operator (altitude-only, for KNN ordering)
+CREATE FUNCTION tle_alt_distance(tle, tle) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE OPERATOR && (
+ LEFTARG = tle,
+ RIGHTARG = tle,
+ FUNCTION = tle_overlap,
+ COMMUTATOR = &&,
+ RESTRICT = areasel,
+ JOIN = areajoinsel
+);
+
+COMMENT ON OPERATOR && (tle, tle) IS 'Orbital key overlap (altitude band AND inclination range) — necessary condition for conjunction';
+
+CREATE OPERATOR <-> (
+ LEFTARG = tle,
+ RIGHTARG = tle,
+ FUNCTION = tle_alt_distance,
+ COMMUTATOR = <->
+);
+
+COMMENT ON OPERATOR <-> (tle, tle) IS '2-D orbital distance in km: L2 norm of altitude-band gap and inclination gap (radians × Earth radius). Returns 0 when both dimensions overlap.';
+
+
+-- ============================================================
+-- 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);
+
+
+-- ============================================================
+-- 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.';
+
+
+-- ============================================================
+-- 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.';
+
+
+-- ============================================================
+-- 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.';
+
+
+-- ============================================================
+-- 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.';
+
+
+-- ============================================================
+-- 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.';
+
+
+-- ============================================================
+-- DE ephemeris functions (optional high-precision)
+-- ============================================================
+
+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.';
+
+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.';
+
+CREATE FUNCTION galilean_observe_de(int4, observer, timestamptz) RETURNS topocentric
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION galilean_observe_de(int4, observer, timestamptz) IS
+ 'Observe Galilean moon with JPL DE parent position. L1.2 moon offsets. Body IDs: 0-3 (Io-Callisto).';
+
+CREATE FUNCTION saturn_moon_observe_de(int4, observer, timestamptz) RETURNS topocentric
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION saturn_moon_observe_de(int4, observer, timestamptz) IS
+ 'Observe Saturn moon with JPL DE parent position. TASS 1.7 moon offsets. Body IDs: 0-7 (Mimas-Hyperion).';
+
+CREATE FUNCTION uranus_moon_observe_de(int4, observer, timestamptz) RETURNS topocentric
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION uranus_moon_observe_de(int4, observer, timestamptz) IS
+ 'Observe Uranus moon with JPL DE parent position. GUST86 moon offsets. Body IDs: 0-4 (Miranda-Oberon).';
+
+CREATE FUNCTION mars_moon_observe_de(int4, observer, timestamptz) RETURNS topocentric
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION mars_moon_observe_de(int4, observer, timestamptz) IS
+ 'Observe Mars moon with JPL DE parent position. MarsSat moon offsets. Body IDs: 0-1 (Phobos-Deimos).';
+
+
+-- Diagnostic function
+
+CREATE FUNCTION pg_orrery_ephemeris_info(
+ OUT provider text, OUT file_path text,
+ OUT start_jd float8, OUT end_jd float8,
+ OUT version int4, OUT au_km float8
+) RETURNS RECORD
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE;
+COMMENT ON FUNCTION pg_orrery_ephemeris_info() IS
+ 'Returns current ephemeris provider status: VSOP87 or JPL_DE with file path, JD range, version, and AU value.';
+
+
+-- ============================================================
+-- Orbit determination (TLE fitting from observations)
+-- ============================================================
+
+-- Fit TLE from ECI position/velocity ephemeris
+-- weights: per-observation weighting (NULL = uniform)
+
+CREATE FUNCTION tle_from_eci(
+ positions eci_position[], times timestamptz[],
+ seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false,
+ max_iter int4 DEFAULT 15,
+ weights float8[] DEFAULT NULL,
+ OUT fitted_tle tle, OUT iterations int4,
+ OUT rms_final float8, OUT rms_initial float8, OUT status text,
+ OUT condition_number float8, OUT covariance float8[], OUT nstate int4
+) RETURNS RECORD
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE;
+COMMENT ON FUNCTION tle_from_eci(eci_position[], timestamptz[], tle, boolean, int4, float8[]) IS
+ 'Fit a TLE from ECI position/velocity observations via differential correction. Optional per-observation weights for heterogeneous sensor fusion. Returns fitted TLE, RMS residuals, convergence status, condition number, and formal covariance matrix.';
+
+-- Fit TLE from topocentric observations (az/el/range) — single observer
+-- fit_range_rate: include range_rate as 4th residual component
+-- weights: per-observation weighting (NULL = uniform)
+
+CREATE FUNCTION tle_from_topocentric(
+ observations topocentric[], times timestamptz[],
+ obs observer,
+ seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false,
+ max_iter int4 DEFAULT 15,
+ fit_range_rate boolean DEFAULT false,
+ weights float8[] DEFAULT NULL,
+ OUT fitted_tle tle, OUT iterations int4,
+ OUT rms_final float8, OUT rms_initial float8, OUT status text,
+ OUT condition_number float8, OUT covariance float8[], OUT nstate int4
+) RETURNS RECORD
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE;
+COMMENT ON FUNCTION tle_from_topocentric(topocentric[], timestamptz[], observer, tle, boolean, int4, boolean, float8[]) IS
+ 'Fit a TLE from topocentric (az/el/range) observations via differential correction. Optional range_rate fitting and per-observation weights. Returns fitted TLE, RMS residuals, convergence status, condition number, and formal covariance matrix.';
+
+-- Fit TLE from topocentric observations — multiple observers
+
+CREATE FUNCTION tle_from_topocentric(
+ observations topocentric[], times timestamptz[],
+ observers observer[], observer_ids int4[],
+ seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false,
+ max_iter int4 DEFAULT 15,
+ fit_range_rate boolean DEFAULT false,
+ weights float8[] DEFAULT NULL,
+ OUT fitted_tle tle, OUT iterations int4,
+ OUT rms_final float8, OUT rms_initial float8, OUT status text,
+ OUT condition_number float8, OUT covariance float8[], OUT nstate int4
+) RETURNS RECORD
+ AS 'MODULE_PATHNAME', 'tle_from_topocentric_multi'
+ LANGUAGE C STABLE PARALLEL SAFE;
+COMMENT ON FUNCTION tle_from_topocentric(topocentric[], timestamptz[], observer[], int4[], tle, boolean, int4, boolean, float8[]) IS
+ 'Fit a TLE from topocentric observations collected by multiple ground stations. observer_ids[i] indexes into observers[]. Optional range_rate fitting and per-observation weights.';
+
+-- Per-observation residuals diagnostic
+
+CREATE FUNCTION tle_fit_residuals(
+ fitted tle,
+ positions eci_position[],
+ times timestamptz[]
+) RETURNS TABLE (
+ t timestamptz,
+ dx_km float8,
+ dy_km float8,
+ dz_km float8,
+ pos_err_km float8
+)
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION tle_fit_residuals(tle, eci_position[], timestamptz[]) IS
+ 'Compute per-observation position residuals (km) between a TLE and ECI observations. Useful for fit quality assessment.';
+
+-- Fit TLE from RA/Dec observations — single observer
+-- Uses Gauss method for initial orbit determination when no seed is provided.
+-- RA in hours [0,24), Dec in degrees [-90,90] (matches star_observe convention).
+-- RMS output is in radians for angles-only mode.
+
+CREATE FUNCTION tle_from_angles(
+ ra_hours float8[], dec_degrees float8[],
+ times timestamptz[],
+ obs observer,
+ seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false,
+ max_iter int4 DEFAULT 15,
+ weights float8[] DEFAULT NULL,
+ OUT fitted_tle tle, OUT iterations int4,
+ OUT rms_final float8, OUT rms_initial float8, OUT status text,
+ OUT condition_number float8, OUT covariance float8[], OUT nstate int4
+) RETURNS RECORD AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE;
+COMMENT ON FUNCTION tle_from_angles(float8[], float8[], timestamptz[], observer, tle, boolean, int4, float8[]) IS
+ 'Fit a TLE from angles-only (RA/Dec) observations via Gauss IOD + differential correction. RA in hours [0,24), Dec in degrees [-90,90]. RMS output in radians. Uses Gauss method for seed-free initial guess.';
+
+-- Fit TLE from RA/Dec observations — multiple observers
+
+CREATE FUNCTION tle_from_angles(
+ ra_hours float8[], dec_degrees float8[],
+ times timestamptz[],
+ observers observer[], observer_ids int4[],
+ seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false,
+ max_iter int4 DEFAULT 15,
+ weights float8[] DEFAULT NULL,
+ OUT fitted_tle tle, OUT iterations int4,
+ OUT rms_final float8, OUT rms_initial float8, OUT status text,
+ OUT condition_number float8, OUT covariance float8[], OUT nstate int4
+) RETURNS RECORD
+ AS 'MODULE_PATHNAME', 'tle_from_angles_multi'
+ LANGUAGE C STABLE PARALLEL SAFE;
+COMMENT ON FUNCTION tle_from_angles(float8[], float8[], timestamptz[], observer[], int4[], tle, boolean, int4, float8[]) IS
+ 'Fit a TLE from angles-only (RA/Dec) observations from multiple ground stations via Gauss IOD + differential correction.';
+-- pg_orrery 0.6.0 -> 0.7.0 migration
+--
+-- Adds SP-GiST orbital trie index for satellite pass prediction.
+-- 2-level trie: SMA (L0) + inclination (L1) with query-time RAAN filter.
+-- The &? operator answers "might this satellite be visible?"
+
+-- ============================================================
+-- observer_window composite type (query parameter bundle)
+-- ============================================================
+
+CREATE TYPE observer_window AS (
+ obs observer,
+ t_start timestamptz,
+ t_end timestamptz,
+ min_el float8
+);
+
+COMMENT ON TYPE observer_window IS
+ 'Observation query parameters: observer location, time window, and minimum elevation angle (degrees). Used with the &? visibility cone operator.';
+
+-- ============================================================
+-- Visibility cone operator function
+-- ============================================================
+
+CREATE FUNCTION tle_visibility_possible(tle, observer_window) RETURNS boolean
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+
+COMMENT ON FUNCTION tle_visibility_possible(tle, observer_window) IS
+ 'Could this satellite be visible from the observer during the time window? Combines altitude, inclination, and RAAN checks. Conservative superset — survivors need SGP4 propagation for ground truth.';
+
+-- ============================================================
+-- &? operator (visibility cone check)
+-- ============================================================
+-- The indexed column (tle) MUST be the left argument so PostgreSQL
+-- can form a ScanKey and pass it to inner_consistent for pruning.
+
+CREATE OPERATOR &? (
+ LEFTARG = tle,
+ RIGHTARG = observer_window,
+ FUNCTION = tle_visibility_possible,
+ RESTRICT = contsel,
+ JOIN = contjoinsel
+);
+
+COMMENT ON OPERATOR &? (tle, observer_window) IS
+ 'Visibility cone check: could this satellite be visible from the observer during the time window? Index-accelerated via SP-GiST orbital trie.';
+
+-- ============================================================
+-- SP-GiST support functions
+-- ============================================================
+
+CREATE FUNCTION spgist_tle_config(internal, internal) RETURNS void
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION spgist_tle_choose(internal, internal) RETURNS void
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION spgist_tle_picksplit(internal, internal) RETURNS void
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION spgist_tle_inner_consistent(internal, internal) RETURNS void
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION spgist_tle_leaf_consistent(internal, internal) RETURNS void
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+-- ============================================================
+-- SP-GiST operator class (opt-in, not DEFAULT)
+-- ============================================================
+
+CREATE OPERATOR CLASS tle_spgist_ops
+ FOR TYPE tle USING spgist AS
+ OPERATOR 1 &? (tle, observer_window),
+ FUNCTION 1 spgist_tle_config(internal, internal),
+ FUNCTION 2 spgist_tle_choose(internal, internal),
+ FUNCTION 3 spgist_tle_picksplit(internal, internal),
+ FUNCTION 4 spgist_tle_inner_consistent(internal, internal),
+ FUNCTION 5 spgist_tle_leaf_consistent(internal, internal);
+-- pg_orrery 0.7.0 -> 0.8.0 migration
+--
+-- Adds orbital_elements type for comets/asteroids, MPC MPCORB.DAT parser,
+-- and small_body_observe()/small_body_heliocentric() observation functions.
+
+-- ============================================================
+-- orbital_elements type
+-- ============================================================
+
+CREATE TYPE orbital_elements;
+
+CREATE FUNCTION orbital_elements_in(cstring) RETURNS orbital_elements
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION orbital_elements_out(orbital_elements) RETURNS cstring
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION orbital_elements_recv(internal) RETURNS orbital_elements
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION orbital_elements_send(orbital_elements) RETURNS bytea
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE TYPE orbital_elements (
+ INPUT = orbital_elements_in,
+ OUTPUT = orbital_elements_out,
+ RECEIVE = orbital_elements_recv,
+ SEND = orbital_elements_send,
+ INTERNALLENGTH = 72,
+ ALIGNMENT = double,
+ STORAGE = plain
+);
+
+COMMENT ON TYPE orbital_elements IS
+ 'Classical Keplerian orbital elements for comets and asteroids (epoch, q, e, inc, omega, Omega, tp, H, G). 72 bytes, fixed-size.';
+
+
+-- ============================================================
+-- Accessor functions
+-- ============================================================
+
+CREATE FUNCTION oe_epoch(orbital_elements) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION oe_epoch(orbital_elements) IS 'Osculation epoch (Julian date)';
+
+CREATE FUNCTION oe_perihelion(orbital_elements) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION oe_perihelion(orbital_elements) IS 'Perihelion distance q (AU)';
+
+CREATE FUNCTION oe_eccentricity(orbital_elements) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION oe_eccentricity(orbital_elements) IS 'Eccentricity';
+
+CREATE FUNCTION oe_inclination(orbital_elements) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION oe_inclination(orbital_elements) IS 'Inclination (degrees)';
+
+CREATE FUNCTION oe_arg_perihelion(orbital_elements) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION oe_arg_perihelion(orbital_elements) IS 'Argument of perihelion (degrees)';
+
+CREATE FUNCTION oe_raan(orbital_elements) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION oe_raan(orbital_elements) IS 'Longitude of ascending node (degrees)';
+
+CREATE FUNCTION oe_tp(orbital_elements) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION oe_tp(orbital_elements) IS 'Time of perihelion passage (Julian date)';
+
+CREATE FUNCTION oe_h_mag(orbital_elements) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION oe_h_mag(orbital_elements) IS 'Absolute magnitude H (NaN if unknown)';
+
+CREATE FUNCTION oe_g_slope(orbital_elements) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION oe_g_slope(orbital_elements) IS 'Slope parameter G (NaN if unknown)';
+
+CREATE FUNCTION oe_semi_major_axis(orbital_elements) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION oe_semi_major_axis(orbital_elements) IS 'Semi-major axis a = q/(1-e) in AU. NULL for parabolic/hyperbolic orbits (e >= 1).';
+
+CREATE FUNCTION oe_period_years(orbital_elements) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION oe_period_years(orbital_elements) IS 'Orbital period in years = a^1.5 (Kepler third law). NULL for parabolic/hyperbolic orbits (e >= 1).';
+
+
+-- ============================================================
+-- MPC MPCORB.DAT parser
+-- ============================================================
+
+CREATE FUNCTION oe_from_mpc(text) RETURNS orbital_elements
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION oe_from_mpc(text) IS
+ 'Parse one MPCORB.DAT fixed-width line into orbital_elements. Converts MPC packed epoch, computes perihelion distance and tp from (a, e, M).';
+
+
+-- ============================================================
+-- Observation functions
+-- ============================================================
+
+CREATE FUNCTION small_body_heliocentric(orbital_elements, timestamptz) RETURNS heliocentric
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION small_body_heliocentric(orbital_elements, timestamptz) IS
+ 'Heliocentric ecliptic J2000 position of a comet/asteroid from its orbital elements at a given time.';
+
+CREATE FUNCTION small_body_observe(orbital_elements, observer, timestamptz) RETURNS topocentric
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION small_body_observe(orbital_elements, observer, timestamptz) IS
+ 'Observe a comet/asteroid from orbital elements. Auto-fetches Earth via VSOP87. Returns topocentric az/el with geocentric range in km.';
+-- pg_orrery 0.8.0 -> 0.9.0 migration
+--
+-- Adds equatorial type (apparent RA/Dec of date), atmospheric refraction,
+-- stellar proper motion, and light-time corrected _apparent() functions.
+
+-- ============================================================
+-- equatorial type — apparent RA/Dec of date
+-- ============================================================
+
+CREATE TYPE equatorial;
+
+CREATE FUNCTION equatorial_in(cstring) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION equatorial_out(equatorial) RETURNS cstring
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION equatorial_recv(internal) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION equatorial_send(equatorial) RETURNS bytea
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE TYPE equatorial (
+ INPUT = equatorial_in,
+ OUTPUT = equatorial_out,
+ RECEIVE = equatorial_recv,
+ SEND = equatorial_send,
+ INTERNALLENGTH = 24,
+ ALIGNMENT = double,
+ STORAGE = plain
+);
+
+COMMENT ON TYPE equatorial IS
+ 'Apparent equatorial coordinates of date: RA (hours), Dec (degrees), distance (km). Solar system: J2000 precessed via IAU 1976. Satellites: TEME frame (~of-date to ~arcsecond). 24 bytes, fixed-size.';
+
+
+-- ============================================================
+-- Equatorial accessor functions
+-- ============================================================
+
+CREATE FUNCTION eq_ra(equatorial) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION eq_ra(equatorial) IS 'Right ascension in hours [0, 24)';
+
+CREATE FUNCTION eq_dec(equatorial) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION eq_dec(equatorial) IS 'Declination in degrees [-90, 90]';
+
+CREATE FUNCTION eq_distance(equatorial) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION eq_distance(equatorial) IS 'Distance in km (0 for stars without parallax)';
+
+
+-- ============================================================
+-- Satellite RA/Dec functions
+-- ============================================================
+
+CREATE FUNCTION eci_to_equatorial(eci_position, observer, timestamptz) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION eci_to_equatorial(eci_position, observer, timestamptz) IS
+ 'Topocentric apparent RA/Dec from ECI position. Observer parallax-corrected — LEO parallax is ~1 degree.';
+
+CREATE FUNCTION eci_to_equatorial_geo(eci_position, timestamptz) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION eci_to_equatorial_geo(eci_position, timestamptz) IS
+ 'Geocentric apparent RA/Dec from ECI position. Observer-independent — the direction of the TEME position vector.';
+
+
+-- ============================================================
+-- Solar system equatorial functions (VSOP87)
+-- ============================================================
+
+CREATE FUNCTION planet_equatorial(int4, timestamptz) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION planet_equatorial(int4, timestamptz) IS
+ 'Geocentric apparent RA/Dec of a planet via VSOP87. Body IDs: 1=Mercury through 8=Neptune.';
+
+CREATE FUNCTION sun_equatorial(timestamptz) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_equatorial(timestamptz) IS
+ 'Geocentric apparent RA/Dec of the Sun via VSOP87.';
+
+CREATE FUNCTION moon_equatorial(timestamptz) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION moon_equatorial(timestamptz) IS
+ 'Geocentric apparent RA/Dec of the Moon via ELP2000-82B.';
+
+CREATE FUNCTION small_body_equatorial(orbital_elements, timestamptz) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION small_body_equatorial(orbital_elements, timestamptz) IS
+ 'Geocentric apparent RA/Dec of a comet/asteroid from orbital elements. Earth via VSOP87.';
+
+CREATE FUNCTION star_equatorial(float8, float8, timestamptz) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION star_equatorial(float8, float8, timestamptz) IS
+ 'Apparent RA/Dec of a star at a given time. Precesses J2000 catalog coordinates (RA hours, Dec degrees) to date via IAU 1976.';
+
+
+-- ============================================================
+-- Atmospheric refraction (Bennett 1982)
+-- ============================================================
+
+CREATE FUNCTION atmospheric_refraction(float8) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION atmospheric_refraction(float8) IS
+ 'Atmospheric refraction correction in degrees for a given geometric elevation (degrees). Standard atmosphere: P=1010 mbar, T=10C. Bennett (1982) formula with domain guard at -1 degree.';
+
+CREATE FUNCTION atmospheric_refraction_ext(float8, float8, float8) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION atmospheric_refraction_ext(float8, float8, float8) IS
+ 'Atmospheric refraction with pressure/temperature correction. Args: elevation_deg, pressure_mbar, temperature_celsius. Meeus P/T factor applied to Bennett formula.';
+
+CREATE FUNCTION topo_elevation_apparent(topocentric) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION topo_elevation_apparent(topocentric) IS
+ 'Apparent elevation in degrees — geometric elevation plus atmospheric refraction correction.';
+
+
+-- ============================================================
+-- Refracted pass prediction
+-- ============================================================
+
+CREATE FUNCTION predict_passes_refracted(
+ tle, observer, timestamptz, timestamptz, float8 DEFAULT 0.0
+) RETURNS SETOF pass_event
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE
+ ROWS 20;
+COMMENT ON FUNCTION predict_passes_refracted(tle, observer, timestamptz, timestamptz, float8) IS
+ 'Predict satellite passes using a refracted horizon threshold (-0.569 deg geometric). Atmospheric refraction makes satellites visible ~35 seconds earlier at AOS and later at LOS.';
+
+
+-- ============================================================
+-- Stellar proper motion
+-- ============================================================
+
+CREATE FUNCTION star_observe_pm(
+ float8, float8, float8, float8, float8, float8, observer, timestamptz
+) RETURNS topocentric
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION star_observe_pm(float8, float8, float8, float8, float8, float8, observer, timestamptz) IS
+ 'Observe a star with proper motion. Args: ra_hours, dec_deg, pm_ra_masyr (mu_alpha*cos(delta)), pm_dec_masyr, parallax_mas, rv_kms, observer, time. Hipparcos/Gaia convention for pm_ra.';
+
+CREATE FUNCTION star_equatorial_pm(
+ float8, float8, float8, float8, float8, float8, timestamptz
+) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION star_equatorial_pm(float8, float8, float8, float8, float8, float8, timestamptz) IS
+ 'Apparent RA/Dec of a star with proper motion. Args: ra_hours, dec_deg, pm_ra_masyr, pm_dec_masyr, parallax_mas, rv_kms, time. Distance from parallax if > 0.';
+
+
+-- ============================================================
+-- Light-time corrected observation functions
+-- ============================================================
+
+CREATE FUNCTION planet_observe_apparent(int4, observer, timestamptz) RETURNS topocentric
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION planet_observe_apparent(int4, observer, timestamptz) IS
+ 'Observe a planet with single-iteration light-time correction. Body at retarded time, Earth at observation time. VSOP87.';
+
+CREATE FUNCTION sun_observe_apparent(observer, timestamptz) RETURNS topocentric
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_observe_apparent(observer, timestamptz) IS
+ 'Observe the Sun with light-time correction (~8.3 min). VSOP87.';
+
+CREATE FUNCTION small_body_observe_apparent(orbital_elements, observer, timestamptz) RETURNS topocentric
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION small_body_observe_apparent(orbital_elements, observer, timestamptz) IS
+ 'Observe a comet/asteroid with single-iteration light-time correction. Kepler propagation at retarded time, Earth via VSOP87 at observation time.';
+
+
+-- ============================================================
+-- Light-time corrected equatorial functions
+-- ============================================================
+
+CREATE FUNCTION planet_equatorial_apparent(int4, timestamptz) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION planet_equatorial_apparent(int4, timestamptz) IS
+ 'Geocentric apparent RA/Dec of a planet with light-time correction. VSOP87.';
+
+CREATE FUNCTION moon_equatorial_apparent(timestamptz) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION moon_equatorial_apparent(timestamptz) IS
+ 'Geocentric apparent RA/Dec of the Moon with light-time correction (~1.3 sec). ELP2000-82B.';
+
+CREATE FUNCTION small_body_equatorial_apparent(orbital_elements, timestamptz) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION small_body_equatorial_apparent(orbital_elements, timestamptz) IS
+ 'Geocentric apparent RA/Dec of a comet/asteroid with light-time correction.';
+
+
+-- ============================================================
+-- DE ephemeris equatorial variants (STABLE)
+-- ============================================================
+
+CREATE FUNCTION planet_equatorial_de(int4, timestamptz) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION planet_equatorial_de(int4, timestamptz) IS
+ 'Geocentric apparent RA/Dec of a planet via JPL DE ephemeris (falls back to VSOP87 + equatorial).';
+
+CREATE FUNCTION moon_equatorial_de(timestamptz) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION moon_equatorial_de(timestamptz) IS
+ 'Geocentric apparent RA/Dec of the Moon via JPL DE ephemeris (falls back to ELP2000-82B + equatorial).';
+-- pg_orrery 0.9.0 -> 0.10.0 migration
+--
+-- Adds annual aberration to existing _apparent() functions,
+-- 6 new _apparent_de() variants, equatorial angular separation
+-- operator and cone predicate, and stellar annual parallax.
+
+-- ============================================================
+-- Equatorial angular distance and cone search
+-- ============================================================
+
+CREATE FUNCTION eq_angular_distance(equatorial, equatorial) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION eq_angular_distance(equatorial, equatorial) IS
+ 'Angular separation in degrees between two equatorial positions. Vincenty formula (stable at 0 and 180 degrees).';
+
+CREATE FUNCTION eq_within_cone(equatorial, equatorial, float8) RETURNS bool
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION eq_within_cone(equatorial, equatorial, float8) IS
+ 'True if first position is within radius_deg of second position. Cosine shortcut for fast rejection.';
+
+CREATE OPERATOR <-> (
+ LEFTARG = equatorial,
+ RIGHTARG = equatorial,
+ FUNCTION = eq_angular_distance,
+ COMMUTATOR = <->
+);
+COMMENT ON OPERATOR <-> (equatorial, equatorial) IS
+ 'Angular separation in degrees between two equatorial positions.';
+
+
+-- ============================================================
+-- DE apparent observation functions (STABLE, light-time + aberration)
+-- ============================================================
+
+CREATE FUNCTION planet_observe_apparent_de(int4, observer, timestamptz) RETURNS topocentric
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION planet_observe_apparent_de(int4, observer, timestamptz) IS
+ 'Observe a planet with light-time correction and annual aberration via JPL DE (falls back to VSOP87).';
+
+CREATE FUNCTION sun_observe_apparent_de(observer, timestamptz) RETURNS topocentric
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_observe_apparent_de(observer, timestamptz) IS
+ 'Observe the Sun with aberration via JPL DE (falls back to VSOP87).';
+
+CREATE FUNCTION moon_observe_apparent_de(observer, timestamptz) RETURNS topocentric
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION moon_observe_apparent_de(observer, timestamptz) IS
+ 'Observe the Moon with light-time correction and annual aberration via JPL DE (falls back to ELP2000-82B).';
+
+CREATE FUNCTION planet_equatorial_apparent_de(int4, timestamptz) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION planet_equatorial_apparent_de(int4, timestamptz) IS
+ 'Geocentric apparent RA/Dec of a planet with light-time correction and annual aberration via JPL DE (falls back to VSOP87).';
+
+CREATE FUNCTION moon_equatorial_apparent_de(timestamptz) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION moon_equatorial_apparent_de(timestamptz) IS
+ 'Geocentric apparent RA/Dec of the Moon with light-time correction and annual aberration via JPL DE (falls back to ELP2000-82B).';
+
+CREATE FUNCTION small_body_observe_apparent_de(orbital_elements, observer, timestamptz) RETURNS topocentric
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION small_body_observe_apparent_de(orbital_elements, observer, timestamptz) IS
+ 'Observe a comet/asteroid with light-time correction and annual aberration. Earth position via JPL DE (falls back to VSOP87).';
+-- pg_orrery 0.10.0 -> 0.11.0 migration
+--
+-- Adds make_orbital_elements() constructors and
+-- geocentric equatorial functions for planetary moons.
+
+-- ============================================================
+-- orbital_elements constructors
+-- ============================================================
+
+CREATE FUNCTION make_orbital_elements(
+ epoch_jd float8, q_au float8, e float8,
+ inc_rad float8, omega_rad float8, node_rad float8,
+ tp_jd float8, h_mag float8, g_slope float8
+) RETURNS orbital_elements
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION make_orbital_elements(float8,float8,float8,float8,float8,float8,float8,float8,float8) IS
+ 'Construct orbital_elements from 9 floats (angular elements in radians).';
+
+CREATE FUNCTION make_orbital_elements_deg(
+ epoch_jd float8, q_au float8, e float8,
+ inc_deg float8, omega_deg float8, node_deg float8,
+ tp_jd float8, h_mag float8, g_slope float8
+) RETURNS orbital_elements
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION make_orbital_elements_deg(float8,float8,float8,float8,float8,float8,float8,float8,float8) IS
+ 'Construct orbital_elements from 9 floats (angular elements in degrees). Matches text I/O and most catalog column layouts.';
+
+
+-- ============================================================
+-- Planetary moon equatorial functions
+-- ============================================================
+
+CREATE FUNCTION galilean_equatorial(int4, timestamptz) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION galilean_equatorial(int4, timestamptz) IS
+ 'Geometric geocentric RA/Dec of a Galilean moon (0=Io, 1=Europa, 2=Ganymede, 3=Callisto). L1.2 theory + VSOP87. No light-time or aberration correction.';
+
+CREATE FUNCTION saturn_moon_equatorial(int4, timestamptz) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION saturn_moon_equatorial(int4, timestamptz) IS
+ 'Geometric geocentric RA/Dec of a Saturn moon (0=Mimas..7=Hyperion). TASS17 theory + VSOP87. No light-time or aberration correction.';
+
+CREATE FUNCTION uranus_moon_equatorial(int4, timestamptz) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION uranus_moon_equatorial(int4, timestamptz) IS
+ 'Geometric geocentric RA/Dec of a Uranus moon (0=Miranda..4=Oberon). GUST86 theory + VSOP87. No light-time or aberration correction.';
+
+CREATE FUNCTION mars_moon_equatorial(int4, timestamptz) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION mars_moon_equatorial(int4, timestamptz) IS
+ 'Geometric geocentric RA/Dec of a Mars moon (0=Phobos, 1=Deimos). MarsSat theory + VSOP87. No light-time or aberration correction.';
+-- pg_orrery 0.11.0 -> 0.12.0 migration
+--
+-- Adds equatorial GiST operator class for KNN sky queries
+-- and DE moon equatorial functions for all 4 planetary moon families.
+
+-- ============================================================
+-- GiST support functions for equatorial type
+-- ============================================================
+
+CREATE FUNCTION gist_eq_consistent(internal, equatorial, smallint, oid, internal) RETURNS bool
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION gist_eq_union(internal, internal) RETURNS internal
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION gist_eq_compress(internal) RETURNS internal
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION gist_eq_decompress(internal) RETURNS internal
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION gist_eq_penalty(internal, internal, internal) RETURNS internal
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION gist_eq_picksplit(internal, internal) RETURNS internal
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION gist_eq_same(internal, internal, internal) RETURNS internal
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+CREATE FUNCTION gist_eq_distance(internal, equatorial, smallint, oid, internal) RETURNS float8
+ AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+-- ============================================================
+-- Equatorial GiST operator class (KNN ordering only)
+-- ============================================================
+
+CREATE OPERATOR CLASS eq_gist_ops
+ DEFAULT FOR TYPE equatorial USING gist AS
+ OPERATOR 15 <-> (equatorial, equatorial) FOR ORDER BY pg_catalog.float_ops,
+ FUNCTION 1 gist_eq_consistent(internal, equatorial, smallint, oid, internal),
+ FUNCTION 2 gist_eq_union(internal, internal),
+ FUNCTION 3 gist_eq_compress(internal),
+ FUNCTION 4 gist_eq_decompress(internal),
+ FUNCTION 5 gist_eq_penalty(internal, internal, internal),
+ FUNCTION 6 gist_eq_picksplit(internal, internal),
+ FUNCTION 7 gist_eq_same(internal, internal, internal),
+ FUNCTION 8 gist_eq_distance(internal, equatorial, smallint, oid, internal);
+
+-- ============================================================
+-- DE moon equatorial functions (STABLE, fall back to VSOP87)
+-- ============================================================
+
+CREATE FUNCTION galilean_equatorial_de(int4, timestamptz) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION galilean_equatorial_de(int4, timestamptz) IS
+ 'Geocentric RA/Dec of a Galilean moon via DE parent position (falls back to VSOP87). 0=Io, 1=Europa, 2=Ganymede, 3=Callisto.';
+
+CREATE FUNCTION saturn_moon_equatorial_de(int4, timestamptz) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION saturn_moon_equatorial_de(int4, timestamptz) IS
+ 'Geocentric RA/Dec of a Saturn moon via DE parent position (falls back to VSOP87). 0=Mimas..7=Hyperion.';
+
+CREATE FUNCTION uranus_moon_equatorial_de(int4, timestamptz) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION uranus_moon_equatorial_de(int4, timestamptz) IS
+ 'Geocentric RA/Dec of a Uranus moon via DE parent position (falls back to VSOP87). 0=Miranda..4=Oberon.';
+
+CREATE FUNCTION mars_moon_equatorial_de(int4, timestamptz) RETURNS equatorial
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION mars_moon_equatorial_de(int4, timestamptz) IS
+ 'Geocentric RA/Dec of a Mars moon via DE parent position (falls back to VSOP87). 0=Phobos, 1=Deimos.';
+
+
+-- ============================================================
+-- v0.13.0: make_equatorial() constructor
+-- ============================================================
+
+CREATE FUNCTION make_equatorial(ra_hours float8, dec_deg float8, distance_km float8)
+ RETURNS equatorial
+ AS 'MODULE_PATHNAME', 'make_equatorial'
+ LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+
+COMMENT ON FUNCTION make_equatorial(float8, float8, float8) IS
+ 'Construct equatorial from RA (hours [0,24)), Dec (degrees [-90,90]), distance (km).';
+
+
+-- ============================================================
+-- v0.13.0: Rise/set prediction functions
+-- ============================================================
+
+CREATE FUNCTION planet_next_rise(body_id int4, obs observer, t timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION planet_next_rise(int4, observer, timestamptz) IS
+ 'Next geometric rise time for a planet. Returns NULL if no rise within 7 days. body_id: 1=Mercury..8=Neptune.';
+
+CREATE FUNCTION planet_next_set(body_id int4, obs observer, t timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION planet_next_set(int4, observer, timestamptz) IS
+ 'Next geometric set time for a planet. Returns NULL if no set within 7 days. body_id: 1=Mercury..8=Neptune.';
+
+CREATE FUNCTION sun_next_rise(obs observer, t timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_next_rise(observer, timestamptz) IS
+ 'Next geometric sunrise. Returns NULL if Sun does not rise within 7 days (polar night).';
+
+CREATE FUNCTION sun_next_set(obs observer, t timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_next_set(observer, timestamptz) IS
+ 'Next geometric sunset. Returns NULL if Sun does not set within 7 days (midnight sun).';
+
+CREATE FUNCTION moon_next_rise(obs observer, t timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION moon_next_rise(observer, timestamptz) IS
+ 'Next geometric moonrise. Returns NULL if Moon does not rise within 7 days.';
+
+CREATE FUNCTION moon_next_set(obs observer, t timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION moon_next_set(observer, timestamptz) IS
+ 'Next geometric moonset. Returns NULL if Moon does not set within 7 days.';
+
+CREATE FUNCTION sun_next_rise_refracted(obs observer, t timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_next_rise_refracted(observer, timestamptz) IS
+ 'Next refracted sunrise (-0.833 deg threshold: refraction + semidiameter). Earlier than geometric by ~4 min.';
+
+CREATE FUNCTION sun_next_set_refracted(obs observer, t timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_next_set_refracted(observer, timestamptz) IS
+ 'Next refracted sunset (-0.833 deg threshold: refraction + semidiameter). Later than geometric by ~4 min.';
+
+
+-- ============================================================
+-- v0.14.0: Refracted planet/moon rise/set
+-- ============================================================
+
+CREATE FUNCTION planet_next_rise_refracted(body_id int4, obs observer, t timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION planet_next_rise_refracted(int4, observer, timestamptz) IS
+ 'Next refracted rise time for a planet (-0.569 deg threshold: atmospheric refraction only). Earlier than geometric.';
+
+CREATE FUNCTION planet_next_set_refracted(body_id int4, obs observer, t timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION planet_next_set_refracted(int4, observer, timestamptz) IS
+ 'Next refracted set time for a planet (-0.569 deg threshold: atmospheric refraction only). Later than geometric.';
+
+CREATE FUNCTION moon_next_rise_refracted(obs observer, t timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION moon_next_rise_refracted(observer, timestamptz) IS
+ 'Next refracted moonrise (-0.833 deg threshold: refraction + semidiameter). Earlier than geometric.';
+
+CREATE FUNCTION moon_next_set_refracted(obs observer, t timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION moon_next_set_refracted(observer, timestamptz) IS
+ 'Next refracted moonset (-0.833 deg threshold: refraction + semidiameter). Later than geometric.';
+
+
+-- ============================================================
+-- v0.14.0: Constellation identification (Roman 1987, CDS VI/42)
+-- ============================================================
+
+CREATE FUNCTION constellation(eq equatorial) RETURNS text
+ AS 'MODULE_PATHNAME', 'constellation_from_equatorial'
+ LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION constellation(equatorial) IS
+ 'IAU constellation abbreviation (3 letters) from equatorial coordinates (Roman 1987).';
+
+CREATE FUNCTION constellation(ra_hours float8, dec_deg float8) RETURNS text
+ AS 'MODULE_PATHNAME', 'constellation_from_radec'
+ LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION constellation(float8, float8) IS
+ 'IAU constellation from J2000 RA (hours [0,24)) and Dec (degrees [-90,90]).';
+-- pg_orrery 0.14.0 -> 0.15.0 migration
+--
+-- Adds: constellation_full_name (1 function),
+-- rise/set status diagnostics (3 functions).
+
+-- ============================================================
+-- Constellation full name lookup
+-- ============================================================
+
+CREATE FUNCTION constellation_full_name(abbr text) RETURNS text
+ AS 'MODULE_PATHNAME', 'constellation_full_name_from_abbr'
+ LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION constellation_full_name(text) IS
+ 'Full IAU constellation name from 3-letter abbreviation. Returns NULL for invalid abbreviation.';
+
+-- ============================================================
+-- Rise/set status diagnostics
+-- ============================================================
+
+CREATE FUNCTION sun_rise_set_status(obs observer, t timestamptz) RETURNS text
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_rise_set_status(observer, timestamptz) IS
+ 'Classify Sun visibility: rises_and_sets, circumpolar, or never_rises.';
+
+CREATE FUNCTION moon_rise_set_status(obs observer, t timestamptz) RETURNS text
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION moon_rise_set_status(observer, timestamptz) IS
+ 'Classify Moon visibility: rises_and_sets, circumpolar, or never_rises.';
+
+CREATE FUNCTION planet_rise_set_status(body_id int4, obs observer, t timestamptz) RETURNS text
+ AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION planet_rise_set_status(int4, observer, timestamptz) IS
+ 'Classify planet visibility: rises_and_sets, circumpolar, or never_rises. Body IDs 1-8 (Mercury-Neptune).';
+-- pg_orrery 0.15.0 -> 0.16.0: twilight, lunar phase, planet magnitude
+
+-- ============================================================
+-- Twilight functions (6)
+-- ============================================================
+
+CREATE FUNCTION sun_civil_dawn(observer, timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME', 'sun_civil_dawn'
+ LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_civil_dawn(observer, timestamptz) IS
+ 'Next civil dawn (Sun crosses -6 deg rising). Outdoor activities without artificial light.';
+
+CREATE FUNCTION sun_civil_dusk(observer, timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME', 'sun_civil_dusk'
+ LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_civil_dusk(observer, timestamptz) IS
+ 'Next civil dusk (Sun crosses -6 deg setting). Artificial light needed.';
+
+CREATE FUNCTION sun_nautical_dawn(observer, timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME', 'sun_nautical_dawn'
+ LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_nautical_dawn(observer, timestamptz) IS
+ 'Next nautical dawn (Sun crosses -12 deg rising). Horizon visible at sea.';
+
+CREATE FUNCTION sun_nautical_dusk(observer, timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME', 'sun_nautical_dusk'
+ LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_nautical_dusk(observer, timestamptz) IS
+ 'Next nautical dusk (Sun crosses -12 deg setting). Horizon no longer visible at sea.';
+
+CREATE FUNCTION sun_astronomical_dawn(observer, timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME', 'sun_astronomical_dawn'
+ LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_astronomical_dawn(observer, timestamptz) IS
+ 'Next astronomical dawn (Sun crosses -18 deg rising). Sky was fully dark.';
+
+CREATE FUNCTION sun_astronomical_dusk(observer, timestamptz) RETURNS timestamptz
+ AS 'MODULE_PATHNAME', 'sun_astronomical_dusk'
+ LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_astronomical_dusk(observer, timestamptz) IS
+ 'Next astronomical dusk (Sun crosses -18 deg setting). Sky becomes fully dark.';
+
+-- ============================================================
+-- Lunar phase functions (4)
+-- ============================================================
+
+CREATE FUNCTION moon_phase_angle(timestamptz) RETURNS float8
+ AS 'MODULE_PATHNAME', 'moon_phase_angle'
+ LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION moon_phase_angle(timestamptz) IS
+ 'Sun-Earth-Moon phase angle in degrees [0,360). 0=new, 90=first quarter, 180=full, 270=last quarter.';
+
+CREATE FUNCTION moon_illumination(timestamptz) RETURNS float8
+ AS 'MODULE_PATHNAME', 'moon_illumination'
+ LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION moon_illumination(timestamptz) IS
+ 'Illuminated fraction of the Moon disk [0.0, 1.0].';
+
+CREATE FUNCTION moon_phase_name(timestamptz) RETURNS text
+ AS 'MODULE_PATHNAME', 'moon_phase_name'
+ LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION moon_phase_name(timestamptz) IS
+ 'Moon phase name: new_moon, waxing_crescent, first_quarter, waxing_gibbous, full_moon, waning_gibbous, last_quarter, waning_crescent.';
+
+CREATE FUNCTION moon_age(timestamptz) RETURNS float8
+ AS 'MODULE_PATHNAME', 'moon_age'
+ LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION moon_age(timestamptz) IS
+ 'Days since last new moon [0, ~29.53), approximated from phase angle.';
+
+-- ============================================================
+-- Planet magnitude (1)
+-- ============================================================
+
+CREATE FUNCTION planet_magnitude(int4, timestamptz) RETURNS float8
+ AS 'MODULE_PATHNAME', 'planet_magnitude'
+ LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION planet_magnitude(int4, timestamptz) IS
+ 'Apparent visual magnitude of a planet (Mallama & Hilton 2018). Body IDs 1-8. Saturn ring tilt not modeled.';
diff --git a/src/lunar_phase_funcs.c b/src/lunar_phase_funcs.c
new file mode 100644
index 0000000..e029432
--- /dev/null
+++ b/src/lunar_phase_funcs.c
@@ -0,0 +1,236 @@
+/*
+ * lunar_phase_funcs.c -- Lunar phase calculations
+ *
+ * Phase angle, illumination, phase name, and lunar age from
+ * VSOP87 Sun + ELP2000-82B Moon positions.
+ *
+ * The Sun-Earth-Moon geometry determines phase:
+ * 1. Earth heliocentric (VSOP87) -> negate -> Sun geocentric
+ * 2. Moon geocentric ecliptic (ELP2000-82B)
+ * 3. Elongation = angle between geocentric Sun and Moon
+ * 4. Cross product z-component distinguishes waxing from waning
+ */
+
+#include "postgres.h"
+#include "fmgr.h"
+#include "utils/timestamp.h"
+#include "utils/builtins.h"
+#include "types.h"
+#include "astro_math.h"
+#include "vsop87.h"
+#include "elp82b.h"
+#include
+
+PG_FUNCTION_INFO_V1(moon_phase_angle);
+PG_FUNCTION_INFO_V1(moon_illumination);
+PG_FUNCTION_INFO_V1(moon_phase_name);
+PG_FUNCTION_INFO_V1(moon_age);
+
+/* Mean synodic month in days (Meeus, Astronomical Algorithms) */
+#define SYNODIC_MONTH 29.530588853
+
+
+/*
+ * compute_phase_angle -- Sun-Earth-Moon phase angle in degrees [0, 360)
+ *
+ * 0 = new moon (Sun and Moon same direction from Earth)
+ * 90 = first quarter (waxing)
+ * 180 = full moon
+ * 270 = last quarter (waning)
+ *
+ * Waxing/waning determined by the z-component of (sun_geo x moon_geo)
+ * in the ecliptic plane: positive = Moon east of Sun = waxing.
+ */
+static double
+compute_phase_angle(double jd)
+{
+ double earth_xyz[6];
+ double moon_ecl[3];
+ double sun_x, sun_y, sun_z;
+ double moon_x, moon_y, moon_z;
+ double dot, cross_z;
+ double sun_r, moon_r;
+ double elongation, angle;
+
+ /* Earth heliocentric -> Sun geocentric (negate) */
+ GetVsop87Coor(jd, 2, earth_xyz); /* VSOP87 body 2 = Earth */
+ sun_x = -earth_xyz[0];
+ sun_y = -earth_xyz[1];
+ sun_z = -earth_xyz[2];
+
+ /* Moon geocentric ecliptic (AU) */
+ GetElp82bCoor(jd, moon_ecl);
+ moon_x = moon_ecl[0];
+ moon_y = moon_ecl[1];
+ moon_z = moon_ecl[2];
+
+ /* Magnitudes for unit-vector dot product */
+ sun_r = sqrt(sun_x * sun_x + sun_y * sun_y + sun_z * sun_z);
+ moon_r = sqrt(moon_x * moon_x + moon_y * moon_y + moon_z * moon_z);
+
+ /* Elongation = angle between geocentric Sun and Moon directions */
+ dot = (sun_x * moon_x + sun_y * moon_y + sun_z * moon_z)
+ / (sun_r * moon_r);
+ if (dot > 1.0) dot = 1.0;
+ if (dot < -1.0) dot = -1.0;
+ elongation = acos(dot);
+
+ /* Cross product z-component: sign determines waxing vs waning */
+ cross_z = sun_x * moon_y - sun_y * moon_x;
+
+ /*
+ * Elongation is already the right quantity:
+ * elongation 0 = same direction = new moon
+ * elongation pi = opposite = full moon
+ *
+ * acos gives [0, 180]. Use cross_z to extend to [0, 360):
+ * cross_z > 0: Moon east of Sun (waxing), stays [0, 180)
+ * cross_z < 0: Moon west of Sun (waning), maps to [180, 360)
+ */
+ angle = elongation * RAD_TO_DEG;
+
+ if (cross_z < 0.0)
+ angle = 360.0 - angle;
+
+ return angle;
+}
+
+
+/* ================================================================
+ * moon_phase_angle(timestamptz) -> float8
+ *
+ * Returns the Sun-Earth-Moon phase angle in degrees [0, 360).
+ * 0 = new moon
+ * 90 = first quarter
+ * 180 = full moon
+ * 270 = last quarter
+ * ================================================================
+ */
+Datum
+moon_phase_angle(PG_FUNCTION_ARGS)
+{
+ int64 ts = PG_GETARG_INT64(0);
+ double jd;
+
+ jd = timestamptz_to_jd(ts);
+
+ PG_RETURN_FLOAT8(compute_phase_angle(jd));
+}
+
+
+/* ================================================================
+ * moon_illumination(timestamptz) -> float8
+ *
+ * Returns the illuminated fraction of the Moon's disk [0.0, 1.0].
+ * Uses the elongation between geocentric Sun and Moon directions:
+ * k = (1 - cos(elongation)) / 2
+ * ================================================================
+ */
+Datum
+moon_illumination(PG_FUNCTION_ARGS)
+{
+ int64 ts = PG_GETARG_INT64(0);
+ double jd;
+ double earth_xyz[6];
+ double moon_ecl[3];
+ double sun_x, sun_y, sun_z;
+ double dot;
+ double sun_r, moon_r;
+ double elongation, illum;
+
+ jd = timestamptz_to_jd(ts);
+
+ /* Earth heliocentric -> Sun geocentric (negate) */
+ GetVsop87Coor(jd, 2, earth_xyz);
+ sun_x = -earth_xyz[0];
+ sun_y = -earth_xyz[1];
+ sun_z = -earth_xyz[2];
+
+ /* Moon geocentric ecliptic */
+ GetElp82bCoor(jd, moon_ecl);
+
+ /* Elongation between geocentric Sun and Moon */
+ sun_r = sqrt(sun_x * sun_x + sun_y * sun_y + sun_z * sun_z);
+ moon_r = sqrt(moon_ecl[0] * moon_ecl[0] +
+ moon_ecl[1] * moon_ecl[1] +
+ moon_ecl[2] * moon_ecl[2]);
+
+ dot = (sun_x * moon_ecl[0] + sun_y * moon_ecl[1] + sun_z * moon_ecl[2])
+ / (sun_r * moon_r);
+ if (dot > 1.0) dot = 1.0;
+ if (dot < -1.0) dot = -1.0;
+ elongation = acos(dot);
+
+ illum = (1.0 - cos(elongation)) / 2.0;
+
+ PG_RETURN_FLOAT8(illum);
+}
+
+
+/* ================================================================
+ * moon_phase_name(timestamptz) -> text
+ *
+ * Returns one of 8 phase names based on phase angle:
+ * [0, 22.5) or [337.5, 360) -> 'new_moon'
+ * [22.5, 67.5) -> 'waxing_crescent'
+ * [67.5, 112.5) -> 'first_quarter'
+ * [112.5, 157.5) -> 'waxing_gibbous'
+ * [157.5, 202.5) -> 'full_moon'
+ * [202.5, 247.5) -> 'waning_gibbous'
+ * [247.5, 292.5) -> 'last_quarter'
+ * [292.5, 337.5) -> 'waning_crescent'
+ * ================================================================
+ */
+Datum
+moon_phase_name(PG_FUNCTION_ARGS)
+{
+ int64 ts = PG_GETARG_INT64(0);
+ double jd;
+ double angle;
+ const char *name;
+
+ jd = timestamptz_to_jd(ts);
+ angle = compute_phase_angle(jd);
+
+ if (angle < 22.5 || angle >= 337.5)
+ name = "new_moon";
+ else if (angle < 67.5)
+ name = "waxing_crescent";
+ else if (angle < 112.5)
+ name = "first_quarter";
+ else if (angle < 157.5)
+ name = "waxing_gibbous";
+ else if (angle < 202.5)
+ name = "full_moon";
+ else if (angle < 247.5)
+ name = "waning_gibbous";
+ else if (angle < 292.5)
+ name = "last_quarter";
+ else
+ name = "waning_crescent";
+
+ PG_RETURN_TEXT_P(cstring_to_text(name));
+}
+
+
+/* ================================================================
+ * moon_age(timestamptz) -> float8
+ *
+ * Returns the Moon's age in days since the last new moon [0, ~29.53).
+ * Approximation from phase angle and mean synodic month:
+ * age = phase_angle / 360 * 29.530588853
+ * ================================================================
+ */
+Datum
+moon_age(PG_FUNCTION_ARGS)
+{
+ int64 ts = PG_GETARG_INT64(0);
+ double jd;
+ double angle, age;
+
+ jd = timestamptz_to_jd(ts);
+ angle = compute_phase_angle(jd);
+ age = angle / 360.0 * SYNODIC_MONTH;
+
+ PG_RETURN_FLOAT8(age);
+}
diff --git a/src/magnitude_funcs.c b/src/magnitude_funcs.c
new file mode 100644
index 0000000..77cbc07
--- /dev/null
+++ b/src/magnitude_funcs.c
@@ -0,0 +1,144 @@
+/*
+ * magnitude_funcs.c -- Planet apparent visual magnitude
+ *
+ * Uses the Mallama & Hilton (2018) magnitude model with
+ * VSOP87 positions for distances and phase angles.
+ *
+ * Reference: Mallama & Hilton, "Computing Apparent Planetary
+ * Magnitudes for The Astronomical Almanac", A&C vol. 25, 2018.
+ */
+
+#include "postgres.h"
+#include "fmgr.h"
+#include "utils/timestamp.h"
+#include "types.h"
+#include "astro_math.h"
+#include "vsop87.h"
+#include
+
+PG_FUNCTION_INFO_V1(planet_magnitude);
+
+
+/*
+ * Planet magnitude parameters -- Mallama & Hilton (2018), simplified.
+ *
+ * V(1,0) = absolute magnitude at r=1 AU, delta=1 AU, i=0 deg
+ * Phase corrections are polynomial fits to i (phase angle in degrees).
+ *
+ * We use the linear+quadratic terms which are sufficient for
+ * phase angles encountered from Earth (typically <180 deg).
+ *
+ * Saturn caveat: visual magnitude depends heavily on ring tilt
+ * (can vary by ~1.5 mag). The simplified model here uses a fixed
+ * V(1,0) without ring correction.
+ */
+typedef struct {
+ double v10; /* V(1,0) */
+ double c1; /* coefficient for i */
+ double c2; /* coefficient for i^2 */
+ double c3; /* coefficient for i^3 (0 if unused) */
+} mag_params;
+
+static const mag_params planet_mag[] = {
+ [0] = { 0, 0, 0, 0 }, /* Sun: unused placeholder */
+ [1] = { -0.613, 6.328e-2, -1.6336e-3, 0 }, /* Mercury */
+ [2] = { -4.384, 1.044e-3, 3.687e-4, 0 }, /* Venus */
+ [3] = { 0, 0, 0, 0 }, /* Earth: unused */
+ [4] = { -1.601, 2.267e-2, -1.302e-4, 0 }, /* Mars */
+ [5] = { -9.395, 3.7e-4, 0, 0 }, /* Jupiter */
+ [6] = { -8.95, 0, 0, 0 }, /* Saturn (ring tilt NOT modeled) */
+ [7] = { -7.110, 0, 0, 0 }, /* Uranus */
+ [8] = { -7.00, 0, 0, 0 }, /* Neptune */
+};
+
+
+/*
+ * Compute apparent visual magnitude of a planet from Earth.
+ *
+ * Phase angle is the Sun-Planet-Earth angle, computed via the law
+ * of cosines from three heliocentric/geocentric distances.
+ */
+static double
+compute_planet_magnitude(int body_id, double jd)
+{
+ double earth_xyz[6], planet_xyz[6];
+ double geo[3];
+ double r, delta, R;
+ double cos_i, i_deg;
+ const mag_params *p;
+ double V;
+ int vsop_body = body_id - 1; /* pg_orrery 1-based -> VSOP87 0-based */
+
+ GetVsop87Coor(jd, 2, earth_xyz); /* Earth (VSOP87 body 2) */
+ GetVsop87Coor(jd, vsop_body, planet_xyz); /* target planet */
+
+ /* Heliocentric distance to planet */
+ r = sqrt(planet_xyz[0] * planet_xyz[0] +
+ planet_xyz[1] * planet_xyz[1] +
+ planet_xyz[2] * planet_xyz[2]);
+
+ /* Geocentric vector and distance */
+ geo[0] = planet_xyz[0] - earth_xyz[0];
+ geo[1] = planet_xyz[1] - earth_xyz[1];
+ geo[2] = planet_xyz[2] - earth_xyz[2];
+ delta = sqrt(geo[0] * geo[0] + geo[1] * geo[1] + geo[2] * geo[2]);
+
+ /* Sun-Earth distance */
+ R = sqrt(earth_xyz[0] * earth_xyz[0] +
+ earth_xyz[1] * earth_xyz[1] +
+ earth_xyz[2] * earth_xyz[2]);
+
+ /* Phase angle via law of cosines: triangle Sun-Planet-Earth */
+ cos_i = (r * r + delta * delta - R * R) / (2.0 * r * delta);
+ if (cos_i > 1.0) cos_i = 1.0;
+ if (cos_i < -1.0) cos_i = -1.0;
+ i_deg = acos(cos_i) * RAD_TO_DEG;
+
+ /* Mallama & Hilton (2018) magnitude formula */
+ p = &planet_mag[body_id];
+ V = p->v10
+ + 5.0 * log10(r * delta)
+ + p->c1 * i_deg
+ + p->c2 * i_deg * i_deg
+ + p->c3 * i_deg * i_deg * i_deg;
+
+ return V;
+}
+
+
+/* ================================================================
+ * planet_magnitude(body_id int4, timestamptz) -> float8
+ *
+ * Returns the apparent visual magnitude of a planet as seen from
+ * Earth. Uses Mallama & Hilton (2018) magnitude model.
+ *
+ * Body IDs: 1=Mercury, ..., 8=Neptune (not Sun 0, Earth 3, or Moon 10)
+ *
+ * NOTE: Saturn magnitude does not account for ring tilt, which
+ * can vary the apparent magnitude by ~1.5 mag. The returned value
+ * is approximate for Saturn.
+ * ================================================================
+ */
+Datum
+planet_magnitude(PG_FUNCTION_ARGS)
+{
+ int32 body_id = PG_GETARG_INT32(0);
+ int64 ts = PG_GETARG_INT64(1);
+ double jd, mag;
+
+ if (body_id < BODY_MERCURY || body_id > BODY_NEPTUNE)
+ ereport(ERROR,
+ (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+ errmsg("planet_magnitude: 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 compute magnitude for Earth from Earth")));
+
+ jd = timestamptz_to_jd(ts);
+ mag = compute_planet_magnitude(body_id, jd);
+
+ PG_RETURN_FLOAT8(mag);
+}
diff --git a/src/rise_set_funcs.c b/src/rise_set_funcs.c
index 3ef1ddf..e41d4bd 100644
--- a/src/rise_set_funcs.c
+++ b/src/rise_set_funcs.c
@@ -38,6 +38,12 @@ PG_FUNCTION_INFO_V1(moon_next_set_refracted);
PG_FUNCTION_INFO_V1(sun_rise_set_status);
PG_FUNCTION_INFO_V1(moon_rise_set_status);
PG_FUNCTION_INFO_V1(planet_rise_set_status);
+PG_FUNCTION_INFO_V1(sun_civil_dawn);
+PG_FUNCTION_INFO_V1(sun_civil_dusk);
+PG_FUNCTION_INFO_V1(sun_nautical_dawn);
+PG_FUNCTION_INFO_V1(sun_nautical_dusk);
+PG_FUNCTION_INFO_V1(sun_astronomical_dawn);
+PG_FUNCTION_INFO_V1(sun_astronomical_dusk);
#define COARSE_STEP_JD (60.0 / 86400.0) /* 60 seconds */
#define BISECT_TOL_JD (0.1 / 86400.0) /* 0.1 second */
@@ -65,6 +71,11 @@ PG_FUNCTION_INFO_V1(planet_rise_set_status);
*/
#define REFRACTION_ONLY_HORIZON_RAD (-0.00993) /* -0.569 deg */
+/* Twilight depression angles (geometric Sun center below horizon) */
+#define CIVIL_TWILIGHT_RAD (-0.10472) /* -6.0 deg */
+#define NAUTICAL_TWILIGHT_RAD (-0.20944) /* -12.0 deg */
+#define ASTRONOMICAL_TWILIGHT_RAD (-0.30416) /* -18.0 deg */
+
/* ----------------------------------------------------------------
* elevation_at_jd_body -- compute topocentric elevation for a body
@@ -684,3 +695,177 @@ planet_rise_set_status(PG_FUNCTION_ARGS)
PG_RETURN_TEXT_P(cstring_to_text(status));
}
+
+
+/* ================================================================
+ * sun_civil_dawn(observer, timestamptz) -> timestamptz
+ *
+ * Returns the next time civil twilight begins (Sun crosses -6 deg
+ * heading upward). Civil twilight = enough light for outdoor
+ * activities without artificial lighting.
+ * ================================================================
+ */
+Datum
+sun_civil_dawn(PG_FUNCTION_ARGS)
+{
+ pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
+ int64 ts = PG_GETARG_INT64(1);
+ double start_jd, stop_jd, result_jd;
+
+ start_jd = timestamptz_to_jd(ts);
+ stop_jd = start_jd + DEFAULT_WINDOW_DAYS;
+
+ result_jd = find_next_crossing(BTYPE_SUN, 0, obs,
+ start_jd, stop_jd,
+ CIVIL_TWILIGHT_RAD, true);
+
+ if (result_jd < 0.0)
+ PG_RETURN_NULL();
+
+ PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
+}
+
+
+/* ================================================================
+ * sun_civil_dusk(observer, timestamptz) -> timestamptz
+ *
+ * Returns the next time civil twilight ends (Sun crosses -6 deg
+ * heading downward). After civil dusk, outdoor activities require
+ * artificial lighting.
+ * ================================================================
+ */
+Datum
+sun_civil_dusk(PG_FUNCTION_ARGS)
+{
+ pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
+ int64 ts = PG_GETARG_INT64(1);
+ double start_jd, stop_jd, result_jd;
+
+ start_jd = timestamptz_to_jd(ts);
+ stop_jd = start_jd + DEFAULT_WINDOW_DAYS;
+
+ result_jd = find_next_crossing(BTYPE_SUN, 0, obs,
+ start_jd, stop_jd,
+ CIVIL_TWILIGHT_RAD, false);
+
+ if (result_jd < 0.0)
+ PG_RETURN_NULL();
+
+ PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
+}
+
+
+/* ================================================================
+ * sun_nautical_dawn(observer, timestamptz) -> timestamptz
+ *
+ * Returns the next time nautical twilight begins (Sun crosses -12 deg
+ * heading upward). At nautical dawn the horizon becomes visible at
+ * sea and bright stars are still visible for navigation.
+ * ================================================================
+ */
+Datum
+sun_nautical_dawn(PG_FUNCTION_ARGS)
+{
+ pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
+ int64 ts = PG_GETARG_INT64(1);
+ double start_jd, stop_jd, result_jd;
+
+ start_jd = timestamptz_to_jd(ts);
+ stop_jd = start_jd + DEFAULT_WINDOW_DAYS;
+
+ result_jd = find_next_crossing(BTYPE_SUN, 0, obs,
+ start_jd, stop_jd,
+ NAUTICAL_TWILIGHT_RAD, true);
+
+ if (result_jd < 0.0)
+ PG_RETURN_NULL();
+
+ PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
+}
+
+
+/* ================================================================
+ * sun_nautical_dusk(observer, timestamptz) -> timestamptz
+ *
+ * Returns the next time nautical twilight ends (Sun crosses -12 deg
+ * heading downward). After nautical dusk the horizon is no longer
+ * visible at sea; bright stars remain visible.
+ * ================================================================
+ */
+Datum
+sun_nautical_dusk(PG_FUNCTION_ARGS)
+{
+ pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
+ int64 ts = PG_GETARG_INT64(1);
+ double start_jd, stop_jd, result_jd;
+
+ start_jd = timestamptz_to_jd(ts);
+ stop_jd = start_jd + DEFAULT_WINDOW_DAYS;
+
+ result_jd = find_next_crossing(BTYPE_SUN, 0, obs,
+ start_jd, stop_jd,
+ NAUTICAL_TWILIGHT_RAD, false);
+
+ if (result_jd < 0.0)
+ PG_RETURN_NULL();
+
+ PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
+}
+
+
+/* ================================================================
+ * sun_astronomical_dawn(observer, timestamptz) -> timestamptz
+ *
+ * Returns the next time astronomical twilight begins (Sun crosses
+ * -18 deg heading upward). Before astronomical dawn the sky is
+ * fully dark — faintest objects are observable.
+ * ================================================================
+ */
+Datum
+sun_astronomical_dawn(PG_FUNCTION_ARGS)
+{
+ pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
+ int64 ts = PG_GETARG_INT64(1);
+ double start_jd, stop_jd, result_jd;
+
+ start_jd = timestamptz_to_jd(ts);
+ stop_jd = start_jd + DEFAULT_WINDOW_DAYS;
+
+ result_jd = find_next_crossing(BTYPE_SUN, 0, obs,
+ start_jd, stop_jd,
+ ASTRONOMICAL_TWILIGHT_RAD, true);
+
+ if (result_jd < 0.0)
+ PG_RETURN_NULL();
+
+ PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
+}
+
+
+/* ================================================================
+ * sun_astronomical_dusk(observer, timestamptz) -> timestamptz
+ *
+ * Returns the next time astronomical twilight ends (Sun crosses
+ * -18 deg heading downward). After astronomical dusk the sky is
+ * fully dark — faintest objects become observable.
+ * ================================================================
+ */
+Datum
+sun_astronomical_dusk(PG_FUNCTION_ARGS)
+{
+ pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
+ int64 ts = PG_GETARG_INT64(1);
+ double start_jd, stop_jd, result_jd;
+
+ start_jd = timestamptz_to_jd(ts);
+ stop_jd = start_jd + DEFAULT_WINDOW_DAYS;
+
+ result_jd = find_next_crossing(BTYPE_SUN, 0, obs,
+ start_jd, stop_jd,
+ ASTRONOMICAL_TWILIGHT_RAD, false);
+
+ if (result_jd < 0.0)
+ PG_RETURN_NULL();
+
+ PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
+}
diff --git a/test/expected/v016_features.out b/test/expected/v016_features.out
new file mode 100644
index 0000000..daa2b9b
--- /dev/null
+++ b/test/expected/v016_features.out
@@ -0,0 +1,243 @@
+-- v016_features.sql -- Tests for v0.16.0: twilight, lunar phase, planet magnitude
+--
+-- Verifies twilight dawn/dusk, lunar phase calculations,
+-- and planet apparent magnitude.
+CREATE EXTENSION IF NOT EXISTS pg_orrery;
+NOTICE: extension "pg_orrery" already exists, skipping
+-- ============================================================
+-- Twilight: ordering (astronomical < nautical < civil < sunrise)
+-- Eagle, Idaho on the 2024 summer solstice
+-- ============================================================
+-- Dawn ordering: astronomical dawn < nautical dawn < civil dawn < sunrise
+-- Use midnight MDT (07:00 UTC) so all "next" events land on the same morning
+SELECT sun_astronomical_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+ < sun_nautical_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+ AS astro_before_nautical;
+ astro_before_nautical
+-----------------------
+ t
+(1 row)
+
+SELECT sun_nautical_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+ < sun_civil_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+ AS nautical_before_civil;
+ nautical_before_civil
+-----------------------
+ t
+(1 row)
+
+SELECT sun_civil_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+ < sun_next_rise('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+ AS civil_before_sunrise;
+ civil_before_sunrise
+----------------------
+ t
+(1 row)
+
+-- Dusk ordering: sunset < civil dusk < nautical dusk < astronomical dusk
+-- Noon MDT (18:00 UTC) ensures all dusk events are still ahead
+SELECT sun_next_set('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+ < sun_civil_dusk('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+ AS sunset_before_civil;
+ sunset_before_civil
+---------------------
+ t
+(1 row)
+
+SELECT sun_civil_dusk('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+ < sun_nautical_dusk('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+ AS civil_before_nautical;
+ civil_before_nautical
+-----------------------
+ t
+(1 row)
+
+SELECT sun_nautical_dusk('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+ < sun_astronomical_dusk('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+ AS nautical_before_astro;
+ nautical_before_astro
+-----------------------
+ t
+(1 row)
+
+-- ============================================================
+-- Twilight: civil dawn ~30 min before sunrise at mid-latitude
+-- ============================================================
+SELECT extract(epoch FROM
+ sun_next_rise('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+ - sun_civil_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+) BETWEEN 1200 AND 3600
+ AS civil_dawn_reasonable_offset;
+ civil_dawn_reasonable_offset
+------------------------------
+ t
+(1 row)
+
+-- ============================================================
+-- Twilight: high latitude summer -- no astronomical darkness
+-- At 60N in June, astronomical dusk should be NULL (never gets dark enough)
+-- ============================================================
+SELECT sun_astronomical_dusk('(60.0,25.0,0)'::observer, '2024-06-21 00:00:00+00'::timestamptz) IS NULL
+ AS no_astro_dark_60n_summer;
+ no_astro_dark_60n_summer
+--------------------------
+ t
+(1 row)
+
+-- ============================================================
+-- Lunar phase: known full moon (2024-01-25 ~17:54 UTC)
+-- Phase angle should be near 180 deg, illumination near 1.0
+-- ============================================================
+SELECT round(moon_phase_angle('2024-01-25 18:00:00+00'::timestamptz)::numeric, 0)
+ BETWEEN 170 AND 190
+ AS full_moon_angle_near_180;
+ full_moon_angle_near_180
+--------------------------
+ t
+(1 row)
+
+SELECT round(moon_illumination('2024-01-25 18:00:00+00'::timestamptz)::numeric, 2)
+ >= 0.95
+ AS full_moon_high_illumination;
+ full_moon_high_illumination
+-----------------------------
+ t
+(1 row)
+
+SELECT moon_phase_name('2024-01-25 18:00:00+00'::timestamptz) = 'full_moon'
+ AS full_moon_named;
+ full_moon_named
+-----------------
+ t
+(1 row)
+
+-- ============================================================
+-- Lunar phase: known new moon (2024-01-11 ~11:57 UTC)
+-- Phase angle should be near 0 or 360, illumination near 0
+-- ============================================================
+SELECT moon_illumination('2024-01-11 12:00:00+00'::timestamptz)
+ < 0.05
+ AS new_moon_low_illumination;
+ new_moon_low_illumination
+---------------------------
+ t
+(1 row)
+
+SELECT moon_phase_name('2024-01-11 12:00:00+00'::timestamptz) = 'new_moon'
+ AS new_moon_named;
+ new_moon_named
+----------------
+ t
+(1 row)
+
+-- ============================================================
+-- Lunar phase: first quarter (2024-01-18 ~03:53 UTC)
+-- Phase angle near 90, illumination near 0.5
+-- ============================================================
+SELECT round(moon_phase_angle('2024-01-18 04:00:00+00'::timestamptz)::numeric, 0)
+ BETWEEN 80 AND 100
+ AS first_quarter_angle_near_90;
+ first_quarter_angle_near_90
+-----------------------------
+ t
+(1 row)
+
+SELECT moon_illumination('2024-01-18 04:00:00+00'::timestamptz)
+ BETWEEN 0.4 AND 0.6
+ AS first_quarter_half_illuminated;
+ first_quarter_half_illuminated
+--------------------------------
+ t
+(1 row)
+
+SELECT moon_phase_name('2024-01-18 04:00:00+00'::timestamptz) = 'first_quarter'
+ AS first_quarter_named;
+ first_quarter_named
+---------------------
+ t
+(1 row)
+
+-- ============================================================
+-- Moon age: new moon has age near 0, full moon near 14.7
+-- ============================================================
+SELECT moon_age('2024-01-11 12:00:00+00'::timestamptz) < 2.0
+ AS new_moon_young;
+ new_moon_young
+----------------
+ t
+(1 row)
+
+SELECT moon_age('2024-01-25 18:00:00+00'::timestamptz)
+ BETWEEN 12.0 AND 17.0
+ AS full_moon_age_midcycle;
+ full_moon_age_midcycle
+------------------------
+ t
+(1 row)
+
+-- ============================================================
+-- Illumination range: always [0, 1]
+-- ============================================================
+SELECT moon_illumination('2024-06-01 00:00:00+00'::timestamptz) BETWEEN 0.0 AND 1.0
+ AND moon_illumination('2024-09-01 00:00:00+00'::timestamptz) BETWEEN 0.0 AND 1.0
+ AND moon_illumination('2024-12-01 00:00:00+00'::timestamptz) BETWEEN 0.0 AND 1.0
+ AS illumination_always_valid;
+ illumination_always_valid
+---------------------------
+ t
+(1 row)
+
+-- ============================================================
+-- Planet magnitude: Jupiter should be bright (negative mag)
+-- ============================================================
+SELECT planet_magnitude(5, '2024-01-15 00:00:00+00'::timestamptz) < 0.0
+ AS jupiter_is_bright;
+ jupiter_is_bright
+-------------------
+ t
+(1 row)
+
+-- ============================================================
+-- Planet magnitude: Venus is the brightest planet
+-- ============================================================
+SELECT planet_magnitude(2, '2024-06-01 12:00:00+00'::timestamptz)
+ < planet_magnitude(4, '2024-06-01 12:00:00+00'::timestamptz)
+ AS venus_brighter_than_mars;
+ venus_brighter_than_mars
+--------------------------
+ t
+(1 row)
+
+-- ============================================================
+-- Planet magnitude: Neptune is faint (~+7-8)
+-- ============================================================
+SELECT planet_magnitude(8, '2024-01-15 00:00:00+00'::timestamptz) > 7.0
+ AS neptune_is_faint;
+ neptune_is_faint
+------------------
+ t
+(1 row)
+
+-- ============================================================
+-- Planet magnitude: all planets return finite values
+-- ============================================================
+SELECT bool_and(
+ planet_magnitude(body_id, '2024-01-15 00:00:00+00'::timestamptz) IS NOT NULL
+ AND planet_magnitude(body_id, '2024-01-15 00:00:00+00'::timestamptz) > -30
+ AND planet_magnitude(body_id, '2024-01-15 00:00:00+00'::timestamptz) < 30
+) AS all_magnitudes_finite
+FROM (VALUES (1),(2),(4),(5),(6),(7),(8)) AS t(body_id);
+ all_magnitudes_finite
+-----------------------
+ t
+(1 row)
+
+-- ============================================================
+-- Planet magnitude: error cases
+-- ============================================================
+DO $$ BEGIN PERFORM planet_magnitude(0, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'body_id=0(Sun): %', SQLERRM; END $$;
+NOTICE: body_id=0(Sun): planet_magnitude: body_id 0 must be 1-8 (Mercury-Neptune)
+DO $$ BEGIN PERFORM planet_magnitude(3, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'body_id=3(Earth): %', SQLERRM; END $$;
+NOTICE: body_id=3(Earth): cannot compute magnitude for Earth from Earth
+DO $$ BEGIN PERFORM planet_magnitude(9, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'body_id=9: %', SQLERRM; END $$;
+NOTICE: body_id=9: planet_magnitude: body_id 9 must be 1-8 (Mercury-Neptune)
diff --git a/test/sql/v016_features.sql b/test/sql/v016_features.sql
new file mode 100644
index 0000000..4b40e03
--- /dev/null
+++ b/test/sql/v016_features.sql
@@ -0,0 +1,161 @@
+-- v016_features.sql -- Tests for v0.16.0: twilight, lunar phase, planet magnitude
+--
+-- Verifies twilight dawn/dusk, lunar phase calculations,
+-- and planet apparent magnitude.
+CREATE EXTENSION IF NOT EXISTS pg_orrery;
+
+-- ============================================================
+-- Twilight: ordering (astronomical < nautical < civil < sunrise)
+-- Eagle, Idaho on the 2024 summer solstice
+-- ============================================================
+
+-- Dawn ordering: astronomical dawn < nautical dawn < civil dawn < sunrise
+-- Use midnight MDT (07:00 UTC) so all "next" events land on the same morning
+SELECT sun_astronomical_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+ < sun_nautical_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+ AS astro_before_nautical;
+
+SELECT sun_nautical_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+ < sun_civil_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+ AS nautical_before_civil;
+
+SELECT sun_civil_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+ < sun_next_rise('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+ AS civil_before_sunrise;
+
+-- Dusk ordering: sunset < civil dusk < nautical dusk < astronomical dusk
+-- Noon MDT (18:00 UTC) ensures all dusk events are still ahead
+SELECT sun_next_set('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+ < sun_civil_dusk('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+ AS sunset_before_civil;
+
+SELECT sun_civil_dusk('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+ < sun_nautical_dusk('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+ AS civil_before_nautical;
+
+SELECT sun_nautical_dusk('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+ < sun_astronomical_dusk('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+ AS nautical_before_astro;
+
+-- ============================================================
+-- Twilight: civil dawn ~30 min before sunrise at mid-latitude
+-- ============================================================
+
+SELECT extract(epoch FROM
+ sun_next_rise('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+ - sun_civil_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+) BETWEEN 1200 AND 3600
+ AS civil_dawn_reasonable_offset;
+
+-- ============================================================
+-- Twilight: high latitude summer -- no astronomical darkness
+-- At 60N in June, astronomical dusk should be NULL (never gets dark enough)
+-- ============================================================
+
+SELECT sun_astronomical_dusk('(60.0,25.0,0)'::observer, '2024-06-21 00:00:00+00'::timestamptz) IS NULL
+ AS no_astro_dark_60n_summer;
+
+-- ============================================================
+-- Lunar phase: known full moon (2024-01-25 ~17:54 UTC)
+-- Phase angle should be near 180 deg, illumination near 1.0
+-- ============================================================
+
+SELECT round(moon_phase_angle('2024-01-25 18:00:00+00'::timestamptz)::numeric, 0)
+ BETWEEN 170 AND 190
+ AS full_moon_angle_near_180;
+
+SELECT round(moon_illumination('2024-01-25 18:00:00+00'::timestamptz)::numeric, 2)
+ >= 0.95
+ AS full_moon_high_illumination;
+
+SELECT moon_phase_name('2024-01-25 18:00:00+00'::timestamptz) = 'full_moon'
+ AS full_moon_named;
+
+-- ============================================================
+-- Lunar phase: known new moon (2024-01-11 ~11:57 UTC)
+-- Phase angle should be near 0 or 360, illumination near 0
+-- ============================================================
+
+SELECT moon_illumination('2024-01-11 12:00:00+00'::timestamptz)
+ < 0.05
+ AS new_moon_low_illumination;
+
+SELECT moon_phase_name('2024-01-11 12:00:00+00'::timestamptz) = 'new_moon'
+ AS new_moon_named;
+
+-- ============================================================
+-- Lunar phase: first quarter (2024-01-18 ~03:53 UTC)
+-- Phase angle near 90, illumination near 0.5
+-- ============================================================
+
+SELECT round(moon_phase_angle('2024-01-18 04:00:00+00'::timestamptz)::numeric, 0)
+ BETWEEN 80 AND 100
+ AS first_quarter_angle_near_90;
+
+SELECT moon_illumination('2024-01-18 04:00:00+00'::timestamptz)
+ BETWEEN 0.4 AND 0.6
+ AS first_quarter_half_illuminated;
+
+SELECT moon_phase_name('2024-01-18 04:00:00+00'::timestamptz) = 'first_quarter'
+ AS first_quarter_named;
+
+-- ============================================================
+-- Moon age: new moon has age near 0, full moon near 14.7
+-- ============================================================
+
+SELECT moon_age('2024-01-11 12:00:00+00'::timestamptz) < 2.0
+ AS new_moon_young;
+
+SELECT moon_age('2024-01-25 18:00:00+00'::timestamptz)
+ BETWEEN 12.0 AND 17.0
+ AS full_moon_age_midcycle;
+
+-- ============================================================
+-- Illumination range: always [0, 1]
+-- ============================================================
+
+SELECT moon_illumination('2024-06-01 00:00:00+00'::timestamptz) BETWEEN 0.0 AND 1.0
+ AND moon_illumination('2024-09-01 00:00:00+00'::timestamptz) BETWEEN 0.0 AND 1.0
+ AND moon_illumination('2024-12-01 00:00:00+00'::timestamptz) BETWEEN 0.0 AND 1.0
+ AS illumination_always_valid;
+
+-- ============================================================
+-- Planet magnitude: Jupiter should be bright (negative mag)
+-- ============================================================
+
+SELECT planet_magnitude(5, '2024-01-15 00:00:00+00'::timestamptz) < 0.0
+ AS jupiter_is_bright;
+
+-- ============================================================
+-- Planet magnitude: Venus is the brightest planet
+-- ============================================================
+
+SELECT planet_magnitude(2, '2024-06-01 12:00:00+00'::timestamptz)
+ < planet_magnitude(4, '2024-06-01 12:00:00+00'::timestamptz)
+ AS venus_brighter_than_mars;
+
+-- ============================================================
+-- Planet magnitude: Neptune is faint (~+7-8)
+-- ============================================================
+
+SELECT planet_magnitude(8, '2024-01-15 00:00:00+00'::timestamptz) > 7.0
+ AS neptune_is_faint;
+
+-- ============================================================
+-- Planet magnitude: all planets return finite values
+-- ============================================================
+
+SELECT bool_and(
+ planet_magnitude(body_id, '2024-01-15 00:00:00+00'::timestamptz) IS NOT NULL
+ AND planet_magnitude(body_id, '2024-01-15 00:00:00+00'::timestamptz) > -30
+ AND planet_magnitude(body_id, '2024-01-15 00:00:00+00'::timestamptz) < 30
+) AS all_magnitudes_finite
+FROM (VALUES (1),(2),(4),(5),(6),(7),(8)) AS t(body_id);
+
+-- ============================================================
+-- Planet magnitude: error cases
+-- ============================================================
+
+DO $$ BEGIN PERFORM planet_magnitude(0, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'body_id=0(Sun): %', SQLERRM; END $$;
+DO $$ BEGIN PERFORM planet_magnitude(3, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'body_id=3(Earth): %', SQLERRM; END $$;
+DO $$ BEGIN PERFORM planet_magnitude(9, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'body_id=9: %', SQLERRM; END $$;