CompanionDilutionFactors.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as colors\n",
    "from matplotlib.patches import Circle\n",
    "from matplotlib.collections import PatchCollection\n",
    "\n",
    "from astropy.io import fits\n",
    "import os, sys\n",
    "\n",
    "import photutils\n",
    "\n",
    "import astropy.units as u\n",
    "\n",
    "from astropy.wcs import WCS\n",
    "from astropy.wcs.utils import skycoord_to_pixel\n",
    "from astropy.coordinates import SkyCoord\n",
    "\n",
    "import urllib\n",
    "\n",
    "from SPCA import Photometry_Aperture"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch = 2\n",
    "oversamp = True\n",
    "oversampTimes = 10\n",
    "method = 'exact'\n",
    "\n",
    "# # WASP-103b\n",
    "# aperad = 2.6\n",
    "# dx = 0.13218\n",
    "# dy = 0.15009\n",
    "# tinytim_path = '/home/taylor/Downloads/stinytimv2.0/stinytim/WASP103.fits'\n",
    "\n",
    "# WASP-12b\n",
    "aperad = 2.2 # Old\n",
    "# aperad = 3.2 # New\n",
    "dx = 0.31836\n",
    "dy = 0.831538\n",
    "tinytim_path = '/home/taylor/Downloads/stinytimv2.0/stinytim/WASP12.fits'\n",
    "\n",
    "# # CoRoT-2b\n",
    "# aperad = 3.6\n",
    "# dx = 1.625\n",
    "# dy = 2.9933\n",
    "# tinytim_path = '/home/taylor/Downloads/stinytimv2.0/stinytim/CoRoT2.fits'\n",
    "\n",
    "if oversamp:\n",
    "    aperad *= oversampTimes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 201,
   "metadata": {},
   "outputs": [],
   "source": [
    "file = fits.open(tinytim_path)\n",
    "data = file[0].data\n",
    "data = data[np.newaxis,:]\n",
    "header = file[0].header\n",
    "file.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 202,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "187.44419918062755 187.64919235360014\n",
      "187.12583918062757 186.81765435360015\n"
     ]
    }
   ],
   "source": [
    "if oversamp:\n",
    "    dx *= oversampTimes\n",
    "    dy *= oversampTimes\n",
    "\n",
    "cx_empirical = np.mean(np.arange(data[0].shape[0])*np.sum(data[0], axis=1)/np.mean(np.sum(data[0], axis=1)))\n",
    "cy_empirical = np.mean(np.arange(data[0].shape[1])*np.sum(data[0], axis=0)/np.mean(np.sum(data[0], axis=0)))\n",
    "\n",
    "cx = (cx_empirical-dx)\n",
    "cy = (cy_empirical-dy)\n",
    "\n",
    "print(cx_empirical/oversampTimes, cy_empirical/oversampTimes)\n",
    "print(cx/oversampTimes, cy/oversampTimes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 203,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeEAAAHWCAYAAABANTEbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3da4xk533f+e+/rl3dPT0XzpAazkwsRqYSk9osbREEF1oE3sAbMXqxlBYrgH5hEVkhdAQJsAFnASkBNvK7bBDbgLCQABoWRCVeEwRsQ4xh7YYmbBiLZUSTDiXeRGtsKuKIIw7JmWZf6l713xd1ptUe9dwvT3fP9wMU6tRT55x6Hh70/Hie85znRGYiSZJuvFrpCkiSdLMyhCVJKsQQliSpEENYkqRCDGFJkgoxhCVJKuSGh3BEPBARr0XE8Yj4/I3+fUmStou4kfcJR0Qd+CvgfwROAH8B/GJmvnLDKiFJ0jZxo8+E7wOOZ+bfZOYQeBx48AbXQZKkbeFGh/AR4I1Nn09UZZIk3XQaN/j3Youyn+gPj4hHgEeqjx++rjWSJOn6eyczD51beKND+ARwbNPno8Cb566UmY8CjwJEhJNbS5J2uv+6VeGN7o7+C+DOiLgjIlrAQ8CTN7gOkiRtCzf0TDgzxxHxOeD/AerAVzPz5RtZB0mStosbeovSlbA7WpK0CzyfmfeeW+iMWZIkFWIIS5JUiCEsSVIhhrAkSYUYwpIkFWIIS5JUiCEsSVIhhrAkSYUYwpIkFWIIS5JUiCEsSVIhhrAkSYUYwpIkFWIIS5JUiCEsSVIhhrAkSYUYwpIkFWIIS5JUiCEsSVIhhrAkSYUYwpIkFWIIS5JUiCEsSVIhhrAkSYUYwpIkFWIIS5JUiCEsSVIhhrAkSYUYwpIkFWIIS5JUiCEsSVIhhrAkSYUYwpIkFWIIS5JUiCEsSVIhhrAkSYUYwpIkFWIIS5JUiCEsSVIhhrAkSYUYwpIkFWIIS5JUiCEsSVIhhrAkSYUYwpIkFWIIS5JUiCEsSVIhhrAkSYUYwpIkFXJVIRwR34+IFyPihYh4rio7EBFPRcT3qvf9m9b/QkQcj4jXIuKjV1t5SZJ2smtxJvw/ZOY9mXlv9fnzwNOZeSfwdPWZiLgLeAi4G3gA+HJE1K/B70uStCNdj+7oB4HHquXHgI9vKn88MweZ+TpwHLjvOvy+JEk7wtWGcAL/KSKej4hHqrLbMvMkQPV+a1V+BHhj07YnqjJJkm5Kjavc/iOZ+WZE3Ao8FRHfvcC6sUVZbrniLNAf2eo7SZJ2i6s6E87MN6v3U8AfMutefisiDgNU76eq1U8AxzZtfhR48zz7fTQz7910nVmSpF3nikM4IhYiYs/ZZeAfAy8BTwIPV6s9DHyjWn4SeCgi2hFxB3An8OyV/r4kSTvd1XRH3wb8YUSc3c//lZn/d0T8BfBERHwa+AHwSYDMfDkingBeAcbAZzNzclW1lyRpB4vMLS/LbhsRsb0rKEnSxT2/1SVWZ8ySJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKmQi4ZwRHw1Ik5FxEubyg5ExFMR8b3qff+m774QEccj4rWI+Oim8g9HxIvVd1+KiLj2zZEkaee4lDPhrwEPnFP2eeDpzLwTeLr6TETcBTwE3F1t8+WIqFfbfAV4BLizep27T0mSbioXDeHM/HPg9DnFDwKPVcuPAR/fVP54Zg4y83XgOHBfRBwGljLzmcxM4OubtpEk6aZ0pdeEb8vMkwDV+61V+RHgjU3rnajKjlTL55ZLknTTalzj/W11nTcvUL71TiIeYdZ1LUnSrnWlZ8JvVV3MVO+nqvITwLFN6x0F3qzKj25RvqXMfDQz783Me6+wfpIkbXtXGsJPAg9Xyw8D39hU/lBEtCPiDmYDsJ6tuqxXI+L+alT0pzZtI0nSTemi3dER8XvAzwMHI+IE8K+BfwM8ERGfBn4AfBIgM1+OiCeAV4Ax8NnMnFS7+gyzkdYd4JvVS5Kkm1bMBitvXxGxvSsoSdLFPb/VJVZnzJIkqRBDWJKkQgxhSZIKMYQlSSrEEJYkqRBDWJKkQgxhSZIKMYQlSSrEEJYkqRBDWJKkQgxhSZIKMYQlSSrEEJYkqRBDWJKkQgxhSZIKMYQlSSrEEJYkqRBDWJKkQgxhSZIKMYQlSSrEEJYkqRBDWJKkQgxhSZIKMYQlSSrEEJYkqRBDWJKkQgxhSZIKMYQlSSrEEJYkqRBDWJKkQgxhSZIKMYQlSSrEEJYkqRBDWJKkQgxhSZIKMYQlSSrEEJYkqRBDWJKkQgxhSZIKMYQlSSrEEJYkqRBDWJKkQgxhSZIKMYQlSSrEEJYkqRBDWJKkQgxhSZIKMYQlSSrEEJYkqZCLhnBEfDUiTkXES5vKvhgRP4yIF6rXxzZ994WIOB4Rr0XERzeVfzgiXqy++1JExLVvjiRJO8elnAl/DXhgi/Lfysx7qtcfA0TEXcBDwN3VNl+OiHq1/leAR4A7q9dW+5Qk6aZx0RDOzD8HTl/i/h4EHs/MQWa+DhwH7ouIw8BSZj6TmQl8Hfj4lVZakqTd4GquCX8uIr5TdVfvr8qOAG9sWudEVXakWj63XJKkm9aVhvBXgA8A9wAngd+oyre6zpsXKN9SRDwSEc9FxHNXWD9Jkra9KwrhzHwrMyeZOQV+G7iv+uoEcGzTqkeBN6vyo1uUn2//j2bmvZl575XUT5KkneCKQri6xnvWJ4CzI6efBB6KiHZE3MFsANazmXkSWI2I+6tR0Z8CvnEV9ZYkacdrXGyFiPg94OeBgxFxAvjXwM9HxD3MupS/D/wyQGa+HBFPAK8AY+CzmTmpdvUZZiOtO8A3q5ckSTetmA1W3r4iYntXUJKki3t+q0uszpglSVIhhrAkSYUYwpIkFWIIS5JUiCEsSVIhhrAkSYUYwpIkFWIIS5JUiCEsSVIhhrAkSYUYwpIkFWIIS5JUiCEsSVIhhrAkSYUYwpIuU9BstqnV6qUrIu14hrC0S0TUaDbbQFzX3+l0Fmm15qrfknQ1DGFpl4gIGo0W8/N7iLg+f9qNRotms0291oDM6/Ib0s2kUboCkq6N6XRCRNBuz1Or1en31xmPh9ds/3Nzi3Q6izQaLQaDLpPp5JrtW7pZeSYs7SLT6ZRarUZnbpHFxf0XvW5bq9VpNFqXsOdgbm6BVmuOiCBzek0DXrpZGcLSLjKdjmdBXK8zN7fA0tLBi24TEVzsOvLevQdptzvUavUqhPOi20i6OLujpV1kOOyzsDA7w20229TrdcbjIWtry8DfvoYbUaPZaFFvNAEYjQZb7rPVmmN+folGo0VOJ0wmY0bD/k/sT9Ll80xY2mWGwx4Rs9uI2u15Op09LO058LfWabXmZqOc2x3m55cueMvRXHuBubnZq9FsMxoNGAx7N6Ip0q5nCEu7TL+/znDYJ6JGuz3PwsJe9iwdYGFh78Y6zWZ74xajs+E7O3NubCxDsLR0kIXFfXTmFoFZwM/OmO2Klq4Fu6OlXSNotWbBOhz06PdWaTUPMj+/BCwRUWc6ndLvr1OvNWi1Okwmo9ngrHqTRr25Mdiq0WjRmVvkwIHDtFpz1BtN+oN1RsMB4/EIu6Kla8MQlnaJdhWW0+mUqM06uSbTMbVanbn2Anv27KffXwOgVm/Qbs0xnjSYTie02p3qFqfZhB+NRpM9SwdotztE1JhMxgCMJyP6/XVgdk05c1qmsdIuYQhLu8Rg2GexNUej3iSnU2r1xmxSDYAI2u0OS0u3MD+/tHE/8WQ8YjgaMJmMmEzG7N17kGZzjl5vlbm5Rer1Jo1Gk1qtRre7Qre7Qr1ep15vMBz2yzZY2gUMYWnXSHq9VRYW9kEmve4q7fY88/N7aLXaTCYTarUGi4uLNJstWq0Ovd4a8xEsL5+iXm+ysLCXer3JcNirBnXNrgX3uqucOfMjJpMxc3MLGyOpa7U6tVrdLmrpChnC0o53dpBUkpmMRgPq9QaT6Zheb9b93OkssTC/QK1WYzqdMp1OqNfq1Gq12X2/0wlEjVarQ7PZ5tZbf4q9ew8B8N7yKVZWTzMeD+l0Fomobbym1axZzWbrvLc4STo/Q1jawc6eiZ6dN7rZaNGZ3wPAYDC7jWg47DEc9mi1ZrcsZc7OmKNWZ2npIN3uCodu/SlW3nubyWS8cctSRNDvrzMYdGk0mtWMWZ1q392N94ggHC0tXRFDWNrBarU69VodImg2WhBBvd6g09nD3NwCk8mEhYV9tJptRqMhOZ2wtHRwNuK53mB+fomlpVtmA68yqdXrzM8v8e67b9Kq7jPef+Awy8unmJ/fS61WYzDokZk0Gi1qUWMw6DL0LFi6IoawtIONx0NotJibW4AIGo0mg0GXWtRozy3QbjdoNmdzQ8/NLTCdjGk0W9TqDebmFmi2m7Q7c/TWuyzu2U8taswv7GUw6NKurv0Oq4k5ms0Ww0GP6XRCs9lmPB4yncxGTDtKWroyTtYh7XDj8ZDRsE/mlHqtQb3eZDDoMZmMZyOlM1lZeZe1tTPs238b+w/eQkTQ662yvrrCdDyhXm8wGvZpttosLu3lwIHbIZPl5VOsr79HqzVHuz0P/Hhyj8lkzHgyKtl0acfzTFja0WZnv0kynU4YVKOaa/U6rdYcSTIc9je6rdfWztCaa2/cCzy30GZuYQ4iGA57dLurDHp9er1Vmq02zWaLXm+VzCnNZptWu8Okt8Z0OqHRaDKdTqqHOUi6EoawtKMlmVMajTnm5haoVfcFt1odGvUm9UaTbneF97//Q8zP76XbfY+1lffodBa5/advZ9+t+xgNR5BJe67F6R+dobu2Tr+/xuLi0dkZMWw8OanT2UNmEhEMh/2N0dGSrowhLO1wk8mYbneVzJw9ZKHe3Piu319n795DtNvz1OsNBoMee/bcwsGjB9l/+h1+5r/8f6wtLvHa37ub8XjCgfftZ89wD+vry6yvr9Bstti//zCnT7/JeDSk3mhSi9ltTmcn62g222ROZzN1VWHts4alS2MIS7tAo9Gsngs8k9MJvf4atVqd/ftvo9XqMBh02b/vNt7/offzkRee4YE/+Y+8dtc/4LYffJdPPPsn/Nv/6Z8y3neQpbkm4+FPc+L1v2FhYR/j8YDl5R8xmY5hDPVGc/aEplaHabaYTCbVIxNHTKcTRqMh0+nEs2TpEhjC0o4yuwbcbndoNudoteao1eobczvXajX27j20MWvWaDyk212lXmswnow4cPTvMN9d52NPPclv/fL/xs/9r/8zB24/TPPX/gX/y6vf4kt/514OHj3Ivlv3cvzVLuvrKwDs2XMLnc4S/d4ak94Kk8mIfftvm4U/wWDYo9tdYTIZ025PmE73zEZWD3oMRwNHT0vnYQhL29TGJBz1Ju32PI1mi0ajRb3eICKo1Rp0OouzgVi1GuPxiEF/nXZ7fnb7UbPFZDJmYWEfmVM6nUUOHj3Isddf441Dh1m4/2e5be9e7jh0K6cf/hQ/8yu/yt//5C/x1n99i/Z8myPHPsCZd07RaLbITFZXT9NstdnXvpWlpYNMJuONrudmTllc3D+7bWk62eiqrtXqNJqzs+XxeFhtYyBLZxnC0jbUbnVotuY2ZsGKWo16vUGtVmOuvUCj2aJea5Ako1GfubkFFhb2sn//bVXYjRiNBuzbd4hDhw+T0ynTSTIZT3lzfi+HT5/iZ+54H3cfOQLAgZdfoffBD7Lv1r2snl7hzKllmu0mBw7dxnQ6pbc2mzmrXqsz19lDTid0e6usr79H5nR2e1R99kSms+Ec1X3LABFjsuqidjC19GOGsLQNTXO6cV11Mh3TqLU2bgUaT0ZMc0otatQbs+cATyYTRqM+g0F3Njir0aquE9dYfufd2dlyo0aj1WD1fYd56Y6/z/3//J/T/Wf/jPapU0y//GXOPPEE7f6YTLj9p2/njVd/QHelS6+3RrPZ5sCBw6yunqbbfY9ud5XpZAwRTCZjRqMB3e7KRvjO5qeeMhoNGI0GG9NcSvrbDGFpGzobXgCNRmvj2m+tVofMjdmx5jt7aLbmyMyNe3b377+NTmd2C9Hqymn27ruVwaDLZDLh0LFDtDtt/uMv/lPW3vwrPvRHf8TawYO8/R/+A9+ewNqZNbqrXerNOq25Fmtry6yvL2/UaX39PSaTMb3eKpPJpArc8exsubdWPRJxsvFoRAdnSRcW2/1G+4jY3hWUCoiYdU83Gi1qtRoL83upN5q0WnOzkcvteTKnjMcj3ve+OzYexnDs2N/j/f/NHUwnUxb3LVBvNnjnxDs0Wg3IpNmeXb9dfusMZ95a5t13fwjMHtRw9jao4bBHr7e2cRvScNhjNBzQH6zP5qfOxMcaSj/h+cy899xCz4SlHWgWsEPG4yGNRovJeDZ9ZLs9z8LiPmA2nWVEjeGwx4EDt9NotHjnnTc5tHyIg0cPbdzSNBmNWTuzxtLBJcajMYPegDNvLfPOOyeICOr1Js3mHLBKvV6n0WiR0wmrq6c3zr5nj1AcOuhKukyGsLTDnX2IQ6ezSKvdod9fZzjsMZlMaLc7dLur7N07ptls0++v89q3X2LY/yAHDt9CZ3GO5lyL8XjCsDegt9Zn+e0znD59kmazzXQ6ZWXlnY3rurPgHzEaD6lFjdHE673S1TCEpR0uokatVmM46LG2trzRTb24uI/MZDIesb6+TKvVmY2int/L6pk11pbXGfT6rK8vb0y40eksMRoNGI9H1OtNcjohoka3u0qrNcd4PJwNzlp/j8l0dk04ouYZsHSFDGFph8v88RSSwOxBDoMu9Xp9NuPVZMTbb7/Brbf+FIcOHWNh7wLr762ztnaGZrNNt7tKo9FkZeUdMpOlpYNkTun31njn3R+S0ylZ7fPs2TARTk0pXQOGsLRrBHNzC/T768BsHum9ew9tzOncas3R76+TmfT76xu3HgGMxyOGwz5vv/0GtVp99tSk1myA1+nTJ4HZHNXD4Y9vRWq15hgOBzgIS7pyhrC0ayT9/howu60posZ7771NrVbnlltu35ixajweMhh0WVl5h35/nTNnfkS/v8Zg0OPAgcOsrS0zmYzozC3SarbZs+cA7777Q/r9dcbj2S1IkAyH3n4kXa1a6QpIuj5WV08zGg2o1Wo0m7NQjqgxGg1ZWXmH9bVl3n77B6ytnWE8HjE3twiwcd/vaDxkVAV2o9GaTYs56NJstgq3TNo9LhrCEXEsIv40Il6NiJcj4leq8gMR8VREfK96379pmy9ExPGIeC0iPrqp/MMR8WL13Zdi82NfJF0z4/FwY+as0Wh27Taixng8pNt9j15vjXpjNtNWrVan0WhSr9erCTlWNrqtm41WNenGlGazReZ0YxIRSVfvUs6Ex8CvZebPAPcDn42Iu4DPA09n5p3A09Vnqu8eAu4GHgC+HBH1al9fAR4B7qxeD1zDtkiqtFud2exawC233E6ns8R0OmZ5+RQrK+8yHPaZjEcMhz3W199jfX2FXm+N1dXTrK6+y/KZH9HtzkZALy7uY35+D7WoE2HnmXQtXfQvKjNPZuZfVsurwKvAEeBB4LFqtceAj1fLDwKPZ+YgM18HjgP3RcRhYCkzn8nZlDpf37SNpGsoqscbLizsrQZQ9eh2Vzce7DAc9hgMuvR6q6ysvMPy8lsb14zX1s5UD4YYVo8nnNBqdZhfWDKEpWvssv6iIuL9wM8C3wJuy8yTMAtq4NZqtSPAG5s2O1GVHamWzy2XdB10Oos0m+1Njxyc3f/b660yHPSZTCeMq5m2gNncz+MRa6tnGAx6DPrr5HTKYNBlNOrTqCYEkXTtXPLo6IhYBH4f+NXMXLnA5dytvsgLlG/1W48w67aWdJlm80rPru9OJ2NyOmUyGTGuJu2Y3cIEjGAy+XEIz6a5DAbDPrG2TESNA9UgrFqtzmTSmz1asdn2urB0jVxSCEdEk1kA/25m/kFV/FZEHM7Mk1VX86mq/ARwbNPmR4E3q/KjW5T/hMx8FHi0+m1vQpQuQ6PRrO7j7dBotqB6eMNg0GNt9QyjTZNsbJ7kY/OTmzKndLvv0Wg0aTbb1OtL1Gt16vWGISxdQ5cyOjqA3wFezczf3PTVk8DD1fLDwDc2lT8UEe2IuIPZAKxnqy7r1Yi4v9rnpzZtI+kaGY9HtFodWq05plV3cr+/zvraMsPRgPF4WE3WsXrefZydhWt9bZled5XBoEvU6rTbHdrt+RvYGml3u5Qz4Y8AvwS8GBEvVGX/Evg3wBMR8WngB8AnATLz5Yh4AniF2cjqz2bm2bv6PwN8DegA36xekq6hdnueZqO18Yzhsw9f6PZWNp68dDHT6Wwu6P5gnWk1e1az2WYyHtFstj0blq4Rnycs7TKdzh7m5/fMwrjZptdbo99fp9tdueT5ns/e3jSdTqnVauzde2hjUNZ0OmV19V3W19+7bm2QdiGfJyzdDBqNJrVaY2OO6G53hW53hen00qeZ3LzudDphbe0MMDvLzpxSr/lPh3Qt+Jck7TIRNabTMYPB7KEM/f76ZQXwVkajAf3+OhFBvd6k3mgyu+HBjirpahjC0i5y9qlIw2F/42EN1+qRg73eKqPRgIWFvQC0Wu2/Nbpa0uUzhKVdJCIYDLobAXytjcdDhsN+Ndd0EzCEpathCEu7yNmHLWROr9tv9HprtFtzTK6yi1uSISztKpPJ+Ab8SjIY9vF6sHT1nI1d0hUwgKVrwRCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIuGsIRcSwi/jQiXo2IlyPiV6ryL0bEDyPiher1sU3bfCEijkfEaxHx0U3lH46IF6vvvhQRcX2aJUnS9te4hHXGwK9l5l9GxB7g+Yh4qvrutzLz321eOSLuAh4C7gZuB/4kIj6YmRPgK8AjwH8G/hh4APjmtWmKJEk7y0XPhDPzZGb+ZbW8CrwKHLnAJg8Cj2fmIDNfB44D90XEYWApM5/JzAS+Dnz8qlsgSdIOdVnXhCPi/cDPAt+qij4XEd+JiK9GxP6q7AjwxqbNTlRlR6rlc8slSbopXXIIR8Qi8PvAr2bmCrOu5Q8A9wAngd84u+oWm+cFyrf6rUci4rmIeO5S6ydJ0k5zSSEcEU1mAfy7mfkHAJn5VmZOMnMK/DZwX7X6CeDYps2PAm9W5Ue3KP8JmfloZt6bmfdeTmMkSdpJLmV0dAC/A7yamb+5qfzwptU+AbxULT8JPBQR7Yi4A7gTeDYzTwKrEXF/tc9PAd+4Ru2QJGnHuZTR0R8Bfgl4MSJeqMr+JfCLEXEPsy7l7wO/DJCZL0fEE8ArzEZWf7YaGQ3wGeBrQIfZqGhHRkuSbloxG6i8fUXE9q6gJEkX9/xWl1idMUuSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIMYUmSCjGEJUkqxBCWJKkQQ1iSpEIuGsIRMRcRz0bEtyPi5Yj49ar8QEQ8FRHfq973b9rmCxFxPCJei4iPbir/cES8WH33pYiI69MsSZK2v0s5Ex4A/ygz/1vgHuCBiLgf+DzwdGbeCTxdfSYi7gIeAu4GHgC+HBH1al9fAR4B7qxeD1zDtkiStKNcNIRzZq362KxeCTwIPFaVPwZ8vFp+EHg8MweZ+TpwHLgvIg4DS5n5TGYm8PVN20iSdNO5pGvCEVGPiBeAU8BTmfkt4LbMPAlQvd9arX4EeGPT5ieqsiPV8rnlkiTdlC4phDNzkpn3AEeZndV+6AKrb3WdNy9Q/pM7iHgkIp6LiOcupX6SJO1ElzU6OjOXgT9jdi33raqLmer9VLXaCeDYps2OAm9W5Ue3KN/qdx7NzHsz897LqZ8kSTvJpYyOPhQR+6rlDvALwHeBJ4GHq9UeBr5RLT8JPBQR7Yi4g9kArGerLuvViLi/GhX9qU3bSJJ002lcwjqHgceqEc414InM/KOIeAZ4IiI+DfwA+CRAZr4cEU8ArwBj4LOZOan29Rnga0AH+Gb1kiTpphSzgcrbV0Rs7wpKknRxz291idUZsyRJKsQQliSpEENYkqRCDGFJkgoxhCVJKsQQliSpEENYkqRCDGFJkgoxhCVJKsQQliSpEENYkqRCDGFJkgoxhCVJKsQQliSpEENYkqRCDGFJkgoxhCVJKsQQliSpEENYkqRCDGFJkgoxhCVJKsQQliSpEENYkqRCDGFJkgoxhCVJKsQQliSpEENYkqRCDGFJkgoxhCVJKsQQliSpEENYkqRCDGFJkgoxhCVJKsQQliSpEENYkqRCDGFJkgoxhCVJKsQQliSpEENYkqRCDGFJkgoxhIp6QuoAAAbQSURBVCVJKsQQliSpEENYkqRCDGFJkgoxhCVJKsQQliSpEENYkqRCDGFJkgoxhCVJKuSiIRwRcxHxbER8OyJejohfr8q/GBE/jIgXqtfHNm3zhYg4HhGvRcRHN5V/OCJerL77UkTE9WmWJEnbX+MS1hkA/ygz1yKiCfy/EfHN6rvfysx/t3nliLgLeAi4G7gd+JOI+GBmToCvAI8A/xn4Y+AB4JtIknQTuuiZcM6sVR+b1SsvsMmDwOOZOcjM14HjwH0RcRhYysxnMjOBrwMfv7rqS5K0c13SNeGIqEfEC8Ap4KnM/Fb11eci4jsR8dWI2F+VHQHe2LT5iarsSLV8brkkSTelSwrhzJxk5j3AUWZntR9i1rX8AeAe4CTwG9XqW13nzQuU/4SIeCQinouI5y6lfpIk7USXNTo6M5eBPwMeyMy3qnCeAr8N3FetdgI4tmmzo8CbVfnRLcq3+p1HM/PezLz3cuonSdJOcimjow9FxL5quQP8AvDd6hrvWZ8AXqqWnwQeioh2RNwB3Ak8m5kngdWIuL8aFf0p4BvXsC2SJO0olzI6+jDwWETUmYX2E5n5RxHx7yPiHmZdyt8HfhkgM1+OiCeAV4Ax8NlqZDTAZ4CvAR1mo6IdGS1JumnFbKDy9hUR27uCkiRd3PNbXWJ1xixJkgoxhCVJKsQQliSpEENYkqRCDGFJkgoxhCVJKsQQliSpEENYkqRCDGFJkgoxhCVJKsQQliSpEENYkqRCDGFJkgoxhCVJKsQQliSpEENYkqRCDGFJkgoxhCVJKsQQliSpEENYkqRCDGFJkgoxhCVJKsQQliSpEENYkqRCGqUrcAnWgNdKV+IGOAi8U7oSN4Dt3F1s5+5iO6+fn9qqcCeE8GuZeW/pSlxvEfGc7dw9bOfuYjt3l+3UTrujJUkqxBCWJKmQnRDCj5auwA1iO3cX27m72M7dZdu0MzKzdB0kSbop7YQzYUmSdqVtG8IR8UBEvBYRxyPi86Xrc7Ui4vsR8WJEvBARz1VlByLiqYj4XvW+f9P6X6ja/lpEfLRczS8sIr4aEaci4qVNZZfdroj4cPXf53hEfCki4ka35ULO084vRsQPq2P6QkR8bNN3O7WdxyLiTyPi1Yh4OSJ+pSrfVcf0Au3cVcc0IuYi4tmI+HbVzl+vynfb8TxfO7f/8czMbfcC6sBfA38XaAHfBu4qXa+rbNP3gYPnlP1b4PPV8ueB/6Navqtqcxu4o/pvUS/dhvO06x8CPwe8dDXtAp4F/jsggG8C/6R02y6hnV8E/sUW6+7kdh4Gfq5a3gP8VdWeXXVML9DOXXVMqzotVstN4FvA/bvweJ6vndv+eG7XM+H7gOOZ+TeZOQQeBx4sXKfr4UHgsWr5MeDjm8ofz8xBZr4OHGf232Tbycw/B06fU3xZ7YqIw8BSZj6Ts7+Cr2/aZls4TzvPZye382Rm/mW1vAq8Chxhlx3TC7TzfHZqOzMz16qPzeqV7L7jeb52ns+2aed2DeEjwBubPp/gwn8gO0EC/ykino+IR6qy2zLzJMz+UQBurcp3evsvt11HquVzy3eCz0XEd6ru6rNderuinRHxfuBnmZ1V7Npjek47YZcd04ioR8QLwCngqczclcfzPO2EbX48t2sIb9UHv9OHcX8kM38O+CfAZyPiH15g3d3Yfjh/u3Zqe78CfAC4BzgJ/EZVvuPbGRGLwO8Dv5qZKxdadYuyHdPWLdq5645pZk4y8x7gKLOzvQ9dYPXd1s5tfzy3awifAI5t+nwUeLNQXa6JzHyzej8F/CGz7uW3qu4PqvdT1eo7vf2X264T1fK55dtaZr5V/eFPgd/mx5cMdnQ7I6LJLJh+NzP/oCredcd0q3bu1mMKkJnLwJ8BD7ALj+dZm9u5E47ndg3hvwDujIg7IqIFPAQ8WbhOVywiFiJiz9ll4B8DLzFr08PVag8D36iWnwQeioh2RNwB3MlssMBOcVntqrrDViPi/mok4qc2bbNtnf1HrPIJZscUdnA7q3r9DvBqZv7mpq921TE9Xzt32zGNiEMRsa9a7gC/AHyX3Xc8t2znjjie13PU19W8gI8xG7H418C/Kl2fq2zL32U2Eu/bwMtn2wPcAjwNfK96P7Bpm39Vtf01ttEoxC3a9nvMunlGzP4v8tNX0i7gXmZ/IH8N/J9UE8lsl9d52vnvgReB7zD7oz68C9r53zPrfvsO8EL1+thuO6YXaOeuOqbAPwD+S9Wel4D/vSrfbcfzfO3c9sfTGbMkSSpku3ZHS5K06xnCkiQVYghLklSIISxJUiGGsCRJhRjCkiQVYghLklSIISxJUiH/P1yVkxgdVkS5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,8))\n",
    "plt.imshow(data[0], cmap='bone', norm=colors.LogNorm(vmin=data[0].max()/50000, vmax=data[0].max()/10))\n",
    "c = plt.Circle((cx, cy), aperad, ec='red', fc='None')\n",
    "plt.gca().add_artist(c)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 204,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7826989890559137\n"
     ]
    }
   ],
   "source": [
    "flux = Photometry_Aperture.A_photometry([0], cx=cx, cy=cy, r=[aperad], methods=[method], img_data=data)\n",
    "contam = flux[0]/np.sum(data)\n",
    "print(contam)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### WASP-12b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 205,
   "metadata": {},
   "outputs": [],
   "source": [
    "R_comp = 0.56 # Stevenson+2014a\n",
    "Teff_comp = 3660\n",
    "logg_comp = 5.0 # median value from https://arxiv.org/pdf/1810.13252.pdf\n",
    "feh_comp = 0 # median value from https://arxiv.org/pdf/1810.13252.pdf\n",
    "N_comp = 2\n",
    "\n",
    "R_host = 1.657\n",
    "Teff_host = 6360\n",
    "logg_host = 4.17\n",
    "feh_host = 0.30"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### WASP-103b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {},
   "outputs": [],
   "source": [
    "R_comp = 0.665 # https://en.wikipedia.org/wiki/61_Cygni\n",
    "Teff_comp = 4400\n",
    "logg_comp = 4.40 # https://en.wikipedia.org/wiki/61_Cygni\n",
    "feh_comp = -0.20 # https://en.wikipedia.org/wiki/61_Cygni\n",
    "N_comp = 1\n",
    "\n",
    "R_host = 1.436\n",
    "Teff_host = 6110\n",
    "logg_host = 4.22\n",
    "feh_host = 0.06"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CoRoT-2b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 191,
   "metadata": {},
   "outputs": [],
   "source": [
    "R_comp = 0.65 # https://en.wikipedia.org/wiki/HD_113538\n",
    "Teff_comp = 4000\n",
    "logg_comp = 4.28 # median value from https://arxiv.org/pdf/1810.13252.pdf\n",
    "feh_comp = -0.17 # median value from https://arxiv.org/pdf/1810.13252.pdf\n",
    "N_comp = 1\n",
    "\n",
    "R_host = 0.902\n",
    "Teff_host = 5620\n",
    "logg_host = 4.66\n",
    "feh_host = -0.03"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "metadata": {},
   "outputs": [],
   "source": [
    "phoenixPath = '/home/taylor/Documents/Research/PHOENIX/HiResFITS/'\n",
    "phoenixWavFile = '/home/taylor/Documents/Research/PHOENIX/HiResFITS/WAVE_PHOENIX-ACES-AGSS-COND-2011.fits'\n",
    "\n",
    "with fits.open(phoenixWavFile) as f:\n",
    "    wavStar = f[0].data*1e-4 # angstrom to micron"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "metadata": {},
   "outputs": [],
   "source": [
    "fStars = []\n",
    "\n",
    "for Rstar, teff, logg, feh, Nstar in [[R_host, Teff_host, logg_host, feh_host, 1], [R_comp, Teff_comp, logg_comp, feh_comp, N_comp]]:\n",
    "\n",
    "    # Get the phoenix file ready to compute the stellar brightness temperature\n",
    "    teffStr = teff\n",
    "    if teffStr <= 7000:\n",
    "        teffStr = teffStr - (teffStr%100) + np.rint((teffStr%100)/100)*100\n",
    "    elif teffStr > 7000:\n",
    "        teffStr = teffStr - (teffStr%200) + np.rint((teffStr%200)/200)*200\n",
    "    elif teffStr > 12000:\n",
    "        teffStr = 12000\n",
    "    teffStr = str(int(teffStr)).zfill(5)\n",
    "\n",
    "    if np.isnan(logg):\n",
    "        logg = 4.5\n",
    "    logg = logg - (logg%0.5) + np.rint((logg%0.5)*2)/2.\n",
    "    logg = -logg\n",
    "    \n",
    "    if np.isnan(feh):\n",
    "        feh = 0.\n",
    "    feh = (feh - (feh%0.5) + np.rint((feh%0.5)*2)/2.)\n",
    "    if feh<-2.:\n",
    "        feh = (feh - (feh%1) + np.rint((feh%1)))\n",
    "        \n",
    "    webfile = ('lte'+teffStr\n",
    "           +(\"{0:+.02f}\".format(logg) if logg!=0 else '-0.00')\n",
    "           +(\"{0:+.01f}\".format(feh) if feh!=0 else '-0.0')\n",
    "           +'.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits')\n",
    "    \n",
    "    webfolder = 'ftp://phoenix.astro.physik.uni-goettingen.de/HiResFITS/'\n",
    "    webfolder += 'PHOENIX-ACES-AGSS-COND-2011/Z'+(\"{0:+.01f}\".format(feh) if feh!=0 else '-0.0')+'/'\n",
    "    \n",
    "    phoenixSpectraFile = phoenixPath+webfile\n",
    "    \n",
    "    if not os.path.exists(phoenixSpectraFile):\n",
    "        print('Downloading relevant PHOENIX spectra...', flush=True)\n",
    "        try:\n",
    "            _ = urllib.request.urlretrieve(webfolder+webfile, phoenixSpectraFile)\n",
    "        except:\n",
    "            # Fix: throw a proper error\n",
    "            print('ERROR: No previously downloaded PHOENIX data - try again when you are connected to the internet.',\n",
    "                  flush=True)\n",
    "            break\n",
    "        print('Done download.', flush=True)\n",
    "    \n",
    "    with fits.open(phoenixSpectraFile) as f:\n",
    "        fStars.append(f[0].data*1e-1*Rstar**2*Nstar) # 'erg/s/cm^2/cm' to kg/s^3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 208,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('SpitzerResponseCurves/Spitzer_Ch2_ResponseCurve_Subarray.txt', 'r') as f:\n",
    "# with open('SpitzerResponseCurves/Spitzer_Ch2_ResponseCurve_FullFrame.txt', 'r') as f:\n",
    "    f.readline()\n",
    "    f.readline()\n",
    "    f.readline()\n",
    "    lines = f.readlines()\n",
    "\n",
    "lines = np.array([line.strip().split('\\t') for line in lines]).astype(float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 209,
   "metadata": {},
   "outputs": [],
   "source": [
    "fluxes = []\n",
    "\n",
    "wavTemp = wavStar[np.logical_and(wavStar>=lines[0][0], wavStar<=lines[-1][0])]\n",
    "for f in fStars:\n",
    "    fTemp = f[np.logical_and(wavStar>=lines[0][0], wavStar<=lines[-1][0])]\n",
    "    fTemp *= lines[[np.argmin(np.abs(lines[:,0]-wavTemp[i])) for i in range(len(wavTemp))],1]\n",
    "    \n",
    "    fluxes.append(fTemp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 210,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha_comp = np.mean(fluxes[1]/fluxes[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 211,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7826989890559137\n",
      "0.11607887\n",
      "0.09085481303223122\n"
     ]
    }
   ],
   "source": [
    "print(contam)\n",
    "print(alpha_comp)\n",
    "print(contam*alpha_comp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:python3] *",
   "language": "python",
   "name": "conda-env-python3-py"
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}