Comets and Asteroids

pg_orrery propagates comets and asteroids using two-body Keplerian mechanics. You provide six classical orbital elements from the Minor Planet Center (MPC) or any other source, and pg_orrery computes the body’s heliocentric position at any time. Combined with Earth’s position from VSOP87, you can observe the body from any location on Earth.

How you do it today

Tracking comets and asteroids typically involves:

• JPL Small-Body Database (SBDB): Look up an object, get its orbital elements, request an ephemeris. One object at a time, web-based.
• Find_Orb: Fit orbits from observations and propagate forward. Powerful but desktop-only, primarily for orbit determination rather than ephemeris computation.
• Skyfield: Can propagate comets from MPC elements, but you need to load a planetary ephemeris for the Earth’s position and write the observation pipeline yourself.
• Minor Planet Center (MPC): Publishes orbital elements for over 1.3 million objects. Getting batch ephemerides means downloading elements and running them through your own propagation code.

The pattern is familiar: download elements, propagate in Python or C, transform to observer coordinates, and import results into your database.

What changes with pg_orrery

Five functions handle comet/asteroid computation:

Function	What it does
kepler_propagate(q, e, i, omega, Omega, T, time)	Propagates orbital elements to a heliocentric position (AU)
comet_observe(q, e, i, omega, Omega, T, ex, ey, ez, observer, time)	Full observation pipeline: propagate + geocentric transform + topocentric
oe_from_mpc(line)	Parses one MPCORB.DAT line into an orbital_elements type
small_body_heliocentric(oe, time)	Heliocentric position from bundled elements
small_body_observe(oe, observer, time)	Topocentric observation — auto-fetches Earth via VSOP87

kepler_propagate() solves Kepler’s equation for elliptic (e < 1), parabolic (e = 1), and hyperbolic (e > 1) orbits. The solver handles all three cases with appropriate numerical methods.

comet_observe() wraps the full chain: propagate the comet’s position, subtract Earth’s heliocentric position, and transform to horizon coordinates. You supply Earth’s position as three floats (ecliptic J2000, AU) because you might want to compute it once and reuse it across many comets.

small_body_observe() (v0.8.0) does the same thing but fetches Earth’s position automatically — you just pass the orbital_elements type and an observer. See the orbital_elements type section below.

The parameters map directly to MPC orbital element format:

Parameter	MPC field	Units
q	Perihelion distance	AU
e	Eccentricity	dimensionless
i	Inclination	degrees
omega	Argument of perihelion	degrees
Omega	Longitude of ascending node	degrees
T	Perihelion time	Julian date

What pg_orrery does not replace

Two-body limitations

Keplerian propagation assumes the body is influenced only by the Sun. Real small bodies experience planetary perturbations that accumulate over time.

• No perturbations. Jupiter alone can shift a comet’s position by degrees over a few years. Two-body propagation is most accurate near perihelion, within a few months of the elements’ epoch.
• No non-gravitational forces. Comet outgassing produces accelerations not captured by Keplerian mechanics. For long-period comets far from the Sun, this is negligible. For short-period comets near perihelion, it matters.
• No magnitude estimation. pg_orrery returns position only. Comet brightness depends on heliocentric distance, geocentric distance, and a magnitude slope parameter that varies per comet.
• No orbit determination. pg_orrery propagates known orbits. It does not fit orbits from observations.

For MPC elements less than a few months old, two-body propagation is typically accurate to a few arcminutes for asteroids and tens of arcminutes for comets. Fresh elements give better results.

The orbital_elements type

The raw-parameter functions (kepler_propagate, comet_observe) work well when you have elements in variables or a table with individual columns. But they require passing 6–11 float8 arguments per call, and comet_observe requires you to manually fetch Earth’s position.

The orbital_elements type (v0.8.0) bundles all nine classical elements into a single 72-byte PostgreSQL datum:

-- Construct from a literal
SELECT '(2460605.5,2.5478,0.0789126,10.58664,73.42937,80.2686,2460319.0,3.33,0.12)'::orbital_elements;

-- Or parse directly from the MPC catalog
SELECT oe_from_mpc(
  '00001    3.33  0.12 K24AM  60.07966  73.42937   80.26860  10.58664  0.0789126  0.21406048   2.7660961  0 MPO838504  8738 115 1801-2024 0.65 M-v 30k MPCLINUX   0000      (1) Ceres              20240825'
);

With bundled elements, observation becomes a single function call:

-- Before (comet_observe): 11 arguments, manual Earth fetch
WITH earth AS (SELECT planet_heliocentric(3, now()) AS h)
SELECT topo_elevation(comet_observe(
    2.5478, 0.0789, 10.59, 73.43, 80.27, 2460319.0,
    helio_x(h), helio_y(h), helio_z(h),
    '40.0N 105.3W 1655m'::observer, now()))
FROM earth;

-- After (small_body_observe): 3 arguments, Earth auto-fetched
SELECT topo_elevation(small_body_observe(
    oe, '40.0N 105.3W 1655m'::observer, now()))
FROM asteroid_catalog;

Load and query the MPC catalog

The MPC publishes MPCORB.DAT — orbital elements for every numbered asteroid. Here’s how to load it into PostgreSQL:

1. Create a table with an orbital_elements column:

   CREATE TABLE asteroids (
       number      int PRIMARY KEY,
       name        text,
       oe          orbital_elements NOT NULL
   );

2. Load via a staging table:

   -- Stage the raw text lines
   CREATE TEMP TABLE mpc_raw (line text);
   \copy mpc_raw FROM 'MPCORB.DAT'
   
   -- Parse into orbital_elements, extract number and name
   INSERT INTO asteroids (number, name, oe)
   SELECT substring(line FROM 1 FOR 7)::int,
          trim(substring(line FROM 167 FOR 30)),
          oe_from_mpc(line)
   FROM mpc_raw
   WHERE length(line) >= 103
     AND substring(line FROM 1 FOR 7) ~ '^\s*\d+$';
   
   DROP TABLE mpc_raw;

3. Query: what asteroids are above 20 degrees tonight?

   SELECT name, number,
          round(topo_elevation(t)::numeric, 1) AS el,
          round(topo_azimuth(t)::numeric, 1)   AS az,
          round((topo_range(t) / 149597870.7)::numeric, 2) AS dist_au
   FROM asteroids,
        small_body_observe(oe, '40.0N 105.3W 1655m'::observer, now()) AS t
   WHERE topo_elevation(t) > 20
   ORDER BY topo_elevation(t) DESC
   LIMIT 20;

4. Query: heliocentric distance of Ceres over 6 months:

   SELECT t::date AS date,
          round(helio_distance(
              small_body_heliocentric(oe, t))::numeric, 4) AS dist_au
   FROM asteroids,
        generate_series(
            '2025-01-01'::timestamptz,
            '2025-07-01'::timestamptz,
            interval '15 days'
        ) AS t
   WHERE number = 1;

Tip

For batch observation at a single time, comet_observe() is still more efficient — it lets you compute Earth’s VSOP87 position once with planet_heliocentric(3, t) and reuse it across all objects. small_body_observe() re-fetches Earth on every call. For interactive single-object queries, small_body_observe() is simpler.

Try it

Circular orbit sanity check

A body in a circular orbit at 1 AU with all angles zero should return to its starting position after one year:

-- At perihelion (T=0), position should be (1, 0, 0) AU
SELECT round(helio_x(kepler_propagate(
           1.0, 0.0, 0.0, 0.0, 0.0,
           2451545.0,  -- J2000.0
           '2000-01-01 12:00:00+00'))::numeric, 6) AS x,
       round(helio_y(kepler_propagate(
           1.0, 0.0, 0.0, 0.0, 0.0,
           2451545.0,
           '2000-01-01 12:00:00+00'))::numeric, 6) AS y;

At time = perihelion, the position is exactly (q, 0, 0) in the orbital plane. After a quarter orbit (~91 days), it moves to approximately (0, 1, 0).

Eccentric elliptic orbit

An orbit with e=0.5 and q=0.5 AU has a semi-major axis of 1.0 AU and the same period as Earth, but a very different shape:

-- Position over one orbit
SELECT t::date AS date,
       round(helio_distance(kepler_propagate(
           0.5, 0.5, 0.0, 0.0, 0.0,
           2451545.0, t))::numeric, 4) AS dist_au
FROM generate_series(
    '2000-01-01 12:00:00+00'::timestamptz,
    '2001-01-01 12:00:00+00'::timestamptz,
    interval '30 days'
) AS t;

The distance ranges from 0.5 AU (perihelion) to 1.5 AU (aphelion). This is the classic comet behavior: fast and close to the Sun at perihelion, slow and distant at aphelion.

Inclined orbit

Orbital inclination rotates the orbital plane out of the ecliptic:

-- A polar orbit (i=90 deg) at 1 AU
SELECT round(helio_x(kepler_propagate(
           1.0, 0.0, 90.0, 0.0, 0.0,
           2451545.0,
           '2000-01-01 12:00:00+00'))::numeric, 6) AS x,
       round(helio_z(kepler_propagate(
           1.0, 0.0, 90.0, 0.0, 0.0,
           2451545.0,
           '2000-01-01 12:00:00+00'))::numeric, 6) AS z;

At perihelion, the position is still along the node line (x-axis) regardless of inclination. The inclination only shows up when the body moves away from the node.

Hyperbolic orbit

Interstellar objects like ‘Oumuamua travel on hyperbolic trajectories (e > 1):

-- Hyperbolic orbit: e=1.5, q=1.0 AU
-- At perihelion
SELECT round(helio_x(kepler_propagate(
           1.0, 1.5, 0.0, 0.0, 0.0,
           2451545.0,
           '2000-01-01 12:00:00+00'))::numeric, 6) AS x_at_perihelion;

-- 6 months later: body is receding rapidly
SELECT round(helio_distance(kepler_propagate(
           1.0, 1.5, 0.0, 0.0, 0.0,
           2451545.0,
           '2000-07-01 12:00:00+00'))::numeric, 2) AS dist_6mo;

The body approaches from infinity, swings past the Sun at perihelion distance, and departs on a hyperbola.

Near-parabolic comet

Many long-period comets have eccentricities very close to 1.0. pg_orrery handles the parabolic case (e=1.0 exactly) with a dedicated Barker equation solver:

SELECT round(helio_x(kepler_propagate(
           1.0, 1.0, 0.0, 0.0, 0.0,
           2451545.0,
           '2000-01-01 12:00:00+00'))::numeric, 6) AS x_parabolic;

Track a comet with real MPC elements

Here is how you would observe a comet using elements from the Minor Planet Center. This example uses hypothetical elements for a Halley-type orbit:

1. Get Earth’s heliocentric position at the observation time:

   SELECT helio_x(planet_heliocentric(3, '2024-06-15 04:00:00+00')) AS ex,
          helio_y(planet_heliocentric(3, '2024-06-15 04:00:00+00')) AS ey,
          helio_z(planet_heliocentric(3, '2024-06-15 04:00:00+00')) AS ez;

2. Observe the comet using all parameters together:

   WITH earth AS (
       SELECT planet_heliocentric(3, '2024-06-15 04:00:00+00') AS pos
   )
   SELECT round(topo_azimuth(comet_observe(
              0.587, 0.967, 162.3, 111.3, 58.4, 2446467.4,
              helio_x(pos), helio_y(pos), helio_z(pos),
              '40.0N 105.3W 1655m'::observer,
              '2024-06-15 04:00:00+00'))::numeric, 1) AS az,
          round(topo_elevation(comet_observe(
              0.587, 0.967, 162.3, 111.3, 58.4, 2446467.4,
              helio_x(pos), helio_y(pos), helio_z(pos),
              '40.0N 105.3W 1655m'::observer,
              '2024-06-15 04:00:00+00'))::numeric, 1) AS el,
          round(topo_range(comet_observe(
              0.587, 0.967, 162.3, 111.3, 58.4, 2446467.4,
              helio_x(pos), helio_y(pos), helio_z(pos),
              '40.0N 105.3W 1655m'::observer,
              '2024-06-15 04:00:00+00'))::numeric, 0) AS range_km
   FROM earth;

   The orbital elements are: q=0.587 AU, e=0.967, i=162.3 deg, omega=111.3 deg, Omega=58.4 deg, T=JD 2446467.4.

Store a comet catalog

For batch operations, store orbital elements in a table:

CREATE TABLE comets (
    designation     text PRIMARY KEY,
    name            text,
    q_au            float8 NOT NULL,
    eccentricity    float8 NOT NULL,
    inclination_deg float8 NOT NULL,
    arg_peri_deg    float8 NOT NULL,
    long_node_deg   float8 NOT NULL,
    perihelion_jd   float8 NOT NULL,
    epoch_jd        float8
);

-- Insert a few examples
INSERT INTO comets VALUES
    ('1P',   'Halley',    0.587,  0.967, 162.3, 111.3, 58.4,  2446467.4, 2446480.0),
    ('2P',   'Encke',     0.336,  0.847,  11.8, 186.5, 334.6, 2460585.0, 2460600.0),
    ('67P',  'C-G',       1.243,  0.641,   7.0,  12.8,  50.1, 2457257.0, 2457260.0);

Batch observe all comets

WITH earth AS (
    SELECT planet_heliocentric(3, '2024-06-15 04:00:00+00') AS pos
)
SELECT c.name,
       round(topo_azimuth(obs)::numeric, 1) AS az,
       round(topo_elevation(obs)::numeric, 1) AS el,
       round(topo_range(obs)::numeric, 0) AS range_km
FROM comets c, earth,
     LATERAL comet_observe(
         c.q_au, c.eccentricity, c.inclination_deg,
         c.arg_peri_deg, c.long_node_deg, c.perihelion_jd,
         helio_x(earth.pos), helio_y(earth.pos), helio_z(earth.pos),
         '40.0N 105.3W 1655m'::observer,
         '2024-06-15 04:00:00+00') obs
WHERE topo_elevation(obs) > 0
ORDER BY topo_elevation(obs) DESC;

This observes every comet in the catalog and filters to those above the horizon. The Earth position is computed once and reused for all comets.

Heliocentric distance over time

Track how a comet’s distance from the Sun changes through its orbit:

SELECT t::date AS date,
       round(helio_distance(kepler_propagate(
           0.336, 0.847, 11.8, 186.5, 334.6,
           2460585.0, t))::numeric, 3) AS encke_dist_au
FROM generate_series(
    '2024-01-01'::timestamptz,
    '2024-12-31'::timestamptz,
    interval '15 days'
) AS t;

Comet Encke (q=0.336 AU, e=0.847) ranges from 0.336 AU at perihelion to about 4.1 AU at aphelion. Its 3.3-year period means it passes through the inner solar system frequently.

Verify: distance conservation for circular orbit

A useful sanity check. A circular orbit should maintain constant heliocentric distance:

SELECT t::date AS date,
       round(helio_distance(kepler_propagate(
           1.0, 0.0, 0.0, 0.0, 0.0,
           2451545.0, t))::numeric, 6) AS dist_au
FROM generate_series(
    '2000-01-01'::timestamptz,
    '2001-01-01'::timestamptz,
    interval '60 days'
) AS t;

Every row should read 1.000000 AU. If it does, the Kepler solver is working correctly.