{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling Basics\n",
    "\n",
    "Here are a few examples of how to get started modeling transit light curves."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## `batman` for simulating transits\n",
    "\n",
    "The [batman](https://www.cfa.harvard.edu/~lkreidberg/batman/) package is a Python toolkit for generating model transit light curves for planets of given parameters. Laura Kreidberg wrote it to make transit-modeling with Python easier, and she did a great job documenting how it works.\n",
    "\n",
    "We wrote some quick wrapper so you can immediately start playing, by looking at transit models by themselves:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import henrietta as hsl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can make a plot of a model transit light curve for some set of parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hsl.example_transit_model(radius=0.1, a=4, b=0.9,  period=0.47, t0=0.05);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If I want to compare two or more different models, I can keep track of the plotting axes in which the first is plotted, and keep adding new curves to the same plot, on the same scale. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = hsl.example_transit_model(radius=0.1, a=4, b=0.9,  period=0.47, t0=0.05);\n",
    "ax = hsl.example_transit_model(radius=0.1, a=4, ax=ax);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to access the modeling tool directly, we can use the `BATMAN` function to make an array of model flux values. The inputs to this code are an array of times and the transit model parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np, matplotlib.pyplot as plt\n",
    "\n",
    "# make some times that span 0 to 5 days, in 2-minute increments\n",
    "times = np.arange(0, 1, 2.0/60.0/24.0)\n",
    "\n",
    "# use the BATMAN function to make a transit model\n",
    "modelflux = hsl.BATMAN(times, period=4.2, t0=0.3, radius=0.11, a=9.0, b=0.2)\n",
    "\n",
    "# plot the times and model fluxes\n",
    "plt.scatter(times, modelflux, s=3)\n",
    "plt.xlabel('Time (days)')\n",
    "plt.ylabel('Relative Flux');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All of the other visualization tools described on this page are calling this `BATMAN` function to make their models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## making a simulated transit dataset \n",
    "Sometimes it can be really helpful to be able to make a simulated noisy dataset, starting from some model parameters. The `simulate_transit_data` function can do that for you!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lc = hsl.simulate_transit_data(period=1.23, t0=0.1, radius=0.1, a=10, b=0.5)\n",
    "ax = lc.scatter()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## visually comparing a model to data\n",
    "\n",
    "We want to be able to compare our transit model to our data. Here's a little wrapper called `plot_with_transit_model`, which can plot a light curve, and then overplot a transit model on top. As inputs, this function takes a light curve, and the transit parameters `period`, `t0`, `radius`, `a`, `b` as keyword arguments. The top panel shows the light curve data; the bottom panel shows the \"residuals\", which are defined as the data minus the model. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lc = hsl.simulate_transit_data(period=1.23, t0=0.1, radius=0.1, a=10, b=0.5)\n",
    "hsl.plot_with_transit_model(lc, period=1.23, t0=0.1, radius=0.1, a=10, b=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we plot the wrong model for our data, our residuals show big bumps in them!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lc = hsl.simulate_transit_data(period=1.23, t0=0.1, radius=0.1, a=10, b=0.5)\n",
    "hsl.plot_with_transit_model(lc, period=1.23, t0=0.1, radius=0.1, a=10, b=0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we feed this function a folded light curve, it will plot everything in terms of phase instead of time. The x-axis will span only one orbital period, with data repeating back on it many times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lc = hsl.simulate_transit_data(period=1.23, t0=0.1, radius=0.1, a=10, b=0.5)\n",
    "folded = lc.fold(period=1.23, phase=0.1/1.23)\n",
    "hsl.plot_with_transit_model(folded, period=1.23, t0=0.1, radius=0.1, a=10, b=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## keeping times straight\n",
    "\n",
    "One tricky thing might be how to match up a mid-transit time you take off the NASA Exoplanet Archive (in BJD) with the times in your light curves. Remember that most Kepler light curves store their time as BKJD (Barycentric Kepler Julian Date), which is offset from actual BJD (Barycentric Julian Date) by 2454833 days. We've added a couple of small tools to help with this conversion:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print the original t0, in BJD\n",
    "t0 = 2455678.1234\n",
    "print(t0, \"is the BJD.\")\n",
    "\n",
    "# convert to BKJD (to match Kepler light curves)\n",
    "t0_bkjd = hsl.bjd2bkjd(t0)\n",
    "print(t0_bkjd, \"is the BKJD.\")\n",
    "\n",
    "# convert back again\n",
    "t0_bjd = hsl.bkjd2bjd(t0_bkjd)\n",
    "print(t0_bjd, \"is the BJD, after converting back and forth to BKJD.\")"
   ]
  }
 ],
 "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.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}