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 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=: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
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"
└ @ DataDeps ~/.julia/packages/DataDeps/Y2lje/src/verification.jl:44
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)
Extracting archive: /home/runner/.julia/scratchspaces/124859b0-ceae-595e-8997-d05f6a7a8dfe/datadeps/Whereistheplanet/whereistheplanet-master.zip
--
Path = /home/runner/.julia/scratchspaces/124859b0-ceae-595e-8997-d05f6a7a8dfe/datadeps/Whereistheplanet/whereistheplanet-master.zip
Type = zip
Physical Size = 11870993
Comment = d393657c7b6224ed832209e1218a4cd00c6ba0ec
Everything is Ok
Folders: 8
Files: 99
Size: 15005637
Compressed: 11870993We 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)
chainChains 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.