Merge pull request #356 from sblunt/fit-arbitrary-abs-astrom

Fit arbitrary abs astrom

docs/index.rst

@@ -55,13 +55,15 @@ User Guide:
 Changelog:
 ++++++++++
 
-**3.0.0 (TBD)**
+**3.0.0 (2024-4-11)**
 
+- implementation of Hipparcos-Gaia catalog of accelerations fitting! (@semaphoreP)
+- fit arbitrary absolute astrometry (@sblunt)
+- implement O'Neil observation-based priors (@sblunt/@clarissardoo)
 - discuss MCMC autocorrelation in MCMC tutorial (@michaelkmpoon)
 - add time warning if OFTI doesn't accept an orbit in first 60 s (@michaelkmpoon)
 - add first parts of orbitize! manual (@sofiacovarrubias/@sblunt)
 - bugfix for rebound MCMC fits (issue #357; @sblunt)
-- implementation of Hipparcos-Gaia catalog of accelerations fitting! (@semaphoreP)
 - implementation of residual plotting method for orbit plots (@Saanikachoudhary and @semaphoreP)
 - plot companion RVs (@chihchunhsu)
 - add documentation about referencing issues when modifying priors to tutorial (@wcroberson)

docs/tutorials.rst

@@ -56,6 +56,7 @@ us if you are still confused).
    tutorials/Changing_bases_tutorial.ipynb
    tutorials/Hipparcos_IAD.ipynb
    tutorials/HGCA_tutorial.ipynb
+   tutorials/abs_astrometry.ipynb
 
 
 

docs/tutorials/Hipparcos_IAD.ipynb

@@ -49,17 +49,14 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Go get a coffee. This will take a few mins! :)\n",
-      "Starting burn-in!\n",
-      "Starting production chain!\n",
-      "Done! This fit took 3.1 mins on my machine.\n"
+      "Go get a coffee. This will take a few mins! :)\n"
      ]
     },
     {
@@ -4281,7 +4278,6 @@
     "#   although I'd probably run it for 5,000-10,000 steps if I wanted to publish it.\n",
     "burn_steps = 100\n",
     "mcmc_steps = 100\n",
-    "\n",
     "start = datetime.now()\n",
     "\n",
     "# run the fit\n",

docs/tutorials/abs_astrometry.ipynb

@@ -0,0 +1,275 @@
+{
+ "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",
+    "deoffs = 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, deoff = 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",
+    "    deoffs[i] = deoff.to(u.mas).value\n",
+    "\n",
+    "print(raoffs, deoffs)"
+   ]
+  },
+  {
+   "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": [
+       "<div>\n",
+       "<style scoped>\n",
+       "    .dataframe tbody tr th:only-of-type {\n",
+       "        vertical-align: middle;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe tbody tr th {\n",
+       "        vertical-align: top;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr style=\"text-align: right;\">\n",
+       "      <th></th>\n",
+       "      <th>epoch</th>\n",
+       "      <th>object</th>\n",
+       "      <th>raoff</th>\n",
+       "      <th>deoff</th>\n",
+       "      <th>deoff_err</th>\n",
+       "      <th>raoff_err</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>0</th>\n",
+       "      <td>58894.1644</td>\n",
+       "      <td>0</td>\n",
+       "      <td>377.813056</td>\n",
+       "      <td>1044.819572</td>\n",
+       "      <td>123.4</td>\n",
+       "      <td>123.4</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1</th>\n",
+       "      <td>58934.8270</td>\n",
+       "      <td>0</td>\n",
+       "      <td>1425.176093</td>\n",
+       "      <td>933.702905</td>\n",
+       "      <td>123.4</td>\n",
+       "      <td>123.4</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "        epoch  object        raoff        deoff  deoff_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[\"deoff\"] = deoffs\n",
+    "\n",
+    "df_orbitize[\"deoff_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
+}