# 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_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 with fields .ℓπcallback and .∇ℓπcallback


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

# Typical depth is ~4, so limiting the max_depth down from default 12 speeds warm-up
chain = octofit(model, iterations=4_000, max_depth=6)
Chains MCMC chain (4000×30×1 Array{Float64, 3}):

Iterations        = 1:1:4000
Number of chains  = 1
Samples per chain = 4000
Wall duration     = 5.87 seconds
Compute duration  = 5.87 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         = 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  ⋯
Symbol      Float64   Float64   Float64     Float64     Float64  ⋯

plx      33.9872    0.3383    0.0048   4885.6451   3056.9226  ⋯
pmra      44.2135    0.3368    0.0050   4482.2881   3276.2329  ⋯
pmdec     -64.3921    0.2739    0.0041   4374.9667   3311.7790  ⋯
ra_hip_offset_mas      -0.0009    0.2016    0.0032   3960.6553   3091.6610  ⋯
dec_hip_offset_mas      -0.0053    0.2909    0.0043   4584.8065   2948.7832  ⋯
M       1.0000    0.0000       NaN         NaN         NaN  ⋯
rv       0.0000    0.0000       NaN         NaN         NaN  ⋯
dec      -2.4734    0.0000    0.0000   3960.6443   3091.6610  ⋯
ra      69.4004    0.0000    0.0000   4584.8065   2948.7832  ⋯
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      33.2964      33.7609      33.9909      34.2121       ⋯
pmra      43.5613      43.9853      44.2145      44.4371       ⋯
pmdec     -64.9221     -64.5771     -64.3945     -64.2052     - ⋯
ra_hip_offset_mas      -0.4043      -0.1350      -0.0000       0.1353       ⋯
dec_hip_offset_mas      -0.5870      -0.2030      -0.0043       0.1923       ⋯
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


Plot the posterior values:

octoplot(model,chain,show_astrom=false,show_astrom_time=false)

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=: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

## Constrain Planet Mass

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

astrom_like1,astrom_like2 = Octofitter.Whereistheplanet_astrom("51erib",object=1)
┌ Warning: Checksum not provided, add to the Datadep Registration the following hash line
│   hash = "344043fa845f1a096f0a67f32bce9eeb8c94b1ae210a6b91a7fddfe0602c30ec"

7-Zip (a) [64] 17.04 : Copyright (c) 1999-2021 Igor Pavlov : 2017-08-28
p7zip Version 17.04 (locale=C.UTF-8,Utf16=on,HugeFiles=on,64 bits,4 CPUs AMD EPYC 7763 64-Core Processor                 (A00F11),ASM,AES-NI)

Scanning the drive for archives:
1 file, 11870993 bytes (12 MiB)

--
Type = zip
Physical Size = 11870993
Comment = d393657c7b6224ed832209e1218a4cd00c6ba0ec

Everything is Ok

Folders: 8
Files: 99
Size:       15005637
Compressed: 11870993

We specify our full model:

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

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

rv =  12.60e3 # m/s
plx ~ Uniform(0,100)
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, 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/cos(system.dec)

ref_epoch = Octofitter.hipparcos_catalog_epoch_mjd
end hip_like b

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


Now we sample:

chain = octofit(model)
chain
Chains MCMC chain (1000×39×1 Array{Float64, 3}):

Iterations        = 1:1:1000
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 14.98 seconds
Compute duration  = 14.98 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_ωy, b_ωx, b_i, b_Ωy, b_Ωx, b_θy, b_θx, b_ω, b_Ω, b_θ, b_tp, b_mass
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    ⋯
Symbol      Float64   Float64   Float64    Float64    Float64    ⋯

M_pri       0.0300    0.0000    0.0000     2.7478        NaN    ⋯
M_sec     100.0000    0.0000    0.0000     2.4460    12.2549    ⋯
plx      99.8193    0.0063    0.0034     2.3835    12.2549    ⋯
pmra     -76.8476    0.8557    0.4492     2.3837    12.2549    ⋯
pmdec     -86.4913    0.6983    0.3647     2.3837    12.2549    ⋯
ra_hip_offset_mas    -203.9307    0.0011    0.0006     2.3975    12.2549    ⋯
dec_hip_offset_mas    4871.0452    0.0053    0.0028     2.3836    12.0444    ⋯
M       0.1255    0.0000    0.0000     2.4430    11.9953    ⋯
rv   12600.0000    0.0000       NaN        NaN        NaN    ⋯
dec      -2.4734    0.0000    0.0000     2.3975    12.2549    ⋯
ra      69.3987    0.0000    0.0000     2.3836    12.2549    ⋯
ref_epoch   48348.5625    0.0000       NaN        NaN        NaN    ⋯
b_a     477.8678   11.1742    5.3018     2.4035    12.2555    ⋯
b_e       0.9851    0.0005    0.0003     2.3835    12.2549    ⋯
b_ωy       2.3556    0.0298    0.0155     2.3845    12.2549    ⋯
b_ωx      -0.9565    0.0185    0.0093     2.4412    12.2549    ⋯
b_i       1.5198    0.0040    0.0021     2.4453    12.2549    ⋯
⋮                ⋮           ⋮         ⋮         ⋮          ⋮        ⋱
2 columns and 9 rows omitted

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

M_pri       0.0300       0.0300       0.0300       0.0300       ⋯
M_sec     100.0000     100.0000     100.0000     100.0000     1 ⋯
plx      99.8032      99.8208      99.8222      99.8228       ⋯
pmra     -77.3292     -77.2923     -77.2220     -77.0826     - ⋯
pmdec     -86.8785     -86.8511     -86.7953     -86.6866     - ⋯
ra_hip_offset_mas    -203.9339    -203.9304    -203.9302    -203.9301    -2 ⋯
dec_hip_offset_mas    4871.0421    4871.0424    4871.0429    4871.0438    48 ⋯
M       0.1255       0.1255       0.1255       0.1255       ⋯
rv   12600.0000   12600.0000   12600.0000   12600.0000   126 ⋯
dec      -2.4734      -2.4734      -2.4734      -2.4734       ⋯
ra      69.3987      69.3987      69.3987      69.3987       ⋯
ref_epoch   48348.5625   48348.5625   48348.5625   48348.5625   483 ⋯
b_a     472.8102     473.0719     473.2738     473.5860     5 ⋯
b_e       0.9836       0.9853       0.9854       0.9854       ⋯
b_ωy       2.2701       2.3639       2.3684       2.3709       ⋯
b_ωx      -1.0115      -0.9509      -0.9488      -0.9475       ⋯
b_i       1.5082       1.5208       1.5215       1.5219       ⋯
⋮                ⋮            ⋮            ⋮            ⋮            ⋱
1 column and 9 rows omitted

octoplot(model, chain, show_mass=true)

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