Add v0.16.0: twilight dawn/dusk, lunar phase, planet apparent magnitude

Twilight: 6 functions (civil/nautical/astronomical × dawn/dusk) reusing the existing find_next_crossing() bisection search with Sun depression angle thresholds (-6°, -12°, -18°). Returns NULL for polar regions where the threshold is never reached. Lunar phase: 4 functions computing Sun-Earth-Moon geometry from VSOP87 + ELP2000-82B. Phase angle [0,360) via elongation + cross product z-component for waxing/waning discrimination. 8 named phases in 45° bins. Moon age approximated from phase angle and mean synodic month. Planet magnitude: Mallama & Hilton (2018) polynomial model with VSOP87 heliocentric distances and phase angle via law of cosines. All 7 planets (Mercury-Neptune, excluding Earth). Saturn ring tilt not modeled. 151 → 162 SQL objects. 26 → 27 test suites, all passing.
2026-02-26 12:42:01 -07:00 · 2026-02-26 12:42:01 -07:00 · 46c8a30575
commit 46c8a30575
parent e670ed7ed1
10 changed files with 2747 additions and 14 deletions
--- a/CLAUDE.md
+++ b/CLAUDE.md
@ -1,9 +1,9 @@
 # pg_orrery — A Database Orrery for PostgreSQL

 ## What This Is
-A database orrery — celestial mechanics types and functions for PostgreSQL. Native C extension using PGXS, 151 SQL objects (135 user-visible functions + 16 GiST support), 9 custom types, covering satellites (SGP4/SDP4), planets (VSOP87 + optional JPL DE441), Moon (ELP2000-82B), 19 planetary moons (L1.2/TASS17/GUST86/MarsSat), stars (with proper motion and annual parallax), comets, asteroids (MPC catalog), Jupiter radio bursts, interplanetary Lambert transfers, equatorial RA/Dec coordinates with GiST-indexed angular separation, atmospheric refraction, annual stellar aberration, light-time correction, rise/set prediction (geometric + refracted) with status diagnostics, and IAU constellation identification with full name lookup (Roman 1987).
+A database orrery — celestial mechanics types and functions for PostgreSQL. Native C extension using PGXS, 162 SQL objects (146 user-visible functions + 16 GiST support), 9 custom types, covering satellites (SGP4/SDP4), planets (VSOP87 + optional JPL DE441), Moon (ELP2000-82B), 19 planetary moons (L1.2/TASS17/GUST86/MarsSat), stars (with proper motion and annual parallax), comets, asteroids (MPC catalog), Jupiter radio bursts, interplanetary Lambert transfers, equatorial RA/Dec coordinates with GiST-indexed angular separation, atmospheric refraction, annual stellar aberration, light-time correction, rise/set prediction (geometric + refracted) with status diagnostics, IAU constellation identification with full name lookup (Roman 1987), twilight dawn/dusk (civil/nautical/astronomical), lunar phase (angle, illumination, name, age), and planet apparent magnitude (Mallama & Hilton 2018).

-**Current version:** 0.15.0
+**Current version:** 0.16.0
 **Repository:** https://git.supported.systems/warehack.ing/pg_orrery
 **Documentation:** https://pg-orrery.warehack.ing

@ -11,7 +11,7 @@ A database orrery — celestial mechanics types and functions for PostgreSQL. Na
 ```bash
 make PG_CONFIG=/usr/bin/pg_config           # Compile with PGXS
 sudo make install PG_CONFIG=/usr/bin/pg_config  # Install extension
-make installcheck PG_CONFIG=/usr/bin/pg_config  # Run 26 regression test suites
+make installcheck PG_CONFIG=/usr/bin/pg_config  # Run 27 regression test suites
 ```

 Requires: PostgreSQL 17 development headers, GCC, Make.
@ -27,7 +27,7 @@ Image: `git.supported.systems/warehack.ing/pg_orrery:pg17`

 ## Project Layout
 ```
-pg_orrery.control                # Extension metadata (version 0.15.0)
+pg_orrery.control                # Extension metadata (version 0.16.0)
 Makefile                        # PGXS build + Docker targets
 sql/
  pg_orrery--0.1.0.sql           # v0.1.0: satellite types/functions/operators
@ -45,6 +45,7 @@ sql/
  pg_orrery--0.13.0.sql          # v0.13.0: nutation, make_equatorial, rise/set (141 objects)
  pg_orrery--0.14.0.sql          # v0.14.0: refracted rise/set, constellation ID (147 objects)
  pg_orrery--0.15.0.sql          # v0.15.0: constellation full name, rise/set status (151 objects)
+  pg_orrery--0.16.0.sql          # v0.16.0: twilight, lunar phase, planet magnitude (162 objects)
  pg_orrery--0.1.0--0.2.0.sql    # Migration: v0.1.0 → v0.2.0 (adds solar system)
  pg_orrery--0.2.0--0.3.0.sql    # Migration: v0.2.0 → v0.3.0 (adds DE ephemeris)
  pg_orrery--0.3.0--0.4.0.sql    # Migration: v0.3.0 → v0.4.0
@ -59,6 +60,7 @@ sql/
  pg_orrery--0.12.0--0.13.0.sql  # Migration: v0.12.0 → v0.13.0 (nutation, make_equatorial, rise/set)
  pg_orrery--0.13.0--0.14.0.sql  # Migration: v0.13.0 → v0.14.0 (refracted rise/set, constellation ID)
  pg_orrery--0.14.0--0.15.0.sql  # Migration: v0.14.0 → v0.15.0 (constellation full name, rise/set status)
+  pg_orrery--0.15.0--0.16.0.sql  # Migration: v0.15.0 → v0.16.0 (twilight, lunar phase, planet magnitude)
 src/
  pg_orrery.c                    # PG_MODULE_MAGIC + _PG_init() (GUC registration)
  types.h                       # All struct definitions + constants + DE body ID mapping
@ -85,9 +87,11 @@ src/
  orbital_elements_type.c       # orbital_elements type, MPC parser, small_body_observe/equatorial/apparent()
  equatorial_funcs.c            # equatorial type I/O, accessors, satellite/planet/sun/moon RA/Dec
  refraction_funcs.c            # atmospheric_refraction(), _ext(), topo_elevation_apparent()
-  rise_set_funcs.c              # planet/sun/moon rise/set (geometric + refracted)
+  rise_set_funcs.c              # planet/sun/moon rise/set (geometric + refracted) + twilight dawn/dusk
  constellation_data.h / .c     # Roman (1987) IAU boundary table (CDS VI/42, 357 segments)
  constellation_funcs.c         # constellation() from equatorial or RA/Dec
+  lunar_phase_funcs.c           # moon_phase_angle(), moon_illumination(), moon_phase_name(), moon_age()
+  magnitude_funcs.c             # planet_magnitude() (Mallama & Hilton 2018)
  l12.c / l12.h                 # L1.2 Galilean moon theory (Lieske 1998)
  tass17.c / tass17.h           # TASS 1.7 Saturn moon theory (Vienne & Duriez 1995)
  gust86.c / gust86.h           # GUST86 Uranus moon theory (Laskar & Jacobson 1987)
@ -112,7 +116,7 @@ src/
    PROVENANCE.md               # Vendoring decision, modifications, verification
    LICENSE                     # MIT license (Bill Gray / Project Pluto)
 test/
-  sql/                          # 26 regression test suites
+  sql/                          # 27 regression test suites
  expected/                     # Expected output
  data/vallado_518.json         # 518 Vallado test vectors (AIAA 2006-6753-Rev1)
 docs/
@ -139,7 +143,7 @@ All types are fixed-size, `STORAGE = plain`, `ALIGNMENT = double`. No TOAST over
 | `orbital_elements` | 72 | Classical Keplerian elements for comets/asteroids (epoch, q, e, inc, omega, Omega, tp, H, G) |
 | `equatorial` | 24 | Apparent RA (hours), Dec (degrees), distance (km) — of date |

-## Function Domains (151 SQL objects)
+## Function Domains (162 SQL objects)

 | Domain | Theory | Key Functions | Count |
 |--------|--------|---------------|-------|
@ -157,6 +161,9 @@ All types are fixed-size, `STORAGE = plain`, `ALIGNMENT = double`. No TOAST over
 | GiST index (TLE) | Altitude-band approximation | `&&` (overlap), `<->` (distance) | 8 |
 | GiST index (equatorial) | Spherical bounding box | `<->` (KNN ordering) | 8 |
 | Rise/set | Bisection (60s scan) | `planet_next_rise()`, `sun_next_rise_refracted()`, `*_rise_set_status()` | 15 |
+| Twilight | Sun depression angles | `sun_civil_dawn()`, `sun_nautical_dusk()`, `sun_astronomical_dawn()` | 6 |
+| Lunar phase | VSOP87 + ELP2000-82B geometry | `moon_phase_angle()`, `moon_illumination()`, `moon_phase_name()`, `moon_age()` | 4 |
+| Planet magnitude | Mallama & Hilton (2018) | `planet_magnitude()` | 1 |
 | Constellation | Roman (1987) CDS VI/42 | `constellation()`, `constellation_full_name()` | 3 |
 | Diagnostics | -- | `pg_orrery_ephemeris_info()` | 1 |

@ -291,7 +298,7 @@ All numerical logic is byte-identical to upstream. Verified against 518 Vallado

 ## Testing

-26 regression test suites via `make installcheck`:
+27 regression test suites via `make installcheck`:

 | Suite | What it tests |
 |-------|--------------|
@ -321,10 +328,11 @@ All numerical logic is byte-identical to upstream. Verified against 518 Vallado
 | rise_set | Planet/Sun/Moon rise/set (geometric + refracted), circumpolar, polar night |
 | constellation | Roman (1987) boundary lookup, known stars, solar system objects, edge cases |
 | v015_features | constellation_full_name lookup, rise_set_status diagnostics (circumpolar/never_rises) |
+| v016_features | Twilight ordering/offset/polar, lunar phase at known events, planet magnitude ranges/errors |

 ### PG Version Matrix

-Test all 26 regression suites + DE reader unit test across PostgreSQL 14-18 using Docker:
+Test all 27 regression suites + DE reader unit test across PostgreSQL 14-18 using Docker:

 ```bash
 make test-matrix                              # Full matrix (PG 14-18)
@ -350,7 +358,7 @@ Logs saved to `test/matrix-logs/pg${ver}.log`. The script reuses the Dockerfile

 Starlight docs at `docs/` — 44+ MDX pages covering all domains.

-Sections: Getting Started, Guides (9 domain walkthroughs incl. DE ephemeris), Workflow Translation (Skyfield/Horizons/GMAT/Radio Jupiter Pro comparisons), Reference (all 151 SQL objects incl. DE variants, equatorial GiST, refraction, rise/set, constellation), Architecture (Hamilton's principles, constant custody, observation pipeline), Performance (benchmarks).
+Sections: Getting Started, Guides (9 domain walkthroughs incl. DE ephemeris), Workflow Translation (Skyfield/Horizons/GMAT/Radio Jupiter Pro comparisons), Reference (all 162 SQL objects incl. DE variants, equatorial GiST, refraction, rise/set, constellation, twilight, lunar phase, planet magnitude), Architecture (Hamilton's principles, constant custody, observation pipeline), Performance (benchmarks).

 ### Local Development
 ```bash
--- a/9
+++ b/9
@ -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_).
--- a/pg_orrery.control
+++ b/pg_orrery.control
@ -1,4 +1,4 @@
 comment = 'A database orrery — celestial mechanics types and functions for PostgreSQL'
-default_version = '0.15.0'
+default_version = '0.16.0'
 module_pathname = '$libdir/pg_orrery'
 relocatable = true
--- a/sql/pg_orrery--0.15.0--0.16.0.sql
+++ b/sql/pg_orrery--0.15.0--0.16.0.sql
@ -0,0 +1,79 @@
+-- pg_orrery 0.15.0 -> 0.16.0: twilight, lunar phase, planet magnitude
+
+-- ============================================================
+-- Twilight functions (6)
+-- ============================================================
+
+CREATE FUNCTION sun_civil_dawn(observer, timestamptz) RETURNS timestamptz
+    AS 'MODULE_PATHNAME', 'sun_civil_dawn'
+    LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_civil_dawn(observer, timestamptz) IS
+    'Next civil dawn (Sun crosses -6 deg rising). Outdoor activities without artificial light.';
+
+CREATE FUNCTION sun_civil_dusk(observer, timestamptz) RETURNS timestamptz
+    AS 'MODULE_PATHNAME', 'sun_civil_dusk'
+    LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_civil_dusk(observer, timestamptz) IS
+    'Next civil dusk (Sun crosses -6 deg setting). Artificial light needed.';
+
+CREATE FUNCTION sun_nautical_dawn(observer, timestamptz) RETURNS timestamptz
+    AS 'MODULE_PATHNAME', 'sun_nautical_dawn'
+    LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_nautical_dawn(observer, timestamptz) IS
+    'Next nautical dawn (Sun crosses -12 deg rising). Horizon visible at sea.';
+
+CREATE FUNCTION sun_nautical_dusk(observer, timestamptz) RETURNS timestamptz
+    AS 'MODULE_PATHNAME', 'sun_nautical_dusk'
+    LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_nautical_dusk(observer, timestamptz) IS
+    'Next nautical dusk (Sun crosses -12 deg setting). Horizon no longer visible at sea.';
+
+CREATE FUNCTION sun_astronomical_dawn(observer, timestamptz) RETURNS timestamptz
+    AS 'MODULE_PATHNAME', 'sun_astronomical_dawn'
+    LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_astronomical_dawn(observer, timestamptz) IS
+    'Next astronomical dawn (Sun crosses -18 deg rising). Sky was fully dark.';
+
+CREATE FUNCTION sun_astronomical_dusk(observer, timestamptz) RETURNS timestamptz
+    AS 'MODULE_PATHNAME', 'sun_astronomical_dusk'
+    LANGUAGE C STABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION sun_astronomical_dusk(observer, timestamptz) IS
+    'Next astronomical dusk (Sun crosses -18 deg setting). Sky becomes fully dark.';
+
+-- ============================================================
+-- Lunar phase functions (4)
+-- ============================================================
+
+CREATE FUNCTION moon_phase_angle(timestamptz) RETURNS float8
+    AS 'MODULE_PATHNAME', 'moon_phase_angle'
+    LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION moon_phase_angle(timestamptz) IS
+    'Sun-Earth-Moon phase angle in degrees [0,360). 0=new, 90=first quarter, 180=full, 270=last quarter.';
+
+CREATE FUNCTION moon_illumination(timestamptz) RETURNS float8
+    AS 'MODULE_PATHNAME', 'moon_illumination'
+    LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION moon_illumination(timestamptz) IS
+    'Illuminated fraction of the Moon disk [0.0, 1.0].';
+
+CREATE FUNCTION moon_phase_name(timestamptz) RETURNS text
+    AS 'MODULE_PATHNAME', 'moon_phase_name'
+    LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION moon_phase_name(timestamptz) IS
+    'Moon phase name: new_moon, waxing_crescent, first_quarter, waxing_gibbous, full_moon, waning_gibbous, last_quarter, waning_crescent.';
+
+CREATE FUNCTION moon_age(timestamptz) RETURNS float8
+    AS 'MODULE_PATHNAME', 'moon_age'
+    LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION moon_age(timestamptz) IS
+    'Days since last new moon [0, ~29.53), approximated from phase angle.';
+
+-- ============================================================
+-- Planet magnitude (1)
+-- ============================================================
+
+CREATE FUNCTION planet_magnitude(int4, timestamptz) RETURNS float8
+    AS 'MODULE_PATHNAME', 'planet_magnitude'
+    LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
+COMMENT ON FUNCTION planet_magnitude(int4, timestamptz) IS
+    'Apparent visual magnitude of a planet (Mallama & Hilton 2018). Body IDs 1-8. Saturn ring tilt not modeled.';
--- a/sql/pg_orrery--0.16.0.sql
+++ b/sql/pg_orrery--0.16.0.sql
--- a/src/lunar_phase_funcs.c
+++ b/src/lunar_phase_funcs.c
@ -0,0 +1,236 @@
+/*
+ * lunar_phase_funcs.c -- Lunar phase calculations
+ *
+ * Phase angle, illumination, phase name, and lunar age from
+ * VSOP87 Sun + ELP2000-82B Moon positions.
+ *
+ * The Sun-Earth-Moon geometry determines phase:
+ *   1. Earth heliocentric (VSOP87) -> negate -> Sun geocentric
+ *   2. Moon geocentric ecliptic (ELP2000-82B)
+ *   3. Elongation = angle between geocentric Sun and Moon
+ *   4. Cross product z-component distinguishes waxing from waning
+ */
+
+#include "postgres.h"
+#include "fmgr.h"
+#include "utils/timestamp.h"
+#include "utils/builtins.h"
+#include "types.h"
+#include "astro_math.h"
+#include "vsop87.h"
+#include "elp82b.h"
+#include <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);
+}
--- a/src/magnitude_funcs.c
+++ b/src/magnitude_funcs.c
@ -0,0 +1,144 @@
+/*
+ * magnitude_funcs.c -- Planet apparent visual magnitude
+ *
+ * Uses the Mallama & Hilton (2018) magnitude model with
+ * VSOP87 positions for distances and phase angles.
+ *
+ * Reference: Mallama & Hilton, "Computing Apparent Planetary
+ * Magnitudes for The Astronomical Almanac", A&C vol. 25, 2018.
+ */
+
+#include "postgres.h"
+#include "fmgr.h"
+#include "utils/timestamp.h"
+#include "types.h"
+#include "astro_math.h"
+#include "vsop87.h"
+#include <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);
+}
--- a/src/rise_set_funcs.c
+++ b/src/rise_set_funcs.c
@ -38,6 +38,12 @@ PG_FUNCTION_INFO_V1(moon_next_set_refracted);
 PG_FUNCTION_INFO_V1(sun_rise_set_status);
 PG_FUNCTION_INFO_V1(moon_rise_set_status);
 PG_FUNCTION_INFO_V1(planet_rise_set_status);
+PG_FUNCTION_INFO_V1(sun_civil_dawn);
+PG_FUNCTION_INFO_V1(sun_civil_dusk);
+PG_FUNCTION_INFO_V1(sun_nautical_dawn);
+PG_FUNCTION_INFO_V1(sun_nautical_dusk);
+PG_FUNCTION_INFO_V1(sun_astronomical_dawn);
+PG_FUNCTION_INFO_V1(sun_astronomical_dusk);

 #define COARSE_STEP_JD      (60.0 / 86400.0)     /* 60 seconds */
 #define BISECT_TOL_JD       (0.1 / 86400.0)      /* 0.1 second */
@ -65,6 +71,11 @@ PG_FUNCTION_INFO_V1(planet_rise_set_status);
 */
 #define REFRACTION_ONLY_HORIZON_RAD     (-0.00993)  /* -0.569 deg */

+/* Twilight depression angles (geometric Sun center below horizon) */
+#define CIVIL_TWILIGHT_RAD          (-0.10472)  /* -6.0 deg */
+#define NAUTICAL_TWILIGHT_RAD       (-0.20944)  /* -12.0 deg */
+#define ASTRONOMICAL_TWILIGHT_RAD   (-0.30416)  /* -18.0 deg */
+

 /* ----------------------------------------------------------------
 * elevation_at_jd_body -- compute topocentric elevation for a body
@ -684,3 +695,177 @@ planet_rise_set_status(PG_FUNCTION_ARGS)

    PG_RETURN_TEXT_P(cstring_to_text(status));
 }
+
+
+/* ================================================================
+ * sun_civil_dawn(observer, timestamptz) -> timestamptz
+ *
+ * Returns the next time civil twilight begins (Sun crosses -6 deg
+ * heading upward).  Civil twilight = enough light for outdoor
+ * activities without artificial lighting.
+ * ================================================================
+ */
+Datum
+sun_civil_dawn(PG_FUNCTION_ARGS)
+{
+    pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
+    int64        ts  = PG_GETARG_INT64(1);
+    double       start_jd, stop_jd, result_jd;
+
+    start_jd = timestamptz_to_jd(ts);
+    stop_jd  = start_jd + DEFAULT_WINDOW_DAYS;
+
+    result_jd = find_next_crossing(BTYPE_SUN, 0, obs,
+                                   start_jd, stop_jd,
+                                   CIVIL_TWILIGHT_RAD, true);
+
+    if (result_jd < 0.0)
+        PG_RETURN_NULL();
+
+    PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
+}
+
+
+/* ================================================================
+ * sun_civil_dusk(observer, timestamptz) -> timestamptz
+ *
+ * Returns the next time civil twilight ends (Sun crosses -6 deg
+ * heading downward).  After civil dusk, outdoor activities require
+ * artificial lighting.
+ * ================================================================
+ */
+Datum
+sun_civil_dusk(PG_FUNCTION_ARGS)
+{
+    pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
+    int64        ts  = PG_GETARG_INT64(1);
+    double       start_jd, stop_jd, result_jd;
+
+    start_jd = timestamptz_to_jd(ts);
+    stop_jd  = start_jd + DEFAULT_WINDOW_DAYS;
+
+    result_jd = find_next_crossing(BTYPE_SUN, 0, obs,
+                                   start_jd, stop_jd,
+                                   CIVIL_TWILIGHT_RAD, false);
+
+    if (result_jd < 0.0)
+        PG_RETURN_NULL();
+
+    PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
+}
+
+
+/* ================================================================
+ * sun_nautical_dawn(observer, timestamptz) -> timestamptz
+ *
+ * Returns the next time nautical twilight begins (Sun crosses -12 deg
+ * heading upward).  At nautical dawn the horizon becomes visible at
+ * sea and bright stars are still visible for navigation.
+ * ================================================================
+ */
+Datum
+sun_nautical_dawn(PG_FUNCTION_ARGS)
+{
+    pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
+    int64        ts  = PG_GETARG_INT64(1);
+    double       start_jd, stop_jd, result_jd;
+
+    start_jd = timestamptz_to_jd(ts);
+    stop_jd  = start_jd + DEFAULT_WINDOW_DAYS;
+
+    result_jd = find_next_crossing(BTYPE_SUN, 0, obs,
+                                   start_jd, stop_jd,
+                                   NAUTICAL_TWILIGHT_RAD, true);
+
+    if (result_jd < 0.0)
+        PG_RETURN_NULL();
+
+    PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
+}
+
+
+/* ================================================================
+ * sun_nautical_dusk(observer, timestamptz) -> timestamptz
+ *
+ * Returns the next time nautical twilight ends (Sun crosses -12 deg
+ * heading downward).  After nautical dusk the horizon is no longer
+ * visible at sea; bright stars remain visible.
+ * ================================================================
+ */
+Datum
+sun_nautical_dusk(PG_FUNCTION_ARGS)
+{
+    pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
+    int64        ts  = PG_GETARG_INT64(1);
+    double       start_jd, stop_jd, result_jd;
+
+    start_jd = timestamptz_to_jd(ts);
+    stop_jd  = start_jd + DEFAULT_WINDOW_DAYS;
+
+    result_jd = find_next_crossing(BTYPE_SUN, 0, obs,
+                                   start_jd, stop_jd,
+                                   NAUTICAL_TWILIGHT_RAD, false);
+
+    if (result_jd < 0.0)
+        PG_RETURN_NULL();
+
+    PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
+}
+
+
+/* ================================================================
+ * sun_astronomical_dawn(observer, timestamptz) -> timestamptz
+ *
+ * Returns the next time astronomical twilight begins (Sun crosses
+ * -18 deg heading upward).  Before astronomical dawn the sky is
+ * fully dark — faintest objects are observable.
+ * ================================================================
+ */
+Datum
+sun_astronomical_dawn(PG_FUNCTION_ARGS)
+{
+    pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
+    int64        ts  = PG_GETARG_INT64(1);
+    double       start_jd, stop_jd, result_jd;
+
+    start_jd = timestamptz_to_jd(ts);
+    stop_jd  = start_jd + DEFAULT_WINDOW_DAYS;
+
+    result_jd = find_next_crossing(BTYPE_SUN, 0, obs,
+                                   start_jd, stop_jd,
+                                   ASTRONOMICAL_TWILIGHT_RAD, true);
+
+    if (result_jd < 0.0)
+        PG_RETURN_NULL();
+
+    PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
+}
+
+
+/* ================================================================
+ * sun_astronomical_dusk(observer, timestamptz) -> timestamptz
+ *
+ * Returns the next time astronomical twilight ends (Sun crosses
+ * -18 deg heading downward).  After astronomical dusk the sky is
+ * fully dark — faintest objects become observable.
+ * ================================================================
+ */
+Datum
+sun_astronomical_dusk(PG_FUNCTION_ARGS)
+{
+    pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
+    int64        ts  = PG_GETARG_INT64(1);
+    double       start_jd, stop_jd, result_jd;
+
+    start_jd = timestamptz_to_jd(ts);
+    stop_jd  = start_jd + DEFAULT_WINDOW_DAYS;
+
+    result_jd = find_next_crossing(BTYPE_SUN, 0, obs,
+                                   start_jd, stop_jd,
+                                   ASTRONOMICAL_TWILIGHT_RAD, false);
+
+    if (result_jd < 0.0)
+        PG_RETURN_NULL();
+
+    PG_RETURN_TIMESTAMPTZ(jd_to_timestamptz(result_jd));
+}
--- a/test/expected/v016_features.out
+++ b/test/expected/v016_features.out
@ -0,0 +1,243 @@
+-- v016_features.sql -- Tests for v0.16.0: twilight, lunar phase, planet magnitude
+--
+-- Verifies twilight dawn/dusk, lunar phase calculations,
+-- and planet apparent magnitude.
+CREATE EXTENSION IF NOT EXISTS pg_orrery;
+NOTICE:  extension "pg_orrery" already exists, skipping
+-- ============================================================
+-- Twilight: ordering (astronomical < nautical < civil < sunrise)
+-- Eagle, Idaho on the 2024 summer solstice
+-- ============================================================
+-- Dawn ordering: astronomical dawn < nautical dawn < civil dawn < sunrise
+-- Use midnight MDT (07:00 UTC) so all "next" events land on the same morning
+SELECT sun_astronomical_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+     < sun_nautical_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+    AS astro_before_nautical;
+ astro_before_nautical 
+-----------------------
+ t
+(1 row)
+
+SELECT sun_nautical_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+     < sun_civil_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+    AS nautical_before_civil;
+ nautical_before_civil 
+-----------------------
+ t
+(1 row)
+
+SELECT sun_civil_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+     < sun_next_rise('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+    AS civil_before_sunrise;
+ civil_before_sunrise 
+----------------------
+ t
+(1 row)
+
+-- Dusk ordering: sunset < civil dusk < nautical dusk < astronomical dusk
+-- Noon MDT (18:00 UTC) ensures all dusk events are still ahead
+SELECT sun_next_set('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+     < sun_civil_dusk('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+    AS sunset_before_civil;
+ sunset_before_civil 
+---------------------
+ t
+(1 row)
+
+SELECT sun_civil_dusk('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+     < sun_nautical_dusk('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+    AS civil_before_nautical;
+ civil_before_nautical 
+-----------------------
+ t
+(1 row)
+
+SELECT sun_nautical_dusk('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+     < sun_astronomical_dusk('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+    AS nautical_before_astro;
+ nautical_before_astro 
+-----------------------
+ t
+(1 row)
+
+-- ============================================================
+-- Twilight: civil dawn ~30 min before sunrise at mid-latitude
+-- ============================================================
+SELECT extract(epoch FROM
+    sun_next_rise('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+    - sun_civil_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+) BETWEEN 1200 AND 3600
+    AS civil_dawn_reasonable_offset;
+ civil_dawn_reasonable_offset 
+------------------------------
+ t
+(1 row)
+
+-- ============================================================
+-- Twilight: high latitude summer -- no astronomical darkness
+-- At 60N in June, astronomical dusk should be NULL (never gets dark enough)
+-- ============================================================
+SELECT sun_astronomical_dusk('(60.0,25.0,0)'::observer, '2024-06-21 00:00:00+00'::timestamptz) IS NULL
+    AS no_astro_dark_60n_summer;
+ no_astro_dark_60n_summer 
+--------------------------
+ t
+(1 row)
+
+-- ============================================================
+-- Lunar phase: known full moon (2024-01-25 ~17:54 UTC)
+-- Phase angle should be near 180 deg, illumination near 1.0
+-- ============================================================
+SELECT round(moon_phase_angle('2024-01-25 18:00:00+00'::timestamptz)::numeric, 0)
+    BETWEEN 170 AND 190
+    AS full_moon_angle_near_180;
+ full_moon_angle_near_180 
+--------------------------
+ t
+(1 row)
+
+SELECT round(moon_illumination('2024-01-25 18:00:00+00'::timestamptz)::numeric, 2)
+    >= 0.95
+    AS full_moon_high_illumination;
+ full_moon_high_illumination 
+-----------------------------
+ t
+(1 row)
+
+SELECT moon_phase_name('2024-01-25 18:00:00+00'::timestamptz) = 'full_moon'
+    AS full_moon_named;
+ full_moon_named 
+-----------------
+ t
+(1 row)
+
+-- ============================================================
+-- Lunar phase: known new moon (2024-01-11 ~11:57 UTC)
+-- Phase angle should be near 0 or 360, illumination near 0
+-- ============================================================
+SELECT moon_illumination('2024-01-11 12:00:00+00'::timestamptz)
+    < 0.05
+    AS new_moon_low_illumination;
+ new_moon_low_illumination 
+---------------------------
+ t
+(1 row)
+
+SELECT moon_phase_name('2024-01-11 12:00:00+00'::timestamptz) = 'new_moon'
+    AS new_moon_named;
+ new_moon_named 
+----------------
+ t
+(1 row)
+
+-- ============================================================
+-- Lunar phase: first quarter (2024-01-18 ~03:53 UTC)
+-- Phase angle near 90, illumination near 0.5
+-- ============================================================
+SELECT round(moon_phase_angle('2024-01-18 04:00:00+00'::timestamptz)::numeric, 0)
+    BETWEEN 80 AND 100
+    AS first_quarter_angle_near_90;
+ first_quarter_angle_near_90 
+-----------------------------
+ t
+(1 row)
+
+SELECT moon_illumination('2024-01-18 04:00:00+00'::timestamptz)
+    BETWEEN 0.4 AND 0.6
+    AS first_quarter_half_illuminated;
+ first_quarter_half_illuminated 
+--------------------------------
+ t
+(1 row)
+
+SELECT moon_phase_name('2024-01-18 04:00:00+00'::timestamptz) = 'first_quarter'
+    AS first_quarter_named;
+ first_quarter_named 
+---------------------
+ t
+(1 row)
+
+-- ============================================================
+-- Moon age: new moon has age near 0, full moon near 14.7
+-- ============================================================
+SELECT moon_age('2024-01-11 12:00:00+00'::timestamptz) < 2.0
+    AS new_moon_young;
+ new_moon_young 
+----------------
+ t
+(1 row)
+
+SELECT moon_age('2024-01-25 18:00:00+00'::timestamptz)
+    BETWEEN 12.0 AND 17.0
+    AS full_moon_age_midcycle;
+ full_moon_age_midcycle 
+------------------------
+ t
+(1 row)
+
+-- ============================================================
+-- Illumination range: always [0, 1]
+-- ============================================================
+SELECT moon_illumination('2024-06-01 00:00:00+00'::timestamptz) BETWEEN 0.0 AND 1.0
+   AND moon_illumination('2024-09-01 00:00:00+00'::timestamptz) BETWEEN 0.0 AND 1.0
+   AND moon_illumination('2024-12-01 00:00:00+00'::timestamptz) BETWEEN 0.0 AND 1.0
+    AS illumination_always_valid;
+ illumination_always_valid 
+---------------------------
+ t
+(1 row)
+
+-- ============================================================
+-- Planet magnitude: Jupiter should be bright (negative mag)
+-- ============================================================
+SELECT planet_magnitude(5, '2024-01-15 00:00:00+00'::timestamptz) < 0.0
+    AS jupiter_is_bright;
+ jupiter_is_bright 
+-------------------
+ t
+(1 row)
+
+-- ============================================================
+-- Planet magnitude: Venus is the brightest planet
+-- ============================================================
+SELECT planet_magnitude(2, '2024-06-01 12:00:00+00'::timestamptz)
+     < planet_magnitude(4, '2024-06-01 12:00:00+00'::timestamptz)
+    AS venus_brighter_than_mars;
+ venus_brighter_than_mars 
+--------------------------
+ t
+(1 row)
+
+-- ============================================================
+-- Planet magnitude: Neptune is faint (~+7-8)
+-- ============================================================
+SELECT planet_magnitude(8, '2024-01-15 00:00:00+00'::timestamptz) > 7.0
+    AS neptune_is_faint;
+ neptune_is_faint 
+------------------
+ t
+(1 row)
+
+-- ============================================================
+-- Planet magnitude: all planets return finite values
+-- ============================================================
+SELECT bool_and(
+    planet_magnitude(body_id, '2024-01-15 00:00:00+00'::timestamptz) IS NOT NULL
+    AND planet_magnitude(body_id, '2024-01-15 00:00:00+00'::timestamptz) > -30
+    AND planet_magnitude(body_id, '2024-01-15 00:00:00+00'::timestamptz) < 30
+) AS all_magnitudes_finite
+FROM (VALUES (1),(2),(4),(5),(6),(7),(8)) AS t(body_id);
+ all_magnitudes_finite 
+-----------------------
+ t
+(1 row)
+
+-- ============================================================
+-- Planet magnitude: error cases
+-- ============================================================
+DO $$ BEGIN PERFORM planet_magnitude(0, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'body_id=0(Sun): %', SQLERRM; END $$;
+NOTICE:  body_id=0(Sun): planet_magnitude: body_id 0 must be 1-8 (Mercury-Neptune)
+DO $$ BEGIN PERFORM planet_magnitude(3, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'body_id=3(Earth): %', SQLERRM; END $$;
+NOTICE:  body_id=3(Earth): cannot compute magnitude for Earth from Earth
+DO $$ BEGIN PERFORM planet_magnitude(9, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'body_id=9: %', SQLERRM; END $$;
+NOTICE:  body_id=9: planet_magnitude: body_id 9 must be 1-8 (Mercury-Neptune)
--- a/test/sql/v016_features.sql
+++ b/test/sql/v016_features.sql
@ -0,0 +1,161 @@
+-- v016_features.sql -- Tests for v0.16.0: twilight, lunar phase, planet magnitude
+--
+-- Verifies twilight dawn/dusk, lunar phase calculations,
+-- and planet apparent magnitude.
+CREATE EXTENSION IF NOT EXISTS pg_orrery;
+
+-- ============================================================
+-- Twilight: ordering (astronomical < nautical < civil < sunrise)
+-- Eagle, Idaho on the 2024 summer solstice
+-- ============================================================
+
+-- Dawn ordering: astronomical dawn < nautical dawn < civil dawn < sunrise
+-- Use midnight MDT (07:00 UTC) so all "next" events land on the same morning
+SELECT sun_astronomical_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+     < sun_nautical_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+    AS astro_before_nautical;
+
+SELECT sun_nautical_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+     < sun_civil_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+    AS nautical_before_civil;
+
+SELECT sun_civil_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+     < sun_next_rise('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+    AS civil_before_sunrise;
+
+-- Dusk ordering: sunset < civil dusk < nautical dusk < astronomical dusk
+-- Noon MDT (18:00 UTC) ensures all dusk events are still ahead
+SELECT sun_next_set('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+     < sun_civil_dusk('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+    AS sunset_before_civil;
+
+SELECT sun_civil_dusk('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+     < sun_nautical_dusk('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+    AS civil_before_nautical;
+
+SELECT sun_nautical_dusk('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+     < sun_astronomical_dusk('(43.7,-116.4,800)'::observer, '2024-06-21 18:00:00+00'::timestamptz)
+    AS nautical_before_astro;
+
+-- ============================================================
+-- Twilight: civil dawn ~30 min before sunrise at mid-latitude
+-- ============================================================
+
+SELECT extract(epoch FROM
+    sun_next_rise('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+    - sun_civil_dawn('(43.7,-116.4,800)'::observer, '2024-06-21 07:00:00+00'::timestamptz)
+) BETWEEN 1200 AND 3600
+    AS civil_dawn_reasonable_offset;
+
+-- ============================================================
+-- Twilight: high latitude summer -- no astronomical darkness
+-- At 60N in June, astronomical dusk should be NULL (never gets dark enough)
+-- ============================================================
+
+SELECT sun_astronomical_dusk('(60.0,25.0,0)'::observer, '2024-06-21 00:00:00+00'::timestamptz) IS NULL
+    AS no_astro_dark_60n_summer;
+
+-- ============================================================
+-- Lunar phase: known full moon (2024-01-25 ~17:54 UTC)
+-- Phase angle should be near 180 deg, illumination near 1.0
+-- ============================================================
+
+SELECT round(moon_phase_angle('2024-01-25 18:00:00+00'::timestamptz)::numeric, 0)
+    BETWEEN 170 AND 190
+    AS full_moon_angle_near_180;
+
+SELECT round(moon_illumination('2024-01-25 18:00:00+00'::timestamptz)::numeric, 2)
+    >= 0.95
+    AS full_moon_high_illumination;
+
+SELECT moon_phase_name('2024-01-25 18:00:00+00'::timestamptz) = 'full_moon'
+    AS full_moon_named;
+
+-- ============================================================
+-- Lunar phase: known new moon (2024-01-11 ~11:57 UTC)
+-- Phase angle should be near 0 or 360, illumination near 0
+-- ============================================================
+
+SELECT moon_illumination('2024-01-11 12:00:00+00'::timestamptz)
+    < 0.05
+    AS new_moon_low_illumination;
+
+SELECT moon_phase_name('2024-01-11 12:00:00+00'::timestamptz) = 'new_moon'
+    AS new_moon_named;
+
+-- ============================================================
+-- Lunar phase: first quarter (2024-01-18 ~03:53 UTC)
+-- Phase angle near 90, illumination near 0.5
+-- ============================================================
+
+SELECT round(moon_phase_angle('2024-01-18 04:00:00+00'::timestamptz)::numeric, 0)
+    BETWEEN 80 AND 100
+    AS first_quarter_angle_near_90;
+
+SELECT moon_illumination('2024-01-18 04:00:00+00'::timestamptz)
+    BETWEEN 0.4 AND 0.6
+    AS first_quarter_half_illuminated;
+
+SELECT moon_phase_name('2024-01-18 04:00:00+00'::timestamptz) = 'first_quarter'
+    AS first_quarter_named;
+
+-- ============================================================
+-- Moon age: new moon has age near 0, full moon near 14.7
+-- ============================================================
+
+SELECT moon_age('2024-01-11 12:00:00+00'::timestamptz) < 2.0
+    AS new_moon_young;
+
+SELECT moon_age('2024-01-25 18:00:00+00'::timestamptz)
+    BETWEEN 12.0 AND 17.0
+    AS full_moon_age_midcycle;
+
+-- ============================================================
+-- Illumination range: always [0, 1]
+-- ============================================================
+
+SELECT moon_illumination('2024-06-01 00:00:00+00'::timestamptz) BETWEEN 0.0 AND 1.0
+   AND moon_illumination('2024-09-01 00:00:00+00'::timestamptz) BETWEEN 0.0 AND 1.0
+   AND moon_illumination('2024-12-01 00:00:00+00'::timestamptz) BETWEEN 0.0 AND 1.0
+    AS illumination_always_valid;
+
+-- ============================================================
+-- Planet magnitude: Jupiter should be bright (negative mag)
+-- ============================================================
+
+SELECT planet_magnitude(5, '2024-01-15 00:00:00+00'::timestamptz) < 0.0
+    AS jupiter_is_bright;
+
+-- ============================================================
+-- Planet magnitude: Venus is the brightest planet
+-- ============================================================
+
+SELECT planet_magnitude(2, '2024-06-01 12:00:00+00'::timestamptz)
+     < planet_magnitude(4, '2024-06-01 12:00:00+00'::timestamptz)
+    AS venus_brighter_than_mars;
+
+-- ============================================================
+-- Planet magnitude: Neptune is faint (~+7-8)
+-- ============================================================
+
+SELECT planet_magnitude(8, '2024-01-15 00:00:00+00'::timestamptz) > 7.0
+    AS neptune_is_faint;
+
+-- ============================================================
+-- Planet magnitude: all planets return finite values
+-- ============================================================
+
+SELECT bool_and(
+    planet_magnitude(body_id, '2024-01-15 00:00:00+00'::timestamptz) IS NOT NULL
+    AND planet_magnitude(body_id, '2024-01-15 00:00:00+00'::timestamptz) > -30
+    AND planet_magnitude(body_id, '2024-01-15 00:00:00+00'::timestamptz) < 30
+) AS all_magnitudes_finite
+FROM (VALUES (1),(2),(4),(5),(6),(7),(8)) AS t(body_id);
+
+-- ============================================================
+-- Planet magnitude: error cases
+-- ============================================================
+
+DO $$ BEGIN PERFORM planet_magnitude(0, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'body_id=0(Sun): %', SQLERRM; END $$;
+DO $$ BEGIN PERFORM planet_magnitude(3, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'body_id=3(Earth): %', SQLERRM; END $$;
+DO $$ BEGIN PERFORM planet_magnitude(9, '2024-01-15 00:00:00+00'::timestamptz); EXCEPTION WHEN OTHERS THEN RAISE NOTICE 'body_id=9: %', SQLERRM; END $$;