CBP_stability_interpolation_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interpolation Functions for Circumbinary Planet Stability\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The go-to source for determining whether a circumbinary planet or CBP (i.e., planet that two stars completely) 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 discoveries made by the Kepler Space Telescope, these assumptions have not held up to reality.  It is true that we may expect CBPs 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.  Secondly, the parameter estimation of CBPs (e.g., planetary mass and eccentricity) **depend** on how the planet perturbs the binary orbit enough to affect the eclipse times.  \n",
    "\n",
    "Quarles et al. (2018) provides 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 on 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. (2018) using Kepler-16 as an example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Physical observables\n",
    "For eclipsing 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": [
    "# Kepler-16\n",
    "From Doyle et al. (2011) and Bender et al. (2012), we have measurements of the Kepler-16 binary for $q=0.2994 \\pm 0.0031$, $e_{bin} = 0.15894 \\pm 0.00079$, and $P_{bin} = 41.077580 \\pm 0.000008$ days.  Now to the preliminaries of the code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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,P_bin):\n",
    "    return f(mu,e_bin)[0]**(1.5)*P_bin"
   ]
  },
  {
   "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.genfromtxt(\"a_crit.txt\",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 (this ratio is being converted to period in our stability fuction)\n",
    "\n",
    "xi = np.linspace(0,0.5,51)\n",
    "yi = np.linspace(0,0.8,81)\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": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using Doyle, P_c = 181.8849 days\n",
      "Using Bender, P_c = 182.1355 days\n"
     ]
    }
   ],
   "source": [
    "mu_1 = 0.2937/(1 + 0.2937)  #converting from q --> mu  Doyle et al. (2011)\n",
    "mu_2 = 0.2994/(1 + 0.2994)  #converting from q --> mu  Bender et al. (2012)\n",
    "e_bin_1 = 0.15944\n",
    "e_bin_2 = 0.15894\n",
    "P_bin = 41.077580\n",
    "\n",
    "print(\"Using Doyle, P_c = %3.4f days\" % get_stability_limit(f,mu_1,e_bin_1,P_bin))\n",
    "print(\"Using Bender, P_c = %3.4f days\" % get_stability_limit(f,mu_2,e_bin_2,P_bin))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare the critical period $P_c$ to the measured orbital period of the CBP $P_b = 228.776$ days.  The ratio of the CBP orbital period to the critical period (using Bender) is $P_b/P_c = 1.25607$, which tells us that Kepler-16b is ~25% beyond the critical period, or ~15% beyond the critical semimajor axis ratio $a_c/a_{bin}$."
   ]
  },
  {
   "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
}