Fitting

Orbit fitting and uncertainty estimation from observations.

This module provides tools for determining and refining orbits from astronomical observations:

Initial Orbit Determination (IOD) – find candidate orbits from a small number of observations using statistical ranging.
Orbit Fitting (fit_orbit()) – refine an orbit guess to best match the data, with automatic outlier rejection. Produces a best-fit state and Gaussian uncertainty estimate (covariance).
MCMC Uncertainty Estimation (fit_orbit_mcmc()) – for short arcs where the Gaussian approximation is unreliable, sample the full range of plausible orbits consistent with the data.
Lambert Solver – compute the transfer orbit between two positions.

class kete.fitting.Observation

Astronomical observation for orbit determination.

Observations can be optical (RA/Dec), radar range, or radar range-rate. Each carries the observer state, measured value(s), and 1-sigma uncertainties.

Use the static methods to construct instances:

Observation.optical()
Observation.radar_range()
Observation.radar_rate()

dec: Declination in degrees (optical only, None otherwise).

epoch: The observation epoch (from the observer state).

observer: The observer state (Sun-centered, Ecliptic).

static optical(observer, ra, dec, sigma_ra, sigma_dec)

Create an optical (RA/Dec) observation.

The observer state is automatically re-centered to the solar system barycenter if needed (the fitting engine works internally in SSB-centered coordinates).

Parameters:

observer (State) – Observer state (any center / frame - will be converted to SSB-centered Equatorial internally). The observation epoch is taken from the observer’s epoch.
ra (float) – Right ascension in degrees.
dec (float) – Declination in degrees.
sigma_ra (float) – 1-sigma RA uncertainty in arcseconds (should include cos(dec) factor).
sigma_dec (float) – 1-sigma Dec uncertainty in arcseconds.

ra: Right ascension in degrees (optical only, None otherwise).

static radar_range(observer, range, sigma_range)

Create a radar range observation.

The observer state is automatically re-centered to SSB.

Parameters:

observer (State) – Observer state (any center / frame – will be converted).
range (float) – Measured range in AU.
sigma_range (float) – 1-sigma range uncertainty in AU.

static radar_rate(observer, range_rate, sigma_range_rate)

Create a radar range-rate (Doppler) observation.

The observer state is automatically re-centered to SSB.

Parameters:

observer (State) – Observer state (any center / frame – will be converted).
range_rate (float) – Measured range-rate in AU/day (positive = receding).
sigma_range_rate (float) – 1-sigma range-rate uncertainty in AU/day.

range: Measured range in AU (radar range only, None otherwise).

range_rate: Measured range-rate in AU/day (radar rate only, None otherwise).

sigma_dec: 1-sigma Dec uncertainty in arcseconds (optical only, None otherwise).

sigma_ra: 1-sigma RA uncertainty in arcseconds (optical only, None otherwise).

sigma_range: 1-sigma range uncertainty in AU (radar range only, None otherwise).

sigma_range_rate: 1-sigma range-rate uncertainty in AU/day (radar rate only, None otherwise).

class kete.fitting.OrbitFit

Result of fitting an orbit to observations.

Returned by fit_orbit(). Contains the best-fit orbital state, its uncertainty (covariance), and diagnostic information about the fit.

all_observations: All input observations (time-sorted), including rejected outliers.

converged

Whether the solver achieved strict convergence.

When False the fit is the best found within the iteration limit but the correction norm did not drop below tol.

included

Per-observation inclusion mask (time-sorted).

True means the observation was used in the final fit; False means it was rejected as an outlier.

non_grav

Fitted non-gravitational model, or None if not fitted.

Convenience shortcut for self.uncertain_state.non_grav.

observations

Observations included in the final fit (time-sorted).

Rejected outliers are not present in this list. To see all observations (including rejected ones), use all_observations.

residuals

Post-fit residuals as a list of lists (included observations only, time-sorted order).

For optical observations the two elements are (DeltaRA, DeltaDec) in arcseconds. Radar residuals remain in AU or AU/day.

rms: Reduced weighted RMS of post-fit residuals.

state

Best-fit state at the reference epoch (Sun-centered, Ecliptic).

Convenience shortcut for self.uncertain_state.state.

uncertain_state: The uncertain orbit state (state + covariance + non-grav model).

class kete.fitting.OrbitSamples

Collection of plausible orbits from MCMC uncertainty estimation.

Returned by fit_orbit_mcmc(). Each orbit (draw) is statistically consistent with the observations; the spread of the collection represents the uncertainty in the orbit.

chain_id: Seed index (0-based) that generated each draw.

desig: Designator of the fitted object.

divergent

Per-draw divergence flag.

A divergent sample indicates the sampler had difficulty exploring that region of orbit space. A small fraction of divergences is normal; many divergences suggest the model or data are problematic.

draws

Sampled orbits as a list of State objects.

Each state is Sun-centered Ecliptic at the reference epoch.

Non-gravitational parameters (if fitted) are available via raw_draws.

epoch: Common reference epoch (JD, TDB).

raw_draws

Raw orbit samples as a list of lists.

Each inner list is [x, y, z, vx, vy, vz, ng_params...] in the Equatorial frame at the reference epoch. Use np.array(samples.raw_draws) to convert to a 2-D array.

class kete.fitting.UncertainState

Uncertain orbit state: a best-fit Cartesian state together with a covariance matrix that may span the 6 position/velocity components and any fitted non-gravitational parameters.

This is the canonical representation of orbit uncertainty in kete, providing correct handling of non-gravitational model templates (preserving fixed physical parameters such as alpha, r_0, etc.).

Construction

from_state() – from a State with isotropic uncertainties.
from_cometary() – from cometary orbital elements and an element-space covariance (e.g. from JPL Horizons).
Returned as part of OrbitFit from orbit fitting.

cov_matrix: Covariance matrix as a list of lists (use np.array() to convert).

desig: Object designator (shortcut for self.state.desig).

epoch: Reference epoch as a Time (shortcut for self.state.epoch).

static from_cometary(elements, cov_matrix, non_grav=None)

Build an UncertainState from cometary orbital elements and a covariance expressed in element space.

The element-space covariance is transformed to a Cartesian covariance via a numerically evaluated Jacobian.

Parameters:

elements (CometElements) – Cometary orbital elements (with desig and epoch).
cov_matrix (list[list[float]]) – Covariance matrix in element space, (6+Np)x(6+Np). Element order: [e, q, tp, node, w, i, <nongrav...>].
non_grav (NonGravModel, optional) – Non-gravitational model template.

static from_state(state, pos_sigma=0.01, vel_sigma=0.0001, non_grav=None)

Build an UncertainState from a state with isotropic diagonal uncertainties.

The covariance is initialised to a diagonal matrix with the given pos_sigma (AU) and vel_sigma (AU/day) on the diagonal. Useful for seeding MCMC from an IOD candidate.

The input state is automatically re-centered to the solar system barycenter if needed.

Parameters:

state (State) – Object state (any center / frame – will be converted to SSB-centered Equatorial internally).
pos_sigma (float) – 1-sigma position uncertainty in AU (default 0.01).
vel_sigma (float) – 1-sigma velocity uncertainty in AU/day (default 0.0001).
non_grav (NonGravModel, optional) – Non-gravitational model template. If provided, the covariance is extended to (6+Np)x(6+Np) with tiny diagonal entries for the non-grav parameters.

non_grav: Non-gravitational model template, or None.

param_names

Names of all parameters in the covariance matrix, in row/column order.

Always starts with ["x", "y", "z", "vx", "vy", "vz"], followed by any non-gravitational parameter names.

sample(n_samples, seed=None)

Draw random samples from the covariance distribution.

Returns a tuple (states, non_gravs) where states is a list of State objects and non_gravs is a list of NonGravModel or None.

Parameters:

n_samples (int) – Number of samples to draw.
seed (int) – Random seed for reproducibility (optional).

state: Best-fit state at the reference epoch (Sun-centered, Ecliptic).

kete.fitting.fit_orbit(initial_state, observations, include_asteroids=False, non_grav=None, max_iter=50, tol=1e-08, chi2_threshold=9.0, max_reject_passes=3, auto_sigma=False)

Fit an orbit to observations using iterative least squares.

Given an initial guess and a set of observations, this function refines the orbital state until it best matches the data, and estimates the uncertainty of the result via a covariance matrix. It can also automatically identify and reject outlier observations.

This is the standard approach for orbit determination when you have a reasonable initial guess (e.g. from initial_orbit_determination()). It works well for arcs of any length, but the Gaussian uncertainty estimate is most reliable for long, well-sampled arcs. For short arcs where the uncertainty is non-Gaussian, consider fit_orbit_mcmc() instead.

For arcs longer than 180 days, progressively wider time windows are fitted around the reference epoch so that each stage bootstraps from the previous converged solution. The final pass fits the full arc and re-evaluates all observations for outlier rejection (if enabled).

The input state is automatically re-centered to the solar system barycenter internally.

Parameters:

initial_state (State) – Initial guess for the object state (any center / frame).
observations (list) – List of Observation to fit.
include_asteroids (bool) – If True, include asteroid masses in the force model (slower but more accurate for near-Earth objects). Default is False.
non_grav (NonGravModel, optional) – Non-gravitational force model.
max_iter (int) – Maximum number of iterations per convergence pass. Default is 50.
tol (float) – Convergence tolerance on the state correction norm. Default is 1e-8.
chi2_threshold (float) – Chi-squared threshold for outlier rejection. Default is 9.0 (roughly 3-sigma per component). Only used when max_reject_passes > 0.
max_reject_passes (int) – Maximum number of outlier-rejection cycles. Set to 0 to disable rejection entirely. Default is 3.
auto_sigma (bool) – If True, adaptively rescale the rejection threshold based on the actual residual scatter rather than the stated uncertainties. Useful when observation uncertainties are unreliable. Default is False.

Returns:

The fitted orbit, including the best-fit state, covariance, residuals, and convergence diagnostics.

Return type:

OrbitFit

kete.fitting.fit_orbit_mcmc(seeds, observations, include_asteroids=False, num_draws=1000, num_tune=500, non_grav=None, maxdepth=10)

Estimate orbit uncertainty from observations using Markov Chain Monte Carlo.

Given one or more candidate orbital states (seeds) and a set of observations, this function produces a collection of plausible orbits that are statistically consistent with the data. The spread of returned orbits represents the uncertainty in the orbit determination – wider spread means less certainty about the true orbit.

This is most useful for short-arc observations (a few nights) where the usual least-squares approach (fit_orbit()) underestimates the true uncertainty. For well-observed objects with long arcs, fit_orbit() alone is usually sufficient and far cheaper.

Under the hood this uses the No-U-Turn Sampler (NUTS), an adaptive variant of Hamiltonian Monte Carlo that efficiently explores the space of possible orbits. Each draw requires a full numerical propagation, so this is orders of magnitude more expensive than least squares – but the result captures non-Gaussian and multi-modal uncertainty that least squares cannot represent.

Seeds are raw State objects (typically from initial_orbit_determination()). No prior orbit fit is required – the sampler builds its own internal mass matrix from a linearization at each seed.

Sampling is parallelized automatically across available CPU cores. When there are fewer seeds than cores, each seed spawns multiple independent sub-chains. The chain_id in the returned OrbitSamples identifies the seed (orbital mode), not the sub-chain.

num_draws is the total number of orbit samples returned across all seeds. Each seed receives roughly num_draws / len(seeds) draws.

Parameters:

seeds (list) – List of State candidate states (e.g. from initial_orbit_determination()), one per orbital mode. Seeds at different epochs are automatically propagated to the first seed’s epoch. The input states are re-centered to SSB Equatorial internally.
observations (list) – List of Observation to evaluate against.
include_asteroids (bool) – If True, include asteroid masses in the force model. Default is False.
num_draws (int) – Total orbit samples across all seeds (after tuning). Default is 1000.
num_tune (int) – Number of warmup steps per sub-chain used to adapt internal sampling parameters. These draws are discarded. Default is 500.
non_grav (NonGravModel, optional) – Shared non-gravitational force model applied to all chains.
maxdepth (int) – Maximum tree depth for the sampler. Higher values allow more thorough exploration at greater computational cost. Default is 10.

Returns:

Collection of plausible orbits sampled from the posterior.

Return type:

OrbitSamples

Raises:

ValueError – If seeds is empty or two-body epoch propagation fails.

kete.fitting.initial_orbit_determination(observations, epoch=None)

Compute an initial orbit from observations.

Parameters:

observations (list[Observation]) – At least 2 optical observations.
epoch (float) – Reference epoch (JD, TDB) for returned states (optional). Defaults to the last observation epoch (for forward prediction).

Returns:

One or more candidate initial states at the reference epoch.

Return type:

list[State]

kete.fitting.lambert(r1, r2, dt, prograde=True)

Solve Lambert’s problem for a single-revolution Keplerian transfer.

Given two heliocentric position vectors and a transfer time, compute the velocity vectors at departure and arrival that connect them via two-body (Keplerian) motion.

Parameters:

r1 (Vector) – Heliocentric position at departure (AU).
r2 (Vector) – Heliocentric position at arrival (AU).
dt (float) – Transfer time in days. Must be positive.
prograde (bool) – If True (default), selects the short-way transfer (transfer angle less than 180 degrees for prograde orbits). If False, selects the long-way transfer.

Returns:

(v1, v2) – velocity at r1 and r2 respectively (AU/day).

Return type:

tuple[Vector, Vector]

Raises:

ValueError – If dt <= 0, positions are zero-length, or positions are nearly collinear (transfer angle near 0 or 180 degrees).
RuntimeError – If the iterative solver fails to converge.

kete.fitting.mpc_obs_to_observations(mpc_obs, sigma_ra=0.1, sigma_dec=0.1)

Convert a list of MPCObservation objects to fitting Observations.

Only optical (RA/Dec) observations are supported. Each MPCObservation is converted to an Observation.optical with the observer state computed from the MPC observatory code (ground-based) or from the stored spacecraft position.

Parameters:

mpc_obs (list) – List of MPCObservation objects (see kete.mpc).
sigma_ra (float) – Default 1-sigma RA uncertainty in arcseconds. The cos(dec) factor is applied automatically.
sigma_dec (float) – Default 1-sigma Dec uncertainty in arcseconds.

Returns:

One Observation.optical per input observation.

Return type:

list[Observation]

Raises:

ValueError – If an observation has a spacecraft flag but no valid position.

Examples

import kete

lines = [...]  # 80-char MPC observation lines
mpc_obs = kete.mpc.MPCObservation.from_lines(lines)
observations = kete.fitting.mpc_obs_to_observations(mpc_obs)
fit = kete.fitting.fit_orbit(initial_state, observations)