GitHub

Fit Relative RV Data

Octofitter includes support for fitting relative radial velocity data — the velocity of a companion measured against its host, rather than the reflex velocity of the host itself.

The convention we adopt is that positive relative radial velocity is the velocity of the companion (exoplanets) minus the velocity of the host (star).

For both stellar reflex RV and relative RVs, you use the same RadialVelocityObs type and just change which ref points to (Barycentre for stellar reflex, another body for relative).

RadialVelocityObs(tab; target=A, ref=Barycentre)   # the star's reflex motion
RadialVelocityObs(tab; target=b, ref=A)            # the companion, relative to the star

To fit relative RV data, start by building the bodies and then the observation:

using Octofitter
using CairoMakie
using Distributions

rv_dat_1 = Table(
    epoch=55000:100:57400,
    rv = [
         -24022.74
        -18571.33
        14221.56
        26076.89
        -459.26
        -26319.26
        -13430.96
        19230.96
        23580.26
        -6786.28
        -27161.78
        -7548.58
        23177.95
        19780.94
        -12738.39
        -26503.74
        -1249.19
        25844.47
        14888.83
        -17986.76
        -24381.49
        5119.22
        27083.2
        9174.18
        -22241.45
    ],
    # Hint! Type as \sigma + <TAB>
    σ_rv= fill(15000.0, 25),
)

See the Basic RV Fit tutorial for examples on how this data can be loaded from a CSV file.

A = Body(
    name="A",
    variables=@variables begin
        # The companion is a test particle here, so the star carries the whole
        # gravitating mass of the orbit.
        mass ~ truncated(Normal(1.2, 0.1), lower=0.1)   # M⊙
    end
)

b = Body(
    name="b",
    about=A,
    variables=@variables begin
        mass = 0.0            # M⊙; a test particle
        a ~ Uniform(0, 10)    # AU
        e ~ Uniform(0.0, 0.5)
        i ~ Sine()
        ω ~ Uniform(0, 2pi)
        Ω ~ Uniform(0, 2pi)
        M0 ~ Uniform(0, 2pi)  # mean anomaly at `epoch`
        epoch = 60000.0       # choose an MJD date near your data
    end
)

rel_rv_obs = RadialVelocityObs(
    rv_dat_1;
    target=b, ref=A,           # <-- this is what makes it *relative* RV
    name="simulated data",
    variables=@variables begin
        jitter ~ LogUniform(0.1, 1000) # m/s
    end
)

The relative RV likelihood does not need an instrument-specific zero point — the two stellar spectra are differenced against each other, so there is nothing to offset. You may still declare an offset variable if your reduction has one. A jitter parameter can be specified in the observation's @variables block, as can parameters for a Gaussian process model of correlated noise (see Fit Gaussian Process). Create one RadialVelocityObs per instrument if you have several, each with its own jitter.

Next, assemble the system:

sys = System(
    name="Example_System",
    bodies=[A, b],
    observations=[rel_rv_obs],
)

model = Octofitter.LogDensityModel(sys)
LogDensityModel for System Example_System of dimension 8 and 25 epochs with fields .ℓπcallback and .∇ℓπcallback

Initialize the model and verify starting point

init_chain = initialize!(model)

octoplot(model, init_chain)
Example block output
using Random
rng = Random.Xoshiro(123)
chain = octofit(rng, model)
Chains MCMC chain (1000×24×1 Array{Float64, 3}):

Iterations        = 1:1:1000
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 3.81 seconds
Compute duration  = 3.81 seconds
parameters        = A_mass, b_a, b_e, b_i, b_ω, b_Ω, b_M0, b_mass, b_epoch, simulated_data_jitter
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, logprior, logpost, tree_depth, numerical_error

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

octoplot(model, chain)
Example block output