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

Merge v0.16.0: twilight, lunar phase, planet magnitude

Browse Source
This commit is contained in:
Ryan Malloy 2026-02-26 12:42:13 -07:00
commit 84d13bd7d6
14 changed files with 3264 additions and 152 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
CLAUDE.md
View File

@ -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

9
Makefile
View File

@ -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_).

138
docs/TODO-v0.10.0.md
View File

@ -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

2
docs/astro.config.mjs
View File

@ -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" },
],
},
{

206
docs/src/content/docs/guides/constellation-identification.mdx Normal file
View File

@ -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
<Aside type="caution" title="Boundary precision">
The Roman boundary table defines constellation boundaries as line segments in RA and Dec at epoch B1875.0. For positions very close to a boundary (within a few arcseconds), the linear scan may disagree with higher-precision boundary implementations that account for the curved nature of precession over short segments.
</Aside>
- **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.

309
docs/src/content/docs/guides/rise-set-prediction.mdx Normal file
View File

@ -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
<Aside type="caution" title="Analytic ephemeris, not almanac precision">
Rise/set times are computed from the VSOP87 and ELP2000-82B analytic theories (~1 arcsecond for planets, ~10 arcseconds for the Moon). The bisection precision is 0.1 seconds, but the underlying positional accuracy limits the practical accuracy to roughly 1-2 seconds for planets and up to 30 seconds for the Moon near extreme latitudes.
</Aside>
- **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:
<Tabs>
<TabItem label="Sun">
```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.
</TabItem>
<TabItem label="Moon">
```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.
</TabItem>
<TabItem label="Jupiter">
```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.
</TabItem>
</Tabs>
### 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.
<Aside type="note" title="The NULL contract">
All rise/set functions return NULL when the body does not cross the horizon within the 7-day search window. Three scenarios produce NULL:
- **Circumpolar** -- the body is always above the horizon (midnight sun, circumpolar stars).
- **Never rises** -- the body is always below the horizon (polar night at high latitudes in winter).
- **Slow-moving body near threshold** -- rare, but a planet near the celestial pole can hover within a degree of the horizon for days without cleanly crossing.
Use `sun_rise_set_status()`, `moon_rise_set_status()`, or `planet_rise_set_status()` to distinguish these cases. Each returns one of three strings: `'rises_and_sets'`, `'circumpolar'`, or `'never_rises'`.
</Aside>
### 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.

2
pg_orrery.control
View File

@ -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

79
sql/pg_orrery--0.15.0--0.16.0.sql Normal file
View File

@ -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.';

1674
sql/pg_orrery--0.16.0.sql Normal file
View File

File diff suppressed because it is too large Load Diff

236
src/lunar_phase_funcs.c Normal file
View File

@ -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 <math.h>
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);
}

144
src/magnitude_funcs.c Normal file
View File

@ -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 <math.h>
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);
}

185
src/rise_set_funcs.c
View File

@ -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));
}

243
test/expected/v016_features.out Normal file
View File

@ -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)

161
test/sql/v016_features.sql Normal file
View File

@ -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 $$;