Fitting Interferometric Observables

In this tutorial, we fit a planet & orbit model to a sequence of interferometric observations. Closure phases and squared visibilities are supported.

We load the observations in OI-FITS format and model them as a point source orbiting a star.

Note

Interferometer modelling is supported in Octofitter via the extension package OctofitterInterferometry. To install it, run pkg> add http://github.com/sefffal/Octofitter.jl:OctofitterInterferometry

using Octofitter
using OctofitterInterferometry
using Distributions
using CairoMakie
using PairPlots
┌ Warning: Module Octofitter with build ID ffffffff-ffff-ffff-0000-008764614b64 is missing from the cache.
│ This may mean Octofitter [daf3887e-d01a-44a1-9d7e-98f15c5d69c9] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:1948

Download simulated JWST AMI observations from our examples folder on GitHub:

download("https://github.com/sefffal/Octofitter.jl/raw/main/examples/AMI_data/Sim_data_2023_1_.oifits", "Sim_data_2023_1_.oifits")
download("https://github.com/sefffal/Octofitter.jl/raw/main/examples/AMI_data/Sim_data_2023_2_.oifits", "Sim_data_2023_2_.oifits")
download("https://github.com/sefffal/Octofitter.jl/raw/main/examples/AMI_data/Sim_data_2024_1_.oifits", "Sim_data_2024_1_.oifits")
"Sim_data_2024_1_.oifits"

Create the likelihood object:

vis_like = InterferometryLikelihood(
    (; filename="Sim_data_2023_1_.oifits", epoch=mjd("2023-06-01"), band=:F480M, use_vis2=false),
    (; filename="Sim_data_2023_2_.oifits", epoch=mjd("2023-08-15"), band=:F480M, use_vis2=false),
    (; filename="Sim_data_2024_1_.oifits", epoch=mjd("2024-06-01"), band=:F480M, use_vis2=false),
)
OctofitterInterferometry.InterferometryLikelihood Table with 14 columns and 3 rows:
     filename              epoch    band   use_vis2  u                     ⋯
   ┌────────────────────────────────────────────────────────────────────────
 1 │ Sim_data_2023_1_.oi…  60096.0  F480M  false     [-3.46641e5; 6.7596…  ⋯
 2 │ Sim_data_2023_2_.oi…  60171.0  F480M  false     [-3.46641e5; 6.7596…  ⋯
 3 │ Sim_data_2024_1_.oi…  60462.0  F480M  false     [-3.46641e5; 6.7596…  ⋯

Plot the closure phases:

fig = Makie.Figure()
ax = Axis(
    fig[1,1],
    xlabel="index",
    ylabel="closure phase",
)
Makie.stem!(
    vis_like.table.cps_data[1][:],
    label="epoch 1",
)
Makie.stem!(
    vis_like.table.cps_data[2][:],
    label="epoch 2"
)
Makie.stem!(
    vis_like.table.cps_data[3][:],
    label="epoch 3"
)
Makie.Legend(fig[1,2], ax)
fig
Example block output
@planet b Visual{KepOrbit} begin
    a ~ truncated(Normal(2,0.1), lower=0)
    e ~ truncated(Normal(0, 0.05),lower=0, upper=1.0)
    i ~ Sine()
    ω ~ UniformCircular()
    Ω ~ UniformCircular()

    # Our prior on the planet's photometry
    # 0 +- 10% of stars brightness (assuming this is unit of data files)
    F480M ~ truncated(Normal(0, 0.1),lower=0)

    θ ~ UniformCircular()
    tp = θ_at_epoch_to_tperi(system,b,60171)  # reference epoch for θ. Choose an MJD date near your data.
end

@system Tutoria begin
    M ~ truncated(Normal(1.5, 0.01), lower=0)
    plx ~ truncated(Normal(100., 0.1), lower=0)
end vis_like b
System model Tutoria
Derived:
  
Priors:
  M	Truncated(Distributions.Normal{Float64}(μ=1.5, σ=0.01); lower=0.0)
  plx	Truncated(Distributions.Normal{Float64}(μ=100.0, σ=0.1); lower=0.0)
Planet b
Derived:
  ω, Ω, θ, tp, 
Priors:
  a	Truncated(Distributions.Normal{Float64}(μ=2.0, σ=0.1); lower=0.0)
  e	Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.05); lower=0.0, upper=1.0)
  i	Sine()
  ωy	Distributions.Normal{Float64}(μ=0.0, σ=1.0)
  ωx	Distributions.Normal{Float64}(μ=0.0, σ=1.0)
  Ωy	Distributions.Normal{Float64}(μ=0.0, σ=1.0)
  Ωx	Distributions.Normal{Float64}(μ=0.0, σ=1.0)
  F480M	Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.1); lower=0.0)
  θy	Distributions.Normal{Float64}(μ=0.0, σ=1.0)
  θx	Distributions.Normal{Float64}(μ=0.0, σ=1.0)
Octofitter.UnitLengthPrior{:ωy, :ωx}: √(ωy^2+ωx^2) ~ LogNormal(log(1), 0.02)
Octofitter.UnitLengthPrior{:Ωy, :Ωx}: √(Ωy^2+Ωx^2) ~ LogNormal(log(1), 0.02)
Octofitter.UnitLengthPrior{:θy, :θx}: √(θy^2+θx^2) ~ LogNormal(log(1), 0.02)


OctofitterInterferometry.InterferometryLikelihood Table with 14 columns and 3 rows:
     filename              epoch    band   use_vis2  u                     ⋯
   ┌────────────────────────────────────────────────────────────────────────
 1 │ Sim_data_2023_1_.oi…  60096.0  F480M  false     [-3.46641e5; 6.7596…  ⋯
 2 │ Sim_data_2023_2_.oi…  60171.0  F480M  false     [-3.46641e5; 6.7596…  ⋯
 3 │ Sim_data_2024_1_.oi…  60462.0  F480M  false     [-3.46641e5; 6.7596…  ⋯

Create the model object and run octofit:

model = Octofitter.LogDensityModel(Tutoria)

using Random
rng = Xoshiro(0) # Optional seeded random number generator.
results = octofit(rng, model,  adaptation=2000,iterations=5000)
Chains MCMC chain (5000×29×1 Array{Float64, 3}):

Iterations        = 1:1:5000
Number of chains  = 1
Samples per chain = 5000
Wall duration     = 13.99 seconds
Compute duration  = 13.99 seconds
parameters        = M, plx, b_a, b_e, b_i, b_ωy, b_ωx, b_Ωy, b_Ωx, b_F480M, 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      rh ⋯
      Symbol      Float64    Float64   Float64     Float64     Float64   Float ⋯

           M       1.4994     0.0099    0.0002   4241.1553   3392.6981    0.99 ⋯
         plx      99.9993     0.1005    0.0011   8447.9477   3367.6732    1.00 ⋯
         b_a       2.0730     0.0555    0.0013   1895.7368   3411.5071    0.99 ⋯
         b_e       0.0352     0.0273    0.0006   1847.3337   1950.4820    1.00 ⋯
         b_i       0.5956     0.1093    0.0022   2509.8985   2703.9002    0.99 ⋯
        b_ωy       0.0435     0.6671    0.0181   1438.6176   3424.8590    1.00 ⋯
        b_ωx      -0.2196     0.7260    0.0565    200.4483   2260.6931    1.00 ⋯
        b_Ωy      -0.4455     0.6021    0.1305     41.8775    105.4821    1.00 ⋯
        b_Ωx      -0.4000     0.5449    0.1144     39.0518     98.3379    1.00 ⋯
     b_F480M       0.0005     0.0001    0.0000   3518.9199   2665.9701    1.00 ⋯
        b_θy       0.4997     0.0610    0.0011   3250.6819   2670.3905    0.99 ⋯
        b_θx      -0.8722     0.0885    0.0014   4081.2970   3318.0875    0.99 ⋯
         b_ω      -0.3730     1.6922    0.0992    500.8647   3204.2509    1.00 ⋯
         b_Ω      -1.7499     1.3540    0.2740     48.9449    223.1968    1.00 ⋯
         b_θ      -1.0506     0.0382    0.0007   3081.4607   2789.4043    1.00 ⋯
        b_tp   59759.5298   225.8954    5.8285   1643.6982   2088.3700    0.99 ⋯
                                                               2 columns omitted

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

           M       1.4804       1.4928       1.4990       1.5059       1.5193
         plx      99.8078      99.9291      99.9983     100.0699     100.1918
         b_a       1.9605       2.0370       2.0743       2.1085       2.1801
         b_e       0.0012       0.0136       0.0298       0.0506       0.1012
         b_i       0.3820       0.5214       0.5966       0.6664       0.8133
        b_ωy      -1.0480      -0.5606       0.0838       0.6482       1.0531
        b_ωx      -1.0996      -0.8750      -0.4570       0.4721       1.0515
        b_Ωy      -1.0424      -0.8243      -0.6801      -0.4513       0.9448
        b_Ωx      -0.9973      -0.7608      -0.6033      -0.3071       0.8675
     b_F480M       0.0003       0.0004       0.0005       0.0005       0.0006
        b_θy       0.3897       0.4565       0.4974       0.5395       0.6269
        b_θx      -1.0556      -0.9283      -0.8687      -0.8112      -0.7109
         b_ω      -2.9368      -1.6655      -0.7787       0.9998       2.9057
         b_Ω      -2.9186      -2.5282      -2.3180      -2.0157       1.1170
         b_θ      -1.1281      -1.0757      -1.0502      -1.0249      -0.9773
        b_tp   59293.0295   59605.3737   59794.5015   59923.1404   60134.5723

Examine the recovered photometry posterior:

hist(results[:b_F480M][:], axis=(;xlabel="F480M"))
Example block output

Determine the significance of the detection:

using Statistics
phot = results[:b_F480M][:]
snr = mean(phot)/std(phot)
7.236867536007492

Plot the resulting orbit:

octoplot(model, results)
Example block output

Plot position at each epoch:

using PlanetOrbits
els = Octofitter.construct_elements(results,:b,:);
fig = Makie.Figure()
ax = Makie.Axis(
    fig[1,1],
    autolimitaspect = 1,
    xreversed=true,
    xlabel="ΔR.A. (mas)",
    ylabel="ΔDec. (mas)",
)
for epoch in vis_like.table.epoch
    Makie.scatter!(
        ax,
        raoff.(els, epoch)[:],
        decoff.(els, epoch)[:],
        label=string(mjd2date(epoch)),
        markersize=1.5,
    )
end
Makie.Legend(fig[1,2], ax, "date")
fig
Example block output

We can use PairPlots.jl to create a contour plot of positions at all three epochs:

els = Octofitter.construct_elements(results,:b,:);
fig = pairplot(
    [
        (;
                ra=raoff.(els, epoch)[:],
                dec=decoff.(els, epoch)[:],
        )=>(PairPlots.Contourf(),)
        for epoch in vis_like.table.epoch
    ]...,
    bodyaxis=(;width=400,height=400),
    axis=(;
        ra=(;reversed=true, lims=(;low=250,high=-250,)),
        dec=(;lims=(;low=-250,high=250,)),
    ),
    labels=Dict(:ra=>"ra offset [mas]", :dec=>"dec offset [mas]"),
)
Makie.scatter!(fig.content[1], [0],[0],marker='⭐', markersize=30, color=:black)
fig
Example block output

Finally we can examine the joint photometry and orbit posterior as a corner plot:

using PairPlots
using CairoMakie: Makie
octocorner(model, results)
Example block output