Skip to content

Quickstart

This quickstart shows the smallest end-to-end example: create an environment, define (or fetch) an observer, ingest three astrometric observations for a single trajectory, and run a Gauss Initial Orbit Determination (IOD) to obtain a preliminary orbit.

If you have not installed the package yet, see the Installation page first.

1. Import and environment

The environment holds ephemerides and the astrometric error model. The first call may trigger a download of planetary ephemerides (JPL DE440) depending on cache state.

Environment initialization
1
2
3
from py_outfit import PyOutfit

env = PyOutfit("horizon:DE440", "FCCT14")  # (ephemerides_id, error_model)

2. Define or fetch an observer

You can supply your own geodetic coordinates or use an MPC code if available. Elevation is expected in kilometers.

Observer creation
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
from py_outfit import Observer

observer = Observer(
    longitude=12.345,  # degrees East
    latitude=-5.0,  # degrees North
    elevation=1.5,  # kilometers
    name="DemoSite",
    ra_accuracy=None,
    dec_accuracy=None,
)
env.add_observer(observer)

# Alternatively using MPC code:
# observer = env.get_observer_from_mpc_code("I41")

3. Prepare a minimal observation set

We create three observations (the minimum for Gauss IOD) at distinct times. Use either the degree ingestion path (shown) or radians. Times are MJD(TT). Uncertainties are 1-sigma values (arcseconds in the degree path).

Setup small data
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
import numpy as np

trajectory_id = np.array([0, 0, 0], dtype=np.uint32)  # single trajectory with ID 0
ra_deg = np.array(
    [
        20.9191548,
        19.8927380,
        18.2218784,
    ]
)  # Right Ascension in degrees
dec_deg = np.array(
    [
        20.0550441,
        20.2977473,
        20.7096409,
    ]
)  # Declination in degrees
mjd_tt = np.array(
    [
        58789.13709704,
        58793.19047664,
        58801.28714334,
    ]
)  # Observation epochs (TT)

err_ra_arcsec = 0.5  # uniform RA uncertainty
err_dec_arcsec = 0.5  # uniform Dec uncertainty

4. Build a TrajectorySet and extract the single trajectory

The TrajectorySet groups observations by ID. Even for one object it is the simplest way to obtain an Observations handle that exposes estimate_best_orbit.

Push data into a trajectory container
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
from py_outfit import TrajectorySet

traj_set = TrajectorySet.from_numpy_degrees(
    env,
    trajectory_id,
    ra_deg.astype(float),
    dec_deg.astype(float),
    float(err_ra_arcsec),
    float(err_dec_arcsec),
    mjd_tt.astype(float),
    observer,
)

obs = traj_set[0]  # Observations wrapper for trajectory ID 0
print("Number of observations:", len(obs))

5. Configure IOD parameters

Use the builder to tweak only what you need. Here we disable noise realizations for determinism and cap the search space for speed.

IOD Parameter configuration
1
2
3
4
5
6
7
8
from py_outfit import IODParams

params = (
    IODParams.builder()
    .n_noise_realizations(0)  # purely geometric Gauss solve
    .max_triplets(50)  # limit combinatorial expansion
    .build()
)

6. Run Gauss IOD for the single trajectory

estimate_best_orbit returns a (GaussResult, rms) pair. The RMS is an internal quality metric (angular residual scale).

Initial orbit determination
1
2
3
gauss_result, rms = obs.estimate_best_orbit(env, params, seed=42)
print("RMS:", rms)
print("Elements family:", gauss_result.elements_type())

7. Inspect the resulting orbital elements

Depending on internal selection, a Keplerian, Equinoctial, or Cometary set is produced. Access helpers return None when the family does not match.

inspect results
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
kep = gauss_result.keplerian()
if kep is not None:
    print("Semi-major axis (AU):", kep.semi_major_axis)
    print("Eccentricity:", kep.eccentricity)
    print("Inclination (rad):", kep.inclination)
else:
    eq = gauss_result.equinoctial()
    if eq is not None:
        print("Equinoctial h,k:", eq.h, eq.k)
    else:
        com = gauss_result.cometary()
        print("Cometary perihelion distance (AU):", com.perihelion_distance)

8. Full minimal script (copy & run)

Full quickstart script (copy and past)
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
from py_outfit import PyOutfit

env = PyOutfit("horizon:DE440", "FCCT14")  # (ephemerides_id, error_model)

from py_outfit import Observer

observer = Observer(
    longitude=12.345,  # degrees East
    latitude=-5.0,  # degrees North
    elevation=1.5,  # kilometers
    name="DemoSite",
    ra_accuracy=None,
    dec_accuracy=None,
)
env.add_observer(observer)

# Alternatively using MPC code:
# observer = env.get_observer_from_mpc_code("I41")

import numpy as np

trajectory_id = np.array([0, 0, 0], dtype=np.uint32)  # single trajectory with ID 0
ra_deg = np.array(
    [
        20.9191548,
        19.8927380,
        18.2218784,
    ]
)  # Right Ascension in degrees
dec_deg = np.array(
    [
        20.0550441,
        20.2977473,
        20.7096409,
    ]
)  # Declination in degrees
mjd_tt = np.array(
    [
        58789.13709704,
        58793.19047664,
        58801.28714334,
    ]
)  # Observation epochs (TT)

err_ra_arcsec = 0.5  # uniform RA uncertainty
err_dec_arcsec = 0.5  # uniform Dec uncertainty

from py_outfit import TrajectorySet

traj_set = TrajectorySet.from_numpy_degrees(
    env,
    trajectory_id,
    ra_deg.astype(float),
    dec_deg.astype(float),
    float(err_ra_arcsec),
    float(err_dec_arcsec),
    mjd_tt.astype(float),
    observer,
)

obs = traj_set[0]  # Observations wrapper for trajectory ID 0
print("Number of observations:", len(obs))

from py_outfit import IODParams

params = (
    IODParams.builder()
    .n_noise_realizations(0)  # purely geometric Gauss solve
    .max_triplets(50)  # limit combinatorial expansion
    .build()
)

gauss_result, rms = obs.estimate_best_orbit(env, params, seed=42)
print("RMS:", rms)
print("Elements family:", gauss_result.elements_type())

kep = gauss_result.keplerian()
if kep is not None:
    print("Semi-major axis (AU):", kep.semi_major_axis)
    print("Eccentricity:", kep.eccentricity)
    print("Inclination (rad):", kep.inclination)
else:
    eq = gauss_result.equinoctial()
    if eq is not None:
        print("Equinoctial h,k:", eq.h, eq.k)
    else:
        com = gauss_result.cometary()
        print("Cometary perihelion distance (AU):", com.perihelion_distance)

At the end of the run, typical console output looks like: 

Number of observations: 3
RMS: 0.62            # angular residual scale (approx.)
Elements family: keplerian
Semi-major axis (AU): 2.72084815
Eccentricity: 0.27511014
Inclination (rad): 0.27433785 (~15.7 deg)

Note: Exact numbers can vary slightly (last decimals) depending on platform, floating‑point rounding, and random seed (if noise realizations were enabled). Differences at the 1e-9–1e-12 level are normal.

9. Next steps

Proceed to:

  • Tutorials for batch processing and parameter tuning.
  • API reference (IODParams, GaussResult, orbital element classes) for deeper control.

If your first run downloads ephemerides, subsequent runs should start much faster.