Hierarchical Co-Planar, Resonant Model
This example shows how you can fit a two planet model to relative astrometry data. This functionality would work equally well with RV, images, etc.
We will restrict the two planets to being exactly co-planar, and have a near 2:1 period resonance, and zero-eccentricity.
For this example, we will use astrometry from the HR8799 system collated by Jason Wang and retrieved from the website Whereistheplanet.com.
using Octofitter
using CairoMakie
using PairPlots
using Distributions
using PlanetOrbits┌ Warning: Module Octofitter with build ID ffffffff-ffff-ffff-0000-004287710a13 is missing from the cache.
│ This may mean Octofitter [daf3887e-d01a-44a1-9d7e-98f15c5d69c9] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:1948astrom_b = PlanetRelAstromLikelihood(
(epoch=53200.0, object=1, ra=1471.0, σ_ra=6.0, dec=887.0, σ_dec=6.0, cor=0, ),
(epoch=54314.0, object=1, ra=1504.0, σ_ra=3.0, dec=837.0, σ_dec=3.0, cor=0, ),
(epoch=54398.0, object=1, ra=1500.0, σ_ra=7.0, dec=836.0, σ_dec=7.0, cor=0, ),
(epoch=54727.0, object=1, ra=1516.0, σ_ra=4.0, dec=818.0, σ_dec=4.0, cor=0, ),
(epoch=55042.0, object=1, ra=1526.0, σ_ra=4.0, dec=797.0, σ_dec=4.0, cor=0, ),
(epoch=55044.0, object=1, ra=1531.0, σ_ra=7.0, dec=794.0, σ_dec=7.0, cor=0, ),
(epoch=55136.0, object=1, ra=1524.0, σ_ra=10.0, dec=795.0, σ_dec=10.0, cor=0, ),
(epoch=55390.0, object=1, ra=1532.0, σ_ra=5.0, dec=783.0, σ_dec=5.0, cor=0, ),
(epoch=55499.0, object=1, ra=1535.0, σ_ra=15.0, dec=766.0, σ_dec=15.0, cor=0, ),
(epoch=55763.0, object=1, ra=1541.0, σ_ra=5.0, dec=762.0, σ_dec=5.0, cor=0, ),
(epoch=56130.0, object=1, ra=1545.0, σ_ra=5.0, dec=747.0, σ_dec=5.0, cor=0, ),
(epoch=56226.0, object=1, ra=1549.0, σ_ra=4.0, dec=743.0, σ_dec=4.0, cor=0, ),
(epoch=56581.0, object=1, ra=1545.0, σ_ra=22.0, dec=724.0, σ_dec=22.0, cor=0, ),
(epoch=56855.0, object=1, ra=1560.0, σ_ra=13.0, dec=725.0, σ_dec=13.0, cor=0, ),
(epoch=58798.03906, object=1, ra=1611.002, σ_ra=0.133, dec=604.893, σ_dec=0.199, cor=-0.406, ),
(epoch=59453.245, object=1, ra=1622.924, σ_ra=0.32, dec=570.534, σ_dec=0.296, cor=-0.905, ),
(epoch=59454.231, object=1, ra=1622.872, σ_ra=0.204, dec=571.296, σ_dec=0.446, cor=-0.79, ),
)
astrom_c = PlanetRelAstromLikelihood(
(epoch=53200.0, object=2, ra=-739.0, σ_ra=6.0, dec=612.0, σ_dec=6.0, cor=0, ),
(epoch=54314.0, object=2, ra=-683.0, σ_ra=4.0, dec=671.0, σ_dec=4.0, cor=0, ),
(epoch=54398.0, object=2, ra=-678.0, σ_ra=7.0, dec=678.0, σ_dec=7.0, cor=0, ),
(epoch=54727.0, object=2, ra=-663.0, σ_ra=3.0, dec=693.0, σ_dec=3.0, cor=0, ),
(epoch=55042.0, object=2, ra=-639.0, σ_ra=4.0, dec=712.0, σ_dec=4.0, cor=0, ),
(epoch=55136.0, object=2, ra=-636.0, σ_ra=9.0, dec=720.0, σ_dec=9.0, cor=0, ),
(epoch=55390.0, object=2, ra=-619.0, σ_ra=4.0, dec=728.0, σ_dec=4.0, cor=0, ),
(epoch=55499.0, object=2, ra=-607.0, σ_ra=12.0, dec=744.0, σ_dec=12.0, cor=0, ),
(epoch=55763.0, object=2, ra=-595.0, σ_ra=4.0, dec=747.0, σ_dec=4.0, cor=0, ),
(epoch=56130.0, object=2, ra=-578.0, σ_ra=5.0, dec=761.0, σ_dec=5.0, cor=0, ),
(epoch=56226.0, object=2, ra=-572.0, σ_ra=3.0, dec=768.0, σ_dec=3.0, cor=0, ),
(epoch=56581.0, object=2, ra=-542.0, σ_ra=22.0, dec=784.0, σ_dec=22.0, cor=0, ),
(epoch=56855.0, object=2, ra=-540.0, σ_ra=12.0, dec=799.0, σ_dec=12.0, cor=0, ),
)
fig = Makie.scatter(astrom_b.table.ra, astrom_b.table.dec, axis=(;autolimitaspect=1))
Makie.scatter!(astrom_c.table.ra, astrom_c.table.dec)
Makie.scatter!([0], [0], marker='⋆', markersize=50, color=:black)
fig
We now specify our two planet model for planets b & c.
@planet b Visual{KepOrbit} begin
e = 0.0
ω = 0.0
# Use the system inclination and longitude of ascending node
# variables
i = system.i
Ω = system.Ω
# Specify the period as ~ 1% around 2X the P_nominal variable
P_mul ~ Normal(1, 0.1)
P = 2*system.P_nominal * b.P_mul
a = cbrt(system.M * b.P^2)
θ ~ UniformCircular()
tp = θ_at_epoch_to_tperi(system,b,59454.231) # reference epoch for θ. Choose an MJD date near your data.
end astrom_b
@planet c Visual{KepOrbit} begin
e = 0.0
ω = 0.0
# Use the system inclination and longitude of ascending node
# variables
i = system.i
Ω = system.Ω
# Specify the period as ~ 1% the P_nominal variable
P_mul ~ truncated(Normal(1, 0.1), lower=0)
P = system.P_nominal * c.P_mul
a = cbrt(system.M * c.P^2)
θ ~ UniformCircular()
tp = θ_at_epoch_to_tperi(system,c,59454.231) # reference epoch for θ. Choose an MJD date near your data.
end astrom_c
@system HR8799_res_co begin
plx ~ gaia_plx(;gaia_id=2832463659640297472)
M ~ truncated(Normal(1.5, 0.02), lower=0)
# We create inclination and longitude of ascending node variables at the
# system level.
i ~ Sine()
Ω ~ UniformCircular()
# We create a nominal period of planet c variable.
P_nominal ~ Uniform(50, 300) # years
end b c
model = Octofitter.LogDensityModel(HR8799_res_co)LogDensityModel for System HR8799_res_co of dimension 12 with fields .ℓπcallback and .∇ℓπcallback
Let's now sample from the model:
using Pigeons
results,pt = octofit_pigeons(model, n_rounds=10);┌ Warning: Module Octofitter with build ID ffffffff-ffff-ffff-0000-004287710a13 is missing from the cache.
│ This may mean Octofitter [daf3887e-d01a-44a1-9d7e-98f15c5d69c9] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:1948
[ Info: Determining initial positions and metric using pathfinder
┌ Info: Found a sample of initial positions
└ initial_logpost_range = (-181.12810535746516, -170.0900294661454)
┌ Warning: Invalid log likelihood encountered. (maxlog=1)
│ θ = (plx = Inf, M = Inf, i = 0.0005726099324290644, Ωy = 4.8690096939477945, Ωx = -3.2264241007218293, P_nominal = 299.89735119519906, Ω = -0.5852130842848838, planets = (b = (P_mul = 1795.599574523399, θy = -973.792566996004, θx = 341.02361331659904, e = 0.0, ω = 0.0, i = 0.0005726099324290644, Ω = -0.5852130842848838, P = 1.0769911124135877e6, a = Inf, θ = 2.804738339425028, tp = NaN), c = (P_mul = 2.8475596141009845e35, θy = 44.311589394218956, θx = 112.03157268236912, e = 0.0, ω = 0.0, i = 0.0005726099324290644, Ω = -0.5852130842848838, P = 8.539755856393085e37, a = Inf, θ = 1.1941513232117695, tp = NaN)))
│ llike = NaN
│ θ_transformed =
│ 12-element Vector{Float64}:
│ 2583.6006322056305
│ ⋮
└ @ Octofitter ~/work/Octofitter.jl/Octofitter.jl/src/logdensitymodel.jl:105
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 3.66 3.17 7.53 3.96e+08 -1.56e+07 0 0.78 0.939 0.966
4 0 4.08 5.2 1.21 9.23e+06 -2.29e+06 0 0.701 0.905 0.939
8 0 4.89 5.02 1.81 1.02e+06 -3.73e+05 0 0.68 0.901 0.927
16 0 7.61 6.96 3.63 1.91e+06 -1.02e+05 0 0.53 0.913 0.931
32 0 9.21 7.91 7.18 3.62e+06 -3.61e+03 0 0.447 0.909 0.926
64 0 9.86 5.76 15.9 2.82e+08 -210 1.67e-38 0.496 0.907 0.924
128 1 10 6.75 30.9 4.85e+08 -212 8.13e-14 0.46 0.908 0.926
256 4 10.2 6.94 62 9.67e+08 -211 2.13e-27 0.448 0.912 0.927
512 4 10.7 8.28 121 1.9e+09 -209 0.00093 0.388 0.913 0.927
1.02e+03 14 11 8.16 253 4.01e+09 -205 0.00108 0.381 0.916 0.928
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────Plots the orbits:
octoplot(model, results)
Corner plot:
octocorner(model, results, small=true)
Now examine the period ratio:
hist(
results[:b_P][:] ./ results[:c_P][:],
axis=(;
xlabel="period ratio",
ylabel="counts",
)
)