Propagation

Propagation of objects using orbital mechanics, this includes a simplified 2 body model as well as a N body model which includes some general relativistic effects.

class kete.propagation.NonGravModel

Non-gravitational force models.

This is used optionally by the N-Body propagation methods to compute orbits including non-gravitational forces, such as solar radiation pressure, or poynting-robertson force.

There are two generic non-gravitational models available, one is specifically intended for dust modeling, and includes the solar radiation pressure, the other implements the functional form documented for the JPL Horizons comet model (see NonGravModel.new_comet() for the formula).

See NonGravModel.new_dust() and NonGravModel.new_comet() for more details. Note that the Comet model can also represent asteroids which undergo the Yarkovsky effect, see NonGravModel.new_asteroid(), which is a convenience function over the NonGravModel.new_comet() method, but with 1/r^2 falloff.

a_over_m

Stored area-to-mass ratio A/M (m^2 / kg) for a FarnocchiaModel (Eq. 6).

Returns NaN unless this is a FarnocchiaModel.

beta: Get the beta value for this dust model.

bulk_density(diameter)

Recover bulk density (kg / m^3) of this FarnocchiaModel given the auxiliary inputs that were collapsed away at construction time.

Returns NaN unless this is a FarnocchiaModel.

diameter(density=1000.0, c_pr=1.19, q_pr=1.0)

Estimate the diameter of the dust particle in meters.

Only works for dust models, returns NaN for asteroid/comet models.

Parameters:

density – Density in kg/m^3, defaults to 1000 kg/m^3
c_pr – Radiation pressure coefficient, defaults to 1.19 kg/m^2
q_pr – Scattering efficiency for radiation pressure, defaults to 1.0 1.0 is a good estimate for particles larger than 1um (Burns, Lamy & Soter 1979)

items: Return a dictionary of the values used in this non-grav model.

lambda_0

Stored thermal parameter lambda_0 (dimensionless, Eq. 12) for a FarnocchiaModel.

Returns NaN unless this is a FarnocchiaModel.

static new_asteroid(a1, a2, a3, alpha=1.0, r_0=1.0, m=2.0, n=1.0, k=0.0, dt=0.0)

This is the same as NonGravModel.new_comet(), but with default values set so that \(g(r) = 1/r^2\).

See NonGravModel.new_comet() for more details.

Setting an A term to float('nan') excludes it from fitting and treats it as zero in the force model.

static new_comet(a1=0.0, a2=0.0, a3=0.0, alpha=0.1112620426, r_0=2.808, m=2.15, n=5.093, k=4.6142, dt=0.0)

JPL’s non-gravitational forces are modeled as defined on page 139 of the Comets II textbook.

This model adds 3 “A” terms to the acceleration which the object feels. These A terms represent additional radial, tangential, and normal forces on the object.

The defaults of this method are the defaults that JPL Horizons uses for comets when they are not otherwise specified.

\[\text{accel} = A_1 g(r) \vec{r} + A_2 g(r) \vec{t} + A_3 g(r) \vec{n}\]

Where \(\vec{r}\), \(\vec{t}\), \(\vec{n}\) are the radial, tangential, and normal unit vectors for the object.

The \(g(r)\) function is defined by the equation:

\[g(r) = \alpha \big(\frac{r}{r_0}\big) ^ {-m} \bigg(1 + \big(\frac{r}{r_0}\big) ^ n\bigg) ^ {-k}\]

When alpha=1.0, n=0.0, k=0.0, r0=1.0, and m=2.0, this is equivalent to a \(1/r^2\) correction.

This includes an optional time delay, which the non-gravitational forces are time delayed.

Setting an A term to float('nan') excludes it from fitting and treats it as zero in the force model. This allows fitting any subset of A1, A2, A3. For example, to fit only A2:

NonGravModel.new_comet(a1=float('nan'), a2=0.0, a3=float('nan'))

static new_dust(beta=None, diameter=None, density=1000.0, c_pr=0.00119, q_pr=1.0)

Create a new non-gravitational forces Dust model.

This implements the radiative force model presented in: “Radiation forces on small particles in the solar system” Icarus, Vol 40, Issue 1, Pages 1-48, 1979 Oct https://doi.org/10.1016/0019-1035(79)90050-2

The model calculated has the acceleration of the form:

\[\text{accel} = \frac{L_0 A Q_{pr}}{r^2 c m} \bigg((1 - \frac{\dot{r}}{c}) \vec{S} - \vec{v} / c \bigg)\]

Where \(L_0\) is the luminosity of the Sun, A is the effective cross sectional area of the dust, \(Q_{pr}\) is a scattering coefficient (~1 for dust larger than about 0.1 micron), m mass, c speed of light, and r heliocentric distance.

The vectors on the right are \(\vec{S}\) the position with respect to the Sun. \(\vec{v}\) the velocity with respect to the Sun. \(\dot{r}\) is the radial velocity toward the sun.

This equation includes both the effects from solar radiation pressure in addition to the Poynting-Robertson effect. By neglecting the Poynting-Robertson components of the above formula, it is possible to find a mapping from the standard \(\beta\) formalism to the above coefficient:

\[\beta = \frac{L_0 A Q_{pr}}{c m G}\]

Where G is the solar standard gravitational parameter (GM). Making the above equation equivalent to:

\[\text{accel} = \frac{\beta G}{r^2} \bigg((1 - \frac{\dot{r}}{c}) \vec{S} - \vec{v} / c \bigg)\]

Parameters:

beta – Beta value of the dust, if this is specified, all other inputs are ignored. If this value is specified, diameter cannot be specified.
diameter – Diameter of the dust particle in meters, this uses the following parameters to estimate the beta value. If beta is specified, this cannot be specified.
density – Density in kg/m^3, defaults to 1000 kg/m^3
c_pr – Radiation pressure coefficient, defaults to 1.19e-3 kg/m^2
q_pr – Scattering efficiency for radiation pressure, defaults to 1.0 1.0 is a good estimate for particles larger than 1um (Burns, Lamy & Soter 1979)

static new_farnocchia(a_over_m, lambda_0, albedo, absorptivity, flattening, spin_pole)

Construct a physical radiation force model from Farnocchia et al. 2025.

Models the body as an oblate spheroid with a fixed spin pole and computes solar radiation pressure plus thermal recoil (Yarkovsky) acceleration.

The two fittable parameters are taken in the form used by the paper: a_over_m (Eq. 6) and lambda_0 (Eq. 12). Use the helpers kete.propagation.a_over_m_from_physical() and kete.propagation.lambda_0_from_physical() to compute them from physical surface inputs (density, thermal inertia, diameter, rotation period, etc.).

Parameters:

a_over_m – Area-to-mass ratio in m^2 / kg (Eq. 6: A/M = 3 / (4 * rho * R_P)).
lambda_0 – Dimensionless thermal lag parameter at 1 AU (Eq. 12). Set to 0 to disable the thermal (Yarkovsky) component.
albedo – Geometric (Lambert) albedo, enters SRP only.
absorptivity – alpha = 1 - A_B where A_B is the Bond albedo. Multiplies the thermal terms.
flattening – Axis ratio e = R_P / R_E. Use 1.0 for a sphere.
spin_pole – Spin pole unit vector (any Vector or length-3 sequence). Must be fixed in inertial space.

thermal_inertia(emissivity, rotation_period)

Recover surface thermal inertia Gamma (SI units) of this FarnocchiaModel given the auxiliary inputs that were collapsed away at construction time. rotation_period is in hours.

Returns NaN unless this is a FarnocchiaModel.

kete.propagation.a_over_m_from_physical(density, diameter, flattening)

Compute A/M (m^2 / kg) from physical surface inputs (Farnocchia 2025 Eq. 6).

Parameters:

density – Bulk density in kg / m^3.
diameter – Volume-equivalent diameter in km.
flattening – Axis ratio e = R_P / R_E. Use 1.0 for a sphere.

kete.propagation.closest_approach(state_a, state_b, jd_start, jd_end, include_asteroids=False)

Find the epoch and distance of closest approach between two objects.

Both objects are propagated using full N-body mechanics over the search window. A coarse grid scan followed by golden-section refinement locates the minimum separation.

Parameters:

state_a – Orbital state of the first object.
state_b – Orbital state of the second object.
jd_start – Start of the search window (JD, TDB).
jd_end – End of the search window (JD, TDB).
include_asteroids – Whether to include asteroid masses in the force model.

Returns:

(epoch, distance) where epoch is a float JD (TDB) and distance is in AU.

Return type:

tuple

kete.propagation.density_from_a_over_m(a_over_m, diameter, flattening): Inverse of a_over_m_from_physical(): solve for bulk density (kg / m^3) given A/M, diameter (km), and flattening.

kete.propagation.lambda_0_from_physical(thermal_inertia, emissivity, absorptivity, flattening, rotation_period)

Compute lambda_0 (dimensionless, Farnocchia 2025 Eq. 12) from physical surface inputs.

Parameters:

thermal_inertia – Surface thermal inertia Gamma in SI units (J m^-2 s^-1/2 K^-1).
emissivity – Thermal emissivity.
absorptivity – alpha = 1 - A_B where A_B is the Bond albedo.
flattening – Axis ratio e = R_P / R_E.
rotation_period – Rotation period in hours.

kete.propagation.moid(state_a, state_b=None)

Compute the MOID between the input state and an optional second state. If the second state is not provided, default to Earth.

Returns the MOID in units of au.

Parameters:

state_a – State of the first object.
state_b – Optional state of the second object, defaults to Earth.

kete.propagation.propagate_n_body(states, jd, include_asteroids=False, non_gravs=None, suppress_errors=True, suppress_impact_errors=False)

Propagate the provided State using N body mechanics to the specified times, no approximations are made, this can be very CPU intensive.

This does not compute light delay, however it does include corrections for general relativity due to the Sun.

Parameters:

states – The initial states, this is a list of multiple State objects.
jd – A JD to propagate the initial states to.
include_asteroids – If this is true, the computation will include the largest 5 asteroids. The asteroids are: Ceres, Pallas, Interamnia, Hygiea, and Vesta.
non_gravs – A list of non-gravitational terms for each object. If provided, then every object must have an associated NonGravModel or None.
suppress_errors – If True, errors during propagation will return NaN for the relevant state vectors, but propagation will continue.
suppress_impact_errors – If True, impacts will be printed to stderr, but states will still return filled with NaN. If False, impacts are not printed.

Returns:

A State at the new time.

Return type:

Iterable

kete.propagation.propagate_n_body_long(states, jd_final, planet_states=None, non_gravs=None, batch_size=10)

It is STRONGLY recommended to use propagate_n_body instead of this function wherever possible. This function is specifically meant for kilo-year or longer simulations, it is slower and less accurate than propagate_n_body, but that function only works for as long as there are SPICE kernels for the planets available.

This is designed to not require SPICE kernels to function, and is meant for long term simulations on the other of less than a mega-year. It is not recommended for orbits longer than this.

Propagation using this will treat the Earth and Moon as a single object for performance reasons.

Propagate the provided State using N body mechanics to the specified times, very few approximations are made, this can be very CPU intensive.

This does not compute light delay, however it does include corrections for general relativity due to the Sun.

This returns two lists of states: - First one contains the states of the objects at the end of the integration - Second contains the states of the planets at the end of the integration.

The second set of states may be used as input for continuing the integration.

Parameters:

states – The initial states, this is a list of multiple State objects.
jd – A JD to propagate the initial states to.
planet_states – Optional list of the planet’s states at the same time as the provided states. If this is not provided, the planets positions are loaded from the Spice kernels if that information is available.
non_gravs – A list of non-gravitational terms for each object. If provided, then every object must have an associated NonGravModel.
batch_size – Number of objects to propagate at once with the planets. This is used to break up the simulation for multi-core support. It additionally has effects on the integrator stepsize which is difficult to predict before running. This can be manually tuned for increased performance, it should have no other effects than performance.

Returns:

A State at the new time.

Return type:

Iterable

kete.propagation.propagate_two_body(states, epoch, observer_pos=None)

Propagate the State for all the objects to the specified time. This assumes 2 body interactions.

This is a multi-core operation.

Parameters:

state – List of states, which are in units of AU from the Sun and velocity is in AU/Day.
epoch – Time to integrate to in JD days with TDB scaling.
observer_pos – A vector of length 3 describing the position of an observer. If this is provided then the estimated states will be returned as a result of light propagation delay.

Returns:

Final states after propagating to the target time.

Return type:

State

kete.propagation.thermal_inertia_from_lambda_0(lambda_0, emissivity, absorptivity, flattening, rotation_period): Inverse of lambda_0_from_physical(): solve for thermal inertia (SI units) given lambda_0 and the auxiliary surface inputs. rotation_period is in hours.