Add Cosmic Queries Cookbook guide — 9 cross-domain SQL recipes

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.
2026-02-18 16:45:25 -07:00 · 2026-02-18 16:45:25 -07:00 · b6c5149cd7
commit b6c5149cd7
parent 53733daeba
2 changed files with 436 additions and 0 deletions
--- 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" },
--- a/docs/src/content/docs/guides/cosmic-queries.mdx
+++ 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.
 <Aside type="tip" title="Why R² = 1 exactly">
 `oe_period_years()` is defined as `a^1.5` where `a = q/(1-e)`. The regression isn't discovering a physical law — it's confirming that the accessor functions implement Kepler's third law without floating-point drift across the entire catalog. If you ever see R² < 1.0, something is wrong with your data (likely a corrupted MPC record).
 </Aside>
 ### 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;
 ```
 <Aside type="caution" title="Approximation accuracy">
 This treats the satellite's nadir point as the shadow boundary, which is geometrically simplified — it ignores the satellite's altitude above the surface and Earth's atmospheric refraction. For the ISS at ~400 km altitude, the shadow entry/exit times are accurate to roughly 10–20 seconds. For precise eclipse predictions, you'd need a cylindrical or conical shadow model. But for observation planning — knowing *approximately* when the ISS goes dark — this is very usable.
 </Aside>
 ### 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.
 <Aside type="note" title="PostGIS setup" id="postgis-setup">
 Download [Natural Earth 110m countries](https://www.naturalearthdata.com/downloads/110m-cultural-vectors/) and load them:
 ```bash
 wget https://naciscdn.org/naturalearth/110m/cultural/ne_110m_admin_0_countries.zip
 unzip ne_110m_admin_0_countries.zip
 shp2pgsql -s 4326 ne_110m_admin_0_countries.shp countries | psql -d your_database
 ```
 This creates a `countries` table with `name` (text) and `geom` (geometry) columns. Add a spatial index for faster lookups:
 ```sql
 CREATE EXTENSION IF NOT EXISTS postgis;
 CREATE INDEX ON countries USING gist (geom);
 ```
 </Aside>
 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.
 <Aside type="tip" title="Parameterized version">
 For multiple observers, replace the literal with a function parameter or a lookup table:
 ```sql
 SELECT s.* FROM observers o,
       LATERAL (SELECT * FROM solar_system_at(o.loc, now())) s;
 ```
 That requires wrapping the view logic in a `CREATE FUNCTION`, but the pattern is the same.
 </Aside>
 ### 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).