Comets and Asteroids

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

Two 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

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.

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_orbit 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_orbit returns position only. Comet brightness depends on heliocentric distance, geocentric distance, and a magnitude slope parameter that varies per comet.
• No orbit determination. pg_orbit 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.

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