PyAstrOFit_Orbit.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PyAstrOFit: the orbit module"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We present here methods related to the **orbit** module, part of the PyAstrOFit package. Futhermore, this notebook constitutes a how-to-use illustrated with several examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import\n",
    "\n",
    "Whatever part of the notebook you plan to execute, you better first run this cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "\n",
    "import PyAstrOFit.Orbit as Orbit\n",
    "import numpy as np\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instanciate an orbit object"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will create a planet orbit, in this case $\\beta$ Pic b. To define a planet orbit, we need 6 parameters, for example the Kepler elements. One also need the mass of the star and its distance to the Earth: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 8.42        #[A.U.]\n",
    "e = 0.07\n",
    "i = 89.04       #[deg]\n",
    "node = 31.87    #[deg]\n",
    "w = 338.25      #[deg]\n",
    "tp = 2452353.32 #[JD]\n",
    "\n",
    "starMass = 1.75     #[solar mass]\n",
    "starDistance = 19.3 #[pc]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "betapicb = Orbit(semiMajorAxis = a,\n",
    "                   eccentricity = e,\n",
    "                   inclinaison = i,\n",
    "                   longitudeAscendingNode = node,\n",
    "                   periastron = w,\n",
    "                   periastronTime = tp,\n",
    "                   dStar = starDistance,\n",
    "                   starMass = starMass)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize the orbit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The method -showOrbit- allows to vizualize the full planet orbit. Several options are available. Let's discover them step by step."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simple representation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we can just display the full orbit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "betapicb.showOrbit(pointOfView = 'earth', \n",
    "                   unit = 'arcsec', \n",
    "                   figsize = (10,10), \n",
    "                   title = [r'$\\beta$ Pic b: a probable orbit',20])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Position at specific observation dates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can depict particular points on the orbit, corresponding to specific observation dates. This information can be retreived from the file 'betaPicb.txt' (included in the package resource folder) which contains the most recent observational data. For that, we use the method getInfo( ) from the module PyAstrOFit.FileHandler. We then display the planet orbit and zoom on a small region of interest. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import PyAstrOFit.FileHandler as FileHandler\n",
    "fileh = FileHandler('PyAstrOFit/res/exo/betaPicb.txt')\n",
    "observation_dates = fileh.getInfo()[0]['timePosition']\n",
    "\n",
    "betapicb.showOrbit(pointOfView = 'earth', \n",
    "                   unit = 'arcsec', \n",
    "                   figsize = (10,10), \n",
    "                   title = [r'$\\beta$ Pic b: a probable orbit',20],\n",
    "                   addPosition = observation_dates,\n",
    "                   lim = [[-0.30,-0.10],[-0.42,-0.22]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Adding observational data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To this figure, we can represent observational data with their associated error bars. This information can also be retreived from the file 'betaPicb.txt'. Furthermore, by setting the attribut -link- to True, we link the observed positions with the modeled positions at specific observation dates.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import PyAstrOFit.FileHandler as FileHandler\n",
    "fileh = FileHandler('PyAstrOFit/res/exo/betaPicb.txt')\n",
    "data_RA = fileh.getInfo()[0]['raPosition']\n",
    "err_RA = fileh.getInfo()[1]['raError']\n",
    "data_DEC = fileh.getInfo()[0]['decPosition']\n",
    "err_DEC = fileh.getInfo()[1]['decError']\n",
    "\n",
    "betapicb.showOrbit(pointOfView = 'earth', \n",
    "                   unit = 'arcsec', \n",
    "                   figsize = (10,10), \n",
    "                   title = [r'$\\beta$ Pic b: a probable orbit',20],\n",
    "                   addPosition = observation_dates,\n",
    "                   lim = [[-0.30,-0.10],[-0.42,-0.22]],\n",
    "                   addPoints = [data_RA,data_DEC,err_RA,err_DEC],\n",
    "                   link = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Saving the figure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, it's time to save our figure into a PDF file. For that, we modify the attribut save = (True,file_name). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name = 'betaPicb_test'\n",
    "betapicb.showOrbit(pointOfView = 'earth', \n",
    "                   unit = 'arcsec', \n",
    "                   figsize = (10,10), \n",
    "                   title = [r'$\\beta$ Pic b: a probable orbit',20],\n",
    "                   addPosition = observation_dates,\n",
    "                   lim = [[-0.30,-0.10],[-0.42,-0.22]],\n",
    "                   addPoints = [data_RA,data_DEC,err_RA,err_DEC],\n",
    "                   link = True,\n",
    "                   save = (True,file_name))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model Parameters and other quantities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can get the Kepler elements by using the method whichParameters( ). If it was not initially give, the period is now returns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = betapicb.whichParameters()\n",
    "for key in p:\n",
    "    print('{} = {}'.format(key,p[key]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The True anomaly $\\nu$ of an orbit point, given the time at position."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trueAnomalies = [betapicb.trueAnomaly(date) for date in observation_dates ]\n",
    "for k,nu in enumerate(trueAnomalies):\n",
    "    print('True anomaly corresponding to position {} = {}'.format(k,nu))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The time at position $t$, given the true anomaly $\\nu$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time_at_position = [betapicb.timeFromTrueAnomaly(nu,'iso') for nu in trueAnomalies]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determination of the $\\chi^2$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "chi2 = betapicb.chiSquare(obs = fileh.getInfo()[0],\n",
    "                          error = fileh.getInfo()[1])\n",
    "print('The non-reduced chi^2 equals {}'.format(chi2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Example: a more complete code for HR8799e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import PyAstrOFit.Orbit as Orbit\n",
    "import PyAstrOFit.FileHandler as FileHandler\n",
    "from PyAstrOFit.Planet_data import get_planet_data\n",
    "%matplotlib inline\n",
    "\n",
    "#path = 'PyAstrOFit/res/exo/hr8799c.txt'\n",
    "a = 3.80502142e+01        #[A.U.]\n",
    "e = 4.06273243e-02\n",
    "i = 2.50047414e+01       #[deg]\n",
    "node = 3.17132192e+02    #[deg]\n",
    "w = 8.15415302e+01      #[deg]\n",
    "tp = 2.46965054e+06 #[JD]\n",
    "\n",
    "starMass = 1.51     #[solar mass]\n",
    "starDistance = 39.4 #[pc]\n",
    "\n",
    "planet = Orbit(semiMajorAxis = a,\n",
    "                   eccentricity = e,\n",
    "                   inclinaison = i,\n",
    "                   longitudeAscendingNode = node,\n",
    "                   periastron = w,\n",
    "                   periastronTime = tp,\n",
    "                   dStar = starDistance,\n",
    "                   starMass = starMass)\n",
    "\n",
    "\n",
    "planet_data = get_planet_data('hr8799c')\n",
    "#fileh = FileHandler(path)\n",
    "observation_dates = planet_data[0]['timePosition']\n",
    "data_RA = planet_data[0]['raPosition']\n",
    "err_RA = planet_data[1]['raError']\n",
    "data_DEC = planet_data[0]['decPosition']\n",
    "err_DEC = planet_data[1]['decError']\n",
    "\n",
    "planet.showOrbit(pointOfView = 'earth', \n",
    "                   unit = 'arcsec', \n",
    "                   figsize = (10,10), \n",
    "                   title = ['Planet c',20],\n",
    "                   addPosition = observation_dates,\n",
    "                   lim = [[-0.9,-0.5],[0.45,0.85]],\n",
    "                   addPoints = [data_RA,data_DEC,err_RA,err_DEC],\n",
    "                   link = True,\n",
    "                   save = (False,'test_hr8799c'),\n",
    "                   color =(0,0.38,1.),\n",
    "                   link_color = 'r',\n",
    "                   addPosition_options={'marker':'o',\n",
    "                                        'markerfacecolor':(0,0.38,1.),\n",
    "                                        'markeredgecolor':(0,0.38,1.),\n",
    "                                        'markersize':8},\n",
    "                   invert_xaxis=True,\n",
    "                   invert_yaxis=False,\n",
    "                   cardinal=(False,'br'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}