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 = 8.94 seconds
Compute duration = 8.94 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.0382 0.0023 286.0594 213.0921 1.0213 ⋯
jitter1 0.4587 0.3340 0.0177 262.2884 177.6767 0.9981 ⋯
plx 100.0000 0.0000 NaN NaN NaN NaN ⋯
b_e 0.2466 0.1835 0.0259 51.4972 108.2819 1.0094 ⋯
b_a 1.7363 0.5336 0.0746 48.2024 79.5997 1.0034 ⋯
b_mass 0.9302 0.5699 0.0509 106.3104 112.7010 1.0014 ⋯
b_i 1.1415 0.0902 0.0135 47.3333 80.0045 1.0050 ⋯
b_Ωy -0.6885 0.1580 0.0199 66.8446 208.6899 1.0003 ⋯
b_Ωx -0.7112 0.1413 0.0177 64.7448 135.5576 1.0009 ⋯
b_ωy -0.5615 0.5735 0.0742 74.8850 95.2298 0.9999 ⋯
b_ωx 0.3224 0.5141 0.0835 41.7226 74.5113 0.9979 ⋯
b_θy -0.9849 0.0882 0.0055 263.4601 244.5979 0.9994 ⋯
b_θx -0.1657 0.0157 0.0009 337.1740 208.3242 0.9997 ⋯
b_Ω -2.3379 0.1908 0.0261 55.0234 144.8706 0.9985 ⋯
b_ω 1.2052 2.0692 0.2505 142.2506 202.3846 1.0008 ⋯
b_θ -2.9749 0.0060 0.0003 421.6103 161.9141 1.0052 ⋯
b_tp 58504.9357 273.0580 33.8677 57.3736 92.7018 0.9986 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
M 0.9281 0.9738 0.9966 1.0251 1.0742
jitter1 0.0171 0.1991 0.4134 0.6220 1.2334
plx 100.0000 100.0000 100.0000 100.0000 100.0000
b_e 0.0112 0.0894 0.1996 0.3986 0.6085
b_a 1.0505 1.3372 1.6076 2.0500 2.8748
b_mass 0.0881 0.4510 0.9011 1.2973 2.2032
b_i 0.9467 1.0895 1.1564 1.2072 1.2812
b_Ωy -0.9918 -0.8150 -0.6829 -0.5755 -0.3978
b_Ωx -0.9453 -0.8120 -0.7127 -0.6096 -0.4279
b_ωy -1.1430 -0.9221 -0.7945 -0.4817 0.9620
b_ωx -0.9210 0.0193 0.4703 0.6946 0.9804
b_θy -1.1501 -1.0447 -0.9893 -0.9193 -0.8231
b_θx -0.1990 -0.1751 -0.1658 -0.1559 -0.1358
b_Ω -2.6734 -2.4911 -2.3416 -2.1886 -2.0043
b_ω -3.0596 0.3082 2.3689 2.6297 3.0519
b_θ -2.9869 -2.9790 -2.9743 -2.9705 -2.9637
b_tp 57722.6058 58415.9632 58580.3389 58694.7337 58773.0250
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)