Orbit class overview

Orbit class overview#

Orbit is the front door to MeepMeep. It is built for the workflow that dominates exoplanet light-curve and radial-velocity work: instantiate the class once with the time array of your observations, and then — inside a modelling or fitting loop — just update the orbital parameters and read out the quantities you need. The expensive per-orbit setup happens once; each iteration costs one small Taylor solve plus the requested observable.

The available observables are planet position, velocity, sky-projected separation, phase angle, radial velocity, reflected-light and thermal phase curves, ellipsoidal variation, and light travel time. With derivatives=True every observable also returns analytic gradients w.r.t. the seven orbital parameters and any method-specific extras, which is what gradient-based optimisers and HMC samplers want.

Quickstart #

import numpy as np
from meepmeep import Orbit

o = Orbit(npt=15)
o.set_pars(tc=0.0, p=3.0, a=8.5, i=np.radians(89.0),
           e=0.1, w=np.radians(90.0))
o.set_data(np.linspace(-0.05, 0.05, 1001))

x, y, z = o.xyz()                       # planet position
d = o.star_planet_distance()            # 3D separation

For analytic gradients in addition to values, switch to derivative mode at construction time and unpack the extra arrays:

from meepmeep import Orbit

o = Orbit(npt=15, derivatives=True)
o.set_pars(tc=0.0, p=3.0, a=8.5, i=np.radians(89.0),
           e=0.1, w=np.radians(90.0))
o.set_data(np.linspace(-0.05, 0.05, 1001))

x, y, z, dx, dy, dz = o.xyz()           # values plus (N, 6) gradients

The shape contract is the same for every observable — see Derivative mode below.

Typical modelling loop #

The class is designed so that the same instance can be reused across many parameter updates without rebuilding internal state from scratch. A transit-fitting or RV-fitting loop typically looks like this:

import numpy as np
from meepmeep import Orbit

# One-time setup: expansion-point grid (constant) and the observation times.
o = Orbit(npt=15)
o.set_data(times)                                  # observation times

def log_likelihood(theta):
    tc, p, a, i, e, w, k = theta
    o.set_pars(tc=tc, p=p, a=a, i=i, e=e, w=w)     # update the orbit
    rv_model = o.radial_velocity(k=k)              # read out observable
    return -0.5 * np.sum(((rv_obs - rv_model) / rv_err)**2)

# ... drop log_likelihood into your sampler / optimiser ...

Only the very first call after construction pays the numba JIT-compile cost; every subsequent iteration is a fast Taylor solve followed by a Horner-scheme evaluation per requested observable.

Construction #

Orbit is constructed with three optional arguments:

Argument	Meaning
`npt`	Number of expansion points used by the multi-expansion-point Taylor expansion (default 15). Includes the periodic-image slot. Raise this when the orbit is eccentric enough that 15 expansion points no longer cover periastron with adequate per-expansion-point accuracy.
`ep_placement`	Expansion-point placement strategy: `'mm'` (uniform in mean motion), `'ea'` (uniform in eccentric anomaly; default and preferred for eccentric orbits), or `'ta'` (uniform in true anomaly). See Multi-expansion-point orbit-spanning evaluation for the placement strategies’ implications.
`derivatives`	If `True`, every observable also returns analytic parameter derivatives. Costs roughly 2–4× the value-only runtime.

The expansion-point grid is built once at construction time; subsequent set_pars() calls reuse it.

Binding orbital elements #

set_pars() accepts the orbital elements as keyword-only arguments. The time anchor is specified by exactly one of tc (transit center) or tp (periastron passage); the two are related by

\[t_p \;=\; t_c \;-\; \frac{M_\mathrm{tr}(e, w)}{2\pi}\,p,\]

where \(M_\mathrm{tr}\) is the mean anomaly at transit. Whichever you pass, the other is derived and stored, so both anchors are available afterwards regardless of which form you used.

# Same orbit, two equivalent calls:
o.set_pars(tc=0.0,   p=3.0, a=8.5, i=np.radians(89.0), e=0.1, w=np.radians(90.0))
# ... or ...
o.set_pars(tp=-0.4786815, p=3.0, a=8.5, i=np.radians(89.0), e=0.1, w=np.radians(90.0))

# After either call:
o._tc   # 0.0
o._tp   # -0.4786815...

Passing both tc and tp, or neither, raises TypeError — the call site is always explicit about which convention it uses. The remaining elements (p, a, i, e, w) are keyword-only too; this keeps the call sites self-documenting and matches the rest of the package.

Binding times #

set_data() binds a 1-D time array to the instance. The bound grid is used as the default times argument by every evaluator that exposes a times=None parameter (xyz(), star_planet_distance(), lambert_phase_curve(), emission_phase_curve(), ellipsoidal_variation()), and used unconditionally by methods that do not expose times at all (vxyz(), cos_phase(), phase(), theta(), radial_velocity(), light_travel_time()).

You can rebind the grid as often as you like without recomputing the Taylor coefficients — set_data() only stores the array.

Observables #

The class groups its observables into a small set of conceptual families. Per-method detail (parameters, return shapes, units, edge cases) lives in the docstrings, surfaced on the API page; this section is a tour.

Geometry. xyz() returns the planet position in the sky frame; vxyz() returns the velocity; star_planet_distance() returns the 3D separation \(r = \sqrt{x^2 + y^2 + z^2}\). cos_phase() returns \(\cos\alpha = -z/r\) (the cosine of the planet-observer-star phase angle), with phase() and theta() returning \(\alpha\) and \(\pi-\alpha\) respectively. The latter two clamp slightly inside the \(\arccos\) domain in derivative mode so the gradient stays finite at the conjunction extrema; see the phase() docstring for the caveat.

Photometry. lambert_phase_curve() evaluates the reflected-light phase curve assuming Lambertian scattering. emission_phase_curve() evaluates a simple cosine thermal-emission phase curve with a hotspot offset (amplitude \(k^2 f_\mathrm{ratio}\)). ellipsoidal_variation() evaluates the stellar ellipsoidal-distortion signal of Lillo-Box et al. (2014).

Spectroscopy. radial_velocity() scales the line-of-sight planet velocity by the closed-form factor \(K / [(2\pi/p)\,(a\sin i)/\sqrt{1-e^2}]\) so the result is the stellar RV signal in m/s.

Timing. mean_anomaly() and true_anomaly() return the corresponding anomalies at the bound times. true_anomaly(exact=True) switches to the Newton-Raphson reference solver for validation runs; the exact path does not provide parameter derivatives, so it raises NotImplementedError if combined with derivatives=True. light_travel_time() returns the LTT correction referenced to primary transit, with the rstar derivative intentionally omitted from the gradient.

Diagnostics. plot() renders a three-panel sky-frame plot of the orbit, with an optional show_exact=True overlay against the Newton-Raphson reference. The private helpers _xyz_error and _cos_phase_error are available for quick sanity checks against the same reference.

Derivative mode #

Construct with derivatives=True to switch every observable into gradient-returning form. The shape contract is uniform:

Multi-coordinate value returns are extended with one gradient array per coordinate. For example xyz() returns (xs, ys, zs, dxs, dys, dzs) instead of (xs, ys, zs), with each d*s of shape (N, ndp).
Single-value returns become two-tuples (value, dvalue) with dvalue of shape (N, ndp).

The trailing axis ndp is the number of differentiable parameters. The first seven slots are always the orbital block (tc, p, a, i, e, w, lan) — with tp replacing tc when the orbit is bound by periastron time (see Convention bridge below); some methods append physical extras:

Method	Trailing extras (in order)
`radial_velocity()`	`k`
`lambert_phase_curve()`	`ag`, `k`
`emission_phase_curve()`	`k`, `fratio`, `offset`
`ellipsoidal_variation()`	`alpha`, `mass_ratio`, `inc`

The rstar derivative of light_travel_time() is intentionally not returned — only the 7 orbital derivatives.

For the gradient math (Kepler implicit-differentiation step, orbital-plane derivative chain, evaluator-level chain rules) see Analytic parameter derivatives.

Convention bridge #

This class accepts either tc (transit-center time) or tp (periastron-passage time) via set_pars(). The underlying Taylor backend anchors its expansion-point grid at periastron, so set_pars() converts once and stores both values internally: self._tc and self._tp. Every Taylor-backend dispatcher call uses self._tp; self._tc is used only by the Newton-Raphson diagnostic paths and by plot(show_exact=True)().

If you ever drop down to the Taylor backend directly, note that its multi-expansion-point dispatchers (pos_o() and friends) take a tpa argument that is the periastron-anchored time — not the transit-center time. See Taylor-series backend overview for the low-level convention and Two ways to use the Taylor backend for when to drop down at all.

Orbit class overview

Contents

Orbit class overview#