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.41 seconds
Compute duration = 8.41 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.9944 0.0428 0.0030 205.4638 250.2706 0.9991 ⋯
jitter1 0.4965 0.3468 0.0182 285.0751 187.3742 0.9989 ⋯
plx 100.0000 0.0000 NaN NaN NaN NaN ⋯
b_e 0.2057 0.1669 0.0302 21.7839 46.9137 1.0443 ⋯
b_a 1.9406 0.5528 0.0690 52.8409 81.3475 1.0168 ⋯
b_mass 0.9753 0.6064 0.0356 233.2062 121.7101 1.0006 ⋯
b_i 1.1748 0.0800 0.0145 43.9264 39.1719 1.0242 ⋯
b_Ωy -0.7082 0.1516 0.0204 54.4391 80.5919 1.0164 ⋯
b_Ωx -0.6923 0.1626 0.0229 55.5687 88.3022 1.0102 ⋯
b_ωy -0.1896 0.7656 0.1709 24.4334 169.4649 1.0136 ⋯
b_ωx 0.2077 0.5917 0.0828 63.6769 219.8237 1.0143 ⋯
b_θy -0.9928 0.0995 0.0060 271.7965 224.2204 0.9977 ⋯
b_θx -0.1676 0.0179 0.0010 305.7911 272.8527 0.9989 ⋯
b_Ω -2.3685 0.1953 0.0300 45.2494 52.3727 1.0157 ⋯
b_ω 0.8517 1.8521 0.2530 77.4492 157.5793 1.0090 ⋯
b_θ -2.9744 0.0054 0.0003 392.1820 292.6372 1.0008 ⋯
b_tp 58322.1720 454.0348 86.4111 24.5771 109.5364 1.0285 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
M 0.9119 0.9664 0.9929 1.0249 1.0778
jitter1 0.0333 0.2180 0.4438 0.7140 1.3285
plx 100.0000 100.0000 100.0000 100.0000 100.0000
b_e 0.0056 0.0673 0.1571 0.3040 0.5965
b_a 1.0997 1.5357 1.8474 2.2865 3.0816
b_mass 0.0711 0.5349 0.8690 1.3005 2.3824
b_i 0.9891 1.1386 1.1888 1.2302 1.2958
b_Ωy -1.0021 -0.8144 -0.7162 -0.5960 -0.4345
b_Ωx -0.9753 -0.7965 -0.7041 -0.6084 -0.3213
b_ωy -1.0705 -0.8533 -0.5783 0.6322 1.0849
b_ωx -1.0127 -0.1858 0.3683 0.6390 1.0108
b_θy -1.1905 -1.0601 -0.9893 -0.9218 -0.7984
b_θx -0.2046 -0.1796 -0.1680 -0.1545 -0.1322
b_Ω -2.8032 -2.4817 -2.3622 -2.2216 -2.0219
b_ω -2.9460 -0.3132 1.3414 2.5283 2.9325
b_θ -2.9854 -2.9782 -2.9742 -2.9708 -2.9643
b_tp 57249.2554 58112.3950 58467.9316 58658.5236 58799.1431
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)