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(0,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

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

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     = 377.55 seconds
Compute duration  = 377.55 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.7515     0.0509    0.0034   224.6770   187.7860   ⋯
               M_sec       9.6336    12.2401    3.7490     8.0793    22.2639   ⋯
                 plx      33.9395     0.3376    0.0329   111.7684   116.3528   ⋯
                pmra      44.4192     0.4373    0.0590    50.4479   155.6490   ⋯
               pmdec     -64.3261     0.2765    0.0218   159.0450   132.9624   ⋯
   ra_hip_offset_mas      -2.2463     2.8974    0.9343    17.4322    22.2639   ⋯
  dec_hip_offset_mas      -0.3274     0.6096    0.0491   172.1545   123.0522   ⋯
                   M       1.7607     0.0510    0.0035   209.3220   225.4685   ⋯
                  rv   12600.0000     0.0000       NaN        NaN        NaN   ⋯
                 dec      -2.4734     0.0000    0.0000    17.4322    22.2639   ⋯
                  ra      69.4004     0.0000    0.0000   172.1545   123.0522   ⋯
           ref_epoch   48348.5625     0.0000       NaN        NaN        NaN   ⋯
                 b_a      10.2443     0.5444    0.1756    10.6619    85.3130   ⋯
                 b_e       0.5975     0.0560    0.0198    10.8520    98.1379   ⋯
                 b_ω       2.5531     1.5358    0.2058   192.4059    88.4120   ⋯
                 b_i       2.4749     0.0696    0.0247     6.2835    32.1471   ⋯
                 b_Ω       2.3924     1.6482    0.1813   211.4972    29.7334   ⋯
          ⋮                ⋮           ⋮          ⋮         ⋮          ⋮       ⋱
                                                    2 columns and 3 rows omitted

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

               M_pri       1.6631       1.7143       1.7509       1.7873       ⋯
               M_sec       0.1187       0.5229       4.6815      16.3014       ⋯
                 plx      33.3014      33.7282      33.9428      34.1771       ⋯
                pmra      43.6068      44.1447      44.4122      44.6593       ⋯
               pmdec     -64.8232     -64.5007     -64.3496     -64.1563     - ⋯
   ra_hip_offset_mas      -7.3031      -3.9142      -0.9929      -0.1693       ⋯
  dec_hip_offset_mas      -1.7490      -0.5585      -0.1780       0.0434       ⋯
                   M       1.6715       1.7248       1.7605       1.7962       ⋯
                  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.4228       9.8504      10.1217      10.5998       ⋯
                 b_e       0.4674       0.5646       0.6153       0.6372       ⋯
                 b_ω       0.9617       1.6069       1.7821       4.8210       ⋯
                 b_i       2.3053       2.4348       2.4910       2.5287       ⋯
                 b_Ω       0.2361       1.4130       1.7655       4.7276       ⋯
          ⋮                ⋮            ⋮            ⋮            ⋮            ⋱
                                                     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.