Hipparcos Modelling

This tutorial explains how to model Hipparcos IAD data. The first example reproduces the catalog values of position, parallax, and proper motion. The second uses Hipparcos to constrain the mass of a directly imaged planet.

Reproduce Catalog Values

This is the so-called "Nielsen" test from Nielsen et al (2020) and available in Orbitize!.

We start by using a system with a planet with zero mass to fit the straight line motion.

using Octofitter
using Distributions
using CairoMakie

hip_like = Octofitter.HipparcosIADLikelihood(;
    hip_id=21547,
    renormalize=true, # default: true
)

@planet b AbsoluteVisual{KepOrbit} begin
    mass = 0.
    e = 0.
    ω = 0.
    a = 1.
    i = 0
    Ω = 0.
    tp = 0.
end
@system c_Eri_straight_line begin
    M = 1.0 # Host mass not important for this example
    rv = 0.0 # system RV not significant for this example
    plx ~ Uniform(10,100)
    pmra ~ Uniform(-100, 100)
    pmdec ~  Uniform(-100, 100)


    # It is convenient to put a prior of the catalog value +- 10,000 mas on position
    ra_hip_offset_mas ~  Normal(0, 10000)
    dec_hip_offset_mas ~ Normal(0, 10000)
    dec = $hip_like.hip_sol.dedeg + system.ra_hip_offset_mas/60/60/1000
    ra = $hip_like.hip_sol.radeg + system.dec_hip_offset_mas/60/60/1000/cosd(system.dec)

    ref_epoch = Octofitter.hipparcos_catalog_epoch_mjd

end hip_like b

model = Octofitter.LogDensityModel(c_Eri_straight_line)
LogDensityModel for System c_Eri_straight_line of dimension 5 and 104 epochs with fields .ℓπcallback and .∇ℓπcallback

Let's initialize the starting point for the chains to reasonable values

initialize!(model, (;
    plx=34.,
    pmra=44.25,
    pmdec=-64.5,
    ra_hip_offset_mas=0.,
    dec_hip_offset_mas=0.,
))
Chains MCMC chain (1000×18×1 Array{Float64, 3}):

Iterations        = 1:1:1000
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 2.91 seconds
Compute duration  = 2.91 seconds
parameters        = plx, pmra, pmdec, ra_hip_offset_mas, dec_hip_offset_mas, M, rv, dec, ra, ref_epoch, b_mass, b_e, b_ω, b_a, b_i, b_Ω, b_tp
internals         = logpost

Summary Statistics
          parameters         mean       std      mcse   ess_bulk   ess_tail    ⋯
              Symbol      Float64   Float64   Float64    Float64    Float64    ⋯

                 plx      34.0000    0.0000       NaN        NaN        NaN    ⋯
                pmra      44.2500    0.0000       NaN        NaN        NaN    ⋯
               pmdec     -64.5000    0.0000       NaN        NaN        NaN    ⋯
   ra_hip_offset_mas       0.0000    0.0000       NaN        NaN        NaN    ⋯
  dec_hip_offset_mas       0.0000    0.0000       NaN        NaN        NaN    ⋯
                   M       1.0000    0.0000       NaN        NaN        NaN    ⋯
                  rv       0.0000    0.0000       NaN        NaN        NaN    ⋯
                 dec      -2.4734    0.0000    0.0000        NaN        NaN    ⋯
                  ra      69.4004    0.0000       NaN        NaN        NaN    ⋯
           ref_epoch   48348.5625    0.0000       NaN        NaN        NaN    ⋯
              b_mass       0.0000    0.0000       NaN        NaN        NaN    ⋯
                 b_e       0.0000    0.0000       NaN        NaN        NaN    ⋯
                 b_ω       0.0000    0.0000       NaN        NaN        NaN    ⋯
                 b_a       1.0000    0.0000       NaN        NaN        NaN    ⋯
                 b_i       0.0000    0.0000       NaN        NaN        NaN    ⋯
                 b_Ω       0.0000    0.0000       NaN        NaN        NaN    ⋯
                b_tp       0.0000    0.0000       NaN        NaN        NaN    ⋯
                                                               2 columns omitted

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

                 plx      34.0000      34.0000      34.0000      34.0000       ⋯
                pmra      44.2500      44.2500      44.2500      44.2500       ⋯
               pmdec     -64.5000     -64.5000     -64.5000     -64.5000     - ⋯
   ra_hip_offset_mas       0.0000       0.0000       0.0000       0.0000       ⋯
  dec_hip_offset_mas       0.0000       0.0000       0.0000       0.0000       ⋯
                   M       1.0000       1.0000       1.0000       1.0000       ⋯
                  rv       0.0000       0.0000       0.0000       0.0000       ⋯
                 dec      -2.4734      -2.4734      -2.4734      -2.4734       ⋯
                  ra      69.4004      69.4004      69.4004      69.4004       ⋯
           ref_epoch   48348.5625   48348.5625   48348.5625   48348.5625   483 ⋯
              b_mass       0.0000       0.0000       0.0000       0.0000       ⋯
                 b_e       0.0000       0.0000       0.0000       0.0000       ⋯
                 b_ω       0.0000       0.0000       0.0000       0.0000       ⋯
                 b_a       1.0000       1.0000       1.0000       1.0000       ⋯
                 b_i       0.0000       0.0000       0.0000       0.0000       ⋯
                 b_Ω       0.0000       0.0000       0.0000       0.0000       ⋯
                b_tp       0.0000       0.0000       0.0000       0.0000       ⋯
                                                                1 column omitted

We can now sample from the model using Hamiltonian Monte Carlo. This should only take about 15 seconds.

using Pigeons
chain,pt = octofit_pigeons(model, n_rounds=6)
(chain = MCMC chain (64×21×1 Array{Float64, 3}), pt = PT(checkpoint = false, ...))

Plot the posterior values:

octoplot(model,chain,show_astrom=false,show_astrom_time=false)
Example block output

We now visualize the model fit compared to the Hipparcos catalog values:

using LinearAlgebra, StatsBase
fig = Figure(size=(1080,720))
j = i = 1
for prop in (
    (;chain=:ra, hip=:radeg, hip_err=:e_ra),
    (;chain=:dec, hip=:dedeg, hip_err=:e_de),
    (;chain=:plx, hip=:plx, hip_err=:e_plx),
    (;chain=:pmra, hip=:pm_ra, hip_err=:e_pmra),
    (;chain=:pmdec, hip=:pm_de, hip_err=:e_pmde)
)
    global i, j, ax
    ax = Axis(
        fig[j,i],
        xlabel=string(prop.chain),
    )
    i+=1
    if i > 3
        j+=1
        i = 1
    end
    unc = hip_like.hip_sol[prop.hip_err]
    if prop.chain == :ra
        unc /= 60*60*1000 * cosd(hip_like.hip_sol.dedeg)
    end
    if prop.chain == :dec
        unc /= 60*60*1000
    end
    if prop.hip == :zero
        n = Normal(0, unc)
    else
        mu = hip_like.hip_sol[prop.hip]
        n = Normal(mu, unc)
    end
    n0,n1=quantile.(n,(1e-4, 1-1e-4))
    nxs = range(n0,n1,length=200)
    h = fit(Histogram, chain[prop.chain][:], nbins=55)
    h = normalize(h, mode=:pdf)
    barplot!(ax, (h.edges[1][1:end-1] .+ h.edges[1][2:end])./2, h.weights, gap=0, color=:red, label="posterior")
    lines!(ax, nxs, pdf.(n,nxs), label="Hipparcos Catalog", color=:black, linewidth=2)
end
Legend(fig[i-1,j+1],ax,tellwidth=false)
fig
Example block output

Constrain Planet Mass

We now allow the planet to have a non zero mass and have free orbit. We start by specifying relative astrometry data on the planet, collated by Jason Wang and co. on whereistheplanet.com.

astrom_like1 = PlanetRelAstromLikelihood(
    (;epoch=57009.1, sep=454.24,  σ_sep=1.88, pa=2.98835, σ_pa=0.00401426),
    (;epoch=57052.1, sep=451.81,  σ_sep=2.06, pa=2.96723, σ_pa=0.00453786),
    (;epoch=57053.1, sep=456.8 ,  σ_sep=2.57, pa=2.97038, σ_pa=0.00523599),
    (;epoch=57054.3, sep=461.5 ,  σ_sep=23.9 ,pa=2.97404, σ_pa=0.0523599 ,),
    (;epoch=57266.4, sep=455.1 ,  σ_sep=2.23, pa=2.91994, σ_pa=0.00453786),
    (;epoch=57332.2, sep=452.88,  σ_sep=5.41, pa=2.89934, σ_pa=0.00994838),
    (;epoch=57374.2, sep=455.91,  σ_sep=6.23, pa=2.89131, σ_pa=0.00994838),
    (;epoch=57376.2, sep=455.01,  σ_sep=3.03, pa=2.89184, σ_pa=0.00750492),
    (;epoch=57415.0, sep=454.46,  σ_sep=6.03, pa=2.8962 , σ_pa=0.00890118),
    (;epoch=57649.4, sep=454.81,  σ_sep=2.02, pa=2.82394, σ_pa=0.00453786),
    (;epoch=57652.4, sep=451.43,  σ_sep=2.67, pa=2.82272, σ_pa=0.00541052),
    (;epoch=57739.1, sep=449.39,  σ_sep=2.15, pa=2.79357, σ_pa=0.00471239),
    (;epoch=58068.3, sep=447.54,  σ_sep=3.02, pa=2.70927, σ_pa=0.00680678),
    (;epoch=58442.2, sep=434.22,  σ_sep=2.01, pa=2.61171, σ_pa=0.00401426),
)
PlanetRelAstromLikelihood Table with 5 columns and 14 rows:
      epoch    sep     σ_sep  pa       σ_pa
    ┌────────────────────────────────────────────
 1  │ 57009.1  454.24  1.88   2.98835  0.00401426
 2  │ 57052.1  451.81  2.06   2.96723  0.00453786
 3  │ 57053.1  456.8   2.57   2.97038  0.00523599
 4  │ 57054.3  461.5   23.9   2.97404  0.0523599
 5  │ 57266.4  455.1   2.23   2.91994  0.00453786
 6  │ 57332.2  452.88  5.41   2.89934  0.00994838
 7  │ 57374.2  455.91  6.23   2.89131  0.00994838
 8  │ 57376.2  455.01  3.03   2.89184  0.00750492
 9  │ 57415.0  454.46  6.03   2.8962   0.00890118
 10 │ 57649.4  454.81  2.02   2.82394  0.00453786
 11 │ 57652.4  451.43  2.67   2.82272  0.00541052
 12 │ 57739.1  449.39  2.15   2.79357  0.00471239
 13 │ 58068.3  447.54  3.02   2.70927  0.00680678
 14 │ 58442.2  434.22  2.01   2.61171  0.00401426

We specify our full model:

@planet b AbsoluteVisual{KepOrbit} begin
    a ~ truncated(Normal(10,1),lower=0.1)
    e ~ Uniform(0,0.99)
    ω ~ Uniform(0, 2pi)
    i ~ Sine()
    Ω ~ Uniform(0, 2pi)
    θ ~ Uniform(0, 2pi)
    tp = θ_at_epoch_to_tperi(system,b,58442.2)
    mass = system.M_sec
end astrom_like1

@system cEri begin
    M_pri ~ truncated(Normal(1.75,0.05), lower=0.03) # Msol
    M_sec ~ LogUniform(0.1, 100) # MJup
    M = system.M_pri + system.M_sec*Octofitter.mjup2msol # Msol

    rv =  12.60e3 # m/s
    plx ~ Uniform(20,40)
    pmra ~ Uniform(-100, 100)
    pmdec ~  Uniform(-100, 100)

    # It is convenient to put a prior of the catalog value +- 1000 mas on position
    ra_hip_offset_mas ~  Normal(0, 1000)
    dec_hip_offset_mas ~ Normal(0, 1000)
    dec = $hip_like.hip_sol.dedeg + system.ra_hip_offset_mas/60/60/1000
    ra = $hip_like.hip_sol.radeg + system.dec_hip_offset_mas/60/60/1000/cos(system.dec)

    ref_epoch = Octofitter.hipparcos_catalog_epoch_mjd
end hip_like b

model = Octofitter.LogDensityModel(cEri)
LogDensityModel for System cEri of dimension 13 and 118 epochs with fields .ℓπcallback and .∇ℓπcallback

Initialize the starting points, and confirm the data are entered correcly:

init_chain = initialize!(model, (;
    plx=34.,
    pmra=44.25,
    pmdec=-64.5,
    ra_hip_offset_mas=0.,
    dec_hip_offset_mas=0.,
))
octoplot(model, init_chain)
Example block output

Now we sample:

using Pigeons
chain,pt = octofit_pigeons(model, n_rounds=8, explorer=SliceSampler())
chain
Chains MCMC chain (256×24×1 Array{Float64, 3}):

Iterations        = 1:1:256
Number of chains  = 1
Samples per chain = 256
Wall duration     = 463.58 seconds
Compute duration  = 463.58 seconds
parameters        = M_pri, M_sec, plx, pmra, pmdec, ra_hip_offset_mas, dec_hip_offset_mas, M, rv, dec, ra, ref_epoch, b_a, b_e, b_ω, b_i, b_Ω, b_θ, b_tp, b_mass
internals         = loglike, logpost, logprior, pigeons_logpotential

Summary Statistics
          parameters         mean        std      mcse   ess_bulk   ess_tail   ⋯
              Symbol      Float64    Float64   Float64    Float64    Float64   ⋯

               M_pri       1.7390     0.0473    0.0039   148.7037   137.3135   ⋯
               M_sec       9.8555    17.7061    1.3836   193.9058   108.8616   ⋯
                 plx      33.9609     0.3415    0.0265   166.7147   165.4535   ⋯
                pmra      44.4331     0.5267    0.0406   186.8091   167.0159   ⋯
               pmdec     -64.4064     0.3122    0.0225   210.1332   184.4831   ⋯
   ra_hip_offset_mas      -2.2972     4.1344    0.3209   171.2353   108.8616   ⋯
  dec_hip_offset_mas      -0.2278     0.5002    0.0979    78.7606    89.7199   ⋯
                   M       1.7484     0.0511    0.0041   160.1364   113.6174   ⋯
                  rv   12600.0000     0.0000       NaN        NaN        NaN   ⋯
                 dec      -2.4734     0.0000    0.0000   171.2353   108.8616   ⋯
                  ra      69.4004     0.0000    0.0000    78.7606    89.7199   ⋯
           ref_epoch   48348.5625     0.0000       NaN        NaN        NaN   ⋯
                 b_a      10.3245     0.5027    0.0658    55.8552    62.2894   ⋯
                 b_e       0.5764     0.0551    0.0052   160.6194    86.2594   ⋯
                 b_ω       1.8610     1.1776    0.1432   169.1137    76.4464   ⋯
                 b_i       2.4499     0.0695    0.0082    62.8152    67.3279   ⋯
                 b_Ω       1.5590     1.3227    0.1418   166.6806    98.9104   ⋯
          ⋮                ⋮           ⋮          ⋮         ⋮          ⋮       ⋱
                                                    2 columns and 3 rows omitted

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

               M_pri       1.6551       1.7026       1.7413       1.7704       ⋯
               M_sec       0.1186       0.5503       1.7445       7.9406       ⋯
                 plx      33.3021      33.7332      33.9505      34.1827       ⋯
                pmra      43.7083      44.0545      44.3348      44.6960       ⋯
               pmdec     -65.1464     -64.6099     -64.3859     -64.1823     - ⋯
   ra_hip_offset_mas     -15.3784      -1.8820      -0.4571      -0.1770       ⋯
  dec_hip_offset_mas      -1.5949      -0.3369      -0.1007       0.0641       ⋯
                   M       1.6629       1.7113       1.7468       1.7813       ⋯
                  rv   12600.0000   12600.0000   12600.0000   12600.0000   126 ⋯
                 dec      -2.4734      -2.4734      -2.4734      -2.4734       ⋯
                  ra      69.4004      69.4004      69.4004      69.4004       ⋯
           ref_epoch   48348.5625   48348.5625   48348.5625   48348.5625   483 ⋯
                 b_a       9.3849       9.9990      10.3327      10.6806       ⋯
                 b_e       0.4518       0.5564       0.5873       0.6073       ⋯
                 b_ω       0.9434       1.0863       1.5106       1.8898       ⋯
                 b_i       2.3117       2.4119       2.4477       2.4859       ⋯
                 b_Ω       0.1785       0.5082       1.2699       2.0023       ⋯
          ⋮                ⋮            ⋮            ⋮            ⋮            ⋱
                                                     1 column and 3 rows omitted
octoplot(model, chain, show_mass=true)
Example block output

We see that we constrained both the orbit and the parallax. The mass is not strongly constrained by Hipparcos.