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 = 10.67 seconds
Compute duration = 10.67 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 rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
M 0.9987 0.0405 0.0042 89.5802 86.3339 1.0000 ⋯
jitter1 0.4720 0.4121 0.0246 196.3376 118.1601 1.0061 ⋯
plx 100.0000 0.0000 NaN NaN NaN NaN ⋯
b_e 0.6125 0.2042 0.0496 20.1891 89.9835 1.0093 ⋯
b_a 1.2901 0.4209 0.1087 18.4180 52.8151 1.0617 ⋯
b_mass 0.6965 0.5193 0.0576 67.1341 67.3337 1.0041 ⋯
b_i 0.8491 0.2614 0.0508 30.2233 113.1292 1.0088 ⋯
b_Ωy -0.4265 0.6047 0.1533 15.9817 40.2910 1.0698 ⋯
b_Ωx -0.2448 0.6577 0.4176 2.1711 15.2546 1.4728 ⋯
b_ωy -0.3695 0.6343 0.2042 12.6916 45.1658 1.1257 ⋯
b_ωx -0.0341 0.6763 0.4920 2.2693 30.6540 1.4653 ⋯
b_θy -1.0036 0.1011 0.0084 146.0911 140.8038 1.0042 ⋯
b_θx -0.1688 0.0182 0.0014 177.5505 220.8101 0.9978 ⋯
b_Ω -1.2444 1.8490 0.5994 20.6788 24.1018 1.1569 ⋯
b_ω 0.3669 2.1830 1.2574 7.8361 110.5506 1.2513 ⋯
b_θ -2.9749 0.0062 0.0004 199.6859 172.1110 1.0198 ⋯
b_tp 58745.8956 52.9554 12.1529 21.1729 75.5356 1.0048 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
M 0.9205 0.9702 0.9983 1.0245 1.0786
jitter1 0.0086 0.1643 0.3446 0.7316 1.4851
plx 100.0000 100.0000 100.0000 100.0000 100.0000
b_e 0.1782 0.4497 0.6560 0.7992 0.8589
b_a 0.8136 0.9704 1.2254 1.4867 2.3750
b_mass 0.0221 0.2842 0.6272 1.0042 1.9895
b_i 0.2249 0.7341 0.9185 1.0473 1.1539
b_Ωy -1.0637 -0.8758 -0.7055 0.1167 0.8691
b_Ωx -1.0165 -0.7023 -0.4946 -0.0090 1.1479
b_ωy -1.1148 -0.9118 -0.6445 0.1502 0.8017
b_ωx -1.0372 -0.6630 0.0236 0.5583 1.0463
b_θy -1.2087 -1.0720 -0.9970 -0.9350 -0.8138
b_θx -0.2047 -0.1830 -0.1682 -0.1560 -0.1345
b_Ω -2.9909 -2.5845 -2.3138 -0.0845 2.9290
b_ω -3.0376 -1.6795 1.0611 2.5758 3.0629
b_θ -2.9865 -2.9793 -2.9754 -2.9705 -2.9631
b_tp 58617.9434 58717.2253 58767.6680 58786.6890 58795.8871
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)