{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Fitting Arbitrary Absolute Astrometry\n", "\n", "by Sarah Blunt (2023)\n", "\n", "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](https://orbitize.readthedocs.io/en/latest/tutorials/Hipparcos_IAD.html)). 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. \n", "\n", "This tutorial will take you through:\n", "- formatting absolute astrometry measurements for input into orbitize!\n", "- setting up an orbit fit incorporating these measurements\n", "\n", "This tutorial assumes the following prerequities:\n", "- [Using the Hipparcos IAD](https://orbitize.readthedocs.io/en/latest/tutorials/Hipparcos_IAD.html)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Input Data Format\n", "\n", "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](https://orbitize.readthedocs.io/en/latest/tutorials/Hipparcos_IAD.html) for info on how to download the data for your object." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from orbitize import hipparcos, DATADIR\n", "\n", "hip_num = \"027321\" # beta Pic\n", "\n", "# Location of the Hipparcos IAD file.\n", "IAD_file = \"{}H{}.d\".format(DATADIR, hip_num)\n", "\n", "# The HipparcosLogProb object needs to know how many companions are in your fit\n", "# in order to compute likelihood. There are 2 known planets around beta Pic, but let's\n", "# keep it simple for the tutorial\n", "num_secondary_bodies = 1\n", "\n", "betaPicHipObject = hipparcos.HipparcosLogProb(IAD_file, hip_num, num_secondary_bodies)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 377.81305615 1425.17609269] [1044.81957171 933.70290544]\n" ] } ], "source": [ "from astropy.coordinates import SkyCoord\n", "from astropy import units as u\n", "import numpy as np\n", "\n", "# let's imagine our data look like this:\n", "datapoints = [\"05 47 17.123456 -51 03 59.123456\", \"05 47 17.234567 -51 03 59.234567\"]\n", "data_epochs = [\"2020.1234\", \"2020.2345\"]\n", "num_datapoints = len(datapoints)\n", "\n", "hipparcos_coordinate = SkyCoord(\n", " betaPicHipObject.alpha0, betaPicHipObject.delta0, unit=(u.deg, u.deg)\n", ")\n", "\n", "raoffs = np.zeros(num_datapoints)\n", "decoffs = np.zeros(num_datapoints)\n", "for i in range(num_datapoints):\n", " my_data_coordinate = SkyCoord(datapoints[i], unit=(u.hourangle, u.deg))\n", "\n", " # take difference between reported Hipparcos position and convert to mas\n", " raoff, decoff = hipparcos_coordinate.spherical_offsets_to(my_data_coordinate)\n", "\n", " # n.b. orbitize! expects raw ra offsets, NOT multiplied by cos(delta0). Don't\n", " # multiply by cos(delta0) here.\n", " raoffs[i] = raoff.to(u.mas).value\n", " decoffs[i] = decoff.to(u.mas).value\n", "\n", "print(raoffs, decoffs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
epochobjectraoffdecoffdecoff_errraoff_err
058894.16440377.8130561044.819572123.4123.4
158934.827001425.176093933.702905123.4123.4
\n", "
" ], "text/plain": [ " epoch object raoff decoff decoff_err raoff_err\n", "0 58894.1644 0 377.813056 1044.819572 123.4 123.4\n", "1 58934.8270 0 1425.176093 933.702905 123.4 123.4" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from pandas import DataFrame\n", "from astropy.time import Time\n", "\n", "df_orbitize = DataFrame(Time(data_epochs, format=\"decimalyear\").mjd, columns=[\"epoch\"])\n", "\n", "# this line tells orbitize! \"these measurements are astrometry of the primary\"\n", "df_orbitize[\"object\"] = 0\n", "\n", "df_orbitize[\"raoff\"] = raoffs\n", "df_orbitize[\"decoff\"] = decoffs\n", "\n", "df_orbitize[\"decoff_err\"] = 123.4 # error on the declination measurement, in mas\n", "df_orbitize[\"raoff_err\"] = 123.4 # error on the RA measurement, in mas\n", "\n", "df_orbitize.to_csv(\"data_for_orbit_fit.csv\", index=False)\n", "df_orbitize" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Setting up & Running Your Fit\n", "\n", "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:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from orbitize import read_input, system, priors, sampler\n", "import os\n", "\n", "data_table = read_input.read_file(\"data_for_orbit_fit.csv\")\n", "\n", "fit_secondary_mass = True # tell orbitize! we want to get dynamical masses\n", "m0 = 1\n", "plx = 1\n", "\n", "# this sets up a joint fit of Hipparcos time series data and the absolute astrometry\n", "# from the data table we just created.\n", "betaPicSystem = system.System(\n", " num_secondary_bodies,\n", " data_table,\n", " m0,\n", " plx,\n", " hipparcos_IAD=betaPicHipObject,\n", " fit_secondary_mass=fit_secondary_mass,\n", ")\n", "\n", "# change any priors you want to:\n", "plx_idx = betaPicSystem.param_idx[\"plx\"]\n", "betaPicSystem.sys_priors[plx_idx] = priors.UniformPrior(10, 15)\n", "\n", "# run the fit!\n", "tutorialSampler = sampler.MCMC(betaPicSystem)\n", "# tutorialSampler.run_sampler(you_choose, burn_steps=you_choose)\n", "\n", "# clean up\n", "os.system(\"rm data_for_orbit_fit.csv\")" ] } ], "metadata": { "kernelspec": { "display_name": "python3.12", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.2" } }, "nbformat": 4, "nbformat_minor": 2 }