API Documentation

Body(; mass, flux=(;), name=:auto)

A point mass participating in a System. mass is in solar masses (see the msun, mjup, mearth constants). flux is an optional NamedTuple of per-band fluxes (arbitrary but consistent units within a band), used to compute photocentres.

Body values are construction-time inputs, typically rebuilt every MCMC sample. A named Body doubles as a handle: observables resolve it by name against the solved system (raoff(sol, b, A) with the same variables used to build it), reading only its name — a value left over from a previous sample resolves identically. Unnamed bodies (and hot loops) use the persistent references from bodies(sys) instead.

Orbit(exterior; about, <size>, e=0, i=0, ω=0, Ω=0, tp=0, M=nothing)

One Keplerian relationship in a system: the orbit of exterior about about. Both endpoints take the same grammar — a Body, or a tuple of Bodys denoting their barycentre:

Orbit(b, about=A;      a=5.2, e=0.05)   # astrocentric
Orbit(c, about=(A, b); a=9.6)           # Jacobi
Orbit((Ba, Bb), about=(Aa, Ab); a=120)  # 2+2 quadruple

The convention is always explicit: there is no default about= and no inference from semi-major axis. See Jacobi and Astrocentric for the two standard chains.

Elements

Supply exactly one alternative from each group. Orientation is always i [rad] and Ω [rad]; every angle is in radians and every epoch is MJD.

secosω = √e·cosω and ecosω = e·cosω sample the eccentricity disc rather than the half-plane, removing the ω degeneracy at e → 0. M0 is the mean anomaly at epoch [rad]; θ is the sky-plane position angle at epoch [rad]. τ is deliberately not accepted — it needs hidden period and reference-epoch state and has no clean meaning under N-body integration.

P is in days, matching period(sys) so the two round-trip. Users who think in years get a plausible-looking 365× error rather than a crash, so show prints the period in both units.

Cartesian initial conditions

Alternatively give the full relative state, which determines every element and so replaces all of the groups above:

Orbit(b, about=A; x=1.2, y=0.3, z=-0.1,       # [AU]
                  vx=-0.9, vy=2.1, vz=0.2,    # [AU / julian yr]
                  epoch=59000.0)              # [MJD]

Positions and velocities are of exterior relative to about, in the same frame and units as posx/velx, so a state read back out of a solution reconstructs the same orbit. Unbound states are fine — a comes out negative and e > 1 — which makes this the natural way to specify hyperbolic orbits.

Eccentricity

e > 1 is supported (unbound); a is negative by convention there and a positive value is taken as |a|. e == 1 exactly is rejected: the elements are degenerate for parabolae. Use Cartesian initial conditions instead.

M (compatibility escape hatch)

The row's gravitating mass is normally the total mass of every body the row binds, taken from the Body values themselves. Passing M [M⊙] overrides it. This is physically inconsistent — it decouples Kepler's third law from the masses that drive the reflex amplitudes — and exists only to reproduce published fits bit-for-bit and to match orbitize!/RadVel conventions. It is a compatibility switch, not a modelling choice.

System(bodies, orbits; plx=…, ra=…, dec=…, pmra=…, pmdec=…, rv=…, ref_epoch=…)
System(orbits; …)

A hierarchical system. bodies is a tuple of Body values defining the body order (and hence the indices of bodies(sys)); orbits is a tuple of Orbits, one per degree of freedom — exactly length(bodies) - 1 of them. Given only orbits, the bodies are collected in order of first appearance.

A single body with no orbits is legal: System((A,), ()). The lone body is the system barycentre, so its motion is purely the frame's (position, proper motion, parallax, and the observing-geometry corrections); all displacement observables of the body relative to the barycentre are identically zero. This is the natural model for fitting the absolute astrometry of an isolated star.

A = Body(mass=1.1, name=:A)
b = Body(mass=8mjup, name=:b)
c = Body(mass=2mjup, name=:c)
System((A, b, c), (Orbit(b, about=A;      a=2.5),
                   Orbit(c, about=(A, b); a=8.0)))

The topology is whatever the orbits say — Jacobi chains, astrocentric sets, moons, 2+2 quadruples, and mixtures of these are all expressible, and the convention is never inferred. See Jacobi and Astrocentric for the two standard chains, and show(sys) for the convention actually resolved.

The keyword arguments define the system's frame, attached to the system barycentre:

• none → physical-unit observables only,
• plx [mas] → angular observables in mas,
• all of them → full absolute frame (see AbsoluteFrame).

System values are cheap, immutable, isbits structures rebuilt every sample. Use bodies(sys) for the persistent per-body references.

reframe(sys, frame)
reframe(sys; plx=…, ra=…, dec=…, pmra=…, pmdec=…, rv=…, ref_epoch=…)

sys with its frame replaced and nothing else changed. The keyword form takes exactly the frame arguments System does, so reframe(sys; plx=p) and reframe(sys) (→ NoFrame) are spelled the same way as construction.

Masses, rows, row specs, A⁻¹ and fluxes are all frame-independent by construction — the frame enters only in frame_pass!/observe_pass!, never in the propagation — so they are carried over verbatim rather than rebuilt. That is what makes this exact rather than merely equivalent: no element is recomputed, so there is no room for a recomputation to differ in the last bit.

The motivating use is a two-stage frame: the barycentric AbsoluteFrame inputs are not known until the model's own bodies have been solved, but body–barycentre kinematics are fully defined at Parallax level (they need a distance, not a sky position). So a system can be built and solved at Parallax level, that solution used to derive the absolute frame, and the frame then installed here. Octofitter's anchored-frame parameterization is exactly this: sample the anchor source's observed catalog quantities, and subtract the model's own motion of the anchor body about the system barycentre to get the barycentric ones.

The scalar type only ever widens: promote_type(T, …) of the current type and the new frame's. A frame narrower than the system (a Float64 frame onto a Dual system) leaves the system Dual, where fresh construction would have given Float64 — same values, different type. Widening cannot lose a derivative; narrowing would.

See also System.

OrbitDomainError(msg)

The parameter values — not the system's structure — put an orbit outside the domain where it is defined: a non-finite mass, e == 1 exactly, a semi-major axis that does not size the conic (a ≤ 0 or a non-finite for an ellipse, a == 0 or non-finite for a hyperbola) or a period that is not positive and finite, a radial or zero-separation Cartesian state.

This is deliberately a distinct type from the ErrorExceptions the structural checks throw, and the distinction is the useful one for a sampler. A structural problem (a malformed hierarchy, a body with no orbit) is a property of the model: it is true of every draw, so it deserves to be loud and to stop the run. A domain problem is a property of one proposal: a sampler exploring a wide prior — a parallel-tempering scheme's hottest chains especially — will reach parameters that overflow to Inf under the unconstrained-space transform, and the right response is to score that proposal -Inf and move on, not to warn.

So callers that build a system per likelihood evaluation should catch OrbitDomainError quietly and let anything else through.

Jacobi(inner, b => (; a=…, e=…), c => (; a=…), …)

Build a Jacobi chain: each body orbits the barycentre of everything interior to it. Equivalent to spelling out

Orbit(b, about=inner;       …)
Orbit(c, about=(inner…, b); …)

Returns a tuple of Orbits for System. inner is a Body or a tuple of Bodys (e.g. a tight binary at the centre of a circumbinary chain).

Astrocentric(centre, b => (; a=…), c => (; a=…), …)

Build an astrocentric set: every body orbits centre directly. Equivalent to Orbit(b, about=centre; …), Orbit(c, about=centre; …), ….

Note this is a materially different model from Jacobi under KeplerianApprox (not a relabelling): the rows are the approximation. Under AHL21 the two agree, since rows only set initial conditions.

ThieleInnes(; A, B, F, G, plx=nothing)

Convert Thiele-Innes constants to the size and orientation elements, as a NamedTuple to splat into Orbit:

Orbit(b, about=A; ThieleInnes(A=…, B=…, F=…, G=…, plx=24.5)...,
                  e=0.3, tp=59000.0)

Thiele-Innes replaces a, i, ω and Ω jointly; supply the shape and phase elements as usual. This parametrization has no coordinate singularity at e → 0 or i → 0, which is why Gaia-only astrometric fits use it.

A, B, F, G are in AU, or in mas if plx [mas] is given.

The inverse is thieleinnes(sys[, k]; plx=nothing).

thieleinnes(sys, k=1; plx=nothing)

Thiele-Innes constants (A, B, F, G) of hierarchy row k, in AU — or in mas if plx [mas] is given. Inverse of ThieleInnes.

orbit(; M, a, e=0, i=0, ω=0, Ω=0, tp=0, plx=…, ra=…, …)

Construct a trivial two-body System: a primary of mass M in solar masses, named :A, with a massless secondary named :b on the given orbit. This reproduces the classic PlanetOrbits v0.11 KepOrbit/Visual/AbsoluteVisual semantics, where M is the total mass and one-argument observables (raoff(sol), …) give the secondary relative to the primary.

Frame keywords are passed through to System (none, plx, or the full absolute frame set).

This v0.11-compatibility convenience is deliberately not exported — opt in with using PlanetOrbits: orbit. For a secondary with real mass — and for anything hierarchical — construct System(Orbit(…)) explicitly.

rvorbit(; M, msini=0, a=…|P=…, e=0, ω=0, tp=0)

Construct a two-body System in the radial-velocity-only convention: no parallax (so angular observables are unavailable rather than silently wrong), i = π/2 and Ω = 0.

Radial velocities constrain only m·sin i, never m and i separately, so msini is exactly what the data measure — under the i = π/2 convention the secondary's mass is its m·sin i. Reading it as a true mass is a lower bound.

Like orbit, this is an opt-in convenience: using PlanetOrbits: rvorbit. For a full 3D fit, build the System directly and let the astrometry constrain i.

bodies(sys)

NamedTuple of persistent BodyRefs for the bodies of sys, keyed by the names given at construction (Body(… , name=:b)), in the order the bodies were listed. These are the resolved form of what observables accept — guaranteed cheap in hot loops, and the only handles to unnamed bodies. (Named Body values and Symbols resolve to them automatically.)

(; A, b) = bodies(sys)
raoff(sol, b, A)

barycentre(sys)
barycentre(sys, members...)

The mass-weighted barycentre of the whole system, or of the subsystem spanned by the given members, as a WeightedPoint. Members can be given as BodyRefs, named Body values, or Symbols, e.g. barycentre(sys, jup, gan).

photocentre(sys; band=nothing)
photocentre(sys, members...; band=nothing)
photocentre(weights::StaticVector)

The flux-weighted photocentre of the system — the point astrometric instruments observe for blended sources — as a WeightedPoint. With more than one band defined on the system's bodies, pass band to select one.

Given members, the photocentre is over that subset only: weights f_j / Σ_members f_k for members and zero for everything else. Members can be given as BodyRefs, named Body values, or Symbols, exactly as for barycentre:

photocentre(sys, Aa, Ab; band=:G)

A subset whose total flux is zero is an error: a structural membership declaration over bodies that are all dark has no meaning. Membership that varies per draw or per epoch is the likelihood's business — build the weight vector there and pass it to the third method, which normalizes it:

photocentre(fluxes(sys, :G) .* member)

WeightedPoints are isbits, so constructing one per epoch is free.

fluxes(sys)
fluxes(sys, band)

Per-band fluxes of sys's bodies, as they were declared on the Body values: a NamedTuple of SVector{NB}s keyed by band, or the SVector{NB} for one band. Bodies that declared no flux in a band read zero.

This is the entry point for likelihoods that need to build their own weights — a per-epoch resolution taper, a sampled membership indicator — and hand the result to photocentre:

f = fluxes(sys, :G)
wp = photocentre(f .* member)     # a WeightedPoint, normalized

Vector order matches bodies(sys).

framedirection

The system barycentre's apparent direction at each epoch — where frame_ra and frame_dec point, and the origin of the tangent plane every angular observable is measured in. Accepted anywhere a reference is.

It is a direction, not a place, and the name says so because that is the entire distinction it exists to make. Against the zero-argument observables it is nearly the same reference as barycentre(sys) — they differ only by the barycentre's own retardation residual, far below a µas — so the two look interchangeable right up until you pass an observer position, at which point they differ by the whole parallax ellipse:

raoff(sol, A, barycentre(sys), obs_pos)  # relative: the differential part
raoff(sol, A, framedirection, obs_pos)   # absolute: parallax + orbit

A barycentre sits at a finite distance, so displacing the observer shifts it and the target together and only the differential (Kopeikin) term survives. A direction has no parallax of its own, so the target keeps its parallax factor in full. Relative versus absolute astrometry is therefore expressed by which reference you name, not by a flag — and both are designed uses, so neither is refused.

Nothing emits from a direction, so it carries no radvel Einstein term, and it does not move within the tangent plane: the frame's own drift is frame_pmra/frame_pmdec, not a pairwise observable.

BodyRef

Persistent, isbits reference to a body of a System (an index into its body list). Obtained from bodies(sys); valid across samples for any system sharing the same topology. This is the resolved form observables work with — interactively, named Body values and Symbols are accepted anywhere a BodyRef is and resolve to one by name.

WeightedPoint(w::SVector{NB}; emits::Bool)

A normalized weighted combination of body states — the generalization of a single body used to represent barycentres (mass weights) and photocentres (flux weights). Valid only for the sample whose masses/fluxes produced it.

Every observable accepts a WeightedPoint anywhere a reference goes, in either the target or the reference position, and evaluates it as one dot product over absolute body states. That is what makes a blended source containing a hierarchical pair correct by construction: the point carries the intra-pair wobble and the pair's wide-orbit motion at once, with no per-level bookkeeping and under either propagator.

Normally obtained from barycentre or photocentre, which carry the semantics in their names. A likelihood whose weights are its own business — a per-epoch resolution taper, a sampled membership indicator — may build one directly; the weights are expected to sum to one, and photocentre(w) is the normalizing constructor. WeightedPoint is isbits, so constructing one per epoch inside a scan loop costs nothing.

emits has no default, deliberately. It says whether light comes from this point, which changes what radvel means for it and nothing else — so a wrong value is invisible in every other observable, and would be invisible here too if it could be omitted. A photocentre emits; a barycentre does not. Build a photocentre-like point with emits=true and a mass-weighted or geometric one with emits=false.

orbitsolve(sys::System, epochs; method=KeplerianApprox())

Solve sys at the sorted epochs [MJD], returning a Trajectory. Index it (traj[k]) for per-epoch solutions to pass to observables.

Allocates the trajectory storage; for allocation-free hot loops see orbitsolve!.

orbitsolve(sys::System, t::Real; method=KeplerianApprox())

Single-epoch convenience: solve at one epoch [MJD] and return the solution directly. Allocates; batch epochs into a vector for performance.

orbitsolve!(traj::Trajectory, sys::System; method=KeplerianApprox())

Fill a caller-allocated Trajectory with per-body barycentric states of sys at traj.epochs. Performs no allocation itself (unless threads > 1): with caller-provided (e.g. bump-allocated) column storage the whole construct → solve → query path is allocation-free.

The two precision opt-outs are independent, and gate different corrections — see the "Precision opt-outs" page in the manual before setting either.

• observing_geometry=false skips the observing-geometry pass, whose terms all scale with the system's angular extent ρ. See observe_pass!.
• barycentric_lighttime=false skips the barycentric light-travel solve, a whole-system timing correction that scales with proximity and proper motion, not with ρ. See frame_pass!.

threads=n (default 1) splits the epochs into up to n contiguous chunks solved on concurrent tasks. Every pass is epoch-local, each chunk writes a disjoint epoch range of the same storage, and epochs keep their identity — so the result is identical to the serial solve, bit for bit. Only the KeplerianApprox propagator supports this (AHL21 marches through time sequentially); with any other method, or too few epochs for the task overhead to amortize (512 per task), the solve silently runs serial.

Trajectory(sys, epochs)
Trajectory{T}(sys, epochs)

Allocate storage for solving sys at epochs (sorted, MJD). The element type defaults to the system's scalar type; pass T explicitly when solving with a different element type (e.g. ForwardDiff Duals).

For allocation-free hot loops the columns can come from a caller-owned allocator instead; see the Trajectory(alloc, T, sys, epochs) method and orbitsolve!.

Trajectory(alloc, T, sys, epochs)

Build a Trajectory whose columns come from alloc rather than the heap: alloc(S, n) must return an AbstractVector{S} of length n and alloc(S, n, m) an n × m AbstractMatrix{S}.

This is the entry point for callers that own their scratch storage — e.g. Octofitter's per-sample Bumper buffers:

@no_escape begin
    traj = Trajectory((S, dims...) -> @alloc(S, dims...), T, sys, epochs)
    orbitsolve!(traj, sys)
    …
end

The column set is deliberately not part of the public interface (it has grown as passes were added); go through this constructor rather than the inner one so callers do not have to track it.

sol = traj[k]

Immutable zero-cost view of a Trajectory at epoch index k. Observables (raoff(sol, b, A), radvel(sol, A, barycentre(sys)), …) read from it.

soltime(sol)

The observation epoch [MJD] this solution corresponds to — identically the value passed into orbitsolve, preserving the epoch-indexing contract (soltime(traj[k]) === epochs[k]).

KeplerianApprox(; solver=Auto(), simd=true)

Propagator in which every hierarchy row evolves on an independent Keplerian orbit (exact for two bodies; the classic approximation for hierarchical systems). solver selects the Kepler-equation algorithm (see PlanetOrbits.Markley, PlanetOrbits.RootsMethod).

With simd=true (the default), solves with the Markley/Auto solver batch across epochs through branch-free vectorizable kernels (≈4× on AVX2; agrees with the scalar solver to ≤4e-15). This covers both Float64 elements and first-order ForwardDiff Duals, which solve their primal roots through the same kernel and attach partials analytically — so a gradient evaluation carries a value bit-identical to a plain Float64 evaluation. Other solvers, nested Duals (Hessians) and hyperbolic orbits use the scalar path, which applies the same implicit rule one epoch at a time.

AHL21(; h, t0=nothing)

N-body symplectic propagator of Agol, Hernandez & Langford (2021), MNRAS 507, 1582 (arXiv:2106.02188) — please cite this paper when publishing results computed with it. Pass as orbitsolve(sys, epochs; method=AHL21(h=…)).

• h: fixed timestep [days]. Guidance: h ≲ P_min/20 for the shortest period in the system; a warning is emitted (once) beyond that. Tight moons make likelihood evaluations expensive — the cost is one map evaluation per h of timespan covered, there is no adaptivity.
• t0: osculating epoch [MJD] — the orbital elements are interpreted as osculating elements at t0. Defaults to the frame ref_epoch for systems with an absolute frame; must be given explicitly otherwise.

Unlike KeplerianApprox, every body gravitates: hierarchical approximation error disappears, and phenomena it cannot express (TTVs, resonant interactions, moon-planet-star coupling) are captured. For a two-body system the map is exact and the two propagators agree to roundoff.

posx(sol, target, reference)

Position offset [AU] of target relative to reference along the right-ascension (east) direction. target/reference are BodyRefs or WeightedPoints — or named Body values / Symbols, resolved by name.

Resolved on the local East/North/line-of-sight triad of the system barycentre's apparent direction at the observation epoch (not at ref_epoch), with each body taken at its own light-travel-retarded time. See observe.jl.

posy(sol, target, reference)

Position offset [AU] of target relative to reference along the declination (north) direction, on the same triad as posx.

posz(sol, target, reference)

Position offset [AU] of target relative to reference along the line of sight, positive away from the observer, on the same triad as posx.

velx(sol, target, reference)

Velocity [AU / julian year] of target relative to reference along the right-ascension (east) direction, on the same triad as posx.

vely(sol, target, reference)

Velocity [AU / julian year] of target relative to reference along the declination (north) direction, on the same triad as posx.

velz(sol, target, reference)

Line-of-sight velocity [AU / julian year] of target relative to reference — the kinematic quantity, for dynamics.

With the observing-geometry pass on, this already carries the projection onto each body's own apparent direction (see observe.jl); what it does not carry is the relativistic Einstein term. For the quantity a spectrograph reports, in m/s, use radvel. The distinction is kinematic vs. spectroscopic, not coordinate vs. apparent.

radvel(sol, target, reference)

Radial velocity [m/s] of target relative to reference along the line of sight, positive receding — the spectroscopic quantity, i.e. what a spectrograph reports. E.g. radvel(sol, b, A) for a relative RV, or radvel(sol, A, barycentre(sys)) for the stellar reflex.

Two pieces: the kinematic projected line-of-sight velocity (velz, in physical units) and the Einstein term — the second-order Doppler and gravitational-redshift difference between the two references,

Ein_i = ( ½|v_tot,i|² + Σ_{j≠i} G·mⱼ / r_ij ) / c

with v_tot the body's total barycentric velocity (orbital, plus the frame's space velocity when the system has an AbsoluteFrame). Nothing emits from a barycentre, so its Einstein term is zero and the stellar-reflex case carries the star's own in full.

There is no keyword to decline this. The orbit-varying part depends on the sampled orbit (e, masses, r(t)), so no reduction pipeline can have removed it, and the constant part is absorbed by the instrument offset either way. Its size, and which of the two uses of radvel it matters for, are tabulated on the "Precision opt-outs" page — briefly, sub-cm/s for a stellar reflex with a planetary companion, but several m/s of variation for the relative RV of a close-in eccentric one.

Use velz instead when you want the kinematic velocity — for dynamics, or to compare against a coordinate-velocity reference.

raoff(sol, target, reference)

Right-ascension offset [mas] of target relative to reference: the difference of their gnomonic (tangent-plane) coordinates about the system barycentre's apparent direction at the observation epoch, with each body taken at its own light-travel-retarded time.

Note this is not posx * cart2angle for a single shared scale factor: the two references are divided by their own d + z, which differs from the shared-scale answer by ρ² in radians (≈ 4.85·ρ[″]² µas).

raoff(sol, target, reference, obs_pos)
decoff(sol, target, reference, obs_pos)

Right-ascension / declination offset [mas] of target relative to reference as seen from obs_pos, the observer's barycentric position in ICRS Cartesian coordinates [AU] at this epoch (e.g. the Earth's, Gaia's at L2, or (0, 0, 0) for the solar-system barycentre).

This is the seam for annual–orbital parallax (Kopeikin 1995) and exact per-body parallax factors: the apparent direction of each body is computed from the observer's actual position by the same exact geometry the zero-argument forms use from the SSB, so the full coupling falls out with no series expansion.

Which reference you name is what decides relative versus absolute, and there is no flag. A body and a barycentre both sit at the system's distance, so displacing the observer shifts target and reference together: the first-order parallax cancels and only the differential (Kopeikin) part survives — that is relative astrometry. framedirection is a direction, not a place, and has no parallax of its own, so against it the target keeps its parallax factor in full — that is absolute astrometry. Same code path, no special case. (An absolute-astrometry likelihood that instead supplies its own parallax term, as Gaia's published parallax_factor_al does, uses the zero-argument forms against a barycentre and never passes obs_pos at all.)

Conventions: obs_pos is ICRS, in AU, and the epochs passed to orbitsolve are barycentric (BJD_TDB-like MJD). Requires an AbsoluteFrame — an ICRS observer position is meaningless without the target's ICRS direction — and a trajectory solved with observing_geometry=true.

Ephemerides are deliberately not PlanetOrbits' business: the caller supplies the position. A likelihood that needs one and has no ephemeris source should say so when it is constructed, not degrade silently here.

raoff(sol, b, A, earth_pos_au)              # relative: differential part only
raoff(sol, A, framedirection, earth_pos_au) # absolute: the full parallax ellipse
raoff(sol, A, barycentre(sys), (0,0,0)) == raoff(sol, A, barycentre(sys))

decoff(sol, target, reference)

Declination offset [mas] of target relative to reference. See raoff.

pmra(sol, target, reference)

Instantaneous relative proper motion [mas/julian year] of target with respect to reference in right ascension — the exact time derivative of raoff, so it carries the same per-body depth scaling.

pmdec(sol, target, reference)

Instantaneous relative proper motion [mas/julian year] of target with respect to reference in declination — the exact time derivative of decoff. See pmra.

projectedseparation(sol, target, reference, obs_pos)
posangle(sol, target, reference, obs_pos)

Separation [mas] and position angle [rad] as seen from obs_pos. See the four-argument raoff.

projectedseparation(sol, target, reference)

Projected separation [mas] between target and reference.

posangle(sol, target, reference)

Position angle [rad] of target about reference, measured from north through east.

Nearly — but no longer exactly — parallax-free. For a body-vs-body pair the per-body depth factor 1/(d + z) is common to both components of the separation and cancels, so the distance drops out as it always did. When the two references sit at different line-of-sight depths (a body versus a barycentre or photocentre) it does not cancel exactly, and the position angle acquires a dependence on distance at the ~1e-6 rad level. On systems built without plx the physical-unit fallback is used and no parallax is required.

frame_ra(sol)     [deg]

Apparent right ascension of the system barycentre frame at this epoch, from rigorous 3D space-motion propagation of the AbsoluteFrame catalog values.

This and its four siblings — frame_dec, frame_pmra, frame_pmdec, frame_rv — describe the frame, not any one body. Absolute quantities of a body compose as frame plus the pairwise observable taken against barycentre(sys), e.g.

frame_pmra(sol) + pmra(sol, A, barycentre(sys))   # the star's absolute pmra

Requires a system built with a full absolute frame (ra, dec, plx, pmra, pmdec, rv, ref_epoch).

Convention: with barycentric_lighttime=true the angular quantities here (frame_ra/frame_dec/frame_pmra/frame_pmdec) are rigorous apparent quantities — the light-time-affected direction and its rates against the observation epoch, d/dt_obs — from the de-Dopplered true worldline. With barycentric_lighttime=false they are the light-time-free standard model the astrometric catalogs are reduced with. frame_rv is in the spectroscopic convention either way (see its docstring). All five reproduce the catalog values exactly at ref_epoch; away from it the two settings differ only by the genuine second-order light-time terms. Compare catalog-convention data against the light-time-free setting (Butkevich & Lindegren 2014, Sects. 5.5, 6.1).

frame_ra and frame_dec are computed on demand from the solved emission epoch and the trajectory's frame rather than stored per epoch: they are the only frame quantities requiring a transcendental, and nothing inside the solver consumes them. Reading both costs ~33 ns/epoch; not storing them takes frame_pass! from 76 to 40 ns/epoch for every model that never asks.

frame_dec(sol)    [deg]

Apparent declination of the system barycentre frame at this epoch. Computed on demand; see frame_ra.

frame_pmra(sol)   [mas / julian yr]

Apparent proper motion of the system barycentre frame in right ascension at this epoch (already including the cos δ factor). Propagated, not frozen at the catalog value, so perspective acceleration is present. See frame_ra.

frame_pmdec(sol)  [mas / julian yr]

Apparent proper motion of the system barycentre frame in declination at this epoch. See frame_ra.

frame_rv(sol)     [m/s]

Radial velocity of the system barycentre frame at this epoch, positive receding — the propagated frame quantity, not a body's reflex. See frame_ra.

Unlike the angular frame quantities, this one carries no apparent-rate conversion under barycentric_lighttime=true: it is the coordinate radial rate at the emission event, which is what the Doppler shift of light leaving that event corresponds to — not d(distance)/dt_obs, which is an astrometric quantity (the two differ by v_r²/c, 40 m/s at Barnard's star). The astrometric catalogs make the same split, publishing apparent proper motions beside a spectroscopic radial velocity. Consequently frame_rv equals the frame's rv at ref_epoch under either setting of the flag, so code composing it against a catalog radial velocity needs no knowledge of which was used.

PlanetOrbits.Auto()

Automatic choice of Kepler solver algorithm. Currently defaults to PlanetOrbits.Markley()

PlanetOrbits.Markley()

Kepler solver implementation from AstroLib, based on Markley (1995) Celestial Mechanics and Dynamical Astronomy, 63, 101 (DOI:10.1007/BF00691917).

PlanetOrbits.HyperbolicHalley()

Allocation-free Halley solver for the hyperbolic Kepler equation (M = e·sinh(H) − H, e > 1). Selected automatically by Auto() when e > 1; the elliptical solvers do not apply there.

PlanetOrbits.RootsMethod(method::Roots.AbstractUnivariateZeroMethod, kwargs...)

Wraps a root finding method from Roots.jl. You can also pass keyword arguments that will be forwarded to Roots to control the tolerance.

Examples:

method = PlanetOrbits.RootsMethod(Roots.Newton())
method = PlanetOrbits.RootsMethod(Roots.Thukral5B())
method = PlanetOrbits.RootsMethod(Roots.Bisection())
method = PlanetOrbits.RootsMethod(Roots.A42())
method = PlanetOrbits.RootsMethod(Roots.Newton(), rtol=1e-3, atol=1e-3)

period(sys)
period(sys, k)

Period [days] of the system's only orbit, or of hierarchy row k.

totalmass(sys)

Total mass of every body in the system [M⊙].

semimajoraxis(sys)
semimajoraxis(sys, k)

Semi-major axis [AU] of the system's only orbit, or of hierarchy row k. Negative for unbound (hyperbolic) orbits.

eccentricity(sys)
eccentricity(sys, k)

Eccentricity of the system's only orbit, or of hierarchy row k. Greater than 1 for unbound orbits.

inclination(sys)
inclination(sys, k)

Inclination [rad] of the system's only orbit, or of hierarchy row k.

meanmotion(sys)
meanmotion(sys, k)

Mean motion [rad / julian year] of the system's only orbit, or of hierarchy row k.

periastron(sys)
periastron(sys, k)

Epoch of periastron passage tp [MJD] of the system's only orbit, or of hierarchy row k.

distance(sys)

Distance to the system [pc]. Requires a parallax or absolute frame.

mjd("2020-01-01")

The modified julian day of a calendar timestamp, for writing epochs down readably: mjd("2020-01-01") == 58849.0.

The timestamp is read on the TT (Terrestrial Time) scale, which is the scale epochs are wanted on — TT and TDB differ by periodic terms below 2 ms, which no orbit fit can see. The Conventions page of the manual has the full statement.

mjd(Date("2020-01-01"))

The modified julian day of a Date or DateTime, read on the TT scale. Same caveats as the string method above: no timescale conversion is performed.

mjd()

The modified julian day of right now, on the TT scale.

mjd2date(modified_julian)

Get a Date value from a modfied julian day, rounded to closest day

Examples

julia> mjd2date(59160.8)
2020-11-08

years2mjd()

Convert from decimal years (e.g. 1995.25) into modified julian date, rounded to closest second

msun

One IAU solar mass, in the mass unit used throughout PlanetOrbits (which is the solar mass, so msun == 1.0). Provided so masses read naturally: Body(mass=1.2msun).

mjup

One jupiter mass in solar masses (exact ratio of the IAU nominal values of GMjup and GMsun): Body(mass=5.3mjup).

mearth

One earth mass in solar masses (exact ratio of the IAU nominal values of GMearth and GMsun): Body(mass=23mearth).

year2day_julian

Days per Julian year, 365.25 exactly (IAU). Exported because the P orbital element is in days while periods are usually quoted, prior'd and plotted in years, so user code converts constantly:

P = P_yr * year2day_julian     # sampled a period in years

Deliberately not named year2day: this package carries two years that differ by 1.9e-5 — the Julian year here and kepler_year_to_julian_day_conversion_factor (365.2568983840419 days), the period of a 1 M⊙ / 1 AU orbit under the IAU nominal constants. Mixing them produces results that look right and quietly drift. This is the one you want for a period expressed in years; the other only appears inside Kepler's third law.

day2year_julian

Julian years per day, 1/year2day_julian. See year2day_julian.

PlanetOrbits.kepler_year_to_julian_day_conversion_factor

The "Kepler year" in julian days, 365.2568983840419 — the period of a hypothetical planet orbiting a 1 M⊙ star at exactly 1 AU under the IAU nominal constants. This is what appears inside Kepler's third law:

period_days = √(a³ / M) * kepler_year_to_julian_day_conversion_factor

It exists because the IAU definitions of the AU and the solar mass do not combine to give an orbital period of exactly one julian year. It therefore differs from year2day_julian (365.25) by 1.9e-5 — small enough that substituting one for the other looks right and quietly drifts. Use this one only in Kepler's third law; use year2day_julian for a period quoted in years.

Unexported, and deliberately verbose, so that the distinction has to be made on purpose.

Derivation, from Gilles Otten (thank you for tracking this down): with G*M_sun = 1.3271244e20 m³/s² and 1 AU = 149_597_870_700 m by definition, sqrt(4π²/1.3271244e20 * 149597870700³) / 86400 = 365.2568983840419 julian days.

plotinfo(f) -> (; label, unit, flip, wrap)

Axis metadata for an observable function f (e.g. raoff, radvel):

• label — human-readable quantity name.
• unit — unit string, or "" for dimensionless.
• flip — whether the axis should increase leftward/downward (RA-like axes: east is to the left on the sky).
• wrap — the period the quantity wraps with (2π for position angle), or nothing.

plotlabel(f) renders these as "label [unit]".

plotlabel(f) -> String

Axis label "label [unit]" for an observable function, from plotinfo.

paraminfo(name::Symbol) -> (; label, unit, angle) | nothing

Display metadata for a parameter by its conventional name — the orbital element keywords Orbit accepts, plus frame variables and common fit parameters. angle marks radian-valued parameters conventionally displayed in degrees. Corner plots and summary tables share this table instead of each keeping a private label dictionary; nothing means "not a name this table knows".

orbit_track_epochs(sys, k=only; n=150, tstart=nothing) -> Vector{Float64}

Epochs [MJD] tracing one full cycle of hierarchy row k, spaced uniformly in eccentric anomaly — dense at periastron, sparse at apastron — so the plotted track is smooth with few points. The first and last epoch describe the same orbital phase one period apart, closing the track.

The cycle drawn is the one containing tstart (default: the cycle starting at the row's periastron epoch). Under KeplerianApprox every cycle is identical; under time-dependent effects (absolute frames, N-body) the choice of cycle matters, so pass tstart near your data.

Errors for hyperbolic rows — an unbound orbit has no period; sample a time span instead.

plot_epochs(sys, tmin, tmax; points_per_period=30, min_points=200, max_points=1000)

Epochs [MJD] for plotting model curves over the span [tmin, tmax], allocated per hierarchy row: each bound row contributes points at period/points_per_period spacing over the whole span, so a short-period inner orbit stays smooth without inflating the grid for a long-period outer one. A uniform base grid of min_points guarantees coverage; the union is clamped to max_points by thinning the densest contributions first (the uniform base and the data span are never dropped).

This replaces sampling the whole span at the minimum period over all orbits, which over-resolves every wide orbit in a system with a broad period spread.

orbit_phase(row_or_sys[, k], t) -> Float64

Mean-anomaly phase in [0, 2π) of hierarchy row k at epoch t [MJD], computed from the row's elements (2π·(t − tp)/P, wrapped). This is what plots colour orbit tracks by; it is exact under KeplerianApprox and the osculating phase under AHL21. Returns NaN for hyperbolic rows.

orbitlines!(ax, sys, target, ref; n=150, kwargs...)

Draw the closed orbit track(s) of target relative to ref on ax, sampled uniformly in eccentric anomaly and coloured by orbital phase. Requires Makie to be loaded.

orbit_theme()

A Makie theme carrying the PlanetOrbits/Octofitter plot look: no axis grid lines, and no colorbar tick marks. Deliberately minimal — everything else is left at Makie's defaults so that a user theme composes with it rather than being overridden. Requires Makie to be loaded.

with_theme(orbit_theme()) do
    octoplot(model, chain)
end

MJDConversion()

A Makie dimension conversion for epoch axes: data are Modified Julian Days (plain Float64), tick labels are calendar dates that adapt as you zoom. Date/DateTime values plot into the same axis and convert automatically (for single-argument plots like vlines!, which bypass Makie's dim conversions, use mjd):

ax = Axis(fig[1,1]; dim1_conversion=PlanetOrbits.MJDConversion())
lines!(ax, epochs_mjd, rvs)
scatter!(ax, [DateTime("2035-06-01")], [0.0])
vlines!(ax, mjd("2035-06-01"))

Requires Makie to be loaded (the method is added by the PlanetOrbitsMakie extension).

add_mjd_axis!(axisparent, ax; position=:bottom, label="MJD") -> Axis

Add a companion x-axis displaying raw MJD values below/above an epoch axis built with MJDConversion. Calendar dates say when; the MJD number is what you type back into a script, so publication figures generally want both. Requires Makie to be loaded.