Fitting Arbitrary Absolute Astrometry

by Sarah Blunt (2023)

This tutorial walks you through using orbitize! to perform a fit on arbitary absolute astrometry. By “arbitrary,” I mean astrometry not taken by Gaia or Hipparcos (which orbitize! has dedicated modules for; see the HGCA and Hipparcos IAD tutorials). Let’s imagine we have astrometry for a single star derived from wide-field images taken over several years, and we want to combine these data with measurements from Hipparcos. We are going to perform a fit to jointly constrain astrometric parameters (parallax and proper motion) and orbital parameters of a secondary companion.

This tutorial will take you through: - formatting absolute astrometry measurements for input into orbitize! - setting up an orbit fit incorporating these measurements

This tutorial assumes the following prerequities: - Using the Hipparcos IAD

Input Data Format

Following Nielsen et al 2020 (see the Hipparcos IAD tutorial), orbitize! defines astrometric data points as offset from the reported Hipparcos position at the reported Hipparcos epoch. Let’s start by defining an orbitize.hipparcos.Hipparcos object, which holds onto information from the Hipparcos mission observations of our object of interest. I’m going to use beta Pictoris as an example since you already have that IAD file in your orbitize! distribution. See the IAD tutorial for info on how to download the data for your object.

from orbitize import hipparcos, DATADIR

hip_num = "027321"  # beta Pic

# Location of the Hipparcos IAD file.
IAD_file = "{}H{}.d".format(DATADIR, hip_num)

# The HipparcosLogProb object needs to know how many companions are in your fit
# in order to compute likelihood. There are 2 known planets around beta Pic, but let's
# keep it simple for the tutorial
num_secondary_bodies = 1

betaPicHipObject = hipparcos.HipparcosLogProb(IAD_file, hip_num, num_secondary_bodies)

Generally, when you’re deriving (or using published) absolute astrometry, it will be in the form 82 02 14.35787 (J2000). However, orbitize! expects astrometry to be input relative to the Hipparcos position. Our friends at astropy have made these calculations very easy to do! Here’s an example:

from astropy.coordinates import SkyCoord
from astropy import units as u
import numpy as np

# let's imagine our data look like this:
datapoints = ["05 47 17.123456 -51 03 59.123456", "05 47 17.234567 -51 03 59.234567"]
data_epochs = ["2020.1234", "2020.2345"]
num_datapoints = len(datapoints)

hipparcos_coordinate = SkyCoord(
    betaPicHipObject.alpha0, betaPicHipObject.delta0, unit=(u.deg, u.deg)
)

raoffs = np.zeros(num_datapoints)
deoffs = np.zeros(num_datapoints)
for i in range(num_datapoints):
    my_data_coordinate = SkyCoord(datapoints[i], unit=(u.hourangle, u.deg))

    # take difference between reported Hipparcos position and convert to mas
    raoff, deoff = hipparcos_coordinate.spherical_offsets_to(my_data_coordinate)

    # n.b. orbitize! expects raw ra offsets, NOT multiplied by cos(delta0). Don't
    # multiply by cos(delta0) here.
    raoffs[i] = raoff.to(u.mas).value
    deoffs[i] = deoff.to(u.mas).value

print(raoffs, deoffs)

[ 377.81305615 1425.17609269] [1044.81957171  933.70290544]

Sweet! These absolute astrometry points are now suitable for an orbitize! input file. You can add them to an existing file with other types of data (relative astrometry and RVs) and/or fit them on their own. Here’s what the data file for our two points would look like:

from pandas import DataFrame
from astropy.time import Time

df_orbitize = DataFrame(Time(data_epochs, format="decimalyear").mjd, columns=["epoch"])

# this line tells orbitize! "these measurements are astrometry of the primary"
df_orbitize["object"] = 0

df_orbitize["raoff"] = raoffs
df_orbitize["deoff"] = deoffs

df_orbitize["deoff_err"] = 123.4  # error on the declination measurement, in mas
df_orbitize["raoff_err"] = 123.4  # error on the RA measurement, in mas

df_orbitize.to_csv("data_for_orbit_fit.csv", index=False)
df_orbitize

	epoch	object	raoff	deoff	deoff_err	raoff_err
0	58894.1644	0	377.813056	1044.819572	123.4	123.4
1	58934.8270	0	1425.176093	933.702905	123.4	123.4

Setting up & Running Your Fit

The hard part is over– we have formatted our input data! orbitize! will now function the same as any other fit. Behind the scenes, orbitize! will automatically recognize that you have inputted absolute astrometry, and set up a fit that includes position, parallax, and proper motion terms as free parameters. Observe:

from orbitize import read_input, system, priors, sampler
import os

data_table = read_input.read_file("data_for_orbit_fit.csv")

fit_secondary_mass = True  # tell orbitize! we want to get dynamical masses
m0 = 1
plx = 1

# this sets up a joint fit of Hipparcos time series data and the absolute astrometry
# from the data table we just created.
betaPicSystem = system.System(
    num_secondary_bodies,
    data_table,
    m0,
    plx,
    hipparcos_IAD=betaPicHipObject,
    fit_secondary_mass=fit_secondary_mass,
)

# change any priors you want to:
plx_idx = betaPicSystem.param_idx["plx"]
betaPicSystem.sys_priors[plx_idx] = priors.UniformPrior(10, 15)

# run the fit!
tutorialSampler = sampler.MCMC(betaPicSystem)
# tutorialSampler.run_sampler(you_choose, burn_steps=you_choose)

# clean up
os.system("rm data_for_orbit_fit.csv")

0