Add OD documentation and workflow translation pages for v0.4.0-v0.6.0

New pages: - OD function reference (tle_from_eci, tle_from_topocentric, tle_from_angles, tle_fit_residuals) - OD guide (ECI, topocentric, angles-only, range rate, weights, multi-observer, covariance interpretation) - From find_orb to SQL (OD workflow comparison) - From Poliastro to SQL (Lambert/Kepler comparison) Updated pages: - Corrected stale "No orbit determination" claim - Updated function counts and test suite counts - Added v0.4.0-v0.6.0 upgrade paths - Added OD to capabilities table, theory-to-code mapping, constants/accuracy reference - Added OD examples to Skyfield comparison and SQL Advantage - Fixed stale version references across workflow pages
2026-02-17 18:43:13 -07:00 · 2026-02-17 18:43:13 -07:00 · 3ea4e2bf53
commit 3ea4e2bf53
parent adfb6949e1
14 changed files with 1645 additions and 9 deletions
--- a/docs/astro.config.mjs
+++ b/docs/astro.config.mjs
@ -68,6 +68,7 @@ export default defineConfig({
            { label: "Interplanetary Trajectories", slug: "guides/interplanetary-trajectories" },
            { label: "Conjunction Screening", slug: "guides/conjunction-screening" },
            { label: "JPL DE Ephemeris", slug: "guides/de-ephemeris" },
            { label: "Orbit Determination", slug: "guides/orbit-determination" },
          ],
        },
        {
@ -77,6 +78,8 @@ export default defineConfig({
            { label: "From JPL Horizons to SQL", slug: "workflow/from-jpl-horizons" },
            { label: "From GMAT to SQL", slug: "workflow/from-gmat" },
            { label: "From Radio Jupiter Pro to SQL", slug: "workflow/from-radio-jupiter-pro" },
            { label: "From find_orb to SQL", slug: "workflow/from-find-orb" },
            { label: "From Poliastro to SQL", slug: "workflow/from-poliastro" },
            { label: "The SQL Advantage", slug: "workflow/sql-advantage" },
          ],
        },
@ -91,6 +94,7 @@ export default defineConfig({
            { label: "Functions: Radio", slug: "reference/functions-radio" },
            { label: "Functions: Transfers", slug: "reference/functions-transfers" },
            { label: "Functions: DE Ephemeris", slug: "reference/functions-de" },
            { label: "Functions: Orbit Determination", slug: "reference/functions-od" },
            { label: "Operators & GiST Index", slug: "reference/operators-gist" },
            { label: "Body ID Reference", slug: "reference/body-ids" },
            { label: "Constants & Accuracy", slug: "reference/constants-accuracy" },
--- a/docs/src/content/docs/architecture/theory-to-code.mdx
+++ b/docs/src/content/docs/architecture/theory-to-code.mdx
@ -103,6 +103,21 @@ The 82B revision is the version implemented. It provides geocentric ecliptic coo
 The Izzo solver uses Householder iterations for fast convergence and handles both short-way and long-way transfers. pg_orrery uses the prograde (short-way) solution by default.
 ## Orbit determination
 | Theory | Source | What it computes | Code location |
 |--------|--------|------------------|---------------|
 | Differential correction | Vallado (2013) Ch. 10 | Equinoctial element refinement via SVD least-squares | `src/od_solver.c` |
 | Gibbs method | Vallado Algorithm 54 | Initial velocity from three position vectors | `src/od_iod.c` |
 | Herrick-Gibbs | Vallado Algorithm 55 | Short-arc initial velocity (closely-spaced obs) | `src/od_iod.c` |
 | Gauss method | Vallado Algorithm 52 | Initial orbit from three angles-only (RA/Dec) observations | `src/od_iod.c` |
 ### References
 - Vallado, D.A. (2013). *Fundamentals of Astrodynamics and Applications*, 4th ed. Microcosm Press.
 The OD solver uses equinoctial elements to avoid singularities at zero eccentricity and inclination. LAPACK's `dgelss_` provides the SVD solve, with `dpotrf_`/`dpotri_` for formal covariance estimation.
 ## Radio emission
 | Theory | Source | What it computes | Code location |
--- a/docs/src/content/docs/getting-started/installation.mdx
+++ b/docs/src/content/docs/getting-started/installation.mdx
@ -97,7 +97,7 @@ import { Tabs, TabItem, Steps, Aside } from "@astrojs/starlight/components";
 ## Running the test suite
-If building from source, the regression tests verify all 68 functions across 12 test suites:
+If building from source, the regression tests verify all functions across 14 test suites:
 ```bash
 make installcheck PG_CONFIG=/usr/bin/pg_config
@ -115,6 +115,15 @@ ALTER EXTENSION pg_orrery UPDATE TO '0.2.0';
 -- From v0.2.0 to v0.3.0 (DE ephemeris support)
 ALTER EXTENSION pg_orrery UPDATE TO '0.3.0';
 -- From v0.3.0 to v0.4.0 (orbit determination)
 ALTER EXTENSION pg_orrery UPDATE TO '0.4.0';
 -- From v0.4.0 to v0.5.0 (OD enhancements: multi-observer, Gibbs IOD, covariance)
 ALTER EXTENSION pg_orrery UPDATE TO '0.5.0';
 -- From v0.5.0 to v0.6.0 (range rate, weighted obs, Gauss angles-only IOD)
 ALTER EXTENSION pg_orrery UPDATE TO '0.6.0';
 ```
 Each migration adds new functions while preserving existing data and functions.
--- a/docs/src/content/docs/getting-started/quick-start.mdx
+++ b/docs/src/content/docs/getting-started/quick-start.mdx
@ -112,5 +112,6 @@ You've seen the five domains pg_orrery covers. For deeper dives:
 - **[Tracking Satellites](/guides/tracking-satellites/)** — batch observation, conjunction screening, pass prediction workflows
 - **[Observing the Solar System](/guides/observing-solar-system/)** — "what's up tonight?" queries, rise/set times, conjunctions
 - **[Orbit Determination](/guides/orbit-determination/)** — fit TLEs from ECI positions, ground station observations, or angles-only RA/Dec data
 - **[JPL DE Ephemeris](/guides/de-ephemeris/)** — opt-in sub-milliarcsecond accuracy using JPL DE440/441 files
 - **[The SQL Advantage](/workflow/sql-advantage/)** — why doing this in SQL changes what's possible
--- a/docs/src/content/docs/getting-started/what-is-pg-orrery.mdx
+++ b/docs/src/content/docs/getting-started/what-is-pg-orrery.mdx
@ -26,6 +26,7 @@ PostGIS added spatial awareness to PostgreSQL — suddenly your database underst
 | Jupiter radio | Carr et al. (1983) sources | `jupiter_burst_probability()` | Empirical probability |
 | Transfers | Lambert (Izzo, 2015) | `lambert_transfer()`, `lambert_c3()` | Ballistic two-body |
 | DE ephemeris (optional) | JPL DE440/441 | `planet_observe_de()`, `moon_observe_de()` | ~0.1 milliarcsecond |
 | Orbit determination | Differential correction (v0.4.0+) | `tle_from_eci()`, `tle_from_topocentric()`, `tle_from_angles()` | Depends on observation quality |
 ## Who it's for
@ -57,7 +58,7 @@ pg_orrery is a computation engine, not a complete application. Understanding wha
 **Not sub-arcsecond by default.** The built-in VSOP87 pipeline is accurate to about 1 arcsecond — sufficient for observation planning and visual astronomy. For precision work (dish pointing, occultation timing, astrometry), pg_orrery v0.3.0 supports [optional JPL DE440/441 ephemeris files](/guides/de-ephemeris/) that bring accuracy to ~0.1 milliarcsecond. DE is opt-in and requires a one-time GUC configuration.
-**Not a TLE source.** Bring your own TLEs from Space-Track, CelesTrak, or any other provider. pg_orrery parses and propagates them; it doesn't fetch them.
+**Not a TLE source.** Bring your own TLEs from Space-Track, CelesTrak, or any other provider. pg_orrery parses, propagates, and — since v0.4.0 — can [fit new TLEs from observations](/guides/orbit-determination/). But it doesn't fetch TLE catalogs.
 **Not a replacement for SPICE.** No BSP kernel support, no light-time iteration, no aberration corrections at the IAU 2000A level. With DE enabled, pg_orrery matches SPICE on raw planet position accuracy — the remaining gap is in apparent-position corrections (aberration, light-time, nutation) that matter for sub-arcsecond apparent coordinates.
--- a/docs/src/content/docs/guides/orbit-determination.mdx
+++ b/docs/src/content/docs/guides/orbit-determination.mdx
@ -0,0 +1,226 @@
 ---
 title: Orbit Determination
 sidebar:
  order: 10
 ---
 import { Steps, Aside, Tabs, TabItem } from "@astrojs/starlight/components";
 Orbit determination (OD) is the inverse of propagation: instead of computing where a satellite will be given its orbital elements, you compute the orbital elements given where the satellite was observed. pg_orrery's OD solver fits SGP4/SDP4 mean elements (as a TLE) from observations using differential correction with an SVD least-squares solve.
 Three observation types are supported:
 - **ECI** — position/velocity vectors in the TEME frame
 - **Topocentric** — azimuth, elevation, range (and optionally range rate) from a ground station
 - **Angles-only** — right ascension and declination (RA/Dec) from optical observations
 ## How differential correction works
 The solver starts with an initial orbit estimate (the "seed") and iteratively refines it to minimize the residuals between observed and computed positions. Each iteration:
 1. Propagates the current TLE to each observation time
 2. Computes residuals (observed minus computed)
 3. Builds the Jacobian matrix (partial derivatives of residuals with respect to orbital elements)
 4. Solves the normal equations via SVD to get a correction vector
 5. Applies the correction to the equinoctial elements and converts back to a TLE
 The solver uses equinoctial elements rather than classical Keplerian elements to avoid singularities at zero eccentricity and zero inclination — both common for real satellites.
 ## ECI fitting
 The simplest case: you have ECI position/velocity vectors (perhaps from another propagator, a simulation, or a precise ephemeris) and want to fit a TLE.
 ```sql
 -- Generate synthetic observations by propagating a known TLE
 WITH iss AS (
    SELECT '1 25544U 98067A   24001.50000000  .00016717  00000-0  10270-3 0  9025
 2 25544  51.6400 208.9163 0006703  30.1694  61.7520 15.50100486 00001'::tle AS tle
 ),
 obs AS (
    SELECT array_agg(sgp4_propagate(iss.tle, t) ORDER BY t) AS positions,
           array_agg(t ORDER BY t) AS times
    FROM iss,
         generate_series(
             '2024-01-01 12:00:00+00'::timestamptz,
             '2024-01-01 13:30:00+00'::timestamptz,
             interval '5 minutes') t
 )
 SELECT iterations,
       round(rms_final::numeric, 6) AS rms_km,
       status,
       round(condition_number::numeric, 1) AS cond
 FROM obs, tle_from_eci(obs.positions, obs.times);
 ```
 With clean synthetic data (propagated from the same TLE model), this converges in a few iterations with sub-meter RMS. With real observations containing measurement noise, expect RMS in the 0.1-10 km range depending on data quality.
 <Aside type="tip" title="No seed needed">
 `tle_from_eci()` can bootstrap its own seed using the Gibbs method (three position vectors to derive an initial velocity). You only need to provide a seed TLE when the Gibbs geometry is degenerate (nearly colinear positions) or when you want to start from a known approximate orbit.
 </Aside>
 ## Topocentric fitting
 When observations come from a ground station (radar tracks, optical telescopes with range data), use `tle_from_topocentric()`. The solver accounts for the observer's position on the rotating Earth when computing the observation geometry.
 ```sql
 -- Observe a satellite, then fit from the topocentric data
 WITH iss AS (
    SELECT '1 25544U 98067A   24001.50000000  .00016717  00000-0  10270-3 0  9025
 2 25544  51.6400 208.9163 0006703  30.1694  61.7520 15.50100486 00001'::tle AS tle
 ),
 obs AS (
    SELECT array_agg(
               observe(iss.tle, '40.0N 105.3W 1655m'::observer, t) ORDER BY t
           ) AS observations,
           array_agg(t ORDER BY t) AS times
    FROM iss,
         generate_series(
             '2024-01-01 12:00:00+00'::timestamptz,
             '2024-01-01 12:10:00+00'::timestamptz,
             interval '30 seconds') t
 )
 SELECT iterations,
       round(rms_final::numeric, 4) AS rms_km,
       status
 FROM iss, obs,
     tle_from_topocentric(
         obs.observations, obs.times,
         '40.0N 105.3W 1655m'::observer,
         seed := iss.tle
     );
 ```
 <Aside type="caution" title="Seed usually required">
 Unlike ECI fitting, topocentric fitting typically needs a seed TLE. The azimuth/elevation/range → orbit mapping is more nonlinear than the ECI case, and the solver may not converge from a poor initial guess.
 </Aside>
 ## Range rate fitting
 When your observations include Doppler or radar range-rate data, enable `fit_range_rate` to use it as an additional constraint:
 ```sql
 SELECT iterations, round(rms_final::numeric, 4) AS rms_km, status
 FROM tle_from_topocentric(
    observations := topo_obs,
    times := obs_times,
    obs := '40.0N 105.3W 1655m'::observer,
    seed := seed_tle,
    fit_range_rate := true
 );
 ```
 Range rate adds a 4th residual component per observation (alongside azimuth, elevation, and range). This is particularly valuable for short observation arcs where position-only data doesn't constrain the velocity well. The range-rate residuals are internally scaled by a factor of 10 (1 km/s maps to 10 km equivalent) to balance their magnitude against position residuals.
 ## Multi-observer fitting
 When observations come from multiple ground stations, use the multi-observer variant. Each observation is tagged with its originating station via the `observer_ids` array:
 ```sql
 -- Define two ground stations
 -- Station 1: Boulder, CO
 -- Station 2: Los Angeles, CA
 SELECT iterations, round(rms_final::numeric, 4) AS rms_km, status
 FROM tle_from_topocentric(
    observations := all_topo_obs,
    times := all_times,
    observers := ARRAY[
        '40.0N 105.3W 1655m',
        '34.1N 118.3W 100m'
    ]::observer[],
    observer_ids := station_ids,    -- e.g., ARRAY[1,1,1,2,2,2]
    seed := seed_tle
 );
 ```
 The `observer_ids` array uses 1-based indexing into `observers[]`. Observations from different stations can be interleaved in any order — the solver uses `observer_ids[i]` to look up the correct station geometry for each observation.
 Multi-observer data improves orbit determination geometry, especially for short arcs. Two well-separated stations observing the same pass can constrain an orbit that a single station cannot.
 ## Weighted observations
 When observations have different accuracies (e.g., radar range vs. optical angles, or a high-quality station vs. a noisy one), use the `weights` parameter:
 ```sql
 SELECT iterations, round(rms_final::numeric, 4) AS rms_km, status
 FROM tle_from_eci(
    positions := obs_positions,
    times := obs_times,
    weights := ARRAY[1.0, 1.0, 0.5, 1.0, 0.2, 1.0]
 );
 ```
 Higher weights give observations more influence on the solution. A weight of 0.5 means "half as trustworthy as weight 1.0." This is useful for:
 - Down-weighting observations near the horizon (higher atmospheric refraction)
 - Down-weighting observations from less accurate sensors
 - Implementing iterative reweighting after identifying outliers via `tle_fit_residuals()`
 ## Angles-only (Gauss method)
 Optical observations often provide only RA/Dec — no range information. The `tle_from_angles()` function handles this case, using the Gauss method for initial orbit determination when no seed TLE is available:
 ```sql
 SELECT iterations,
       round(rms_final::numeric, 6) AS rms_rad,
       status,
       round(condition_number::numeric, 1) AS cond
 FROM tle_from_angles(
    ra_hours := ARRAY[12.345, 12.567, 12.789, 13.012, 13.234, 13.456],
    dec_degrees := ARRAY[45.1, 44.8, 44.3, 43.6, 42.8, 41.9],
    times := ARRAY[
        '2024-01-01 12:00+00', '2024-01-01 12:01+00',
        '2024-01-01 12:02+00', '2024-01-01 12:03+00',
        '2024-01-01 12:04+00', '2024-01-01 12:05+00'
    ]::timestamptz[],
    obs := '40.0N 105.3W 1655m'::observer
 );
 ```
 <Aside type="note" title="RMS units">
 For angles-only fitting, `rms_final` and `rms_initial` are in **radians**, not km. A typical good fit produces RMS around 0.001 radians (~3.4 arcminutes).
 </Aside>
 The Gauss method requires at least 3 observations with sufficient angular separation and time spread. The method works by solving for the geocentric range at the middle observation time, then deriving a full state vector. This bootstrap orbit is then refined by the standard DC solver.
 <Aside type="caution" title="Gauss limitations">
 The Gauss method can fail or produce spurious results when:
 - Observations span too short an arc (insufficient angular motion)
 - The satellite is near zenith for all observations (poor parallax geometry)
 - The time gaps between observations are too large (the linear approximation breaks down)
 When Gauss fails, provide a seed TLE from a catalog or from a separate IOD method.
 </Aside>
 RA is in hours [0, 24) and declination in degrees [-90, 90], matching the output convention of `star_observe()`. This lets you use the same coordinate system for both stellar calibration and satellite observations.
 ## Interpreting results
 Every OD function returns the same 8-column record. Here's how to use each field:
 | Field | What to check |
 |-------|---------------|
 | `fitted_tle` | NULL means the solver failed completely. Check `status` for the reason. |
 | `iterations` | Should be well below `max_iter`. If it equals `max_iter`, the solver didn't converge — increase `max_iter` or improve the seed. |
 | `rms_final` | The fit quality. For ECI: sub-km is good. For topocentric: 1-10 km is typical. For angles-only: < 0.01 rad is good. |
 | `rms_initial` | Compare with `rms_final`. A large ratio (rms_initial / rms_final) means the solver made significant improvement. If they're similar, the seed was already good or the solver made little progress. |
 | `status` | `'converged'` is success. `'max_iterations'` means more iterations needed. Other strings describe specific failures. |
 | `condition_number` | Measures how well-conditioned the geometry is. Values below ~1e6 are well-conditioned. Above ~1e10 suggests the observations don't constrain all orbital elements. |
 | `covariance` | Formal uncertainty. Extract diagonal elements for per-parameter variance. Off-diagonal elements show correlations between parameters. |
 | `nstate` | 6 (equinoctial elements) or 7 (with B*). The covariance matrix is `nstate x nstate`. |
 ## Tips
 **Seed selection.** For ECI fitting, the Gibbs IOD bootstrap usually works. For topocentric and angles-only fitting, start with a TLE from a catalog (Space-Track, CelesTrak) for the target object. The seed doesn't need to be precise — it just needs to be in the right ballpark (correct orbit regime, approximate inclination).
 **Minimum observations.** The solver needs at least 6 observations for a 6-parameter fit (or 7 for B*). In practice, more is better — 10-20 well-distributed observations across at least one orbit give the solver enough leverage to constrain all elements.
 **Observation spacing.** Spread observations across the orbit. Clustered observations from a short arc constrain some elements well (position) but others poorly (period, eccentricity). A full orbit of data is ideal.
 **Convergence troubleshooting.** If `status` shows `'max_iterations'`:
 1. Try a better seed TLE (closer to the true orbit)
 2. Increase `max_iter` to 30 or 50
 3. Check if observations contain outliers using `tle_fit_residuals()`
 4. Verify that the observation timestamps are correct (off-by-one-hour timezone errors are common)
 **B* fitting.** Only enable `fit_bstar` when you have observations spanning multiple orbits. B* captures atmospheric drag, which manifests as a slow period decay — a short arc doesn't have enough signal to constrain it. With insufficient data, fitting B* degrades the solution by introducing an under-determined 7th parameter.
--- a/docs/src/content/docs/guides/tracking-satellites.mdx
+++ b/docs/src/content/docs/guides/tracking-satellites.mdx
@ -44,7 +44,7 @@ pg_orrery propagates TLEs and computes look angles. It does not replace the full
 - **No real-time GUI.** GPredict and STK provide map displays, polar plots, and Doppler displays. pg_orrery returns numbers. Use any visualization tool to render its output.
 - **No rotator control.** Hamlib drives antenna rotators. pg_orrery computes the azimuth and elevation values Hamlib would consume, but it has no hardware interface.
 - **No TLE fetching.** Bring your own TLEs from Space-Track, CelesTrak, or any provider. pg_orrery parses and propagates them.
- **No orbit determination.** pg_orrery propagates existing TLEs. It does not fit orbits from observations.
+- **Orbit determination available.** Since v0.4.0, pg_orrery can fit TLEs from ECI, topocentric, or angles-only observations via differential correction. See the [Orbit Determination guide](/guides/orbit-determination/).
 - **No high-precision propagation.** SGP4/SDP4 accuracy degrades with TLE age. For operational conjunction assessment, use SP ephemerides or owner/operator-provided state vectors. pg_orrery's GiST screening finds candidates; you verify with better data.
 ## Try it
--- a/docs/src/content/docs/index.mdx
+++ b/docs/src/content/docs/index.mdx
@ -43,6 +43,12 @@ import { Card, CardGrid, LinkCard } from "@astrojs/starlight/components";
    Pork chop plots as SQL CROSS JOINs — 22,500 transfer solutions in
    8.3 seconds. Departure C3, arrival C3, time of flight, transfer SMA.
  </Card>
  <Card title="Determine orbits" icon="seti:satellite">
    Fit TLEs from ECI ephemeris, ground station observations (az/el/range with
    optional range rate), or angles-only RA/Dec. Gauss and Gibbs methods for
    seed-free initial orbit determination. Weighted least-squares with
    covariance output.
  </Card>
 </CardGrid>
 ## Explore the docs
@ -55,12 +61,12 @@ import { Card, CardGrid, LinkCard } from "@astrojs/starlight/components";
  />
  <LinkCard
    title="Guides"
-    description="Domain-specific walkthroughs for satellites, planets, moons, stars, comets, radio, and trajectories"
+    description="Domain-specific walkthroughs for satellites, planets, moons, stars, comets, radio, trajectories, and orbit determination"
    href="/guides/tracking-satellites/"
  />
  <LinkCard
    title="Workflow Translation"
-    description="Side-by-side comparisons: how you do it today vs. how pg_orrery changes the game"
+    description="Side-by-side comparisons with Skyfield, Horizons, GMAT, find_orb, and Poliastro"
    href="/workflow/from-skyfield/"
  />
  <LinkCard
--- a/docs/src/content/docs/reference/constants-accuracy.mdx
+++ b/docs/src/content/docs/reference/constants-accuracy.mdx
@ -207,6 +207,18 @@ Keplerian propagation ignores gravitational perturbations from planets, non-grav
 For preliminary mission design and pork chop plot generation, these limitations are standard and expected.
 ### Orbit Determination
 | Quantity | Notes |
 |----------|-------|
 | ECI fitting RMS | Sub-km typical with clean observations (round-trip from propagated state) |
 | Topocentric fitting RMS | ~1-10 km depending on arc length and observation spacing |
 | Angles-only fitting RMS | Output in radians. ~0.001 rad with good geometry |
 | Condition number | Formal indicator of solution quality. Values above ~1e10 suggest poorly-constrained geometry |
 | Covariance | Formal (H^T H)^{-1} from final Jacobian. Optimistic; does not include systematic errors |
 **Limitations:** The DC solver fits SGP4/SDP4 mean elements (6 equinoctial + optional B*). Accuracy is bounded by the TLE/SGP4 accuracy floor (~1 km at epoch for LEO). Range rate uses a fixed scale factor (OD_RR_SCALE = 10.0, mapping 1 km/s to 10 km equivalent). Gauss IOD requires at least 3 well-spaced observations with sufficient angular separation.
 ---
 ## Reference Publications
@ -222,3 +234,4 @@ For preliminary mission design and pork chop plot generation, these limitations
 | MarsSat | Jacobson, R.A. "The orbits and masses of the Martian satellites and the libration of Phobos." The Astronomical Journal, 139, 668-679, 2010. |
 | Carr source regions | Carr, T.D., Desch, M.D., Alexander, J.K. "Phenomenology of magnetospheric radio emissions." In Physics of the Jovian Magnetosphere, Cambridge Univ. Press, 1983. |
 | Lambert solver | Battin, R.H. "An Introduction to the Mathematics and Methods of Astrodynamics." AIAA Education Series, Revised Edition, 1999. |
 | Orbit determination | Vallado, D.A. "Fundamentals of Astrodynamics and Applications." 4th ed., Microcosm Press, 2013. |
--- a/docs/src/content/docs/reference/functions-od.mdx
+++ b/docs/src/content/docs/reference/functions-od.mdx
@ -0,0 +1,434 @@
 ---
 title: "Functions: Orbit Determination"
 sidebar:
  order: 8
 ---
 import { Aside, Tabs, TabItem } from "@astrojs/starlight/components";
 Functions for fitting TLE orbital elements from observations via differential correction (DC). The solver uses equinoctial elements internally to avoid singularities at zero eccentricity and inclination, with LAPACK SVD (`dgelss_`) for the least-squares solve.
 Three observation types are supported: ECI position/velocity, topocentric (az/el/range with optional range rate), and angles-only (RA/Dec). Each has single-observer and multi-observer variants where applicable.
 All OD functions share the same 8-column output record.
 ---
 ## Shared Output Record
 Every OD function returns a `RECORD` with these fields:
 | Field | Type | Description |
 |-------|------|-------------|
 | `fitted_tle` | `tle` | The fitted TLE. NULL if the solver did not converge. |
 | `iterations` | `int4` | Number of DC iterations performed. |
 | `rms_final` | `float8` | RMS residual after final iteration. Units depend on the observation type (km for ECI/topocentric, radians for angles-only). |
 | `rms_initial` | `float8` | RMS residual before the first iteration. Compare with `rms_final` to assess improvement. |
 | `status` | `text` | Convergence status: `'converged'`, `'max_iterations'`, or an error description. |
 | `condition_number` | `float8` | Condition number of the normal equations matrix. Values above ~1e10 suggest poorly-constrained geometry. |
 | `covariance` | `float8[]` | Formal covariance matrix `(H^T H)^{-1}` from the final Jacobian, stored as a flat array in row-major order. Length is `nstate * nstate`. |
 | `nstate` | `int4` | Number of estimated parameters (6 for equinoctial elements, 7 if `fit_bstar` is true). |
 <Aside type="note" title="Covariance interpretation">
 The covariance matrix is formal — it reflects the linear approximation at the solution and does not account for systematic errors in the observation model. It is optimistic. Use it for relative comparisons between solutions, not as an absolute uncertainty bound.
 </Aside>
 ---
 ## tle_from_eci
 Fit a TLE from ECI position/velocity observations. This is the simplest OD mode — the observations are already in the SGP4 propagation frame (TEME), so no observer geometry is involved.
 ### Signature
 ```sql
 tle_from_eci(
    positions eci_position[],
    times timestamptz[],
    seed tle DEFAULT NULL,
    fit_bstar boolean DEFAULT false,
    max_iter int4 DEFAULT 15,
    weights float8[] DEFAULT NULL,
    OUT fitted_tle tle,
    OUT iterations int4,
    OUT rms_final float8,
    OUT rms_initial float8,
    OUT status text,
    OUT condition_number float8,
    OUT covariance float8[],
    OUT nstate int4
 ) RETURNS RECORD
 ```
 ### Parameters
 | Parameter | Type | Default | Description |
 |-----------|------|---------|-------------|
 | `positions` | `eci_position[]` | — | Array of ECI position/velocity observations. Requires >= 6 observations. |
 | `times` | `timestamptz[]` | — | Observation timestamps. Must be same length as `positions`. |
 | `seed` | `tle` | `NULL` | Initial TLE estimate. When NULL, the solver uses the Gibbs or Herrick-Gibbs method to bootstrap an initial orbit from three position vectors. |
 | `fit_bstar` | `boolean` | `false` | When true, estimates B* drag coefficient as a 7th parameter. Requires a well-distributed observation arc. |
 | `max_iter` | `int4` | `15` | Maximum differential correction iterations. |
 | `weights` | `float8[]` | `NULL` | Per-observation weights. NULL means uniform weighting. Length must equal length of `positions`. Higher weight = more influence on the solution. |
 <Aside type="tip" title="Seed-free fitting">
 When no seed TLE is provided, the solver uses the Gibbs method (or Herrick-Gibbs for closely-spaced observations) to derive an initial velocity estimate from three position vectors. This makes `tle_from_eci()` fully seed-free for most use cases.
 </Aside>
 ### Example
 ```sql
 -- Round-trip test: propagate a TLE, then fit it back
 WITH iss AS (
    SELECT '1 25544U 98067A   24001.50000000  .00016717  00000-0  10270-3 0  9025
 2 25544  51.6400 208.9163 0006703  30.1694  61.7520 15.50100486 00001'::tle AS tle
 ),
 obs AS (
    SELECT array_agg(sgp4_propagate(iss.tle, t) ORDER BY t) AS positions,
           array_agg(t ORDER BY t) AS times
    FROM iss,
         generate_series(
             '2024-01-01 12:00:00+00'::timestamptz,
             '2024-01-01 13:30:00+00'::timestamptz,
             interval '5 minutes') t
 )
 SELECT iterations,
       round(rms_final::numeric, 6) AS rms_km,
       round(rms_initial::numeric, 3) AS rms_init_km,
       status,
       round(condition_number::numeric, 1) AS cond
 FROM obs, tle_from_eci(obs.positions, obs.times);
 ```
 ---
 ## tle_from_topocentric (single observer)
 Fit a TLE from topocentric (azimuth/elevation/range) observations collected by a single ground station.
 ### Signature
 ```sql
 tle_from_topocentric(
    observations topocentric[],
    times timestamptz[],
    obs observer,
    seed tle DEFAULT NULL,
    fit_bstar boolean DEFAULT false,
    max_iter int4 DEFAULT 15,
    fit_range_rate boolean DEFAULT false,
    weights float8[] DEFAULT NULL,
    OUT fitted_tle tle,
    OUT iterations int4,
    OUT rms_final float8,
    OUT rms_initial float8,
    OUT status text,
    OUT condition_number float8,
    OUT covariance float8[],
    OUT nstate int4
 ) RETURNS RECORD
 ```
 ### Parameters
 | Parameter | Type | Default | Description |
 |-----------|------|---------|-------------|
 | `observations` | `topocentric[]` | — | Array of topocentric observations (az, el, range, range_rate). Requires >= 6 observations. |
 | `times` | `timestamptz[]` | — | Observation timestamps. |
 | `obs` | `observer` | — | Ground station location. |
 | `seed` | `tle` | `NULL` | Initial TLE estimate. Unlike ECI fitting, topocentric fitting typically needs a seed for convergence. |
 | `fit_bstar` | `boolean` | `false` | Estimate B* drag coefficient. |
 | `max_iter` | `int4` | `15` | Maximum iterations. |
 | `fit_range_rate` | `boolean` | `false` | When true, includes range rate as a 4th residual component per observation. Use when you have Doppler or radar range-rate data. |
 | `weights` | `float8[]` | `NULL` | Per-observation weights. NULL = uniform. |
 <Aside type="note" title="Range rate scaling">
 Range rate residuals are internally scaled by `OD_RR_SCALE = 10.0`, mapping 1 km/s of range-rate error to 10 km equivalent in the residual vector. This prevents the typically small range-rate values from being swamped by position residuals.
 </Aside>
 ### Example
 ```sql
 -- Observe a satellite, then fit back from topocentric data
 WITH iss AS (
    SELECT '1 25544U 98067A   24001.50000000  .00016717  00000-0  10270-3 0  9025
 2 25544  51.6400 208.9163 0006703  30.1694  61.7520 15.50100486 00001'::tle AS tle
 ),
 obs AS (
    SELECT array_agg(observe(iss.tle, '40.0N 105.3W 1655m'::observer, t) ORDER BY t)
               AS observations,
           array_agg(t ORDER BY t) AS times
    FROM iss,
         generate_series(
             '2024-01-01 12:00:00+00'::timestamptz,
             '2024-01-01 12:10:00+00'::timestamptz,
             interval '30 seconds') t
 )
 SELECT iterations,
       round(rms_final::numeric, 4) AS rms_km,
       status
 FROM obs,
     tle_from_topocentric(
         obs.observations, obs.times,
         '40.0N 105.3W 1655m'::observer,
         seed := iss.tle
     );
 ```
 ---
 ## tle_from_topocentric (multi-observer)
 Fit a TLE from topocentric observations collected by multiple ground stations. The `observer_ids` array maps each observation to its originating station.
 ### Signature
 ```sql
 tle_from_topocentric(
    observations topocentric[],
    times timestamptz[],
    observers observer[],
    observer_ids int4[],
    seed tle DEFAULT NULL,
    fit_bstar boolean DEFAULT false,
    max_iter int4 DEFAULT 15,
    fit_range_rate boolean DEFAULT false,
    weights float8[] DEFAULT NULL,
    OUT fitted_tle tle,
    OUT iterations int4,
    OUT rms_final float8,
    OUT rms_initial float8,
    OUT status text,
    OUT condition_number float8,
    OUT covariance float8[],
    OUT nstate int4
 ) RETURNS RECORD
 ```
 ### Parameters
 | Parameter | Type | Default | Description |
 |-----------|------|---------|-------------|
 | `observations` | `topocentric[]` | — | All observations from all stations, interleaved in time order. |
 | `times` | `timestamptz[]` | — | Observation timestamps. |
 | `observers` | `observer[]` | — | Array of ground station locations (1-indexed). |
 | `observer_ids` | `int4[]` | — | Per-observation index into `observers[]`. `observer_ids[i]` identifies which station produced `observations[i]`. |
 | `seed` | `tle` | `NULL` | Initial TLE estimate. |
 | `fit_bstar` | `boolean` | `false` | Estimate B* drag coefficient. |
 | `max_iter` | `int4` | `15` | Maximum iterations. |
 | `fit_range_rate` | `boolean` | `false` | Include range rate in residuals. |
 | `weights` | `float8[]` | `NULL` | Per-observation weights. Useful when stations have different measurement accuracies. |
 ### Example
 ```sql
 -- Two ground stations observe the same satellite
 SELECT iterations, round(rms_final::numeric, 4) AS rms_km, status
 FROM tle_from_topocentric(
    observations := ARRAY[obs1_t1, obs1_t2, obs2_t1, obs2_t2]::topocentric[],
    times := ARRAY['2024-01-01 12:00+00', '2024-01-01 12:05+00',
                    '2024-01-01 12:02+00', '2024-01-01 12:07+00']::timestamptz[],
    observers := ARRAY['40.0N 105.3W 1655m', '34.1N 118.3W 100m']::observer[],
    observer_ids := ARRAY[1, 1, 2, 2],
    seed := seed_tle
 );
 ```
 <Aside type="tip" title="Observer ID indexing">
 `observer_ids` uses 1-based indexing into the `observers` array. Station 1 is `observers[1]`, station 2 is `observers[2]`, etc. This matches PostgreSQL's native array indexing convention.
 </Aside>
 ---
 ## tle_from_angles (single observer)
 Fit a TLE from angles-only (RA/Dec) observations collected by a single ground station. When no seed TLE is provided, the Gauss method derives an initial orbit from three observations — making this function fully seed-free.
 ### Signature
 ```sql
 tle_from_angles(
    ra_hours float8[],
    dec_degrees float8[],
    times timestamptz[],
    obs observer,
    seed tle DEFAULT NULL,
    fit_bstar boolean DEFAULT false,
    max_iter int4 DEFAULT 15,
    weights float8[] DEFAULT NULL,
    OUT fitted_tle tle,
    OUT iterations int4,
    OUT rms_final float8,
    OUT rms_initial float8,
    OUT status text,
    OUT condition_number float8,
    OUT covariance float8[],
    OUT nstate int4
 ) RETURNS RECORD
 ```
 ### Parameters
 | Parameter | Type | Default | Description |
 |-----------|------|---------|-------------|
 | `ra_hours` | `float8[]` | — | Right ascension values in hours, range [0, 24). Matches the `star_observe()` convention. |
 | `dec_degrees` | `float8[]` | — | Declination values in degrees, range [-90, 90]. |
 | `times` | `timestamptz[]` | — | Observation timestamps. Requires >= 3 observations. |
 | `obs` | `observer` | — | Ground station location. |
 | `seed` | `tle` | `NULL` | Initial TLE estimate. When NULL, the Gauss method bootstraps an initial orbit from 3 observations. |
 | `fit_bstar` | `boolean` | `false` | Estimate B* drag coefficient. |
 | `max_iter` | `int4` | `15` | Maximum iterations. |
 | `weights` | `float8[]` | `NULL` | Per-observation weights. |
 <Aside type="caution" title="RMS units">
 For angles-only fitting, `rms_final` and `rms_initial` are in **radians**, not km. A typical good fit produces RMS values around 0.001 radians (~0.06 degrees).
 </Aside>
 <Aside type="note" title="RA/Dec conventions">
 Right ascension is in hours [0, 24), matching the `star_observe()` output convention. Declination is in degrees [-90, 90]. Internally, both are converted to radians for the residual computation.
 </Aside>
 ### Example
 ```sql
 -- Angles-only OD: RA/Dec observations of a satellite
 SELECT iterations,
       round(rms_final::numeric, 6) AS rms_rad,
       status,
       round(condition_number::numeric, 1) AS cond
 FROM tle_from_angles(
    ra_hours := ARRAY[12.345, 12.567, 12.789, 13.012, 13.234, 13.456],
    dec_degrees := ARRAY[45.1, 44.8, 44.3, 43.6, 42.8, 41.9],
    times := ARRAY[
        '2024-01-01 12:00+00', '2024-01-01 12:01+00',
        '2024-01-01 12:02+00', '2024-01-01 12:03+00',
        '2024-01-01 12:04+00', '2024-01-01 12:05+00'
    ]::timestamptz[],
    obs := '40.0N 105.3W 1655m'::observer
 );
 ```
 ---
 ## tle_from_angles (multi-observer)
 Fit a TLE from angles-only (RA/Dec) observations collected by multiple ground stations. Uses the same Gauss IOD bootstrap as the single-observer variant when no seed is provided.
 ### Signature
 ```sql
 tle_from_angles(
    ra_hours float8[],
    dec_degrees float8[],
    times timestamptz[],
    observers observer[],
    observer_ids int4[],
    seed tle DEFAULT NULL,
    fit_bstar boolean DEFAULT false,
    max_iter int4 DEFAULT 15,
    weights float8[] DEFAULT NULL,
    OUT fitted_tle tle,
    OUT iterations int4,
    OUT rms_final float8,
    OUT rms_initial float8,
    OUT status text,
    OUT condition_number float8,
    OUT covariance float8[],
    OUT nstate int4
 ) RETURNS RECORD
 ```
 ### Parameters
 | Parameter | Type | Default | Description |
 |-----------|------|---------|-------------|
 | `ra_hours` | `float8[]` | — | Right ascension values in hours [0, 24). |
 | `dec_degrees` | `float8[]` | — | Declination values in degrees [-90, 90]. |
 | `times` | `timestamptz[]` | — | Observation timestamps from all stations, interleaved in time order. |
 | `observers` | `observer[]` | — | Array of ground station locations (1-indexed). |
 | `observer_ids` | `int4[]` | — | Per-observation index into `observers[]`. |
 | `seed` | `tle` | `NULL` | Initial TLE estimate. NULL = Gauss IOD bootstrap. |
 | `fit_bstar` | `boolean` | `false` | Estimate B* drag coefficient. |
 | `max_iter` | `int4` | `15` | Maximum iterations. |
 | `weights` | `float8[]` | `NULL` | Per-observation weights. Useful when stations have different apertures or sky conditions. |
 ### Example
 ```sql
 -- Two optical stations observe the same satellite
 SELECT iterations, round(rms_final::numeric, 6) AS rms_rad, status
 FROM tle_from_angles(
    ra_hours := ARRAY[12.3, 12.5, 12.7, 12.4, 12.6, 12.8],
    dec_degrees := ARRAY[45.0, 44.5, 44.0, 44.8, 44.3, 43.7],
    times := ARRAY[
        '2024-01-01 12:00+00', '2024-01-01 12:02+00',
        '2024-01-01 12:04+00', '2024-01-01 12:01+00',
        '2024-01-01 12:03+00', '2024-01-01 12:05+00'
    ]::timestamptz[],
    observers := ARRAY['40.0N 105.3W 1655m', '34.1N 118.3W 100m']::observer[],
    observer_ids := ARRAY[1, 1, 1, 2, 2, 2]
 );
 ```
 ---
 ## tle_fit_residuals
 Compute per-observation position residuals between a fitted TLE and the original ECI observations. Returns one row per observation with the XYZ and total position error in km. Use this to identify outlier observations or assess spatial error distribution.
 ### Signature
 ```sql
 tle_fit_residuals(
    fitted tle,
    positions eci_position[],
    times timestamptz[]
 ) RETURNS TABLE (
    t timestamptz,
    dx_km float8,
    dy_km float8,
    dz_km float8,
    pos_err_km float8
 )
 ```
 ### Parameters
 | Parameter | Type | Description |
 |-----------|------|-------------|
 | `fitted` | `tle` | The fitted TLE (typically the `fitted_tle` output from `tle_from_eci()`). |
 | `positions` | `eci_position[]` | The original ECI observations used for fitting. |
 | `times` | `timestamptz[]` | The original observation timestamps. |
 ### Returns
 One row per observation:
 | Column | Type | Unit | Description |
 |--------|------|------|-------------|
 | `t` | `timestamptz` | — | Observation time |
 | `dx_km` | `float8` | km | X-axis residual (observed - computed) |
 | `dy_km` | `float8` | km | Y-axis residual |
 | `dz_km` | `float8` | km | Z-axis residual |
 | `pos_err_km` | `float8` | km | Total position error: sqrt(dx^2 + dy^2 + dz^2) |
 ### Example
 ```sql
 -- After fitting, inspect per-observation residuals
 WITH fit AS (
    SELECT fitted_tle, positions, times
    FROM obs, tle_from_eci(obs.positions, obs.times)
 )
 SELECT t,
       round(dx_km::numeric, 4) AS dx,
       round(dy_km::numeric, 4) AS dy,
       round(dz_km::numeric, 4) AS dz,
       round(pos_err_km::numeric, 4) AS total_err
 FROM fit, tle_fit_residuals(fit.fitted_tle, fit.positions, fit.times)
 ORDER BY pos_err_km DESC;
 ```
 <Aside type="tip" title="Outlier detection">
 Observations with residuals significantly larger than the RMS are potential outliers. Consider removing them and re-fitting, or using the `weights` parameter to down-weight them.
 </Aside>
--- a/docs/src/content/docs/workflow/from-find-orb.mdx
+++ b/docs/src/content/docs/workflow/from-find-orb.mdx
@ -0,0 +1,373 @@
 ---
 title: From find_orb to SQL
 sidebar:
  order: 5
 description: Comparing find_orb's orbit determination workflow with pg_orrery's SQL-based differential correction.
 ---
 import { Tabs, TabItem, Aside, Steps } from "@astrojs/starlight/components";
 [find_orb](https://github.com/Bill-Gray/find_orb) is Bill Gray's orbit determination software — the same Bill Gray whose [sat_code SGP4 library](https://github.com/Bill-Gray/sat_code) is vendored inside pg_orrery. find_orb takes astrometric observations (RA/Dec positions with timestamps) and fits orbital elements via differential correction. It handles asteroids, comets, and artificial satellites. Amateur astronomers, minor planet observers, and satellite trackers use it worldwide.
 Since v0.4.0, pg_orrery has its own differential correction solver. The domain is narrower — pg_orrery fits SGP4/SDP4 mean elements (TLEs) for Earth-orbiting satellites, not heliocentric orbits for asteroids — but within that domain, the SQL approach offers batch processing, data integration, and automation that a desktop application cannot match.
 <Aside type="note" title="Different orbit models">
 find_orb fits osculating Keplerian elements suitable for numerical integration. pg_orrery fits SGP4 mean elements (TLEs) designed for the SGP4/SDP4 analytical propagator. These are different representations. A TLE absorbs SGP4's internal model biases — it only produces accurate positions when propagated with SGP4. This is a design choice, not a limitation: TLEs are the operational standard for satellite tracking, and pg_orrery's fitted TLEs plug directly into `observe()`, `predict_passes()`, and every other satellite function.
 </Aside>
 ## Fitting an orbit from RA/Dec observations
 The core workflow: you have a series of sky positions (right ascension and declination) for an object, and you want to determine its orbit.
 <Tabs>
  <TabItem label="find_orb">
    find_orb reads MPC 80-column format observation files:
    ```
    ISS      C2024 01 01.50000 12 20 42.0 +45 06 00.0          15.0 V      XXX
    ISS      C2024 01 01.50069 12 34 01.2 +44 48 00.0          15.0 V      XXX
    ISS      C2024 01 01.50139 12 47 18.6 +44 18 00.0          15.0 V      XXX
    ISS      C2024 01 01.50208 13 00 43.2 +43 36 00.0          15.0 V      XXX
    ISS      C2024 01 01.50278 13 14 02.4 +42 48 00.0          15.0 V      XXX
    ISS      C2024 01 01.50347 13 27 21.6 +41 54 00.0          15.0 V      XXX
    ```
    Then either:
    **GUI:** Launch find_orb, open the observation file, select the object, click "Full Step" repeatedly until the residuals converge. Inspect the residual map visually. Copy the fitted elements from the output panel.
    **CLI:** Run `fo obs_file.txt` to process observations non-interactively. Parse the output file for fitted elements.
    For each object, this is a separate run. Processing 50 objects from a night of observing means 50 file preparations and 50 find_orb runs. The results end up in text files that you then parse and load into your database.
  </TabItem>
  <TabItem label="pg_orrery (SQL)">
    ```sql
    -- Angles-only OD: RA/Dec observations → fitted TLE
    SELECT iterations,
           round(rms_final::numeric, 6) AS rms_rad,
           status,
           round(condition_number::numeric, 1) AS cond,
           fitted_tle
    FROM tle_from_angles(
        ra_hours := ARRAY[12.345, 12.567, 12.789, 13.012, 13.234, 13.456],
        dec_degrees := ARRAY[45.1, 44.8, 44.3, 43.6, 42.8, 41.9],
        times := ARRAY[
            '2024-01-01 12:00+00', '2024-01-01 12:01+00',
            '2024-01-01 12:02+00', '2024-01-01 12:03+00',
            '2024-01-01 12:04+00', '2024-01-01 12:05+00'
        ]::timestamptz[],
        obs := '40.0N 105.3W 1655m'::observer
    );
    ```
    No file format to prepare. RA in hours, Dec in degrees — direct numerical input. The Gauss method bootstraps the initial orbit automatically when no seed TLE is provided.
    The fitted TLE is immediately usable with any pg_orrery satellite function:
    ```sql
    -- Fit and immediately predict passes
    WITH od AS (
        SELECT fitted_tle
        FROM tle_from_angles(
            ra_hours := ARRAY[12.345, 12.567, 12.789, 13.012, 13.234, 13.456],
            dec_degrees := ARRAY[45.1, 44.8, 44.3, 43.6, 42.8, 41.9],
            times := ARRAY[
                '2024-01-01 12:00+00', '2024-01-01 12:01+00',
                '2024-01-01 12:02+00', '2024-01-01 12:03+00',
                '2024-01-01 12:04+00', '2024-01-01 12:05+00'
            ]::timestamptz[],
            obs := '40.0N 105.3W 1655m'::observer
        )
    )
    SELECT pass_aos(p) AS rise,
           pass_max_el(p) AS max_el,
           pass_los(p) AS set
    FROM od,
         predict_passes(od.fitted_tle,
                        '40.0N 105.3W 1655m'::observer,
                        now(), now() + interval '24 hours', 10.0) p;
    ```
    Observation → orbit fit → pass prediction in a single query.
  </TabItem>
 </Tabs>
 ## Batch orbit determination
 This is where the architecture diverges most sharply. find_orb processes one object at a time. pg_orrery processes however many your query describes.
 <Tabs>
  <TabItem label="find_orb">
    ```bash
    # Process a night's observations — one object at a time
    for obj in $(ls obs_files/*.txt); do
        fo "$obj" > results/$(basename "$obj" .txt).out 2>&1
    done
    # Parse each output file for fitted elements
    # Build a pipeline to extract orbital elements, residuals, etc.
    # Load results into your database
    ```
    For 50 objects, this means 50 invocations, 50 output files, and a parsing pipeline. Convergence failures need manual attention — find_orb's differential corrector may not converge for all objects, and the failure modes differ per object.
  </TabItem>
  <TabItem label="pg_orrery (SQL)">
    ```sql
    -- Batch OD: fit orbits for all objects observed tonight
    -- Assumes an observation table with RA/Dec per object
    SELECT obj_id,
           (od).iterations,
           round((od).rms_final::numeric, 6) AS rms_rad,
           (od).status,
           (od).fitted_tle
    FROM (
        SELECT obj_id,
               tle_from_angles(
                   ra_hours := array_agg(ra_hours ORDER BY obs_time),
                   dec_degrees := array_agg(dec_deg ORDER BY obs_time),
                   times := array_agg(obs_time ORDER BY obs_time),
                   obs := '40.0N 105.3W 1655m'::observer
               ) AS od
        FROM tonight_observations
        WHERE obs_time > now() - interval '12 hours'
        GROUP BY obj_id
        HAVING count(*) >= 3  -- Need at least 3 for Gauss IOD
    ) sub
    ORDER BY (od).rms_final;
    ```
    Every object observed tonight, fitted in one query. Objects that fail to converge show `status != 'converged'` — they stay in the result set instead of halting the pipeline. Sort by RMS to see the best fits first.
  </TabItem>
 </Tabs>
 ## Multi-observer observations
 Observations from multiple stations improve orbit geometry. Both tools support this, but the workflow differs.
 <Tabs>
  <TabItem label="find_orb">
    In MPC format, each observation line includes a 3-character observatory code (column 78-80). find_orb looks up the observatory coordinates from a `stations.txt` file. To combine data from multiple observers:
    1. Ensure all observers have MPC observatory codes (or define custom ones)
    2. Concatenate observation files, maintaining the 80-column format
    3. Make sure `stations.txt` contains coordinates for all observatory codes
    4. Load the combined file into find_orb
    The station lookup is implicit — based on the observatory code column in the fixed-width format.
  </TabItem>
  <TabItem label="pg_orrery (SQL)">
    ```sql
    -- Multi-observer OD with explicit station geometry
    SELECT (od).iterations,
           round((od).rms_final::numeric, 6) AS rms_rad,
           (od).status,
           round((od).condition_number::numeric, 1) AS cond
    FROM tle_from_angles(
        ra_hours := array_agg(ra ORDER BY t),
        dec_degrees := array_agg(dec ORDER BY t),
        times := array_agg(t ORDER BY t),
        observers := ARRAY[
            '40.0N 105.3W 1655m',    -- Boulder
            '34.1N 118.3W 100m',     -- Los Angeles
            '32.0N 110.9W 730m'      -- Tucson
        ]::observer[],
        observer_ids := array_agg(station_id ORDER BY t)
    ) AS od
    FROM multi_station_observations
    WHERE target = 'UNKNOWN-2024A';
    ```
    Station geometry is explicit — no lookup file, no observatory code registry. Each observation carries its station ID directly. Adding a new station is adding a value to the `observers` array.
  </TabItem>
 </Tabs>
 ## Weighted observations and range rate
 v0.6.0 added per-observation weights and range rate fitting. These features are most useful when combining heterogeneous data — different sensors, different accuracies, or a mix of optical and radar observations.
 <Tabs>
  <TabItem label="find_orb">
    find_orb handles observation weighting internally. It assigns weights based on observatory statistics from the MPC's Observatory Performance page — observatories with historically tighter residuals get higher effective weight. You can override this by editing the observation weights in the GUI, but there's no direct numerical control in the input file format.
    Range rate (Doppler) observations are supported through a separate observation type in the MPC format. Combining angle and range-rate data requires properly formatted input with both observation types interleaved.
  </TabItem>
  <TabItem label="pg_orrery (SQL)">
    ```sql
    -- Weighted multi-sensor fitting with range rate
    -- Station 1: precision radar (high weight, has range rate)
    -- Station 2: optical telescope (lower weight, angles + range only)
    SELECT (od).iterations,
           round((od).rms_final::numeric, 4) AS rms_km,
           (od).status,
           round((od).condition_number::numeric, 1) AS cond,
           (od).nstate AS params_fitted
    FROM tle_from_topocentric(
        observations := all_topo_obs,
        times := all_times,
        observers := ARRAY[
            '40.0N 105.3W 1655m',   -- radar station
            '34.1N 118.3W 100m'     -- optical station
        ]::observer[],
        observer_ids := station_ids,
        seed := seed_tle,
        fit_range_rate := true,      -- use Doppler data from radar
        weights := ARRAY[
            1.0, 1.0, 1.0, 1.0,     -- radar obs: full weight
            0.3, 0.3, 0.3, 0.3      -- optical obs: lower weight
        ]
    ) AS od;
    ```
    The `weights` array gives explicit numerical control per observation. The `fit_range_rate` flag adds range rate as a 4th residual component — particularly valuable when short observation arcs don't constrain velocity well from position data alone.
    ```sql
    -- Extract formal uncertainty from the covariance matrix
    -- Diagonal elements are the per-parameter variances
    WITH fit AS (
        SELECT (od).covariance AS cov,
               (od).nstate AS n,
               (od).condition_number AS cond
        FROM tle_from_topocentric(/* ... */) AS od
    )
    SELECT round(cond::numeric, 1) AS condition_number,
           round(sqrt(cov[1])::numeric, 6) AS sigma_1,
           round(sqrt(cov[n + 2])::numeric, 6) AS sigma_2,
           round(sqrt(cov[2 * n + 3])::numeric, 6) AS sigma_3
    FROM fit;
    ```
    The covariance matrix is formal `(H^T H)^{-1}` — optimistic, but useful for comparing relative solution quality across different observation geometries or weighting schemes.
  </TabItem>
 </Tabs>
 ## Residual inspection
 After fitting, you want to know which observations are good and which are outliers.
 <Tabs>
  <TabItem label="find_orb">
    find_orb displays a residual map in the GUI — a scatter plot of RA and Dec residuals per observation. Outliers are visually obvious. You can click to exclude individual observations and re-fit.
    In CLI mode, residuals appear in the output file, one line per observation. You parse the file to identify outliers programmatically.
  </TabItem>
  <TabItem label="pg_orrery (SQL)">
    For ECI-based fitting, `tle_fit_residuals()` returns per-observation residuals:
    ```sql
    -- Inspect per-observation residuals after fitting
    WITH fit AS (
        SELECT fitted_tle
        FROM tle_from_eci(obs_positions, obs_times)
    )
    SELECT t,
           round(pos_err_km::numeric, 4) AS error_km,
           CASE WHEN pos_err_km > 5.0 THEN 'OUTLIER' ELSE '' END AS flag
    FROM fit, tle_fit_residuals(fit.fitted_tle, obs_positions, obs_times)
    ORDER BY pos_err_km DESC;
    ```
    For iterative outlier rejection:
    ```sql
    -- Fit, identify outliers, remove them, re-fit
    WITH first_fit AS (
        SELECT fitted_tle, positions, times
        FROM tle_from_eci(obs_positions, obs_times)
    ),
    good_obs AS (
        SELECT array_agg(positions[i] ORDER BY i) AS clean_pos,
               array_agg(times[i] ORDER BY i) AS clean_times
        FROM first_fit,
             generate_series(1, array_length(positions, 1)) AS i
        WHERE (tle_fit_residuals(fitted_tle, positions, times)).pos_err_km < 5.0
    )
    SELECT iterations,
           round(rms_final::numeric, 6) AS rms_km,
           status
    FROM good_obs, tle_from_eci(good_obs.clean_pos, good_obs.clean_times);
    ```
    No manual clicking. The rejection criterion is a `WHERE` clause — change the threshold, re-run.
  </TabItem>
 </Tabs>
 ## Closing the loop: fit → propagate → observe
 find_orb produces orbital elements. To then predict passes, compute rise/set times, or screen for conjunctions, you need a separate tool — GPredict, Skyfield, or a custom propagator. The orbit determination result and the operational tracking live in different systems.
 pg_orrery keeps everything in one place:
 ```sql
 -- Full pipeline: observations → fit → catalog → predict passes
 WITH fit AS (
    SELECT tle_from_angles(
        ra_hours := obs_ra,
        dec_degrees := obs_dec,
        times := obs_times,
        obs := '40.0N 105.3W 1655m'::observer
    ) AS od
    FROM tonight_observations
    WHERE target = 'UNKNOWN-2024A'
 )
 -- Insert the fitted TLE into the satellite catalog
 INSERT INTO satellites (norad_id, name, tle, source)
 SELECT 99999, 'UNKNOWN-2024A', (fit.od).fitted_tle, 'local_od'
 FROM fit
 WHERE (fit.od).status = 'converged';
 -- Now predict passes for the newly cataloged object
 SELECT pass_aos(p) AS rise,
       pass_max_el(p) AS max_el,
       pass_los(p) AS set
 FROM satellites s,
     predict_passes(s.tle,
                    '40.0N 105.3W 1655m'::observer,
                    now(), now() + interval '48 hours', 10.0) p
 WHERE s.name = 'UNKNOWN-2024A';
 ```
 The fitted TLE goes into the same catalog table as Space-Track TLEs. Every pg_orrery function — `observe()`, `predict_passes()`, the GiST conjunction index — works on it immediately. No export, no format conversion, no tool switching.
 ## Where find_orb wins
 <Aside type="note" title="find_orb is the deeper tool for its domain">
 find_orb's OD capabilities are far broader than pg_orrery's. The comparison here is specifically about the satellite tracking workflow.
 </Aside>
 **Heliocentric orbits.** find_orb fits orbits for asteroids, comets, and interplanetary objects — heliocentric Keplerian elements with perturbations. pg_orrery's OD solver fits geocentric SGP4 mean elements for Earth-orbiting satellites only.
 **Perturbation models.** find_orb can include planetary perturbations, solar radiation pressure, and non-gravitational forces in the orbit fit. pg_orrery's DC solver fits pure SGP4 mean elements with optional B* drag.
 **Interactive refinement.** find_orb's GUI lets you visually inspect residuals, toggle observations, try different force models, and watch the solution converge. pg_orrery's solver is a function call — you set parameters and read results.
 **MPC ecosystem.** find_orb reads MPC 80-column format natively and integrates with the Minor Planet Center's observation database. If you're reporting asteroid discoveries, find_orb speaks the language.
 **Multiple solutions.** find_orb can identify and present multiple orbit solutions when the data is ambiguous (e.g., short-arc asteroid observations with two equally valid orbits). pg_orrery's solver returns a single solution.
 **Mature solver.** find_orb's differential corrector has decades of refinement and handles edge cases (near-parabolic orbits, objects at the boundary of Earth's sphere of influence, multi-opposition linkage) that pg_orrery's newer solver does not.
 ## Where pg_orrery wins
 **Batch processing.** Fit orbits for every object observed tonight in a single query. find_orb processes objects one at a time.
 **Data integration.** Fitted TLEs land in the same database as your satellite catalog, observation logs, contact schedules, and frequency assignments. `JOIN` the results with anything.
 **Automation.** A SQL query is a complete, repeatable specification. Set it up in a cron job or a materialized view and the pipeline runs itself. find_orb requires either manual GUI interaction or scripted CLI runs with output parsing.
 **Closed-loop tracking.** The fitted TLE immediately works with `observe()`, `predict_passes()`, and the GiST conjunction index. No format conversion, no tool switching, no exporting elements and importing them elsewhere.
 **Multi-observer without MPC codes.** pg_orrery takes observer coordinates directly — no registry, no observatory code, no station lookup file. Define a station as a string and use it.
 **Weighted observations.** pg_orrery's `weights` parameter lets you explicitly down-weight noisy observations or prioritize high-quality data. find_orb handles weighting internally based on observatory statistics.
 ## Migrating gradually
 <Steps>
 1. **Use find_orb for heliocentric objects.** Asteroids, comets, and interplanetary objects belong in find_orb. pg_orrery's OD solver is designed for Earth-orbiting satellites.
 2. **Use pg_orrery for satellite OD.** If your observations are of artificial satellites and you want TLEs that work with SGP4, use `tle_from_angles()` or `tle_from_topocentric()`. The fitted TLE integrates directly with your existing pg_orrery workflow.
 3. **Store observations in PostgreSQL.** Whether you process them with find_orb or pg_orrery, keeping observations in a database table lets you re-process with different parameters, track observation quality over time, and correlate with metadata.
 4. **Automate the pipeline.** Observation ingestion → OD → catalog update → pass prediction can be a scheduled SQL procedure. Manual find_orb runs become the exception for difficult cases, not the default.
 </Steps>
--- a/docs/src/content/docs/workflow/from-poliastro.mdx
+++ b/docs/src/content/docs/workflow/from-poliastro.mdx
@ -0,0 +1,418 @@
 ---
 title: From Poliastro to SQL
 sidebar:
  order: 6
 description: Side-by-side comparison of Poliastro Python workflows and equivalent pg_orrery SQL queries for orbital mechanics.
 ---
 import { Tabs, TabItem, Aside, Steps } from "@astrojs/starlight/components";
 [Poliastro](https://docs.poliastro.space/) is a Python library for interactive astrodynamics and orbital mechanics. Built on Astropy, it provides Keplerian propagation, Lambert transfers, orbit plotting, and maneuver analysis with a clean Pythonic API. If you're doing orbital mechanics in Jupyter notebooks, you've probably used it.
 pg_orrery overlaps with Poliastro on two core computations — Lambert transfers and Keplerian propagation — and extends beyond it with SGP4 satellite tracking, orbit determination, and GiST-indexed catalog operations. This page compares the workflows where both tools operate, and identifies where each is the better fit.
 ## Lambert transfers
 Both tools implement the Izzo (2015) Lambert solver. The same algorithm, different execution environments.
 <Tabs>
  <TabItem label="Poliastro (Python)">
    ```python
    from astropy import units as u
    from astropy.time import Time
    from poliastro.bodies import Earth, Mars, Sun
    from poliastro.ephem import Ephem
    from poliastro.iod import izzo
    from poliastro.util import norm
    # Get planet positions from built-in ephemeris
    dep_time = Time("2028-10-01", scale="tdb")
    arr_time = Time("2029-06-15", scale="tdb")
    earth_ephem = Ephem.from_body(Earth, dep_time)
    mars_ephem = Ephem.from_body(Mars, arr_time)
    r_earth = earth_ephem.rv(dep_time)[0]   # position vector
    r_mars = mars_ephem.rv(arr_time)[0]     # position vector
    tof = arr_time - dep_time               # time of flight
    # Solve Lambert problem
    (v_dep, v_arr), = izzo.lambert(Sun.k, r_earth, r_mars, tof)
    # Compute C3
    v_earth = earth_ephem.rv(dep_time)[1]
    v_inf = v_dep - v_earth
    c3 = norm(v_inf).to(u.km / u.s) ** 2
    print(f"Departure C3: {c3:.2f}")
    ```
    Several objects are involved: `Ephem` for planet positions, `Time` for epochs, Astropy units for dimensional analysis, and the `izzo` module for the actual solve. Each intermediate result has attached units. This is excellent for interactive exploration in a notebook — you can inspect each step, change units, and see what you're working with.
    For a single transfer, this is perfectly fine. The friction appears when you need many transfers.
  </TabItem>
  <TabItem label="pg_orrery (SQL)">
    ```sql
    SELECT round(c3_departure::numeric, 2) AS c3_km2s2,
           round(c3_arrival::numeric, 2) AS c3_arrive,
           round(v_inf_departure::numeric, 3) AS v_inf_dep_kms,
           round(tof_days::numeric, 1) AS flight_days,
           round(transfer_sma::numeric, 4) AS sma_au
    FROM lambert_transfer(3, 4,
           '2028-10-01'::timestamptz,
           '2029-06-15'::timestamptz);
    ```
    One function call. Planet positions, Lambert solve, and C3 computation happen internally. The body IDs (3=Earth, 4=Mars) and timestamps are the only inputs. No unit objects, no intermediate vectors, no ephemeris loading.
  </TabItem>
 </Tabs>
 ## Pork chop plots
 The parameter sweep that reveals launch windows. This is where the architectural difference matters most.
 <Tabs>
  <TabItem label="Poliastro (Python)">
    ```python
    from astropy import units as u
    from astropy.time import Time
    import numpy as np
    from poliastro.bodies import Earth, Mars, Sun
    from poliastro.ephem import Ephem
    from poliastro.iod import izzo
    from poliastro.util import norm
    dep_dates = Time(np.arange(
        Time("2028-08-01").jd,
        Time("2029-01-01").jd,
        1.0  # 1-day steps
    ), format="jd")
    arr_dates = Time(np.arange(
        Time("2029-03-01").jd,
        Time("2029-10-01").jd,
        1.0
    ), format="jd")
    c3_grid = np.full((len(dep_dates), len(arr_dates)), np.nan)
    for i, dep in enumerate(dep_dates):
        earth = Ephem.from_body(Earth, dep)
        r_earth, v_earth = earth.rv(dep)
        for j, arr in enumerate(arr_dates):
            mars = Ephem.from_body(Mars, arr)
            r_mars = mars.rv(arr)[0]
            tof = arr - dep
            if tof.to(u.day).value < 90:
                continue
            try:
                (v_dep, _), = izzo.lambert(Sun.k, r_earth, r_mars, tof)
                v_inf = v_dep - v_earth
                c3_grid[i, j] = norm(v_inf).to(u.km / u.s).value ** 2
            except Exception:
                pass
    # Plot with matplotlib
    import matplotlib.pyplot as plt
    plt.contourf(arr_dates.datetime, dep_dates.datetime, c3_grid, levels=20)
    plt.colorbar(label="C3 (km²/s²)")
    plt.xlabel("Arrival"); plt.ylabel("Departure")
    plt.show()
    ```
    The nested loop iterates over ~23,000 date combinations. Each iteration constructs ephemeris objects, solves Lambert, handles exceptions, and stores the result. With Poliastro's overhead per iteration (unit conversions, object construction), a dense grid can take minutes.
    The plotting integration is the payoff — `matplotlib` is right there, and the result renders inline in a Jupyter notebook.
  </TabItem>
  <TabItem label="pg_orrery (SQL)">
    ```sql
    -- Same grid: 153 departure x 214 arrival = ~32,700 transfers
    SELECT dep::date AS departure,
           arr::date AS arrival,
           round(c3_departure::numeric, 2) AS c3_km2s2,
           round(tof_days::numeric, 0) AS flight_days
    FROM generate_series(
           '2028-08-01'::timestamptz, '2029-01-01'::timestamptz,
           interval '1 day') AS dep,
         generate_series(
           '2029-03-01'::timestamptz, '2029-10-01'::timestamptz,
           interval '1 day') AS arr,
    LATERAL lambert_transfer(3, 4, dep, arr) AS xfer
    WHERE tof_days > 90
      AND c3_departure < 50;
    ```
    ~32,700 Lambert solves. Runs in seconds with automatic parallelism across cores. Export to CSV for plotting:
    ```sql
    COPY (
      SELECT dep::date, arr::date,
             round(c3_departure::numeric, 2) AS c3
      FROM generate_series(
             '2028-08-01'::timestamptz, '2029-01-01'::timestamptz,
             interval '1 day') AS dep,
           generate_series(
             '2029-03-01'::timestamptz, '2029-10-01'::timestamptz,
             interval '1 day') AS arr,
      LATERAL lambert_transfer(3, 4, dep, arr) AS xfer
      WHERE tof_days > 90
    ) TO '/tmp/porkchop.csv' WITH CSV HEADER;
    ```
    Then plot with gnuplot, matplotlib, or any contour tool. The computation and visualization are decoupled — pg_orrery produces data, you render it however you prefer.
  </TabItem>
 </Tabs>
 <Aside type="tip" title="Same algorithm, different throughput">
 Both tools use the Izzo (2015) Lambert solver with Householder iterations. The performance difference comes from execution overhead: Poliastro constructs Astropy unit-annotated objects and ephemeris instances per iteration. pg_orrery calls compiled C functions inline with PostgreSQL's batch execution and parallel query engine.
 </Aside>
 ## Keplerian propagation
 Propagating an orbit forward in time using two-body mechanics. This is how you track comets and asteroids from their osculating orbital elements.
 <Tabs>
  <TabItem label="Poliastro (Python)">
    ```python
    from astropy import units as u
    from astropy.time import Time
    from poliastro.bodies import Sun
    from poliastro.twobody import Orbit
    import numpy as np
    # Define a comet orbit from classical elements
    comet = Orbit.from_classical(
        Sun,
        a=2.7 * u.AU,
        ecc=0.85 * u.one,
        inc=12.0 * u.deg,
        raan=65.0 * u.deg,
        argp=180.0 * u.deg,
        nu=0.0 * u.deg,
        epoch=Time("2025-01-15")
    )
    # Propagate forward 90 days
    future = comet.propagate(90 * u.day)
    r = future.r.to(u.AU)
    print(f"Position: {r}")
    # Time series
    times = np.linspace(0, 365, 100) * u.day
    positions = [comet.propagate(t).r.to(u.AU).value for t in times]
    ```
    Poliastro's `Orbit` object handles the Kepler equation internally. Propagation returns a new `Orbit` with updated elements. The unit-annotated result is self-documenting.
    The time series loop is explicit — each propagation creates a new object.
  </TabItem>
  <TabItem label="pg_orrery (SQL)">
    ```sql
    -- Propagate a comet orbit and observe from Earth
    SELECT t,
           topo_azimuth(obs) AS az,
           topo_elevation(obs) AS el,
           topo_range(obs) / 149597870.7 AS range_au
    FROM generate_series(
           '2025-01-15'::timestamptz,
           '2026-01-15'::timestamptz,
           interval '1 day'
         ) AS t,
    LATERAL comet_observe(
        2.7,           -- semi-major axis (AU)
        0.85,          -- eccentricity
        12.0,          -- inclination (deg)
        65.0,          -- RAAN (deg)
        180.0,         -- argument of perihelion (deg)
        0.0,           -- mean anomaly (deg)
        '2025-01-15'::timestamptz,  -- epoch
        '40.0N 105.3W 1655m'::observer,
        t
    ) AS obs;
    ```
    `comet_observe()` wraps Keplerian propagation with the topocentric observation pipeline. The time series is `generate_series` — one row per timestep, parallel-safe, no loop management.
    For raw heliocentric position without the observation step:
    ```sql
    SELECT t,
           helio_x(kepler_propagate(2.7, 0.85, 12.0, 65.0, 180.0, 0.0,
                                     '2025-01-15'::timestamptz, t)) AS x_au,
           helio_y(kepler_propagate(2.7, 0.85, 12.0, 65.0, 180.0, 0.0,
                                     '2025-01-15'::timestamptz, t)) AS y_au,
           helio_z(kepler_propagate(2.7, 0.85, 12.0, 65.0, 180.0, 0.0,
                                     '2025-01-15'::timestamptz, t)) AS z_au
    FROM generate_series(
           '2025-01-15'::timestamptz,
           '2026-01-15'::timestamptz,
           interval '1 day') AS t;
    ```
  </TabItem>
 </Tabs>
 ## What Poliastro does that pg_orrery doesn't
 Poliastro is a broader orbital mechanics toolkit. Several of its capabilities have no pg_orrery equivalent.
 ### Orbit plotting
 ```python
 from poliastro.plotting import OrbitPlotter3D
 plotter = OrbitPlotter3D()
 plotter.plot(earth_orbit, label="Earth")
 plotter.plot(mars_orbit, label="Mars")
 plotter.plot(transfer, label="Transfer")
 plotter.show()
 ```
 Interactive 3D orbit visualization in Jupyter notebooks. pg_orrery returns coordinates — you build your own visualization externally.
 ### Maneuver planning
 ```python
 from poliastro.maneuver import Maneuver
 hohmann = Maneuver.hohmann(initial_orbit, target_radius)
 bielliptic = Maneuver.bielliptic(initial_orbit, r_b, target_radius)
 print(f"Hohmann delta-v: {hohmann.get_total_cost()}")
 print(f"Bielliptic delta-v: {bielliptic.get_total_cost()}")
 ```
 Poliastro computes delta-v budgets for standard maneuvers. pg_orrery's Lambert solver gives you the transfer orbit but doesn't model the departure/arrival maneuvers.
 ### Perturbation models
 ```python
 from poliastro.twobody.propagation import CowellPropagator
 from poliastro.twobody.perturbations import J2_perturbation
 orbit_with_j2 = orbit.propagate(
    30 * u.day,
    method=CowellPropagator(f=J2_perturbation)
 )
 ```
 Poliastro supports numerical integration with custom perturbation functions — J2 oblateness, atmospheric drag, third-body gravity, solar radiation pressure. pg_orrery's Keplerian propagation is strictly two-body.
 ### Astropy integration
 Poliastro builds on Astropy's unit system, time handling, and coordinate frames. Results carry physical units and can be converted freely. pg_orrery returns bare doubles — you know the units from the documentation, but the database doesn't enforce them.
 ## What pg_orrery does that Poliastro doesn't
 ### SGP4/SDP4 satellite tracking
 ```sql
 -- Observe 12,000 satellites in 17ms
 SELECT s.name,
       topo_elevation(observe(s.tle, '40.0N 105.3W 1655m'::observer, now())) AS el
 FROM satellites s
 WHERE topo_elevation(observe(s.tle, '40.0N 105.3W 1655m'::observer, now())) > 0;
 ```
 Poliastro has no TLE support, no SGP4, no satellite catalog operations. The entire satellite tracking domain — pass prediction, conjunction screening, GiST indexing — is absent.
 ### Orbit determination
 ```sql
 -- Fit an orbit from observations
 SELECT iterations, round(rms_final::numeric, 6) AS rms, status, fitted_tle
 FROM tle_from_angles(
    ra_hours := obs_ra, dec_degrees := obs_dec,
    times := obs_times, obs := station
 );
 ```
 Poliastro does not fit orbits from observations. pg_orrery's differential correction solver handles ECI, topocentric, and angles-only OD with Gauss/Gibbs IOD bootstrap.
 ### Data correlation via SQL
 ```sql
 -- Combine Lambert transfer analysis with mission constraint database
 SELECT t.target, t.dep::date, t.c3,
       m.max_c3_budget, m.launch_vehicle,
       CASE WHEN t.c3 < m.max_c3_budget THEN 'FEASIBLE' ELSE 'OVER BUDGET' END
 FROM transfer_survey t
 JOIN mission_constraints m ON m.target = t.target
 WHERE t.c3 < m.max_c3_budget * 1.1;
 ```
 Poliastro results live in Python memory. Correlating with mission databases, launch vehicle catalogs, or historical data requires loading everything into Python. pg_orrery results are database rows — correlation is a JOIN.
 ### Automatic parallelism
 All pg_orrery functions are `PARALLEL SAFE`. PostgreSQL distributes work across CPU cores automatically. Poliastro runs single-threaded by default — parallelism requires explicit `multiprocessing` or `joblib` setup.
 ## Where Poliastro wins
 <Aside type="note" title="Different tools for different phases">
 Poliastro excels at interactive exploration and visualization. pg_orrery excels at batch computation and data integration. They address different phases of the same workflow.
 </Aside>
 **Interactive exploration.** Poliastro in a Jupyter notebook is the best way to explore orbital mechanics interactively. Plot orbits, change parameters, see results immediately. pg_orrery returns tabular data — useful, but not visual.
 **Unit safety.** Astropy units catch dimension errors at runtime. Passing kilometers where AU are expected raises an exception. pg_orrery uses bare doubles — unit errors are silent.
 **Perturbation modeling.** Poliastro's numerical integrators with custom force models handle real-world orbit analysis that two-body Keplerian propagation cannot.
 **Maneuver analysis.** Delta-v budgets, Hohmann transfers, bielliptic maneuvers — Poliastro models the spacecraft operations side. pg_orrery provides the trajectory geometry.
 **Orbit visualization.** 2D and 3D orbit plots, ground tracks on map projections, interactive widgets in notebooks. pg_orrery produces numbers for external rendering.
 **Broader ephemeris support.** Poliastro accesses any Astropy-compatible ephemeris source, including custom SPK kernels for spacecraft and small bodies.
 ## Where pg_orrery wins
 **Throughput.** 22,500 Lambert solves in 8.3 seconds vs. minutes in Poliastro for the same grid. The difference is compiled C + parallel execution vs. Python object construction per iteration.
 **Satellite operations.** SGP4, TLEs, pass prediction, conjunction screening, orbit determination — Poliastro has none of these. If you work with artificial satellites, pg_orrery covers the domain.
 **Data integration.** Results live in PostgreSQL alongside your other data. JOIN with mission databases, observation logs, or any other table.
 **Reproducibility.** A SQL query is a complete specification. No virtual environment, no package versions, no imports. Share the query and it runs identically on any pg_orrery installation.
 **Deployment simplicity.** pg_orrery is a PostgreSQL extension — install once, available to every client. Poliastro requires a Python environment with Astropy, NumPy, SciPy, and their transitive dependencies.
 ## A combined workflow
 The two tools complement each other naturally:
 <Steps>
 1. **Explore with Poliastro.** Use a Jupyter notebook to understand the problem — plot orbits, visualize transfer geometry, prototype maneuver sequences. This is the design phase.
 2. **Survey with pg_orrery.** Once you know what you're looking for, run the parameter sweep in SQL. 100,000 Lambert solves to identify the best launch windows. Store the results in a table.
 3. **Correlate with pg_orrery.** JOIN transfer results with launch vehicle constraints, budget timelines, or historical mission data. Filter to feasible solutions.
 4. **Refine with Poliastro.** Take the best candidates back to Python for detailed analysis — add perturbations, model the departure spiral, compute the precise delta-v budget.
 5. **Operate with pg_orrery.** Once you have TLEs (from catalog data or from pg_orrery's OD solver), the operational tracking — observation, pass prediction, conjunction screening — stays in SQL.
 </Steps>
 <Aside type="tip" title="Connecting the two">
 Export pg_orrery results to CSV or JSON for import into Poliastro notebooks:
 ```sql
 COPY (
  SELECT dep::date, arr::date,
         c3_departure, c3_arrival, tof_days, transfer_sma
  FROM lambert_transfer(3, 4, dep, arr)
  -- ... your survey grid ...
 ) TO '/tmp/survey.csv' WITH CSV HEADER;
 ```
 ```python
 import pandas as pd
 survey = pd.read_csv('/tmp/survey.csv')
 best = survey.nsmallest(10, 'c3_departure')
 # Now explore the best candidates in Poliastro
 ```
 The boundary between "bulk computation" and "interactive analysis" is a CSV file.
 </Aside>
--- a/docs/src/content/docs/workflow/from-skyfield.mdx
+++ b/docs/src/content/docs/workflow/from-skyfield.mdx
@ -273,7 +273,7 @@ Predict when a satellite will be visible from a location. This is where Skyfield
 pg_orrery does not replace Skyfield for all use cases. Be clear about where the trade-offs fall.
 </Aside>
-**Apparent-position corrections.** Skyfield uses the full IAU 2000A nutation model, polar motion corrections, delta-T from IERS data, and iterates for light-time and stellar aberration. pg_orrery v0.3.0 can [optionally use DE441](/guides/de-ephemeris/) for the same underlying geometric accuracy (~0.1 milliarcsecond), but Skyfield still applies corrections that pg_orrery does not — corrections that matter for precision apparent-coordinate work like occultation timing or sub-arcsecond astrometry.
+**Apparent-position corrections.** Skyfield uses the full IAU 2000A nutation model, polar motion corrections, delta-T from IERS data, and iterates for light-time and stellar aberration. Since v0.3.0, pg_orrery can [optionally use DE441](/guides/de-ephemeris/) for the same underlying geometric accuracy (~0.1 milliarcsecond), but Skyfield still applies corrections that pg_orrery does not — corrections that matter for precision apparent-coordinate work like occultation timing or sub-arcsecond astrometry.
 **Rise/set finding.** `find_events()` uses numerical root-finding to pinpoint the exact moment a body crosses an elevation threshold. pg_orrery's `predict_passes` uses a step-and-refine approach that's fast for batches but less precise for individual events.
@ -295,6 +295,42 @@ pg_orrery does not replace Skyfield for all use cases. Be clear about where the
 **Reproducibility.** A SQL query is a complete, self-contained specification of a computation. No virtual environment, no package versions, no file paths. The same query produces the same result on any PostgreSQL instance with pg_orrery installed.
 **Orbit determination.** Skyfield is a propagation and observation library — it does not fit orbits from observations. Since v0.4.0, pg_orrery can [fit TLEs from ECI, topocentric, or angles-only observations](/guides/orbit-determination/) via differential correction, entirely in SQL. This closes the loop: observe a satellite, fit its orbit, and propagate the result — all within the database.
 ```sql
 -- The closed-loop workflow Skyfield can't do:
 -- observe → fit (with range rate + weighted obs) → catalog → predict
 WITH obs_data AS (
    -- Collect observations from tonight's tracking pass
    SELECT array_agg(observe(seed_tle, station, t) ORDER BY t) AS topo_obs,
           array_agg(t ORDER BY t) AS times
    FROM generate_series(
             '2024-01-01 12:00+00'::timestamptz,
             '2024-01-01 12:10+00'::timestamptz,
             interval '30 seconds') t
 ),
 fit AS (
    SELECT fitted_tle, rms_final, condition_number, status
    FROM obs_data,
         tle_from_topocentric(
             obs_data.topo_obs, obs_data.times,
             '40.0N 105.3W 1655m'::observer,
             seed := seed_tle,
             fit_range_rate := true,   -- v0.6.0: use Doppler data
             weights := obs_weights    -- v0.6.0: per-obs weighting
         )
    WHERE status = 'converged'
 )
 -- Immediately predict passes with the freshly-fitted TLE
 SELECT pass_aos(p) AS rise,
       pass_max_el(p) AS max_el,
       pass_los(p) AS set
 FROM fit,
     predict_passes(fit.fitted_tle,
                    '40.0N 105.3W 1655m'::observer,
                    now(), now() + interval '48 hours', 10.0) p;
 ```
 ## Migrating gradually
 You don't have to choose one or the other. A practical migration path:
@ -308,5 +344,7 @@ You don't have to choose one or the other. A practical migration path:
 4. **Move reporting to SQL.** "What was above 20 degrees from each of our 5 observers last night?" is a single query with a CROSS JOIN, not a Python loop over observers and timestamps.
-5. **Enable DE when accuracy matters.** If you find VSOP87's ~1 arcsecond isn't enough for a specific use case, [configure a DE file](/guides/de-ephemeris/) and switch to `_de()` function variants. Same SQL patterns, same parameters — just add `_de` to the function name.
+5. **Move orbit determination to SQL.** If you're fitting orbits with separate Python tools (find_orb, OREKit bindings, custom scripts), pg_orrery's [OD solver](/guides/orbit-determination/) can do it alongside your existing satellite catalog. Observations in, TLE out — no data pipeline between Python and PostgreSQL.
 6. **Enable DE when accuracy matters.** If you find VSOP87's ~1 arcsecond isn't enough for a specific use case, [configure a DE file](/guides/de-ephemeris/) and switch to `_de()` function variants. Same SQL patterns, same parameters — just add `_de` to the function name.
 </Steps>
--- a/docs/src/content/docs/workflow/sql-advantage.mdx
+++ b/docs/src/content/docs/workflow/sql-advantage.mdx
@ -344,7 +344,7 @@ This is a screening filter, not a precision conjunction analysis. It identifies
 ## Provider switching: accuracy when you need it
-pg_orrery v0.3.0 has two ephemeris providers — the built-in VSOP87 pipeline (~1 arcsecond) and optional [JPL DE440/441](/guides/de-ephemeris/) (~0.1 milliarcsecond). The SQL interface makes switching between them a one-character change.
+pg_orrery has two ephemeris providers (since v0.3.0) — the built-in VSOP87 pipeline (~1 arcsecond) and optional [JPL DE440/441](/guides/de-ephemeris/) (~0.1 milliarcsecond). The SQL interface makes switching between them a one-character change.
 ```sql
 -- VSOP87 (built-in, IMMUTABLE, no setup)
@ -437,6 +437,104 @@ ORDER BY night;
 For each night, this finds which planet reaches the highest elevation — useful for deciding what to observe. The window function handles the ranking within each night without a correlated subquery.
 ## Orbit determination: observations in, TLE out
 Since v0.4.0, pg_orrery can fit TLE orbital elements from observations using differential correction. This closes a loop that previously required external tools — you can observe, fit, and track entirely within SQL.
 ### Seed-free fitting with Gauss IOD
 The Gauss method (v0.6.0) bootstraps an initial orbit from three angles-only observations, eliminating the need for a seed TLE from a catalog:
 ```sql
 -- Unknown satellite: RA/Dec observations from last night
 -- No seed TLE needed — Gauss IOD handles the initial guess
 SELECT iterations,
       round(rms_final::numeric, 6) AS rms_rad,
       status,
       round(condition_number::numeric, 1) AS cond,
       fitted_tle
 FROM tle_from_angles(
    ra_hours := ARRAY[12.345, 12.567, 12.789, 13.012, 13.234, 13.456],
    dec_degrees := ARRAY[45.1, 44.8, 44.3, 43.6, 42.8, 41.9],
    times := ARRAY[
        '2024-01-01 12:00+00', '2024-01-01 12:01+00',
        '2024-01-01 12:02+00', '2024-01-01 12:03+00',
        '2024-01-01 12:04+00', '2024-01-01 12:05+00'
    ]::timestamptz[],
    obs := '40.0N 105.3W 1655m'::observer
 );
 ```
 The fitted TLE is immediately usable with `observe()`, `predict_passes()`, and the GiST conjunction index. No export, no format conversion.
 ### Weighted multi-sensor fusion
 When observations come from sensors with different accuracies, the `weights` parameter (v0.6.0) gives explicit control over each observation's influence:
 ```sql
 -- Radar + optical: weight radar higher, include Doppler data
 SELECT iterations,
       round(rms_final::numeric, 4) AS rms_km,
       status,
       nstate AS params_fitted
 FROM tle_from_topocentric(
    observations := all_topo_obs,
    times := all_times,
    observers := ARRAY[
        '40.0N 105.3W 1655m',   -- radar
        '34.1N 118.3W 100m'     -- optical
    ]::observer[],
    observer_ids := station_ids,
    seed := seed_tle,
    fit_range_rate := true,    -- include Doppler from radar
    weights := ARRAY[
        1.0, 1.0, 1.0, 1.0,   -- radar: full weight
        0.3, 0.3, 0.3, 0.3    -- optical: down-weighted
    ]
 );
 ```
 ### The full pipeline as SQL
 The OD pattern composes with every other pattern on this page:
 ```sql
 -- Full pipeline: batch OD → catalog insert → pass prediction
 -- Process all objects observed tonight
 WITH fits AS (
    SELECT obj_id,
           tle_from_angles(
               ra_hours := array_agg(ra ORDER BY obs_time),
               dec_degrees := array_agg(dec ORDER BY obs_time),
               times := array_agg(obs_time ORDER BY obs_time),
               obs := '40.0N 105.3W 1655m'::observer
           ) AS od
    FROM tonight_observations
    GROUP BY obj_id
    HAVING count(*) >= 3
 ),
 new_tles AS (
    INSERT INTO satellites (norad_id, name, tle, source)
    SELECT 90000 + row_number() OVER (),
           obj_id,
           (od).fitted_tle,
           'local_od'
    FROM fits
    WHERE (od).status = 'converged'
    RETURNING norad_id, tle
 )
 -- Immediately predict passes for every newly cataloged object
 SELECT s.norad_id,
       pass_aos(p) AS rise,
       pass_max_el(p) AS max_el
 FROM new_tles s,
     predict_passes(s.tle, '40.0N 105.3W 1655m'::observer,
                    now(), now() + interval '48 hours', 10.0) p
 ORDER BY pass_aos(p);
 ```
 Observation ingestion, orbit fitting, catalog insertion, and pass prediction — one statement. Schedule it with `pg_cron` and the pipeline runs itself.
 ## Composition: building complex queries from simple parts
 The real power of SQL is that these patterns compose. A single query can use `generate_series` for time steps, `CROSS JOIN` for parameter sweeps, `JOIN` for data correlation, window functions for change detection, and `COPY TO` for export — all in one statement.
@ -499,4 +597,4 @@ This query:
 In a traditional workflow, each of these steps would be a separate script, a separate data file, and a separate tool. In SQL, they compose into a single declarative statement that the database engine optimizes and parallelizes.
-That's the advantage. Not that SQL is a better programming language — it isn't. But for the specific pattern of "evaluate a function over structured parameter spaces and correlate the results with existing data," SQL is exactly the right tool. pg_orrery puts 68 functions inside that tool — from 17ms satellite batch propagation to sub-milliarcsecond DE441 planet positions — and every one of them composes with every SQL pattern on this page.
+That's the advantage. Not that SQL is a better programming language — it isn't. But for the specific pattern of "evaluate a function over structured parameter spaces and correlate the results with existing data," SQL is exactly the right tool. pg_orrery puts over 100 functions inside that tool — from 17ms satellite batch propagation to sub-milliarcsecond DE441 planet positions — and every one of them composes with every SQL pattern on this page.