From b6c5149cd77febf4169f2f94dc3549462dff09c7 Mon Sep 17 00:00:00 2001 From: Ryan Malloy Date: Wed, 18 Feb 2026 16:45:25 -0700 Subject: [PATCH] =?UTF-8?q?Add=20Cosmic=20Queries=20Cookbook=20guide=20?= =?UTF-8?q?=E2=80=94=209=20cross-domain=20SQL=20recipes?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit New guide combining multiple pg_orrery function families with PostgreSQL analytical functions: Kirkwood gap histograms, Kepler regression, asteroid family taxonomy, universal sky report, planetary alignment detection, ISS eclipse timing, PostGIS ground tracks, solar system dashboard view, and satellite shell census. --- docs/astro.config.mjs | 1 + .../content/docs/guides/cosmic-queries.mdx | 435 ++++++++++++++++++ 2 files changed, 436 insertions(+) create mode 100644 docs/src/content/docs/guides/cosmic-queries.mdx 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).