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 PlanetOrbits

We 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)
Example block output

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.0

And plot our simulated astrometry measurments:

fig = Makie.lines(orb_template,)
Makie.scatter!(astrom.table.ra, astrom.table.dec)
fig
Example block output

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[:])
Example block output

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     = 4.75 seconds
Compute duration  = 4.75 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.9999     0.0292    0.0010    785.7538   480.0167    1.000 ⋯
    jitter_1       0.7912     0.5868    0.0164    901.8288   487.6043    1.001 ⋯
       rv0_1       6.7626     0.9467    0.0405    555.8480   549.3671    0.999 ⋯
         plx     100.0000     0.0000       NaN         NaN        NaN       Na ⋯
         b_e       0.7028     0.0135    0.0005    732.1189   680.0086    1.002 ⋯
         b_a       0.9995     0.0115    0.0004    700.9839   560.3749    0.999 ⋯
      b_mass       1.0786     0.1071    0.0042    679.5368   522.9171    1.003 ⋯
         b_i       0.7720     0.0651    0.0027    623.9378   421.6864    1.001 ⋯
        b_Ωy       1.0045     0.0987    0.0039    657.1127   608.4801    0.999 ⋯
        b_Ωx       0.1057     0.0582    0.0023    632.9678   619.7583    1.006 ⋯
        b_ωy       0.7145     0.0969    0.0038    659.7751   497.7267    0.999 ⋯
        b_ωx       0.6961     0.0979    0.0036    733.4854   629.9088    1.001 ⋯
        b_θy       0.7623     0.0747    0.0024   1014.4876   528.4769    0.999 ⋯
        b_θx       0.6575     0.0668    0.0021   1085.3257   787.6426    0.999 ⋯
         b_Ω       0.1051     0.0568    0.0022    663.9083   675.2338    1.005 ⋯
         b_ω       0.7722     0.0929    0.0037    650.3109   715.3534    1.003 ⋯
         b_θ       0.7117     0.0315    0.0010   1016.2578   763.9165    1.000 ⋯
        b_tp   58658.8881   182.4131    6.4223    703.4303   689.3369    0.999 ⋯
                                                               2 columns omitted

Quantiles
  parameters         2.5%        25.0%        50.0%        75.0%        97.5%
      Symbol      Float64      Float64      Float64      Float64      Float64

           M       0.9435       0.9808       1.0000       1.0202       1.0559
    jitter_1       0.0395       0.3310       0.6774       1.1459       2.1501
       rv0_1       5.0059       6.1440       6.7173       7.3722       8.8020
         plx     100.0000     100.0000     100.0000     100.0000     100.0000
         b_e       0.6736       0.6947       0.7034       0.7127       0.7266
         b_a       0.9770       0.9916       0.9998       1.0077       1.0216
      b_mass       0.8779       1.0081       1.0734       1.1392       1.3180
         b_i       0.6359       0.7322       0.7753       0.8171       0.8923
        b_Ωy       0.8266       0.9339       1.0021       1.0644       1.2073
        b_Ωx      -0.0074       0.0687       0.1056       0.1440       0.2219
        b_ωy       0.5464       0.6475       0.7128       0.7760       0.9088
        b_ωx       0.5118       0.6273       0.6951       0.7624       0.8993
        b_θy       0.6274       0.7111       0.7579       0.8108       0.9192
        b_θx       0.5328       0.6109       0.6563       0.6974       0.7939
         b_Ω      -0.0072       0.0696       0.1066       0.1406       0.2193
         b_ω       0.5887       0.7131       0.7718       0.8355       0.9492
         b_θ       0.6487       0.6910       0.7126       0.7339       0.7721
        b_tp   58479.3095   58483.2914   58488.9251   58848.7104   58848.9716

Display results as a corner plot:

octocorner(model, results, small=true)
Example block output

Plot RV curve, phase folded curve, and binned residuals:

OctofitterRadialVelocity.rvpostplot(model, results,)
Example block output

Display RV, PMA, astrometry, relative separation, position angle, and 3D projected views:

octoplot(model, results)
Example block output