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 model for n-body propagation.

The n-body propagation functions accept these models to include forces other than gravity. The constructors are:

NonGravModel.new_dust(): solar radiation pressure and Poynting-Robertson drag on dust.
NonGravModel.new_comet(): the functional form of the JPL Horizons comet model; see that method for the formula.
NonGravModel.new_asteroid(): the comet model with a 1/r^2 falloff, for asteroids with the Yarkovsky effect.
NonGravModel.new_farnocchia(): the Farnocchia et al. radiation model.

The orbit fitting tools read the fittable parameters of each model (for example A1, A2 and A3 of the comet model) with a NaN convention. A NaN value marks the parameter as free, and the fit starts it from 0. A concrete value holds the parameter fixed at that value. NonGravModel.with_free() frees parameters while keeping their values as the starting point of the fit. Propagation treats NaN values as 0.

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=Ellipsis, q_pr=1.0)

Estimate the diameter of the dust particle in meters.

This inverts the beta relation used by NonGravModel.new_dust() and takes the same defaults, so a diameter passed to that constructor is returned unchanged here. Since only beta is stored, the density, c_pr and q_pr used to build the model are not recovered with it and must be supplied again to get back the same diameter.

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

free_parameters: The names of the parameters an orbit fit would fit. These are the NaN values and those freed by with_free().

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, including the NaN convention: pass NaN for any of A1/A2/A3 to leave that parameter free during orbit fitting; concrete values are frozen.

static new_comet(a1=Ellipsis, a2=Ellipsis, a3=Ellipsis, 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.

The A1, A2, and A3 terms follow the NaN convention: a NaN value marks the parameter as free when the model is passed to the orbit fitting tools (the fit starts it at 0), while a concrete value freezes it at that value. Propagation treats NaN as 0.0. The defaults leave all three terms free, so new_comet() requests a fit of A1/A2/A3; pass explicit values to freeze them instead.

static new_dust(beta=None, diameter=None, density=1000.0, c_pr=Ellipsis, 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)\]

NonGravModel.diameter() is the inverse conversion and shares these defaults, so new_dust(diameter=d).diameter() == d.

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. Pass float("nan") to leave beta free during orbit fitting.
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). Pass float("nan") for either to leave it free during orbit fitting; concrete values are frozen. 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). At 0 (zero thermal lag) the transverse Yarkovsky component vanishes but the radial recoil from instantaneous re-emission remains; set absorptivity to 0 to disable the thermal terms entirely.
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, in (0, 1]. Use 1.0 for a sphere.
spin_pole – Spin pole unit vector (any Vector or length-3 sequence). Must be fixed in inertial space.

Raises:

ValueError – If a_over_m or lambda_0 is neither NaN nor a finite value >= 0, if albedo or absorptivity is negative or not finite, if flattening is outside (0, 1], or if spin_pole is zero.

static new_farnocchia_from_h_mag(h_mag, spin_pole, albedo=0.15, density=2500.0, thermal_inertia=200.0, rotation_period=6.0, emissivity=0.9, absorptivity=0.9, flattening=1.0)

Construct a Farnocchia radiation force model from an absolute magnitude H and assumed physical properties.

This is the usual entry point for a collisional family, where H is measured and the rest is assumed. It composes the standard chain:

H and albedo give the diameter, D = 1329 / sqrt(albedo) * 10 ** (-H / 5) km (see compute_diameter()).
the diameter and density give the area-to-mass ratio, A/M = 3 / (4 * density * R), which is the radiation pressure coupling (see a_over_m_from_physical()). Since A/M scales as 1 / (density * D), the drift rate scales as 1 / D: five magnitudes fainter is ten times the drift.
thermal_inertia and rotation_period give the thermal lag (see lambda_0_from_physical()).

The result is an ordinary Farnocchia NonGravModel, usable with both propagate_n_body() and SymplecticSim; the stored a_over_m and lambda_0 are readable so the chain can be checked, and bulk_density() / thermal_inertia() invert it.

Parameters:

h_mag – Absolute magnitude H.
spin_pole – Spin pole (any Vector or length-3 sequence), fixed in inertial space. A collisional family should be given randomly oriented poles rather than one shared pole.
albedo – Geometric albedo. Sets the diameter along with H, and enters radiation pressure.
density – Bulk density in kg / m^3.
thermal_inertia – Thermal inertia in SI units (J m^-2 K^-1 s^-1/2).
rotation_period – Rotation period in hours.
emissivity – Surface emissivity.
absorptivity – alpha = 1 - A_B where A_B is the Bond albedo.
flattening – Axis ratio e = R_P / R_E. Use 1.0 for a sphere.

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.

with_free(*names)

Free parameters for orbit fitting, starting from their current values.

Returns a copy of this model in which the named fittable parameters (every one when none are named) are free, like a NaN, but keep their values. An orbit fit starts those parameters there, and takes initial_state as the matching starting point for the joint fit: the gravity-only pass that otherwise precedes the non-gravitational fit is skipped. This is how a fit is continued from an earlier one:

fit = kete.orbit_fitting.fit_orbit(fit.state, obs, non_grav=fit.non_grav.with_free())

Propagation uses the values as given. Parameters that are NaN stay free and start from 0.

Parameters:: names – Names of fittable parameters, as in free_parameters of a model with every parameter NaN (for example "a1").

class kete.propagation.SymplecticSim(massive_states, masses, test_particles=Ellipsis, dt=4.0, include_gr=True, include_j2=True, use_correctors=True, non_gravs=None)

A long-term symplectic N-body simulation of the solar system.

This is a fixed-step Wisdom-Holman integrator (democratic heliocentric splitting) intended for million-to-billion year orbital evolution of the planets, massive asteroids, and massless test particles. Unlike propagate_n_body(), which follows the SPICE kernels and is limited to their time span, this evolves the entire system self-consistently and has no time limit.

Because the step size is fixed, it should be at most 1/20th of the shortest orbital period in the simulation (about 4 days when Mercury is included). Close encounters are not resolved by this method; the simulation records the closest approach seen so results can be audited (see closest_encounter).

Test particles which fall into the Sun are removed and recorded in lost_particles.

Test particles may optionally carry a NonGravModel, the same per-object non-gravitational description accepted by propagate_n_body(). A Farnocchia model gives the particle radiation pressure and thermal recoil (Yarkovsky), whose along-track component drifts the semi-major axis at a rate set by the angle between the spin pole and the orbit normal; over millions of years this is what spreads a collisional family, so a population should be given randomly oriented poles (see new_farnocchia_from_h_mag()). A dust model makes the particle a dust grain: it moves on the radiation-reduced two-body orbit (1 - beta) * GM_sun and is damped by Poynting-Robertson drag, so its orbit decays secularly. A JPL-style A1/A2/A3 model applies the fitted radial/transverse/normal accelerations with their g(r) scaling – the form JPL orbit solutions use both for comet outgassing and for asteroid Yarkovsky detections via A2, so covariance samples from sample() feed directly into non_gravs. The fitted constants are extrapolated unchanged over the whole integration, and the time-lagged (dt != 0) outgassing variant is rejected.

Total energy of the massive system oscillates within a bounded band rather than trending, which is the usual quality diagnostic for a run; see energy.

This is a reduced dynamical model built for long-term statistical dynamics, not short-term ephemeris prediction: point-mass Newtonian gravity plus the optional GR and J2 terms, a limited asteroid list, the Earth-Moon pair merged into a barycenter carrying its orbit-averaged lunar quadrupole, and no Pluto unless registered. Measured against JPL DE441 over 15 kyr backwards, planet positions drift away roughly linearly in time, dominated by the remaining model truncations rather than by integration error: every planet stays within a few 1e-4 AU over the full 15 kyr except Neptune, which reaches ~1e-2 AU unless the Pluto system is included (kete.register_mass(9) before from_spice() reduces Neptune’s drift by roughly an order of magnitude). For ephemeris-accurate short-term work use propagate_n_body().

Parameters:

massive_states – Barycentric states of all massive bodies. The first entry must be the Sun. For long integrations use the Earth-Moon barycenter (NAIF id 3) rather than the separate Earth and Moon.
masses – Mass of each massive body as a fraction of the Sun’s mass, in the same order as massive_states (the Sun itself is 1.0).
test_particles – States of massless test particles. They feel every massive body but affect nothing.
dt – Fixed step size in days. Negative integrates backwards in time.
include_gr – Apply the general relativity correction of the Sun: the first-order Schwarzschild acceleration shared with propagate_n_body(), which reproduces both the secular apsidal precession (43 arcsec/century for Mercury) and the relativistic mean motion. The term is velocity dependent, so with it enabled the integrator is only approximately symplectic; the energy error remains a bounded band with no secular trend.
include_j2 – Apply the solar J2 oblateness term, with the same coefficient and solar spin pole (the IAU pole, as in DE440) used by propagate_n_body(). This is a position-only potential, so the integrator remains symplectic; it drives the small secular nodal regression and apsidal precession of low semi-major axis orbits.
use_correctors – Wrap each integration call in the order-17 symplectic corrector of Wisdom (2006), which removes the dominant oscillating error of the map (roughly a factor of 100 in the energy band). The overhead is a fixed ~30 steps worth of work per integrate_n_steps() or integrate_to() call, so batch many steps per call rather than looping over step().
non_gravs – Optional list of NonGravModel, one per test particle, with None entries for particles which feel gravity alone. A NaN parameter, left free for orbit fitting, is treated as 0, as in every other propagation.

angular_momentum: Total angular momentum of the massive system multiplied by the gravitational constant, in AU^5 / Day^3, on equatorial axes.

closest_encounter

The closest approach between two bodies seen so far, as a (designation, designation, hill_ratio, jd) tuple, or None.

Only approaches within 3 Hill radii of a massive body are recorded. If this is not None, the trajectories of the involved bodies should be treated with suspicion, since the fixed-step map does not resolve close encounters.

dt: Fixed step size in days.

energy

Total energy of the massive system multiplied by the gravitational constant, in AU^5 / Day^4.

For a symplectic map this oscillates within a narrow bounded band; its relative change is the standard quality diagnostic for a run.

static from_spice(jd, dt=4.0, include_registered=True, test_particles=Ellipsis, include_gr=True, include_j2=True, use_correctors=True, non_gravs=None)

Build a simulation from the loaded SPICE kernels.

The Sun and the 8 planets (Earth and Moon merged into their barycenter) are always included. The currently registered massive bodies (by default Ceres, Pallas, Vesta, Hygiea, and Interamnia; see register_mass() and register_custom_mass() to add more, for example the Pluto system with kete.register_mass(9)) may optionally be added as fully massive bodies.

Parameters:

jd – Epoch at which to pull the initial states from SPICE. Must be within the loaded kernels’ time span.
dt – Fixed step size in days. Negative integrates backwards in time.
include_registered – Include the registered massive bodies beyond the planets.
test_particles – States of massless test particles.
include_gr – Apply the general relativity correction of the Sun.
include_j2 – Apply the solar J2 oblateness term; see the class documentation.
use_correctors – Wrap each integration call in the order-17 symplectic corrector, see the class documentation.
non_gravs – Optional list of NonGravModel, one per test particle, with None entries for particles which feel gravity alone; see the class documentation.

include_gr: Whether the GR correction is enabled.

include_j2: Whether the solar J2 oblateness term is enabled.

integrate_n_steps(n)

Advance the simulation by the given number of steps.

The GIL is released for the duration of the integration. A keyboard interrupt stops it within a few thousand steps, at the end of a step, with the simulation usable from there.

integrate_to(jd)

Advance the simulation to approximately the target time.

The step size is fixed, so the integration lands on the whole step closest to the requested time and no partial step is taken; check jd for the resulting epoch. The GIL is released for the duration of the integration. A keyboard interrupt stops it within a few thousand steps, at the end of a step, with the simulation usable from there.

Parameters:: jd – Target time. Must not be behind the current epoch with respect to the sign of dt.

jd: Current epoch of the simulation as a TDB scaled JD.

lost_particles: Test particles removed from the simulation, as a list of (designation, jd, reason) tuples. reason is either "sun_impact" or "kepler_failure".

massive_states: Current barycentric states of the Sun and all massive bodies, in the order they were provided (the Sun first).

n_test_particles: Number of currently active test particles.

shortest_period: The shortest osculating orbital period among the massive bodies, in days. dt should be at most 1/20th of this.

step()

Advance the simulation by a single step of dt.

With correctors enabled every call pays the fixed corrector overhead; prefer integrate_n_steps() or integrate_to() for more than a handful of steps.

steps_taken: Number of steps taken so far.

test_particle_states: Current barycentric states of the active test particles, including any carrying non-gravitational forces.

use_correctors: Whether the order-17 symplectic corrector is enabled.

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 Minimum Orbital Intersection Distance (MOID).

The MOID is the smallest distance between a point on the orbit of state_a and a point on the orbit of state_b. Each state defines a two-body orbit about the Sun. A state with a center other than the Sun is converted to a Sun centered state with the loaded SPICE kernels.

Parameters:

state_a (State or list of State) – State of the first object, or a list of states.
state_b (State, optional) – State of the second object. The default is the Earth, at the epoch of the first state in state_a.

Returns:

MOID in au, one value for each state in state_a. The value is NaN for a state whose MOID cannot be computed, such as a state with purely radial motion.

Return type:

float or list of float

Raises:

ValueError – If state_a is empty, or if the Earth state or state_b cannot be converted to a Sun centered state.

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, Ceres, Pallas, Interamnia, Hygiea, and Vesta, along with any masses registered in addition.
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 are not printed to stderr. With suppress_errors the state of an object that impacted a massive body is returned filled with NaN either way.

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 two-body motion about the Sun. This is a multi-core operation.

Parameters:

states (State or list of State) – States to propagate, in AU and AU/Day.
epoch (float or Time) – Time to propagate to, in JD with TDB scaling.
observer_pos (Vector, optional) – Sun-centered position of an observer. If it is given, the states are corrected for light travel time to the observer. The delay is iterated until it changes by less than 1e-12 days, at most 3 times. Defaults to no correction.

Returns:

States after propagation. A state that fails to propagate has NaN values.

Return type:

State or list of State

kete.propagation.register_custom_mass(naif_id, mass, radius=0.0)

Register a new massive object to be used in the extended list of objects.

Masses must be provided as a fraction of the Sun’s mass, and radius in AU.

An object already registered with the same NAIF ID is replaced.

kete.propagation.register_mass(naif_id)

Register a massive object from the known-masses table by its NAIF ID.

Use [register_custom_mass] to add a body that is not in the table.

# Errors

Returns an error if the NAIF ID is not present in the built-in mass table.

kete.propagation.register_polyhedron(naif_id, vertices, faces, switch_radius, mass=None, frame_id=None, rotation=None)

Register a massive body whose gravity near it is a constant-density polyhedron.

Registered bodies are part of the force model when include_asteroids=True is passed to propagate_n_body() or orbit fitting; the default force model is a fixed list of planets and the Moon, which a registration does not change. Within switch_radius of the body those functions use the polyhedron field of the shape model in place of the point mass; beyond it, the point mass. The body must have an SPK ephemeris loaded; the shape model’s origin is placed at the body’s SPK position, so vertices are given in the body frame about that point. The point mass beyond switch_radius is also at the SPK position, while the polyhedron’s center of mass is its volume centroid; if the two differ the force changes across the switch by that offset as well. A body already registered with the same NAIF ID is replaced.

Exactly one orientation is required: frame_id, a body frame read at each evaluation inside switch_radius from the loaded PCK files or, when they have no frame with that id, the CK files with their clock, where a time without orientation is an error, or rotation, a fixed matrix taking body-frame vectors to equatorial (J2000) axes.

The polyhedron field and the point mass agree ever more closely with distance, so switch_radius trades evaluation cost against the step in the force where they meet; that step is largest for elongated bodies and falls as the inverse square of switch_radius.

Parameters:

naif_id – NAIF ID of the body.
vertices – Shape model vertices in the body frame, in km, shape (n, 3).
faces – Faces as 0-based vertex indices, counter-clockwise seen from outside, shape (m, 3).
switch_radius – Distance from the body in km inside which the polyhedron field is used.
mass – Mass of the body as a fraction of the Sun’s mass. Defaults to the value in the built-in mass table.
frame_id – SPICE frame ID or frame name of the body frame.
rotation – Fixed rotation matrix from the body frame to equatorial axes, shape (3, 3).

Returns:

The registered field, in AU and AU^3/day^2.

Return type:

Polyhedron

kete.propagation.register_spherical_harmonics(naif_id, c, s, radius, min_radius, switch_radius, mass=None, frame_id=None, rotation=None)

Register a massive body whose gravity near it is a spherical harmonic field.

Registered bodies are part of the force model when include_asteroids=True is passed to propagate_n_body() or orbit fitting; the default force model is a fixed list of planets and the Moon, which a registration does not change. Within switch_radius of the body those functions use the series in place of the point mass; beyond it, the point mass. The body must have an SPK ephemeris loaded; the expansion origin is the body’s SPK position. Inside min_radius (the body’s Brillouin radius, or larger) the series is not valid, and a propagation that reaches it raises as an impact with the body. A body already registered with the same NAIF ID is replaced.

The coefficients are fully normalized, without the Condon-Shortley phase (see SphericalHarmonics); check the convention of a published field before using it. C_00 must be 1, since the series’ mass is the body’s, which is also the point mass beyond switch_radius.

Exactly one orientation is required: frame_id, a body frame read at each evaluation inside switch_radius from the loaded PCK files or, when they have no frame with that id, the CK files with their clock, where a time without orientation is an error, or rotation, a fixed matrix taking body-frame vectors to equatorial (J2000) axes.

Parameters:

naif_id – NAIF ID of the body.
c – Normalized C_nm, one row per degree n with n + 1 entries.
s – Normalized S_nm, the same shape; S_n0 must be 0.
radius – Reference radius of the coefficients, in km.
min_radius – Distance from the body in km inside which the series is not valid.
switch_radius – Distance from the body in km inside which the series is used; at least min_radius.
mass – Mass of the body as a fraction of the Sun’s mass. Defaults to the value in the built-in mass table.
frame_id – SPICE frame ID or frame name of the body frame.
rotation – Fixed rotation matrix from the body frame to equatorial axes, shape (3, 3).

Returns:

The registered field, in AU and AU^3/day^2.

Return type:

SphericalHarmonics

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.