{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exposure time integration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we'll briefly discuss how to account for finite exposure time integration, which tends to smooth out occultation light curves. Unfortunately, there's no analytic way (that I know of) to tackle this, so the thing to do is to oversample the light curve on a fine time grid and numerically integrate over the exposure window."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide_input"
    ]
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide_input"
    ]
   },
   "outputs": [],
   "source": [
    "%run notebook_setup.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import starry\n",
    "\n",
    "starry.config.lazy = False\n",
    "starry.config.quiet = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating a star-planet system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's instantiate a simple `Primary` object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "star = starry.Primary(starry.Map(ydeg=0, udeg=2, amp=1.0), m=1.0, r=1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's give it some limb darkening:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "star.map[1] = 0.40\n",
    "star.map[2] = 0.26"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what that looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "star.map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now create a featureless hot Jupiter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "planet = starry.kepler.Secondary(\n",
    "    starry.Map(2, amp=0), m=0, r=0.1, porb=1.0, ecc=0.3, w=30, t0=0,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have a star and a planet, we can instantiate the planetary system. By default, exposure time integration is **disabled**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys = starry.System(star, planet)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing a transit light curve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're ready to compute a transit light curve. Let's do this over 1000 cadences between $t=-0.1$ and $t=0.1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide_input"
    ]
   },
   "outputs": [],
   "source": [
    "# HACK: run this to pre-compile the flux method\n",
    "sys.flux(0.0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "time = np.linspace(-0.1, 0.1, 1000)\n",
    "flux = sys.flux(time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cool -- `starry` computed 1,000 cadences in just a few ms. Let's check it out:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time, flux)\n",
    "plt.xlabel(\"time [days]\")\n",
    "plt.ylabel(\"system flux\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, now let's enable exposure time integration. We have to instantiate a new ``System`` object for this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys_exp = starry.System(star, planet, texp=0.01, oversample=9, order=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We passed in the three keywords controlling exposure time integration. The first is ``texp``, the length of the exposure window in ``sys_exp.time_unit`` (usually days). The second is ``oversample``, the factor by which the light curve is oversampled. The larger this number, the more accurate the model will be, at the cost of extra computation time. Finally, ``order`` controls the order of the numerical integration: ``0`` for a centered Riemann sum\n",
    "(equivalent to the \"resampling\" procedure suggested by Kipping 2010), ``1`` for the trapezoid rule, or ``2`` for Simpson’s rule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide_input"
    ]
   },
   "outputs": [],
   "source": [
    "# HACK: run this to pre-compile the flux method\n",
    "sys_exp.flux(0.0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute the flux (and time the computation):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "flux_exp = sys_exp.flux(time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That was a factor of a few slower than the original evaluation, but it's not bad. Here's the comparison of the two light curves:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time, flux, label=r\"$t_{exp} = 0$\")\n",
    "plt.plot(time, flux_exp, label=r\"$t_{exp} = 0.01$\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"time [days]\")\n",
    "plt.ylabel(\"system flux\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing a phase curve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Exposure time integration also affects phase curves. Let's give the planet a random map and compute its phase curve with and without exposure time integration. We'll make the planet's rotational period really short so we can clearly tell the difference between the two:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "planet.map.amp = 0.1\n",
    "planet.prot = 0.05\n",
    "planet.map[1:, :] = 0.1 * np.random.randn(planet.map.Ny - 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "planet.map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's grab just the planet flux (by passing ``total=False`` to the ``flux`` method and keeping only the second vector):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "flux = sys.flux(time, total=False)[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "flux_exp = sys_exp.flux(time, total=False)[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the two light curves; it's clear that the finite exposure time has smoothed out the phase curve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time, flux, label=r\"$t_{exp} = 0$\")\n",
    "plt.plot(time, flux_exp, label=r\"$t_{exp} = 0.01$\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"time [days]\")\n",
    "plt.ylabel(\"planet flux\");"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "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
}