Cosmic Queries Cookbook

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 — load the MPC export with oe_from_mpc()
satellites table (TLE catalog)	4, 6, 7, 9	See Building TLE Catalogs — any catalog with a tle column works
countries table (Natural Earth)	7	PostGIS + Natural Earth boundaries — setup below
PostGIS extension	7	CREATE EXTENSION IF NOT EXISTS postgis;
None — built-in functions only	5, 8	Just pg_orrery

The expected table schemas:

-- 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:

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:

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.

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).

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:

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, e0.15), Themis (3.13 AU, e0.15), Koronis (2.87 AU, e0.05), and Eos (3.01 AU, e0.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:

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():

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:

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:

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;

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.

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:

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.

PostGIS setup

Download Natural Earth 110m countries and load them:

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:

CREATE EXTENSION IF NOT EXISTS postgis;
CREATE INDEX ON countries USING gist (geom);

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:

-- 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:

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.

Parameterized version

For multiple observers, replace the literal with a function parameter or a lookup table:

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.

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:

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:

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).