lightcurve_fitting.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from lightcurve_fitting import lightcurve, models, fitting, bolometric"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Read the Light Curve\n",
    "Change the path to the file.\n",
    "It should have at least these columns: MJD, mag, dmag, filter.\n",
    "If the columns are not fixed width, you may need to add the keyword `format` (see [`astropy.table.Table.read`](http://docs.astropy.org/en/stable/io/unified.html#built-in-readers-writers)).\n",
    "Most reasonable filter names should be recognized.\n",
    "Also give the extinction ($E(B-V)$) and the distance modulus to calculate the absolute magnitudes.\n",
    "If necessary, you can give a second extinction correction for the host galaxy with the `host_ebv` keyword, which is applied at the host redshift."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "lc = lightcurve.LC.read('lightcurve_fitting/example/SN2016bkv.txt')\n",
    "lc.meta['dm'] = 30.79  # distance modulus\n",
    "lc.meta['ebv'] = 0.016  # Milky Way extinction\n",
    "lc.meta['host_ebv'] = 0.  # host galaxy extinction\n",
    "lc.meta['redshift'] = 0.002  # redshift\n",
    "lc.show_in_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make a bolometric light curve and/or color curves, skip to the [last section](#Bolometric-Light-Curve)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set Up the Parameters for the Fit\n",
    "If you only want to fit a subset of the data, you can do that here.\n",
    "Choose your model and initialize it with your light curve.\n",
    "(See the package documentation for information on the various choices of model.)\n",
    "You also need to give some kind of guesses for the parameters.\n",
    "Define your priors and starting guesses for each parameter.\n",
    "The walkers start in the boundaries defined by `p_lo` and `p_up` but can walk out to the prior boundaries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lc_early = lc.where(MJD_min=57468., MJD_max=57485.)\n",
    "model = models.ShockCooling2(lc_early)\n",
    "priors = [\n",
    "    models.UniformPrior(0., 100.),\n",
    "    models.UniformPrior(0., 100.),\n",
    "    models.UniformPrior(0., 100.),\n",
    "    models.UniformPrior(57468., 57468.7),\n",
    "]\n",
    "p_lo = [20., 2., 20., 57468.5]\n",
    "p_up = [50., 5., 50., 57468.7]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run the Fit\n",
    "You can modify the number of walkers and the number of steps here.\n",
    "I'm starting them with numbers that are probably **too small** so that you can test that everything works.\n",
    "You can save the results to a file using `save_sampler_as='filename.npy'`.\n",
    "\n",
    "When the fit is done, check the plots to make sure they have converged during the burn-in period."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampler = fitting.lightcurve_mcmc(lc_early, model, priors=priors, p_lo=p_lo, p_up=p_up,\n",
    "                                  nwalkers=10, nsteps=100, nsteps_burnin=100, show=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# View the Results\n",
    "This makes a corner plot with the posterior distributions and the $1\\sigma$ credible intervals, as well as a plot showing the best-fit models compared to the data in the upper right.\n",
    "You can save this plot with `save_plot_as='filename.pdf'`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig, ax_corner, ax_model = fitting.lightcurve_corner(lc_early, model, sampler.flatchain)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Check the Validity Times\n",
    "The shock cooling models are only valid for temperatures above 0.7 eV = 8120 K (Sapir & Waxman 2017), so you should check that you have not included observations where the model goes below that.\n",
    "If you have, you should rerun the fit without those points.\n",
    "If you used the Rabinak & Waxman option, the model fails even earlier, but you will have to check that manually."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_mean = sampler.flatchain.mean(axis=0)\n",
    "t_max = model.t_max(p_mean)\n",
    "print(t_max)\n",
    "if lc_early['MJD'].max() > t_max:\n",
    "    print('Warning: your model is not valid for all your observations')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bolometric Light Curve and Color Curves"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also make a bolometric light curve from the photometry table you defined in the first section.\n",
    "The light curve is divided into epochs (defined by the `bin` argument to `calculate_bolometric`), and processed four different ways:\n",
    "- Fitting the Planck function using `scipy.curve_fit`. This is very fast but may not give reliable uncertainties. The columns `temp`, `radius`, `dtemp`, and `dradius` come from this fit.\n",
    "  - The Stefan-Bolzmann law gives the total bolometric luminosity, `lum` and `dlum`.\n",
    "  - Integrating the Planck function between $U$ and $I$ band (observed) gives `L_opt`.\n",
    "- Fitting the Planck function using an MCMC routine. This is slower, depending on how many walkers (`nwalkers`) and steps (`burnin_steps` and `steps`) you use, but gives more robust uncertainties. The columns `temp_mcmc`, `radius_mcmc`, `dtemp0`, `dtemp1`, `dradius0`, and `dradius1` come from this fit. My convention for non-Gaussian uncertainties is that 0 is the lower uncertainty and 1 is the upper uncertainty.\n",
    "  - Integrating the Planck function between $U$ and $I$ band (observed) gives `L_mcmc`, `dL_mcmc0`, and `dL_mcmc1`.\n",
    "- Directly integrating the observed SED, assuming 0 flux outside of $U$ to $I$. Use this if you do not want to assume the SED is a blackbody. This yields the column `L_int`.\n",
    "\n",
    "The MCMC routine saves a corner plot for each fit in the folder you specify (`outpath`).\n",
    "I highly recommend looking through these to make sure the fits converged.\n",
    "If they didn't, try adjusting the number of burn-in steps (`burnin_steps`).\n",
    "To save the table, give `save_table_as='filename.table'` as an argument to `calculate_bolometric`.\n",
    "To save the plot, give `save_plot_as='filename.pdf'` as an argument to `plot_bolometric_results`.\n",
    "\n",
    "Beware of the units I'm using:\n",
    "- Temperatures are in kilokelvins (kK).\n",
    "- Radii are in thousands of solar radii ($1000R_\\odot$).\n",
    "- Luminosities are in watts (W). $1\\,\\mathrm{W} = 10^7\\,\\mathrm{erg}\\,\\mathrm{s}^{-1}$\n",
    "\n",
    "Optionally, you can calculate one or more colors at each epoch by giving the argument `colors` to `calculate_bolometric`).\n",
    "These get saved in the same output table in four columns per color, e.g., for $B-V$:\n",
    "- the color itself, `B-V`,\n",
    "- the uncertainty on the color, `d(B-V)`,\n",
    "- whether the color is a lower limit, `lolims(B-V)` (i.e., $B$ was an upper limit), and\n",
    "- whether the color is an upper limit, `uplims(B-V)` (i.e., $V$ was an upper limit).\n",
    "\n",
    "The color curves can be plotted conveniently with `plot_color_curves`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = bolometric.calculate_bolometric(lc, outpath='./SN2016bkv_bolometric', colors=['B-V', 'g-r', 'r-i'])\n",
    "t.show_in_notebook()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig1 = bolometric.plot_bolometric_results(t)\n",
    "fig2 = bolometric.plot_color_curves(t)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}