Fit Radial Velocity and Astrometry
You can use Octofitter to jointly fit relative astrometry data and radial velocity data. Below is an example. For more information on these functions, see previous guides.
Import required packages
using Octofitter
using OctofitterRadialVelocity
using CairoMakie
using PairPlots
using Distributions
using PlanetOrbitsWe now use PlanetOrbits.jl to create sample data. We start with a template orbit and record it's positon and velocity at a few epochs.
orb_template = orbit(
a = 1.0,
e = 0.7,
i= pi/4,
Ω = 0.1,
ω = 1π/4,
M = 1.0,
plx=100.0,
m =0,
tp =58829-40
)
Makie.lines(orb_template,axis=(;autolimitaspect=1))
Sample position and store as relative astrometry measurements:
epochs = [58849,58852,58858,58890]
astrom = PlanetRelAstromLikelihood(Table(
epoch=epochs,
ra=raoff.(orb_template, epochs),
dec=decoff.(orb_template, epochs),
σ_ra=fill(1.0, size(epochs)),
σ_dec=fill(1.0, size(epochs)),
))PlanetRelAstromLikelihood Table with 5 columns and 4 rows:
epoch ra dec σ_ra σ_dec
┌───────────────────────────────────────
1 │ 58849 -18.3017 -108.907 1.0 1.0
2 │ 58852 -21.1181 -111.432 1.0 1.0
3 │ 58858 -26.6349 -115.889 1.0 1.0
4 │ 58890 -53.1334 -129.043 1.0 1.0And plot our simulated astrometry measurments:
fig = Makie.lines(orb_template,axis=(;autolimitaspect=1))
Makie.scatter!(astrom.table.ra, astrom.table.dec)
fig
Generate a simulated RV curve from the same orbit:
using Random
Random.seed!(1)
epochs = 58849 .+ range(0,step=1.5, length=20)
planet_sim_mass = 0.001 # solar masses here
rvlike = MarginalizedStarAbsoluteRVLikelihood(
Table(
epoch=epochs,
rv=radvel.(orb_template, epochs, planet_sim_mass) .+ 150,
σ_rv=fill(5.0, size(epochs)),
),
jitter=:jitter1,
instrument_name="inst1"
)
epochs = 58949 .+ range(0,step=1.5, length=20)
rvlike2 = MarginalizedStarAbsoluteRVLikelihood(
Table(
epoch=epochs,
rv=radvel.(orb_template, epochs, planet_sim_mass) .- 150,
σ_rv=fill(5.0, size(epochs)),
),
jitter=:jitter1,
instrument_name="inst2"
)
fap = Makie.scatter(rvlike.table.epoch[:], rvlike.table.rv[:])
Makie.scatter!(rvlike2.table.epoch[:], rvlike2.table.rv[:])
fap
Now specify model and fit it:
@planet b Visual{KepOrbit} begin
e ~ Uniform(0,0.999999)
a ~ truncated(Normal(1, 1),lower=0.1)
mass ~ truncated(Normal(1, 1), lower=0.)
i ~ Sine()
Ω ~ UniformCircular()
ω ~ UniformCircular()
θ ~ UniformCircular()
tp = θ_at_epoch_to_tperi(system,b,58849.0) # reference epoch for θ. Choose an MJD date near your data.
end astrom
@system test begin
M ~ truncated(Normal(1, 0.04),lower=0.1) # (Baines & Armstrong 2011).
plx = 100.0
jitter1 ~ truncated(Normal(0,10),lower=0)
end rvlike rvlike2 b
model = Octofitter.LogDensityModel(test)
using Random
rng = Xoshiro(0) # seed the random number generator for reproducible results
results = octofit(rng, model, max_depth=9, adaptation=300, iterations=400)Chains MCMC chain (400×30×1 Array{Float64, 3}):
Iterations = 1:1:400
Number of chains = 1
Samples per chain = 400
Wall duration = 5.06 seconds
Compute duration = 5.06 seconds
parameters = M, jitter1, plx, b_e, b_a, b_mass, b_i, b_Ωy, b_Ωx, b_ωy, b_ωx, b_θy, b_θx, b_Ω, b_ω, b_θ, b_tp
internals = n_steps, is_accept, acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size, is_adapt, loglike, logpost, tree_depth, numerical_error
Summary Statistics
parameters mean std mcse ess_bulk ess_tail rha ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float6 ⋯
M 0.9980 0.0344 0.0029 140.8501 141.8130 1.060 ⋯
jitter1 0.3978 0.3506 0.0372 122.2376 150.0601 1.129 ⋯
plx 100.0000 0.0000 NaN NaN NaN Na ⋯
b_e 0.2124 0.1683 0.1387 1.7236 21.3676 1.648 ⋯
b_a 1.6877 0.4728 0.2844 2.5675 34.2740 1.346 ⋯
b_mass 1.0804 0.6152 0.1825 8.3551 28.3530 1.101 ⋯
b_i 1.1418 0.0841 0.0617 2.3690 17.9170 1.395 ⋯
b_Ωy -0.7224 0.1113 0.0184 40.0639 61.2438 1.113 ⋯
b_Ωx -0.7043 0.0988 0.0339 8.3159 73.3527 1.152 ⋯
b_ωy -0.2579 0.7626 0.5505 2.4402 46.4968 1.403 ⋯
b_ωx 0.3651 0.5004 0.0763 52.7406 57.7689 1.381 ⋯
b_θy -0.9621 0.0845 0.0105 69.8198 154.9883 1.016 ⋯
b_θx -0.1624 0.0145 0.0013 112.9723 181.9498 1.132 ⋯
b_Ω -2.3678 0.1174 0.0323 13.3347 64.4314 1.121 ⋯
b_ω 1.5183 1.4815 1.0531 2.0058 130.5647 1.512 ⋯
b_θ -2.9742 0.0056 0.0009 38.9323 84.7630 1.023 ⋯
b_tp 58414.5495 344.4624 206.8557 4.4497 28.3333 1.200 ⋯
2 columns omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
M 0.9309 0.9822 0.9950 1.0173 1.0700
jitter1 0.0163 0.1908 0.2674 0.5425 1.3958
plx 100.0000 100.0000 100.0000 100.0000 100.0000
b_e 0.0026 0.0656 0.1509 0.3922 0.5052
b_a 1.1107 1.2274 1.6192 1.9671 2.6967
b_mass 0.0662 0.6536 1.0328 1.4675 2.1889
b_i 1.0157 1.0374 1.1649 1.2068 1.2755
b_Ωy -0.9750 -0.7839 -0.7084 -0.6618 -0.5248
b_Ωx -0.9044 -0.7477 -0.7182 -0.6447 -0.4934
b_ωy -1.0753 -0.8903 -0.7264 0.4727 1.0898
b_ωx -0.9741 0.2861 0.3939 0.6595 1.0637
b_θy -1.1618 -1.0114 -0.9594 -0.9132 -0.8109
b_θx -0.1965 -0.1704 -0.1614 -0.1541 -0.1365
b_Ω -2.6156 -2.4350 -2.3562 -2.2966 -2.1440
b_ω -2.1649 0.4877 2.2542 2.7343 2.9034
b_θ -2.9864 -2.9778 -2.9738 -2.9696 -2.9647
b_tp 57559.8493 58222.8478 58536.6641 58697.0249 58791.5400
Display results as a corner plot:
octocorner(model,results, small=true)
Plot RV curve, phase folded curve, and binned residuals:
Octofitter.rvpostplot(model, results)
Display RV, PMA, astrometry, relative separation, position angle, and 3D projected views:
octoplot(model, results)