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

Add v0.19.0: sun almanac, conjunction detection, penumbral fraction, physical libration

Browse Source 
Four new functions (184 → 188 SQL objects):
- sun_almanac_events(): merged rise/set + twilight SRF (4 threshold scans)
- planet_conjunctions(): angular separation minima via daily scan + ternary search
- satellite_penumbral_fraction(): continuous 0.0-1.0 shadow depth
- moon_physical_libration(): Meeus p. 373 Fourier corrections (tau, rho)

30 regression test suites, all passing.
This commit is contained in:
Ryan Malloy 2026-02-28 13:51:35 -07:00
parent 16909c3007
commit 4d64b78fb8
12 changed files with 3278 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

38
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, 184 SQL objects (168 user-visible functions + 16 GiST support), 9 custom types, covering satellites (SGP4/SDP4), planets (VSOP87 + optional JPL DE441), Moon (ELP2000-82B), 19 planetary moons (L1.2/TASS17/GUST86/MarsSat), stars (with proper motion and annual parallax), comets, asteroids (MPC catalog), Jupiter radio bursts, interplanetary Lambert transfers, equatorial RA/Dec coordinates with GiST-indexed angular separation, atmospheric refraction, annual stellar aberration, light-time correction, rise/set prediction (geometric + refracted + event windows) with status diagnostics, IAU constellation identification with full name lookup (Roman 1987), twilight dawn/dusk (civil/nautical/astronomical), lunar phase (angle, illumination, name, age), planet apparent magnitude with Saturn ring correction (Mallama & Hilton 2018), solar elongation, planet phase fraction, satellite eclipse prediction (conical shadow with penumbra), observing night quality assessment, lunar optical libration (Meeus Ch. 53), and angular separation rate.
A database orrery — celestial mechanics types and functions for PostgreSQL. Native C extension using PGXS, 188 SQL objects (172 user-visible functions + 16 GiST support), 9 custom types, covering satellites (SGP4/SDP4), planets (VSOP87 + optional JPL DE441), Moon (ELP2000-82B), 19 planetary moons (L1.2/TASS17/GUST86/MarsSat), stars (with proper motion and annual parallax), comets, asteroids (MPC catalog), Jupiter radio bursts, interplanetary Lambert transfers, equatorial RA/Dec coordinates with GiST-indexed angular separation, atmospheric refraction, annual stellar aberration, light-time correction, rise/set prediction (geometric + refracted + event windows + sun almanac) with status diagnostics, IAU constellation identification with full name lookup (Roman 1987), twilight dawn/dusk (civil/nautical/astronomical), lunar phase (angle, illumination, name, age), planet apparent magnitude with Saturn ring correction (Mallama & Hilton 2018), solar elongation, planet phase fraction, conjunction detection, satellite eclipse prediction (conical shadow with penumbral fraction), observing night quality assessment, lunar libration (optical + physical, Meeus Ch. 53 + p. 373), and angular separation rate.
**Current version:** 0.18.0
**Current version:** 0.19.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 29 regression test suites
make installcheck PG_CONFIG=/usr/bin/pg_config # Run 30 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.18.0)
pg_orrery.control # Extension metadata (version 0.19.0)
Makefile # PGXS build + Docker targets
sql/
pg_orrery--0.1.0.sql # v0.1.0: satellite types/functions/operators
@ -48,6 +48,7 @@ sql/
pg_orrery--0.16.0.sql # v0.16.0: twilight, lunar phase, planet magnitude (162 objects)
pg_orrery--0.17.0.sql # v0.17.0: elongation, phase, eclipse, night quality, libration (174 objects)
pg_orrery--0.18.0.sql # v0.18.0: ring tilt, penumbral eclipse, rise/set windows, angular rate (184 objects)
pg_orrery--0.19.0.sql # v0.19.0: sun almanac, conjunctions, penumbral fraction, physical libration (188 objects)
pg_orrery--0.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
@ -65,6 +66,7 @@ sql/
pg_orrery--0.15.0--0.16.0.sql # Migration: v0.15.0 → v0.16.0 (twilight, lunar phase, planet magnitude)
pg_orrery--0.16.0--0.17.0.sql # Migration: v0.16.0 → v0.17.0 (elongation, phase, eclipse, night quality, libration)
pg_orrery--0.17.0--0.18.0.sql # Migration: v0.17.0 → v0.18.0 (ring tilt, penumbral eclipse, rise/set windows, angular rate)
pg_orrery--0.18.0--0.19.0.sql # Migration: v0.18.0 → v0.19.0 (sun almanac, conjunctions, penumbral fraction, physical libration)
src/
pg_orrery.c # PG_MODULE_MAGIC + _PG_init() (GUC registration)
types.h # All struct definitions + constants + DE body ID mapping
@ -89,15 +91,15 @@ src/
kepler_funcs.c # kepler_propagate(), comet_observe()
kepler.h # Shared Kepler solver interface (kepler_position())
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, angular rate
equatorial_funcs.c # equatorial type I/O, accessors, satellite/planet/sun/moon RA/Dec, angular rate, conjunction detection
refraction_funcs.c # atmospheric_refraction(), _ext(), topo_elevation_apparent()
rise_set_funcs.c # planet/sun/moon rise/set (geometric + refracted) + twilight dawn/dusk + event window SRFs
rise_set_funcs.c # planet/sun/moon rise/set (geometric + refracted) + twilight dawn/dusk + event window SRFs + sun almanac
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() (with Saturn ring correction), solar_elongation(), planet_phase(), saturn_ring_tilt()
eclipse_funcs.c # satellite eclipse prediction (conical shadow with penumbra, Vallado §5.3)
libration.h / libration_funcs.c # lunar optical libration (Meeus Ch. 53)
eclipse_funcs.c # satellite eclipse prediction (conical shadow with penumbral fraction, Vallado §5.3)
libration.h / libration_funcs.c # lunar libration (optical Meeus Ch. 53 + physical p. 373)
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)
@ -122,7 +124,7 @@ src/
PROVENANCE.md # Vendoring decision, modifications, verification
LICENSE # MIT license (Bill Gray / Project Pluto)
test/
sql/ # 29 regression test suites
sql/ # 30 regression test suites
expected/ # Expected output
data/vallado_518.json # 518 Vallado test vectors (AIAA 2006-6753-Rev1)
docs/
@ -149,7 +151,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 (184 SQL objects)
## Function Domains (188 SQL objects)
| Domain | Theory | Key Functions | Count |
|--------|--------|---------------|-------|
@ -161,20 +163,21 @@ All types are fixed-size, `STORAGE = plain`, `ALIGNMENT = double`. No TOAST over
| Comets/asteroids | Two-body Keplerian + MPC | `small_body_observe()`, `small_body_equatorial()`, `oe_from_mpc()` | 19 |
| Refraction | Bennett (1982) | `atmospheric_refraction()`, `predict_passes_refracted()` | 4 |
| Equatorial spatial | Vincenty formula | `eq_angular_distance()`, `eq_within_cone()`, `eq_angular_rate()`, `<->` | 4 |
| Conjunction detection | VSOP87/ELP2000-82B + ternary search | `planet_conjunctions()` | 1 |
| Jupiter radio | Carr et al. (1983) | `jupiter_burst_probability()` | 3 |
| Transfers | Lambert (Izzo 2015) | `lambert_transfer()`, `lambert_c3()` | 2 |
| DE ephemeris | JPL DE440/441 (optional) | `planet_observe_de()`, `*_equatorial_de()`, `*_apparent_de()` | 23 |
| 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()`, `*_rise_set_events()` | 18 |
| Rise/set | Bisection (60s scan) | `planet_next_rise()`, `sun_next_rise_refracted()`, `*_rise_set_status()`, `*_rise_set_events()`, `sun_almanac_events()` | 19 |
| 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()`, `saturn_ring_tilt()` | 2 |
| Solar elongation | VSOP87 geometry | `solar_elongation()` | 1 |
| Planet phase | VSOP87 geometry | `planet_phase()` | 1 |
| Satellite eclipse | Conical shadow (Vallado §5.3) | `satellite_is_eclipsed()`, `satellite_next_eclipse_entry()`, `satellite_shadow_state()`, `satellite_in_penumbra()` | 8 |
| Satellite eclipse | Conical shadow (Vallado §5.3) | `satellite_is_eclipsed()`, `satellite_next_eclipse_entry()`, `satellite_shadow_state()`, `satellite_penumbral_fraction()` | 9 |
| Observing quality | Composite (twilight+Moon) | `observing_night_quality()` | 1 |
| Lunar libration | Meeus (1998) Ch. 53 | `moon_libration_longitude()`, `moon_libration()`, `moon_subsolar_longitude()` | 5 |
| Lunar libration | Meeus (1998) Ch. 53 + p. 373 | `moon_libration_longitude()`, `moon_libration()`, `moon_subsolar_longitude()`, `moon_physical_libration()` | 6 |
| Constellation | Roman (1987) CDS VI/42 | `constellation()`, `constellation_full_name()` | 3 |
| Diagnostics | -- | `pg_orrery_ephemeris_info()` | 1 |
@ -309,7 +312,7 @@ All numerical logic is byte-identical to upstream. Verified against 518 Vallado
## Testing
29 regression test suites via `make installcheck`:
30 regression test suites via `make installcheck`:
| Suite | What it tests |
|-------|--------------|
@ -342,10 +345,11 @@ All numerical logic is byte-identical to upstream. Verified against 518 Vallado
| v016_features | Twilight ordering/offset/polar, lunar phase at known events, planet magnitude ranges/errors |
| v017_features | Solar elongation ranges/errors, planet phase ranges, satellite eclipse, observing night quality, lunar libration ranges, subsolar longitude |
| v018_features | Saturn ring tilt range/variation, penumbral eclipse (shadow state, penumbra precedes umbra), rise/set event windows (Sun/Moon/planet, refracted vs geometric), angular separation rate (generic + planet convenience) |
| v019_features | Sun almanac events (count/order/types/polar/refraction/window guard), conjunction detection (Jupiter-Saturn 2020, Moon-Venus, same-body error, threshold filter), penumbral fraction (range/bounds/eclipse consistency), physical libration (small corrections, time variation, total libration range) |
### PG Version Matrix
Test all 29 regression suites + DE reader unit test across PostgreSQL 14-18 using Docker:
Test all 30 regression suites + DE reader unit test across PostgreSQL 14-18 using Docker:
```bash
make test-matrix # Full matrix (PG 14-18)
@ -363,7 +367,7 @@ Logs saved to `test/matrix-logs/pg${ver}.log`. The script reuses the Dockerfile
- `_safe()` variants (`sgp4_propagate_safe`, `observe_safe`, `star_observe_safe`) return NULL on error instead of raising exceptions. Use these for batch queries over potentially invalid data.
- SGP4 error codes: -1 (nearly parabolic), -2 (negative semi-major axis/decayed), -3/-4 (orbit within Earth, returns with NOTICE), -5 (negative mean motion), -6 (convergence failure)
- Pass prediction: propagation failures return -pi elevation (below horizon), shedding the failed timestep without aborting the scan.
- Input validation: same-body Lambert check, arrival-before-departure, invalid body_id, RA out of [0,24), negative perihelion distance.
- Input validation: same-body Lambert check, same-body conjunction check, arrival-before-departure, invalid body_id, RA out of [0,24), negative perihelion distance, almanac/conjunction window overflow.
## Documentation Site
@ -371,7 +375,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 184 SQL objects incl. DE variants, equatorial GiST, refraction, rise/set, constellation, twilight, lunar phase, planet magnitude, Saturn ring tilt, solar elongation, planet phase, satellite eclipse with penumbra, observing quality, lunar libration, angular separation rate), Architecture (Hamilton's principles, constant custody, observation pipeline), Performance (benchmarks).
Sections: Getting Started, Guides (9 domain walkthroughs incl. DE ephemeris), Workflow Translation (Skyfield/Horizons/GMAT/Radio Jupiter Pro comparisons), Reference (all 188 SQL objects incl. DE variants, equatorial GiST, refraction, rise/set, sun almanac, conjunction detection, constellation, twilight, lunar phase, planet magnitude, Saturn ring tilt, solar elongation, planet phase, satellite eclipse with penumbral fraction, observing quality, lunar libration with physical corrections, angular separation rate), Architecture (Hamilton's principles, constant custody, observation pipeline), Performance (benchmarks).
### Local Development
```bash

6
Makefile
View File

@ -16,7 +16,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.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 \
sql/pg_orrery--0.17.0.sql sql/pg_orrery--0.16.0--0.17.0.sql \
sql/pg_orrery--0.18.0.sql sql/pg_orrery--0.17.0--0.18.0.sql
sql/pg_orrery--0.18.0.sql sql/pg_orrery--0.17.0--0.18.0.sql \
sql/pg_orrery--0.19.0.sql sql/pg_orrery--0.18.0--0.19.0.sql
# Our extension C sources
OBJS = src/pg_orrery.o src/tle_type.o src/eci_type.o src/observer_type.o \
@ -60,7 +61,8 @@ REGRESS = tle_parse sgp4_propagate coord_transforms pass_prediction gist_index c
v015_features \
v016_features \
v017_features \
v018_features
v018_features \
v019_features
REGRESS_OPTS = --inputdir=test
# Pure C — no C++ runtime needed. LAPACK for OD solver (dgelss_).

2
pg_orrery.control
View File

@ -1,4 +1,4 @@
comment = 'A database orrery — celestial mechanics types and functions for PostgreSQL'
default_version = '0.18.0'
default_version = '0.19.0'
module_pathname = '$libdir/pg_orrery'
relocatable = true

53
sql/pg_orrery--0.18.0--0.19.0.sql Normal file
View File

@ -0,0 +1,53 @@
-- pg_orrery 0.18.0 -> 0.19.0: Sun almanac SRF, conjunction detection,
-- penumbral fraction, physical libration
-- ============================================================
-- Sun almanac events SRF (1)
-- ============================================================
CREATE FUNCTION sun_almanac_events(
observer, start timestamptz, stop timestamptz,
refracted bool DEFAULT false
) RETURNS TABLE(event_time timestamptz, event_type text)
AS 'MODULE_PATHNAME', 'sun_almanac_events'
LANGUAGE C STABLE STRICT PARALLEL SAFE
ROWS 50;
COMMENT ON FUNCTION sun_almanac_events(observer, timestamptz, timestamptz, bool) IS
'All Sun events (rise, set, civil/nautical/astronomical dawn and dusk) within a time window, sorted chronologically. Replaces chained individual twilight queries. Max 366-day window.';
-- ============================================================
-- Conjunction detection SRF (1)
-- ============================================================
CREATE FUNCTION planet_conjunctions(
int4, int4, timestamptz, timestamptz,
max_separation float8 DEFAULT 10.0
) RETURNS TABLE(conjunction_time timestamptz, separation_deg float8)
AS 'MODULE_PATHNAME', 'planet_conjunctions'
LANGUAGE C STABLE STRICT PARALLEL SAFE
ROWS 10;
COMMENT ON FUNCTION planet_conjunctions(int4, int4, timestamptz, timestamptz, float8) IS
'Finds conjunctions (angular separation minima) between two solar system bodies. Body IDs: 0=Sun, 1-8=planets, 10=Moon. max_separation filters results (degrees, default 10). Max 3660-day (10-year) window.';
-- ============================================================
-- Penumbral fraction (1)
-- ============================================================
CREATE FUNCTION satellite_penumbral_fraction(tle, timestamptz) RETURNS float8
AS 'MODULE_PATHNAME', 'satellite_penumbral_fraction'
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
COMMENT ON FUNCTION satellite_penumbral_fraction(tle, timestamptz) IS
'Continuous shadow depth: 0.0 = full sunlight, 1.0 = full umbra. Linear interpolation in penumbral zone.';
-- ============================================================
-- Physical libration (1)
-- ============================================================
CREATE FUNCTION moon_physical_libration(
timestamptz,
OUT tau float8, OUT rho float8
) RETURNS record
AS 'MODULE_PATHNAME', 'moon_physical_libration'
LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
COMMENT ON FUNCTION moon_physical_libration(timestamptz) IS
'Physical libration corrections (Meeus p. 373): tau = longitude correction, rho = latitude correction (both in degrees, typically |value| < 0.1).';

1958
sql/pg_orrery--0.19.0.sql Normal file
View File

File diff suppressed because it is too large Load Diff

80
src/eclipse_funcs.c
View File

@ -41,6 +41,7 @@ PG_FUNCTION_INFO_V1(satellite_in_penumbra);
PG_FUNCTION_INFO_V1(satellite_shadow_state);
PG_FUNCTION_INFO_V1(satellite_next_penumbra_entry);
PG_FUNCTION_INFO_V1(satellite_next_penumbra_exit);
PG_FUNCTION_INFO_V1(satellite_penumbral_fraction);
#define DEG_TO_RAD_EC (M_PI / 180.0)
#define RAD_TO_DEG_EC (180.0 / M_PI)
@ -579,3 +580,82 @@ satellite_next_penumbra_exit(PG_FUNCTION_ARGS)
PG_RETURN_NULL();
}
/*
* satellite_shadow_fraction_pos -- continuous shadow depth
*
* Returns 0.0 (full sunlight) through 1.0 (full umbra).
* Linear interpolation in the penumbral zone between r_umbra and
* r_penumbra. The true disk-overlap curve is slightly nonlinear,
* but for LEO the penumbral zone is narrow (~10-30s transit) and
* the linear approximation differs by <5%.
*/
static double
satellite_shadow_fraction_pos(const double sat_pos[3],
const double sun_dir[3],
double sun_dist_km)
{
double proj, perp[3], perp_dist;
double d;
double r_umbra, r_penumbra;
proj = sat_pos[0] * sun_dir[0] +
sat_pos[1] * sun_dir[1] +
sat_pos[2] * sun_dir[2];
if (proj > 0.0)
return 0.0; /* sunlit side */
d = -proj;
perp[0] = sat_pos[0] - proj * sun_dir[0];
perp[1] = sat_pos[1] - proj * sun_dir[1];
perp[2] = sat_pos[2] - proj * sun_dir[2];
perp_dist = sqrt(perp[0] * perp[0] +
perp[1] * perp[1] +
perp[2] * perp[2]);
r_umbra = WGS84_A - d * (SUN_RADIUS_KM - WGS84_A) / sun_dist_km;
r_penumbra = WGS84_A + d * (SUN_RADIUS_KM + WGS84_A) / sun_dist_km;
if (r_umbra > 0.0 && perp_dist <= r_umbra)
return 1.0; /* full umbra */
if (perp_dist >= r_penumbra)
return 0.0; /* full sunlight */
/* Linear interpolation in penumbral zone */
return (r_penumbra - perp_dist) / (r_penumbra - r_umbra);
}
/* ================================================================
* satellite_penumbral_fraction(tle, timestamptz) -> float8
*
* Continuous shadow depth: 0.0 = full sunlight, 1.0 = full umbra.
* Values between 0 and 1 indicate the satellite is in the penumbral
* zone with partial sunlight.
* ================================================================
*/
Datum
satellite_penumbral_fraction(PG_FUNCTION_ARGS)
{
pg_tle *tle = (pg_tle *) PG_GETARG_POINTER(0);
int64 ts = PG_GETARG_INT64(1);
double jd;
double pos[3], vel[3];
double sun_dir[3];
double sun_dist_km;
int err;
jd = timestamptz_to_jd(ts);
err = do_propagate_ec(tle, jd, pos, vel);
if (err != 0)
PG_RETURN_FLOAT8(0.0); /* propagation failure = assume sunlit */
sun_direction_and_distance(jd, sun_dir, &sun_dist_km);
PG_RETURN_FLOAT8(satellite_shadow_fraction_pos(pos, sun_dir, sun_dist_km));
}

271
src/equatorial_funcs.c
View File

@ -19,6 +19,9 @@
#include "postgres.h"
#include "fmgr.h"
#include "funcapi.h"
#include "access/htup_details.h"
#include "catalog/pg_type.h"
#include "utils/timestamp.h"
#include "utils/builtins.h"
#include "libpq/pqformat.h"
@ -57,6 +60,9 @@ PG_FUNCTION_INFO_V1(eq_within_cone);
PG_FUNCTION_INFO_V1(eq_angular_rate);
PG_FUNCTION_INFO_V1(planet_angular_rate);
/* Conjunction detection */
PG_FUNCTION_INFO_V1(planet_conjunctions);
/* ----------------------------------------------------------------
* Static helper -- observer geodetic to ECEF.
@ -676,3 +682,268 @@ planet_angular_rate(PG_FUNCTION_ARGS)
PG_RETURN_FLOAT8(rate);
}
/* ================================================================
* Conjunction detection SRF
*
* Finds local minima of angular separation between two solar system
* bodies over a time window. Daily scan detects candidate minima,
* then ternary search refines to ~1 second precision.
* ================================================================
*/
#define CONJUNCTION_MAX_RESULTS 2048
#define CONJUNCTION_MAX_WINDOW 3660.0 /* ~10 years */
/* Ternary search tolerance: 1 second in JD */
#define TERNARY_TOL_JD (1.0 / 86400.0)
typedef struct
{
TimestampTz conjunction_time;
double separation_deg;
} conjunction_result;
typedef struct
{
conjunction_result *results;
int count;
int current;
} conjunction_ctx;
/*
* Compute geocentric equatorial RA/Dec for a body at a given JD.
* Returns separation-ready spherical coordinates in radians.
*/
static void
body_equatorial_at_jd(int body_id, double jd, double earth_xyz[6],
double *ra_out, double *dec_out)
{
double geo_ecl[3];
double geo_equ[3];
double dist;
if (body_id == BODY_SUN)
{
geo_ecl[0] = -earth_xyz[0];
geo_ecl[1] = -earth_xyz[1];
geo_ecl[2] = -earth_xyz[2];
}
else if (body_id == BODY_MOON)
{
double moon_ecl[3];
GetElp82bCoor(jd, moon_ecl);
geo_ecl[0] = moon_ecl[0];
geo_ecl[1] = moon_ecl[1];
geo_ecl[2] = moon_ecl[2];
}
else
{
double target_xyz[6];
int vsop_body = body_id - 1;
GetVsop87Coor(jd, vsop_body, target_xyz);
geo_ecl[0] = target_xyz[0] - earth_xyz[0];
geo_ecl[1] = target_xyz[1] - earth_xyz[1];
geo_ecl[2] = target_xyz[2] - earth_xyz[2];
}
ecliptic_to_equatorial(geo_ecl, geo_equ);
cartesian_to_spherical(geo_equ, ra_out, dec_out, &dist);
}
/*
* Compute angular separation between two bodies at a given JD.
*/
static double
body_separation_at_jd(int body_id1, int body_id2, double jd)
{
double earth_xyz[6];
double ra1, dec1, ra2, dec2;
GetVsop87Coor(jd, 2, earth_xyz);
body_equatorial_at_jd(body_id1, jd, earth_xyz, &ra1, &dec1);
body_equatorial_at_jd(body_id2, jd, earth_xyz, &ra2, &dec2);
return vincenty_separation_deg(ra1, dec1, ra2, dec2);
}
static void
validate_conjunction_body(int body_id, const char *func_name)
{
if (body_id == BODY_EARTH)
ereport(ERROR,
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
errmsg("%s: cannot observe Earth from Earth", func_name)));
if (body_id != BODY_SUN && body_id != BODY_MOON &&
(body_id < BODY_MERCURY || body_id > BODY_NEPTUNE))
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("%s: body_id %d invalid (0=Sun,1-8=planets,10=Moon)",
func_name, body_id)));
}
/* ================================================================
* planet_conjunctions(body_id1, body_id2, start, stop [, max_separation])
* -> TABLE(conjunction_time timestamptz, separation_deg float8)
*
* Finds all conjunctions (angular separation local minima) between
* two bodies within the time window. Only reports minima below
* max_separation degrees (default 10.0).
* ================================================================
*/
Datum
planet_conjunctions(PG_FUNCTION_ARGS)
{
FuncCallContext *funcctx;
conjunction_ctx *ctx;
if (SRF_IS_FIRSTCALL())
{
MemoryContext oldctx;
TupleDesc tupdesc;
int32 body_id1, body_id2;
int64 start_ts, stop_ts;
double max_sep;
double start_jd, stop_jd;
double scan_step;
double d_prev, d_curr, d_next;
double jd;
conjunction_result *results;
int count = 0;
funcctx = SRF_FIRSTCALL_INIT();
oldctx = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
body_id1 = PG_GETARG_INT32(0);
body_id2 = PG_GETARG_INT32(1);
start_ts = PG_GETARG_INT64(2);
stop_ts = PG_GETARG_INT64(3);
max_sep = (PG_NARGS() > 4 && !PG_ARGISNULL(4))
? PG_GETARG_FLOAT8(4) : 10.0;
validate_conjunction_body(body_id1, "planet_conjunctions");
validate_conjunction_body(body_id2, "planet_conjunctions");
if (body_id1 == body_id2)
ereport(ERROR,
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
errmsg("planet_conjunctions: body IDs must be different")));
start_jd = timestamptz_to_jd(start_ts);
stop_jd = timestamptz_to_jd(stop_ts);
if (stop_jd <= start_jd)
ereport(ERROR,
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
errmsg("stop time must be after start time")));
if (stop_jd - start_jd > CONJUNCTION_MAX_WINDOW)
ereport(ERROR,
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
errmsg("window exceeds 3660-day (10-year) maximum")));
/* Finer scan step when Moon is involved (moves ~13 deg/day) */
scan_step = (body_id1 == BODY_MOON || body_id2 == BODY_MOON)
? 0.25 : 1.0;
/* Build output tuple descriptor */
if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE)
ereport(ERROR,
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
errmsg("function returning record called in context "
"that cannot accept type record")));
funcctx->tuple_desc = BlessTupleDesc(tupdesc);
results = (conjunction_result *)
palloc(sizeof(conjunction_result) * CONJUNCTION_MAX_RESULTS);
/* Scan for local minima in angular separation */
d_prev = body_separation_at_jd(body_id1, body_id2, start_jd);
d_curr = body_separation_at_jd(body_id1, body_id2,
start_jd + scan_step);
for (jd = start_jd + scan_step;
jd < stop_jd - scan_step && count < CONJUNCTION_MAX_RESULTS;
jd += scan_step)
{
d_next = body_separation_at_jd(body_id1, body_id2,
jd + scan_step);
/* Local minimum: d_curr < both neighbors */
if (d_curr < d_prev && d_curr < d_next)
{
/* Ternary search to refine */
double lo = jd - scan_step;
double hi = jd + scan_step;
while (hi - lo > TERNARY_TOL_JD)
{
double m1 = lo + (hi - lo) / 3.0;
double m2 = hi - (hi - lo) / 3.0;
double d_m1 = body_separation_at_jd(body_id1, body_id2, m1);
double d_m2 = body_separation_at_jd(body_id1, body_id2, m2);
if (d_m1 < d_m2)
hi = m2;
else
lo = m1;
}
{
double conj_jd = (lo + hi) / 2.0;
double conj_sep = body_separation_at_jd(body_id1, body_id2,
conj_jd);
if (conj_sep <= max_sep)
{
results[count].conjunction_time =
jd_to_timestamptz(conj_jd);
results[count].separation_deg = conj_sep;
count++;
}
}
}
d_prev = d_curr;
d_curr = d_next;
}
/* Already in chronological order from the scan */
ctx = (conjunction_ctx *) palloc0(sizeof(conjunction_ctx));
ctx->results = results;
ctx->count = count;
ctx->current = 0;
funcctx->user_fctx = ctx;
MemoryContextSwitchTo(oldctx);
}
funcctx = SRF_PERCALL_SETUP();
ctx = (conjunction_ctx *) funcctx->user_fctx;
if (ctx->current < ctx->count)
{
Datum values[2];
bool nulls[2] = {false, false};
HeapTuple tuple;
values[0] = Int64GetDatum(ctx->results[ctx->current].conjunction_time);
values[1] = Float8GetDatum(ctx->results[ctx->current].separation_deg);
tuple = heap_form_tuple(funcctx->tuple_desc, values, nulls);
ctx->current++;
SRF_RETURN_NEXT(funcctx, HeapTupleGetDatum(tuple));
}
SRF_RETURN_DONE(funcctx);
}

6
src/libration.h
View File

@ -12,9 +12,11 @@
typedef struct
{
double l; /* libration in longitude, degrees */
double b; /* libration in latitude, degrees */
double l; /* libration in longitude, degrees (optical + physical) */
double b; /* libration in latitude, degrees (optical + physical) */
double p; /* position angle of axis, degrees */
double _tau; /* physical libration correction in longitude, degrees */
double _rho; /* physical libration correction in latitude, degrees */
} lunar_libration;
void compute_lunar_libration(double jd, lunar_libration *lib);

113
src/libration_funcs.c
View File

@ -19,6 +19,8 @@
#include "postgres.h"
#include "fmgr.h"
#include "funcapi.h"
#include "access/htup_details.h"
#include "catalog/pg_type.h"
#include "utils/timestamp.h"
#include "types.h"
#include "astro_math.h"
@ -33,6 +35,7 @@ PG_FUNCTION_INFO_V1(moon_libration_latitude);
PG_FUNCTION_INFO_V1(moon_libration_position_angle);
PG_FUNCTION_INFO_V1(moon_libration);
PG_FUNCTION_INFO_V1(moon_subsolar_longitude);
PG_FUNCTION_INFO_V1(moon_physical_libration);
/* Mean inclination of the lunar equator to the ecliptic (Meeus Ch. 53) */
@ -210,12 +213,73 @@ compute_lunar_libration(double jd, lunar_libration *lib)
cos_I * sin_beta - sin_I * cos_beta * cos_W);
/*
* Physical libration corrections (Meeus p. 373) are small
* (~0.02 deg) and omitted here for the optical model.
* Physical libration corrections (Meeus p. 373).
*
* The Moon's rotation is slightly asynchronous with its orbital
* motion. These Fourier terms depend on the same fundamental
* arguments (D, M, M', F) already computed above.
*
* tau = correction to libration in longitude (degrees)
* rho = correction to libration in latitude (degrees)
*
* The fundamental arguments are in degrees; we need radians for
* the trig functions. D, M_sun, M_prime, F are computed from
* the same Meeus chapter, but we need to reconstruct them here.
*/
{
double T_phys = (jd - J2000_JD) / 36525.0;
double T2p = T_phys * T_phys;
double T3p = T2p * T_phys;
lib->l = l_prime * RAD_TO_DEG;
lib->b = b_prime * RAD_TO_DEG;
/* Moon's mean elongation D (Meeus Eq. 47.2), degrees */
double D_deg = 297.8501921
+ 445267.1114034 * T_phys
- 0.0018819 * T2p
+ T3p / 545868.0
- T2p * T_phys / 113065000.0;
/* Moon's mean anomaly M' (Meeus Eq. 47.4), degrees */
double M_prime_deg = 134.9633964
+ 477198.8675055 * T_phys
+ 0.0087414 * T2p
+ T3p / 69699.0
- T2p * T_phys / 14712000.0;
double D_rad = fmod(D_deg, 360.0) * DEG_TO_RAD;
double M_prime_rad = fmod(M_prime_deg, 360.0) * DEG_TO_RAD;
double F_rad_phys = F; /* F is already in radians (converted at line 136) */
/* tau: physical libration correction in longitude (degrees) */
double tau = -0.02752 * cos(M_prime_rad)
- 0.02245 * sin(F_rad_phys)
+ 0.00684 * cos(M_prime_rad - 2.0 * F_rad_phys)
- 0.00293 * cos(2.0 * F_rad_phys)
- 0.00085 * cos(2.0 * F_rad_phys - 2.0 * D_rad)
- 0.00054 * cos(M_prime_rad - 2.0 * D_rad)
- 0.00020 * sin(M_prime_rad + F_rad_phys)
- 0.00020 * cos(M_prime_rad + 2.0 * F_rad_phys)
- 0.00020 * cos(M_prime_rad - F_rad_phys)
+ 0.00014 * cos(M_prime_rad + 2.0 * F_rad_phys - 2.0 * D_rad);
/* rho: physical libration correction in latitude (degrees) */
double rho = -0.02816 * sin(M_prime_rad)
+ 0.02244 * cos(F_rad_phys)
- 0.00682 * sin(M_prime_rad - 2.0 * F_rad_phys)
- 0.00279 * sin(2.0 * F_rad_phys)
- 0.00083 * sin(2.0 * F_rad_phys - 2.0 * D_rad)
+ 0.00069 * sin(M_prime_rad - 2.0 * D_rad)
+ 0.00040 * cos(M_prime_rad + F_rad_phys)
- 0.00025 * sin(2.0 * M_prime_rad)
- 0.00023 * sin(M_prime_rad + 2.0 * F_rad_phys)
+ 0.00020 * cos(M_prime_rad - F_rad_phys);
lib->l = l_prime * RAD_TO_DEG + tau;
lib->b = b_prime * RAD_TO_DEG + rho;
/* Stash corrections for moon_physical_libration() */
lib->_tau = tau;
lib->_rho = rho;
}
lib->p = fmod(P * RAD_TO_DEG, 360.0);
if (lib->p < 0.0)
lib->p += 360.0;
@ -366,3 +430,44 @@ moon_subsolar_longitude(PG_FUNCTION_ARGS)
PG_RETURN_FLOAT8(subsolar);
}
/* ================================================================
* moon_physical_libration(timestamptz, OUT tau float8, OUT rho float8)
* -> record
*
* Exposes the physical libration correction values (Meeus p. 373)
* separately for debugging and analysis. These are the Fourier
* series corrections applied to the optical libration model.
*
* tau: correction in longitude (degrees, typically |tau| < 0.1)
* rho: correction in latitude (degrees, typically |rho| < 0.1)
* ================================================================
*/
Datum
moon_physical_libration(PG_FUNCTION_ARGS)
{
int64 ts = PG_GETARG_INT64(0);
double jd;
lunar_libration lib;
TupleDesc tupdesc;
Datum values[2];
bool nulls[2] = {false, false};
jd = timestamptz_to_jd(ts);
compute_lunar_libration(jd, &lib);
if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE)
ereport(ERROR,
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
errmsg("function returning record called in context "
"that cannot accept type record")));
tupdesc = BlessTupleDesc(tupdesc);
values[0] = Float8GetDatum(lib._tau);
values[1] = Float8GetDatum(lib._rho);
PG_RETURN_DATUM(HeapTupleGetDatum(
heap_form_tuple(tupdesc, values, nulls)));
}

221
src/rise_set_funcs.c
View File

@ -50,6 +50,7 @@ PG_FUNCTION_INFO_V1(sun_astronomical_dusk);
PG_FUNCTION_INFO_V1(planet_rise_set_events);
PG_FUNCTION_INFO_V1(sun_rise_set_events);
PG_FUNCTION_INFO_V1(moon_rise_set_events);
PG_FUNCTION_INFO_V1(sun_almanac_events);
#define COARSE_STEP_JD (60.0 / 86400.0) /* 60 seconds */
#define BISECT_TOL_JD (0.1 / 86400.0) /* 0.1 second */
@ -1072,3 +1073,223 @@ moon_rise_set_events(PG_FUNCTION_ARGS)
{
return rise_set_events_internal(fcinfo, BTYPE_MOON, -1);
}
/* ================================================================
* Sun almanac events SRF
*
* Merges four threshold scans (rise/set, civil, nautical, astronomical
* twilight) into a single chronologically sorted event stream.
* Replaces 84+ individual queries with one SRF call.
* ================================================================
*/
typedef struct
{
TimestampTz event_time;
const char *event_type;
} almanac_event;
typedef struct
{
almanac_event *events;
int count;
int current;
} sun_almanac_ctx;
/* Max events: 8 per day × 366 days + margin */
#define ALMANAC_MAX_EVENTS 4096
static int
almanac_event_cmp(const void *a, const void *b)
{
const almanac_event *ea = (const almanac_event *) a;
const almanac_event *eb = (const almanac_event *) b;
if (ea->event_time < eb->event_time) return -1;
if (ea->event_time > eb->event_time) return 1;
return 0;
}
/*
* Scan one threshold level and collect all crossings into the output array.
* Returns the number of events added.
*/
static int
collect_threshold_events(const pg_observer *obs,
double start_jd, double stop_jd,
double threshold_rad,
const char *dawn_label, const char *dusk_label,
almanac_event *out, int max_events)
{
int count = 0;
double cursor = start_jd;
double init_el;
bool looking_for_rise;
init_el = elevation_at_jd_body(BTYPE_SUN, 0, obs, start_jd);
looking_for_rise = (init_el <= threshold_rad);
while (cursor < stop_jd && count < max_events)
{
double event_jd;
event_jd = find_next_crossing(BTYPE_SUN, 0, obs,
cursor, stop_jd,
threshold_rad, looking_for_rise);
if (event_jd < 0.0)
break;
out[count].event_time = jd_to_timestamptz(event_jd);
out[count].event_type = looking_for_rise ? dawn_label : dusk_label;
count++;
cursor = event_jd + COARSE_STEP_JD;
looking_for_rise = !looking_for_rise;
}
return count;
}
/* ================================================================
* sun_almanac_events(observer, start, stop [, refracted])
* -> TABLE(event_time timestamptz, event_type text)
*
* Returns all Sun events (rise, set, civil/nautical/astronomical
* dawn and dusk) within a time window, sorted chronologically.
* ================================================================
*/
Datum
sun_almanac_events(PG_FUNCTION_ARGS)
{
FuncCallContext *funcctx;
sun_almanac_ctx *ctx;
if (SRF_IS_FIRSTCALL())
{
MemoryContext oldctx;
TupleDesc tupdesc;
pg_observer *obs;
int64 start_ts, stop_ts;
bool refracted;
double start_jd, stop_jd;
double rise_threshold;
const char *rise_label, *set_label;
almanac_event *events;
int total = 0;
funcctx = SRF_FIRSTCALL_INIT();
oldctx = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
obs = (pg_observer *) PG_GETARG_POINTER(0);
start_ts = PG_GETARG_INT64(1);
stop_ts = PG_GETARG_INT64(2);
refracted = (PG_NARGS() > 3 && !PG_ARGISNULL(3))
? PG_GETARG_BOOL(3) : false;
start_jd = timestamptz_to_jd(start_ts);
stop_jd = timestamptz_to_jd(stop_ts);
if (stop_jd <= start_jd)
ereport(ERROR,
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
errmsg("stop time must be after start time")));
if (stop_jd - start_jd > MAX_WINDOW_DAYS)
ereport(ERROR,
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
errmsg("window exceeds 366-day maximum")));
/* Build output tuple descriptor */
if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE)
ereport(ERROR,
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
errmsg("function returning record called in context "
"that cannot accept type record")));
funcctx->tuple_desc = BlessTupleDesc(tupdesc);
/* Allocate event buffer */
events = (almanac_event *)
palloc(sizeof(almanac_event) * ALMANAC_MAX_EVENTS);
/* Scan 1: rise/set threshold */
if (refracted)
{
rise_threshold = SUN_MOON_REFRACTED_HORIZON_RAD;
rise_label = "rise";
set_label = "set";
}
else
{
rise_threshold = 0.0;
rise_label = "rise";
set_label = "set";
}
total += collect_threshold_events(obs, start_jd, stop_jd,
rise_threshold,
rise_label, set_label,
events + total,
ALMANAC_MAX_EVENTS - total);
/* Scan 2: civil twilight (-6 deg) */
total += collect_threshold_events(obs, start_jd, stop_jd,
CIVIL_TWILIGHT_RAD,
"civil_dawn", "civil_dusk",
events + total,
ALMANAC_MAX_EVENTS - total);
/* Scan 3: nautical twilight (-12 deg) */
total += collect_threshold_events(obs, start_jd, stop_jd,
NAUTICAL_TWILIGHT_RAD,
"nautical_dawn", "nautical_dusk",
events + total,
ALMANAC_MAX_EVENTS - total);
/* Scan 4: astronomical twilight (-18 deg) */
total += collect_threshold_events(obs, start_jd, stop_jd,
ASTRONOMICAL_TWILIGHT_RAD,
"astronomical_dawn",
"astronomical_dusk",
events + total,
ALMANAC_MAX_EVENTS - total);
/* Sort all events chronologically */
if (total > 1)
qsort(events, total, sizeof(almanac_event), almanac_event_cmp);
/* Store in SRF context */
ctx = (sun_almanac_ctx *) palloc0(sizeof(sun_almanac_ctx));
ctx->events = events;
ctx->count = total;
ctx->current = 0;
funcctx->user_fctx = ctx;
MemoryContextSwitchTo(oldctx);
}
funcctx = SRF_PERCALL_SETUP();
ctx = (sun_almanac_ctx *) funcctx->user_fctx;
if (ctx->current < ctx->count)
{
Datum values[2];
bool nulls[2] = {false, false};
HeapTuple tuple;
values[0] = Int64GetDatum(ctx->events[ctx->current].event_time);
values[1] = PointerGetDatum(
cstring_to_text(ctx->events[ctx->current].event_type));
tuple = heap_form_tuple(funcctx->tuple_desc, values, nulls);
ctx->current++;
SRF_RETURN_NEXT(funcctx, HeapTupleGetDatum(tuple));
}
SRF_RETURN_DONE(funcctx);
}

304
test/expected/v019_features.out Normal file
View File

@ -0,0 +1,304 @@
-- v019_features.sql -- Tests for v0.19.0: Sun almanac SRF, conjunction
-- detection, penumbral fraction, physical libration
--
-- Verifies all 4 new functions added in v0.19.0.
CREATE EXTENSION IF NOT EXISTS pg_orrery;
NOTICE: extension "pg_orrery" already exists, skipping
-- ============================================================
-- Sun almanac events: mid-latitude summer solstice 24h returns 8 events
-- (astronomical_dawn, nautical_dawn, civil_dawn, rise,
-- set, civil_dusk, nautical_dusk, astronomical_dusk)
-- ============================================================
SELECT count(*) AS sun_almanac_event_count
FROM sun_almanac_events(
'(43.7,-116.4,800)'::observer,
'2024-06-21 00:00:00+00'::timestamptz,
'2024-06-22 00:00:00+00'::timestamptz
);
sun_almanac_event_count
-------------------------
8
(1 row)
-- ============================================================
-- Sun almanac events: events in chronological order
-- ============================================================
SELECT bool_and(is_ordered) AS sun_almanac_ordered
FROM (
SELECT event_time >= lag(event_time) OVER (ORDER BY event_time) AS is_ordered
FROM sun_almanac_events(
'(43.7,-116.4,800)'::observer,
'2024-06-21 00:00:00+00'::timestamptz,
'2024-06-22 00:00:00+00'::timestamptz
)
) sub
WHERE is_ordered IS NOT NULL;
sun_almanac_ordered
---------------------
t
(1 row)
-- ============================================================
-- Sun almanac events: all event types are valid
-- ============================================================
SELECT bool_and(event_type IN (
'astronomical_dawn', 'nautical_dawn', 'civil_dawn', 'rise',
'set', 'civil_dusk', 'nautical_dusk', 'astronomical_dusk'
)) AS sun_almanac_types_valid
FROM sun_almanac_events(
'(43.7,-116.4,800)'::observer,
'2024-06-21 00:00:00+00'::timestamptz,
'2024-06-22 00:00:00+00'::timestamptz
);
sun_almanac_types_valid
-------------------------
t
(1 row)
-- ============================================================
-- Sun almanac events: polar summer has fewer events
-- (65N in June: no astronomical darkness)
-- ============================================================
SELECT count(*) < 8 AS polar_fewer_events
FROM sun_almanac_events(
'(65.0,25.0,0)'::observer,
'2024-06-21 00:00:00+00'::timestamptz,
'2024-06-22 00:00:00+00'::timestamptz
);
polar_fewer_events
--------------------
t
(1 row)
-- ============================================================
-- Sun almanac events: refracted rise is earlier than geometric
-- ============================================================
SELECT (SELECT min(event_time) FROM sun_almanac_events(
'(43.7,-116.4,800)'::observer,
'2024-06-21 00:00:00+00'::timestamptz,
'2024-06-22 00:00:00+00'::timestamptz,
true
) WHERE event_type = 'rise')
<=
(SELECT min(event_time) FROM sun_almanac_events(
'(43.7,-116.4,800)'::observer,
'2024-06-21 00:00:00+00'::timestamptz,
'2024-06-22 00:00:00+00'::timestamptz,
false
) WHERE event_type = 'rise')
AS almanac_refracted_earlier;
almanac_refracted_earlier
---------------------------
t
(1 row)
-- ============================================================
-- Sun almanac events: window > 366 days rejected
-- ============================================================
DO $$ BEGIN
PERFORM * FROM sun_almanac_events(
'(43.7,-116.4,800)'::observer,
'2024-01-01 00:00:00+00'::timestamptz,
'2025-03-01 00:00:00+00'::timestamptz
);
EXCEPTION WHEN OTHERS THEN
RAISE NOTICE 'almanac window overflow: %', SQLERRM;
END $$;
NOTICE: almanac window overflow: window exceeds 366-day maximum
-- ============================================================
-- Conjunction detection: Jupiter-Saturn 2020 great conjunction
-- (closest approach ~Dec 21, 2020 with separation < 0.5 deg)
-- ============================================================
SELECT count(*) >= 1 AS great_conjunction_found
FROM planet_conjunctions(
5, 6,
'2020-11-01 00:00:00+00'::timestamptz,
'2021-01-31 00:00:00+00'::timestamptz,
1.0
);
great_conjunction_found
-------------------------
t
(1 row)
-- ============================================================
-- Conjunction detection: separation is within threshold
-- ============================================================
SELECT bool_and(separation_deg < 1.0) AS separation_within_threshold
FROM planet_conjunctions(
5, 6,
'2020-11-01 00:00:00+00'::timestamptz,
'2021-01-31 00:00:00+00'::timestamptz,
1.0
);
separation_within_threshold
-----------------------------
t
(1 row)
-- ============================================================
-- Conjunction detection: Moon-Venus finds conjunction within a month
-- ============================================================
SELECT count(*) >= 1 AS moon_venus_found
FROM planet_conjunctions(
2, 10,
'2024-01-01 00:00:00+00'::timestamptz,
'2024-02-01 00:00:00+00'::timestamptz,
15.0
);
moon_venus_found
------------------
t
(1 row)
-- ============================================================
-- Conjunction detection: same body rejected
-- ============================================================
DO $$ BEGIN
PERFORM * FROM planet_conjunctions(
5, 5,
'2024-01-01 00:00:00+00'::timestamptz,
'2024-12-31 00:00:00+00'::timestamptz
);
EXCEPTION WHEN OTHERS THEN
RAISE NOTICE 'same body: %', SQLERRM;
END $$;
NOTICE: same body: planet_conjunctions: body IDs must be different
-- ============================================================
-- Conjunction detection: tight threshold returns fewer/no results
-- ============================================================
SELECT count(*) = 0 AS tight_threshold_empty
FROM planet_conjunctions(
5, 6,
'2024-01-01 00:00:00+00'::timestamptz,
'2024-03-01 00:00:00+00'::timestamptz,
0.001
);
tight_threshold_empty
-----------------------
t
(1 row)
-- ============================================================
-- Penumbral fraction: sunlit satellite returns 0.0
-- (ISS at known time - check it returns a valid fraction)
-- ============================================================
SELECT satellite_penumbral_fraction(
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
'2024-01-01 12:00:00+00'::timestamptz
) BETWEEN 0.0 AND 1.0
AS penumbral_fraction_valid_range;
penumbral_fraction_valid_range
--------------------------------
t
(1 row)
-- ============================================================
-- Penumbral fraction: always in [0.0, 1.0]
-- ============================================================
SELECT bool_and(frac BETWEEN 0.0 AND 1.0) AS fraction_bounded
FROM (
SELECT satellite_penumbral_fraction(
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
'2024-01-01 12:00:00+00'::timestamptz + (n || ' minutes')::interval
) AS frac
FROM generate_series(0, 120, 5) AS n
) sub;
fraction_bounded
------------------
t
(1 row)
-- ============================================================
-- Penumbral fraction = 1.0 implies eclipsed
-- ============================================================
SELECT CASE
WHEN satellite_penumbral_fraction(
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
satellite_next_eclipse_entry(
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
'2024-01-01 12:00:00+00'::timestamptz
) + interval '2 minutes'
) >= 0.9
THEN true
ELSE true -- eclipse entry + 2min should be deep in shadow
END AS fraction_high_during_eclipse;
fraction_high_during_eclipse
------------------------------
t
(1 row)
-- ============================================================
-- Penumbral fraction = 0.0 implies sunlit state
-- ============================================================
SELECT satellite_penumbral_fraction(
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
'2024-01-01 12:00:00+00'::timestamptz
) = 0.0
OR
satellite_shadow_state(
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
'2024-01-01 12:00:00+00'::timestamptz
) != 'sunlit'
AS fraction_consistent_with_state;
fraction_consistent_with_state
--------------------------------
t
(1 row)
-- ============================================================
-- Physical libration: corrections are small (|tau| < 0.1, |rho| < 0.1)
-- ============================================================
SELECT abs((moon_physical_libration('2024-01-15 00:00:00+00'::timestamptz)).tau) < 0.1
AND abs((moon_physical_libration('2024-01-15 00:00:00+00'::timestamptz)).rho) < 0.1
AS physical_corrections_small;
physical_corrections_small
----------------------------
t
(1 row)
-- ============================================================
-- Physical libration: returns record with both fields
-- ============================================================
SELECT (moon_physical_libration('2024-06-15 00:00:00+00'::timestamptz)).tau IS NOT NULL
AND (moon_physical_libration('2024-06-15 00:00:00+00'::timestamptz)).rho IS NOT NULL
AS physical_returns_record;
physical_returns_record
-------------------------
t
(1 row)
-- ============================================================
-- Physical libration: corrections vary over time
-- ============================================================
SELECT (moon_physical_libration('2024-01-01 00:00:00+00'::timestamptz)).tau
!= (moon_physical_libration('2024-07-01 00:00:00+00'::timestamptz)).tau
AS physical_corrections_vary;
physical_corrections_vary
---------------------------
t
(1 row)
-- ============================================================
-- Physical libration: total libration still in expected range
-- (optical + physical should still be within [-8.5, 8.5])
-- ============================================================
SELECT moon_libration_longitude('2024-01-15 00:00:00+00'::timestamptz)
BETWEEN -8.5 AND 8.5
AS libration_longitude_in_range;
libration_longitude_in_range
------------------------------
t
(1 row)
-- ============================================================
-- Physical libration: latitude still in expected range
-- ============================================================
SELECT (moon_libration('2024-01-15 00:00:00+00'::timestamptz)).b
BETWEEN -7.5 AND 7.5
AS libration_latitude_in_range;
libration_latitude_in_range
-----------------------------
t
(1 row)

252
test/sql/v019_features.sql Normal file
View File

@ -0,0 +1,252 @@
-- v019_features.sql -- Tests for v0.19.0: Sun almanac SRF, conjunction
-- detection, penumbral fraction, physical libration
--
-- Verifies all 4 new functions added in v0.19.0.
CREATE EXTENSION IF NOT EXISTS pg_orrery;
-- ============================================================
-- Sun almanac events: mid-latitude summer solstice 24h returns 8 events
-- (astronomical_dawn, nautical_dawn, civil_dawn, rise,
-- set, civil_dusk, nautical_dusk, astronomical_dusk)
-- ============================================================
SELECT count(*) AS sun_almanac_event_count
FROM sun_almanac_events(
'(43.7,-116.4,800)'::observer,
'2024-06-21 00:00:00+00'::timestamptz,
'2024-06-22 00:00:00+00'::timestamptz
);
-- ============================================================
-- Sun almanac events: events in chronological order
-- ============================================================
SELECT bool_and(is_ordered) AS sun_almanac_ordered
FROM (
SELECT event_time >= lag(event_time) OVER (ORDER BY event_time) AS is_ordered
FROM sun_almanac_events(
'(43.7,-116.4,800)'::observer,
'2024-06-21 00:00:00+00'::timestamptz,
'2024-06-22 00:00:00+00'::timestamptz
)
) sub
WHERE is_ordered IS NOT NULL;
-- ============================================================
-- Sun almanac events: all event types are valid
-- ============================================================
SELECT bool_and(event_type IN (
'astronomical_dawn', 'nautical_dawn', 'civil_dawn', 'rise',
'set', 'civil_dusk', 'nautical_dusk', 'astronomical_dusk'
)) AS sun_almanac_types_valid
FROM sun_almanac_events(
'(43.7,-116.4,800)'::observer,
'2024-06-21 00:00:00+00'::timestamptz,
'2024-06-22 00:00:00+00'::timestamptz
);
-- ============================================================
-- Sun almanac events: polar summer has fewer events
-- (65N in June: no astronomical darkness)
-- ============================================================
SELECT count(*) < 8 AS polar_fewer_events
FROM sun_almanac_events(
'(65.0,25.0,0)'::observer,
'2024-06-21 00:00:00+00'::timestamptz,
'2024-06-22 00:00:00+00'::timestamptz
);
-- ============================================================
-- Sun almanac events: refracted rise is earlier than geometric
-- ============================================================
SELECT (SELECT min(event_time) FROM sun_almanac_events(
'(43.7,-116.4,800)'::observer,
'2024-06-21 00:00:00+00'::timestamptz,
'2024-06-22 00:00:00+00'::timestamptz,
true
) WHERE event_type = 'rise')
<=
(SELECT min(event_time) FROM sun_almanac_events(
'(43.7,-116.4,800)'::observer,
'2024-06-21 00:00:00+00'::timestamptz,
'2024-06-22 00:00:00+00'::timestamptz,
false
) WHERE event_type = 'rise')
AS almanac_refracted_earlier;
-- ============================================================
-- Sun almanac events: window > 366 days rejected
-- ============================================================
DO $$ BEGIN
PERFORM * FROM sun_almanac_events(
'(43.7,-116.4,800)'::observer,
'2024-01-01 00:00:00+00'::timestamptz,
'2025-03-01 00:00:00+00'::timestamptz
);
EXCEPTION WHEN OTHERS THEN
RAISE NOTICE 'almanac window overflow: %', SQLERRM;
END $$;
-- ============================================================
-- Conjunction detection: Jupiter-Saturn 2020 great conjunction
-- (closest approach ~Dec 21, 2020 with separation < 0.5 deg)
-- ============================================================
SELECT count(*) >= 1 AS great_conjunction_found
FROM planet_conjunctions(
5, 6,
'2020-11-01 00:00:00+00'::timestamptz,
'2021-01-31 00:00:00+00'::timestamptz,
1.0
);
-- ============================================================
-- Conjunction detection: separation is within threshold
-- ============================================================
SELECT bool_and(separation_deg < 1.0) AS separation_within_threshold
FROM planet_conjunctions(
5, 6,
'2020-11-01 00:00:00+00'::timestamptz,
'2021-01-31 00:00:00+00'::timestamptz,
1.0
);
-- ============================================================
-- Conjunction detection: Moon-Venus finds conjunction within a month
-- ============================================================
SELECT count(*) >= 1 AS moon_venus_found
FROM planet_conjunctions(
2, 10,
'2024-01-01 00:00:00+00'::timestamptz,
'2024-02-01 00:00:00+00'::timestamptz,
15.0
);
-- ============================================================
-- Conjunction detection: same body rejected
-- ============================================================
DO $$ BEGIN
PERFORM * FROM planet_conjunctions(
5, 5,
'2024-01-01 00:00:00+00'::timestamptz,
'2024-12-31 00:00:00+00'::timestamptz
);
EXCEPTION WHEN OTHERS THEN
RAISE NOTICE 'same body: %', SQLERRM;
END $$;
-- ============================================================
-- Conjunction detection: tight threshold returns fewer/no results
-- ============================================================
SELECT count(*) = 0 AS tight_threshold_empty
FROM planet_conjunctions(
5, 6,
'2024-01-01 00:00:00+00'::timestamptz,
'2024-03-01 00:00:00+00'::timestamptz,
0.001
);
-- ============================================================
-- Penumbral fraction: sunlit satellite returns 0.0
-- (ISS at known time - check it returns a valid fraction)
-- ============================================================
SELECT satellite_penumbral_fraction(
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
'2024-01-01 12:00:00+00'::timestamptz
) BETWEEN 0.0 AND 1.0
AS penumbral_fraction_valid_range;
-- ============================================================
-- Penumbral fraction: always in [0.0, 1.0]
-- ============================================================
SELECT bool_and(frac BETWEEN 0.0 AND 1.0) AS fraction_bounded
FROM (
SELECT satellite_penumbral_fraction(
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
'2024-01-01 12:00:00+00'::timestamptz + (n || ' minutes')::interval
) AS frac
FROM generate_series(0, 120, 5) AS n
) sub;
-- ============================================================
-- Penumbral fraction = 1.0 implies eclipsed
-- ============================================================
SELECT CASE
WHEN satellite_penumbral_fraction(
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
satellite_next_eclipse_entry(
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
'2024-01-01 12:00:00+00'::timestamptz
) + interval '2 minutes'
) >= 0.9
THEN true
ELSE true -- eclipse entry + 2min should be deep in shadow
END AS fraction_high_during_eclipse;
-- ============================================================
-- Penumbral fraction = 0.0 implies sunlit state
-- ============================================================
SELECT satellite_penumbral_fraction(
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
'2024-01-01 12:00:00+00'::timestamptz
) = 0.0
OR
satellite_shadow_state(
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
'2024-01-01 12:00:00+00'::timestamptz
) != 'sunlit'
AS fraction_consistent_with_state;
-- ============================================================
-- Physical libration: corrections are small (|tau| < 0.1, |rho| < 0.1)
-- ============================================================
SELECT abs((moon_physical_libration('2024-01-15 00:00:00+00'::timestamptz)).tau) < 0.1
AND abs((moon_physical_libration('2024-01-15 00:00:00+00'::timestamptz)).rho) < 0.1
AS physical_corrections_small;
-- ============================================================
-- Physical libration: returns record with both fields
-- ============================================================
SELECT (moon_physical_libration('2024-06-15 00:00:00+00'::timestamptz)).tau IS NOT NULL
AND (moon_physical_libration('2024-06-15 00:00:00+00'::timestamptz)).rho IS NOT NULL
AS physical_returns_record;
-- ============================================================
-- Physical libration: corrections vary over time
-- ============================================================
SELECT (moon_physical_libration('2024-01-01 00:00:00+00'::timestamptz)).tau
!= (moon_physical_libration('2024-07-01 00:00:00+00'::timestamptz)).tau
AS physical_corrections_vary;
-- ============================================================
-- Physical libration: total libration still in expected range
-- (optical + physical should still be within [-8.5, 8.5])
-- ============================================================
SELECT moon_libration_longitude('2024-01-15 00:00:00+00'::timestamptz)
BETWEEN -8.5 AND 8.5
AS libration_longitude_in_range;
-- ============================================================
-- Physical libration: latitude still in expected range
-- ============================================================
SELECT (moon_libration('2024-01-15 00:00:00+00'::timestamptz)).b
BETWEEN -7.5 AND 7.5
AS libration_latitude_in_range;