harv.kepler package

Keplerian orbit mechanics and related utilities.

class harv.kepler.KeplerianBody

Bases: Module

Orbital parameters of a Keplerian body (companion, i.e. body 2).

This class represents the orbital parameters of a companion or second body (relative to some barycenter). So, all parameters represent the orbital elements of a specific body relative to the barycenter.

The primary parameterization uses: - Period (P) - Eccentricity (e) - Semi-major axis (a) or semi_major_axis - Time of pericenter or t_peri

Alternative constructors support different parameterizations.

Parameters:

  • period (Real[Quantity[PhysicalType('time')], '']) – Orbital period.

  • eccentricity (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '']) – Orbital eccentricity.

  • semi_major_axis (Real[Quantity[PhysicalType('length')], '']) – The semi-major axis of the body relative to its system barycenter.

  • t_peri (Real[Quantity[PhysicalType('time')], '']) – Time of pericenter passage.

  • orientation (KeplerianOrientation) – Optional: Orientation of the orbit.

  • ecc_zero_tol (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '']) – Optional: Tolerance for treating eccentricity as zero in circular orbit shortcuts. This is used to avoid numerical issues in the Kepler solver when eccentricity is very small, and can be set to a small multiple of machine epsilon for the data type used.

Examples

>>> from unxt import Q
>>> from harv.kepler.body import KeplerianBody
>>> body = KeplerianBody(
...     period=Q(365.25, "day"),
...     eccentricity=0.1,
...     semi_major_axis=Q(1.0, "AU"),
...     t_peri=Q(0.0, "day"),
... )
__init__(period, eccentricity, semi_major_axis, t_peri, orientation=KeplerianOrientation(), *, ecc_zero_tol=np.float64(2.220446049250313e-15))
Parameters:

  • period (Real[Quantity[PhysicalType('time')], ''])

  • eccentricity (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • semi_major_axis (Real[Quantity[PhysicalType('length')], ''])

  • t_peri (Real[Quantity[PhysicalType('time')], ''])

  • orientation (KeplerianOrientation)

  • ecc_zero_tol (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

Return type:

None

static __new__(cls, *args, **kwargs)
Parameters:
Return type:

TypeVar(_ModuleT, bound= Module)

ecc_zero_tol: Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''] = np.float64(2.220446049250313e-15)
classmethod from_masses(period, eccentricity, m_total, m_body, t_peri, **kwargs)

Construct body’s barycentric orbit from masses and period.

Computes this body’s barycentric semi-major axis from Kepler’s 3rd law: 1. Compute relative orbit: a_rel = (G m_total P^2 / 4 pi^2)^(1/3) 2. Convert to barycentric: a_body = a_rel * (1 - m_body / m_total)

Parameters:

  • period (Real[Quantity[PhysicalType('time')], '']) – Orbital period.

  • eccentricity (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '']) – Orbital eccentricity.

  • m_total (Real[Quantity[PhysicalType('mass')], '']) – Total system mass.

  • m_body (Real[Quantity[PhysicalType('mass')], '']) – Mass of this body.

  • t_peri (Real[Quantity[PhysicalType('time')], '']) – Time of pericenter passage.

  • orientation – Optional: Orientation of the orbit.

  • kwargs (Any) – Additional keyword arguments to pass to the main constructor.

Return type:

KeplerianBody

Returns:

The body’s orbit about the system barycenter.

get_mass(m_total)

Recover this body’s mass from the total system mass.

This inverts from_masses: given the total system mass and the stored barycentric semi-major axis, it recovers the body mass via

\[ \begin{align}\begin{aligned}a_\mathrm{rel} = \left( \frac{G\, m_\mathrm{total}\, P^2}{4\pi^2}\right)^{1/3}\\m_\mathrm{body} = m_\mathrm{total}\left( 1 - \frac{a_\mathrm{body}}{a_\mathrm{rel}}\right)\end{aligned}\end{align} \]
Parameters:

m_total (Real[Quantity[PhysicalType('mass')], '']) – Total system mass (sum of all bodies).

Return type:

Real[Quantity[PhysicalType('mass')], '']

Returns:

The mass of this body.

get_position(time, orientation=None)

Get 3D position of the body in its orbit at given time(s).

By definition and convention of this class, this is the position of the body relative to the system barycenter, accounting for the orbit orientation.

Parameters:
Return type:

Real[Quantity[PhysicalType('length')], '3 *batch']

get_velocity(time, orientation=None)

Get 3D velocity of the body relative to the system barycenter.

Parameters:
Return type:

Real[Quantity[PhysicalType({'speed', 'velocity'})], '3 *batch']

orientation: KeplerianOrientation = KeplerianOrientation()
period: Real[Quantity[PhysicalType('time')], '']
eccentricity: Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '']
semi_major_axis: Real[Quantity[PhysicalType('length')], '']
t_peri: Real[Quantity[PhysicalType('time')], '']
class harv.kepler.KeplerianOrientation

Bases: Module

Orientation of a Keplerian orbit in 3D space.

Stores the three Euler angles that define how the orbital plane is oriented relative to the observer’s reference frame: - Inclination (i): tilt of orbital plane from sky plane - Longitude of ascending node (Omega): where orbit crosses sky plane - Argument of pericenter (omega): orientation of ellipse within orbital plane

Angles are stored as sin/cos pairs for numerical stability.

Examples

>>> from unxt import Q
>>> from harv.kepler.orientation import KeplerianOrientation
>>> orient = KeplerianOrientation.from_angles(
...     arg_peri=Q(0.5, "rad"),
...     lon_asc_node=Q(1.0, "rad"),
...     inclination=Q(0.3, "rad"),
... )
Parameters:

  • sin_arg_peri (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • cos_arg_peri (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • sin_lon_asc_node (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • cos_lon_asc_node (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • sin_i (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • cos_i (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

__init__(sin_arg_peri=0.0, cos_arg_peri=1.0, sin_lon_asc_node=0.0, cos_lon_asc_node=1.0, sin_i=0.0, cos_i=1.0)
Parameters:

  • sin_arg_peri (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • cos_arg_peri (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • sin_lon_asc_node (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • cos_lon_asc_node (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • sin_i (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • cos_i (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

Return type:

None

static __new__(cls, *args, **kwargs)
Parameters:
Return type:

TypeVar(_ModuleT, bound= Module)

property arg_peri: Real[Quantity[PhysicalType('angle')], '']

Argument of pericenter (omega).

cos_arg_peri: Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''] = 1.0
cos_i: Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''] = 1.0
cos_lon_asc_node: Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''] = 1.0
classmethod from_angles(arg_peri=Quantity(Array(0., dtype=float32, weak_type=True), unit='rad'), lon_asc_node=Quantity(Array(0., dtype=float32, weak_type=True), unit='rad'), inclination=Quantity(Array(0., dtype=float32, weak_type=True), unit='rad'))

Construct from angle values.

Parameters:

  • arg_peri (Real[Quantity[PhysicalType('angle')], ''])

  • lon_asc_node (Real[Quantity[PhysicalType('angle')], ''])

  • inclination (Real[Quantity[PhysicalType('angle')], ''])

Return type:

KeplerianOrientation

classmethod from_thiele_innes(A, B, F, G)

Construct orientation from Thiele-Innes constants.

Inverts the Thiele-Innes constants to recover (omega, Omega, i, a). This loosely follows Appendix A of https://arxiv.org/abs/2206.05726 but my convention for angles is different. For this implementation:

0 < i < pi/2 0 < omega < 2pi 0 < Omega < 2pi

whereas the paper linked above assumes:

0 < i < pi 0 < omega < 2pi 0 < Omega < pi

Return type:

tuple[KeplerianOrientation, Real[Quantity[PhysicalType('length')], ''] | Real[Quantity[PhysicalType('angle')], '']]

Returns:

  • orientation – KeplerianOrientation object

  • semi_major_axis – Recovered semi-major axis

Parameters:

  • A (Real[Quantity[PhysicalType('length')], ''] | Real[Quantity[PhysicalType('angle')], ''])

  • B (Real[Quantity[PhysicalType('length')], ''] | Real[Quantity[PhysicalType('angle')], ''])

  • F (Real[Quantity[PhysicalType('length')], ''] | Real[Quantity[PhysicalType('angle')], ''])

  • G (Real[Quantity[PhysicalType('length')], ''] | Real[Quantity[PhysicalType('angle')], ''])

property inclination: Real[Quantity[PhysicalType('angle')], '']

Inclination (i).

property lon_asc_node: Real[Quantity[PhysicalType('angle')], '']

Longitude of ascending node (Omega).

property rotation_matrix: Array

Compute rotation matrix from orbital plane to observer frame.

Returns the rotation matrix R such that: r_observer_frame = R @ r_orbital_frame

The rotation is composed of three sequential rotations: 1. R_z(omega): Rotate by argument of pericenter, omega, in orbital plane 2. R_x(i): Rotate by inclination, i, to tilt orbital plane 3. R_z(Omega): Rotate by longitude of ascending node, Omega, on sky plane

The full rotation matrix is therefore: R = R_z(Omega) @ R_x(i) @ R_z(omega)

We build the matrix directly from the sin/cos pairs for numerical stability and speed, but using the notation below, it is equivalent to:

R1 = jnp.array([[c_w, -s_w, 0], [s_w, c_w, 0], [0, 0, 1]]) R2 = jnp.array([[1., 0, 0], [0, c_i, -s_i], [0, s_i, c_i]]) R3 = jnp.array([[c_W, -s_W, 0], [s_W, c_W, 0], [0, 0, 1]]) R = R3 @ R2 @ R1

Or, alternatively: omega = arg_peri.to_value(“rad”) Omega = lon_asc_node.to_value(“rad”) i = inclination.to_value(“rad”) R = Rotation.from_euler(‘ZXZ’, [Omega, i, omega]).as_matrix()

sin_arg_peri: Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''] = 0.0
sin_i: Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''] = 0.0
sin_lon_asc_node: Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''] = 0.0
thiele_innes_constants(semi_major_axis=None)

Compute Thiele-Innes constants (A, B, F, G).

These constants linearize the relationship between orbital position and sky-plane projection. See Appendix A of: https://arxiv.org/abs/2206.05726

Parameters:

semi_major_axis (Real[Quantity[PhysicalType('length')], ''] | Real[Quantity[PhysicalType('angle')], ''] | None) – Semi-major axis of the orbit

Returns:

The four Thiele-Innes constants.

Return type:

A, B, F, G

class harv.kepler.AbstractNBodySystem

Bases: Module

Abstract base class for Keplerian N-body systems.

__init__()
Return type:

None

static __new__(cls, *args, **kwargs)
Parameters:
Return type:

TypeVar(_ModuleT, bound= Module)

property m_total: Real[Quantity[PhysicalType('mass')], '']

Total system mass.

property n_bodies: int

Total number of bodies (primary + companions).

class harv.kepler.TwoBodySystem

Bases: AbstractNBodySystem

A system with a primary body and one companion.

Bodies are indexed as: - 0: Primary body - 1: Companion

Parameters:

  • m_primary (Real[Quantity[PhysicalType('mass')], ''])

  • companion (KeplerianBody)

__init__(m_primary, companion)
Parameters:

  • m_primary (Real[Quantity[PhysicalType('mass')], ''])

  • companion (KeplerianBody)

Return type:

None

static __new__(cls, *args, **kwargs)
Parameters:
Return type:

TypeVar(_ModuleT, bound= Module)

property m_companion: Real[Quantity[PhysicalType('mass')], '']

Companion mass.

property m_total: Real[Quantity[PhysicalType('mass')], '']

Total system mass, derived from Kepler’s 3rd law.

The companion’s semi-major axis is defined relative to the system barycenter, so we can use Kepler’s 3rd law to derive the total mass from the companion’s orbital parameters and the primary mass: a_body = a_rel * m_primary / m_total a_rel^3 = G * m_total * P^2 / 4 pi^2 a_body^3 = G * m_primary^3 * P^2 / (4 pi^2 * m_total^2) (solve for m_total)

property n_bodies: int

Total number of bodies (primary + companions).

position_barycentric(time, body_idx)

Get barycentric position of specified body at given time(s).

Parameters:

  • time (Real[Quantity[PhysicalType('time')], '*batch']) – Time(s) to evaluate

  • body_idx (int) – Index of body (0=primary, 1=companion)

Return type:

Real[Quantity[PhysicalType('length')], '3 *batch']

Returns:

3D position vector(s) of specified body relative to barycenter

position_relative(time)

Position of companion relative to primary.

Parameters:

time (Real[Quantity[PhysicalType('time')], '*batch']) – Time(s) to evaluate

Return type:

Real[Quantity[PhysicalType('length')], '3 *batch']

Returns:

3D position vector(s) of companion relative to primary

velocity_barycentric(time, body_idx)

Get barycentric velocity of specified body at given time(s).

Parameters:

  • time (Real[Quantity[PhysicalType('time')], '*batch']) – Time(s) to evaluate

  • body_idx (int) – Index of body (0=primary, 1=companion)

Return type:

Real[Quantity[PhysicalType({'speed', 'velocity'})], '3 *batch']

Returns:

3D velocity vector(s) of specified body relative to barycenter

velocity_relative(time)

Velocity of companion relative to primary.

Parameters:

time (Real[Quantity[PhysicalType('time')], '*batch']) – Time(s) to evaluate

Return type:

Real[Quantity[PhysicalType({'speed', 'velocity'})], '3 *batch']

Returns:

3D velocity vector(s) of companion relative to primary

m_primary: Real[Quantity[PhysicalType('mass')], '']
companion: KeplerianBody
harv.kepler.astrometric_mass_function(a_physical, period)

Astrometric mass function from a physical orbit size and period.

\[f(m) = \frac{4\pi^2\,a^3}{G\,P^2}\]

With a the semi-major axis of the primary’s barycentric orbit – the photocentre orbit when the companion is dark or faint – this equals \(m_2^3 / (m_1 + m_2)^2\). See docs/spec.md (“Mass functions”).

Parameters:

  • a_physical (Real[Quantity[PhysicalType('length')], '*batch']) – Physical semi-major axis (a length), e.g. from semi_major_axis_physical().

  • period (Real[Quantity[PhysicalType('time')], '*batch']) – Orbital period.

Return type:

Real[Quantity[PhysicalType('mass')], '*batch']

Returns:

The astrometric mass function, in solar masses.

Examples

>>> from unxt import Q, ustrip
>>> from harv.kepler.masses import astrometric_mass_function
>>> mf = astrometric_mass_function(Q(1.0, "AU"), Q(1.0, "yr"))
>>> round(float(ustrip("Msun", mf)), 4)
1.0
harv.kepler.astrometric_orbit_at_times(times, period, eccentricity, t_peri, arg_peri, cos_i, lon_asc_node, semi_major_axis)

Compute sky-plane astrometric orbit (Deltara, Deltadec) at given times.

Uses the Thiele-Innes parameterization following the Gaia local plane coordinate (LPC) convention (Lindegren & Bastian, GAIA-C3-TN-LU-LL-061-08, Eq. 4):

Deltara  = (B*cos f + G*sin f) * a      (RA / ``a`` direction)
Deltadec = (A*cos f + F*sin f) * a      (Dec / ``d`` direction)

where A, B, F, G are the unit Thiele-Innes constants and a is the photocentric semi-major axis.

Parameters:

  • times (Real[Quantity[PhysicalType('time')], '*batch']) – Observation times.

  • period (Real[Quantity[PhysicalType('time')], '']) – Orbital period.

  • eccentricity (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '']) – Orbital eccentricity.

  • t_peri (Real[Quantity[PhysicalType('time')], '']) – Time of periastron passage. In the likelihood layer this is derived from the dimensionless phase_peri as t_peri = phase_peri * period (see _solve_kepler).

  • arg_peri (Real[Quantity[PhysicalType('angle')], '']) – Argument of periastron omega.

  • cos_i (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '']) – Cosine of orbital inclination.

  • lon_asc_node (Real[Quantity[PhysicalType('angle')], '']) – Longitude of the ascending node Omega.

  • semi_major_axis (Real[Quantity[PhysicalType('angle')], '']) – Photocentric semi-major axis.

Return type:

tuple[Real[Quantity[PhysicalType('angle')], '*batch'], Real[Quantity[PhysicalType('angle')], '*batch']]

Returns:

(delta_ra, delta_dec) — sky-plane offsets in the same unit as semi_major_axis.

Examples

>>> from unxt import Q
>>> from harv.kepler.orbits import astrometric_orbit_at_times
>>> times = Q([0.0, 100.0, 200.0], "day")
>>> dra, ddec = astrometric_orbit_at_times(
...     times,
...     period=Q(300.0, "day"),
...     eccentricity=0.3,
...     t_peri=Q(0.0, "day"),
...     arg_peri=Q(1.2, "rad"),
...     cos_i=0.5,
...     lon_asc_node=Q(0.8, "rad"),
...     semi_major_axis=Q(3.0, "mas"),
... )
>>> dra.unit
Unit("mas")
harv.kepler.binary_mass_function(period, rv_semiamp, eccentricity)

Binary mass function from radial-velocity orbital elements.

\[f(m) = \frac{m_2^3 \sin^3 i}{(m_1 + m_2)^2} = \frac{P\,K^3\,(1 - e^2)^{3/2}}{2\pi G}\]
Parameters:

  • period (Real[Quantity[PhysicalType('time')], '*batch']) – Orbital period.

  • rv_semiamp (Real[Quantity[PhysicalType({'speed', 'velocity'})], '*batch']) – Radial-velocity semi-amplitude K.

  • eccentricity (Float[Array, '*batch'] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '*batch']) – Orbital eccentricity (dimensionless).

Return type:

Real[Quantity[PhysicalType('mass')], '*batch']

Returns:

The binary mass function, in solar masses.

Examples

>>> from unxt import Q, ustrip
>>> from harv.kepler.masses import binary_mass_function
>>> mf = binary_mass_function(Q(100.0, "day"), Q(10.0, "km/s"), 0.1)
>>> round(float(ustrip("Msun", mf)), 4)
0.0102
harv.kepler.companion_mass_from_mass_function(mass_function, m1, sini=1.0)

Solve the mass function for the companion mass \(m_2\).

Inverts \(m_2^3 \sin^3 i / (m_1 + m_2)^2 = f\) for \(m_2\), given the mass function f, primary mass m1, and sin i. sini = 1 (the default) yields the minimum companion mass.

The cubic has a single positive root, so it is solved by bisection – which is robust and jax.jit / jax.vmap friendly. The root is bracketed by [0, max(4 f_eff, (4 f_eff m_1^2)^{1/3})] with \(f_\mathrm{eff} = f / \sin^3 i\).

Parameters:

  • mass_function (Real[Quantity[PhysicalType('mass')], '*batch']) – Binary or astrometric mass function (a mass).

  • m1 (Quantity[PhysicalType('mass')]) – Primary mass.

  • sini (Float[Array, '*batch'] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '*batch']) – Sine of the orbital inclination. Default 1 (edge-on -> minimum mass).

Return type:

Real[Quantity[PhysicalType('mass')], '*batch']

Returns:

The companion mass \(m_2\), in solar masses.

Examples

>>> from unxt import Q, ustrip
>>> from harv.kepler.masses import companion_mass_from_mass_function
>>> m2 = companion_mass_from_mass_function(Q(0.25, "Msun"), Q(1.0, "Msun"))
>>> round(float(ustrip("Msun", m2)), 3)
1.0
harv.kepler.compute_true_anomaly_components(time, period, eccentricity, t_peri)

Compute true anomaly at given times.

Parameters:

  • time (Real[Quantity[PhysicalType('time')], '*batch']) – Observation times, shape (n,)

  • period (Real[Quantity[PhysicalType('time')], '']) – Orbital period

  • eccentricity (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '']) – Orbital eccentricity

  • t_peri (Real[Quantity[PhysicalType('time')], '']) – Time of pericenter passage

Return type:

tuple[Float[Array, '*batch'] | Real[Quantity[PhysicalType('dimensionless')], '*batch'], Float[Array, '*batch'] | Real[Quantity[PhysicalType('dimensionless')], '*batch']]

Returns:

(sin_f, cos_f) — true anomaly components, each shape (n,).

Examples

>>> from unxt import Q
>>> from harv.kepler.orbits import compute_true_anomaly_components
>>> sin_f, cos_f = compute_true_anomaly_components(
...     time=Q([0.0, 25.0, 50.0], "day"),
...     period=Q(100.0, "day"),
...     eccentricity=0.3,
...     t_peri=Q(0.0, "day"),
... )
harv.kepler.rv_at_times(times, period, eccentricity, t_peri, arg_peri, rv_semiamp, v_sys)

Compute the RV model: K*[cos(omega + f(t)) + e*cos(omega)] + v_0.

Parameters:

  • times (Real[Quantity[PhysicalType('time')], '*batch']) – Observation times.

  • period (Real[Quantity[PhysicalType('time')], '']) – Orbital period.

  • eccentricity (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '']) – Orbital eccentricity.

  • t_peri (Real[Quantity[PhysicalType('time')], '']) – Time of periastron passage. In the likelihood layer this is derived from the dimensionless phase_peri as t_peri = phase_peri * period (see _solve_kepler).

  • arg_peri (Real[Quantity[PhysicalType('angle')], '']) – Argument of periastron omega.

  • rv_semiamp (Real[Quantity[PhysicalType({'speed', 'velocity'})], '']) – RV semi-amplitude.

  • v_sys (Real[Quantity[PhysicalType({'speed', 'velocity'})], '']) – Systemic velocity.

Return type:

Real[Quantity[PhysicalType({'speed', 'velocity'})], '*batch']

Returns:

Radial velocities in the same unit as rv_semiamp and v_sys.

Examples

>>> from unxt import Q
>>> from harv.kepler.orbits import rv_at_times
>>> times = Q([0.0, 50.0, 100.0], "day")
>>> rv = rv_at_times(
...     times,
...     period=Q(200.0, "day"),
...     eccentricity=0.3,
...     t_peri=Q(50.0, "day"),
...     arg_peri=Q(1.2, "rad"),
...     rv_semiamp=Q(8.0, "km/s"),
...     v_sys=Q(-5.0, "km/s"),
... )
>>> rv.unit
Unit("km / s")
harv.kepler.mean_anomaly(dt, period)

Compute mean anomaly from elapsed time and period.

M = 2pi * dt / period, returned as a Q with angle units (radians).

Examples

>>> from unxt import Q
>>> from harv.kepler.orbits import mean_anomaly
>>> M = mean_anomaly(Q(50.0, "day"), Q(100.0, "day"))
>>> M.unit
Unit("rad")
Parameters:

  • dt (Real[Quantity[PhysicalType('time')], '*batch'])

  • period (Real[Quantity[PhysicalType('time')], ''])

Return type:

Real[Quantity[PhysicalType('angle')], '*batch']

harv.kepler.rv_shape(sin_f, cos_f, eccentricity, arg_peri)

RV shape function: cos(omega + f) + e*cos(omega).

Returns the dimensionless RV amplitude factor for each observation. arg_peri may be a plain float/array in radians or a Q with angle units; jnp.cos handles both via quax dispatch. eccentricity may likewise be a dimensionless Q or a plain scalar.

Examples

>>> from unxt import Q
>>> from harv.kepler.orbits import true_anomaly_from_mean, rv_shape
>>> sin_f, cos_f = true_anomaly_from_mean(Q(1.0, "rad"), 0.3)
>>> shape = rv_shape(sin_f, cos_f, 0.3, Q(0.5, "rad"))
Parameters:

  • sin_f (Float[Array, '*batch'] | Real[Quantity[PhysicalType('dimensionless')], '*batch'])

  • cos_f (Float[Array, '*batch'] | Real[Quantity[PhysicalType('dimensionless')], '*batch'])

  • eccentricity (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • arg_peri (Real[Quantity[PhysicalType('angle')], ''] | Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

Return type:

Float[Array, '*batch'] | Real[Quantity[PhysicalType('dimensionless')], '*batch']

harv.kepler.semi_major_axis_physical(a_angular, parallax)

Physical semi-major axis (AU) from an angular size and parallax.

For an orbit at heliocentric distance \(d\), the angular semi-major axis is \(\alpha = a / d\) and the parallax is \(\varpi = 1\,\mathrm{AU} / d\), so \(a = (\alpha / \varpi)\; \mathrm{AU}\).

Parameters:

  • a_angular (Real[Quantity[PhysicalType('angle')], '*batch']) – Angular semi-major axis (an angle, e.g. mas).

  • parallax (Real[Quantity[PhysicalType('angle')], '*batch']) – Parallax (an angle, e.g. mas).

Return type:

Real[Quantity[PhysicalType('length')], '*batch']

Returns:

The physical semi-major axis, in AU.

Examples

>>> from unxt import Q
>>> from harv.kepler.masses import semi_major_axis_physical
>>> a = semi_major_axis_physical(Q(10.0, "mas"), Q(5.0, "mas"))
>>> float(a.value)
2.0
harv.kepler.thiele_innes_ABFG(cos_arg_peri, sin_arg_peri, cos_lon_asc_node, sin_lon_asc_node, cos_i)

Compute unit Thiele-Innes constants (A, B, F, G).

Returns the constants with an implicit semi-major axis of 1. Multiply each by a to recover the physical Thiele-Innes constants.

Inputs are dimensionless and may be scalar or batched (plain scalars, JAX arrays, or dimensionless Q); the computation broadcasts elementwise.

The sky-plane orbital displacement uses the Thiele-Innes coordinates X = (cos E - e) and Y = sqrt(1-e^2) sin E, or equivalently in terms of true anomaly: X = r/a cos f, Y = r/a sin f where r/a = (1-e^2)/(1+e cos f).

In the Gaia LPC convention (Lindegren & Bastian, GAIA-C3-TN-LU-LL-061-08, Eq. 4), B and G project into the RA (a) direction, while A and F project into the Dec (d) direction.

See also Eq. A.1 of https://arxiv.org/abs/2206.05726.

Examples

>>> import quaxed.numpy as jnp
>>> from unxt import Q
>>> from harv.kepler.orbits import thiele_innes_ABFG
>>> A, B, F, G = thiele_innes_ABFG(
...     cos_arg_peri=jnp.cos(Q(0.5, "rad")),
...     sin_arg_peri=jnp.sin(Q(0.5, "rad")),
...     cos_lon_asc_node=jnp.cos(Q(1.0, "rad")),
...     sin_lon_asc_node=jnp.sin(Q(1.0, "rad")),
...     cos_i=jnp.cos(Q(0.3, "rad")),
... )
Parameters:

  • cos_arg_peri (Float[Array, '*batch'] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '*batch'])

  • sin_arg_peri (Float[Array, '*batch'] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '*batch'])

  • cos_lon_asc_node (Float[Array, '*batch'] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '*batch'])

  • sin_lon_asc_node (Float[Array, '*batch'] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '*batch'])

  • cos_i (Float[Array, '*batch'] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '*batch'])

Return type:

tuple[Float[Array, '*batch'] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '*batch'], Float[Array, '*batch'] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '*batch'], Float[Array, '*batch'] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '*batch'], Float[Array, '*batch'] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '*batch']]

harv.kepler.true_anomaly_from_mean(M, eccentricity)

Solve Kepler’s equation: mean anomaly -> (sin f, cos f).

Wraps jaxoplanet.core.kepler.kepler. The mean anomaly is stripped to radians internally.

Examples

>>> from unxt import Q
>>> from harv.kepler.orbits import true_anomaly_from_mean
>>> sin_f, cos_f = true_anomaly_from_mean(Q(1.0, "rad"), 0.3)
Parameters:

  • M (Real[Quantity[PhysicalType('angle')], '*batch'])

  • eccentricity (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

Return type:

tuple[Float[Array, '*batch'] | Real[Quantity[PhysicalType('dimensionless')], '*batch'], Float[Array, '*batch'] | Real[Quantity[PhysicalType('dimensionless')], '*batch']]

Submodules

harv.kepler.body module

Keplerian orbit implementation with units support and JAX compatibility.

class harv.kepler.body.KeplerianBody

Bases: Module

Orbital parameters of a Keplerian body (companion, i.e. body 2).

This class represents the orbital parameters of a companion or second body (relative to some barycenter). So, all parameters represent the orbital elements of a specific body relative to the barycenter.

The primary parameterization uses: - Period (P) - Eccentricity (e) - Semi-major axis (a) or semi_major_axis - Time of pericenter or t_peri

Alternative constructors support different parameterizations.

Parameters:

  • period (Real[Quantity[PhysicalType('time')], '']) – Orbital period.

  • eccentricity (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '']) – Orbital eccentricity.

  • semi_major_axis (Real[Quantity[PhysicalType('length')], '']) – The semi-major axis of the body relative to its system barycenter.

  • t_peri (Real[Quantity[PhysicalType('time')], '']) – Time of pericenter passage.

  • orientation (KeplerianOrientation) – Optional: Orientation of the orbit.

  • ecc_zero_tol (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '']) – Optional: Tolerance for treating eccentricity as zero in circular orbit shortcuts. This is used to avoid numerical issues in the Kepler solver when eccentricity is very small, and can be set to a small multiple of machine epsilon for the data type used.

Examples

>>> from unxt import Q
>>> from harv.kepler.body import KeplerianBody
>>> body = KeplerianBody(
...     period=Q(365.25, "day"),
...     eccentricity=0.1,
...     semi_major_axis=Q(1.0, "AU"),
...     t_peri=Q(0.0, "day"),
... )
period: Real[Quantity[PhysicalType('time')], '']
eccentricity: Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '']
semi_major_axis: Real[Quantity[PhysicalType('length')], '']
t_peri: Real[Quantity[PhysicalType('time')], '']
orientation: KeplerianOrientation = KeplerianOrientation()
ecc_zero_tol: Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''] = np.float64(2.220446049250313e-15)
classmethod from_masses(period, eccentricity, m_total, m_body, t_peri, **kwargs)

Construct body’s barycentric orbit from masses and period.

Computes this body’s barycentric semi-major axis from Kepler’s 3rd law: 1. Compute relative orbit: a_rel = (G m_total P^2 / 4 pi^2)^(1/3) 2. Convert to barycentric: a_body = a_rel * (1 - m_body / m_total)

Parameters:

  • period (Real[Quantity[PhysicalType('time')], '']) – Orbital period.

  • eccentricity (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '']) – Orbital eccentricity.

  • m_total (Real[Quantity[PhysicalType('mass')], '']) – Total system mass.

  • m_body (Real[Quantity[PhysicalType('mass')], '']) – Mass of this body.

  • t_peri (Real[Quantity[PhysicalType('time')], '']) – Time of pericenter passage.

  • orientation – Optional: Orientation of the orbit.

  • kwargs (Any) – Additional keyword arguments to pass to the main constructor.

Return type:

KeplerianBody

Returns:

The body’s orbit about the system barycenter.

get_mass(m_total)

Recover this body’s mass from the total system mass.

This inverts from_masses: given the total system mass and the stored barycentric semi-major axis, it recovers the body mass via

\[ \begin{align}\begin{aligned}a_\mathrm{rel} = \left( \frac{G\, m_\mathrm{total}\, P^2}{4\pi^2}\right)^{1/3}\\m_\mathrm{body} = m_\mathrm{total}\left( 1 - \frac{a_\mathrm{body}}{a_\mathrm{rel}}\right)\end{aligned}\end{align} \]
Parameters:

m_total (Real[Quantity[PhysicalType('mass')], '']) – Total system mass (sum of all bodies).

Return type:

Real[Quantity[PhysicalType('mass')], '']

Returns:

The mass of this body.

get_position(time, orientation=None)

Get 3D position of the body in its orbit at given time(s).

By definition and convention of this class, this is the position of the body relative to the system barycenter, accounting for the orbit orientation.

Parameters:
Return type:

Real[Quantity[PhysicalType('length')], '3 *batch']

get_velocity(time, orientation=None)

Get 3D velocity of the body relative to the system barycenter.

Parameters:
Return type:

Real[Quantity[PhysicalType({'speed', 'velocity'})], '3 *batch']

__init__(period, eccentricity, semi_major_axis, t_peri, orientation=KeplerianOrientation(), *, ecc_zero_tol=np.float64(2.220446049250313e-15))
Parameters:

  • period (Real[Quantity[PhysicalType('time')], ''])

  • eccentricity (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • semi_major_axis (Real[Quantity[PhysicalType('length')], ''])

  • t_peri (Real[Quantity[PhysicalType('time')], ''])

  • orientation (KeplerianOrientation)

  • ecc_zero_tol (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

Return type:

None

static __new__(cls, *args, **kwargs)
Parameters:
Return type:

TypeVar(_ModuleT, bound= Module)

harv.kepler.constants module

Constants used in Keplerian calculations.

harv.kepler.masses module

Mass functions and physical-quantity helpers for Keplerian orbits.

These functions turn posterior orbital elements into physical masses and physical (AU) orbit sizes. They are pure, unit-aware, and broadcast over batched inputs (one value per posterior sample), so callers can apply them directly to arrays of samples. See docs/spec.md (“Mass functions”).

harv.kepler.masses.binary_mass_function(period, rv_semiamp, eccentricity)

Binary mass function from radial-velocity orbital elements.

\[f(m) = \frac{m_2^3 \sin^3 i}{(m_1 + m_2)^2} = \frac{P\,K^3\,(1 - e^2)^{3/2}}{2\pi G}\]
Parameters:

  • period (Real[Quantity[PhysicalType('time')], '*batch']) – Orbital period.

  • rv_semiamp (Real[Quantity[PhysicalType({'speed', 'velocity'})], '*batch']) – Radial-velocity semi-amplitude K.

  • eccentricity (Float[Array, '*batch'] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '*batch']) – Orbital eccentricity (dimensionless).

Return type:

Real[Quantity[PhysicalType('mass')], '*batch']

Returns:

The binary mass function, in solar masses.

Examples

>>> from unxt import Q, ustrip
>>> from harv.kepler.masses import binary_mass_function
>>> mf = binary_mass_function(Q(100.0, "day"), Q(10.0, "km/s"), 0.1)
>>> round(float(ustrip("Msun", mf)), 4)
0.0102
harv.kepler.masses.astrometric_mass_function(a_physical, period)

Astrometric mass function from a physical orbit size and period.

\[f(m) = \frac{4\pi^2\,a^3}{G\,P^2}\]

With a the semi-major axis of the primary’s barycentric orbit – the photocentre orbit when the companion is dark or faint – this equals \(m_2^3 / (m_1 + m_2)^2\). See docs/spec.md (“Mass functions”).

Parameters:

  • a_physical (Real[Quantity[PhysicalType('length')], '*batch']) – Physical semi-major axis (a length), e.g. from semi_major_axis_physical().

  • period (Real[Quantity[PhysicalType('time')], '*batch']) – Orbital period.

Return type:

Real[Quantity[PhysicalType('mass')], '*batch']

Returns:

The astrometric mass function, in solar masses.

Examples

>>> from unxt import Q, ustrip
>>> from harv.kepler.masses import astrometric_mass_function
>>> mf = astrometric_mass_function(Q(1.0, "AU"), Q(1.0, "yr"))
>>> round(float(ustrip("Msun", mf)), 4)
1.0
harv.kepler.masses.companion_mass_from_mass_function(mass_function, m1, sini=1.0)

Solve the mass function for the companion mass \(m_2\).

Inverts \(m_2^3 \sin^3 i / (m_1 + m_2)^2 = f\) for \(m_2\), given the mass function f, primary mass m1, and sin i. sini = 1 (the default) yields the minimum companion mass.

The cubic has a single positive root, so it is solved by bisection – which is robust and jax.jit / jax.vmap friendly. The root is bracketed by [0, max(4 f_eff, (4 f_eff m_1^2)^{1/3})] with \(f_\mathrm{eff} = f / \sin^3 i\).

Parameters:

  • mass_function (Real[Quantity[PhysicalType('mass')], '*batch']) – Binary or astrometric mass function (a mass).

  • m1 (Quantity[PhysicalType('mass')]) – Primary mass.

  • sini (Float[Array, '*batch'] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '*batch']) – Sine of the orbital inclination. Default 1 (edge-on -> minimum mass).

Return type:

Real[Quantity[PhysicalType('mass')], '*batch']

Returns:

The companion mass \(m_2\), in solar masses.

Examples

>>> from unxt import Q, ustrip
>>> from harv.kepler.masses import companion_mass_from_mass_function
>>> m2 = companion_mass_from_mass_function(Q(0.25, "Msun"), Q(1.0, "Msun"))
>>> round(float(ustrip("Msun", m2)), 3)
1.0
harv.kepler.masses.semi_major_axis_physical(a_angular, parallax)

Physical semi-major axis (AU) from an angular size and parallax.

For an orbit at heliocentric distance \(d\), the angular semi-major axis is \(\alpha = a / d\) and the parallax is \(\varpi = 1\,\mathrm{AU} / d\), so \(a = (\alpha / \varpi)\; \mathrm{AU}\).

Parameters:

  • a_angular (Real[Quantity[PhysicalType('angle')], '*batch']) – Angular semi-major axis (an angle, e.g. mas).

  • parallax (Real[Quantity[PhysicalType('angle')], '*batch']) – Parallax (an angle, e.g. mas).

Return type:

Real[Quantity[PhysicalType('length')], '*batch']

Returns:

The physical semi-major axis, in AU.

Examples

>>> from unxt import Q
>>> from harv.kepler.masses import semi_major_axis_physical
>>> a = semi_major_axis_physical(Q(10.0, "mas"), Q(5.0, "mas"))
>>> float(a.value)
2.0

harv.kepler.nbody_system module

Keplerian orbit implementation with units support and JAX compatibility.

class harv.kepler.nbody_system.AbstractNBodySystem

Bases: Module

Abstract base class for Keplerian N-body systems.

property n_bodies: int

Total number of bodies (primary + companions).

property m_total: Real[Quantity[PhysicalType('mass')], '']

Total system mass.

__init__()
Return type:

None

static __new__(cls, *args, **kwargs)
Parameters:
Return type:

TypeVar(_ModuleT, bound= Module)

class harv.kepler.nbody_system.TwoBodySystem

Bases: AbstractNBodySystem

A system with a primary body and one companion.

Bodies are indexed as: - 0: Primary body - 1: Companion

Parameters:

  • m_primary (Real[Quantity[PhysicalType('mass')], ''])

  • companion (KeplerianBody)

m_primary: Real[Quantity[PhysicalType('mass')], '']
companion: KeplerianBody
property n_bodies: int

Total number of bodies (primary + companions).

property m_total: Real[Quantity[PhysicalType('mass')], '']

Total system mass, derived from Kepler’s 3rd law.

The companion’s semi-major axis is defined relative to the system barycenter, so we can use Kepler’s 3rd law to derive the total mass from the companion’s orbital parameters and the primary mass: a_body = a_rel * m_primary / m_total a_rel^3 = G * m_total * P^2 / 4 pi^2 a_body^3 = G * m_primary^3 * P^2 / (4 pi^2 * m_total^2) (solve for m_total)

property m_companion: Real[Quantity[PhysicalType('mass')], '']

Companion mass.

position_barycentric(time, body_idx)

Get barycentric position of specified body at given time(s).

Parameters:

  • time (Real[Quantity[PhysicalType('time')], '*batch']) – Time(s) to evaluate

  • body_idx (int) – Index of body (0=primary, 1=companion)

Return type:

Real[Quantity[PhysicalType('length')], '3 *batch']

Returns:

3D position vector(s) of specified body relative to barycenter

position_relative(time)

Position of companion relative to primary.

Parameters:

time (Real[Quantity[PhysicalType('time')], '*batch']) – Time(s) to evaluate

Return type:

Real[Quantity[PhysicalType('length')], '3 *batch']

Returns:

3D position vector(s) of companion relative to primary

velocity_barycentric(time, body_idx)

Get barycentric velocity of specified body at given time(s).

Parameters:

  • time (Real[Quantity[PhysicalType('time')], '*batch']) – Time(s) to evaluate

  • body_idx (int) – Index of body (0=primary, 1=companion)

Return type:

Real[Quantity[PhysicalType({'speed', 'velocity'})], '3 *batch']

Returns:

3D velocity vector(s) of specified body relative to barycenter

velocity_relative(time)

Velocity of companion relative to primary.

Parameters:

time (Real[Quantity[PhysicalType('time')], '*batch']) – Time(s) to evaluate

Return type:

Real[Quantity[PhysicalType({'speed', 'velocity'})], '3 *batch']

Returns:

3D velocity vector(s) of companion relative to primary

__init__(m_primary, companion)
Parameters:

  • m_primary (Real[Quantity[PhysicalType('mass')], ''])

  • companion (KeplerianBody)

Return type:

None

static __new__(cls, *args, **kwargs)
Parameters:
Return type:

TypeVar(_ModuleT, bound= Module)

harv.kepler.orbits module

Helper functions for Keplerian orbits.

Implementations of core orbit computations used by harv.kepler, harv.likelihood, and harv.simulate.

All functions accept Q objects (including dimensionless ones) as well as plain JAX arrays and Python scalars.

harv.kepler.orbits.mean_anomaly(dt, period)

Compute mean anomaly from elapsed time and period.

M = 2pi * dt / period, returned as a Q with angle units (radians).

Examples

>>> from unxt import Q
>>> from harv.kepler.orbits import mean_anomaly
>>> M = mean_anomaly(Q(50.0, "day"), Q(100.0, "day"))
>>> M.unit
Unit("rad")
Parameters:

  • dt (Real[Quantity[PhysicalType('time')], '*batch'])

  • period (Real[Quantity[PhysicalType('time')], ''])

Return type:

Real[Quantity[PhysicalType('angle')], '*batch']

harv.kepler.orbits.rv_shape(sin_f, cos_f, eccentricity, arg_peri)

RV shape function: cos(omega + f) + e*cos(omega).

Returns the dimensionless RV amplitude factor for each observation. arg_peri may be a plain float/array in radians or a Q with angle units; jnp.cos handles both via quax dispatch. eccentricity may likewise be a dimensionless Q or a plain scalar.

Examples

>>> from unxt import Q
>>> from harv.kepler.orbits import true_anomaly_from_mean, rv_shape
>>> sin_f, cos_f = true_anomaly_from_mean(Q(1.0, "rad"), 0.3)
>>> shape = rv_shape(sin_f, cos_f, 0.3, Q(0.5, "rad"))
Parameters:

  • sin_f (Float[Array, '*batch'] | Real[Quantity[PhysicalType('dimensionless')], '*batch'])

  • cos_f (Float[Array, '*batch'] | Real[Quantity[PhysicalType('dimensionless')], '*batch'])

  • eccentricity (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • arg_peri (Real[Quantity[PhysicalType('angle')], ''] | Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

Return type:

Float[Array, '*batch'] | Real[Quantity[PhysicalType('dimensionless')], '*batch']

harv.kepler.orbits.thiele_innes_ABFG(cos_arg_peri, sin_arg_peri, cos_lon_asc_node, sin_lon_asc_node, cos_i)

Compute unit Thiele-Innes constants (A, B, F, G).

Returns the constants with an implicit semi-major axis of 1. Multiply each by a to recover the physical Thiele-Innes constants.

Inputs are dimensionless and may be scalar or batched (plain scalars, JAX arrays, or dimensionless Q); the computation broadcasts elementwise.

The sky-plane orbital displacement uses the Thiele-Innes coordinates X = (cos E - e) and Y = sqrt(1-e^2) sin E, or equivalently in terms of true anomaly: X = r/a cos f, Y = r/a sin f where r/a = (1-e^2)/(1+e cos f).

In the Gaia LPC convention (Lindegren & Bastian, GAIA-C3-TN-LU-LL-061-08, Eq. 4), B and G project into the RA (a) direction, while A and F project into the Dec (d) direction.

See also Eq. A.1 of https://arxiv.org/abs/2206.05726.

Examples

>>> import quaxed.numpy as jnp
>>> from unxt import Q
>>> from harv.kepler.orbits import thiele_innes_ABFG
>>> A, B, F, G = thiele_innes_ABFG(
...     cos_arg_peri=jnp.cos(Q(0.5, "rad")),
...     sin_arg_peri=jnp.sin(Q(0.5, "rad")),
...     cos_lon_asc_node=jnp.cos(Q(1.0, "rad")),
...     sin_lon_asc_node=jnp.sin(Q(1.0, "rad")),
...     cos_i=jnp.cos(Q(0.3, "rad")),
... )
Parameters:

  • cos_arg_peri (Float[Array, '*batch'] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '*batch'])

  • sin_arg_peri (Float[Array, '*batch'] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '*batch'])

  • cos_lon_asc_node (Float[Array, '*batch'] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '*batch'])

  • sin_lon_asc_node (Float[Array, '*batch'] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '*batch'])

  • cos_i (Float[Array, '*batch'] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '*batch'])

Return type:

tuple[Float[Array, '*batch'] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '*batch'], Float[Array, '*batch'] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '*batch'], Float[Array, '*batch'] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '*batch'], Float[Array, '*batch'] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '*batch']]

harv.kepler.orbits.true_anomaly_from_mean(M, eccentricity)

Solve Kepler’s equation: mean anomaly -> (sin f, cos f).

Wraps jaxoplanet.core.kepler.kepler. The mean anomaly is stripped to radians internally.

Examples

>>> from unxt import Q
>>> from harv.kepler.orbits import true_anomaly_from_mean
>>> sin_f, cos_f = true_anomaly_from_mean(Q(1.0, "rad"), 0.3)
Parameters:

  • M (Real[Quantity[PhysicalType('angle')], '*batch'])

  • eccentricity (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

Return type:

tuple[Float[Array, '*batch'] | Real[Quantity[PhysicalType('dimensionless')], '*batch'], Float[Array, '*batch'] | Real[Quantity[PhysicalType('dimensionless')], '*batch']]

harv.kepler.orbits.campbell_from_thiele_innes(A, B, F, G)

Invert Thiele-Innes constants to Campbell orbital elements.

This follows Halbwachs, Pourbaix, et al. 2023 (see the appendix):

\[\begin{split}u &= \frac{1}{2}(A^2+B^2+F^2+G^2) \\ v &= A\,G - B\,F \\ a_0 &= \sqrt{u + \sqrt{\max(u^2 - v^2, 0)}} \\ \omega + \Omega &= \mathrm{atan2}(B - F, A + G) \\ \omega - \Omega &= \mathrm{atan2}(-B - F, A - G) \\ \cos i &= v / a_0^2\end{split}\]

The returned cos_i is signed: it preserves the sign of \(v = AG - BF\), which determines whether the orbit is prograde (cos_i > 0) or retrograde (cos_i < 0) under the LPC convention. This is required for the round-trip thiele_innes_from_campbell(*campbell_from_thiele_innes(...)) to recover the original TI constants – a positive-only convention would silently flip the sky orbit for retrograde fits.

Parameters:

  • A (Real[Quantity[PhysicalType('angle')], '*batch'])

  • B (Real[Quantity[PhysicalType('angle')], '*batch'])

  • F (Real[Quantity[PhysicalType('angle')], '*batch'])

  • G (Real[Quantity[PhysicalType('angle')], '*batch'])

Return type:

dict[str, Quantity]

harv.kepler.orbits.thiele_innes_from_campbell(semi_major_axis, arg_peri, lon_asc_node, cos_i)

Convert Campbell orbital elements to physical Thiele-Innes constants.

The forward direction of the change of variables inverted by campbell_from_thiele_innes(). The unit Thiele-Innes constants from thiele_innes_ABFG() are scaled by the semi-major axis:

\[(A, B, F, G) = a_0 \cdot \mathrm{thiele\_innes\_ABFG}(\cos\omega, \sin\omega, \cos\Omega, \sin\Omega, \cos i)\]

Examples

>>> from unxt import Q
>>> from harv.kepler.orbits import thiele_innes_from_campbell
>>> A, B, F, G = thiele_innes_from_campbell(
...     Q(2.0, "mas"), Q(0.5, "rad"), Q(1.0, "rad"), Q(0.3, ""),
... )
>>> A.unit
Unit("mas")
Parameters:

  • semi_major_axis (Real[Quantity[PhysicalType('angle')], '*batch'])

  • arg_peri (Real[Quantity[PhysicalType('angle')], '*batch'])

  • lon_asc_node (Real[Quantity[PhysicalType('angle')], '*batch'])

  • cos_i (Real[Quantity[PhysicalType('dimensionless')], '*batch'])

Return type:

tuple[Real[Quantity[PhysicalType('angle')], '*batch'], Real[Quantity[PhysicalType('angle')], '*batch'], Real[Quantity[PhysicalType('angle')], '*batch'], Real[Quantity[PhysicalType('angle')], '*batch']]

harv.kepler.orbits.ecc_omega_from_ecosw_esinw(ecosw, esinw)

Convert (e cos omega, e sin omega) to eccentricity and arg. of periastron.

The inverse of ecosw_esinw_from_ecc_omega():

\[\begin{split}e &= \sqrt{\mathrm{ecosw}^2 + \mathrm{esinw}^2} \\ \omega &= \mathrm{atan2}(\mathrm{esinw}, \mathrm{ecosw})\end{split}\]

The returned arg_peri lies in (-pi, pi] (the range of atan2); use harv.samplers.Samples.wrap_angles() to wrap it into [0, 2*pi) if a prior requires that range.

Examples

>>> from unxt import Q
>>> from harv.kepler.orbits import ecc_omega_from_ecosw_esinw
>>> ecc, arg_peri = ecc_omega_from_ecosw_esinw(Q(0.6, ""), Q(0.8, ""))
>>> float(ecc.value)
1.0
>>> arg_peri.unit
Unit("rad")
Parameters:

  • ecosw (Real[Quantity[PhysicalType('dimensionless')], '*batch'])

  • esinw (Real[Quantity[PhysicalType('dimensionless')], '*batch'])

Return type:

tuple[Real[Quantity[PhysicalType('dimensionless')], '*batch'], Real[Quantity[PhysicalType('angle')], '*batch']]

harv.kepler.orbits.ecosw_esinw_from_ecc_omega(eccentricity, arg_peri)

Convert eccentricity and arg. of periastron to (e cos omega, e sin omega).

The inverse of ecc_omega_from_ecosw_esinw():

\[\begin{split}\mathrm{ecosw} &= e \cos\omega \\ \mathrm{esinw} &= e \sin\omega\end{split}\]

Examples

>>> from unxt import Q
>>> from harv.kepler.orbits import ecosw_esinw_from_ecc_omega
>>> ecosw, esinw = ecosw_esinw_from_ecc_omega(Q(0.5, ""), Q(0.0, "rad"))
>>> float(ecosw.value)
0.5
>>> float(esinw.value)
0.0
Parameters:

  • eccentricity (Real[Quantity[PhysicalType('dimensionless')], '*batch'])

  • arg_peri (Real[Quantity[PhysicalType('angle')], '*batch'])

Return type:

tuple[Real[Quantity[PhysicalType('dimensionless')], '*batch'], Real[Quantity[PhysicalType('dimensionless')], '*batch']]

harv.kepler.orbits.compute_true_anomaly_components(time, period, eccentricity, t_peri)

Compute true anomaly at given times.

Parameters:

  • time (Real[Quantity[PhysicalType('time')], '*batch']) – Observation times, shape (n,)

  • period (Real[Quantity[PhysicalType('time')], '']) – Orbital period

  • eccentricity (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '']) – Orbital eccentricity

  • t_peri (Real[Quantity[PhysicalType('time')], '']) – Time of pericenter passage

Return type:

tuple[Float[Array, '*batch'] | Real[Quantity[PhysicalType('dimensionless')], '*batch'], Float[Array, '*batch'] | Real[Quantity[PhysicalType('dimensionless')], '*batch']]

Returns:

(sin_f, cos_f) — true anomaly components, each shape (n,).

Examples

>>> from unxt import Q
>>> from harv.kepler.orbits import compute_true_anomaly_components
>>> sin_f, cos_f = compute_true_anomaly_components(
...     time=Q([0.0, 25.0, 50.0], "day"),
...     period=Q(100.0, "day"),
...     eccentricity=0.3,
...     t_peri=Q(0.0, "day"),
... )
harv.kepler.orbits.rv_at_times(times, period, eccentricity, t_peri, arg_peri, rv_semiamp, v_sys)

Compute the RV model: K*[cos(omega + f(t)) + e*cos(omega)] + v_0.

Parameters:

  • times (Real[Quantity[PhysicalType('time')], '*batch']) – Observation times.

  • period (Real[Quantity[PhysicalType('time')], '']) – Orbital period.

  • eccentricity (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '']) – Orbital eccentricity.

  • t_peri (Real[Quantity[PhysicalType('time')], '']) – Time of periastron passage. In the likelihood layer this is derived from the dimensionless phase_peri as t_peri = phase_peri * period (see _solve_kepler).

  • arg_peri (Real[Quantity[PhysicalType('angle')], '']) – Argument of periastron omega.

  • rv_semiamp (Real[Quantity[PhysicalType({'speed', 'velocity'})], '']) – RV semi-amplitude.

  • v_sys (Real[Quantity[PhysicalType({'speed', 'velocity'})], '']) – Systemic velocity.

Return type:

Real[Quantity[PhysicalType({'speed', 'velocity'})], '*batch']

Returns:

Radial velocities in the same unit as rv_semiamp and v_sys.

Examples

>>> from unxt import Q
>>> from harv.kepler.orbits import rv_at_times
>>> times = Q([0.0, 50.0, 100.0], "day")
>>> rv = rv_at_times(
...     times,
...     period=Q(200.0, "day"),
...     eccentricity=0.3,
...     t_peri=Q(50.0, "day"),
...     arg_peri=Q(1.2, "rad"),
...     rv_semiamp=Q(8.0, "km/s"),
...     v_sys=Q(-5.0, "km/s"),
... )
>>> rv.unit
Unit("km / s")
harv.kepler.orbits.astrometric_orbit_at_times(times, period, eccentricity, t_peri, arg_peri, cos_i, lon_asc_node, semi_major_axis)

Compute sky-plane astrometric orbit (Deltara, Deltadec) at given times.

Uses the Thiele-Innes parameterization following the Gaia local plane coordinate (LPC) convention (Lindegren & Bastian, GAIA-C3-TN-LU-LL-061-08, Eq. 4):

Deltara  = (B*cos f + G*sin f) * a      (RA / ``a`` direction)
Deltadec = (A*cos f + F*sin f) * a      (Dec / ``d`` direction)

where A, B, F, G are the unit Thiele-Innes constants and a is the photocentric semi-major axis.

Parameters:

  • times (Real[Quantity[PhysicalType('time')], '*batch']) – Observation times.

  • period (Real[Quantity[PhysicalType('time')], '']) – Orbital period.

  • eccentricity (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '']) – Orbital eccentricity.

  • t_peri (Real[Quantity[PhysicalType('time')], '']) – Time of periastron passage. In the likelihood layer this is derived from the dimensionless phase_peri as t_peri = phase_peri * period (see _solve_kepler).

  • arg_peri (Real[Quantity[PhysicalType('angle')], '']) – Argument of periastron omega.

  • cos_i (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], '']) – Cosine of orbital inclination.

  • lon_asc_node (Real[Quantity[PhysicalType('angle')], '']) – Longitude of the ascending node Omega.

  • semi_major_axis (Real[Quantity[PhysicalType('angle')], '']) – Photocentric semi-major axis.

Return type:

tuple[Real[Quantity[PhysicalType('angle')], '*batch'], Real[Quantity[PhysicalType('angle')], '*batch']]

Returns:

(delta_ra, delta_dec) — sky-plane offsets in the same unit as semi_major_axis.

Examples

>>> from unxt import Q
>>> from harv.kepler.orbits import astrometric_orbit_at_times
>>> times = Q([0.0, 100.0, 200.0], "day")
>>> dra, ddec = astrometric_orbit_at_times(
...     times,
...     period=Q(300.0, "day"),
...     eccentricity=0.3,
...     t_peri=Q(0.0, "day"),
...     arg_peri=Q(1.2, "rad"),
...     cos_i=0.5,
...     lon_asc_node=Q(0.8, "rad"),
...     semi_major_axis=Q(3.0, "mas"),
... )
>>> dra.unit
Unit("mas")

harv.kepler.orientation module

Keplerian orbit implementation with units support and JAX compatibility.

class harv.kepler.orientation.KeplerianOrientation

Bases: Module

Orientation of a Keplerian orbit in 3D space.

Stores the three Euler angles that define how the orbital plane is oriented relative to the observer’s reference frame: - Inclination (i): tilt of orbital plane from sky plane - Longitude of ascending node (Omega): where orbit crosses sky plane - Argument of pericenter (omega): orientation of ellipse within orbital plane

Angles are stored as sin/cos pairs for numerical stability.

Examples

>>> from unxt import Q
>>> from harv.kepler.orientation import KeplerianOrientation
>>> orient = KeplerianOrientation.from_angles(
...     arg_peri=Q(0.5, "rad"),
...     lon_asc_node=Q(1.0, "rad"),
...     inclination=Q(0.3, "rad"),
... )
Parameters:

  • sin_arg_peri (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • cos_arg_peri (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • sin_lon_asc_node (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • cos_lon_asc_node (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • sin_i (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • cos_i (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

sin_arg_peri: Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''] = 0.0
cos_arg_peri: Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''] = 1.0
sin_lon_asc_node: Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''] = 0.0
cos_lon_asc_node: Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''] = 1.0
sin_i: Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''] = 0.0
cos_i: Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''] = 1.0
classmethod from_angles(arg_peri=Quantity(Array(0., dtype=float32, weak_type=True), unit='rad'), lon_asc_node=Quantity(Array(0., dtype=float32, weak_type=True), unit='rad'), inclination=Quantity(Array(0., dtype=float32, weak_type=True), unit='rad'))

Construct from angle values.

Parameters:

  • arg_peri (Real[Quantity[PhysicalType('angle')], ''])

  • lon_asc_node (Real[Quantity[PhysicalType('angle')], ''])

  • inclination (Real[Quantity[PhysicalType('angle')], ''])

Return type:

KeplerianOrientation

classmethod from_thiele_innes(A, B, F, G)

Construct orientation from Thiele-Innes constants.

Inverts the Thiele-Innes constants to recover (omega, Omega, i, a). This loosely follows Appendix A of https://arxiv.org/abs/2206.05726 but my convention for angles is different. For this implementation:

0 < i < pi/2 0 < omega < 2pi 0 < Omega < 2pi

whereas the paper linked above assumes:

0 < i < pi 0 < omega < 2pi 0 < Omega < pi

Return type:

tuple[KeplerianOrientation, Real[Quantity[PhysicalType('length')], ''] | Real[Quantity[PhysicalType('angle')], '']]

Returns:

  • orientation – KeplerianOrientation object

  • semi_major_axis – Recovered semi-major axis

Parameters:

  • A (Real[Quantity[PhysicalType('length')], ''] | Real[Quantity[PhysicalType('angle')], ''])

  • B (Real[Quantity[PhysicalType('length')], ''] | Real[Quantity[PhysicalType('angle')], ''])

  • F (Real[Quantity[PhysicalType('length')], ''] | Real[Quantity[PhysicalType('angle')], ''])

  • G (Real[Quantity[PhysicalType('length')], ''] | Real[Quantity[PhysicalType('angle')], ''])

property arg_peri: Real[Quantity[PhysicalType('angle')], '']

Argument of pericenter (omega).

property lon_asc_node: Real[Quantity[PhysicalType('angle')], '']

Longitude of ascending node (Omega).

property inclination: Real[Quantity[PhysicalType('angle')], '']

Inclination (i).

property rotation_matrix: Array

Compute rotation matrix from orbital plane to observer frame.

Returns the rotation matrix R such that: r_observer_frame = R @ r_orbital_frame

The rotation is composed of three sequential rotations: 1. R_z(omega): Rotate by argument of pericenter, omega, in orbital plane 2. R_x(i): Rotate by inclination, i, to tilt orbital plane 3. R_z(Omega): Rotate by longitude of ascending node, Omega, on sky plane

The full rotation matrix is therefore: R = R_z(Omega) @ R_x(i) @ R_z(omega)

We build the matrix directly from the sin/cos pairs for numerical stability and speed, but using the notation below, it is equivalent to:

R1 = jnp.array([[c_w, -s_w, 0], [s_w, c_w, 0], [0, 0, 1]]) R2 = jnp.array([[1., 0, 0], [0, c_i, -s_i], [0, s_i, c_i]]) R3 = jnp.array([[c_W, -s_W, 0], [s_W, c_W, 0], [0, 0, 1]]) R = R3 @ R2 @ R1

Or, alternatively: omega = arg_peri.to_value(“rad”) Omega = lon_asc_node.to_value(“rad”) i = inclination.to_value(“rad”) R = Rotation.from_euler(‘ZXZ’, [Omega, i, omega]).as_matrix()

thiele_innes_constants(semi_major_axis=None)

Compute Thiele-Innes constants (A, B, F, G).

These constants linearize the relationship between orbital position and sky-plane projection. See Appendix A of: https://arxiv.org/abs/2206.05726

Parameters:

semi_major_axis (Real[Quantity[PhysicalType('length')], ''] | Real[Quantity[PhysicalType('angle')], ''] | None) – Semi-major axis of the orbit

Returns:

The four Thiele-Innes constants.

Return type:

A, B, F, G

__init__(sin_arg_peri=0.0, cos_arg_peri=1.0, sin_lon_asc_node=0.0, cos_lon_asc_node=1.0, sin_i=0.0, cos_i=1.0)
Parameters:

  • sin_arg_peri (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • cos_arg_peri (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • sin_lon_asc_node (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • cos_lon_asc_node (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • sin_i (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

  • cos_i (Float[Array, ''] | floating[Any] | float | int | Real[Quantity[PhysicalType('dimensionless')], ''])

Return type:

None

static __new__(cls, *args, **kwargs)
Parameters:
Return type:

TypeVar(_ModuleT, bound= Module)