{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analytic maximum likelihood scaled transit fit\n",
    "\n",
    "We can use a factor $s$ to modifiy the transit depth in a nominal model ${\\bf m_0}$ calculated with approximately the correct depth that is scaled as follows ...\n",
    "\n",
    "$ {\\bf m}(s) = 1 + s\\times({\\bf m_0}- 1)$\n",
    "\n",
    "The data are normalised fluxes ${\\bf f} = f_1, \\dots, f_N$ with nominal errors ${\\bf \\sigma} = \\sigma_1,\\dots, \\sigma_N$. \n",
    "\n",
    "Assume that the actual standard errors are underestimated by some factor $b$, and that these are normally distributed and uncorrelated, so that the log-likelihood is\n",
    "\n",
    "$\\ln {\\cal L}  = -\\frac{1}{2b^2}\\chi^2 - \\frac{1}{2}\\sum_{i=1}^N \\ln \\sigma_i^2  - N\\ln b - \\frac{N}{2}\\ln(2\\pi)$\n",
    "\n",
    "where\n",
    "\n",
    "$\\chi^2 = \\sum_i^N \\frac{(f_i - 1 - s(m_{0,i}-1)^2}{\\sigma_i^2} $\n",
    "\n",
    "\n",
    "The maximum likelihood occurs for parameter values $\\hat{s}$, and $\\hat{b}$ such that\n",
    "\n",
    "$\\left. \\frac{\\partial  \\ln {\\cal L}}{\\partial s}\\right|_{\\hat{s},\\hat{b}} = 0$\n",
    "and \n",
    "$\\left. \\frac{\\partial  \\ln {\\cal L}}{\\partial b}\\right|_{\\hat{s},\\hat{b}} = 0$,\n",
    "\n",
    "from which we obtain \n",
    "\n",
    " $\\hat{s} = \\sum_{i=1}^N \\frac{(f_i - 1)(m_{0,i}-1) }{\\sigma_i^2} \\left[ \\sum_{i=1}^N \\frac{(m_{0,i}-1)^2}{\\sigma_i^2}\\right]^{-1}$\n",
    " \n",
    " and \n",
    " \n",
    " $\\hat{b} = \\sqrt{\\chi^2/N}$\n",
    " \n",
    " For the standard errors on these parameters we use \n",
    " \n",
    " $\\sigma_s^{-2}  = -\\frac{\\partial^2\\ln{\\cal L}}{\\partial^2 s^2} \\left.\\right|_{\\hat{s},\\hat{b}}$ \n",
    " \n",
    " and\n",
    " \n",
    " $\\sigma_b^{-2}  = -\\frac{\\partial^2\\ln{\\cal L}}{\\partial^2 b^2} \\left.\\right|_{\\hat{s},\\hat{b}}$ \n",
    "\n",
    "to derive\n",
    "\n",
    "$\\sigma_s = b\\left[\\sum_{i=1}^N \\frac{(m_i-1)^2}{\\sigma_i^2}\\right]^{-1/2}$\n",
    " \n",
    " and\n",
    " \n",
    " $\\sigma_b  = \\left[3\\chi^2/b^4 -N/b^2\\right]^{-1/2} .$\n",
    " \n",
    "**N.B.** Whether or how much of the data outside transit to include depends on whether these data can be assumed to have the same noise characteristics as the data in transit. Note that including these data has no effect on $s$ or $\\sigma_s$, because of the factors $(m_{0,i}-1)$ in their calculation, but will affect the estimates of $b$ and $\\sigma_b$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "ename": "ImportError",
     "evalue": "cannot import name 'scaled_transit_fit'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mImportError\u001b[0m                               Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-1-adf5e034a187>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0mget_ipython\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun_line_magic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'pylab'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'inline'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mpycheops\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmodels\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mTransitModel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mscaled_transit_fit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mpycheops\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minstrument\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mtransit_noise\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mImportError\u001b[0m: cannot import name 'scaled_transit_fit'"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "from pycheops.models import TransitModel, scaled_transit_fit\n",
    "from pycheops.instrument import transit_noise\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load simulated data and clean out CRs\n",
    "\n",
    " As a heuristic approach to setting a limit to reject cosmic rays we use the most negative value in the observed data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Read 2615 observations\n",
      "C_cut = 3515101, N(C > C_cut = 85)\n",
      "Mean counts = 3507779.2\n",
      "Median counts = 3507211.2\n",
      "RMS counts = 3790.9 [1081 ppm]\n",
      "Median standard error = 1954.5 [557 ppm]\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 44372 characters]\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11dd62dd8>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "T = Table.read('CH_PR300005_TG000201_TU2020-01-02T23-44-04_SCI_COR_Lightcurve-OPTIMAL_V0000.fits')\n",
    "N_obs = len(T)\n",
    "print('Read {} observations'.format(N_obs))\n",
    "ok = (T['EVENT'] == 0)  |  (T['EVENT'] == 100)\n",
    "BJD = np.array(T['BJD_TIME'][ok])\n",
    "C = np.array(T['FLUX'][ok])\n",
    "Cerr = np.array(T['FLUXERR'][ok])\n",
    "\n",
    "C_cut = (2*median(C)-min(C))\n",
    "Ecr = C[C > C_cut] - median(C)\n",
    "E_cut = C_cut - median(C)\n",
    "N_cut = len(Ecr)\n",
    "Ecr.sort()\n",
    "print('C_cut = {:0.0f}, N(C > C_cut = {})'.format(C_cut, N_cut))\n",
    "\n",
    "print('Mean counts = {:0.1f}'.format(C.mean()))\n",
    "print('Median counts = {:0.1f}'.format(median(C)))\n",
    "print('RMS counts = {:0.1f} [{:0.0f} ppm]'.format(std(C), 1e6*std(C)/np.median(C)))\n",
    "print('Median standard error = {:0.1f} [{:0.0f} ppm]'.format(np.median(Cerr), 1e6*np.median(Cerr)/np.median(C)))\n",
    "\n",
    "time_raw = BJD[C < C_cut]-int(min(BJD))\n",
    "flux_raw = C[C < C_cut]/median(C)\n",
    "flux_err = Cerr[C < C_cut]/median(C)\n",
    "\n",
    "fig,ax=subplots(1,2,figsize=(8,4))\n",
    "ax[0].plot(time_raw, flux_raw,'b.',ms=1)\n",
    "ax[0].set_xlabel(\"BJD-{}\".format(int(min(BJD))))\n",
    "ax[0].set_ylabel(\"Flux\".format(median(C)));\n",
    "\n",
    "n, bins, patches = ax[1].hist(flux_raw, 50, normed=True)\n",
    "ax[1].set_xlabel(\"Flux\".format(median(C)));\n",
    "v  = np.var(flux_raw)\n",
    "ax[1].plot(bins,exp(-0.5*(bins-1)**2/v)/sqrt(2*np.pi*v));\n",
    "ax[1].plot(bins,exp(-0.5*(bins-1)**2/v)/sqrt(2*np.pi*v));\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**N.B.** we have normalized the data using the median of the data here. There is a noticable offset between the mean and median values due to the tail in the distribution of the observed counts values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Inject a transit\n",
    "\n",
    "This will be the reference model ${\\bf m_0}$. Use a transit depth $D = 85\\,$ppm and a transit width $W=3\\,$hours based on the parameters for SciReq 1.2. For simplicity using impact parameter $b=0$. Use limb darkening appropriate for a K-type dwarf and an orbital period of 13 days star since this is the case under consideration for SciReq 1.2.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R_p/R_* = 0.0092, S = 0.96\n",
      "h_1 = 0.680, h_2 = 0.432\n"
     ]
    }
   ],
   "source": [
    "WIDTH_HOURS = 3\n",
    "DEPTH_PPM = 85\n",
    "PERIOD = 13\n",
    "\n",
    "\n",
    "T_0 = 1.5\n",
    "D = DEPTH_PPM*1e-6  # Set transit depth\n",
    "W = WIDTH_HOURS/24/PERIOD  # Transit Width in phase units\n",
    "b = 0 # Impact parameter S is the ratio (W_23 / W_14)**2 = ((1-k**2)-b**2)/((1+k)**2-b**2) \n",
    "k = sqrt(D)  # Ratio of the radii\n",
    "S = ((1-k)**2-b**2)/((1+k)**2-b**2)  # Transit Shape\n",
    "print('R_p/R_* = {:0.4f}, S = {:0.2f}'.format(k, S))\n",
    "c = 1.0    # Flux outside transit\n",
    "\n",
    "T_eff = 4800\n",
    "log_g = 4.5\n",
    "Fe_H = 0\n",
    "p2C = stagger_power2_interpolator(passband='CHEOPS')\n",
    "c2,a2,h1,h2 = p2C(T_eff, log_g, Fe_H)\n",
    "print('h_1 = {:0.3f}, h_2 = {:0.3f}'.format(h1, h2))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 25124 characters]\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1236930f0>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "model = TransitModel()\n",
    "pars0 = model.make_params(T_0=T_0, P=PERIOD, D=D, W=W, S=S, h_1=h1, h_2=h2)\n",
    "model_0 = model.eval(params=pars0, t=time_raw) \n",
    "\n",
    "fig,ax=subplots(1,1,figsize=(5,4))\n",
    "ax.plot(time_raw, model_0)\n",
    "ax.set_xlim(min(time_raw),max(time_raw))\n",
    "\n",
    "flux_sim = flux_raw * model_0\n",
    "ax.plot(time_raw, flux_sim,'b.',ms=1)\n",
    "ax.set_xlabel(\"BJD-{}\".format(int(min(BJD))))\n",
    "ax.set_ylabel(\"Flux\".format(median(C)));\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit the transit using only in-transit data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Transit scale factor s = 0.78 +/- 0.62\n",
      "Error scale factor  b = 1.35 +/- 0.05\n",
      "Depth = 66 +/- 53 ppm\n"
     ]
    }
   ],
   "source": [
    "j = (np.abs(time_raw-T_0) < (WIDTH_HOURS/24)).nonzero()[0]\n",
    "\n",
    "s0, b0, sigma_s0, sigma_b0 = scaled_transit_fit(flux_sim[j],flux_err[j],model_0[j])\n",
    "print('Transit scale factor s = {:0.2f} +/- {:0.2f}'.format(s0, sigma_s0))\n",
    "print('Error scale factor  b = {:0.2f} +/- {:0.2f}'.format(b0, sigma_b0))\n",
    "depth0 = s0*DEPTH_PPM*1e-6\n",
    "e_depth0 = sigma_s0*DEPTH_PPM*1e-6\n",
    "print('Depth = {:0.0f} +/- {:0.0f} ppm'.format(depth0*1e6, e_depth0*1e6))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Iteration step\n",
    "\n",
    "Now find depth that gives S/N = 1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Noise in 3 hour window at t=1.5 from (S/N=1) transit depth =  52.6 ppm\n",
      "Noise scaling factor = 1.35\n"
     ]
    }
   ],
   "source": [
    "depth_tol = 0.5e-6\n",
    "depth_in = 0\n",
    "e_depth = e_depth0\n",
    "while abs(e_depth-depth_in) > depth_tol:\n",
    "    depth_in = e_depth\n",
    "    k = sqrt(depth_in)\n",
    "    S = ((1-k)**2)/((1+k)**2)  \n",
    "    pars = model.make_params(T_0=T_0, P=PERIOD, D=depth_in, W=W, S=S, h_1=h1, h_2=h2, c=c)\n",
    "    model_1 = model.eval(params=pars, t=time_raw)\n",
    "    s1, b1, sigma_s1, sigma_b1 = scaled_transit_fit(flux_sim[j],flux_err[j],model_1[j])\n",
    "    depth = s1*depth_in\n",
    "    e_depth = sigma_s1*depth_in\n",
    "    \n",
    "print('Noise in {} hour window at t={:0.1f} from (S/N=1) transit depth =  {:0.1f} ppm'\n",
    "      .format(WIDTH_HOURS, T_0,1e6*depth_in))\n",
    "print('Noise scaling factor = {:0.2f}'.format(b1))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Noise estimate in moving window"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a function to do all this in one call"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def transit_noise(time, flux, flux_err, T_0=None, width=3, \n",
    "                  h_1=0.7224, h_2=0.6713, \n",
    "                  tol=0.1, verbose=None):\n",
    "    \"\"\" \n",
    "    Transit noise estimate\n",
    "    \n",
    "    The noise is calculated in a window of width duration in hours \n",
    "    centered at time T_0 by finding the depth of a transit that gives\n",
    "    S/N = 1. \n",
    "    \n",
    "    The transit is calculated from an impact parameter b=0 and orbital period\n",
    "    P days using power-2 limb darkening parameters h_1 and h_2.\n",
    "    Default values for h_1 and h_2 are solar values.\n",
    "    \n",
    "    If T_0 is not specifed that the median value of time is used.\n",
    "    \n",
    "    Returns\n",
    "     ppm_1 = transit noise estimate in ppm\n",
    "     sig_1 - estimate of the noise scaling factor\n",
    "        \n",
    "    The tolerance on ppm_1, tol, is specified in ppm\n",
    "    \"\"\"\n",
    "\n",
    "    if T_0 is None:\n",
    "        T_0 = median(time)\n",
    "\n",
    "    # Use orbital period = 10* data duration so there is certainly 1 transit\n",
    "    P = 10*(max(time)-min(time))\n",
    "    j = (np.abs(time-T_0) < (width/24)).nonzero()[0]\n",
    "    e_depth = median(flux_err[j])/sqrt(len(j))\n",
    "    if verbose:\n",
    "        print('N = {}, initial depth = {:0.1f} ppm'.format(len(j), 1e6*e_depth))\n",
    "        \n",
    "    W = width/24/P   # Transit Width in phase units\n",
    "    model = TransitModel()\n",
    "\n",
    "    depth_tol = tol*1e-6\n",
    "    depth_in = 0 \n",
    "    while abs(e_depth-depth_in) > depth_tol:\n",
    "        depth_in = e_depth\n",
    "        k = sqrt(depth_in)\n",
    "        S = ((1-k)**2)/((1+k)**2)  \n",
    "        pars = model.make_params(T_0=T_0, P=P, D=depth_in, W=W, S=S, h_1=h_1, h_2=h_2)\n",
    "        model_1 = model.eval(params=pars, t=time)\n",
    "        s1, b1, sigma_s1, sigma_b1 = scaled_transit_fit(flux,flux_err,model_1)\n",
    "        e_depth = sigma_s1*depth_in\n",
    "\n",
    "    return 1e6*depth_in, b1\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Noise calculation for K-dwarf limb darkening"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K-dwarf, P=13\n",
      "Noise in 3 hour window from (S/N=1) transit depth = 53.6 +/- 2.1 ppm\n",
      "Noise scaling factor = 1.3158 +/- 0.0004\n"
     ]
    }
   ],
   "source": [
    "T = linspace(0.75, 2.25, 200)\n",
    "N = empty_like(T)\n",
    "F = empty_like(T)\n",
    "\n",
    "for i,_t in enumerate(T):\n",
    "    _n,_f = transit_noise(time_raw, flux_raw, flux_err,\n",
    "                        T_0=_t, width=WIDTH_HOURS,h_1=h1, h_2=h2)\n",
    "    N[i] = _n\n",
    "    F[i] = _f\n",
    "\n",
    "print('K-dwarf, P=13')\n",
    "print('Noise in {} hour window from (S/N=1) transit depth = {:0.1f} +/- {:0.1f} ppm'\n",
    "      .format(WIDTH_HOURS,mean(N),np.std(N)))\n",
    "print('Noise scaling factor = {:0.4f} +/- {:0.4f}'.format(mean(F),np.std(F)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Noise calculation for sun-like star,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sun-like\n",
      "Noise in 3 hour window from (S/N=1) transit depth = 54.4 +/- 2.1 ppm\n",
      "Noise scaling factor = 1.3158 +/- 0.0004\n"
     ]
    }
   ],
   "source": [
    "Nsun = empty_like(T)\n",
    "Fsun = empty_like(T)\n",
    "\n",
    "for i,_t in enumerate(T):\n",
    "    _n,_f = transit_noise(time_raw, flux_raw, flux_err, \n",
    "                        T_0=_t, width=WIDTH_HOURS)\n",
    "    Nsun[i] = _n\n",
    "    Fsun[i] = _f\n",
    "\n",
    "print('Sun-like')\n",
    "print('Noise in {} hour window from (S/N=1) transit depth = {:0.1f} +/- {:0.1f} ppm'\n",
    "      .format(WIDTH_HOURS,mean(Nsun),np.std(Nsun)))\n",
    "print('Noise scaling factor = {:0.4f} +/- {:0.4f}'.format(mean(Fsun),np.std(Fsun)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 28920 characters]\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11dfe0b70>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax=subplots(2,1,figsize=(6,4),sharex=True)\n",
    "ax[0].set_xlim(min(time_raw),max(time_raw))\n",
    "ax[0].plot(time_raw, flux_raw,'b.',ms=1)\n",
    "ax[0].set_ylabel(\"Flux \".format(median(C)))\n",
    "ax[0].set_ylim(0.9975,1.0035)\n",
    "ax[0].errorbar(median(T),1.003,xerr=WIDTH_HOURS/24,\n",
    "           capsize=5,color='b',ecolor='b')\n",
    "ax[1].plot(T,N,'.',ms=1)\n",
    "ax[1].set_ylabel(\"Noise [ppm] \".format(median(C)))\n",
    "ax[1].set_xlabel(\"BJD-{}\".format(int(min(BJD))));\n",
    "ax[1].plot(T,Nsun,'.',ms=1)\n",
    "ax[1].set_ylabel(\"Noise [ppm] \".format(median(C)))\n",
    "ax[1].set_xlabel(\"BJD-{}\".format(int(min(BJD))));\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Same test but with a gap\n",
    " Gaps are are shorter than width because I have not made the function transit_noise robust to empty input data yet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 23380 characters]\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x123618588>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "T = linspace(0.75, 2.25, 200)\n",
    "Ngap = empty_like(T)\n",
    "Fgap = empty_like(T)\n",
    "\n",
    "i_gap = ((abs(time_raw-1.5) > 0.05).nonzero())[0]\n",
    "time_gap = time_raw[i_gap]\n",
    "flux_gap = flux_raw[i_gap]\n",
    "error_gap = flux_err[i_gap]\n",
    "\n",
    "for i,_t in enumerate(T):\n",
    "    _n,_f = transit_noise(time_gap, flux_gap, error_gap,\n",
    "                        T_0=_t, width=WIDTH_HOURS,h_1=h1, h_2=h2)\n",
    "    Ngap[i] = _n\n",
    "    Fgap[i] = _f\n",
    "\n",
    "fig,ax=subplots(2,1,figsize=(6,4),sharex=True)\n",
    "ax[0].set_xlim(min(time_raw),max(time_raw))\n",
    "ax[0].plot(time_gap, flux_gap,'b.',ms=1)\n",
    "ax[0].set_ylabel(\"Flux \".format(median(C)))\n",
    "ax[0].set_ylim(0.9975,1.0035)\n",
    "ax[0].errorbar(median(T),1.003,xerr=WIDTH_HOURS/24,\n",
    "           capsize=5,color='b',ecolor='b')\n",
    "ax[1].plot(T,Ngap,'b.',ms=1)\n",
    "ax[1].set_ylabel(\"Noise [ppm] \".format(median(C)))\n",
    "ax[1].set_xlabel(\"BJD-{}\".format(int(min(BJD))));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&copy; Pierre Maxted, Keele University, 2019"
   ]
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
