diff --git a/docs/astro.config.mjs b/docs/astro.config.mjs
index 256e409..e15802c 100644
--- a/docs/astro.config.mjs
+++ b/docs/astro.config.mjs
@@ -61,6 +61,7 @@ export default defineConfig({
items: [
{ label: "Tracking Satellites", slug: "guides/tracking-satellites" },
{ label: "Observing the Solar System", slug: "guides/observing-solar-system" },
+ { label: "Cosmic Queries Cookbook", slug: "guides/cosmic-queries" },
{ label: "Planetary Moon Tracking", slug: "guides/planetary-moons" },
{ label: "Star Catalogs in SQL", slug: "guides/star-catalogs" },
{ label: "Comet & Asteroid Tracking", slug: "guides/comets-asteroids" },
diff --git a/docs/src/content/docs/guides/cosmic-queries.mdx b/docs/src/content/docs/guides/cosmic-queries.mdx
new file mode 100644
index 0000000..c26117c
--- /dev/null
+++ b/docs/src/content/docs/guides/cosmic-queries.mdx
@@ -0,0 +1,435 @@
+---
+title: Cosmic Queries Cookbook
+sidebar:
+ order: 3
+---
+
+import { Steps, Aside, Tabs, TabItem } from "@astrojs/starlight/components";
+
+Each pg_orrery guide covers a single domain — satellites, planets, comets, radio bursts. This page is different. These nine queries combine multiple pg_orrery function families with PostgreSQL's analytical engine to ask questions that span physical theories, orbital regimes, and even external extensions like PostGIS. They range from statistical analysis of 600,000+ asteroids to real-time cross-domain sky surveys.
+
+Every query here is copy-paste ready. Swap the observer coordinates for your location and the timestamps for your session.
+
+## Prerequisites
+
+Not every query needs the same data. Here's what to load before you start:
+
+| Data | Queries that use it | Setup |
+|------|---------------------|-------|
+| `asteroids` table (MPC catalog) | 1, 2, 3, 4 | See [Comet & Asteroid Tracking](/guides/comets-asteroids) — load the MPC export with `oe_from_mpc()` |
+| `satellites` table (TLE catalog) | 4, 6, 7, 9 | See [Building TLE Catalogs](/guides/catalog-management) — any catalog with a `tle` column works |
+| `countries` table (Natural Earth) | 7 | PostGIS + Natural Earth boundaries — [setup below](#postgis-setup) |
+| PostGIS extension | 7 | `CREATE EXTENSION IF NOT EXISTS postgis;` |
+| None — built-in functions only | 5, 8 | Just pg_orrery |
+
+The expected table schemas:
+
+```sql
+-- Asteroids: name + orbital elements from MPC
+CREATE TABLE asteroids (
+ name text PRIMARY KEY,
+ oe orbital_elements
+);
+
+-- Satellites: NORAD ID + parsed TLE
+CREATE TABLE satellites (
+ norad_id integer PRIMARY KEY,
+ name text,
+ tle tle
+);
+```
+
+---
+
+## The Asteroid Belt as a Dataset
+
+The MPC catalog isn't just a list of orbits — it's a dataset with 600,000+ rows and rich statistical structure. PostgreSQL's aggregate functions turn it into an orbital mechanics laboratory.
+
+### 1. Kirkwood Gaps — Jupiter's Gravitational Fingerprint
+
+In 1866, Daniel Kirkwood noticed that asteroids avoid certain orbital distances. The gaps correspond to mean-motion resonances with Jupiter: an asteroid at 2.50 AU completes exactly 3 orbits for every 1 of Jupiter's (the 3:1 resonance). Over millions of years, Jupiter's periodic gravitational nudges clear these orbits out.
+
+`width_bucket()` bins the semi-major axes into a 200-bin histogram across the main belt. The depletions at 2.50, 2.82, 2.96, and 3.28 AU are unmistakable:
+
+```sql
+WITH belt AS (
+ SELECT oe_semi_major_axis(oe) AS a_au
+ FROM asteroids
+ WHERE oe_eccentricity(oe) < 0.4
+ AND oe_semi_major_axis(oe) IS NOT NULL
+ AND oe_semi_major_axis(oe) BETWEEN 1.5 AND 5.5
+)
+SELECT round((1.5 + (bucket - 1) * 0.02)::numeric, 2) AS a_au,
+ count AS n_asteroids
+FROM (
+ SELECT width_bucket(a_au, 1.5, 5.5, 200) AS bucket, count(*) AS count
+ FROM belt GROUP BY bucket
+) sub
+ORDER BY a_au;
+```
+
+The output is a classic Kirkwood gap diagram. Plot `a_au` vs `n_asteroids` and the resonance depletions jump out — the 3:1 gap at 2.50 AU is the deepest, with the 5:2 (2.82 AU), 7:3 (2.96 AU), and 2:1 (3.28 AU) gaps clearly visible.
+
+### 2. Kepler's Third Law as a Regression
+
+Kepler published his third law in 1619: the square of a planet's orbital period is proportional to the cube of its semi-major axis, or equivalently $\log P = 1.5 \cdot \log a$. With `regr_slope()` and `regr_r2()`, you can verify this 400-year-old relationship against every bound asteroid in the MPC catalog:
+
+```sql
+WITH bounded AS (
+ SELECT oe_semi_major_axis(oe) AS a_au, oe_period_years(oe) AS p_yr
+ FROM asteroids
+ WHERE oe_semi_major_axis(oe) IS NOT NULL
+ AND oe_semi_major_axis(oe) BETWEEN 0.5 AND 100.0
+)
+SELECT round(regr_slope(ln(p_yr), ln(a_au))::numeric, 6) AS slope,
+ round(exp(regr_intercept(ln(p_yr), ln(a_au)))::numeric, 6) AS intercept_years,
+ regr_count(ln(p_yr), ln(a_au)) AS n_objects,
+ round(regr_r2(ln(p_yr), ln(a_au))::numeric, 12) AS r_squared
+FROM bounded;
+```
+
+The slope will be exactly 1.500000 (Kepler's 3/2 power law). The intercept will be 1.000000 years (because for $a = 1$ AU, $P = 1$ year — Earth). The $R^2$ will be 1.000000000000. Not approximately. Exactly. This isn't a statistical correlation — it's a mathematical identity baked into `oe_period_years()`, which computes $a^{3/2}$. The query is a 600,000-row proof that the accessor functions are self-consistent.
+
+
+
+### 3. Asteroid Family Taxonomy
+
+Collisional families — groups of asteroids created by a single catastrophic impact — cluster tightly in (semi-major axis, eccentricity) space. A 2D `width_bucket()` grid reveals these density peaks as hot spots:
+
+```sql
+WITH belt AS (
+ SELECT oe_semi_major_axis(oe) AS a,
+ oe_eccentricity(oe) AS e
+ FROM asteroids
+ WHERE oe_eccentricity(oe) < 0.4
+ AND oe_semi_major_axis(oe) IS NOT NULL
+ AND oe_semi_major_axis(oe) BETWEEN 2.0 AND 3.5
+)
+SELECT round((2.0 + (a_bin - 1) * 0.03)::numeric, 2) AS a_au,
+ round((0.0 + (e_bin - 1) * 0.01)::numeric, 2) AS ecc,
+ count(*) AS n
+FROM (
+ SELECT width_bucket(a, 2.0, 3.5, 50) AS a_bin,
+ width_bucket(e, 0.0, 0.4, 40) AS e_bin
+ FROM belt
+) sub
+GROUP BY a_bin, e_bin
+HAVING count(*) >= 10
+ORDER BY n DESC;
+```
+
+The highest-density cells correspond to known collisional families: Flora (~2.2 AU, e~0.15), Themis (~3.13 AU, e~0.15), Koronis (~2.87 AU, e~0.05), and Eos (~3.01 AU, e~0.07). The `HAVING count(*) >= 10` filter suppresses noise in sparsely populated cells. Increase the threshold to isolate only the major families; decrease it to reveal smaller groupings.
+
+---
+
+## Cross-Domain Observation
+
+These queries combine satellite tracking, planetary ephemerides, and solar observation — functions backed by different physical theories, unified through pg_orrery's common `topocentric` return type.
+
+### 4. Universal Sky Report — Everything at Once
+
+Four gravitational theories in one query. `planet_observe()` uses VSOP87, `moon_observe()` uses ELP2000-82B, `observe_safe()` uses SGP4/SDP4, and `small_body_observe()` uses two-body Keplerian propagation. They all return `topocentric`, so `UNION ALL` works:
+
+```sql
+WITH obs AS (SELECT '40.0N 105.3W 1655m'::observer AS o),
+sky AS (
+ -- Planets (VSOP87)
+ SELECT 'Mercury' AS body, planet_observe(1, o, now()) AS topo FROM obs
+ UNION ALL SELECT 'Venus', planet_observe(2, o, now()) FROM obs
+ UNION ALL SELECT 'Mars', planet_observe(4, o, now()) FROM obs
+ UNION ALL SELECT 'Jupiter', planet_observe(5, o, now()) FROM obs
+ UNION ALL SELECT 'Saturn', planet_observe(6, o, now()) FROM obs
+ UNION ALL SELECT 'Uranus', planet_observe(7, o, now()) FROM obs
+ UNION ALL SELECT 'Neptune', planet_observe(8, o, now()) FROM obs
+ -- Sun and Moon
+ UNION ALL SELECT 'Sun', sun_observe(o, now()) FROM obs
+ UNION ALL SELECT 'Moon', moon_observe(o, now()) FROM obs
+ -- Satellites (SGP4/SDP4) — observe_safe returns NULL for decayed TLEs
+ UNION ALL
+ SELECT s.name, observe_safe(s.tle, obs.o, now())
+ FROM satellites s, obs
+ WHERE s.norad_id IN (25544, 20580, 48274) -- ISS, HST, Tiangong
+ -- Asteroids (two-body Keplerian)
+ UNION ALL
+ SELECT a.name, small_body_observe(a.oe, obs.o, now())
+ FROM asteroids a, obs
+ WHERE a.name IN ('Ceres', 'Vesta', 'Pallas')
+)
+SELECT body,
+ round(topo_azimuth(topo)::numeric, 1) AS az,
+ round(topo_elevation(topo)::numeric, 1) AS el,
+ CASE WHEN topo_elevation(topo) > 0 THEN 'visible' ELSE 'below horizon' END AS status
+FROM sky
+WHERE topo IS NOT NULL
+ORDER BY topo_elevation(topo) DESC;
+```
+
+Replace the NORAD IDs and asteroid names with whatever interests you. The `observe_safe` call is important for satellites — a decayed TLE will return NULL instead of raising an error, and the `WHERE topo IS NOT NULL` filter drops it cleanly.
+
+### 5. Planetary Alignment Detector
+
+How close are any two planets in the sky right now? The angular separation between two objects at (az₁, el₁) and (az₂, el₂) is the spherical law of cosines. PostgreSQL's built-in `sind()`, `cosd()`, and `acosd()` work in degrees — matching the degree output of `topo_azimuth()` and `topo_elevation()`:
+
+```sql
+WITH obs AS (SELECT '40.0N 105.3W 1655m'::observer AS o),
+planets AS (
+ SELECT body_id, name,
+ planet_observe(body_id, o, now()) AS topo
+ FROM obs,
+ (VALUES (1,'Mercury'),(2,'Venus'),(4,'Mars'),
+ (5,'Jupiter'),(6,'Saturn')) AS p(body_id, name)
+)
+SELECT a.name AS body_a, b.name AS body_b,
+ round(acosd(
+ sind(topo_elevation(a.topo)) * sind(topo_elevation(b.topo)) +
+ cosd(topo_elevation(a.topo)) * cosd(topo_elevation(b.topo)) *
+ cosd(topo_azimuth(a.topo) - topo_azimuth(b.topo))
+ )::numeric, 1) AS separation_deg
+FROM planets a
+JOIN planets b ON a.body_id < b.body_id
+WHERE topo_elevation(a.topo) > 0
+ AND topo_elevation(b.topo) > 0
+ORDER BY separation_deg;
+```
+
+The `a.body_id < b.body_id` join condition gives each pair exactly once (Venus–Jupiter, not also Jupiter–Venus). Only above-horizon planets are included — no point measuring the angular separation of objects you can't see.
+
+To find the closest approach over a year, sweep with `generate_series` and pick the tightest dates:
+
+```sql
+WITH obs AS (SELECT '40.0N 105.3W 1655m'::observer AS o),
+sweep AS (
+ SELECT t,
+ planet_observe(5, o, t) AS jupiter,
+ planet_observe(6, o, t) AS saturn
+ FROM obs,
+ generate_series(
+ '2026-01-01'::timestamptz,
+ '2026-12-31'::timestamptz,
+ interval '1 day'
+ ) AS t
+)
+SELECT t::date AS date,
+ round(acosd(
+ sind(topo_elevation(jupiter)) * sind(topo_elevation(saturn)) +
+ cosd(topo_elevation(jupiter)) * cosd(topo_elevation(saturn)) *
+ cosd(topo_azimuth(jupiter) - topo_azimuth(saturn))
+ )::numeric, 1) AS separation_deg
+FROM sweep
+WHERE topo_elevation(jupiter) > 0
+ AND topo_elevation(saturn) > 0
+ORDER BY separation_deg
+LIMIT 10;
+```
+
+### 6. ISS Eclipse Timing — Shadow Entry and Exit
+
+A satellite enters Earth's shadow when the Sun is below the horizon at the satellite's nadir point. This chains three pg_orrery domains together: `subsatellite_point()` (SGP4 → geodetic), `observer_from_geodetic()` (geodetic → observer), and `sun_observe()` (VSOP87 → topocentric). The `lag()` window function then detects the sunlit/shadow transitions:
+
+```sql
+WITH orbit AS (
+ SELECT t,
+ subsatellite_point(s.tle, t) AS geo
+ FROM satellites s,
+ generate_series(now(), now() + interval '93 minutes', interval '30 seconds') AS t
+ WHERE s.norad_id = 25544
+),
+shadow AS (
+ SELECT t,
+ geodetic_lat(geo) AS lat,
+ geodetic_lon(geo) AS lon,
+ topo_elevation(
+ sun_observe(
+ observer_from_geodetic(geodetic_lat(geo), geodetic_lon(geo)),
+ t
+ )
+ ) AS sun_el_at_nadir
+ FROM orbit
+)
+SELECT t,
+ round(lat::numeric, 2) AS lat,
+ round(lon::numeric, 2) AS lon,
+ round(sun_el_at_nadir::numeric, 1) AS sun_el,
+ CASE WHEN sun_el_at_nadir < 0 THEN 'SHADOW' ELSE 'SUNLIT' END AS state,
+ CASE
+ WHEN sun_el_at_nadir < 0 AND lag(sun_el_at_nadir) OVER (ORDER BY t) >= 0
+ THEN '>>> ECLIPSE ENTRY'
+ WHEN sun_el_at_nadir >= 0 AND lag(sun_el_at_nadir) OVER (ORDER BY t) < 0
+ THEN '<<< ECLIPSE EXIT'
+ END AS transition
+FROM shadow
+ORDER BY t;
+```
+
+
+
+### 7. Ground Track Geography with PostGIS
+
+Where on Earth is the ISS flying over? Combine `ground_track()` with PostGIS spatial joins against Natural Earth country boundaries.
+
+The simpler approach: a point-in-polygon test at each time step. Each (lat, lon) from the ground track becomes a PostGIS point, joined against country polygons:
+
+```sql
+WITH track AS (
+ SELECT t, lat, lon, alt
+ FROM satellites s,
+ ground_track(s.tle, now(), now() + interval '93 minutes', interval '30 seconds')
+ WHERE s.norad_id = 25544
+)
+SELECT track.t,
+ round(track.lat::numeric, 2) AS lat,
+ round(track.lon::numeric, 2) AS lon,
+ round(track.alt::numeric, 0) AS alt_km,
+ c.name AS country
+FROM track
+LEFT JOIN countries c
+ ON ST_Contains(c.geom, ST_SetSRID(ST_MakePoint(track.lon, track.lat), 4326));
+```
+
+The `LEFT JOIN` keeps rows over oceans (where `country` is NULL). The `ST_MakePoint()` argument order is (longitude, latitude) — x before y, the PostGIS convention.
+
+
+
+For a communications footprint, buffer the subsatellite point by the satellite's horizon radius. At ISS altitude (~400 km), the radio horizon is approximately 2,300 km:
+
+```sql
+-- Countries within ISS radio line-of-sight right now
+WITH nadir AS (
+ SELECT subsatellite_point(s.tle, now()) AS geo
+ FROM satellites s
+ WHERE s.norad_id = 25544
+)
+SELECT c.name AS country
+FROM nadir, countries c
+WHERE ST_DWithin(
+ c.geom::geography,
+ ST_SetSRID(ST_MakePoint(geodetic_lon(geo), geodetic_lat(geo)), 4326)::geography,
+ 2300000 -- 2,300 km horizon radius in meters
+);
+```
+
+---
+
+## The Solar System as Data
+
+SQL views and aggregation functions turn pg_orrery's observation pipeline into data products — live dashboards and statistical breakdowns that update every time you query them.
+
+### 8. Celestial Clock — Solar System Dashboard
+
+One row. Every planet's elevation, the Moon, the Sun, Io's phase angle, Jupiter's central meridian longitude, and the decametric burst probability. Wrap it in a `CREATE VIEW` and `SELECT * FROM solar_system_now` becomes a real-time dashboard:
+
+```sql
+CREATE VIEW solar_system_now AS
+SELECT now() AS computed_at,
+ round(topo_elevation(sun_observe(o, now()))::numeric, 1) AS sun_el,
+ round(topo_elevation(moon_observe(o, now()))::numeric, 1) AS moon_el,
+ round(topo_elevation(planet_observe(1, o, now()))::numeric, 1) AS mercury_el,
+ round(topo_elevation(planet_observe(2, o, now()))::numeric, 1) AS venus_el,
+ round(topo_elevation(planet_observe(4, o, now()))::numeric, 1) AS mars_el,
+ round(topo_elevation(planet_observe(5, o, now()))::numeric, 1) AS jupiter_el,
+ round(topo_elevation(planet_observe(6, o, now()))::numeric, 1) AS saturn_el,
+ round(topo_elevation(planet_observe(7, o, now()))::numeric, 1) AS uranus_el,
+ round(topo_elevation(planet_observe(8, o, now()))::numeric, 1) AS neptune_el,
+ round(io_phase_angle(now())::numeric, 1) AS io_phase,
+ round(jupiter_cml(o, now())::numeric, 1) AS jupiter_cml,
+ round(jupiter_burst_probability(
+ io_phase_angle(now()), jupiter_cml(o, now()))::numeric, 2) AS burst_prob
+FROM (SELECT '40.0N 105.3W 1655m'::observer) AS cfg(o);
+```
+
+Because the view uses `now()`, every `SELECT` recomputes against the current time — no refresh needed. The conditional aggregation approach (one column per planet) avoids `tablefunc`/`crosstab` entirely. Change the observer literal to your coordinates.
+
+
+
+### 9. Satellite Shell Census
+
+How many satellites occupy each orbital shell? Compute altitude from the TLE's mean motion using Kepler's third law ($a = (\mu / n^2)^{1/3}$, altitude $= a - R_\oplus$), then classify into LEO/MEO/GEO/HEO:
+
+```sql
+WITH altitudes AS (
+ SELECT norad_id, name,
+ power(
+ 398600.8 / power(tle_mean_motion(tle) * 2 * pi() / 86400.0, 2),
+ 1.0 / 3.0
+ ) - 6378.135 AS alt_km
+ FROM satellites
+ WHERE tle_mean_motion(tle) > 0
+)
+SELECT
+ CASE
+ WHEN alt_km < 2000 THEN 'LEO'
+ WHEN alt_km < 35786 THEN 'MEO'
+ WHEN alt_km < 35800 THEN 'GEO'
+ ELSE 'HEO/Other'
+ END AS shell,
+ count(*) AS n_satellites,
+ round(100.0 * count(*) / sum(count(*)) OVER (), 1) AS pct,
+ round(min(alt_km)::numeric, 0) AS min_alt_km,
+ round(percentile_cont(0.5) WITHIN GROUP (ORDER BY alt_km)::numeric, 0) AS median_alt_km,
+ round(max(alt_km)::numeric, 0) AS max_alt_km
+FROM altitudes
+WHERE alt_km > 100 -- filter decayed objects
+GROUP BY
+ CASE
+ WHEN alt_km < 2000 THEN 'LEO'
+ WHEN alt_km < 35786 THEN 'MEO'
+ WHEN alt_km < 35800 THEN 'GEO'
+ ELSE 'HEO/Other'
+ END
+ORDER BY min_alt_km;
+```
+
+The 398600.8 is WGS-72 $\mu$ (km³/s²) and 6378.135 is WGS-72 $a_e$ (km) — the same constants SGP4 uses internally. The `percentile_cont(0.5)` gives the median altitude per shell, which is more informative than the mean when distributions are skewed (LEO has a long tail from Molniya-type parking orbits).
+
+For a finer-grained altitude histogram within LEO — revealing the Starlink, ISS, sun-synchronous, and Iridium clusters:
+
+```sql
+WITH altitudes AS (
+ SELECT power(
+ 398600.8 / power(tle_mean_motion(tle) * 2 * pi() / 86400.0, 2),
+ 1.0 / 3.0
+ ) - 6378.135 AS alt_km
+ FROM satellites
+ WHERE tle_mean_motion(tle) > 0
+)
+SELECT round((150 + (bucket - 1) * 10)::numeric, 0) AS alt_km,
+ count(*) AS n_satellites
+FROM (
+ SELECT width_bucket(alt_km, 150, 2050, 190) AS bucket
+ FROM altitudes
+ WHERE alt_km BETWEEN 150 AND 2050
+) sub
+GROUP BY bucket
+ORDER BY alt_km;
+```
+
+Plot `alt_km` vs `n_satellites` and you'll see pronounced peaks: a massive spike near 550 km (Starlink's operational shell), a cluster around 780 km (Iridium NEXT), concentrations at 500–600 km (sun-synchronous polar orbits), and smaller peaks near 400 km (crewed missions) and 1200 km (older constellations).