notebooks/2_exoplanet_example_orbit.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model the motion of a two-body system and plot its orbit\n",
    "\n",
    "In this notebook we create an OrbitSystem object and use it to plot parallax, proper, and orbital motion.\n",
    "\n",
    "Dependencies:\n",
    "    astropy, pystrometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import OrderedDict\n",
    "import glob\n",
    "import os\n",
    "import pprint \n",
    "\n",
    "import astropy.constants as constants\n",
    "from astropy.table import Table, vstack\n",
    "import astropy.units as u\n",
    "from IPython.display import display\n",
    "import matplotlib as mp\n",
    "import matplotlib.pylab as pl\n",
    "\n",
    "from pystrometry.pystrometry import OrbitSystem, get_cpsi_spsi_for_2Dastrometry "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Setting some defaults"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pp = pprint.PrettyPrinter(indent=4)     \n",
    "%matplotlib inline\n",
    "mp.rcParams['figure.figsize'] = (18, 9)\n",
    "mp.rcParams['font.size'] = 20\n",
    "mp.rcParams['ytick.direction'] = 'in'\n",
    "mp.rcParams['xtick.direction'] = 'in'\n",
    "\n",
    "save_plot = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Create temporary directory to save the figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LOCAL_PATH = os.path.dirname(os.getcwd())\n",
    "TEMPORARY_DIR = os.path.join(LOCAL_PATH, 'temp')\n",
    "if os.path.isdir(TEMPORARY_DIR) is False:\n",
    "    os.makedirs(TEMPORARY_DIR)\n",
    "print(LOCAL_PATH)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define the system parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# GJ 676 A b, see Sahlmann et al. 2016\n",
    "attribute_dict = OrderedDict([  ('RA_deg', 262.54667), \n",
    "                                  ('DE_deg', -51.6369),\n",
    "                                  ('absolute_plx_mas', 59.3), \n",
    "                                  ('Tp_day', 55410.4), \n",
    "                                  ('omega_deg', 87.4),\n",
    "                                  ('P_day', 1052.1), \n",
    "                                  ('ecc', 0.323), \n",
    "                                  ('OMEGA_deg', 208.),\n",
    "                                  ('i_deg', 45.), \n",
    "                                  ('m1_MS', 0.71),\n",
    "                                  ('m2_MJ', 6.7), \n",
    "                                  ('muRA_mas', -253.4), \n",
    "                                  ('muDE_mas', -177.9)])\n",
    "     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute astrometric motion components (parallax, proper motion, Keplerian orbit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# instantiate the pystrometry object                             \n",
    "orb = OrbitSystem(attribute_dict)\n",
    "\n",
    "# get default times for plotting\n",
    "t_mjd = orb.get_t_plot(format='mjd')\n",
    "\n",
    "# compute parallax and proper motion \n",
    "delta_ra_star_mas, delta_dec_mas = orb.ppm(t_mjd)\n",
    "\n",
    "# Generate coefficients to be able to model the data like one-dimensional measurements along the scan angle psi.\n",
    "# This allows us to use the same formalism for 1D (e.g. from Hipparcos) and 2D (e.g. from an imaging camera or catalog) astrometry\n",
    "# cpsi == cos(psi); spsi == sin(psi), see Sahlmann et al. 2011\n",
    "# xi and yi are the array indices corresponding to the measurements along perpendicular axes\n",
    "timestamps_curve_1D, cpsi_curve, spsi_curve, xi_curve, yi_curve = get_cpsi_spsi_for_2Dastrometry( t_mjd )\n",
    "\n",
    "# compute photocentre orbit (here == barycentre orbit)\n",
    "# phi0 == along-scan displacement\n",
    "phi0_curve = orb.photocenter_orbit(timestamps_curve_1D, spsi_curve, cpsi_curve)\n",
    "\n",
    "# compute relative orbit\n",
    "phi0_curve_relative = orb.relative_orbit_fast(timestamps_curve_1D, spsi_curve, cpsi_curve, shift_omega_by_pi = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot default orbit figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "orb.plot_orbits()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate the customised orbit figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_figure_columns = 3\n",
    "n_figure_rows = 1\n",
    "share_axes = False\n",
    "line_width = 3\n",
    "line_color = 'b'\n",
    "line_style = '-'\n",
    "orbit_colour = '0.5'\n",
    "\n",
    "fig, axes = pl.subplots(n_figure_rows, n_figure_columns,\n",
    "                        figsize=(n_figure_columns * 8, n_figure_rows * 6),\n",
    "                        facecolor='w', edgecolor='k', sharex=share_axes,\n",
    "                        sharey=share_axes)\n",
    "\n",
    "# plot parallax and proper motion\n",
    "axes[0].plot(delta_ra_star_mas, delta_dec_mas, 'k-', lw=line_width, color=line_color, ls=line_style, label='Parallax + Proper motion')    \n",
    "axes[0].set_ylabel('Offset in Declination (mas)')\n",
    "# plot photocentre orbit\n",
    "axes[0].plot(phi0_curve[xi_curve], phi0_curve[yi_curve],'k--', lw=line_width, color='k', ls=line_style, label='Photocentre orbit')\n",
    "# plot relative orbit\n",
    "axes[0].plot(phi0_curve_relative[xi_curve], phi0_curve_relative[yi_curve],'k--', lw=line_width, color=orbit_colour, ls='--', label='Relative orbit')\n",
    "\n",
    "# plot photocentre orbit in separate panel\n",
    "axes[1].plot(phi0_curve[xi_curve], phi0_curve[yi_curve],'k--', lw=line_width, color=orbit_colour, ls=line_style, label='Photocentre orbit')\n",
    "axes[1].plot(0,0,'kx', ms=15, mew=3, label='System barycentre')\n",
    "\n",
    "# plot relative orbit\n",
    "axes[2].plot(phi0_curve_relative[xi_curve], phi0_curve_relative[yi_curve],'k--', lw=line_width, color=orbit_colour, ls='--', label='Rel. orbit of companion')\n",
    "axes[2].plot(0,0,'k*', ms=15, mew=3, label='Primary')\n",
    "\n",
    "for ax in axes:\n",
    "    ax.legend(loc=3, fancybox=True, framealpha=0.)\n",
    "    ax.axis('equal')\n",
    "    ax.set_xlabel('Offset in Right Ascension (mas)')\n",
    "    ax.invert_xaxis()\n",
    "\n",
    "for ax in axes[1:]:    \n",
    "    ax.axhline(y=0,color='0.7',ls='--',zorder=-50)\n",
    "    ax.axvline(x=0,color='0.7',ls='--',zorder=-50)\n",
    "    \n",
    "pl.show()    \n",
    "if save_plot:\n",
    "    figure_file_name = os.path.join(TEMPORARY_DIR, 'example_orbit.png')\n",
    "    fig.savefig(figure_file_name, transparent=True, bbox_inches='tight', pad_inches=0.05)        "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "astrom",
   "language": "python",
   "name": "astrom"
  },
  "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.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}