Parametrizations

An orbit needs six numbers. Which six is up to you: Orbit accepts several standard parametrizations, organised into groups. Supply exactly one alternative from each group.

Or replace all of them at once with Cartesian initial conditions, or replace size and orientation together with Thiele-Innes constants. Supplying two alternatives from the same group is an error at construction.

using PlanetOrbits
import PlanetOrbits as PO
A = PO.Body(mass=1.1, name=:A); b = PO.Body(mass=5mjup, name=:b)

Size: a or P

sysa = PO.System((A, b), (PO.Orbit(b, about=A; a=8.0, e=0.1),))
sysp = PO.System((A, b), (PO.Orbit(b, about=A; P=period(sysa), e=0.1),))
(semimajoraxis(sysp), period(sysp))

(7.999999999999998, 7863.155915195313)

P is converted with the orbit's own gravitating mass, so under different conventions the same P gives a different a — as it should.

For a bound orbit both must be positive and finite. A non-positive or non-finite value is rejected at construction with a PlanetOrbits.OrbitDomainError, which a likelihood catches as a quiet -Inf — worth knowing because an ordinary prior like a ~ Uniform(0, 100) can propose exactly 0.0 during sampling. (Hyperbolic orbits are the exception: there a < 0 is the convention — see below.)

Shape

secosω = √e·cosω and ecosω = e·cosω sample the eccentricity disc rather than the half-plane, which removes the ω degeneracy as e → 0. That is usually what you want when sampling.

e, ω = 0.3, 1.1
s1 = PO.System((A, b), (PO.Orbit(b, about=A; a=5.0, e=e, ω=ω),))
s2 = PO.System((A, b), (PO.Orbit(b, about=A; a=5.0, secosω=√e*cos(ω), sesinω=√e*sin(ω)),))
(eccentricity(s1), eccentricity(s2))

(0.3, 0.3)

Phase

tp is the canonical element — the epoch of periastron passage, in MJD. The alternatives are measured at an epoch you supply:

• M0 — mean anomaly at epoch [rad].
• θ — sky-plane position angle at epoch [rad].

base = PO.System((A, b), (PO.Orbit(b, about=A; a=5.0, e=0.4, ω=1.1, i=0.5, Ω=2.2, tp=59000.0),); plx=25.0)
ep = 59211.0
θ_obs = posangle(orbitsolve(base, ep), :b, :A)

viaθ = PO.System((A, b),
    (PO.Orbit(b, about=A; a=5.0, e=0.4, ω=1.1, i=0.5, Ω=2.2, θ=θ_obs, epoch=ep),); plx=25.0)
(periastron(base), periastron(viaθ))

(59000.0, 59000.00851692485)

θ is often the best-constrained quantity for a directly imaged planet with a short observational arc, which is why it is worth parametrizing on directly.

Cartesian initial conditions

Give the full relative state instead of elements. This determines everything, so it replaces every group above:

sol = orbitsolve(base, 59123.0)
cart = PO.System((A, b), (PO.Orbit(b, about=A;
    x  = posx(sol, :b, :A), y  = posy(sol, :b, :A), z  = posz(sol, :b, :A),
    vx = velx(sol, :b, :A), vy = vely(sol, :b, :A), vz = velz(sol, :b, :A),
    epoch = 59123.0),); plx=25.0)
(semimajoraxis(cart), eccentricity(cart), periastron(cart))

(5.000000051529844, 0.40000001184633494, 59000.008510690204)

Positions are AU and velocities AU / julian year — the same frame and units as posx and velx — so a state read out of a solution reconstructs the same orbit exactly, as above.

Because a comes from vis-viva and e from the eccentricity vector, unbound states need no special handling: they simply come out with a < 0 and e > 1.

Hyperbolic orbits

e > 1 is supported, with velocities. The semi-major axis is negative by convention; a positive value is taken as |a|.

hyp = PO.System((A, PO.Body(mass=0.0, name=:b)),
    (PO.Orbit(PO.Body(mass=0.0, name=:b), about=A;
              a=-5.0, e=1.4, i=0.5, ω=1.1, Ω=2.2, tp=59000.0),); plx=25.0)
(period(hyp), semimajoraxis(hyp))

(Inf, -5.0)

The period is Inf, and P= is rejected for unbound orbits. e == 1 exactly is also rejected — the elements are degenerate for parabolae, so use Cartesian initial conditions there. (For why e = 1 is a genuinely singular case rather than an implementation gap, see Parabolic Trajectories; the bound and unbound cases are Elliptical and Hyperbolic.)

Thiele-Innes constants

ThieleInnes returns the size and orientation elements to splat in, so it composes with the shape and phase groups:

ti = thieleinnes(base)

(A = -4.496366802697549, B = -0.46770310673372756, F = 1.0131997159495527, G = -4.774007178304462)

viaTI = PO.System((A, b),
    (PO.Orbit(b, about=A; ThieleInnes(; ti...)..., e=0.4, tp=59000.0),); plx=25.0)
(semimajoraxis(viaTI), inclination(viaTI))

(5.0, 0.4999999999999999)

Pass plx if your constants are in mas rather than AU. This parametrization has no coordinate singularity at e → 0 or i → 0, which is why Gaia-only astrometric fits use it.

Radial-velocity-only fits

Radial velocities constrain m·sin i, never m and i separately, and give no handle on Ω or the parallax. The convention is i = π/2, Ω = 0, and no parallax — so angular observables are unavailable rather than silently wrong:

using PlanetOrbits: rvorbit
rv = rvorbit(M=1.1, msini=2mjup, P=400.0, e=0.2, ω=1.0, tp=59000.0)
radvel(orbitsolve(rv, 59100.0), :A, barycentre(rv))

46.2186019978907

Under this convention the secondary's mass is its m·sin i; reading it as a true mass gives a lower bound.

The M= compatibility override

An orbit's gravitating mass normally comes from the bodies it binds. M= overrides it:

ov = PO.System((A, b), (PO.Orbit(b, about=A; a=8.0, M=2.5),))
(ov.rows[1].M, period(ov))

(2.5, 5227.131225839685)

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. show labels it as a compatibility override. It is not a modelling choice.