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 =58849
)
Makie.lines(orb_template)
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 17.0428 19.6097 1.0 1.0
2 │ 58852 21.3842 9.55631 1.0 1.0
3 │ 58858 24.6966 -12.2283 1.0 1.0
4 │ 58890 0.213769 -87.2649 1.0 1.0And plot our simulated astrometry measurments:
fig = Makie.lines(orb_template,)
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 .+ (20:20:660)
planet_sim_mass = 0.001 # solar masses here
rvlike = StarAbsoluteRVLikelihood(Table(
epoch=epochs,
rv=radvel.(orb_template, epochs, planet_sim_mass) .+ randn.() .+ 100randn(),
σ_rv=fill(5.0, size(epochs)),
inst_idx=ones(Int64, size(epochs)),
))
Makie.scatter(rvlike.table.epoch[:], rvlike.table.rv[:])
Now specify model and fit it:
@planet b Visual{KepOrbit} begin
e ~ Uniform(0,0.999999)
a ~ truncated(Normal(1, 1),lower=0)
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) # (Baines & Armstrong 2011).
plx = 100.0
jitter_1 ~ truncated(Normal(0,10),lower=0)
rv0_1 ~ Normal(0,10)
end rvlike b
model = Octofitter.LogDensityModel(test; autodiff=:ForwardDiff,verbosity=4)
using Random
rng = Xoshiro(0) # seed the random number generator for reproducible results
results = octofit(rng, model)Chains MCMC chain (1000×31×1 Array{Float64, 3}):
Iterations = 1:1:1000
Number of chains = 1
Samples per chain = 1000
Wall duration = 6.32 seconds
Compute duration = 6.32 seconds
parameters = M, jitter_1, rv0_1, 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.9995 0.0279 0.0010 844.7235 701.1969 0.999 ⋯
jitter_1 0.8072 0.6215 0.0178 903.1270 465.6722 0.999 ⋯
rv0_1 -6.7536 0.8934 0.0274 1058.6477 710.9066 1.002 ⋯
plx 100.0000 0.0000 NaN NaN NaN Na ⋯
b_e 0.7019 0.0127 0.0004 918.2507 682.9616 1.000 ⋯
b_a 0.9996 0.0109 0.0004 912.8668 695.7784 0.999 ⋯
b_mass 1.0722 0.0972 0.0031 964.2479 617.7269 1.003 ⋯
b_i 0.7764 0.0599 0.0021 848.6405 720.8847 0.999 ⋯
b_Ωy 0.9986 0.0954 0.0034 799.7514 754.3124 1.001 ⋯
b_Ωx 0.1048 0.0540 0.0022 617.8936 558.0809 1.000 ⋯
b_ωy 0.7178 0.0987 0.0038 685.7096 450.5228 1.003 ⋯
b_ωx 0.7014 0.0911 0.0035 688.3612 539.2990 1.002 ⋯
b_θy 0.7637 0.0814 0.0032 691.6338 600.6544 1.000 ⋯
b_θx 0.6613 0.0705 0.0028 690.0779 402.7450 1.000 ⋯
b_Ω 0.1046 0.0528 0.0022 590.9141 558.0809 1.001 ⋯
b_ω 0.7745 0.0911 0.0037 631.8757 523.9169 1.002 ⋯
b_θ 0.7138 0.0297 0.0009 1054.9739 675.8711 1.012 ⋯
b_tp 58660.1879 182.6232 6.5246 850.6161 828.0038 1.000 ⋯
2 columns omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
M 0.9468 0.9805 0.9989 1.0184 1.0522
jitter_1 0.0479 0.3223 0.6548 1.1492 2.2803
rv0_1 -8.4373 -7.3484 -6.7778 -6.1574 -4.9858
plx 100.0000 100.0000 100.0000 100.0000 100.0000
b_e 0.6760 0.6935 0.7021 0.7109 0.7256
b_a 0.9796 0.9920 0.9995 1.0067 1.0211
b_mass 0.9033 1.0052 1.0664 1.1320 1.2773
b_i 0.6539 0.7377 0.7791 0.8177 0.8858
b_Ωy 0.8229 0.9298 0.9985 1.0633 1.1912
b_Ωx -0.0011 0.0671 0.1055 0.1410 0.2115
b_ωy 0.5362 0.6482 0.7088 0.7834 0.9188
b_ωx 0.5439 0.6401 0.6980 0.7611 0.9050
b_θy 0.6246 0.7088 0.7531 0.8184 0.9407
b_θx 0.5314 0.6145 0.6558 0.7050 0.8059
b_Ω -0.0011 0.0680 0.1052 0.1424 0.2045
b_ω 0.5994 0.7095 0.7763 0.8355 0.9593
b_θ 0.6535 0.6952 0.7146 0.7341 0.7698
b_tp 58478.9173 58483.1345 58488.7358 58848.6754 58848.9740
Display results as a corner plot:
octocorner(model, results, small=true)
Plot RV curve, phase folded curve, and binned residuals:
OctofitterRadialVelocity.rvpostplot(model, results,)
Display RV, PMA, astrometry, relative separation, position angle, and 3D projected views:
octoplot(model, results)