Introduction

This page walks through the whole path: build a system, solve it at a set of epochs, and query observables.

Bodies

A PlanetOrbits.Body is a point mass. Masses are in solar masses; the constants msun, mjup and mearth are exported so you can write the unit you mean.

using PlanetOrbits
import PlanetOrbits as PO

A = PO.Body(mass=1.1, name=:A)
b = PO.Body(mass=5mjup, name=:b)

The name is a label, and it is what observables use to find a body later. Give every body one.

Orbits

An PlanetOrbits.Orbit is one Keplerian relationship: the orbit of something about something else. The reference is always explicit.

o = PO.Orbit(b, about=A; a=8.0, e=0.1, i=0.5, ω=1.1, Ω=2.2, tp=58849.0)

Angles are radians, a is in AU, and tp — the epoch of periastron passage — is an MJD. There is no M keyword: the gravitating mass of an orbit is the total mass of the bodies it binds, taken from the Body values themselves.

You can give P [days] instead of a, and there are alternative parametrizations for the eccentricity and the orbital phase — see Parametrizations.

Systems

A PlanetOrbits.System is a list of bodies plus the orbits relating them. There must be exactly one orbit fewer than there are bodies: each orbit supplies one relative coordinate, and the remaining degree of freedom is the system barycentre.

sys = PO.System((A, b), (o,); plx=24.5)

System{2 bodies, 1 orbits, Float64} — two-body
  frame: PlanetOrbits.Parallax{Float64}(24.5, 24.5)
  bodies:
    A          mass = 1.1 M⊙
    b          mass = 0.00477297 M⊙
  orbits:
    1: :b about :A
       a = 8.0         AU   P = 7863.16     d (= 21.5281 yr)   e = 0.1
       i = 0.5         rad  ω = 1.1         rad  Ω = 2.2         rad  tp = 58849.0 MJD
       M = 1.104773 M⊙

The keyword arguments define the system's frame:

The absolute frame adds 3D-motion compensation and light-travel-time correction, which matter for nearby, high-proper-motion stars.

Solving

orbitsolve takes the system and a sorted vector of epochs, and returns a Trajectory — per-body barycentric states at every epoch.

epochs = collect(range(58800.0, 59600.0, length=5))
traj = orbitsolve(sys, epochs)
sol = traj[1]

Indexing a trajectory gives a per-epoch solution, which is a cheap view rather than a copy. Solving all epochs at once is the primary entry point: it is what N-body integration requires, and it lets the Kepler solves batch across epochs under SIMD. There is a single-epoch convenience, orbitsolve(sys, t), for interactive use.

Observables

Every observable is a query between two references. A reference can be a body (by value, by Symbol, or by a handle from bodies), a barycentre, or a photocentre.

raoff(sol, b, A)          # RA offset of b relative to A [mas]

-8.382494812239646

radvel(sol, A, barycentre(sys))   # reflex radial velocity of the star [m/s]

-12.46997163641763

The argument order reads as "of, relative to". These are all equivalent:

refs = bodies(sys)
(raoff(sol, b, A), raoff(sol, :b, :A), raoff(sol, refs.b, refs.A))

(-8.382494812239646, -8.382494812239646, -8.382494812239646)

bodies(sys) returns persistent, isbits handles. Named Body values and Symbols resolve to them by name at compile time, so all three spellings cost the same in a hot loop.

The available observables:

Asking for an angular observable on a system with no parallax is an error rather than a silently wrong number.

Over a whole trajectory

Every observable broadcasts over a trajectory, giving one value per epoch. References broadcast as scalars, whichever spelling you use:

radvel.(traj, A, barycentre(sys))   # the star's reflex curve, one value per epoch

5-element Vector{Float64}:
 -12.46997163641763
  -8.350126181078021
  -3.973285596260631
   0.4617992370782518
   4.764746639920423

raoff.(traj, :b, :A)

5-element Vector{Float64}:
   -8.382494812239646
  -41.21360355228618
  -72.61210898583201
 -101.51925680358156
 -127.02008247682802

Broadcasting materializes a vector of per-epoch views, which is convenient but not free; the allocation-free hot loop indexes traj[k] directly instead.

The exception is the observer position taken by the four-argument forms (see Observer-aware observables (opt-in, per read)), which is not a scalar: it broadcasts element-wise, so a vector of positions is read one per epoch. To use a single position at every epoch, wrap it — raoff.(traj, :b, :A, Ref(obs_pos)).

Barycentres and photocentres

barycentre(sys) is the whole-system barycentre; barycentre(sys, members...) is the barycentre of a subsystem. photocentre(sys; band) is the flux-weighted point — what a blended astrometric source actually measures.

Al = PO.Body(mass=1.0, flux=(G=1.0,), name=:A)
bl = PO.Body(mass=0.2, flux=(G=1.0,), name=:b)
lum = PO.System((Al, bl),
    (PO.Orbit(bl, about=Al; a=4.0, e=0.1, i=0.5, tp=58849.0),); plx=30.0)
lsol = orbitsolve(lum, 59000.0)

# equal brightness ⇒ the photocentre sits midway between the two bodies
mid = (raoff(lsol, :A, barycentre(lum)) + raoff(lsol, :b, barycentre(lum))) / 2
(raoff(lsol, photocentre(lum), barycentre(lum)), mid)

(13.4102684647594, 13.4102684647594)

Bodies carry per-band fluxes, so photocentre(sys; band=:G) selects a weight set. Setting the host's flux to 1.0 makes every other body's flux a contrast ratio.

Orbit properties

(period(sys), semimajoraxis(sys), eccentricity(sys), totalmass(sys), distance(sys))

(7863.155915195313, 8.0, 0.1, 1.1047729711698466, 40.816326530612244)

For a system with more than one orbit, pass a row index: period(sys, 2).

Allocation-free use

For hot loops, preallocate the trajectory once and reuse it. orbitsolve! performs no allocation of its own, so with caller-owned storage the whole construct → solve → query path is allocation-free — including under ForwardDiff Duals.

traj_buf = Trajectory(sys, epochs)

function total_sep(θ, buf)
    A = PO.Body(mass=θ[1], name=:A)
    b = PO.Body(mass=θ[2], name=:b)
    s = PO.System((A, b),
        (PO.Orbit(b, about=A; a=θ[3], e=θ[4], i=θ[5], tp=58849.0),); plx=24.5)
    orbitsolve!(buf, s)
    sum(raoff(x, :b, :A) for x in buf)
end

θ = [1.1, 5mjup, 8.0, 0.1, 0.5]
total_sep(θ, traj_buf)          # warm up
@allocated total_sep(θ, traj_buf)

0

Note the system is rebuilt from scratch on every call: System values are cheap, immutable and isbits, and that is the intended usage under MCMC.

Gradients flow through the same path. The only extra step is allocating the buffer with the Dual element type, Trajectory{eltype(θ)}:

import ForwardDiff

function total_sep_grad(θ)
    A = PO.Body(mass=θ[1], name=:A)
    b = PO.Body(mass=θ[2], name=:b)
    s = PO.System((A, b),
        (PO.Orbit(b, about=A; a=θ[3], e=θ[4], i=θ[5], tp=58849.0),); plx=24.5)
    buf = Trajectory{eltype(θ)}(s, epochs)
    orbitsolve!(buf, s)
    sum(raoff(x, :b, :A) for x in buf)
end

ForwardDiff.gradient(total_sep_grad, θ)

5-element Vector{Float64}:
  102.64592836325107
  102.65367122626425
  -10.89449875290758
  214.16336367447332
 -138.24944401846284