From find_orb to SQL

find_orb is Bill Gray’s orbit determination software — the same Bill Gray whose sat_code SGP4 library 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.

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.

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.

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.

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

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

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.

find_orb

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

pg_orrery (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.

Multi-observer observations

Observations from multiple stations improve orbit geometry. Both tools support this, but the workflow differs.

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.

pg_orrery (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.

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.

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.

pg_orrery (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.

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

Residual inspection

After fitting, you want to know which observations are good and which are outliers.

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.

pg_orrery (SQL)

For ECI-based fitting, tle_fit_residuals() returns per-observation residuals:

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

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

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:

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

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.

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

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.