CSP_stability_interpolation_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interpolation Functions for Circumstellar Planet Stability\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The go-to source for determining whether a circumstellar planet (i.e., planet that orbits only one of two stars in a binary) has been an empirically determined fitting formula by Holman & Wiegert (1999).  In their study, the key assumptions are that the planet begins on a circular, coplanar orbit and does not have sufficient mass to influence the orbits of the massive binary stars (i.e., massless test particle).  Through the many exoplanet discoveries made over the last few decades, these assumptions have not held up to reality.  It is true that we may expect circumstellar planets to begin with on nearly circular and coplanar orbits relative to the host binary, but disk interactions or planet-planet interactions over the lifetime of the system can perturb the discovered planets away from our ideal assumption.  \n",
    "\n",
    "Quarles et al. (2020) provide a somewhat different approach.  Instead of using a fitting function based on some physical parameters (e.g., stellar dynamical mass ratio $\\mu$ and binary eccentricity $e_{bin}$), they use an interpolation scheme that is built upon a much larger set of N-body simulations.  In these simulations, all three bodies interact with one another gravitationally and more accurately represent what we expect to observe.  This Jupyter notebook will show you how to utilize the results of Quarles et al. (2020) using $\\alpha$ Centauri as an example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Physical observables\n",
    "For wide binary systems, one can usually obtain an estimate of the mass quotient $q=M_B/M_A$, binary orbital period $P_{bin}$, and the binary eccentricity $e_{bin}$ from modeling the observed photometric data (i.e., light curve) and radial velocity measurements (RV curve) together.  More intense modeling can help determine the individual masses of the stars, but we can convert between the observational mass quotient $q$ to the dynamical mass ratio $\\mu$ through the following:\n",
    "\n",
    "$\\mu = \\frac{M_B}{M_A+M_B} = \\frac{q}{1+q}$\n",
    "\n",
    "or\n",
    "\n",
    "$q = \\frac{\\mu}{1-\\mu}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Alpha Centauri \n",
    "From Pourbaix & Boffin (2016), we have measurements of the Alpha Centauri AB binary for $M_A=1.133 \\pm 0.005$, $M_B=0.972 \\pm 0.005$, $e_{bin} = 0.524 \\pm 0.001$, and $P_{bin} = 79.91 \\pm 0.013$ years.  Now to the preliminaries of the code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np  #We need numpy for the convenience function genfromtxt\n",
    "from scipy.interpolate import griddata,interp2d  #We need this to map our data to a grid and interpolate between grid points\n",
    "\n",
    "def get_stability_limit(f,mu,e_bin):\n",
    "    return f(mu,e_bin)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The import statements above include the standard packages from numpy and scipy.  The function **get_stability_limit** is something for convenience so that we may generalize later and is passed a set of parameters.  The interpolation function $f$ that we create using scipy and the binary parameters.  Next, let's create the interpolation function for a coplanar planet ($i_p = 0^\\circ$) and determine the stability limit for a planet orbiting star A."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "ip = 0\n",
    "data = np.genfromtxt(\"a_crit_Incl[%i].txt\" % ip,delimiter=',',comments='#')  #The data contained in this repository\n",
    "\n",
    "X = data[:,0] #mu\n",
    "Y = data[:,1] #e_bin\n",
    "Z = data[:,2] #a_c/a_bin\n",
    "\n",
    "xi = np.concatenate(([0.001],np.arange(0.01,1,0.01),[0.999]))\n",
    "yi = np.arange(0,0.81,0.01)\n",
    "zi = griddata((X,Y),Z,(xi[None,:],yi[:,None]),method = 'linear',fill_value=0)  #make the grid\n",
    "\n",
    "f = interp2d(xi, yi, zi, kind='linear') # make the 2d interpolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a_c = 2.773 AU\n"
     ]
    }
   ],
   "source": [
    "M_tot = 0.972 + 1.133 # Total star mass in M_sun\n",
    "mu = 0.972/(0.972 + 1.133)  #converting from M_A, M_B --> mu  Pourbaix & Boffin (2016)\n",
    "e_bin = 0.524\n",
    "P_bin = 79.91\n",
    "a_bin = (P_bin**2*M_tot)**(1./3)\n",
    "\n",
    "print(\"a_c = %1.3f AU\" % (get_stability_limit(f,mu,e_bin)*a_bin))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare the critical semimajor axis $a_c$ to the long-term (10 Gyr) simulations from Quarles, Lissauer, and Kaib (2018).  They determine the stability limit for a initially ciruclar planet orbiting star A at around 2.7 AU, where this can be extended to 2.77 AU if the planet begins with an appropriate forced eccentricity."
   ]
  },
  {
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}