Hierarchical systems

A system is a flat list of bodies plus a flat list of orbits. Each orbit says what orbits what, and nothing is inferred. That one decision is what lets the same API describe a lone planet, a moon, a circumbinary planet, or a 2+2 quadruple.

The reference grammar

Both endpoints of an orbit — the thing orbiting and the about= it orbits — take the same grammar:

That is all of it. A body, or a set of bodies meaning their barycentre.

using PlanetOrbits
import PlanetOrbits as PO

A = PO.Body(mass=1.1, name=:A)
b = PO.Body(mass=8mjup, name=:b)
c = PO.Body(mass=2mjup, name=:c)

Jacobi vs. astrocentric

For a star with two planets there are two natural conventions, and they are different models, not different spellings:

els1 = (; a=2.5, e=0.1, i=0.5, ω=1.1, Ω=2.2, tp=58849.0)
els2 = (; a=8.0, e=0.3, i=0.6, ω=0.4, Ω=2.0, tp=57000.0)

jacobi = PO.System((A, b, c), (
    PO.Orbit(b, about=A;      els1...),
    PO.Orbit(c, about=(A, b); els2...),   # c orbits the A+b barycentre
); plx=50.0)

astro = PO.System((A, b, c), (
    PO.Orbit(b, about=A; els1...),
    PO.Orbit(c, about=A; els2...),        # c orbits the star itself
); plx=50.0)

Under the default KeplerianApprox propagator, each orbit is solved independently — so the orbits are the model. The two systems above give observably different predictions from identical element values:

ts = [58000.0, 59000.0, 60000.0]
tj = orbitsolve(jacobi, ts); ta = orbitsolve(astro, ts)
maximum(abs(raoff(tj[k], :c, :A) - raoff(ta[k], :c, :A)) for k in 1:3)  # mas

1.644869778016215

They also imply different gravitating masses per orbit — Jacobi's outer orbit is bound by all three bodies, astrocentric's by only two:

(jacobi.rows[2].M, astro.rows[2].M)

(1.1095459423396934, 1.1019091884679388)

There is deliberately no default and no inference from semi-major axis. Pick one and say so.

Convention constructors

Jacobi and Astrocentric expand to the explicit spellings above, which is easier to read for long chains:

jacobi2 = PO.System((A, b, c), Jacobi(A, b => els1, c => els2); plx=50.0)
jacobi2.Ainv ≈ jacobi.Ainv

true

show reports the convention it resolved:

astro

System{3 bodies, 2 orbits, Float64} — astrocentric
  frame: PlanetOrbits.Parallax{Float64}(50.0, 50.0)
  bodies:
    A          mass = 1.1 M⊙
    b          mass = 0.00763675 M⊙
    c          mass = 0.00190919 M⊙
  orbits:
    1: :b about :A
       a = 2.5         AU   P = 1371.86     d (= 3.75595 yr)   e = 0.1
       i = 0.5         rad  ω = 1.1         rad  Ω = 2.2         rad  tp = 58849.0 MJD
       M = 1.1076368 M⊙
    2: :c about :A
       a = 8.0         AU   P = 7873.37     d (= 21.5561 yr)   e = 0.3
       i = 0.6         rad  ω = 0.4         rad  Ω = 2.0         rad  tp = 57000.0 MJD
       M = 1.1019092 M⊙

Moons

A moon orbits its host, not the system barycentre. This is the case that tree-shaped representations cannot reach, because the host would have to appear in two places at once — as the thing the moon orbits, and as the thing that orbits the star. In a flat list it is just another about=:

star = PO.Body(mass=1.0, name=:A)
planet = PO.Body(mass=10mjup, name=:b)
moon = PO.Body(mass=1mearth, name=:m)

sysm = PO.System((star, planet, moon), (
    PO.Orbit(planet, about=star;   a=5.2,   e=0.05, i=0.5, ω=1.1, Ω=2.2, tp=58849.0),
    PO.Orbit(moon,   about=planet; a=0.007, e=0.0,  i=0.4, ω=0.0, Ω=1.0, tp=58849.0),
); plx=50.0)

solm = orbitsolve(sysm, 58900.0)
(projectedseparation(solm, :m, :b), projectedseparation(solm, :b, :A))

(0.32427401233059794, 221.39240575114542)

Circumbinary planets

A planet orbiting a tight binary is the barycentre spelling:

Aa = PO.Body(mass=0.6, name=:Aa)
Ab = PO.Body(mass=0.4, name=:Ab)
p  = PO.Body(mass=1mjup, name=:b)

cb = PO.System((Aa, Ab, p), (
    PO.Orbit(Ab, about=Aa;        a=0.2, e=0.05, i=0.3, ω=0.2, Ω=0.1, tp=58849.0),
    PO.Orbit(p,  about=(Aa, Ab);  a=3.0, e=0.1,  i=0.4, ω=1.0, Ω=2.0, tp=58000.0),
); plx=80.0)

Set exteriors: 2+2 quadruples

Both endpoints take the grammar, so the orbiting side can be a barycentre too. A 2+2 quadruple is two tight pairs orbiting each other:

Ba = PO.Body(mass=0.8, name=:Ba)
Bb = PO.Body(mass=0.7, name=:Bb)

quad = PO.System((Aa, Ab, Ba, Bb), (
    PO.Orbit(Ab, about=Aa; a=0.5, e=0.1, i=0.3, ω=0.2, Ω=0.1, tp=58849.0),
    PO.Orbit(Bb, about=Ba; a=0.6, e=0.2, i=0.4, ω=0.3, Ω=0.2, tp=58849.0),
    PO.Orbit((Ba, Bb), about=(Aa, Ab); a=50.0, e=0.3, i=0.5, ω=0.4, Ω=0.3, tp=58000.0),
); plx=20.0)

System{4 bodies, 3 orbits, Float64} — Jacobi (hierarchical)
  frame: PlanetOrbits.Parallax{Float64}(20.0, 20.0)
  bodies:
    Aa         mass = 0.6 M⊙
    Ab         mass = 0.4 M⊙
    Ba         mass = 0.8 M⊙
    Bb         mass = 0.7 M⊙
  orbits:
    1: :Ab about :Aa
       a = 0.5         AU   P = 129.138     d (= 0.35356 yr)   e = 0.1
       i = 0.3         rad  ω = 0.2         rad  Ω = 0.1         rad  tp = 58849.0 MJD
       M = 1.0 M⊙
    2: :Bb about :Ba
       a = 0.6         AU   P = 138.605     d (= 0.37948 yr)   e = 0.2
       i = 0.4         rad  ω = 0.3         rad  Ω = 0.2         rad  tp = 58849.0 MJD
       M = 1.5 M⊙
    3: barycentre(:Ba, :Bb) about barycentre(:Aa, :Ab)
       a = 50.0        AU   P = 81673.9     d (= 223.611 yr)   e = 0.3
       i = 0.5         rad  ω = 0.4         rad  Ω = 0.3         rad  tp = 58000.0 MJD
       M = 2.5 M⊙

Blended sources & photocentres

A catalog astrometric source is not a body. It is whatever flux fell into the instrument's aperture, reduced to a point — so the thing to model is a flux-weighted point over the bodies that were blended into it. That is a WeightedPoint, exactly as a barycentre is, and observables take it anywhere a body goes.

PlanetOrbits owns the two generic linear reductions — the mass-weighted barycentre and the flux-weighted photocentre. Anything instrument-specific — a grating response, a per-epoch resolution taper, a scan-angle-dependent window — belongs to the observation layer, which builds on these primitives.

Bodies carry per-band fluxes; the band selects the weight set. Units are arbitrary but must be consistent within a band, so setting the host to 1.0 makes every other body's number a contrast ratio.

The worked example: two sources in a 2+2 quadruple

The quadruple above, with fluxes, observed as two catalog sources 50 AU apart — each of which blends only its own pair:

Aaf = PO.Body(mass=0.6, flux=(G=1.0,),  name=:Aa)
Abf = PO.Body(mass=0.4, flux=(G=0.25,), name=:Ab)
Baf = PO.Body(mass=0.8, flux=(G=0.8,),  name=:Ba)
Bbf = PO.Body(mass=0.7, flux=(G=0.5,),  name=:Bb)

quadf = PO.System((Aaf, Abf, Baf, Bbf), (
    PO.Orbit(Abf, about=Aaf; a=0.5, e=0.1, i=0.3, ω=0.2, Ω=0.1, tp=58849.0),
    PO.Orbit(Bbf, about=Baf; a=0.6, e=0.2, i=0.4, ω=0.3, Ω=0.2, tp=58849.0),
    PO.Orbit((Baf, Bbf), about=(Aaf, Abf); a=50.0, e=0.3, i=0.5, ω=0.4, Ω=0.3, tp=58000.0),
); plx=20.0)

srcA = photocentre(quadf, :Aa, :Ab; band=:G)
srcB = photocentre(quadf, :Ba, :Bb; band=:G)
sol = orbitsolve(quadf, 59000.0)
raoff(sol, srcB, srcA)   # source-to-source separation [mas]

485.2467237637502

Each source's motion is one dot product over absolute body states, which is why it needs no per-level bookkeeping: srcA carries both the Aa–Ab photocentric wobble and the A-pair's motion on the wide orbit, under either propagator. Compare the source against the pair's own barycentre to see the wobble alone:

raoff(sol, srcA, barycentre(quadf, :Aa, :Ab))

-1.8344509613526725

A single-member subset degrades to that body, so a dark companion needs no special case:

raoff(sol, photocentre(quadf, :Aa; band=:G), :Aa)

0.0

Membership that changes per draw or per epoch

Structural membership — "these two can never be resolved apart" — is the photocentre(sys, members...) form above. When membership is itself part of the model (a sampled resolved-flag, or a taper in separation that only the instrument knows about), read the fluxes and build the point yourself:

f = fluxes(quadf, :G)                  # SVector, in `bodies(quadf)` order
member = PlanetOrbits.SVector(1.0, 1.0, 0.0, 0.0)   # e.g. a per-epoch gate
wp = photocentre(f .* member)
raoff(sol, wp, barycentre(quadf)) ≈ raoff(sol, srcA, barycentre(quadf))

true

photocentre(w) just normalizes; WeightedPoint is isbits, so building one per epoch inside a scan loop allocates nothing. Use the structural form when membership is fixed — it is validated at model-build time — and the per-draw form when it is not.

Mixed conventions

Mixing is legal, and show says so rather than pretending the system has one convention. A moon spelled about its host, inside a chain that is otherwise Jacobi, is a mixture — and that is a real modelling statement, not a formatting detail.

Validation

The orbit set has to determine every body's position. When it does not, the error names the rows involved:

try
    PO.System((A, b, c), (PO.Orbit(b, about=A; a=1.0),))
catch err
    println(sprint(showerror, err))
end

a system of 3 bodies needs exactly 2 orbits to determine every body's motion; got 1. (Each orbit supplies one relative coordinate; the remaining degree of freedom is the system barycentre.)

Zero-mass bodies

Test particles are supported throughout. A massless set has no mass-weighted barycentre, so its limit — the members' geometric centre, which for a single particle is just its own position — is used instead. The star sits exactly at the barycentre of a system whose companions are massless:

z = PO.System(
    (PO.Body(mass=1.0, name=:A), PO.Body(mass=0.0, name=:b)),
    (PO.Orbit(PO.Body(mass=0.0, name=:b), about=PO.Body(mass=1.0, name=:A);
              a=5.0, e=0.2, i=0.4, tp=58849.0),); plx=40.0)
abs(raoff(orbitsolve(z, 58900.0), :A, barycentre(z)))

0.0

The same limit applies to the system barycentre, so a system in which every body is a test particle is built rather than rejected. Masses may be any finite value and any scale — the hierarchy matrix depends only on their ratios, and is formed so that neither their sum overflowing nor an extreme ratio can turn a well-posed hierarchy into a singular one. A mass that is not finite is a different matter, and raises a PlanetOrbits.OrbitDomainError naming the body.

How it works

Internally each orbit contributes one row to a hierarchy matrix H, which maps absolute barycentric positions to the relative coordinates the orbits describe. Inverting H and dropping the barycentre column gives sys.Ainv, which turns per-orbit relative states back into per-body absolute states.

This is one code path for every convention. It is a little slower than the closed form a pure Jacobi chain admits — around 20 ns against a 25 µs likelihood evaluation — and in exchange there is no second path in which an astrocentric system can be silently evaluated with the Jacobi formula.