Working with orbit results (Gauss IOD)

This tutorial shows how to consume the successful results returned by TrajectorySet.estimate_all_orbits(...) and how to navigate the different orbital element families produced by the Gauss solver.

You will learn how to:

• iterate over the successful results map and inspect GaussResult objects,
• check whether the result is a preliminary or corrected orbit,
• extract the concrete orbital elements (Keplerian, Equinoctial, or Cometary),
• convert between element families when supported,
• serialize results to dictionaries for logging or downstream processing.

The examples assume you already ran batch IOD and obtained (ok, errors) from estimate_all_orbits(env, params, ...). See the Trajectories and IODParams tutorials for how to configure and run the solver.

Reference: run_iod() used by the snippets

The code examples below call a shared helper that builds a small dataset, runs batch IOD, and returns (ok, errors). For completeness, here is the helper once:

from py_outfit import PyOutfit, IODParams, TrajectorySet, Observer
import numpy as np
from astropy.time import Time


def run_iod():
    env = PyOutfit("horizon:DE440", "FCCT14")
    obs = Observer(
        0.0, 0.0, 1.0, "DemoSite", np.deg2rad(0.3 / 3600.0), np.deg2rad(0.3 / 3600.0)
    )
    env.add_observer(obs)

    trajectory_id = np.array(
        [
            0,
            1,
            2,
            1,
            2,
            1,
            0,
            0,
            0,
            1,
            2,
            1,
            1,
            0,
            2,
            2,
            0,
            2,
            2,
            10,
            10,
            10,
            11,
            11,
            11,
        ],
        dtype=np.uint32,
    )
    ra_deg = np.array(
        [
            20.9191548,
            33.4247141,
            32.1435128,
            33.4159091,
            32.1347282,
            33.3829299,
            20.6388309,
            20.6187259,
            20.6137886,
            32.7525147,
            31.4874917,
            32.4518231,
            32.4495403,
            19.8927380,
            30.6416348,
            30.0938936,
            18.2218784,
            28.3859403,
            28.3818327,
            10.0,
            10.01,
            10.02,
            180.0,
            180.02,
            180.05,
        ]
    )
    dec_deg = np.array(
        [
            20.0550441,
            23.5516817,
            26.5139615,
            23.5525348,
            26.5160622,
            23.5555991,
            20.1218532,
            20.1264229,
            20.1275173,
            23.6064063,
            26.6622284,
            23.6270392,
            23.6272157,
            20.2977473,
            26.8303010,
            26.9256271,
            20.7096409,
            27.1602652,
            27.1606420,
            5.0,
            5.01,
            5.015,
            -10.0,
            -10.02,
            -10.03,
        ]
    )
    times_jd_utc = np.array(
        [
            2458789.6362963,
            2458789.6381250,
            2458789.6381250,
            2458789.6663773,
            2458789.6663773,
            2458789.7706481,
            2458790.6995023,
            2458790.7733333,
            2458790.7914120,
            2458791.8445602,
            2458791.8445602,
            2458792.8514699,
            2458792.8590741,
            2458793.6896759,
            2458794.7996759,
            2458796.7965162,
            2458801.7863426,
            2458803.7699537,
            2458803.7875231,
            2458800.0,
            2458800.01,
            2458800.03,
            2458800.0,
            2458800.02,
            2458800.05,
        ]
    )

    # Convert times to MJD (TT) using astropy
    t_utc = Time(times_jd_utc, format="jd", scale="utc")
    mjd_tt = t_utc.tt.mjd.astype(np.float64)

    # Degree path performs a single conversion to radians at ingestion
    ts = TrajectorySet.from_numpy_degrees(
        env,
        trajectory_id,
        ra_deg,
        dec_deg,
        error_ra_arcsec=0.3,
        error_dec_arcsec=0.3,
        mjd_tt=mjd_tt,
        observer=obs,
    )

    # Configure IOD and run batch estimation
    params = IODParams.builder().max_triplets(200).do_sequential().build()
    ok, errors = ts.estimate_all_orbits(env, params, seed=42)

    return ok, errors

---

Iterate over successful results

# Iterate over successful Gauss results and print stage + RMS
from common_tuto import run_iod

ok, errors = run_iod()

if not ok:
    print("No successful results; errors:", errors)
else:
    for obj_id, (g_res, rms) in ok.items():
        stage = "corrected" if g_res.is_corrected() else "preliminary"
        print(f"Object {obj_id}: stage={stage}, RMS={rms:.6e} rad")

• obj_id is the same identifier you used when ingesting trajectories (int or str).
• rms is the post-fit residual RMS (radians) computed over the chosen time window.

---

Determine the element family

# Determine element family and stage from one successful result
from common_tuto import run_iod

ok, errors = run_iod()

if ok:
    obj_id, (g_res, rms) = next(iter(ok.items()))
    fam = g_res.elements_type()
    stage = "corrected" if g_res.is_corrected() else "preliminary"
    print(obj_id, fam, stage)
else:
    print("No successful results; cannot show family/stage example.")

---

Extract concrete elements

Use the typed accessors; they return None if the stored family differs.

# Extract concrete orbital elements from a GaussResult
from common_tuto import run_iod

ok, _ = run_iod()

if ok:
    _, (g_res, rms) = next(iter(ok.items()))
    k = g_res.keplerian()
    q = g_res.equinoctial()
    c = g_res.cometary()
    if k is not None:
        print("K a,e:", k.semi_major_axis, k.eccentricity)
    if q is not None:
        print("Q a,h,k:", q.semi_major_axis, q.eccentricity_sin_lon, q.eccentricity_cos_lon)
    if c is not None:
        print("C q,e:", c.perihelion_distance, c.eccentricity)
else:
    print("No successful results; cannot extract element examples.")

Units reminder: - Epochs are MJD (TDB). Angles are radians. Distances are AU.

---

Convert between element families

Conversions are provided by the element classes themselves.

• Keplerian → Equinoctial:

# Convert between orbital element families
from common_tuto import run_iod

ok, _ = run_iod()

if ok:
    _, (g_res, _) = next(iter(ok.items()))
    k = g_res.keplerian()
    q = g_res.equinoctial()
    c = g_res.cometary()

    if k is not None:
        q2 = k.to_equinoctial()
        print("K→Q a,λ:", q2.semi_major_axis, q2.mean_longitude)
    if q is not None:
        k2 = q.to_keplerian()
        print("Q→K a,e:", k2.semi_major_axis, k2.eccentricity)
    if c is not None and c.eccentricity > 1.0:
        k_h = c.to_keplerian()
        print("C(hyperbolic)→K a,e:", k_h.semi_major_axis, k_h.eccentricity)
else:
    print("No successful results; cannot demonstrate conversions.")

Note: parabolic cometary elements (e = 1) cannot be converted by these helpers and will raise a ValueError.

---

Structured dict serialization

Every GaussResult can be converted to a plain dictionary for easy logging and JSON export:

# Serialize a GaussResult to a structured dict
from common_tuto import run_iod

ok, _ = run_iod()

if ok:
    _, (g_res, _) = next(iter(ok.items()))
    d = g_res.to_dict()
    print(d["stage"], d["type"], sorted(d["elements"].keys()))
else:
    print("No successful results; cannot show to_dict().")

Example for a Keplerian result:

{
    'stage': 'corrected', 
    'type': 'keplerian', 
    'elements': {
        'reference_epoch': 58794.29503864708, 
        'semi_major_axis': 2.618543557694562, 
        'eccentricity': 0.2917924222538649, 
        'inclination': 0.23168624097364912, 
        'ascending_node_longitude': 0.20856161706357348, 
        'periapsis_argument': 6.264575557486691, 
        'mean_anomaly': 0.29001350766154466
    }
}

---

Putting it together: filter, convert, export

Below is a compact pattern you can adapt to your pipeline:

# Filter, convert, and build a plain export from successful results
from common_tuto import run_iod

ok, _ = run_iod()

export = []
for obj_id, (g_res, rms) in ok.items():
    k = g_res.keplerian()
    if k is None:
        q = g_res.equinoctial()
        if q is not None:
            k = q.to_keplerian()
        else:
            c = g_res.cometary()
            if c is not None and c.eccentricity > 1.0:
                k = c.to_keplerian()

    if k is None:
        d = g_res.to_dict()
        export.append({"id": obj_id, "rms": rms, **d})
    else:
        export.append({
            "id": obj_id,
            "rms": rms,
            "stage": "corrected" if g_res.is_corrected() else "preliminary",
            "type": "keplerian",
            "a_au": k.semi_major_axis,
            "e": k.eccentricity,
            "i_rad": k.inclination,
            "Omega_rad": k.ascending_node_longitude,
            "omega_rad": k.periapsis_argument,
            "M_rad": k.mean_anomaly,
            "epoch_mjd_tdb": k.reference_epoch,
        })

print(len(export))

---

Tips

• Always check the element family via elements_type() before calling accessors; the typed helpers return None when mismatched.
• When you need a single canonical representation, prefer converting to Keplerian where defined, but keep native cometary elements for e = 1.
• Store the RMS alongside the elements; it’s a useful quality metric for ranking and filtering.
• If you run estimate_best_orbit repeatedly on the same Observations instance, be aware of the in-place uncertainty scaling caveat described in the Observations tutorial; recreate the object for bitwise reproducibility.