{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Asteroid Lightcurve and Its Rotation Period\n",
    "![Astropy: A Community Python Library for Astronomy](http://www.astropy.org/images/astropy_banner_96.png)\n",
    "========\n",
    "\n",
    "***\n",
    "by Chan-Kao Chang (2015 Big Data Summer School @ IANCU)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import astroML.time_series as ts\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as pl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Showing Lightcurve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read lightcurve data of asteroid (36035) from ./36035_test.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with open('./36035_test.txt','r') as ptf:\n",
    "  data = np.array([i.split () for i in ptf.readlines()])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Filter out bad detection\n",
    "sextractor flags = 0 (the last 2nd column) and ipac flags <= 2 (the last 1st column)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "idx = (data[:,-2] == '0') * (data[:,-1] <= '2')\n",
    "t = data[idx,0].astype(float)\n",
    "mag = data[idx,3].astype(float)\n",
    "magerr = data[idx,4].astype(float)\n",
    "magerr[magerr==0] = 0.01\n",
    "zp = data[idx,5].astype(float)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot lightcurve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.plot(t, mag, 'ro', label='no zp corrected')\n",
    "pl.plot(t, mag+zp, 'bo', label='zp corrected')\n",
    "pl.gca().invert_yaxis()\n",
    "pl.xlabel('MJD')\n",
    "pl.ylabel('R (mag)')\n",
    "pl.legend(loc=0, frameon=False, numpoints=1)\n",
    "pl.title('36035')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Zoom in on 2nd day"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.plot(t, mag, 'ro', label='no zp corrected')\n",
    "pl.plot(t, mag+zp, 'bo', label='zp corrected')\n",
    "pl.gca().invert_yaxis()\n",
    "pl.xlabel('MJD')\n",
    "pl.ylabel('R (mag)')\n",
    "pl.xlim(56708+1., 56708+1.5)\n",
    "pl.legend(loc=0, frameon=False, numpoints=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Periodogram (period analysis)\n",
    "use astroML.time_series.multiterm_periodogram(t, mag, magerr, omega, n_terms=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fres = np.arange(24/0.025)*.025 + 0.025\n",
    "pows = ts.multiterm_periodogram(t, mag+zp, magerr, fres*2*np.pi, n_terms=2)\n",
    "pl.plot(fres, pows)\n",
    "pl.xlabel('Frequency (rev/day)')\n",
    "pl.ylabel('Power')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pick out the best fit frequency and generate best_fit lightcurve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fre_pick = 3.025\n",
    "lc = ts.MultiTermFit(fre_pick*2*np.pi, 2)\n",
    "lc.fit(t, mag+zp, magerr)\n",
    "pre_t, pre_y = lc.predict(1000, adjust_offset=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot folded (phased) lightcurve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.plot(pre_t, pre_y, 'g-')\n",
    "phases = (t*fre_pick)%1\n",
    "pl.errorbar(phases, mag+zp, yerr=magerr, fmt='o', mec='none',capsize=0., elinewidth=0.5)\n",
    "pl.gca().invert_yaxis()\n",
    "pl.xlabel('Phase')\n",
    "pl.ylabel('R (mag)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def lc_ana(fn, fre_pick=0., unzp=True):\n",
    "    with open(fn) as ptf:\n",
    "        data = np.array([i.split () for i in ptf.readlines()])\n",
    "    idx = (data[:,-2] == '0') * (data[:,-1] <= '2')\n",
    "    t = data[idx,0].astype(float)\n",
    "    mag = data[idx,3].astype(float)\n",
    "    magerr = data[idx,4].astype(float)\n",
    "    magerr[magerr==0] = 0.01\n",
    "    zp = data[idx,5].astype(float)\n",
    "    fres = np.arange(24/0.025)*.025 + 0.025\n",
    "    pows = ts.multiterm_periodogram(t, mag+zp, magerr, fres*2*np.pi, n_terms=2)\n",
    "    if fre_pick == 0.:\n",
    "        pow_pick = pows[pows.argmax()]\n",
    "        fre_pick = fres[pows.argmax()]\n",
    "    else:\n",
    "        pow_pick = pows[abs(fres-fre_pick).argmin()]\n",
    "    lc = ts.MultiTermFit(fre_pick*2*np.pi, 2)\n",
    "    lc.fit(t, mag+zp, magerr)\n",
    "    pre_t, pre_y = lc.predict(1000, adjust_offset=False)\n",
    "    \n",
    "    pl.figure(figsize=(8,16))\n",
    "    pl.subplot(3,1,1)\n",
    "    if unzp:\n",
    "        pl.plot(t, mag, 'ro', label='no zp corrected')\n",
    "    pl.plot(t, mag+zp, 'bo', label='zp corrected')\n",
    "    pl.gca().invert_yaxis()\n",
    "    pl.xlabel('MJD')\n",
    "    pl.ylabel('R (mag)')\n",
    "    pl.legend(loc=0, frameon=False, numpoints=1)\n",
    "    pl.title(fn.replace('./','').replace('_test.txt',''))\n",
    "\n",
    "    pl.subplot(3,1,2)\n",
    "    pl.plot(fres, pows)\n",
    "    pl.scatter(fre_pick, pow_pick, s=30, c='r')\n",
    "    pl.xlabel('Frequency (rev/day)')\n",
    "    pl.ylabel('Power')\n",
    "    pl.xlim(0,25)\n",
    "    \n",
    "    pl.subplot(3,1,3)\n",
    "    phases = (t*fre_pick)%1\n",
    "    pl.plot(pre_t, pre_y, 'g-')\n",
    "    pl.errorbar(phases, mag+zp, yerr=magerr, fmt='o', mec='none',capsize=0., elinewidth=0.5)\n",
    "    pl.gca().invert_yaxis()\n",
    "    pl.xlabel('Phase')\n",
    "    pl.ylabel('R (mag)')\n",
    "    pl.title('F=%6.2f (rev/day)' % fre_pick)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lc_ana('./12755_test.txt')#, fre_pick=3.64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lc_ana('30567_test.txt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lc_ana('68855_test.txt')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
