GitHub

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 =58829-40
)
Makie.lines(orb_template,axis=(;autolimitaspect=1))
Example block output

Sample position and store as relative astrometry measurements:

epochs = [58849,58852,58858,58890]
astrom_dat = 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)),
    cor=fill(0.0, size(epochs))
)

astrom = PlanetRelAstromObs(
    astrom_dat,
    name = "simulated",
    variables = @variables begin
        # Fixed values for this example - could be free variables:
        jitter = 0        # mas [could use: jitter ~ Uniform(0, 10)]
        northangle = 0    # radians [could use: northangle ~ Normal(0, deg2rad(1))]
        platescale = 1    # relative [could use: platescale ~ truncated(Normal(1, 0.01), lower=0)]
    end
)
PlanetRelAstromObs Table with 6 columns and 4 rows:
     epoch  ra        dec       σ_ra  σ_dec  cor
   ┌────────────────────────────────────────────
 1 │ 58849  -18.3017  -108.907  1.0   1.0    0.0
 2 │ 58852  -21.1181  -111.432  1.0   1.0    0.0
 3 │ 58858  -26.6349  -115.889  1.0   1.0    0.0
 4 │ 58890  -53.1334  -129.043  1.0   1.0    0.0

And plot our simulated astrometry measurments:

fig = Makie.lines(orb_template,axis=(;autolimitaspect=1))
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 .+ range(0,step=1.5, length=20)
planet_sim_mass = 0.001 # solar masses here


rvlike = MarginalizedStarAbsoluteRVObs(
    Table(
        epoch=epochs,
        rv=radvel.(orb_template, epochs, planet_sim_mass) .+ 150,
        σ_rv=fill(5.0, size(epochs)),
    ),
    name="inst1",
    variables=@variables begin
        jitter ~ LogUniform(0.1, 100) # m/s
    end
)

epochs = 58949 .+ range(0,step=1.5, length=20)

rvlike2 = MarginalizedStarAbsoluteRVObs(
    Table(
        epoch=epochs,
        rv=radvel.(orb_template, epochs, planet_sim_mass) .- 150,
        σ_rv=fill(5.0, size(epochs)),
    ),
    name="inst2",
    variables=@variables begin
        jitter ~ LogUniform(0.1, 100) # m/s
    end
)

fap = Makie.scatter(rvlike.table.epoch[:], rvlike.table.rv[:])
Makie.scatter!(rvlike2.table.epoch[:], rvlike2.table.rv[:])
fap
Example block output

Now specify model and fit it:

planet_b = Planet(
    name="b",
    basis=Visual{KepOrbit},
    observations=[astrom],
    variables=@variables begin
        e ~ Uniform(0,0.999999)
        a ~ truncated(Normal(1, 1),lower=0.1)
        mass ~ truncated(Normal(1, 1), lower=0.)
        i ~ Sine()
        M = system.M
        Ω ~ UniformCircular()
        ω ~ UniformCircular()
        θ ~ UniformCircular()
        tp = θ_at_epoch_to_tperi(θ, 58849.0; M, e, a, i, ω, Ω)  # reference epoch for θ. Choose an MJD date near your data.
    end
)

sys = System(
    name="test",
    companions=[planet_b],
    observations=[rvlike, rvlike2],
    variables=@variables begin
        M ~ truncated(Normal(1, 0.04),lower=0.1) # (Baines & Armstrong 2011).
        plx = 100.0
    end
)

model = Octofitter.LogDensityModel(sys)

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×35×1 Array{Float64, 3}):

Iterations        = 1:1:400
Number of chains  = 1
Samples per chain = 400
Wall duration     = 12.08 seconds
Compute duration  = 12.08 seconds
parameters        = M, plx, inst1_jitter, inst2_jitter, b_e, b_a, b_mass, b_i, b_Ωx, b_Ωy, b_ωx, b_ωy, b_θx, b_θy, b_Ω, b_ω, b_θ, b_M, b_tp, b_simulated_jitter, b_simulated_northangle, b_simulated_platescale
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

Use `describe(chains)` for summary statistics and quantiles.

Display results as a corner plot:

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

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

Octofitter.rvpostplot(model, results)
Example block output

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

octoplot(model, results)
Example block output