{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following features of impulse response simulator have been implemented in this notebook.\n",
    "\n",
    "1- Find lag-frequency spectrum of a simple delta impulse response.\n",
    "\n",
    "2- Find lag-frequency spectrum of a _more_ realistic impulse response based on real physical principles.\n",
    "\n",
    "3- Compute lag-frequency spectrum of delta impulse responses with same intensities and varying positions at different energy levels.\n",
    "\n",
    "4- Compute lag-frequency spectrum of delta impulse responses with same positions and varying intensities at different energy levels."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import libraries and obtain data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from stingray import Crossspectrum, Lightcurve, sampledata\n",
    "import numpy as np\n",
    "from scipy import signal\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define variability signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lc = sampledata.sample_data()\n",
    "s = lc.counts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lag-frequency Spectrum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simple Delta Impulse Response"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a delta impulse response with a delay of 10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "delay = int(10/lc.dt)\n",
    "h_zeros = np.zeros(delay)\n",
    "h = np.append(h_zeros, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find output signal by taking convolution of variability signal and impulse response."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "output = signal.fftconvolve(s, h)\n",
    "# To make two counts of equal size, remove last 'delay' entries and avoid first zeros\n",
    "output = output[delay:-delay]\n",
    "s_mod = s[delay:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize input and output signals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEACAYAAABbMHZzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXd4VEX3xz8hgIDSS2ihSG+CIEVECChFELGACBZ81deC\nInbFnwXhlSqCgGJDmgJSlV4EA9JbIEDo0hJa6BBK2vn9Mbtkk+wmu5tNdjc5n+fZZ++dO3fu2U32\ne+eeOXMGFEVRFEVRFEVRFEVRFEVRFEVRFEVRFEVRFEVRFEVRFD/nCBAOhAGbLGXFgOXAfmAZUMSm\nfj/gALAXaGdT3gjYaTn2TaZarCiKomSYwxixt2UY8IFl+0NgiGW7NrAdyANUAg4CAZZjm4Amlu1F\nQIfMMVdRFEXxBIeB4inK9gJBlu3Sln0wvf4PbeotAZoBZYA9NuVPAd973FJFURQlXXI5WU+Av4At\nwH8tZUHAacv2aZJuBGWBSJtzI4FydsqjLOWKoihKFpPbyXr3ASeBkhg//94Ux8XyUhRFUfwAZ8X/\npOU9GpiL8dufxrh7TmFcOmcsdaKAYJtzy2N6/FGWbdvyqJQXqlKlihw6dMhJsxRFURTgEFDVlROc\ncfsUAApatm/HRO/sBOYBvSzlvYA/LNvzMP78vEBloBpmoPcUcBloihkAftbmnKRPcOgQIuLTr88/\n/9zrNqidaqfaqTZaX0AVV4QfnOv5B2F6+9b6v2FCO7cAM4AXMaGgT1rqRFjKI4B4oDdJLqHewEQg\nPybaZ4mrBiuKoigZxxnxPww0sFN+HnjQwTmDLK+UbAXqOWeaoiiKklk4G+2j2BASEuJtE5xC7fQs\naqdn8Qc7/cFGdwlIv0qWIxYflqIoiuIEAQEB4KKea89fURQlB6LiryiKkgNR8VcURcmBqPgriqLk\nQFT8FUVRciAq/oqiKDkQFX9FUZQciIq/oihKDkTFX1EUJQei4q8oipIDUfFXFEXJgaj4K4qi5EBU\n/BVFUXIgKv6Koig5EBV/RVGUHIiKv6IoSg5ExV9RFCUH4qz4BwJhwHzLfn8g0lIWBjxkU7cfcADY\nC7SzKW8E7LQc+8ZtixVFUZQM46z49wUiAOv6igJ8DdxteS22lNcGulveOwDfkbS02DjgRaCa5dUh\ng7YriuINfvwRtm/3thVKBnFG/MsDHYGfSRLyAOyvF9kFmAbEAUeAg0BToAxQENhkqTcZeNRdoxVF\n8RJHjkCfPjB/frpV06R9ezh71iMmKe7hjPiPBN4HEm3KBOgD7ADGA0Us5WUx7iArkUA5O+VRlnJF\nUXyN1atBxP6xTz+FKlVg71732z97FpYtg3373G9DyTDpif/DwBmMX9+2pz8OqAw0AE4CIzLFOkVR\nspZz56BVKxhh5ycdFgZ//QVjxmRM/MPCzPuRI+63oWSY3Okcbw48gnH75AMKYVw2z9nU+ZmkgeAo\nINjmWHlMjz/Ksm1bHuXoov3797+1HRISQkhISDpmKoriESIioFo1+OoruOcesP3tffABfPYZNGoE\n+/ebp4MAe97fdAgLM+ep+LtNaGgooaGhGWrDlb9cK+A9oDPGh3/SUv420BjoiRnonQo0wbh1/gKq\nYtxEG4E3MX7/hcBoYImd64g4euRUFCVz+f572LIFuneHXr1g82YoV864afr0gV27IE8eKF3a1Ctf\nPv02U9KjBxw7BrVrw08/ef4z5EACzE3YpTuxK3H+ASRF+wwDwjE+/1aYGwCYiKAZlvfFQG+bc3pj\nnhIOYAaC7Qm/oijeJCLCiHLbtvDGG9CtG9y4YXr9gwcb4QeoWdN9n31YGDz2GBw96jm7FZdx45kt\n09Gev6J4iwcfhPfegw4dIDHRiPTx45AvH6xdm+TmefVVqFcPXn/dtfavXoVSpWDjRnjiCeM+UjJM\nZvf8FUXJ7kREQJ06ZjtXLpg0CW67zYwB2Pr3a9Rwb9B3xw7TftWqxvWTmJj+OUqmoOKvKIrhwgXT\nM7f14xcpAuvXQ/Pmyeu66/YJC4O774b8+U3bp05lzGbFbVT8FUUxRERArVrORfCk1fPfvdu4juwR\nFgYNG5rtSpU04seLqPgrimKwDvY6Q8WKEB0NMTGpj82bB6NHw8WLqY9Ze/6g4u9lVPwVRTG4Iv6B\ngWY+gL0B29Wr4Y47zE3AlthY87RQr57ZV/H3Kir+iqIYdu9OGux1Bnuun/h4WLcOPv8cZs5Mfiwi\nAipXhgIFzL6Kv1dR8VcUxeBKzx/sD/pu3w7BwfD887BqFVy6lHRs27Yklw+o+HsZFX9FUYxIX7wI\nFSo4f469nv+qVSY3UOHCJjWEbfZP28FeUPH3Mir+iqLAnj2mJ5/LBUmw1/NfvRpatjTbXbvCrFlJ\nx2wHe8HcaDTW32uo+CuKYvz9rrh8AKpXNwO+VvFOTIR//kkS/0cegZUr4fJlcyw8HBo0SDq/QAHz\nhHD6tGc+g+ISKv6KoiSf2esshQqZiVrHj5v9nTuhRAkoU8bsFykC998PCxfCwYNQvDgULZq8DXX9\neA0Vf0VRXB/stWLr+lm92vj7bena1UT9pHT5WFHx9xoq/oqiuC/+toO+q1YluXysdOliFoBZvVrF\n38dQ8VeUnM6VK2a2bqVKrp9r7fmL2O/5Fytm8gL98kvySB8rKv5eQ8VfUXI61kifwEDXz7X2/Pfu\nhdtvtx8qal0TQHv+PoWKv6LkdNx1+YC5aezdmzzEMyWPPgoPP5w0EGyLir/XUPFXlJxORsQ/ONik\ngl6wILXLx0rx4mayl71soRUraqy/l1DxV5Scjjsx/lZy5TLx/osWOe75p0WBAiZkVGP9sxwVf0XJ\n6WSk5w/G9VO6NFSp4t756vrxCs6KfyAQBlgTdRQDlgP7gWVAEZu6/TCLtO8F2tmUNwJ2Wo59477J\nipINOH3aLI14+LB37YiJMatp3Xmn+23UrGl6/c4sAmOPihVV/L2As+LfF4gArCurf4QR/+rACss+\nQG2gu+W9A/AdSYsKjwNeBKpZXh0yaLui+B+bN8NzzxnBHDUKVqzwni03b8IHH0CjRpA7t/vtvPEG\nDB/u/vmVKsHRo+6fr7iFM+JfHugI/EySkD8CTLJsTwIetWx3AaYBccAR4CDQFCgDFAQ2WepNtjlH\nUbI/Z89CmzYm7LFePTh0yNwEvOXr/vdfuO8+iIpKnnnTHUqUSL7ur6uo28crOCP+I4H3Advh+CDA\n+l972rIPUBaItKkXCZSzUx5lKVeU7M+xY9CiBTRrZkT//ffN5KegIO+I/6xZxpZnn4W5c1Pn28lq\nVPy9QnrPeg8DZzD+/hAHdYQkd5BH6N+//63tkJAQQkIcXVpRfJw9e6B9e3jnHXjrreTHgoLMqldZ\nyYIF8O67Jtla48ZZe21HqPi7TGhoKKGhoRlqIz3xb45x8XQE8gGFgCmY3n5p4BTGpXPGUj8KCLY5\nvzymxx9l2bYtj3J0UVvxVxS/ZeNGk9tm+HDTy06JN3r+mzdDr16+I/xgBnyPHjUpItwdNM5hpOwU\nf/HFFy63kZ7b52OMmFcGngJWAs8C84Beljq9gD8s2/Ms9fJazqmG8fOfAi5j/P8Bljas5yhK9uPY\nMTOrdfx4+8IPUKpU1ov/oUPuh2RmFrffDgULaqx/FuNqnL/VvTMEaIsJ9Wxj2QcTETTD8r4Y6G1z\nTm/MoPEBzEDwEretVhRf59NP4ZVXoFMnx3W80fP3RfEHdf14AV98xhIRjw4hKErWsmOH8fPv329m\nrzoiMRFuu83E2ufNmzW2lSpl7LOXZ8ebPPkkPPYY9OjhbUv8kgDjLnNJz3WGr6J4mg8/hP/7v7SF\nH0xqhJIl4cyZtOt5iitXzI2mdOmsuZ4raKx/lqPiryie5K+/zJKFr7ziXP2sdP0cOmRm8vrioKq6\nfbIcFX9F8RSJiWbG7KBBzrtxslr8fdHfDyr+XkDFX1E8xfTpJk1Ct27On6Pib6hUyft5jnIYGUjo\noSjKLW7eNH7+iRNdc6tktfjXr58113KVSpWS8vrn0j5pVqDfsqJ4grlzTa/a0YImjvC0+H//vRFR\ne/hyz1/z+mc5Kv6K4gl+/x2eecb18zwp/gkJ8PHHJnWDPQ4ehKpVPXOtzED9/lmKir+iZJTLl2Hl\nSrNWrat4Uvy3bDFLKm7enPpYbCycPGl/gXVfQcU/S1HxV5SM8uefxt1TpEj6dVPiSfFfutTYsWVL\n6mNHjpi0y3nyeOZamYGKf5ai4q8oGeX336F7d/fO9aT4L1sG771n3DsxMcmP+bK/34pG/GQpKv6K\nkhEuXIB//oFHHnHv/BIl4OJFiI/PmB2XLkF4ODzwANSpA2FhyY/7g/hXrqw9/yxExV9RMsLcufDg\ngyYrpTsEBpqFXaKjM2bHypXQvDnkz2/SNaf0+/uD+KvbJ0tR8VeUjJARl48VT7h+li6Fdu3Mtr+K\nf8WKSbH+Sqaj4q8o7hIdbRZsSSttszNkVPxFjPi3b2/2/VX88+c3g+anTnnbkhyBir+iuMucOdCh\ng1mMJCNkVPwPHoS4OKhd2+zXqmUE9MIFs5+YaAZS77wzY3ZmBer6yTJU/BXFXTzh8oGMi7/V5WNN\nKxEYCHffDVu3mv2TJ83s2TvuyLitmU3lyhrxk0Wo+CuKO5w6ZSJqHnoo421lVPyXLUvy91u5554k\n148/uHysaM8/y1DxVxR3WLDACH++fBlvKyPiHxsLq1ZB27bJy239/ir+ih1U/BXvkphoctL4Gzt2\nQJMmnmkrI+K/fj3UqAHFiycvtxV/X8/pY4uKf5aRnvjnAzYC2zGLsg+2lPcHIoEwy8v22bcfZpH2\nvYDts2gjYKfl2DcZtFvJLowcaWal+hu7dkHdup5pKyPibxviaUuVKmaW76lT2vNX7JJePv8bQGvg\nmqXuGqAFIMDXlpcttYHulvdywF9ANUv9ccCLwCZgEdABWOKJD6H4McuWpU5F4OuIwM6dviH+y5aZ\nG2hKAgKS/P7+JP4VK8Lx4+ZpMDDQ29Zka5xx+1yzvOcFAgFL/JjdleK7ANOAOOAIcBBoCpQBCmKE\nH2Ay4EYKRCVbERdn3Bbh4f41sce64HpQkGfaK1kSzp1z3f116RLs2wdNm9o/3rixSfLmT+KfL5+Z\n8XzypLctyfY4I/65MG6f08DfwG5LeR9gBzAesKYzLItxB1mJxDwBpCyPspQrOZktW4wvulAhOHrU\n29Y4j9Xl46mF0PPkgcKFzQ3AFTZuhEaNHK8X3LixeTKIjzc5hPwFdf1kCc4s45gINAAKA0uBEIwL\nZ4Dl+EBgBMal4xH69+9/azskJISQkBBPNa34EqGhEBJieq87dpgYb3/Ak/5+K1bXT6lSzp+zdq3J\n5+OIxo1hwwZo0MBzN6qswCr+LVp42xKfJTQ0lNDQ0Ay14coavpeAhcA9gO1VfwbmW7ajgGCbY+Ux\nPf4oy7ZteZSjC9mKv5KNCQ2F3r3httuM68edxVC8wa5dpsftSaziX6+e8+esWwdvveX4eNmyULq0\n/7h8rGjPP11Sdoq/+OILl9tIz+1TgiSXTn6gLSa6p7RNnccwUTwA84CnMOMDlTGDvZuAU8BljP8/\nAHgW+MNla5Xsg9Xff//9ZlHxHTu8bZHzZGbP31ni443b5957HdcJCDC9fxV/xQ7p9fzLAJMwN4lc\nwBRgBWbAtgEmiucw8IqlfgQww/IeD/S21MGyPRFzE1mERvrkbLZuNaJUrBjcdRd8+ql37dm/H6pV\nS989IgK7d5uc+Z7EVfHfudOszFWsWNr1+vZ1b4Uxb1KpkkmdoWQq6Yn/TqChnfLn0jhnkOWVkq2A\nC8+0SrbG6u8HqF4doqLg6lXv5J9JSDA9+SlT0s/Vc+yYGaAuWtSzNrgq/uvWpe3vt/LAA+7b5C10\nUZcsQWf4Kt7BVvxz5zaZKHfuTOuMzCM6GnLlMv7z9CJudu3yfK8fXBf/tWvhvvs8b4cvUKECREb6\n58xvP0LFX8l64uJMz/X++5PK6tc3g77e4MQJc/Pp1g3efTftupnh74fM6/n7I7fdZkJTT5zwtiXZ\nGhV/JevZutXklrf1V991l/ODvp5a8NzKiRMmMubLL+Hvv2H5csd1fUH8rS6y6tU9b4evoIO+mY6K\nv5L1rFqV5PKx4mzPX8SI77//es4eq/gXLAjffw+vvOI45YQviL+11+9PsfuuouKf6aj4K5nLzJnG\nf2uLrb/fyl13OZfmISoKzp41KQs8hVX8waRpbt4cPvssdb34eDMhzbpilicpVcqMPTiT5iI7+/ut\n3HmnyUaqZBoq/krmsXs3vPSSiTW3ulLi4ox42fr7waQkdibNg/XpwJPpIGzFH0yitN9+g+3bk9c7\ndAjKlMn4so32uO0206516cW0SG9mb3agTh3zlKVkGir+SubxySemBz11KvTqBQMHmnw+lSunzj8P\nzvn9w8NNdJAnXQIpxb9kSfj449S9/8xy+VhxxvUTEwMRESZjZ3amXj3vRX/lEFT8lcxhwwYj9L17\nQ+vWZnv5cujSJbXLx4ozfv+dO6Fly8wVf4CXXzbLNFoXRAHzJJPZ4m/NGOqIzZuNMObPn3l2+ALV\nq5vUzteupV9XcQsVf8XziMBHH0H//kkiVbYsrFgBffrAs8/aP8/Znn/nzpnr9gGTWjhl798Xev7r\n1mV/fz+YTKfVq5unHCVTUPFXPM/SpUbEevVKXp4nj0nj4MhlkV7P/+ZNMwj40EOe6/nHxcH58/az\nab74IuzZYwQXfEP8c4K/34q6fjIVFX/FsyQmQr9+JmY+tytJY0me5sEee/ea8YIqVYx7JDbWuXaX\nLIHXXrN/7PRp4+O3t2pU3rxJ4xY3b8Lhw2a93MyibFmTPsIRiYkmGV5OEf+77lLxz0RU/BXPMmOG\n6eE/9pjr51rTPDiK8ggPN4KQO7eJukkZQuqINWuMaNrDnsvHll69jOj/+KMJP3S0cIonuPde07N3\nxO7dJqdQmTKZZ4MvoT3/TEXFX/Ec8fGmpzxkiPsTkNLy++/caY6DWevVWb9/eLhxF4mkPpae+OfJ\nY3r+H36YuS4fgGbNzGd09OSzYoV/JmpzFxX/TEXFX/Ec8+cbF0qbNu63kZbf39rzB9dmgO7YATdu\n2Penpyf+AE8/DcHBmZPQzZZ8+cx4iKPef04T/3LljLstvQgoxS1U/BXPMXasiebJCA0bOnbR2Ip/\nxYrOif/Fi2ZAt1Ej+zNGnRH/3Lnhzz8djxt4klatzAzolMTHw+rVGbux+hsBAdr7z0RU/BXPEBFh\nXl27Zqyd5s3NoG9KoT571sR8B1tWCa1UyTm3z86dxl1Tvbr9lBDOiD9AzZrmqSazCQmxL/6bN5vP\nnBU2+BIq/pmGir/iGb791kyMyuiAaO7c5gYyY0by8p07jRBYxxKcdfvs2GGeFqpUcb/nn5U48vvn\nNJePFRX/TEPFX8k4ly7BtGkmG6Yn6N4dpk9PXmbr8gHnB3yt51Wt6h/inz+/cVFZ5xZYUfFXPIyK\nv5JxJk+Gtm09J6ItWpgVtWxnd6YU/+BgI9zx8Wm3FR5uBpH9Rfwhtevn2jXj9mnZ0lsWeY+6dc3/\ngTPZThWXSE/88wEbge2YRdkHW8qLAcuB/cAywHaF6H7AAWAv0M6mvBFmTeADwDcZNVzxERITzUDv\nG294rs1cuUzv33YR7/Bw0wu0kjev8X+ntdpTQoKZM1Cvnn3xv3kTrlyxn2TOm6Qc9F271tzAChb0\nmkleo3Bhs6qXJ9dvUID0xf8G0BpoANxl2W4BfIQR/+rACss+QG2gu+W9A/AdYA34Hge8CFSzvDp4\n6kMoXmTFChOi2KKFZ9u1un5EjIhHRKSOs0/P7//vv0Y4Chc2Ai9iIn+snDwJpUubm40v0ayZudlZ\n/f451eVjRV0/mYIz//XWtHp5gUDgAvAIMMlSPgl41LLdBZgGxAFHgINAU6AMUBDYZKk32eYcxR0u\nXjR5abKalBOlxowx4Z2eXlWqSROTvmHHDhOlExRk8v3bkl64p9XlA8a+lL1/X3T5ABQoYEJerX5/\nFX8V/0zAGfHPhXH7nAb+BnYDQZZ9LO9Blu2ygO2c+0ignJ3yKEu54i5PPJHcLZLZxMXBf/9ronGK\nFTOpDqwC1bOn568XEABPPWV6/yn9/VbSC/e0RvpYSRnx46viD0l+/wsXTE6jZs28bZH3UPHPFJzJ\nvJWIcfsUBpZiXD+2iOXlMfr3739rOyQkhBBH+d9zKmfPGmHIqtS+MTHw5JPGv3/unHHDXLxohKlE\nCdNTzQy6dzc5gvLksS/+FSvCxo2Ozw8PN7NzrfhLzx+M+H/6qVkFrXlzs9JXTqVePfjiC29b4VOE\nhoYSam8+iAu4knbxErAQM3B7GigNnMK4dKzzr6OAYJtzymN6/FGWbdvyKEcXshV/xQ7z5xsx8OQ6\nto6IjoaHHzYJ1376yQgxZM0gaf365nNOngxffZX6eKVKqecD2BIeDsOGJe1XrQp//52078vi36yZ\neXKZNy9nu3zAZFI9ehSuX8/+i9g4ScpO8Rdu3BzTc/uUICmSJz/QFggD5gHWZO29gD8s2/OApzDj\nA5UxA7ubMDeJyxj/fwDwrM05iqvMnQv/+U/mi//hw+bp4sEHYcKEJOHPKqyun2PHHLt9HPn8L182\nuXyqVEkqq1o1+Xd24oTvZsgsUADuvht+/VXFP08eqFZNF3bxMOmJfxlgJcbnvxGYj4nuGYK5EewH\n2lj2wYSDzrC8LwZ6k+QS6g38jAn1PAgs8dSHyFFcvWpcPn372o9b9xQi8MIL5ibz5ZeeH9B1lqee\nMiGOVaumPlahgknrbC8GfOdOk4jNNk+/P7l9wLh+ChaEBg28bYn30dz+Hic9t89OoKGd8vPAgw7O\nGWR5pWQrUM9OueIKS5aYvO/Vqpk49UuXTCijp1m50uTYef99z7ftCjVrmicQe4ut5Mtn8tufPGky\nQNpiG+ljpXRpM35x+bKJHPJ18X/8cXPTtffZcxo66OtxfCzAWUmXOXPMIGhAgIm4yQzXj4jJy//F\nF66vxpUZpDW+4Cjc016EUMrvzNfF/+67YcAAb1vhGzRpYrKaKh5Dxd+fiI2FxYuhSxezX6VK5oj/\nokXGvdS9u+fb9jSOwj1Thnlasbp+YmJMjv+iRTPdRMUD3H+/uck7u4CPki4q/v7EypVQu3bSIGVa\n4r9zpxFvR0siOiIx0YQYDhjgezNf7WGv55+YmHzVL1usg74nT5pev7fGMhTXyJ0bHn0UZs/2tiXZ\nBj/4dSu3mDs3+dq4aYn/smWmh9umjZmcdfKk89cICDA/NH/AXs//yBHTo7fXq7f2/H3d5aOkpmtX\nmDXL21ZkG1T8/YWEBLOalLPiv327Sba2f78Rwbp1Yfjw9K/x2Wfwv//5T4/YXrinI5cPqPj7M23a\nwL59JsJLyTAq/v7Chg1QqlTyuPW0xD8szIQIFiliJjpt3Wqyb65Z4/ga06eb+h38KOdeyrz+cXHw\n229msNQeKv7+S5488MgjJuhByTAq/v5CSpcPmDj3U6dMyKct16+bjJa2C45XqmR69O+9lzo5G5jw\nx48/9m5MvztYxV/EpJzo2NHkv//gA/v1y5UzKSoOHVLx90eeeEJdPx5Cxd8fELEv/rlzm0VNUro9\ndu0ya9amXFLx6adNxNDMmamv8d570K6dmVjkT9x+O9xxh8nxc999Zl7AvHmOc98HBpob4Zo1Kv7+\nSNu2ZjD/1ClvW+L3qPj7Azt3mgiWlJOWwL7rZ/t2+7NCc+UyOXL69Uv+tLBsGSxdCiNGeNburKJS\nJZOC4tVXTYrp9OYmVK1qvlMVf//jttugUyfTGVIyhIq/P/DHHyb6xp47xt7C5Nu3O/Z5t2ljesfj\nxpn9S5fgpZfg559T58v3F1580SR469PHufpVq5qnKRV//0SjfjyCir8/YM/lYyVlsjJIGux1xLBh\nMGiQScn87rvGT962refszWpefdV8Bmex5glS8fdP2rc3AQzR0d62xK9R8fd1jhwxOXYc5e5P6fZJ\nSDAujbTEv04dM0v48cfNKlHphYBmN6pWNWMFOXFN3OxA/vwmIu0PTQycEVT8fZ0//oDOnR0n90op\n/gcPmoXN00v2NmCAiYf/+eecJ4L16kHr1v4V1aQkp2vXtNdyUNJFxd/X+eMPxy4fMInKjhxJSmvs\naLA3JWXKmIiJnJgrvmxZsyCO4r906gR79qS9kpuSJir+vkx0tBHzBx1lz8Ys+lG0qHENQdqDvSlJ\nGQqqKP5C/vzm6fX99+3PW/FnsujzqPj7MgsWmIHYfPnSrmfr+klvsFdRsgu9esH58+Z3kp146CEz\n294RsbEm624GUfH3ZebOdS7Bmq34O+v2URR/JzAQhg6FDz+E+HhvW+MZzpyBf/6Bd94xYdj26N/f\nI4ssqfh7grVrTW54T2JdrrFTp/TrWsX/5EnzIyhf3rO2KIqv0rGjyXk1caK3LfEMCxeaz9SpE3z+\neerja9aY9bT798/wpVT8M0p0tEmJUL06/PKL53ogS5dCs2Ym0Vp6WMXf2uvXKBYlpxAQYOatfP65\nWaDH35k/30T3DR4M06aZFemsXL4Mzz0HP/4IQUEZvpQz4h8M/A3sBnYBb1rK+wORQJjl9ZDNOf0w\nC7XvBdrZlDfCrAt8APgmA3b7DsuXm7v077/DpEkmlfDcuRkftEkvyscW6yzfsDDnB3sVJbvQpAm0\naAGjRnnbkoxx44aZd9OxI5QoYQa0X389SUv69jXBH507Z5lJpQGrE/kOYB9QC/gceMdO/drAdiAP\nUAk4CFi7opuAJpbtRYC93MHiVzz7rMh335ntxESRRYtEKlUSWbDA/TZjYkSKFhWJjHSufnS0SOHC\nIl27ivz6q/vXVRR/5eBBkeLFRc6c8bYl7rN4sch99yXtx8eL3HOPyOTJIrNmiVStKnLlit1TAZd7\nm870/E9ZxBzgKrAHKGfZt+df6AJMA+KAIxjxbwqUAQpibgAAkwE/WS7KAYmJJila+/ZmPyDAjNS/\n/z5Mnep+u8OGmSifcuXSrwtmgXMRWLVKB3uVnEmVKtCzp0lJ7qts3w5vvmnm1sTGpj5udflYCQyE\n774z6cl5+hzTAAAgAElEQVR794YpU0wGWw/hqs+/EnA3sMGy3wfYAYwHrM7pshh3kJVIzM0iZXkU\nSTeRrEEE7rknKSY+o4SHm2Rod96ZvLxrVzNwc+2a620eO2YyU7qSciEgwPzzX7kCNWq4fk1FyQ58\n8okRyMOHvW1JEnFxZhGlhg1NSpVixUz5Dz8kryeSWvwBGjeGZ54x0T/NmnnUtHRy3ybjDmAW0Bfz\nBDAOGGA5NhAYAbzoCaP624xkh4SEEOIox3xiovnSHKU+SMmOHSYh1MaNJq9NRlmyJKnXb0upUuaP\ntmiRuRG4wgcfmOyUFSq4dl7VqiaVcXrpjBUlu1KqlOlZf/op/Pqrt60xTJtmUqgMH24y6gYGGh1q\n397MU7Bm0g0PNyuV1aqVug07HcHQ0FBCQ0Mz13YLeYClwFsOjlfCDOQCfGR5WVmCcfuUxriMrPQA\nvrfTlvM+sqFDRR5/3Pn6AweK3HabyP/9n/PnpEVIiGPf/s8/izzxhGvtrV4tEhxsfP6u0q+fyH//\n6/p5ipKduHxZpHRpkW3bMtzUpRuX5M5v7pQbcTfcb+TVV0VGjUpd/uyzIp99lrQ/cKBI375uXwY3\nfP7OEIDxz49MUV7GZvttwOrktg745gUqA4dIGhvYiLkRBOCJAd9GjUTy5xdZv965+s2aibzxhkjH\nju59w7Zcvixyxx0iV6/aP37unEihQqaeM8THi9x9t8i0ae7Zc+qU8wPEipKdGTtWpH371OWJiebl\nJH8f/lvoj6w6ssp9Wxo0ENmwIXX54cMixYqJnDxp9ps2FVm+XBITE2XgqoFy6Pwhly5DJg343gc8\nA7QmeVjnUCAc4/NvZbkBAEQAMyzvi4HeNob1Bn7GhHoexDwVuEdkpPHtff21WXs2vdDK6GiIiDDh\nUtu2uX3ZW/z9NzRtalID26NYMbj/frOkoDNMmGDa6t7dPXuCgpwfIFaU7Mx//2tCn1esMPtXrhid\nCA52abW6bSe3ERgQyIp/VyQVLl5s5uCcOZN+AzExsH+//SCMSpVMzP7AgXD6NOzbBy1b0m9FP4as\nGcLQNUOdtjM74dyt7rvvRJ55RiQuTqRaNZFly9KuP2mScRElJpqQsBMnXLqzpuK110SGDUu7zpQp\nIg8/nH5bFy+aR9WtWzNmk6IohunTRRo2FPn0U5ESJUS6dzcaUKGC0Qwn6Dm7p/SY1UNa/NLCFNy8\naZ72Q0JEihQRKVdOpHNnE2Zqj9BQ421wRHS00aJ+/US6dZNR60dJjTE1ZPeZ3VJkSBE5G3PW6Y9L\nJvX8fZN58+CRR8wA58CB6ff+Fywwk7ECAsxEqLAw968tYgZ7O9jzWtnwyCOwerVZMSstvv3WTN5o\n2NB9mxRFSaJbNyhd2vSq1683idKee86kPnEynfe2k9t4s+mbbD+1nauxV02gSM2a5qn//HmTg6dk\nSfjpJ/sNbNiQdoROiRLw9tsweDDTW5dk+LrhLH1mKbVL1uaxmo/xw9YfHJ+bTUn/Nnf5skjBgkn+\n9IQE41ubNct+/dhYc6e2+tc++MAMsLjL/v0iZcs65z98/HGR8eMdH792TSQoSGTXLvftURTFOaZO\nFWnTJt1ql29clgJfFpDY+FhpNaGVLNq/SOTzz4122LJ+vUi9evYbefRR8wSSFlevyvIn75GSQ0vI\njlM7bhVvP7ldyo4oKzfjb6Zrq0h26vnXrAk9ehj/3PXrqY8vXQrNmyetQJUrl1mT9pNPzDKGKVmz\nxoRCli5t9jPa81+61IRqOZNDp3t3k/rBERMmmLGDOnXct0dRFOd44gkz9hcRkWa1Had3ULdUXfIE\n5uGByg+w8vBKWLnShGva0rixSah4/HjycpH0e/5AYoH89LjnCDOenMldQXfdKq9fuj41itdgVkTm\nLVTvm+I/Y4Zxqcyfb1w6KZk3z0yYsKVDB/MYNWVK6voLF8LDDyftN2yYsUFfR/H99ujUyTwu2lts\nOj7exPB+9FHqY4qieJ68eeHll42rNQ22ndxGw9LGDdumchtWHFpuNKNFi+QVAwONFixenLz82DHz\nns58nb1n91LotkKEVApJdeytZm8xcsNIJJMWd/FN8b/rLjMBYupUk8Fu376kY/HxZvJUyplw1ux+\nH31kFjC3xervt1K1Kpw7l74v3h7Xrxs/flqra9ly++3m2sOGpT42c6b557j3XtftUBTFPV55xUy+\ncpQvH9h6cisNyxjxb1KuCYfOHeBck3r2o/s6djSaZIu115+Od2BT1CaalGti91inap24cP0C6yPX\np/153MQ3xd9K2bLwf/9nZrxa735r10LFivZz1t97L4wcae7E+/ebsoMHzR/ZdjA1Vy6oX99118+Z\nMybnzuOPm3w6zjJypLkB2d4ARGDIEO31K0pWU7as+R1PnuywyraT22hUthEAeQLz0CK2NKEtg+1X\nbt/eDALfvJlUtn69U+kYNkZupGm5pnaPBeYKpG/TvozckHKKlWfwbfEHeOMN41ObPdvs//lnapeP\nLT16wP/+Z3rmR44kLY6QK8VHbdjQvvhfvmw/hjc83Pjm27QxeftdoVQp+Osvk89j7FhTtsQyxSG9\niCFFUTxPnz7mt5iYmOrQtbhrHDp/iDolk8bh2uyLZUX5OPttFS9uxuz++SepzAl/P8CmE5scij/A\n8w2eZ+XhlRy9eDTdtrIDqYeyV60SKV/epDOtUkUkLCz94e8xY0TuvFOkcWOROXNSH58wQaRnz9Tl\n3bqJ5M0rcv/9IiNHihw5IvLnnyZWeOpUp0beHXL4sEnfMH68SMuWGW9PURT3SEwUqV/fpGBPwfrj\n66XhDw2TCs6fl7DK+aX66GqO2xs4UOTtt832jRsm84Cj2f8WrsVek/z/yy/XYq+lWe+dJe/I+8ve\nT7MOmZTeIaux/+meecaETlWo4PwU7SFDTC4feykWduwQqVkzednevSIlS4qcPWty9rzwghH9smVF\nNm507prpsW+fSJkyIpUrOz3ZRFGUTGDmTDNpq2FDkXffFVm4UOTKFRm7cay89OdLSfXmzpWEdm2l\n+NDicvzScfttbd0qUqOG2V6/3oSep8PaY2ul0Q+N0q136PwhKT60uMTEOs75RbYJ9bTHsGFmuvYj\njzi/TOGHH5pRd2tIqC21asHRo8mXfhs2zKycU7y4GaQdP964nA4dMqsFeYLq1c3avFOnagZORfEm\nXbuawI9vvoHChU3kXb16bDu6/tZgLwArV5LrgQdpXbm1Cfm0R4MGZmzx0CGnXT4bIzc6HOy15c6i\nd3Jv8L38Fv6bs5/MKfxH/MuUgVmzzIw4VyhVyn55njzGT7djh9mPjDTLL/bpk7xe7tyQL5/r9qZF\n9eoez82tKIob5M1rwjc//dQM2nbuzNatC2iUQvxp0yYp3t8euXIlRf14yN9vS58mfRizaYxHwz79\nR/wB2rVLvXBKRrCd7DViBPznP0mLLSiKkuO4MWgA+/Ncpt4SyzygU6fgxAm4+24T7394hWMB7tiR\n/Stnsnvvaucjfco7J/4P3vkgsQmxrD662tmPki7+Jf6exjrZ6+xZs/j6O/aWJFYUJaew6/JBqhWv\nRv5+n8HeveZpoFUrCAykWrFqBBDAnrN77J/84IN8nH8dL9x7BqlaNc3rRMdEc/76eaoXr+6UXbkC\ncvFGkzcYs2mMqx/JcZsea8kfsYr/6NHG/6cpkRUl2zB913RmR8zm4o2LTp+z9cRWGlZubjIL9Oxp\nZu5aUjoEBATwRK0n+H2X/XQtV/MHsrwKnC6ahzWR69K8zqaoTTQu15hcAc5LcK/6vVh5eCXHLx1P\nv7IT5Gzxr1fPzB4eN84sn6goSrZgy4kt9F3Sl5+2/UTwyGDu++U+BqwawLlr59I8b9vJbTQq08jM\nAq5Y0aSLscnn06NeD6btmmbX9bNw/0LuLVibj0p3Y/i6tNfg3hS1iSZlXQsiKXhbQZ656xnGbRnn\n0nmOyNninz+/Wfj8gQdMygdFUVIxfO1wfgv/jRvxN7xtilMkJCbw6oJXGfbgMJY8s4Qz753h81af\nsyFyQ7qifCutQ0CAWXu3b1+oXfvW8cZlG5MoiWw7mTo32MyImXR7sC+9XvuBjVEb2RPtwD0EbIxy\n3t9vyxtN3uDnbT975G+Rs8Uf4PPP4csvvW2Fovgku8/s5qv1XzFpxyQqjKzA+8ve58C5A942K03G\nbRnH7Xlv57n6zwGQP09+2lVpxyctP2HxwcUOz4tNiCUiOoL6QfVNQfHiMGpUstDygIAAnqr7FNN3\nTU92bkxsDMv/Xc6jNR8lf5789L6nN1+v/9rudUQkzZw+aVG9eHUalW2U6vruoOLftavp/SuKkoqv\n1n/Fm03eZNmzy1j34jpyBeTivl/uo8bYGjz+++N89vdn/L7rd05dPZWpdqw7vo53l76bbr0TV07Q\nP7Q/4zqNIyDFfKCm5ZoSeTmSyMuRds8NPx1O5aKVuT2vg6VZLfSo24Ppu6eTKEmpIRYeWEiz8s0o\nXsDk/OrduDez9syy+70cPH+QgrcVpPQdpdP9PPZ4p9k76bqvnEHFX1FyEBeuX+Dg+YNO1Y26HMWf\ne//ktcavAVC1WFWGth1K1DtRzHlyDt3rmPWmp++eTp3v6vDO0nc4ffV0ptg9fdd0vt7wNWuOrUmz\n3jtL3+HlRi9Tu2TtVMcCcwXSvkp7lhy0v3T4nD1z6FStk91jttQpVYei+Yqy9tjaW2UzI2bSrXa3\nW/slby9Jj7o9GLtpbKrzN0U5H99vj7ZV2vJu8/RvhOnhjPgHA38Du4FdwJuW8mLAcmA/sAwoYnNO\nP8wi7XuBdjbljYCdlmPfZMRwRVFc49KNSzww+QEe+u0h4hPj063/zcZveK7+cxTLn3zuS57APNQp\nVYfudbszoPUA5nafy67XdhGfGE+tb2vx4fIPOXvtrEdtX3l4JW83e5u3lryVrMdty7JDy9gYtZFP\nWn7isJ2Hqj7EogOLUpWLCNN2TaNnvZ5O2dOjrhn4BePyWXZoGY/WfDRZnbebvc0PW38wS0DasDHK\nuZm9mY0z4h8HvA3UAZoBrwO1gI8w4l8dWGHZB6gNdLe8dwC+A6zPX+OAF4FqlpemtFT8ggvXL/jN\ngKc9bsTfoMv0Ltxb/l6CCwUzeYfjdMZgbhTjw8bzdjPnZtSXKViG0Q+NZserO7h88zJVRlehy/Qu\nTNs5LZX4ucrpq6eJuhLFsLbDyBOYx67t0THRvLbwNb7t+C0F8hRw2Fb7qu1ZeXglsQmxyco3RG4g\nX+58Sf7+dHiq7lPMiphFXEIciw4somm5ppQoUCJZnWrFq9GyYktGrBvBtbhrt8oz2vP3FM6I/ylg\nu2X7KrAHKAc8AkyylE8CrLe9LsA0zE3jCHAQaAqUAQoCmyz1Jtucoyg+iYgwZccUqoyuQp9FfdI/\nwQeJT4znqVlPUaZgGcZ0HMPA1gMZsGoAN+NvOjznx60/0qFqByoWqejStYILBzPu4XEce+sYT9R6\nginhUyj3dTmenfssl244XjwlLf4+8jetKrYid67cjGo/iv9b+X/JbijRMdG0mdyGnnV70rFaxzTb\nKnV7KaoXr57MZQMwdedUetbtmWqcwBGVi1bmzqJ3suLwCmZEzEjm8rGlf6v+LDywkJLDS3L3D3fz\n6oJX2XlmZ/LcQV7CVZ9/JeBuYCMQBFgdfKct+wBlAdsRlUjMzSJleZSlXFF8kqjLUXSe1pmv1n/F\n7Cdns+DAArac2OJts1xCRHh5/svciL/BpEcnmQHbCvdRq2QtxoeNt3vOzfibjNo4ivebv+/2dQvn\nK8xz9Z9j0dOL+PfNf7kjzx20nNiSE1dOuNzWysMraVPZxNo3Ld+UNpXbMPifwYAR/gcmP0CXGl0Y\n0HqAU+11rNYxmesnPjGeGREz6FGvh0t29ajbg/Fh41l2aBmP1XrMbp16QfXY9N9NnPvgHN93+p7a\nJWvTv1X/dAeVswJX0kreAcwG+gJXUhzzaD7p/v3739oOCQkhJCTEU00rilPM3D2T1xe9zuuNX2dO\n9znkDczLl22+5M3Fb7LmhTUuzcx0xIXrF9hyYgv7z+03r/P7KX1HaV5v/Dr3lL3HA58ChqwZwp6z\ne/jr2b/IG5j3VvnA1gPpMr0L/2nwH/LnyZ/snKk7p1K3VF0alG7gERuKFyjOd52+Y/Cawdz3y30s\nfnoxNUvUdPr8lYdX8mbTN2/tD35gMPW/r89jtR7jhT9foHP1zgxsPdDpXnvHah35z5//YXi74bfa\nr1i4IlWLuTbX58k6T/LOsndoU7lNKpdPSvLlzkfT8k3diu23R2hoKKGhoR5pKz3yAEuBt2zK9gLW\nWKUyln0wvn/btQmXYNw+pTEuIys9gO/tXCvd/NaKkpkkJiZK0PAgWX98fbLyhMQEuefHe2Ty9skZ\nvsbRi0elwsgK0mpCK3ll/isyYt0Imbd3ngxdM1QqjKwg9/58r0wNnyqnrpySvw79JcPXDpees3tK\n1xldJdHJ9SxOXz0txYYWk8MXDts9/uj0R2XEuhHJyqIuR0mNMTVk+aHlGf2IdpkQNkGChgfJumPr\nnKp/5MIRKTW8VKrP/EXoFxL4RaB8/NfHTn8fVhISE6TEsBJy5MIRERF5/o/n5et1X7vUhpWOv3WU\niWET3TrXk5BJi7kEYPzzKReSHAZ8aNn+CBhi2a6NGSPIC1QGDpE04LsRcyMIABZhf8DX29+jksPZ\ncWqHVB1d1e6x9cfXS9kRZeXyDTsLBDnJqSunpNroajJq/Si7x+MS4mROxBxpPbG1FBpcSFr80kL6\nLOojE8ImSM2xNWXVkVVOXafv4r7SZ1Efh8fDT4VLqeGl5MrNK3Ll5hX5/O/PpdjQYm4Jqiss2r9I\nig4pKjXH1pSWE1pK1xld5Y2Fb9hdKGVC2ATpPrN7qvJrsddkdsRst+18Zs4zMm7zOLked12KDCki\nUZej3GonNj7WrfM8DW6IvzPPSS2A1UC4zQX6YQZuZwAVMAO7TwLWDEofAy8A8Rg30VJLeSNgIpAf\nI/5Jz3JJWD5LqkJemvcSnap34vFajzthtmLL/nP7WX98PXkC83B7ntu5Pe/tFM1XlIZlGjr9uJxT\n+GrdVxy+cJhvO31r9/jzfzxP0O1BDG071OW2L1y/QOtJrXm81uN81uozl88fuX4kYafCmPxY2tE6\nRy8epeGPDYnoHUHQHUEO6/WY3YNrcdfYHLWZ1pVb82WbL6lUpJLLdrnK5ZuXibwcyZmYM5yJOcOi\nA4u4mXCTaU9MS1bvubnP0aJCC15u9LJHrz9t5zSm7ZrG8w2eZ+ymsazs5SBPv59g+Q279EP2xV+9\nXfGfvGMyry96nebBzVn6zFI7pym2iAirjq5i/r75zN8/n5i4GFpVbIUgXI29SkxsDMcuHaPk7SUZ\n0W4EzYObe9tkn6HdlHa83vh1utTsYvf4qaunqPtdXda9uM7plLxg4sHbTmlLs/LNGNFuhFs33bPX\nzlJ1dFWOvHWEIvmKOKz34p8vUqZgGf7X5n9ptnfg3AH6rejHh/d9SONyjV22x1Ncjb1KtTHVWPz0\n4ltjDSJC8MhgQp8Pddkfnx7nrp2j8jeVaVWpFV1qdOGlhi95tP2sxh3x90VSPdIcvXhUSg4rKWuP\nrZVCgwtJdEx0Vj5RZTrhp8Ll478+loTEBI+1OXL9SKk6uqp8EfqFbD2x1e7jcUJigkzaPknKf11e\nnvj9CTlw7oDHrm9l0OpBcurKKY+3m1lci70mdwy6Qy7duJRmvTEbx0jRIUXl2TnPytw9c9NcX1VE\nJDomWlpPbC0v/PFChl0qT858Ur7d9K3D43ui90iJYSXkwvULGbpOVjN6w2jp+FvHW/v7zu6T4K+D\nM80Fde/P90regXnl/LXzmdJ+VkJ2WcA9PiH+1odKSEyQ1hNby+B/BouI+cf/YcsP3vqOPc6RC0ek\n/NflpdroavLpyk890ua5a+ek5LCSsvvMbqfqx8TGyJerv5RiQ4vJlB1TPGKDiMjWE1uF/sjI9SM9\n1mZms/TgUrlv/H1O1Y28FCljN46VNpPaSKHBhaTn7J6y6/SuVPX+OfqPBH8dLO8ve1/iEuIybOOy\ng8vk7u/vdni864yuMuSfIRm+TlZzI+6GVBpVSf45+o+IiIzbPE56ze2VadcbtHqQPDb9sUxrPysh\nu4h/uynt5GzMWRERGbV+lDQf3/zWDWHW7lnywKQHvPk9e4zomGipMaaGfLPhGzl15ZRUGFlBft/1\ne4bbfWvxW/Lq/FddPm/3md0SNDxI/tjzR4ZtEDE36gcnPyj3/3K/220kJibK5qjNHrHHGd5b+p58\nEfqFy+dFx0TL0DVDpeSwktJzdk/Zd3afJCQmyKDVgyRoeJAs2LfAYzYmJCZIxZEVZeuJramObYna\nImVHlE33ScRXmRg2UVr80kISExOl24xuMmn7pEy7VnxCvFy9eTXT2s9KyC7i/97S96TSqEoybec0\nKTGshBw8d/DWh7wWe00KDy4sp6+e9uJXnXGu3rwqTX9qKh8u//BW2bYT26TEsBJ2f9TOcuDcASk+\ntLjb38/mqM1SclhJ+evQX27bYGtHdEy0FB5c2G3Xz4J9CySgf0CWuY7uGndXqhBPV7h045IMXDVQ\nig8tLnW+rSP3jb/PbhRLRhkQOkBeW/BasrIL1y9Ik5+ayHebvvP49bKK+IR4qTW2lizYt0BKDCsh\nxy4e87ZJfgHZRfxFRKbvnC4Fviwg32/+PtUH7TGrh1//g8fGx0rH3zpKr7m9UvkzZ+6eKcFfB8vJ\nKyfdavvx3x+XQasHZci+0MOhUnJYSdlwfIPbbbwy/xX5ZMUnIiLSfWZ3+XHLj261EzIxREp/VVrG\nbR7nti3OcuLyCSk6pKhHXDMXrl+QORFzPNKWPY5dPCZFhxS91cOPvBQpdb+rK28uetOjY0feYHbE\nbCnzVRmpNrqat03xG8hO4i8iDmOp5+6ZKyETQ7Lqe/U4Q/4ZIm0nt3UYI/z5359L8/HNHQpHYmKi\nPP774/L474/L/rP7b5WvPrJaKoysINdir2XYxgX7FkjQ8CDZeXqny+daRfTM1TMiIvL7rt/loV8f\ncrmdzVGbJfjrYJm+c3qarr6L1y96JN568vbJ8sTvT2S4nayi428dZdL2SbL7zG6pMLKCDF0zNFPj\n87OKxMREafxjY3l53sveNsVvILuJvyOsEzNOXD6RrHzdsXXS769+WfoDiImNkcqjKssj0x6RqeFT\n5crNK2nWvxF3Q8p8VUbCT4U7rJOQmCBtJrVJNfvSyqzds6TOt3Vk0OpBUnxocem7uK9Ex0RL4x8b\ne3TAdsqOKVJ5VOVb4y/O8uHyD+WNhW/c2r9847IUGlxILl6/6FI73Wd2lxHrRkhMbIwUGlzo1s3E\nlsTERKk/rr6UG1FOBq4amMrdtTd6r3y19iunfMfPzHnG7pOmrzInYo7UHFtTSg0v5ZFZx77EkQtH\n3J54lRMhp4i/iPmhjtk45tb+2mNrpeSwklJrbC2X3B6D/xksC/cvdPnLtjI1fKq0nthaJoZNlId+\nfUgKDS4k3WZ0c+i2Gb9tvLSf0j7ddq0+83/P/5us/OrNqxL8dbCEHg4VETOF//WFr0vBQQXlnh/v\n8fgj//vL3pc2k9o47b64eP2i3ZQCD099WH4L/83p6/57/l8pNrTYrae/bjO6yU9bf0pVb+W/K6XW\n2FoSdjJMXvrzJSkypIg8N/c5eXvJ21JtdDUpO6KsvPDHC1J0SNE0bz7WlA4pv29fJjY+VlpOaClL\nDy71timKlyEnif+8vfNuRZGsO7ZOSg4rKYsPLJbIS5FSdkRZp6Irtp7YKqWGl5Kg4UFuRxV0+LVD\nMlE7G3NWXpn/inSd0TVV3YTEBKk1tpbTg6mDVg+SDr92SPYk89Hyj6Tn7J6p6u4/u99hDpeMEJ8Q\nL+2ntJe3Fr/lVP3B/wyWZ+Y8k6r8l22/2P1OHPHmojflg2Uf3Nr/fdfvdm+anad2Thb6ezbmrAxf\nO1wGhA6QbSe23frueszqIcPXDnd4vbRSOiiKr0NOEv8bcTek6JCiMnP3zFvCb8X6FLA3eq/D8xMT\nE6XVhFby/ebvJeJMhAR/HezQzeKIE5dPSJEhRVKF1V2LvSZVvqki8/fNT1Y+f998ufv7u512S8XG\nx0q97+rJ1PCpImJcGMWHFs/yx+Hz185L1dFV7d4gr8ddl81Rm+XHLT/Kawtek2JDi9l1aUXHREuh\nwYWcGo84d+2cFB1SVCIvRd4qu3LzihQcVFDOXTt3q2z/2f1SYlgJp8Iat0RtkfJfl3c4NjB87XDp\nvaB3uu0oii9CThJ/EZFec3vJbQNvSyb8Vn7a+pPUGFPD4aP+nIg5UufbOrfcGUcvHpWaY2vKR8s/\nclqcR6wbIc//8bzdY8sOLpOKIysmiyNuNaGVS64PEZENxzdI6a9Ky9mYs9J2clu3sw9mlF2nd0mJ\nYSVk3OZx8kXoF9JtRjepNbaW5P9ffrlr3F3Sa24vGbV+VJox+a0ntnZqDsGXq7+0O7nnsemPyYSw\nCbf2+yzqI/3+6uf0ZwiZGCK/7vjV7rG2k9t6bH6DomQ15DTxP3juYJqpYV9b8Jo89OtDqXqGN+Nv\nSpVvqqTylVoHTZ0dM6g/rr6s/Helw+NPz35a3lv6noiIbIrcJBVGVnArKqXPoj5y17i7pO53db2a\nRXDxgcXy6PRHpd9f/eS38N9kx6kdciPuhtPnj9k4Jt0Zm9fjrkvpr0rbjTL6Lfw3eXjqwyJiQimL\nDinqUgz9gn0L7D55rTm6RooNLebygLSi+ArkNPFPj5vxN+Xp2U/LXePuShYSOWLdiGQ5RGyxDrSm\nlxdl+8ntEvx1cJoDrKevnpaSw0pK2MkweXLmk2732i/fuCz1vqvndCpfX+X4peNSbGixNG9gk7dP\nlnZT2tk9dunGJSk4qKBcvH5RRqwbIT1m9XDp+gmJCVJzbE1Z8e+KW2V7o/dK0PAgu0+PiuIvoOKf\nmnAhuzoAAAmGSURBVMTERPlu03dSclhJmR0xW6JjoqXEsBIScSbC4Tm95vaS/n/3T7Pdd5e+Kx//\n9XG61/9p609S+9vaUnxo8QzlgM8uNPmpiSw7uMzh8ZYTWsrsiNkOj3ee2lkmhk2UiiMrysbIjS5f\n/6etP9268Z+8clIqj6osv2z7xeV2FMWXQMXfMZsiN0nFkRWl9re15fWFr6dZ98C5A2lmRYxLiJMy\nX5WRPdF70r1uQmKCtJzQ0qkbRU7g203fSpdpXewe23d2nwQND5Kb8Tcdnj9p+yQpNbyUNB/f3K3r\nX4+7LkHDg2Rj5EZp+ENDt/L4KIqvgRvin/GFSP2ExuUas/XlrbSv0p7+If3TrFu1WFU6VevE6I2j\n7R5f8e8Kyhcq79Q6pLkCcrHsmWUMbDPQHbOzHc83eJ71kevZe3ZvqmPjt43nufrPJVtrNiWdq3fm\nwvULvNX0LYd10iJf7ny83vh1QiaG0LB0Qz5t+alb7SiKv+OLyf8tNzLvcuDcAZr/0pyDfQ5SOF/h\nZMeenvM095a/lzeavOEl6/ybAasGcOzSMX5+5OdbZXEJcQSPDGbV86uoUaJGmudvjtpMwzINCcwV\n6Nb1z18/zzcbvuHTVp+SO1dut9pQFF8iW6/k5Q16/dGLqkWr8mmrpN7htpPbaDOpDQffPEiJAiW8\naJ3/cvbaWaqPqc7u3rspU7AMAHP3zGXkhpGs/s9qL1unKP6HO+LvjNvnF+A0sNOmrD8QCYRZXg/Z\nHOsHHAD2Au1syhtZ2jgAfOOKkd7ik/s/YfSm0Vy6cYlNUZvoPK0znad1ZmT7kSr8GaBEgRI8Xe/p\nZG61n8N+9vul9BTFn3DmTnE/cBWYDNSzlH0OXAG+TlG3NjAVaAyUA/4CqmEGIzYBb1jeFwGjgSV2\nruczPX8wC0ivPb6WuIQ4PmrxES/c/QL5cufztll+z+ELh2n8U2MO9z3MxRsXqf99fSLfiaRAngLe\nNk1R/A53ev7OODz/ASrZu56dsi7ANCAOOAIcBJoCR4GCGOEHcyN5FPvi71MMfmAwKw6voHud7tyW\n+zZvm5NtqFy0Mm2rtOWnbT8RExvDU3WfUuFXlCwkI6NdfYDngC3Au8BFoCywwaZOJOYJIM6ybSXK\nUu7zlCtUjufqP+dtM7Il7zd/ny7TuxAYEMic7nO8bY6i5CjcDfUcB1QGGgAngREes0jJMTQs05Aa\nxWtQLH8xGpZp6G1zFCVH4W7P/4zN9s/AfMt2FBBsc6w8pscfZdm2LY9y1Hj//v1vbYeEhBASEuKm\nmYqvM7bjWC7euOhtMxTFrwgNDSU0NDRDbTg7QFAJI/DWAd8ymB4/wNuYAd6eJA34NiFpwLcqZsB3\nI/Amxu+/ED8Z8FUURfF1MmvAdxrQCigBHMdE+oRgXD4CHAZesdSNAGZY3uOB3iRNO+4NTATyY6J9\nfH6wV1EUJbuik7wURVH8nMya5KUoiqJkM1T8FUVRciAq/oqiKDkQFX9FUZQciIq/oihKDkTFX1EU\nJQei4q8oipIDUfFXFEXJgaj4K4qi5EBU/BVFUXIgKv6Koig5EBV/RVGUHIiKv6IoSg5ExV9RFCUH\nouKvKIqSA1HxVxRFyYGo+CuKouRAVPwVRVFyICr+iqIoORBnxP8X4DSw06asGLAc2A8sA4rYHOsH\nHAD2Au1syhtZ2jgAfOO+yYqiKEpGcUb8JwAdUpR9hBH/6sAKyz5AbaC75b0D8B1JiwqPA14Eqlle\nKdv0G0JDQ71tglOonZ5F7fQs/mCnP9joLs6I/z/AhRRljwCTLNuTgEct212AaUAccAQ4CDQFygAF\ngU2WepNtzvE7/OUfQu30LGqnZ/EHO/3BRndx1+cfhHEFYXkPsmyXBSJt6kUC5eyUR1nKFUVRFC/g\niQFfsbwURVGUbEYlkg/47gVKW7bLWPbB+P4/sqm3BOP2KQ3ssSnvAXzv4FoHSbqh6Etf+tKXvtJ/\nHSSTqERy8R8GfGjZ/ggYYtmuDWwH8gKVgUMkDfhuxNwIAoBF+PGAr6IoSk5gGnACiAWOA//BhHr+\nhf1Qz48xd6G9QHubcmuo50FgdKZbrSiKoiiKoiiKb9IB87RwgCSXki/g6iQ3bxEM/A3sBnYBb1rK\nfcnWfBj333YgAhhsKfclG20JBMKA+ZZ9X7TzCBCOsdMaSu2LdhYBZmHG/iIwLmBfs7MG5nu0vi5h\nfke+ZieYybS7Mbo0FbgN37QzXQIx7qBKQB6MONTypkE23A/cTeoxjw8s2x+SNObhTUoDDSzbdwD7\nMN+hr9lawPKeG9gAtMD3bLTyDvAbMM+y74t2Hsb86G3xRTsnAS9YtnMDhfFNO63kAk5iOlW+Zmcl\n4F+M4AP8DvTC9+x0insxkUFWUkYNeZtKpI52ss5tKE1StJMv8QfwIL5rawFgM1AH37SxPGZcqzVJ\nPX9ftPMwUDxFma/ZWRgjVinxNTttaYeZ4Aq+Z2cxTOeuKOZGOh9oi+/Z6RRdgZ9s9p8BxnjJFntU\nIrn42854DiD1DGhvUwk4iplV7Wu25sI82V3B9FTA92wEmIl54mtFkvj7op3/YlwUW4D/Wsp8zc4G\nGHffBGAb5rd+O75npy2/AL0t275o58uY39AZYIqlzCU7fSWrp3jbgAxgjbP1Fe4AZgN9Mf8ctviC\nrYkYMSgPtMT0rG3xBRsfxvyowkgKVU6JL9gJcB/mJvUQ8DrGTWmLL9iZG2iIyfXVEIgh9ZO9L9hp\nJS/QGdMBSIkv2FkFeAvTySuL+c0/k6JOunb6ivhHYXxrVoJJng7C1zhN8kluZ7xoiy15MMI/BeP2\nAd+19RKwEBMC7Gs2NsfkrzqMCXVug/lOfc1OMH5pgGhgLtAE37Mz0vLabNmfhbkJnMK37LTyELAV\n852C732f9wDrgHNAPDAH4zp36fv0FfHfgsn0WQlz1+1O0iCbLzIPM8CC5f2PNOpmFQHAeEwkxSib\ncl+ytQRJEQj5MX7KMHzLRjBzVYIxExWfAlYCz+J7dhbAuPbAuFHaYdyTvmbnKcwcoeqW/QcxkSrz\n8S07rfTA3PSt+Nr3uRdohvkNBWC+zwh89/tMl4cwgxgHMWFMvoKrk9y8RQuMS2U7SaFqHfAtW+th\nfL7bMeGJ71vKfcnGlLQiqSPia3ZWxnyX2zHhvdbfja/ZCVAf0/PfgempFsY37bwdOEvSTRV8084P\nSAr1nIR56vdFOxVFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURR7/D+dTFae\nkSkZsQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x115123950>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(s_mod[-80:],'r',output[-80:],'g')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make lightcurves using `Lightcurve` class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 598,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "time = lc.time[delay:]\n",
    "lc1 = Lightcurve(time, s_mod)\n",
    "lc2 = Lightcurve(time, output)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute crossspectrum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 599,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cross = Crossspectrum(lc1, lc2)\n",
    "# Rebin the cross spectrum for ease of visualization\n",
    "cross = cross.rebin(0.0075)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate time lag."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 600,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "lag = np.angle(cross.cs)/ (2 * np.pi * cross.freq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot lag."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 601,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG2lJREFUeJzt3XmUFNXdxvHvAKLiAiLrIIohLhBRwQguGNoFBhRfglET\nNDEuUaOJ0SQaI74n06PGJWo8B5fXKO7HKJgjKAiyKD0gilFARASFUQjIpgQUhJGt3z9+3UwzzNLd\nVdXVVfV8zukzPd1dXZei55k7t+79FYiIiIiIiIiIiIiIiIiIiIiIiIiISMA9CawB5mc8FgdWAHNT\nt4GFb5aIiHjlNKAnuwd/OfAHf5ojIiINaeLCe8wA1tfxeIkL7y0iIi5zI/jrcx0wD3gCaOXhfkRE\nxAdd2H2opx3W4y8B7sDCX0REikAzj953bcb9kcC42i/o2rVrsqqqyqPdi4iEVhXwfSdv4NVQT8eM\n+0PZ/a8BAKqqqkgmk7q5dOv3y36+tyFMt/Lyct/bEKabjqd7N6Cr04B2o8f/AtAPaAMsx2b0xIDj\ngSTwOXC1C/uRBlQurfS7CSISEG4E/7A6HnvShfcVEREPeDmrRwqpi98NCJdYLOZ3E0JFx7O4+DnX\nPpkarxIXlFSUkCzX8RQJu5KSEnCY3erxi4hEjII/JMr7lfvdBBEJCA31iIgEiIZ6REQkZwp+EZGI\nUfCLiESMgl9EJGIU/CERT8T9boKIBIRm9YSEFnCJRINm9YiISM4U/CIiEaPgFxGJGAW/iEjEKPhD\nQrV6RCRbmtUjIhIgmtUjIiI5U/CLiESMgl9EJGIU/CIiEaPgDwnV6hGRbGlWT0ioVo9INGhWj4iI\n5EzBLyISMQp+EZGIUfCLiESMgj8kVKtHRLLlxqyeJ4FzgLVAj9RjrYFRwGHAUuBCYEOt7TSrR0Qk\nR8Uyq+cpYGCtx/4MTAGOBN5IfS8iIkXArXn8XYBx1PT4FwH9gDVAByABHF1rG/X4RURyVCw9/rq0\nx0Kf1Nf2Hu1HRERyVIiTu8nUTUREikAzj943PcSzGuiInfjdQzwe33U/FosRi8U8ak74xRNx4rG4\n380QEZclEgkSiYSr7+nVGP/fgHXAPdiJ3VbseYJXY/wuUq0ekWgoljH+F4C3gaOA5cBlwN1Af+BT\n4IzU9yIiUgTcGOoZVs/jZ7nw3pKvTZtgn32gmVejeSISVFq5GwbV1fb1wQfhl7+E7t3hwANhxAh/\n2yUiRUndwaDZvh0WLID33qu5LVoENwPz50PfvvD738PUqbB0qd+tFZEipOD30/bt8O23sHGjDc1k\n3mo/9uWXMHs2zJsHnTvDiSfa7bLL4LjjKH/3Hsic1bN4Mcya5ds/TUSKl67AVShbt8KMGfDaazBh\nAixbZo/tv/+etwMO2POx1q2hZ0844QQbxmnMzJlw003w9tve/9tEpGDcmNXjpyRxu5VPK0/WVj6t\nvN7HA7fdNd2SyZYtk8k+fZLJ225LJmfPTpZPHp4sf/Mv3rXzs8+SycMOK8y/T9tpO21XsO1wYUGs\nevxuSibhgw9g/Hi7ffIJ9O8PgwfDoEHQrl3h2lJdDS1b2teSwHYORKQWN3r8Cv5M1dXw3HOwdi3s\ntRc0b17zNfN+7cfWrYOJE20YZ7/9LOgHD7YTrc2b+/fvOfhg++XTpo1/bRARV7kR/Dq5C7BtGzz9\nNNx2Gxx/PBx7rI2/b9tmXzPv1/W1RQsoK4Mbb4Qjj/T7X1OjY0dYtUrBLyK7iXbw79gBL74I5eXQ\npQu89BKcdJLfrcpLnbV6Skth5Uro0aPObUQkmqK5gCuZhLFjrXf/8MPw+OM27z2goQ9QUVmx54Pp\n4BcRyRCtHn8yCVOmwP/+rw3T3HUXnHNOeE9+KvhFpA7RCf6ZM+HWW23M+/bb4fzzoUnI/+ApLYWF\nC/1uhYgUmZAnHzBnDpx9Nlx8sdWxWbAALrww/KEP6vGLSJ3C2ePftg1efRUee8zq1wwfDmPGwN57\n+92ywlLwi0gdwhX8VVUwciQ89RQcfTRcdRWcd56VJw658n7lez6o4BeROgR/AdfWrfDKK9a7/+AD\nuOQSuPJKC/6oS9cCqq6OxtCWSAREewHXkiXWu3/6aejWzXr3Q4dGoneftebNoVUrq+zZvr3frRGR\nIhGsbuDWrTB6NJx1FpxyipU1rqyEadNg2DCFfl3Sq3dFRFKC1eM/4wxo2hSuucZ691E7WZuP9Dj/\n8cf73RIRKRLBCv7ly62m/aGH+t2S4NAJXhGpJVhDPTt36iRlPeKJeN1PKPhFpJZgpejOneEtr+BQ\nnbV6QMEvInsIXvCrx58bBb+I1BKsFFXw507BLyK1BCtFFfy5U/CLSC3BSlEFf+7at4evvrI1DyIi\nKPhDo85aPQDNmtm1d9euLWyDRKRoBStFFfz12uOyi5k03CMiGYKVogr+/Khsg4hk8Hrl7lLgG2AH\nsA3o7ejdFPz5UY9fRDJ4HfxJIAb815V3U/DnR8EvIhkKkaLuLbVV8OdHwS8iGbxO0SQwFXgfuNL5\nuyUV/PWot1YPKPhFZDdeD/WcCqwC2gJTgEXAjPST8Xh81wtjsRixWKzhd1OPv14VlRX1z+xR8IsE\nViKRIJFIuPqehax4Vg5sAu5PfZ/7pRebNIEdO1SorQ4lFSUky+s5nqtXw3HHwZo1hW2UiLjOjUsv\netl9bgEckLq/HzAAmJ/3uyWTdlPo565tW1i/HrZt87slIlIEvBzqaQ+MydjP88DkvN9NoZ+/pk2h\nXTvr+Xfu7HdrRMRnXgb/54B71/vT+L4z6XF+Bb9I5AUnSRX8Daq3Vk+aVu+KSEpwklTB36AGa/WA\nZvaIyC7BSVIFvzMKfhFJCU6SKvidUfCLSEpwklTB74yCX0RSgpOkCn5nFPwikhKcJFXwN6jBWj2g\n4BeRXYKTpAr+BlVUVjT8goMPho0b4bvvCtMgESlawUlSBb8zTZpAhw6ayy+FsWEDjBsHN94IJ54I\nV1/td4skQ3CSVCWZndNwj3jlq6/g5ZfhhhugZ09bIT5iBLRqBX/5C4weDdu3+91KSfG6LLN71ON3\nTsEvblmzBiora27Ll8Mpp0C/fvDII3DCCdC8ec3rDzsM3n0XTj3VvzbLLsEKfhVpc0ZlGyQfO3bA\nggXwzjswa5Z9XbMG+va1oL/0UuvlN2sgTsrKYPJkBX+RCFbwq8dfr0Zr9YB6/JKdtWutd54O+fff\nt07DySfDSSfB9ddDjx5W9TVbAwbArbdCRSOTEKQgFPwh0WitHrDgnzbN87ZIgGzfDvPm7d6b/+or\n6NPHQv6mm6B3b5sV5kTfvvDxx/Df/0Lr1u60XfKm4I8S9fgFbJjm9ddhwgSYMqWmN3/66XDLLdCt\nm/s/a3vvDaedBlOnwoUXuvvekjMFf5Qo+KNp504brpkwwW6ffgpnnglnnw0PPGCfi0JIj/Mr+H2n\n4I8SBX90rF9vIfvaa9a7b9vWgv6ee+wEa+aMm0IpK4N779XV9IqAgj9KDjoItmyBzZuhRQu/WyNu\nq6qCf/0Lxo+3cfsf/cjCvqICDj/c79bBkUfaz/DChdC9u9+tibTgJKmCv0GN1uoB62WVlmpKZ5is\nXGnDNX362Dj9smU2e2bNGvsFcO21xRH6YJ+/sjKYNMnvlkRecJJUwd+gRmv1pGm4J/jWrYN//MNO\nxh5zDHz4Idx+u/2/PvIIDBwI++7rdyvrlh7nF19pqCdqFPzBtHEjjB0LL74Ib71l4X799fZ1n338\nbl32zjzTFnxVVwer3SGj4I8ard4Nji1bYOJEeOEF6yX/6Edw8cUwahTsv7/frctPq1a2+GvGDOjf\n3+/WRJaCP2rU4y9uq1fbTJxx42yx3Q9/CMOG2dBOWBY+pcf5Ffy+CU7wqzqnO0pLYf58v1shacmk\njdGPG2e3Tz+1YLzgAnjyyfCEfaayMrjqKr9bEWnBCX4VaWtQVrV6QD3+YlBdDYlETdg3bw7nngt3\n3WWrW/fay+8WeuvEE2HFCvscFmrxmOwmWMGvHn+9sqrVAwp+v2zYYCdnX30V3ngDjj3Wwn7SJDj6\n6Gh1apo2tZO8kyfbiV4pOAV/1Cj4C2fbNgu3Z5+1gD/zTBg6FB57DNq08bt1/kqP8yv4faHgj5oD\nD7T66hs3wgEH+N2a8EkmbdXsM8/YbJyuXeGSS+DRR23ltJiyMisIp59rX3h5xAcCi4DFwM2O300f\nEHdo9a43Vq2C++6D446DH//YfqnOmAEzZ9r1ZhX6u+vc2eoHzZnjd0siyaskbQo8hIV/d2AY0M3R\nOyr43aPhHnds3my9+kGDrPbMwoXw4IPw2Wdw221wxBF+t7C4DRig8g0+8SpJewNLgKXANuBFYIij\nd1TwNyirWj1pCn5nli2zGjidOtmQzi9+AV98AU88YZci1Oc0O6rb4xuvPqGdgOUZ369IPZY/BX+D\nsq7VAxrqydeyZTZs06sXtGwJH31kJY8vukjVTvPRrx/MnQvffON3SyLHq5O7yWxeFI/Hd92PxWLE\nYrH6X6zgd0/Hjurx52LpUrjzTit5/OtfwyefaFaOG1q0sMs7TpsGQ5wNCIRZIpEgkUi4+p5eBf8X\nQOeM7ztjvf7dZAZ/oxT87ikthdmz/W5F8fv8cwv8l1+2wF+82Pm1Z2V36XF+BX+9Yn36ENt3XzsR\nPmcOblyu3qvgfx84AugCrAR+ip3gzZ+C3z0a429YZuBfc42VUVDge6OszNY2iKmutstkzp69K+ip\nqrJFfr16Qc+eruzGq+DfDvwWmITN8HkCWOjoHRX87lHw1+2zzyzwx4xR4BdKjx5WhbSqytY8RM23\n38I770BlJUyfboHfrZuVtTjtNCu9/YMf2MXq0377W8e79XIB18TUzR0K/gZlXasHasb4de1Ts2SJ\n1ckZO9Zm6yxeHM7iaMWopKRmuOfaa/1ujfe+/trWdkyfbmE/fz4cf7yd6B4+HE45pSALK4OzclfV\nORuUda0esA9W06Y2m6JlS8/aVPTmzLGLj7/5pvXwFfj+KCuzC8yEMfjXr6/pzVdW2sSA3r0t6O+8\n005u+3C1tOAEv6pzuis93BO14E8mrTLm3XfDggXwhz/AyJEqX+Gns86yk+dbt1ql0iBLJm0h32uv\n2TWP5861ayH36wcjRtj1FTKHbXwSrOBXj9896eDv5mxBdWDs3AmvvGKBv2ED3HyzXc2qCH4II69t\nW1vlPGuWXWUsaKqrrTefDvvt22HwYPuMnX56UV7/WMEfVVE5wbt1Kzz/PPztb3a5wltusamDTZv6\n3TLJlF7FG5TgX7UKJkywoH/zTbvo/eDB1rk45piiH51Q8EdV2IN/0yZ4/HH4+9/tr5qHH7beV5H/\nQEZWWZkNu/31r363pG47dsB779lK7fHjbQbYgAHwk5/Y5yxgC/qCk6QK/gblVKsHwnvR9S1b4I47\n4PDD4e23babO5MlwxhkK/WJ20kl2cv3LL/1uSY2VK+Gpp+BnP4N27eDKK2365f33w5o1dkL65z8P\nXOiDgj80cqrVA+Hr8SeTFvLdu9s1bGfOhJdeghNO8Ltlko3mzSEWg6lT/WvDd9/ZsM3NN1t57WOO\ngYkTrWc/b55Nvbz3XjtRG/DLY2qoJ6rCFPyffGILXZYvtwqZZ5zhd4skH+lx/mHOFvnnpKrKhm8m\nTbITtEcfDQMH2oVzTjwRmgUnInMRnH+Vgt9dYQj+TZtsWGfkSFv8ct11ge+JRVpZmf1/ermwcMsW\nm847caLdNm2y/V50ETz5ZCCHbfKh4I+qIK/eTSZh1Ci46SY7YTt/vv17JNi6drWpj/Pn28Xo3bJ4\ncU3Qv/WW1bsZNAhGj7ZVs0H7/LtAwR9VLVrYD9n69cFarfrRR1arZMMGu/pV375+t0jcVFZmJ+Od\nBP/mzbv36jdvtuGbK66wz0yrVq41N6iCk6QK/gblVKsnLUjDPRs2wA032Pj9hRdaMSuFfvjke1Wu\nlSttZezAgdC+vS3UKy21ayh88YUN45x/vkI/JThJquBvUE61etKCEPw7d8LTT9tc/M2brczCtddq\nAVZYnX66reDdvLnx127bZpVUzz3XKljOmQO/+pWd5J8+3RbrRXQopzEa6omyYg/+r7+269muWgWv\nvmqzLCTcDjzQxuArK20cvi4ff2w9+Oeeg6OOgssvtzn1++1X2LYGWHCSVEXa3FfMwb9woVUxPPRQ\nm5Ov0I+O9Dh/pm++sRWyJ59sRd322gtmzLCe/aWXKvRzFJwev8oyu6+01C42UmzGjIGrrrLFMpde\n6ndrpNDKyuCSS+xn/q23bG3G2LF2fufWW20cP6Tz6wslOEdPQz3u69jRZj8Uix07oLzc/oSfONFK\n2Er09OplpRuOOMJW9F5xhV03oX17v1sWGgr+kIgn4rmf4C2moZ71661M8ubNVgyrXTu/WyR+adIE\nnnnGphn36aMhXg8EJ0kV/A3KuVYPFE/wf/SRjeEfdRRMmaLQFzj7bCvcptD3RHCSVMHvvo4dYfVq\nO7Z+eeklm8IXj8MDD6jkgkgBaKgnyvbe2y45uG6dXQWpkHbssBN1o0bZDI6ePQu7f5EIU/BHXXq4\np5DBv26dFcVKX9wiIoWxRIpFcJJUwe+NQo/zz5tn4/nHHWflcBX6IgWnHn9I5FWrBwob/K+/bvOz\nR4ywqxqJiC+CE/wtW9qYtNQpr1o9ULjgf/ZZ+NOf7GLUJ5/s/f5EpF7BCf7f/MbvFoRTaanVP/dK\nMgn33QcPPQTTplmxNRHxlcZOos7Li67v3Al//KMtxpk5U6EvUiSC0+MXb3g11LN1K1x2GfznP1ZM\n66CD3N+HiOTFqx5/HFgBzE3dBnq0H3HKi+DfuBEGD7byC5MnK/RFioxXwZ8E/g70TN1e92g/khJP\nxPPbsEMHWLvW5tS7Ye1aW4l7+OG2Knfffd15XxFxjZdj/CqyUUB51eoBK5Fw0EFWDdGpqio49VTr\n7T/6qErnihQpL4P/OmAe8ASgC10WMzeGe+bOhdNOs5O58biKa4kUMSddsilAhzoevxX4P+C21Pe3\nA/cDV9R+YTwe33U/FosRi8UcNEfylg7+Xr3y2/6NN2DYMOvln3eeu20TibhEIkHC5etmFKJb1gUY\nB/So9XgymUwWYPfRUFJRQrI8z+N55ZVWRuGqq3LfdtQo+N3vYPRo6Ncvv/2LSNZK7K9pR9nt1VBP\nx4z7QwEPVwiJY/kO9Tz1FNx4I0ydqtAXCRCvzr7dAxyPze75HLjao/1ISt61esCCf/bs3LaZOhX+\n/Gebo3/kkfnvW0QKzqvgv8Sj95V65F2rB3Jfvfvxx1ZW+aWXFPoiAaSSDZLbUM/atTZd8777NLwj\nElAKfsk++LdsgSFD7KLol+iPOpGg8nOytWb1FIvt26FFCyuxUN+iq507bXgH4J//1LURRHzixqwe\nLa0UC/s2bWDNGujUqe7XlJfDsmXw5psKfZGA009wSORdqyetoeGeZ5+F55+3i6io9o5I4Cn4QyLv\nWj1p9QX/9Ok2V3/8eGjXztk+RKQoKPjF1BX8ixfDBRdYb797d3/aJSKuU/CLqR3869bBOefA7bdD\n//7+tUtEXKfgF5MZ/N99Z8XWhgzJr36PiBQ1Bb+YdPAnkxb2rVvD3Xf73SoR8YCmc4aEo1o9UFO2\n4c47YcECqKyEpk3daZyIFBUt4BKzZg0ccoj9Apg1y/4CEJGiU8xlmSVo2raFnj1h3DiFvkjIqccv\nIhIg6vGLiEjOFPwiIhGj4A8Jx7V6RCQyNMYfEo4uti4igaExfhERyZmCX0QkYhT8IiIRo+AXEYkY\nBX9IOK7VIyKRoVk9IiIBolk9IiKSMwW/iEjEKPhFRCJGwS8iEjFOgv8CYAGwA+hV67lbgMXAImCA\ng31IllSrR0Sy5ST45wNDgem1Hu8O/DT1dSDwiMP9SBYqKiv8boKIBISTQF4EfFrH40OAF4BtwFJg\nCdDbwX5ERMRFXvTES4EVGd+vADp5sB8REclDs0aenwJ0qOPx4cC4HPajlVoiIkWiseDvn8d7fgF0\nzvj+kNRje4jH47vux2IxYrFYHrsTEQmvRCJBIpFw9T3dKNkwDbgRmJ36vjvwT2xcvxMwFfg+e/b6\nVbLBRfFEnHgs7nczRMRjbpRscLLxUGAE0Ab4GpgLDEo9Nxy4HNgOXA9MqmN7Bb+ISI78Dn6nFPwi\nIjlSkTYREcmZgl9EJGIU/CIiEaPgDwnV6hGRbOnkbkiUVJSQLNfxFAk7ndwVEZGcKfhFRCJGwS8i\nEjEKfhGRiFHwh0R5v3K/myAiAaFZPSIiAaJZPSIikjMFv4hIxCj4RUQiRsEvIhIxCv6QUK0eEcmW\nZvWEhGr1iESDZvWIiEjOFPwiIhGj4BcRiRgFv4hIxCj4Q0K1ekQkW5rVIyISIJrVIyIiOVPwi4hE\njIJfRCRiFPwiIhGj4A8J1eoRkWw5Cf4LgAXADqBXxuNdgC3A3NTtEQf7kCxVVFb43QQRCQgnwT8f\nGApMr+O5JUDP1O1aB/uQbH3udwPCJZFI+N2EUNHxLC5Ogn8R8KlbDRGHlvrdgHBRULlLx7O4eDXG\nfzg2zJMA+nq0DxERyUOzRp6fAnSo4/HhwLh6tlkJdAbWY2P/Y4EfABvzbKOIiLjIjZIN04A/AnNy\nfH4J0NWF/YuIREkV8H0nb9BYjz9bmb9A2mC9/R3A94AjgM/q2MZRw0VEpPCGAsuxqZurgYmpx38C\nfISN8c8GzvGldSIiIiIi4p2B2HTPxcDN9bxmROr5edh8/1y2jRonx3Mp8CH2F9i/vWtiYDR2LI8G\n3gGqsXNTuWwbRU6O51L02aytseN5MfYz/iEwEzg2h2091RQ7cdsF2Av4AOhW6zVnAxNS9/sAs3LY\nNmqcHE+wpV2tvW1iYGRzLNsCPwTuYPeg0mdzT06OJ+izWVs2x/NkoGXq/kDyzE4v5vH3TjVgKbAN\neBEYUus1/wM8k7r/LtAKmzaazbZRk+/xbJ/xvJ8X3Ckm2RzLL4H3U8/num3UODmeafps1sjmeL4D\nfJ26/y5wSA7b7uJF8HfCTvqmrUg9ls1rSrPYNmqcHE+AJDAV++G70qM2BkU2x9KLbcPK6THRZ3N3\nuR7PK6j5Sz+nbd2azpkp2+sp6jd9dpwez77Yorq22IK8RcAMF9oVRE6u9anrhO7J6TE5FViFPptp\nuRzP04HLsWOY67ae9Pi/wFbupnXGfvs09JpDUq/JZtuoyfd4fpG6vzL19UtgDPYnYVQ5+Xzps7kn\np8dkVeqrPpsm2+N5LPA4NsS7PsdtPdMMW1nWBWhO4ycjT6LmBEU220aNk+PZAjggdX8/bBbAAA/b\nWuxy+XzF2f1kpD6be3JyPPXZ3FM2x/NQbCz/pDy29dwg4BOsgbekHrs6dUt7KPX8PHav51/XtlGX\n7/H8HvYB+ABbVKfj2fix7ICNlX6N9ab+A+zfwLZRl+/x1Gezbo0dz5HAOmqud/LvRrYVERERERER\nEREREREREREREREREREREREREZFC+n+0KptY6+7dFwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x113fce550>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "\n",
    "# Plot lag-frequency spectrum.\n",
    "plt.plot(cross.freq, lag, 'r')\n",
    "\n",
    "# Find cutoff points\n",
    "v_cutoff = 1.0/(2*10.0)\n",
    "h_cutoff = lag[int((v_cutoff-0.0075)*1/0.0075)]\n",
    "\n",
    "plt.axvline(v_cutoff, color='g',linestyle='--')\n",
    "plt.axhline(h_cutoff, color='g', linestyle='-.')\n",
    "\n",
    "# Define axis\n",
    "plt.axis([0,0.2,-15,15])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to Uttley et al, the lag-frequency spectrum shows a constant delay until the frequency (1/2*time_delay) which is represented by the green vertical line in the above figure. After this point, the phase warps and the lag becomes negative. This is given in page 43 of review."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### More realistic impulse response"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The response of refelection from an accretion disk to an instantaneous flash follows the _top-hat function_ to first\n",
    "order approximation. The response shows an initial steep rise some time after the initial flash (slope depending on\n",
    "the light travel time to the disk) and then gradually decays, as parts of the accretion disk farther away from the \n",
    "source receieve radiations at later times.\n",
    "\n",
    "The secondary peak is caused due to the bending of light in strong gravitational field around the black hole. This is the re-emergence of photons reflected from the far side of accretion disk that although would be classically blocked from our view, are lensed by strong gravitational field around black hole into our line of sight. \n",
    "\n",
    "Below, we obtain an impulse response similar to one in Utley et al.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 602,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEACAYAAABMEua6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYXHWd7/F3rb3vWzagOxtJSEJCNARlKTYNOsAdhUHc\nkVFU9M59uI6Ic5UeZ0bFZ7zjVRQQB2QcBK7ICFwRFKQQCSQhQpLOnnQHsna6O72numu9f5zqpOh0\nUtXdp+ucqvN5PU89qeV05WNLf/qX71kKRERERERERERERERERERERERERERywgNAO7A5zXbvBqLA\nh6Y8kYiIjMmdwTYPAqvTbOMB7gKeBVyTDSUiIhOTSam/DHSn2ebLwONAx6QTiYjIhGVS6unMBK4F\n7kk+TpjwniIiMgFmlPoPgK9hlLkLjV9ERCzjNeE9VgCPJu/XAlcBEeCp1I3mzJmT2LNnjwl/nYiI\no+wB5ma6sRkr9dlAU/L2OPAFRhU6wJ49e0gkEra/3XnnnZZncErONWtm0dbWbOuMdroppzNzAnPG\nU8iZrNQfAS7BWIXvA+4EfMnX7hvPXyYyIhYbYnj4AJFIun3wIjIemZT6jeN4v5smGkScZXj4LSBB\nNKpSFzGTGeOXvBIIBKyOkJFczxkKtQIuotGerOYZS65/L+1GOa2VzSNVEsn5kAgHDvyYAwfuxuer\nY/nyP1kdR8S2XC4XjKOrtVIXS4RCrZSWrtD4RcRkKnWxxNBQG2Vl52lHqYjJVOpiiVColbIyrdRF\nzKZSl6xLJBIMDbVSUrKURCJMPB62OpJI3lCpS9ZFIl24XF58viq83kpbHAEjki9U6pJ1Q0OtFBY2\nAeD1VmkEI2IilbpkXSjUSlHRbAC83krtLBUxkUpdsm5oqI3CwpFS10pdxEwqdcm6oaHUlXqVZuoi\nJlKpS9aFQq3HV+o+n1bqImZSqUvWnbxSV6mLmEWlLlkVj0cYHj5IQcGZgHaUiphNpS5ZNTz8NgUF\nM3C7jUvya6UuYi6VumRVKNR2/Bh10I5SEbOp1CWrjBOPZh9/rB2lIuZSqUtWpZ54BBq/iJhNpS5Z\nNXqlrh2lIuZSqUtWaaUuMrVU6pJVJ6/UK4jFBkgkYhamEskfKnXJmkikh0Qiis9Xc/w5l8uN11tO\nNNprYTKR/KFSl6wxLuTVNPJBuscZ11TXCEbEDJmU+gNAO7D5FK9/DNgIbAJeAZaaE03yTerlAVJ5\nvVXaWSpikkxK/UFg9WlebwUuxijzfwJ+akIuyUOpF/JKpZ2lIubJpNRfBk73E/cqMDIQXQvMmmwo\nyU+nW6nrrFIRc5g9U78ZeMbk95Q8caqVus4qFTGP18T3uhT4DPDeU23Q3Nx8/H4gECAQCJj414ud\nJRIJBgc3U1y84KTXtKNU5IRgMEgwGJzw17vSbwJAI/A0sOQUry8FnsCYve8+xTaJRCIxrnCSPwYH\nt7F58wc4//zWk45+eeutbxON9jFnznctSidiX8mfl0y72pTxy5kYhf5xTl3o4nDd3c9TWXn5SYUO\n2lEqYqZMxi+PAJcAtcA+4E7Al3ztPuCbQBVwT/K5CLDS3JiS67q7X6C+/oYxX9OOUhHzZFLqN6Z5\n/W+TN5ExxeNRentf4uyz7xvzde0oFTGPziiVKTcw8BcKCmbh9zeM+bp2lIqYR6UuU667+3mqqq44\n5es6o1TEPCp1mXLd3S9QWXn5KV/XjlIR86jUZUrFYiH6+9dRWXnxKbcxxi+96JBXkclTqcuU6utb\nQ0nJErze8lNu43b78HiKiMX6s5hMJD+p1GVKGfP0U49eRmhnqYg5VOoypbq7XzjtTtIR2lkqYg6V\nukyZSKSHY8e2U16+Ku222lkqYg6VukyZnp4g5eUX4HYXpN1WZ5WKmEOlLlMm03k66KxSEbOo1GXK\n9PS8kHGpa0epiDlU6jIlentfIRYLUVq6LKPttaNUxBwqdTFdIpGgtfUfaGz8Ji6XJ6Ov0Y5SEXOo\n1MV03d3PEw4foqHhkxl/jXaUiphDpS6mSiQStLV9naamb+F2Z/5piZqpi5hDpS6m6ux8kng8Ql3d\n9eP6Oh39ImIOMz94WhwukYixd+83aGr6Di7X+NYL2lEqYg6t1MU0R448isdTRk3NB8f9tZqpi5hD\nK3UxRTweoa3tTs4++/4xP1w6nZGjXxKJxIS+XkQMWqnLpMXjYbZu/QilpedSVXXphN7D4ykEXMTj\nIXPDiTiMSl0mJRYboqXlQyQSMRYt+uWk3ks7S0UmT6UuExaLHaOl5Ro8nlLOOedXGV2463S0s1Rk\n8lTqMiGRSA+bNl2F3z+dRYsexu32Tfo9tbNUZPIyKfUHgHZg82m2+SGwC9gILDchl9hUJNJNW1sz\n69bNo6zsPBYseDDjSwGko0sFiExeJqX+ILD6NK9/AJgLzAM+B9xjQi6xmaGh/bS23sHatXMZHt7H\n8uWvMnfuv437ePTT0VmlIpOXySGNLwONp3n9GuCh5P21QCXQgLG6lxwUjw8zNLSXvr7X6Ol5iZ6e\nPxGN9lBffz0rVmygqKhxSv5e7SgVmTwzjlOfCexLebwfmIVK3dZCoTZCoV0MDe096RaJdFFQMIuy\nshVUVl7CrFm3UVKyyNRV+Vi0o1Rk8sw6+Wj02SKJsTZqbm4+fj8QCBAIBEz662U8hocPs27d2VRU\nXExhYSOFhY1UV19FYWEThYWNFBRMN21OPh5+/3T6+9dn/e8VsZNgMEgwGJzw12d66l4j8DSwZIzX\n7gWCwKPJx9uBSzh5pZ5IJMbsesmyY8d2sWnTVaxatdvqKO8wOLiNTZtWs2rVXp1VKpKU/FnI+AfC\njH9PPwWMXDh7FdCDRi+2lkiEcbv9Vsc4SXHxAhKJGKHQLqujiOSsTMYvj2CsvGsxZud3AiMHJd8H\nPINxBMxuYBC4yfyYYqZ4fHjSJwpNBZfLRXX1++ju/gPFxfOtjiOSkzIp9Rsz2OZLkw0i2ROPD+Ny\n2a/UAaqqruTIkceYOfNWq6OI5CSdUepAxvjFrqV+OT09QeLxiNVRRHKSSt2BjPGL/WbqAH5/PUVF\nTfT3r7M6ikhOUqk7kJ3HL2CMYI4e/YPVMURykkrdgey6o3REVdWVdHf/3uoYIjlJpe5Adp6pA1RU\nXMjg4Gai0V6ro4jkHJW6AxnjF3vO1AE8niLKyy+gu/tFq6OI5ByVugPZffwCIyMYzdVFxkul7kB2\nH7+A5uoiE6VSdyC7j18ASkuXEo32EgrttTqKSE5RqTtQLoxfXC43VVVXaAQjMk4qdQdKJOxf6qAR\njMhEqNQdKB63/0wdoLp6Nd3dz+uSASLjoFJ3oFyYqQMUFEynqGguvb0vWx1FJGeo1B0oV8YvADU1\nV9PV9bTVMURyhkrdgXJl/AJGqXd2Po0+NUskMyp1B7L7Bb1SlZYuIx4f4tixHVZHEckJKnUHMsYv\n9p+pg/FpSDU1f0VX1/+zOopITlCpO1AuHKeeqrZWc3WRTKnUHSgeD+fM+AWgsvIyBgbeIBI5anUU\nEdtTqTtQLo1fwLhqY2VlgKNHn7U6iojtqdQdKNfGL6BDG0UypVJ3oFw6+mVETc0HOXr0OZ1dKpKG\nSt2BcuHSu6MVFMygsHA2vb2vWB1FxNYyKfXVwHZgF3D7GK/XAs8CbwItwKfNCidTwxi/5M5MfYSO\nghFJL12pe4C7MYp9EXAjsHDUNl8C3gCWAQHg+4DX1JRiqlwcv8CJubrOLhU5tXSlvhLYDewFIsCj\nwLWjtjkElCfvlwNdQNS8iGK2XBy/AJSWLiceH2ZwsMXqKCK2la7UZwL7Uh7vTz6X6n7gHOAgsBH4\nO9PSyZTI1fGLy+Wiru46Ojp+ZXUUEdtKNybJ5N+5X8eYpweAOcAfgHOB/tEbNjc3H78fCAQIBAKZ\npRRT5er4BaCu7np27LiJxsZ/xOVyWR1HxHTBYJBgMDjhr0/3U7EKaMaYqQPcAcSBu1K2eQb4F2Dk\nsIQXMHaovj7qvRKahdpDMOjl4otDuN0+q6OMWyKR4LXXzmLJkmcoLV1sdRyRKZdcvGS8gkk3fnkd\nmAc0An7gBuCpUdtsB65I3m8AzgZaMw0g2ZVIxIA4Lldu7svWCEbk9NKVehTj6JbngK3AY8A24Jbk\nDeDbwLsw5unPA18FdJEOmzKu++LP6dFFXd31dHQ8bnUMEVvKZLn2u+Qt1X0p9zuBq01LJFMqFy8R\nMFp5+fnEYn0MDm6lpGSR1XFEbEVnlDpMrh7OmMrlclNb+2GNYETGoFJ3mHxYqQPU11/PkSMqdZHR\nVOoOYxzOmHvHqI9WXn4B0WgPg4PbrI4iYisqdYcxrqWe+yt1l8tNXZ1GMCKjqdQdJh7P/Zn6COMo\nGJW6SCqVusPky/gFoKLiPUQiRxkc3Gp1FBHbUKk7TL6MX8AYwTQ0fIz29l9YHUXENlTqDpNP4xeA\nhoZPcvjwL5JnyoqISt1hcvliXmMpLV2M319Pd/eLVkcRsQWVusMY45f8mKmPaGj4JO3t/2F1DBFb\nUKk7TL6cfJSqoeFGOjufIhodsDqKiOVU6g5jXNArv0rd72+gouJCOjufsDqKiOVU6g6Tj+MXgGnT\nPsnhwxrBiKjUHSYfxy8ANTXXMDDwBkND+9JvLJLHVOoOk4/jFwCPp5C6uutob3/Y6igillKpO0w+\nnXw02rRpxlEw+thEcTKVusMY45f8m6kDlJe/h3h8mP7+0R+PK+IcKnWHybeTj1K5XC6mT/8Mhw7d\nb3UUEcuo1B0mHz756HSmTbuZjo5fEY32Wh1FxBIqdYfJ16NfRhQUTKOq6kra2//T6igillCpO0w+\nXXr3VGbM+DwHD96rHabiSCp1h8n3lTpAZeWlxONh+vrWWB1FJOtU6g6T7zN1MHaYzphxCwcP3mt1\nFJGsy6TUVwPbgV3A7afYJgC8AbQAQTOCydRwwvgFYNq0T9HZ+TThcKfVUUSyKl2pe4C7MYp9EXAj\nsHDUNpXAj4GrgcXAdSZnFBM5YfwC4PPVUFt7LYcP/9zqKCJZla7UVwK7gb1ABHgUuHbUNh8Ffg3s\nTz7W0sjGnDB+GTFjxuc5dOg+Eom41VFEsiZdqc8EUq+QtD/5XKp5QDXwIvA68AnT0onpnLJSBygv\nX4XbXUx39wtWRxHJGm+a1zM5JswHnAdcDhQDrwKvYczg36G5ufn4/UAgQCAQyDCmmMUpM3UwdpjO\nnHkrBw78kOrqK62OI5KRYDBIMBic8Ne70ry+CmjGmKkD3AHEgbtStrkdKEpuB/Az4Fng8VHvldBx\nw9Zbv34JCxc+TGnpUqujZEUsFuK11xpZtuwlSkoWWB1HZNxcLhek7+rj0o1fXscYrzQCfuAG4KlR\n2zwJXIixU7UYOB/YmmkAya543DkzdQCPp4gZM77A/v3/ZnUUkaxIV+pR4EvAcxhF/RiwDbgleQPj\ncMdngU3AWuB+VOq25aTxy4iZM79IR8f/JRw+YnUUkSmX8ZLeBBq/2MCaNdNZsWIDBQUzrI6SVTt2\nfA6/fwZNTc1WRxEZF7PHL5JnnDZ+GTFr1m0cPHgPsVjI6igiU0ql7jD5fD310ykpWUB5+Ura2/Xh\n1JLfVOoOY3ycnbNm6iPOOOMr7Nv3v3UykuQ1lbqDJBJxEokYLpfP6iiWqKi4GI+njK6u31odRWTK\nqNQdJB4P43L5R3a8OI7L5eKMM77C229/V9dal7ylUncQY/TivHl6qvr664lEunTpAMlbKnUHMa77\n4sx5+giXy0Nj4zfYu7dZq3XJSyp1BzHGL85eqQPU13+ESKSDnp4/Wh1FxHQqdQfR+MVgrNa/qdW6\n5CWVuoNo/HJCff1HCIePaLUueUel7iBOPfFoLC6Xh7PO0mxd8o9K3UGc9KlHmTixWn/R6igiplGp\nO4iTPvUoE263l7PO+gZtbd/Ual3yhkrdQZx42d10GhpuJBbrpatr9McEiOQmlbqDaKV+MpfLw+zZ\n32PPntuJxyNWxxGZNJW6g2imPrbq6tUUFMzi0KH7rY4iMmkqdQfRSn1sLpeLOXP+lb17v0U02md1\nHJFJUak7iGbqp1ZWtozq6vfz9tt3pd9YxMZU6g6i8cvpNTX9MwcP3svQ0H6ro4hMmErdQTR+Ob3C\nwjOYMeMW2tr+l9VRRCZMpe4gGr+kd+aZX6O7+zn6+tZZHUVkQlTqDqILeqXn9ZYze/Zd7Nz5BRKJ\nmNVxRMZNpe4g8bhm6ploaPgEHk8ZBw7cY3UUkXHLpNRXA9uBXcDtp9nu3UAU+JAJuWQK6IJemXG5\nXMyf/xPeeusfGR4+ZHUckXFJV+oe4G6MYl8E3AgsPMV2dwHPAs78AMwcYIxfNFPPREnJIqZP/1v2\n7PmK1VFExiVdqa8EdgN7gQjwKHDtGNt9GXgc6DAznJhL45fxOeusb9Db+4o+z1RySrpSnwnsS3m8\nP/nc6G2uBUYGkLrcnU1p/DI+Hk8x8+b9kJ07v0g8Pmx1HJGMeNO8nklB/wD4WnJbF6cZvzQ3Nx+/\nHwgECAQCGby9mEXjl/Grrb2Gw4d/zt6932L27H+xOo44QDAYJBgMTvjr082/VwHNGDN1gDuAOMb8\nfERryvvUAseAzwKjr2Wa0DWrrbVly99QV/dh6utvsDpKTgmH21m//lyWLHmS8vLzrY4jDuNyuWAc\n+yrTjV9eB+YBjYAfuIGTy3o20JS8PQ58YYxtxAbi8bDGLxPg9zcwb94P2bbtU8RiIavjiJxWulKP\nAl8CngO2Ao8B24BbkjfJITr5aOLq6/+G0tJluoSA2F42Dz/U+MVib755GWed9Q9UVV1udZScFIl0\nsX79EhYtepTKyoutjiMOYfb4RfKIxi+T4/PVMH/+vWzffhPR6IDVcUTGpFJ3EI1fJq+29hoqKy9m\n9+4vWx1FZEwqdQfRpXfNMXfuj+jtfZXDhx+yOorISVTqDqJL75rD6y3lnHN+xZ49X2FwcKvVcUTe\nQaXuIPrkI/OUli5h9uy72LLlemKxQavjiBynUncQjV/MNW3aTZSVrWDXLs3XxT5U6g6i8Yu5XC4X\n8+b9hL6+Vzl06OdWxxEBVOqOopW6+Yz5+q9pbf0qvb2vWR1HRKXuJJqpT42SkkUsWPAgW7Z8mKGh\nt62OIw6nUneIRCJOIhHB5fJZHSUv1dR8kDPOuI2Wlmu141QspVJ3CONsUv/IKccyBWbNuo3S0uVs\n2/YJEom41XHEoVTqDqHRy9QzPtv0HiKRDl34SyyjUncI7STNDre7gHPOeYKOjsfZv/9HVscRB1Kp\nO4QOZ8wev7+OpUt/z75936O9/WGr44jDqNQdQhfzyq6iokaWLn2W3btvo6vrd1bHEQdRqTtEPK6Z\neraVlJzD4sW/Yfv2T9Hbu8bqOOIQKnWHMMYvKvVsq6i4gIULf0FLy1/T17fe6jjiACp1hzDGL5qp\nW6G6+v2cffbP2Lz5gzrrVKacSt0hNH6xVm3t1SxY8BAtLdfQ2/uK1XEkj6nUHULjF+vV1FzFwoX/\nSUvLX9PT85LVcSRPqdQdQke/2EN19ftYtOhRtmy5nq6uZ6yOI3lIpe4QxslHmqnbQVXVZSxe/BTb\nt3+GQ4f+3eo4kmdU6g5hXPtFK3W7qKhYxfLlf+Ktt75NW1sziUTC6kiSJzIt9dXAdmAXcPsYr38M\n2AhsAl4BlpqSTkyj8Yv9FBfP57zz1nD06G/ZseNm4vGI1ZEkD2RS6h7gboxiXwTcCCwctU0rcDFG\nmf8T8FMTM4oJNH6xJ7+/gXPPfZFw+AibNr2PcLjD6kiS4zIp9ZXAbmAvEAEeBa4dtc2rQG/y/lpg\nlkn5xCQav9iX11vKkiVPUl5+ARs2vIv+/g1WR5IclkmpzwT2pTzen3zuVG4GtFvfZjR+sTeXy8Ps\n2d9mzpzvs2nTag4f/g+rI0mO8mawzXj24FwKfAZ471gvNjc3H78fCAQIBALjeGuZDF16NzfU119H\nSclCWlr+G319rzFnzvfxeIqsjiVZFAwGCQaDE/76TD4GZxXQjDFTB7gDiAN3jdpuKfBEcrvdY7xP\nQnv4rdPWdifgoqmp2eookoFIpIedOz/PsWNbWLjwEUpLF1sdSSyS/LSyjD+yLJPxy+vAPKAR8AM3\nAE+N2uZMjEL/OGMXulhMn3yUW3y+ShYteoRZs/4nGzdeyoEDP9Fhj5KRTEo9CnwJeA7YCjwGbANu\nSd4AvglUAfcAbwDrTE8qk6LxS+5xuVxMn/5pli9/hUOHHqCl5RqGhw9YHUtsLpufQqzxi4V27ryV\nkpJFzJx5q9VRZALi8TBvvfVtDh78MU1N32H69Jv1IeIOMRXjF8kDiURYH2eXw9xuP01NzZx77gsc\nPHgvGzdeSSjUanUssSGVukNo/JIfSkuXct55r1Fd/X42bHg3e/d+i1gsZHUssRGVukOo1POH2+3l\nzDP/nhUrNjAwsJH168+ho+M32pEqgErdMYzrqWv8kk+KihpZvPjXzJ//U9ravs6mTe9nYGCj1bHE\nYip1h9AhjfmruvoK3vWujdTUXM3Gje9n69aPEQrtsTqWWESl7hAav+Q3t9vHrFlf5vzzd1FcvIAN\nG85n585bdQikA6nUHUIfZ+cMXm8ZjY3fYOXK7bjdRaxfv4QdO27RkTIOolJ3COOCXpqpO4XfX8vc\nuf/KypU78Pnq2LBhJdu2fYKBgc1WR5MpplJ3iHhcM3Un8vvrmD37n1m1ag/FxYvYtGk1b755GZ2d\nT5JIxKyOJ1NAZ5Q6xNq1C1i8+DeUlCywOopYKB4P09HxOPv3/x8ikU5mzPgC06Z9Cr+/zupocgo6\no1TGpPGLgHFmakPDR1mxYi0LFz7M4GALa9fOY8uW6zl69Dmt3vNAJtdTlzyg8YuMVlGxioqKVUSj\nvbS3P0Jr6z8QiXyW+vobaWj4GCUlS3R9mRyk8YtD/PnPtaxcuR2/v9bqKGJjAwObOXLkl7S3/xKP\np4yGho9SV3cdxcXzrY7mWOMdv6jUHeLll8u44IKDeL1lVkeRHJBIxOntXcORI4/Q2flfeL1V1NV9\niNraD1Faukwr+CxSqcuYXnrJz0UX9WsEI+OWSMTp61tHZ+cTdHQ8QTw+RE3NVVRXX0VV1RV4veVW\nR8xrKnU5SSKR4KWX3FxySVwrLJmURCJBKLSLo0d/R1fXM/T1raG0dDmVlZdRVXU55eXna4e8yVTq\ncpJ4fJiXXy7nkkuGrY4ieSYWG6S39890d/+R7u4XCIV2Ul6+ioqKC6mouIjy8pV4PCVWx8xpKnU5\nSTTax6uvzuKii/qsjiJ5LhI5Sm/vK/T2/pne3j8zMLCRkpKFlJWdT3n5SsrKVlJcPB+XS0dTZ0ql\nLicJhztZt24BF17YaXUUcZhYLER//wb6+9fT17eW/v51RCJdlJYuo7R0OWVlyyktXU5x8QKNbU5B\npS4nGR4+wIYNK3nPe3TFPrFeONzJwMAbDAy8mfzzDYaG9lJYOIeSksXJ2yKKixdQVDTX8WU/3lLX\nyUcOoMvuip34/bVUV19JdfWVx5+LxYY4dmw7g4MtDA5u5vDhhzh2bDtDQ29RWHgmRUXzKSqam3Kb\nQ2HhmfrvegwqdQfQpx6J3Xk8hZSVLaOsbNk7no/Hw4RCewiFdhEK7ebYsa10dT1FKLSH4eH9+Hx1\nFBU1UVjYSEHBGRQUnElh4RnJ+zPxeqsdd8RXJqW+GvgB4AF+Btw1xjY/BK4CjgGfBt4wKZ+YQJ96\nJLnK7fZTUrKQkpKFJ70Wj0cJhw8wNLQ3edvHwMBf6Op6kqGhtwmHDxKLhSgomIHfP4OCgun4/dOS\nt+n4fPX4/fX4fHX4/fV5c5ROulL3AHcDVwAHgPXAU8C2lG0+AMwF5gHnA/cAq0xPmiXBYJBAIGB1\njLTGk9PK8UsufD9zISMo52hut5fCwrMoLDwLuGTMbWKxYwwPHyQcPkA4fPj4rafnT6xZs41ly+KE\nw0eIRDoA8Plq8Plq8flq8Xpr8Pmq8Hqr8fmq8Xqr8HorU24VeDwVeL0VuN2+Kf/fm6l0pb4S2A3s\nTT5+FLiWd5b6NcBDyftrgUqgAWg3LWUW5eMPjkr99HIhIyjnRHg8xRQXz6W4eO5Jrz32WDM339wM\nGCdVxePHiES6iEQ6j/8ZjXYTiRxlePggg4MtRKM9o259RKO9uN0+PJ5yvN5yPJ6y4zevtwyPpxS3\nuwSPpxSPpwSPpwS3uzh5vxi3uxi3uyh5vyh5vwi3uxC3u2jc/5vTlfpMYF/K4/0Yq/F028wiR0s9\nHyUSYc3URU7D5XIdL9zCwjPH9bXGL4QQ0WgvsVg/sVg/0Wh/8v4Asdhg8s8BotFewuFDyecGiccH\nicVCxOMh4vFjxGLHiMeHkjfj+fFKV+qZHoM4ek/EmF+3efPVGb6dddrbd7B58warY6Q1npzh8BF8\nPl2dUWQqGL8QivF4ioHpU/E3mLr1KqAZY2cpwB1AnHfuLL0XCGKMZgC2Ywy4Rq/UdwNzxpVORET2\nYOy3NIU3+YaNgB94Exi9G/oDwDPJ+6uA18z6y0VExHxXATswVtp3JJ+7JXkbcXfy9Y3AeVlNJyIi\nIiIiE7MaY86+C7jd4iypHsCY+29Oea4a+AOwE/g9xuGZVjsDeBHYArQA/z35vN2yFmIc0vomsBX4\nTvJ5u+Uc4cE4Se7p5GM75twLbMLIuS75nN1yVgKPYxzmvBXj6Di7ZTwb43s4cuvF+DmyW04wpiFb\nMHrpl0ABNsvpwRjLNAI+xp7JW+UiYDnvLPXvAV9N3r8d+G62Q41hGjBy7nQpxihsIfbMWpz804ux\nb+VC7JkT4DbgYYyT6cCeOdswfqBT2S3nQ8Bnkve9QAX2y5jKDRzCWCzZLWcj0IpR5ACPAZ/CZjkv\nAJ5Nefy15M0uGnlnqW/HOHEKjDLdnu1AGfgNxhm+ds5ajHH28TnYM+cs4HngUk6s1O2Ysw2oGfWc\nnXJWYJTQaHbKONr7gJeT9+2Wsxpj0VaF8QvyaeBKbJbzOuD+lMcfB35kUZaxNPLOUu9Oue8a9dgO\nGoG3gDLsmdWN8a+xfozVBdgz568w/pV2CSdK3Y45WzHGBa8Dn00+Z6ecyzBGbg8Cf8H4WS/BXhlH\newD4YvKEBcTAAAAB+ElEQVS+HXN+DuPn5wjwi+Rz48o51R8/kssXUE9gr/ylwK+Bv8P4Pz2VXbLG\nMX7QZwEXY6yEU9kh519h/MC8wanP07BDToD3YvzyuQq4FWNkmMrqnF6Mo91+kvxzkJP/JW51xlR+\n4GqMX+qj2SHnHOB/YCzeZmD8zH981DZpc051qR/AmF2NOAPjMgJ21Y7xzxswTg07YmGWVD6MQv8F\nxvgF7JsVjB1RvwVWYL+c78G4XlEb8AhwGcb31W45wZj9AnQA/4VxLSY75dyfvK1PPn4co9wPY5+M\nqa4CNmB8P8Fe30uAdwFrgC4gCjyBMcIe1/dzqkv9dYyrNzZi/Ja8gRM7puzoKYwdEyT//M1pts0W\nF/DvGEcW/CDlebtlreXEXvkijFngG9gv59cxFhdNwEeAPwKfwH45izHGbGCMNN6HMSq0U87DGNd9\nmp98fAXGkRtPY5+MqW7E+EU+wk7fSzBm5aswfn5cGN/Prdjw+znWyUt28AhwEAhj/Id5E8aOiuex\nyaFDSRdijDXe5MQhWauxX9YlGHPVNzEOw/v75PN2y5nqEk4sMuyWswnje/kmxqGsIz87dst5LsZK\nfSPGyrIC+2UE4xdjJyd+UYI9c36VE4c0PoTxr3Q75hQRERERERERERERERERERERERERERERERGR\nXPP/AVSAH7dssJeAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x113fc6250>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Primary peak time, secondary peak time, end time\n",
    "t1, t2, t3 = 3, 4, 10\n",
    "# Peaks' values\n",
    "p1, p2 = 1, 1.4\n",
    "# Rise and decay slopes\n",
    "rise, decay = 0.6, 0.1\n",
    "\n",
    "# Append zeros before start time\n",
    "h_primary = np.append(np.zeros(int(t1/lc.dt)), p1)\n",
    "\n",
    "# Create a rising exponential of user-provided slope that ends at secondary peak time and secondary peak\n",
    "# value\n",
    "x = np.linspace(t1/lc.dt, t2/lc.dt, (t2-t1)/lc.dt)\n",
    "h_rise = np.exp(rise*x)\n",
    "# Find a factor for scaling\n",
    "factor = np.max(h_rise)/(p2-p1)\n",
    "h_secondary = (h_rise/factor) + p1\n",
    "\n",
    "# Create a decaying exponential until the end time\n",
    "x = np.linspace(t2/lc.dt, t3/lc.dt, (t3-t2)/lc.dt)\n",
    "h_decay = (np.exp((-decay)*(x-4/lc.dt))) \n",
    "\n",
    "# Add the three responses\n",
    "h = np.append(h_primary, h_secondary)\n",
    "h = np.append(h, h_decay)\n",
    "\n",
    "# Plot\n",
    "plt.plot(h,'y')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Obtain output through convolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 603,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "delay = (int(t3/lc.dt))\n",
    "output = signal.fftconvolve(s, h)\n",
    "output = output[delay:-delay]\n",
    "s_mod = s[delay:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Form light curves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 604,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "time = lc.time[delay:]\n",
    "lc1 = Lightcurve(time, s_mod)\n",
    "lc2 = Lightcurve(time, output)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find cross spectrum and compute lags."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 605,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cross = Crossspectrum(lc1, lc2)\n",
    "cross = cross.rebin(0.0075)\n",
    "lag = np.angle(cross.cs)/ (2 * np.pi * cross.freq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 606,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFG9JREFUeJzt3WuQFOW9x/HvcodFRERAFIPA4gVD8BKEYMJ6gpFYSYxa\nqVRC5Zgck8qlkqqkzgtNzovtzUmqkipNpRLr5EX0nMqpJGoqaoJJvEB0MRFEUUAMcICNKAhBCrkE\ndGFh57x4ZtgLu7Az07M9Pf39VHXNTM/09OPj8Jvep//9DEiSJEmSJEmSJEmSJEmSJCnl/hvYA2zo\nsm4csAzYAjwFjE2gXZKkXgyK4T3+B1jcY91dhOCfCfw5/1iSVEOm0v2IfzMwMX9/Uv6xJKkKxHHE\n35uJhOEf8rcTT/NaSdIAqlTwd5XLL5KkKjCkQu+7hzDE8w/gfOCtni+YPn16rrW1tUK7l6Sa1QrM\nKOcNKnXEvxS4PX//duB3PV/Q2tpKLpdziWlpampKvA21tKSlP5ueSUk7U9KfaViA6eUGdBzB/wCw\nErgE2AF8AfgBcAOhnPNf8o8lxax5RXPSTVAKxTHU85k+1i+K4b0lSTEbiJO7GgCNjY1JN6Gm2J/x\nsj+rS12C+87lx6sklaiuuY5ck/+OsqSurg7KzG6P+CUpYwx+KcWaFjYl3QSlkEM9kpQiDvVIkopm\n8EtSxhj8kpQxBr8kZYzBL6VY1BIl3QSlkFU9Uop5AVf2WNUjSSqawS9JGWPwS1LGGPySlDEGv5Ri\nztWjUljVI0kpYlWPJKloBr8kZYzBL0kZY/BLUsYY/FKKOVePSmFVj5RiztWTPVb1SJKKZvBLUsYY\n/JKUMQa/JGWMwS+lmHP1qBRW9UhSiljVI0kqmsEvSRlj8EtSxhj8kpQxBr+UYs7Vo1JY1SOlmHP1\nZI9VPZKkohn8kpQxBr8kZYzBL0kZY/BLKeZcPSqFVT2SlCJW9UiSimbwS1LGGPySlDEGvyRljMEv\npZhz9agUVvVIKeZcPdkTR1XPkHia0qftwCHgBNAOzK3w/iRJZ1Dp4M8BjcDbFd6PJKmfBmKMP8nh\nJElSD5UO/hywHFgDfKnC+5Ik9UOlh3oWALuB84BlwGbgL4Unoyg6+cLGxkYaGxsr3ByptjhXT+1r\naWmhpaUl1vccyGGYJuAwcE/+sVU9klSkap+rZxRwVv5+PfARYEMF9ydJ6odKDvVMBB7tsp9fAU9V\ncH+SpH7wAi5JSpFqH+qRJFUhg19KMefqUSkc6pFSzLl6ssehHklS0Sp9AZfSKJeDd96Bw4dPXY4c\n6f64vR2GDIGhQ8NtX0vh+cGDwz46Ovq/nDjRuRw/3vv9ro87OqC+HsaOhXPO6bzten/kSKiL8Q/e\njo7QHwcPwqFDYSncL9y2tcH48TBxIkyYEG4nTgxtjbMt0hkkOtRDFO40LWwiaoy6PVkYu+xtffOK\nZrcrrD9+HPbsIVrRTPPWn4ftch8i2nVJCJv2djh2jGjiJug4QbT1wrCusP6S3TS/d1/Y7sVRRCsG\nhdAfMQJGjyZacAyGDSPaMSME1OjRYf3ETTSPXhO2a/8A0TvXhrbkl6j+RejoINo7K+yrsH7SJpov\nbA3b7ZpJtOcyGDTo5BJN2Ah1dUT739d9/TnraB6zNmz37lyi9uvCl0h+iYb8NWw39IawzZEjcOAA\nUd0Kms8Nl480vTqe6LmhsH9/+II45xyiDx6HESOJ/nFp6M+ODsjliN7zGs0Xvx62a51C1Drl5HPk\nckQzdkJ7O9Gq4aGfDx+GUaOIrq+j+ep/hu12TCM6dDWMGROWESOIjj4JRw4TrR8He/aEJZcjunEY\nzXMOhu0OXUk04qPhS+GWW2DKlD4/L3XNnf+Eq/Lz6XaxbxfHUI9j/AOpowO2bIHVq0P4FI7y6ur6\nXgrP5wOeXbtg9+7O2337wlHk+efD5MnhtrCMHRuOtIcN6/22t3WjR8OoUSE8a1lbGxw4EJb9+0Nw\nd+33QYNO/7iuLnwRFkL9rLM6/5op1pEjnV8Cb73Vef+pp2DhQvje9/rc1DH+7EnDfPzZdvAgvPAC\nrFoVltWr4eyz4dprYdKkk0ePfS7QeX/w4LDN3LndQ37ChDCEouKMGBH6c9KkpFsSvkCmTQtLV9Om\nwdKlp93UuXpUCo/4e9q6NRz9FYY16uvDcqZw7eiAzZs7Q/7552H7drjqKpg/Pyzz5lVH0CgdXngB\nvvIVePnlpFuiKuJQT1x27oQHHoBf/Qr27oXzzus8iXnkSFiGDu3+ZdD1tq0NXnwRzj23M+Dnz4fZ\ns8N2Uin274eLLgrnEDz5qzyDvxwHDsDDD8MvfwmvvAK33gpLlsCHPnTq+HYuF8K9a1VL4Qvh8OEw\nDHPNNeFknBSn8ePh1Vf9S1EnOcZfrKNH4U9/CmG/fDksWgTf+AbcdFMY8+1LXV0o/xs5Mvw1IA2U\nhoYw/GjwK0a1H/wdHfDss2EY55FHwvDLkiVw332hnluqZoXg/+AHk26JakhtBv/+/dDSEo7qH3ss\nBPySJbBuHUyZknTrpP4rBH8fopbolNpv6UxqY4z/6FFYuTIE/fLlsHEjLFgQhnIWL4YrrohnP9JA\ne/BB+O1vw9IL6/izJ7tj/B0dsH59Z9CvXAmzZoWg/+EPQ0XN8OFJt1Iq3xmO+KVSpCv4H3kEHnoI\nnn46lE4uWgRf/WpYN3Zs0q2T4tfQANu2hcoySzoVk3QF/7e+FZa773asXtkwZky4VmTXLrjggqRb\noxqRrglZTpyA224z9JUtDvcoZukL/lInwpLS6jTB71w9KoXBL1W7mTP7DH5LOVUKg1+qdg71KGYG\nv1TtDH7FzOCXqt2MGdDaGq5fkWJg8EvVrr4exo2DHTuSbolqhMEvpUEfwz2F32SVimHwS2nQR/AX\nfqhbKobBL6WBJ3gVo/QEf+HEVs9fx5KywOBXjNKToh7tK8sMfsXI4JfSYPp02L4djh9PuiWqAQa/\nlAYjR8KECfDGG91WO1ePSmHwS2nRy3CPc/WoFAa/lBaO8ysmBr+UFga/YmLwS2lh8CsmBr+UFga/\nYmLwS2kxbVqo6mlvP7nKuXpUCoNfSovhw2Hy5FDPn+dcPSqFwS+licM9ioHBL6WJwa8YGPxSmhj8\nioHBL6WJwa8YGPxSmvQIfufqUSkMfilNLr4Y3nwTjh0DnKtHpTH4pTQZOhSmTIG//z3plijFDH4p\nbRznV5kqGfyLgc3AVuDOst/N4JcCg19lqlTwDwbuJYT/5cBngMvKekeDXwoMfpWpUsE/F9gGbAfa\ngQeBm8t6R4NfCroEv3P1qBSVCv4LgB1dHu/MryudwS8FXYLfuXpUiiEVet9cf14URdHJ+42NjTQ2\nNvb9YoNfCt7zHtizB9rakm6JBkBLSwstLS2xvmelgv9NYEqXx1MIR/3ddA3+MzL4pWDIkBD+ra1J\nt0QDoOdBcXNz+X/lVWqoZw3QAEwFhgGfBpaW9Y4Gv9SpoQG2bEm6FUqpSh3xHwe+DjxJqPC5H9hU\n1jsa/FInK3tUhkoFP8Dj+SUeBr/UqaEB1q6laYlz9ah4XrkrpdHMmbB1q3P1qCQGv5RGDvWoDAa/\nlEZTpsDbb8ORI0m3RClk8EtpNGgQTJsG27Yl3RKlkMEvpZXDPSqRwS+lVUMD0aafJd0KpZDBL6VV\nQwPNHU8n3QoNhOPHYft2iGnqhkrW8cfL4Je6a2iA3Uk3QrE5dgw2boR16+C110LQv/56uN29GyZO\nDFN1xMDgl9KqoQGeTboRKsm778KGDfDyy53Lxo3hN5XnzIEZM2DhQpg6NSwXXgjDhoVt6+rK3r3B\nL6XV5Mnh9tAhGDMm2baob4cPw/r1nQH/0kuhGuuSS+Cqq8LyhS/A7NlQXz8gTTL4pbQalD9Ft21b\nCA8l78ABWLu2+5H8G2/ArFlw9dXwgQ/A178OV1wBw4cn1kyDX0qxpj2XhpJOg3/g7d3bPeRfegne\negve977w/+OGG+DOO+Gyy2Do0KRb243BL6VYNOZma/krLZeDN98MId816A8d6hyq+eQn4bvfDedd\nUpBT6Qn+WbPg6NGkWyFVl4YGeNYzvLE5cSJ8kRZCfu3aUGUzaBBceWVYliyBe+4JV07HcKI1CUm2\nOpfL9esXGiX15dln4a67YOXKpFuSPm1t8Oqr3UN+w4ZQNlkI+Tlzwu3551dNyNeFdpTVGINfSrPd\nu0M1yN69Sbekuh07FkJ+zZrOZfPm8BdT14CfMwfOPjvp1p6WwS9lXS4XSjl37ICxY5NuTXVobw81\n8WvWhBOua9aE0J8+Ha65pnOZPRtGjky6tUWLI/jTM8Yv6RTRimaiGTPCuPT73590c5Kxbx888wys\nWBFC/pVX4KKLOgN+yZJwJD9ANfJp4BG/lGJ1zXXk/vapUFXy2c8m3ZyB0dYGzz0Hy5fDsmXhR+ev\nuw6uvx7mzg1DNjV8QZtH/JJqf3rmjo5QWbN8eVhWrQoXQC1aBD/6Ecyb1zmdgfrF4JfSrqEhBGKt\nyOXCJGV//nP473r6aTj33BD0X/sa/OY3ns8ok8EvpV1DA/wsxfPyHz4cxuZXrYLnnw/L4MFh6Gbx\nYrj77vBTk4qNwS+lXZqGenK50NauIb9lS6iwmT8/nIj96U9D0FdJ3XwtMvilFGta2ATnnReuON23\nLwyJVJNcDl54AZ56KoT96tVw1llhXH7ePPj850PFTYITlmWRVT1SLbjmGrj33hCmSTtxAv76V3jk\nkbDU18PHPx5mprz22s7ppFUSq3okBYXhnqSCv7091NI//DD87nch3G+9FZ58Ei6/PJk2qU8Gv1QL\nkhjnb2sLQzgPPwx/+ENow223hXmDpk8f2LaoKAa/VAsaGuCPf6z8fg4dgieeCEM4TzwRxudvuw2+\n//3w84BKBYNfqgWVPOLfsQMeewyWLg1H8wsWhGGcn/wEJkyozD5VUZ7clVIsaomIGqNQ0XPxxXDw\nYPllkLlc+I3Y3/8+hP3rr8NNN8EnPgE33hiqcpQYZ+eUMq6uuY5cU/7f0bhxsGlTmE++WMeOhUnO\nli4Ny9ChcPPNIewXLIAhDg5UC6t6JHUqDPf0N/jb28NJ2YceCtU3l14agv7xx8PvxHoBVc0y+KVa\nUQj+6647/eu2bIH774df/AJmzoTPfQ5+/GOYNGlg2qnEGfxSrTjdCd533glll/fdF3556vbbw9DO\nJZcMbBtVFQx+qVY0NMCjj3Zf9/LLIewfeihc3PXNb8LHPhbG8JVZBr+UYk0LmzofFI74DxyAX/86\nBP7bb8Mdd4T57J3hUnlW9Ui14sCBMGFbfX0ou/ziF+HDH4ZBg5JumWJkOaek7h5/PPz27vjxSbdE\nFWLwS1LGxBH8/g0oSRlj8EtSxhj8UopFLVHSTVAKOcYvpVi3uXqUCY7xS5KKZvBLUsYY/JKUMQa/\nJGVMpYI/AnYCa/PL4grtR8q0bnP1SP1UqaqeJuCfwI9O8xqreiSpSNVe1ePP90hSFapk8H8DWA/c\nD4yt4H4kSUUoZz7+ZUBvv9X2H8DPgO/mH/8ncA9wR88XRlF08n5jYyONjY1lNEeSak9LSwstLS2x\nvudADMdMBR4D3ttjvWP8klSkah7jP7/L/VuADRXaj5RpztWjUlTqiP9/gTlADngN+DKwp8drPOKX\nyuRcPdkTxxF/pX5z918r9L6SpDJ55a4kZYzBL0kZY/BLUsYY/FKKOVePSuEvcElSilRzHb8kqUoZ\n/JKUMQa/JGWMwS9JGWPwSynmXD0qhVU9Uoo5V0/2WNUjSSqawS9JGWPwS1LGGPySlDEGv5RiztWj\nUljVI0kpYlWPJKloBr8kZYzBL0kZY/BLUsYY/FKKOVePSmFVj5RiztWTPVb1SJKKZvBLUsYY/JKU\nMQa/JGWMwS+lmHP1qBRW9UhSiljVI0kqmsEvSRlj8EtSxhj8kpQxBr+UYs7Vo1JY1SOlmHP1ZI9V\nPZKkohn8kpQxBr8kZYzBL0kZY/BLKeZcPSqFVT2SlCJW9UiSimbwS1LGGPySlDEGvyRljMEvpZhz\n9agU5QT/p4C/ASeAq3o8921gK7AZ+EgZ+5B0Gs0rmpNuglKonODfANwCPNtj/eXAp/O3i4H/KnM/\n6oeWlpakm1BT7M942Z/VpZxA3gxs6WX9zcADQDuwHdgGzC1jP+oH/2HFy/6Ml/1ZXSpxJD4Z2Nnl\n8U7gggrsR5JUgiFneH4ZMKmX9d8BHitiP16iK0lVIo4pG54B/h14Of/4rvztD/K3TwBNwOoe220D\npsewf0nKklZgRtKNeAa4usvjy4F1wDDgYkIjk5wTSJIUk1uAHcC7wD+Ax7s89x3CEf1m4MaBb5ok\nSZKkilpMONrfCtzZx2t+kn9+PXBlkdtmTTn9uR14BVgLvFC5JqbGmfryUmAV0EY4d1XMtllUTn9u\nx89mT2fqzyWEf+OvAM8Bs4vYtqIGE4Z5pgJDCeP9l/V4zU3An/L3rwWeL2LbrCmnPwFeA8ZVtomp\n0Z++PA+4Bvge3YPKz+apyulP8LPZU3/6cz5wdv7+YkrMzkrU8c/NN2A74SKuBwkXdXX1CeAX+fur\ngbGEstH+bJs1pfbnxC7Pe3I96E9f7gXW5J8vdtusKac/C/xsdupPf64CDubvrwYuLGLbkyoR/BcQ\nTvoW9HYBV1+vmdyPbbOmnP6EcA3FcsI/vi9VqI1p0Z++rMS2tarcPvGz2V2x/XkHnX/pF7XtmS7g\nKkV/L9bym75/yu3P64BdhD+5lxHGAP8SQ7vSqJwLCb0I8VTl9skCYDd+NguK6c/rgX8j9GGx21bk\niP9NYEqXx1PoPoVDb6+5MP+a/mybNaX255v5+7vyt3uBR8n2vEnlfL78bJ6q3D7Znb/1sxn0tz9n\nAz8nDPHuL3LbihlCuGhrKuEirjOdjJxH5wmK/mybNeX05yjgrPz9ekIVQJanyS7m8xXR/WSkn81T\nldOffjZP1Z/+vIgwlj+vhG0r7qPA/xEa+O38ui/nl4J788+vp/t8/r1tm3Wl9uc0wgdgHfAq9iec\nuS8nEcZKDxKOpt4ARp9m26wrtT/9bPbuTP15H7CPUALbswzWz6ckSZIkSZIkSZIkSZIkSZIkSZIk\nJe3/ATpfKjoz68jQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x113fa1c50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "\n",
    "# Plot lag-frequency spectrum.\n",
    "plt.plot(cross.freq, lag, 'r')\n",
    "\n",
    "# Define the x-position of vertical line\n",
    "v_cutoff = 1.0/(2*t2)\n",
    "h_cutoff = lag[int((v_cutoff-0.0075)*1/0.0075)]\n",
    "\n",
    "plt.axvline(v_cutoff, color='g', linestyle='--')\n",
    "plt.axhline(h_cutoff, color='g', linestyle='-.')\n",
    "\n",
    "# Define axis\n",
    "plt.axis([0,0.2,-10,10])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Energy Dependence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### With same intensity and varying position"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create different lags for different energy channels, we create delta impulses of same intensity at different positions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 607,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "energies = np.array([4.5,8.5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create impulse responses for all energy channels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 608,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "h_zeros = [np.zeros(int(i/lc.dt)) for i in energies]\n",
    "responses = [np.append(h, 1) for h in h_zeros]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 609,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "delays = [int(i/lc.dt) for i in energies]\n",
    "outputs = [signal.fftconvolve(s, h)[d:-d] for h,d in zip(responses,delays)]\n",
    "s_mods = [s[d:] for d in delays]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make light curves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 610,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "t_mods = [lc.time[d:] for d in delays]\n",
    "lc_input = [Lightcurve(t_mod, s_mod) for t_mod, s_mod in zip(t_mods,s_mods)]\n",
    "lc_output = [Lightcurve(t_mod, output) for t_mod, output in zip(t_mods,outputs)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 611,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cross_spectrums = [Crossspectrum(lc1, lc2).rebin(0.0075) for lc1,lc2 in zip(lc_input,lc_output)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute lags and cutoffs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 612,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lags = [np.angle(cross.cs)/ (2 * np.pi * cross.freq) for cross in cross_spectrums]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get cutoff points for all energy channels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 613,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "v_cutoffs = [1.0/(2*energy) for energy in energies]\n",
    "h_cutoffs = [lag[int((v_cutoff-0.0075)*1/0.0075)] for lag, v_cutoff in zip(lags, v_cutoffs)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We plot lag-frequency spectrum for all energy channels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 614,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEZCAYAAABrUHmEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8E3X++PFX2kJb2gKlHD2hXG25DwVBQCqCoIjKzwtE\nBWVd8HZlXV39KmVR111ldRePBW9U8FxXQUUQCAJrQYVyl7MttNzQQimlRzK/Pz6TJmnTNm2SJmnf\nz8djHklmMjOfmUzmPZ9jPgNCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghGpE04ENvJ6KJSUP2\nudAFeDsBjUw2cAEotBn+5c0EOckMnEel9xTwI3CrV1NUPc1Ny0lEbXdN/4E0oAz73/OPblq/P3HX\nPndkOrAbOAccA74Fwj24vkRq/91FDYK8nYBGRgOuA1Z7eD2BgMnNy+wLHATaANcCrwEpwF/cvB5f\nY6hhmgYsAe6qZRkBqBNRY1XTPnLFSOB5YCywFYhE/X8aQk3b5In/V6Mh0bbhTAPWAy8BZ1An6HE2\n01sB7wBHgFxgLtbfZxqwAfgHKicwG3VyXwqcBTYBzwHr9O+/Drxcaf3fAI86kc4zwEfAfcCf9fU4\nm775QAHqynFUHbatpv3SGViLuhJdAbStlN4hwP+AfCADdSKyMKKC3np9/h+AKH3aT/prASoHcZmD\nfWHA8cnlfeBN4DtUDi0ViAW+BE7o2/CQzfdD9XnOADuBx4HDNtPNQJdKy59r8/k6fdvyUfu5j820\nbGAW6qRbAHwCBNtMv0Gf9yywH3WCvgX4tdI2PQb818G2OvI5cFRf31qgp820KKo/LisbBPyspx3U\n9n2I2qeg9sO/Ub/7OdTv2dFm/hRgJXAayNS3yyIUmIfaPwWo3zsE+9/9HOr4mYb9/yuNqkVyidjn\nUIyo32gD6vj5BnVsfmyz7Z2q2W4hKmQBV1UzbRpQisqOG4CZQJ7N9K9QJ6JQoB2wEfi9zbxlwAOo\ngzYEdXJYrL/vARzC+ocYpC/bcsJrCxTpy3Wk8kkLoJm+zrF1SN8jqKu0W1F/ytZOzlvTfvkZFQCb\nASNQf/RF+rQ41J/cEmRG658tgcEI7AO6ofbTGuCv+rROOFc85ags/319+4bqn0OB34D/Q+XeOwMH\ngKv16S+iTq6tgXhgB+r3sqi8/9/DmsMbABxH/aYGVK4nC7U/0N+nA9GoK/VdwAx92mA9nZZjMhZI\nBpqjTrQpNuvcAkx0tBOouh+mAWF6Gl7R57Wo6bisbDiqODcNGIZ9sAO1n8/p32sOvIo1AIWhAu9U\n1G/YHziprxPUhdNqIEafPkRfhqPffRpV/1+zqT1o7EX91i1RFwP7UBdLgcAHwLvVbLcQFbJRVx35\nNsN0fdo01EFl0QJ1ELYHOgAXUQerxWSsxVzTgBybaYGoE213m3Fzsb+i24U6iQI8CCyrId2Oggao\nq8nJTqYvD3sbgTucnLe6/dIR9WcOtZn+Mdag8YTNe4vlWIuT1gBP2Uy7D/hef5+Ic0GjBOtveQZ1\nEnpfHywuw/73AZVLs5w0bAMIwL3UnNOwDRpvUrWIMBMVQEEFjdttpv1NnwdgAepq25E3UbkAgF6o\nbWtWzXfTqL4ivLWe/gicOy4rG4e6Ss9H/XfmYf1N3kcFIIswoBwVeG+jajBaADyrz38B+xyZRSKO\ng0bl3y+NmoPGGtRvbPEyqj7G4jrsg2mjIXUa7qWhigOqq9M4ZvP+gv4ajsoJNEOdpC0CsL8atT3J\ntEP9drbjciutaxHqpP2j/vpK7cm300xfzxnU1Vlt6ascNHJQV7YdnZi3uv3SHnUyKa603AT9fSdU\nkcQEm+lB2O9/22UXU/dK1k+pWqehYb+9nVDbmm8zLhDrSS0W+9/Kdttr00lfv21xVzN9mRaVtzFG\nfx+P/YnM1geoE/L/AXeitrMMmIIqEkJP//hK8wWi6iFuRh0fZtT+aKuvu7bjsrLl+gDqKv1zYA+w\nUF+u7fxFqOMxFrVfLsN+nwehjvso1EXKgVrWbetw7V+p4rjN+4uooknbz56s0PcaCRq+4TDqijaK\n6itUbVuwnERdcSVgvUpPqPT9j4DtQD9UMYSz5dUWN+jr2IT6A9aWvrhKnzsBX+PctlXnKKrIpQXW\nYNIJayXlIdTV4O+rzlorZ1oEaVRfYWo7/yHUFX9SNd89igqeu/XPHStNv4DaRosYrCexQ6iT9AtO\npLeyw6iiOUfSUbmCK1A5v8n6+I/1oTq3A9ejirxyUDmNM6j95MxxWZPV+tBL/2yoNH84qo4tD7Vf\n1mKfg7MIQJ20uwHbKk2r7nevPP489r9JdC1p92QLM58iFeHuV5+WJkdRlX3/QGXzA4CuqD+0Iybg\nP6gsdCgqKNyJ/YGbi6rsXAR8gTpxO5PuNqirzddQZfH5TqavPfAw6ir4Fj1N36GuguuybbZy9G2Y\noy93OPataz5C5TKuRl0Bh6AqpW0DWHW/x0lUEOtaw/qrm7fy+E2oopU/oX6PQKA3cKk+/TNUUYal\nTuMh7H+rDNQ+D0QV19jum7dQ9TyD9fWGoa7+a7qKtaTvHeBu1BV8AGq/JNt870PU71yKakzgjHDU\nsXRGT4ttMHPmuLR1PaqYKVJP82BUQ4Z0m+9ci6rvaI4q6voZFTS+RQXpO1DHRjNUvU8K6nd9F3XM\nxaD261B9Gc787qB+kytQQasV9kVRFoZq3jdqEjTcbyn27fq/1MdrVP3z2H6+C3VQ70L9IT/HenXj\naN4HUQfzMVRRwxLUn9/WB6hyXWduzNqqp3cfcA+qpVWak+kDVYfRHfWnnAvchLXooK7bZvv5dlQx\nxBlUefUHNtNyUTmip1BFA4dQLYls/8BapfeWzxdQV/Ab9HQOpipHaXM03owKZv1RLadOoopXWurT\n56ACYBaqKGZRpTQ+ggp++fr2fmUz7TdUHchrqH2wD7U/a7pitkz7BRU0XkFViK/BPpfzIeqq/qNq\nluVomYv0bclDVej/XCktzhyXFvn6tu1FtTj6EPi7Po9lvYtRldKnUY0C7tCnFaIuFibpaTmKauTQ\nXJ/+R1RO+xd93r+i9rnt734GdWw5+p1/RBXZbdOXsdTBd6o7thxNF26QgDqId6IOvof18W1Qzej2\noq5QWzucW1T2N1QFqq0RVK3g84Rp1FzZKeylUr8ydHcLRbVOqu2q2xWOjktnvYd902PhA7yZ0ygD\n/oC60hmCau7WA3gSFTSSgFX6Z1FVMuqGPEu2/h7sr1CboXILbzV80oSfuA9VtFaXCuPa1HZc1kWT\nKfIR9fNfVBPRTFQzTVBFGJleS5FvuxRVVFGEKhJ5wmZaD1RF3noapgXHVKpviy+qSqVuLag8IRtV\nXNbPzcut6bisK9umx0LYSUQVo0Rg34TOUOmzEEKIJi4cVdl3o/65cpA407DJEUIIUR1v36fRDNW6\n6EOs9xEcRxVLHUM1lztReaauXbtqBw64sxhWCCGahANUf++OU7xZEW5AtSPfhepTxuIbVBk5+muV\nm9IOHDiApmkyuGmYPXt2/eddU/95G+vgyv60G9y1HD8f3LY/ZQA3tJTzZtAYhmpzfSWqj5YtqBub\nXgTGoJrcjtI/Cx81Z+0cbyeh8Zoj+1b4Hm8WT62n+qA1uprxQgghvEjuCBekpqZ6OwmNiuxP95L9\n6Vv89eYZTS+fE15mmGNAmy2/hUcYDCDHuXAjg8EALp73vd16SgjhojZt2pCfL7czCavIyEjOnPHM\n3QoSNIRLZo+c7e0kNF6zndu3+fn5SM5b2NJzFJ5ZtseW7FlSPCWEzmAwSNAQdqo7JtxRPCUV4UII\nIZwmQUMIIYTTJGgIIRrMunXrSElJ8XYy7EyePJmvv/7a7ct9//33GTFihNuXe/z4cXr27ElpaXXP\ntvIsCRpCCI9ITExk1apVduNGjBhBZqbrTzsoLS1l+vTpJCYm0rJlSwYMGMDy5csrphuNRgICAoiI\niCAiIoKEhARuu+02fv31V7vlbNu2jW3btnHDDTdUWcc999xDQEAABw8erDYdiYmJtGjRomI948aN\nc2m70tPTCQ8Pp6ioqMq0AQMG8MYbb9ChQweuvPJKFi5c6NK66kuChnBJmjHN20lovNLSvJ0ClxgM\nBo+14ikvL6djx4789NNPnDt3jueee45bb72VnBzrgyrj4uIoLCyksLCQ9PR0UlJSGDFiBKtXr674\nzoIFC7jjjjuqLH/9+vUcPHiw1vQbDAaWLVtWsR7bwFUfQ4YMIT4+ni+++MJu/I4dO9i9ezeTJ08G\nYMqUKSxYsMClddWXBA3hEul7yoMaYd9TRqORhISEis+JiYnMmzePfv360bp1ayZNmkRJSUnF9GXL\nltG/f38iIyMZNmwY27dvB6BFixbMnj2bjh3VI8/Hjx9P586d2bx5s8P1xsXFMWfOHH73u9/xxBPW\n50ItX76ckSNH2n23vLychx9+mPnz5zvVKs3ZlmuPP/44I0aMoLCwkLNnzzJ9+nRiY2OJj4/nmWee\nwWw2AzB16lQWLVpkN++iRYsYP348kZGRAAwePJiDBw9y+HDDPzVYgoYQwmsMBgOff/45P/zwA1lZ\nWWzbto33338fgC1btjB9+nTeeustzpw5w4wZM7j++usdluUfP36cvXv30qtXrxrXN3HiRDZv3kxx\ncTFFRUVkZWWRnJxs951XXnmFkSNH0qdPH6e2YcqUKbRv356xY8eybdu2KtM1TePee+9lx44drFy5\nkoiICKZNm0bz5s05cOAAW7ZsYcWKFbz99tsA3HHHHfz000/k5uYCYDabWbJkCVOnTq1YZlBQEN26\ndSMjI8OpNLqTBA0hGjuDwT2Dhzz88MNER0cTGRnJhAkTKk6ECxcuZMaMGQwaNAiDwcBdd91FcHAw\n6enpdvOXlZUxZcoUpk2bRlJSUo3rio2NRdM0CgoKKCgoACAiIqJi+uHDh1m4cCF/+YtzT5ldvHgx\nOTk55OTkcOWVVzJ27FjOnj1rl7ZJkyZRUFDA0qVLCQkJ4fjx43z//fe88sorhIaG0q5dOx599FE+\n+eQTABISEkhNTeXDDz8EYNWqVZSUlDB+/Hi7dUdERNitq6FI0BCisdM09wweEh0dXfE+NDSU8+fP\nA5CTk8O8efOIjIysGHJzczl69GjF981mM3feeSchISG89tprta4rLy8Pg8FA69atad26NQCFhYUV\n0x999FGeffZZIiIiKoqdaip+Gjp0KMHBwYSGhvLkk0/SunVr1q1bVzF9//79LF26lGeffZagoKCK\n7SorKyMmJqZiu2bOnMnJkycr5ps6dWpF0Pjwww+ZPHkygYGBdusuLCys2IaGJEFDCOFTLJXPHTt2\n5OmnnyY/P79iOH/+PLfddhugTubTp0/n5MmTfPnll1VOqo589dVXXHLJJYSGhhIWFkbXrl3Zs2dP\nxfTVq1fz+OOPExMTQ2xsLKACgyUX4GzaLXr06MG7777LNddcw969ewGVkwgODub06dMV23X27NmK\n+hpQxWi5ubmsWbOGr776yq5oClS9y/79++nXr59T6XIn6XtKuET6nvIgJ/ue8mWlpaVcvHix4nN5\neXmt81iu7O+9914mTpzI6NGjGTRoEBcuXMBoNDJy5EjCw8O57777yMzM5McffyQ4OLjG5R05coS3\n336bd955h6VLl1ZMu/baa1m7di2XX345APv27auokNY0jZiYGJYtW0bfvn2rLPfw4cMcOnSIQYMG\nYTabmT9/PqdPn2bYsGF235s0aRKlpaWMHj0ao9FIly5duPrqq3nssceYO3cuYWFhZGVlkZeXxxVX\nXAFAWFgYN998M3fffTeJiYkMHDjQbpmbNm0iMTHRrlGBqJkmhFB89f+QmJioGQwGu2H48OFaQkKC\n3XdWrVpV8TktLU278847Kz4vX75cGzRokNa6dWstJiZGu/XWW7XCwkItOztbMxgMWmhoqBYeHl4x\nLF68WNM0TVuzZo0WEBCghYeHa2FhYVpsbKx2yy23aBs3brRL444dO7RevXpVuw0BAQHagQMHKj7P\nnDlTmzlzpqZpmrZz506tb9++WlhYmBYVFaWNHj1a++233yq++/7772sjRoyo+PzWW29pnTp10nJy\ncrSzZ89q9913nxYfH6+1atVKGzBggPbpp5/ardtoNGoGg0H7+9//XiVd999/vzZ//vxq013dMQG4\nXM4oHRYK4eekw0LXTJkyhVtvvdXhDX6+6MSJE6SmppKRkUHz5s0dfseTHRZK0BDCz0nQEJU15l5u\n3wWOA9ttxqUBucAWfXDtvnwhhBBu4+2g8R5Vg4IG/AMYoA+u3ZcvhBDCbbwdNNYBjp5T6a/FZk2O\n9D3lQX7e95RonLwdNKrzELAVeAdo+LtXhNOk7ykPaoR9Twn/54v3abwJWO7hnwvMA6ZX/lKazVVY\namoqqampDZA0IYTwH0ajEaPR6NZl+kIxUCKwFHDUO1h106T1lI8wzDGgzZbfwiMMBqe675DWU6Ky\nxtx6ypEYm/cTsW9ZJYQQwou8HTSWAP8DkoHDwD3A34BtqDqNkcAfvJY6IYRbyeNeXdfUH/c6GYgF\nmgMJqPs27gL6Av2AG1H3cQgfJX1PeZCf9z3lr497XbhwId26daNVq1YMGjSIDRs21LiNTe1xr/6q\n2j5XhGhqfPX/ULlfKXcqKirS0tLStJycHE3TNG3ZsmVaRESElp2drWma6nsqPj6+4vu5ubnas88+\nq4WEhNil6f7779deeOGFis9btmzRwsPDtc2bN2uapmlvvvmm1q5dO81sNjtMR03b+N5772nDhw+v\n87YlJydr77//vt247du3a8HBwdqZM2c0TdO0DRs2aL179652GdUdE7ih7ylv5zSEEE2Irz/uddeu\nXfTs2ZMBAwYAcOedd3Lq1ClOnDhR7TZp8rhXIYRoGL72uNcRI0aQlZXFpk2bMJlMvPvuuwwYMIAO\nHTpUu0x53KsQonFJS3P8+Nbq7jh39H0P3p3uS497TUhI4LnnnmPYsGGEhIQwd+5cFixYUO3y5HGv\nQojGJy3N8eNbawoazn7XDXzpca/ffPMN8+bNY/fu3ZSVlfHhhx9y3XXX2a3TljzuVYg6kr6nPKiJ\n9j3lzce9/vDDD4wfP55u3boBMHbsWGJiYvj555/rlHaLxvi4VwkawiXS95QHNYK+pyyPe7UMdX3c\n67///W82bdqEpmkUFRXx7bffVuRELI97/eabb2p93GteXh5z5szhnXfe4YUXXqiYZnncq0W/fv34\n9ttvycrKQtM0Vq5cyd69e+ndu3eV5R4+fJgNGzZUbONLL71U7eNeX3jhBUaPHs3BgweJiYmpeNxr\nYWEhZrOZAwcO8NNPP1XMI497db86N2MTnkGa/BYe4+Rx7qv/B3983KvJZNIef/xxLT4+XouIiNB6\n9uypffTRRxXT5XGvvtH3VH3o2y+8Tfqe8iDpe6pByONe60aChnCJBA0PkqAh6qmpdVgohBDCR0nQ\nEC6Rvqc8yM/7nhKNkxRPCeHnpHhKVCbFU0IIIXyCBA0hhBBO88VnhAsh6iAyMrLKnciiabP0husJ\n/nqkSZ2GEELUkdRpCK+Tvqc8qIn2PSV8m+Q0hEvk5j4PcvLmPiGc1RhyGu+ingG+3WZcG2AlsBdY\nATR8379CCCEc8nbQeA+o/CT2J1FBIwlYpX8WQgjhA7wdNNYB+ZXGXQ98oL//ALixQVMkhBCiWt4O\nGo50QBVZob9W/3BeIYQQDcrX79Ootv/3NJuWJampqaSmpjZMioQd6XvKg6TvKeEio9GI0Wh06zJ9\nofVUIrAU6KN/zgRSgWNADLAGSKk0j7SeEkKIOmoMracc+QawPBB3KvBfL6ZFCCGEDW/nNJYAI4G2\nqPqLZ4Gvgc+AjkA2cCtQUGk+yWkIIUQdyZP7hBBCOK2xFk8JIYTwURI0hEuk7ykPkr6nhA+S4inh\nEul7yoOk7ynhZlI8JYQQokFJ0BBCCOE0CRpCCCGc5uvdiAgP0jSNgosFZBdkk1WQpV7zs8g+m01x\nWTEtmrWoGEKDQq3vm1nfA3y1+ys6hHfg8oTLvbxFQghPk4rwRsasmblYfrFiKC4rprC0kJyCHLIL\nsu0DREEWmqbRObIzia0T6dxavSa2TiSsWRjF5cVcKLtQMRSXVfpcXszG3I10bdOV5fuXc+yPx2gd\nIo8/cZu0NGlBJdzKHRXh/kojTQ2z18zWKpu9Zna14/1tvmdXP6vtO71P+3zn59rTq57Wxn88Xoub\nF6d1eKmDljQ/Sev4Sket/UvttZZ/bakFzAmomC94brAW83KM1vnVzlrvN3pr4z8erw1aOEgbs2iM\n9sXOL7Tfjvymnb5wWjObzW5JZ/TL0Vr64XSf358yn8zXlOejmg5g68JfI46+/Y3LxfKL7Dixg4xj\nGRXDtuPbiAyNpH90f/p36E+/6H7EhMcQ2iyUkKAQuyE0KJTmgc0tVxMNasp/pjC261ju6ndXg69b\nCOEcd+Q0pE7DDTRNY/uJ7Xyd+TUFFwsICggiMCCQQENgta+W7xSVFrH1+FYyjmVwIP8ASVFJFQHi\nph430S+6H21C23h7E2uVHJVM5qlMbydDCOFhEjTqSdM0fjv6G1/u+pIvdn9BubmciSkTiY2IxWQ2\nYdJMmMwmys3llJhLMJVZx9m+BgcGc1Xnq5g1dBY92/UkOCjY25tWLyltU/h056feToYQwsMkaNSB\nWTOTnpvOl7u+5MvdX9I8sDk397yZT276hIExA71SLOQrJKchRNMgQaMWJrOJ9YfW88WuL/hP5n+I\nDInkph43sXTyUnq3792kAwWovqfSUtPoHtWdg/kHKTeXExQgh5VbSOsp4YP89Yzn0YrwgosF/JTz\nE9/t+47/Zv6XmIgYbu5xMzf1vImUtpUfIti02fY9lfhqIj/e9SPd2nTzcqoaCel7SriZVIS7SVFp\nEesPrWd11mrWZK9h96ndDIkfwtVdrmbDPRvo2qart5PoF5LbJrPn1B4JGkI0Yk0yaFwsv0h6bjpr\nstawOns1W45uYWDMQEZ1HsXLV7/MZXGX+W2FtDelRKWw5/QexjPe20kRQnhIkwka245vY9neZazO\nWs3GvI30bNeTUYmjeOaKZxiWMIyw5mHeTqLfS26bTMaxDG8nQwjhQU0maNz4yY2M6TKGRy57hCs6\nXUGrkFbeTlKjI81uhWj8fDloZAPnABNQBgx2ZWGlplKeGfkM8S3j3ZA0YTF75OyK99Ls1s1mz679\nO0I0MF8OGhqQCpxxx8JMmolAQ6A7FiVspKWmVbyPjYjlQtkF8ovziQyN9F6iGgtpbit8kK8/T8Nt\nTYJNZhOBARI0PMlgMJAclcye03u8nRQhhIf4ctDQgB+BX4F7XV2Y5DQahqXZrRCicfLl4qlhwFGg\nHbASyATWWSam2WTdU1NTSU1NrXFhktNoGJZmt0II7zMajRiNRrcu01/uCJ8NnAfm6Z/rfEd4+Avh\nHJ11lIjgCHenTdj4bOdnfLLjE/5z23+8nRQhRCXuuCPcV4unWgCWs3sYcDWw3ZUFmjTJaXhCmjHN\n7nNKW8lpuI1UhAsf5KvFUx2Ar/T3QcDHwApXFmgym6QjPQ+Ys3aOXQuq7m26c+DMgbp3XGg2Q3k5\nlJWpwfK+8rjycjCZqr6vbZwzg2UdBgMEBkJAgBocvbcdB1BSUnW4eLH68SaTmk/Tqh+2boWvvlLv\nO3eGr7923w8nRD35ak4jC+ivD72Bv1b5hsGgBkdXYw56BzVpJgLnPlfn+SrG+/N8mgYXLsCsWfDA\nA/C//8H338OSJfDmm3DVVdb5eveGceNgzBi49lq48Ubo1Qv69IG774YZM+Chh+Cxx2D4cLV8gwFG\nj4a0NEKfnk10aXOyZ05Sy0pOVssaNgwGDFCfW7Wyri8w0DqEh0O7dtChAyQkqO9ffrlK34QJkJoK\n/fvDJZfAxIlqe556CubMgRdfhAcfhD/8Ad59Fz7+GL74ApYuhXnzYPp0lfb33oP9+yE3F06dgvPn\nYdUqWLcOIiLU+qOioGVL+PlneO45tfxVq6C4GM6dg9On4dgx+PxztQ+zstS4sjIIDVUn+9deg7ff\nVt+fOBGmTlX77emnISVF7c9//Qvmz4fXX1e/wyWXwJYtkJEBl12m9u2iRfDRR7BpEzzyiG8fZzKf\nf8znIn+p06isTnUamqYR8JcAzM+a/aMr8+JidVKqPBw9ql6LiqxX37ZDaWn144OCoHVr+6FVK8ef\nW7WCZs2s81oGB58NJX9GK3vaesJs0YJx5e/zUMsxjI8aAi1aWAd9esUQHKzWExSkgoY//DYNybaX\n2zvvVEF6xgzvpkn4Nenl1klmzUyAIcAzAeP8edi9G3buhMxMdcIH6xWA5b2jcQaDOvkeP24NCMeO\nqWVER1uHmBj1OnCgugqPiFAnW2eH4GBo3tz92w4w58/qatxG8vd5ZLbqyPjLJ3tmnU3RmDHwzTcS\nNITXNYmg4ZZ7NM6dU8Fh1y4VICyvJ0+qIpeePaFHD3WCt5RJg/2ro/dBQTB4sH2QiIz066vulLYp\nbDm2xdvJaFyuugoefVTVhQRKgw7hPU0iaJSby51rOWUyqbLu/fvVsHevNUCcPq3Konv2VGX8M2eq\n9507N+k/sW3fUxbJbZNZsmOJF1LTyNj2PRUXpy5INm+GQYO8lybR5DWJoGEy2+Q0ysshJ8caGGyH\n7GxVCdqtmxq6d4crr1TBITHR2lJGVLBtOWUhXYm4SeVKzDFjYOVKCRrCq5pG0NBMBJaUqSBw+LC6\nYrMEhm7dVKucbt2gSxdVWStcIh0XesiYMfDyy6rFmBBe0jSChtlEUGmZaoo5eLCqGBYeY9tx4ZD4\nId5OTuMxciTcdptqPRcmDw0T3tE0goZmItCMqpOQgNEgUtqmsOeUBA23Cg9XLejWrVP30gifVFRa\nRO653Irh8LnD5J7LpVVwK14c/aJ/NPuvQdMIGmY9aHiq2amoQh7I5CGWeg0JGg1O0zROF5/maOFR\njp4/St65vCqBIfdcLsXlxcS3jK8YElom0LdDX15Nf5XxSeO5otMV3t4UlzSNoKGZCDRrEjQ8IM2Y\n5rAyPKUWfpV9AAAdUUlEQVRtirSgcpWjO3rHjIHf/94bqWm0zJqZY+ePcaTwSEVAsLweO3+s4vPx\nouOENQsjJiKGmPCYiqAwIGYAE5InVASINqFtHOYmAg2BzPt5nt8HDX/NJ9XpjvDs/CxG/qULOS+X\nN+nmsZ5gmGNAm131t9h2fBuTvpjErgd2eSFVjYTtHeEW5eWqq5Pdu9U9PcJpmqZxougE209sZ8eJ\nHRXDzpM7adGsBXERcRUBISY8xvpef+0Q3oGQoJB6r/9C2QUSX01k/T3rSYpKcuOWOU/uCHeSqayE\nQA0JGA2oe5vuHMw/WPeOC0XNgoJUa79Vq2DKFG+nxmcVXCxg54md1uBwUr2aNTN92vehd/veXBp7\nKdP6T6NXu14N0sqvRbMWzLhkBv9M/yevj3/d4+vzlCbxbzZdLCbQbzNV/im0WSgxETFk5WfRPaq7\nt5PTuFjqNZpg0CguK7YWGxVai49sX48UHqGwpJBe7XvRu11verfvzfXJ19O7fW+iw6O9WhH9wOAH\n6Pl6T+aOmkub0DZeS4crmkbQKL1IoCZBo6FZmt1K0HCzMWPg+edV0ZWft8Rx5HzpebYf307GsQwy\njmWw5/SeiqBQUl5CdHg0MREx6jVcvV4Wd5nd+NiIWAIMvnczbnR4NDem3Mi/f/03T43wz/ttnAka\ntwLLgXPAM8BAYC6w2YPpcitTaQlBEjQanKXZ7XVJ13k7KY1Lt26qI8rdu1VvBX7s2PljFcEh41gG\nW45t4fDZw/Rq34v+HfrTP7o/t/a6ldiIWGIiYmgV3Mrvm6z+YcgfGPvRWGYNnUVwkP/dAuBM0HgG\n+AwYDlwFvAy8CVzmwXS5lamkmECffXSIf3PU95RFclQym4/6zbWF75ldzb41GKxFVH4SNApLCtl3\nZh97Tu1h6/GtFUGi3FxO/2gVHCYkTeCZK54huW1yo64H69NB1al8uvNT7up3l7eTU2fO/DL6I8a4\nDngLWIbKafgNU2mJ1Gl4iKPmthbS7NZFNT1AZ/Ro9XCmRx5psOTUpsxUxsH8g+w9vbdi2HN6D3tP\n7+VsyVm6t+lOUlQSfTv05cHBD9I/uj9xEXF+n3Ooj8eGPsaTPz7JnX3v9LvtdyZo5AELgTHAi0AI\nvvvEP4fKSy8S6IPlm41dclu5wc9jrrpK3a9RVqaKqhqQWTOz7/Q+0nPT2Xp8a0WAOHT2EPEt40mK\nSiIpKqmiaCkpKom4lnE+WcfgLWO7jmXWilmszlrNVV2u8nZy6sTZOo1xwEtAARAD/MmTiXI3U8lF\nKZ7ygpjwGC6WX5SOCz2hbVtVt5GeDiNGeHRV+cX5bMrbRHpuOul56WzM3UirkFYMiR/CgOgBjOw0\nkqSoJLpEdvHLMnpvMBgMPDbkMf6R/o9GGTSCgTX6+zZACbDaYylSxgGvAoHA28DfXFmYqaxEgoYX\nGAwGkttKx4UeY6nXcGPQKDeXs+PEDjbmbiQ9L5303HRyz+VyaeylDIkbwsxLZvLeDe8RHS43Frpq\nSt8pPL36aXaf3E2Pdj28nRynORM0NgMdgXz9cyRwTB/uBX5zc5oCgdeA0aiisV+Ab4Dd9V2gSYqn\nvMbSB5UEDQ8YPRqefRb+8pd6L6KotIj03HTW5qxl3aF1/HrkV+JbxjMkfghD4obw6GWP0qt9r0Zd\nMe0tIUEh3HfpfbyS/goLJyz0dnKc5syRsBL4AvhB/3w1cDPwHqoV1WA3p2kwsB/I1j9/AtyAK0Gj\nrIQgVx/3Khyqru8pC0uzW1EPjvqesjV8OGzfDmfPQqtWTi3yXMk5NhzawNqctfyU8xPbjm+jX3Q/\nruh4BX+6/E8MiR8iRYkN6P5B95P0WhLPj3qedmHtvJ0cpzhz+T0Ua8AAWKGP+xnwRA+AccBhm8+5\n+rh6M5WWSk7DQ+asnVPj9OSoZDJPS2V4vcyped8SEgJDh8KaNdV+5fSF03yd+TWP/fAYly68lNh5\nsfz9f38nJCiE50Y9x4nHT7Dhng38dfRfuab7NRIwGli7sHbc0vMW3vz1TW8nxWnO5DSOAk+grvgN\nqIrx46hiJLMH0uRUT4RpNldgqamppKamVvtdU1mJ9XGvokElt02WnIYnWeo1brwRgENnD7H+0Ho2\nHNrAukPryC7IZmjCUEZ2Gsmr415lUOwgqaz2MY8OeZRRH4ziT8P+5FKHiI4YjUaMRqNbl+lM0Lgd\nmA38V/+8AZiMChq3ujU1Sh6QYPM5AZXbsJNWU7a9knIJGl4jHRd6Trm5nO2DEli/9iU2fHGaDYc3\nUGoqZVjCMIZ3HM7U/lMZGDNQ9ruP69muJwNjBvLxto+ZPnC6W5dd+YJ6Tm25Vyc4czSdBB6sZtp+\nl1NQ1a9AdyAROALchgpS9WYqKyUwQIKGN0jHhe5TWFLIxryNbDi0gfWH17MxdyPxLeMZ1qKQca0v\n5blRz9E1sqvf3SwmYNbQWTy8/GHuGXCPz/9+zgSN9qj7MnoCofo4DRjloTSVo4LUD6jczDu4UAkO\nUjzlbSltU6TjwnooD4D0Q+v5du+3rDi4gsxTmQyIHsDwjsN5ePDDXH7T5US1iIL1kyEnEkZ383aS\nRT2N6jyKoIAgVhxYwdhuY52bSdOgoADy8uDIEfV69CicPg35+Y4HN3AmaHwMfIrqRmQGMA2V+/Ck\n7/XBLUzlktPwlJr6nrKwNLuVjgtrd6LoBMv3L+e7fd+x4pkQOn3/EOO7j+fVsa8yOG6w4/qIMWNg\nxQqY7t6iDdFwbG/2G9ttLJSWqgCQl1d1sASIvDz1fJW4OIiNtb5GR0OPHhAZCW3aqFfL4GQru5o4\nEzSiUDfYPQys1YdfXV5zA5LiKc+pqbmtRUrbFH474u7beRoHs2bmtyO/8d2+7/hu/3fsObWHq7pc\nxbXdrmXe1fOIa+lEw8HRo+GJJ8BshgBpJejTzGY4cUKd+C0nf/110tHD/LnParanRNLnYBF06GAN\nBpahb1/7cRERDb4JzgSNUv31GCq3cQR1g5/fMJWXSmWgFyVHJfPx9o+9nQyfUXCxgBUHVvDdvu/4\nfv/3RIVGMb77eF686kWGdRxG88A6tmTv2FFdUW7dCgMGeCbRomaaBufOVQkEVV6PH4fWrdWJ3zZ3\nMGgQwXE38mBRd17pf5p3b/nIZ5806syZ9HmgNTALmA+0BP7gyUS5m6lTRwLzPV2iJqojzW7hYvlF\nlu1dxuLti1mVtYoRHUdwbfdrmT1yNp0jO7u+AkvTWwka7lderk72ublqyMuzvtoGBbDmACwBoVs3\nGDnSGiRiYiC4+ibPMy4Mpdv8brxQfNJnu2pxJmgs1V8LgFT9vV8FjfJePQg8dMrbyWiyLB0Xnik+\n47ePuKwPk9nEmuw1fLz9Y77O/JqBMQOZ0mcK793wHq1CXC9btjNmDLz2GvzJr/oS9T6zWZ30s7Kq\nBgXL+xMnICoK4uNVILC89u5tHyBatnQ5OVEtopjcezJv/PIGf7my/t3DOFJUWuSW5dS3zOYx4BW3\npKABmMwmaT3lRRUdF57aw9CEod5OjkdpmsavR35l8fbFfLLzE+Ii4pjSZwrPj3qe2IhYz604NRXu\nuAOKiyE0tNavNyllZZCdDQcOwP799q9ZWaqCuEsXFQzi4yExUXXRYgkOMTEN2v38o0MeZfi7w3ly\n+JO0aNaiTvNqmsbJCyfZfXI3macy2X3K+nqi6IRb0tckCvpNmgQNT6mt7ykLS7Pbxho09p3ex+Lt\ni1m8YzEms4kpfaZgnGokuW1y/RdaW99Ttlq1UpWk69erXEdTU1CgAkBWFhw8qAKCJTjk5amTf9eu\nqrioa1dVZNS1qwoWYWHeTr2dpKgkhiYM5cOtHzLj0hnVfq/MVEbmqUzr43KPZ7D12FZMmokebXuo\noV0PRncZTY+2PUhsnUjQ/7l+yq/vXSSHsb9ru6FpmuZUbyMAvL7pdXae3Mkb49/wYJKaJsMcA9rs\n2n+L5356jvOl53lx9IsNkKqGcaTwCJ/t/IzF2xdz6Owhbut1G1P6TmFQ7CD33KBlMKgKVmfNng0X\nL8LfXHqSgG+6eBFycuwDg+V9VpbKTXTubB0sAaJbN+jUCZp7ops8z1mbvZYZy2aw64FdBBgCOFdy\njq3HttoFiN0nd9OxVceKx+X2j+5Pvw79iA6Prvb408e7dHDWFHbOU30/UHXLM3mZ5DS8L6VtSqNo\nQXWm+Axf7vqSJTuWsOXYFm5IvoG5V87lqi5Xeb+F3pgx8PDD/hs0NE1VKO/eDZmZ6nX3bti7F06e\nVMVFXbpYA8Oll1rft22rgmwjcUWnKwhrHsaYD8eQU5DD0fNH6dO+D/2j+zMobhD3XnIvfdr3Iax5\nw+eSajrKwxssFR5mMpvkPg0vs9zg54+KSov4Zs83LNmxhLU5axnTZQwPDHqAa7tfS2gzH6o/uOwy\nVSRz8iS08+Futk0mlVOwBAXLkJmpWhb16GEdbrgBkpJUwAhqEqXpgMoRvHv9u+w6uYsBMQPo3qa7\nz5zDmsSvYNJM3r8KbOK6telGVn6W33RcWGoq5Yf9P7BkxxK+2/cdQ+KHcHuf2/no/31Ey2DXW8l4\nRLNmqqx+1SqYNMnbqVFOnlT3j2RkqNetW2HfPutdyykpMGwY/O536nNUlLdT7DP6RfejX3Q/byej\nCt//97pBublciqe8LLRZKLERsT7dcaGmaWw4vIFFWxfx5e4v6dG2B7f3uZ1Xx71K+7D23k6ec0aP\nhh9/bPigYTKpYGAbIDIy4MIF6NcP+veHK6+ERx9VwaGFX5VwCxtNImhI8ZTnONP3lEVyW1VE5WtB\nw6yZ+Trza17c8CL5xfn8buDv2Pz7zXRq3cm7CZvt/L6tMGYMzJun6gc8VcZ/7hxs22bNOWzdCjt2\nqNyDJUDMmKFeO3ZsVHUNwsVadC+qU+up3HO5mDUzHVt19GCSRG3+sPwPxLWM44+X/9HbSQGgpLyE\nj7Z9xEv/e4mWwS15YtgT3Jhyo39fYGgaJCTA6tWqLsDVZWVl2QeHrVvV3dG9eqkAYQkSffu65eY2\n4Vmebj3VaMS3jPd2EgQqp/HrEe/3dXmu5BwLfl3AqxtfpU/7Prw5/k1SE1N9/jkGTjEYrF2K1CVo\nlJWp4qQtW6zBYft2FQgswWHyZHjxRdWM1Uf7RRKe1ySChvAN3m52e+z8Mf6Z/k/e2vwWY7qOYdnk\nZQyIaYR9NY0ZA599Bg88UP13ysrg11/BaFTDzz+rpqsDB6oAcfPNKvcgFdOiEn+9tKpT8ZTwDUcL\nj9L33305+XjDdh65/8x+Xv7fy3y28zMm957MrMtn0SWyS4OmoUEdPw7JyXDqlLWZamkp/PILrF2r\ngkR6urUzvdRUGDFC9ZQrGjUpnhJ+JTo8mpLykgbruHDXyV2kGdNYk72GmZfMJPPBTP9pBeWKDh3U\nXdBvv60Cx9q1KkgkJakA8dBD8Omnqs8lIepIntgiXJJmTHP6uwaDQfVB5eFu0k9dOMWD3z3IyPdH\nMjhuMAcfPsjcUXP9L2A42++UI7ffroJGQQE88ggcPgy//aZaVk2YIAFD1JsUTwmXONv3lMWdX93J\nqMRR3D3gbrenpdRUyhu/vMHz655nUq9JpKWmqWdo+6u69j0lRC0aa/FUGvA7rM8h/zOw3GupEW6V\nEqV6u3UnTdP4dt+3zFqh6irWTltLz3Y93boOIYTii0FDA/6hD6KRSW6bzEfbPnLb8nac2MFjPzzG\n4XOHeXXsq1zT/Rq3LVsIUZWv1mn4a7GZqEVyVLJbchoni05y/7f3M+qDUVyffD3bZm6TgCFEA/DV\noPEQsBV4B/V8ctFIdI/qTlZ+FmWmsnrNX2oqZd7/5tHj9R40C2hG5oOZPDj4QZoFNtyT1YRoyrxV\nPLUScPTU9KeBNwHLw3HnAvOA6ZW/mGbTsiQ1NZXU1FR3p1E4oS59TwGEBIWojgsLskiKqls3F0v3\nLOWxFY+RFJXEurvX0aNdjzrN73fq0/eUEDaMRiNGo9Gty/T1YqBEYCnQp9J4aT3lx679+Fruu/Q+\nJiRPcOr7+cX5PPj9g2zK28T8a+Yzrts4D6dQiMbJHa2nfLF4Ksbm/URgu7cSIjyjLg9k+n7f9/R5\nsw9RoVFkzMiQgCGEl/li66m/Af1RraiygOqfrC78UkrbFH458kuN3yksKWTWilmsOLCCRRMXMarz\nqAZKnRCiJr6Y07gL6Av0A24Ejns3OcLdLM/VqM7a7LX0+3c/TGYTW2dulYAhhA/xxZyGaOSqa3Zb\nXFbM06uf5tOdn7LgugVcl3SdF1InhKiJL+Y0hB+pS99TFtHh0ZSaSjl94XTFuE15mxiwYABHCo+w\nbeY2CRjgWt9TQniIBA3hkjlr59R5HoPBUJHbKDWV8szqZ5iwZAJzUufwyc2f+Hd/Ue40p+77VghP\nk+Ip4RUpbVP4z+7/8MB3DxAXEUfGjAxiImJqn1EI4VWS0xBe0aNtDxb+tpCHBj/E0slLJWAI4Sd8\n/ea+6sjNfT6irl2jW5wvPc+Fsgv+94yLhiRdows3a6xdo4smILx5OOHNw72dDCFEHUnxlHBJXfue\nEnUgfU8JHyTFU0II0UQ01r6nhBBC+CgJGkIIIZwmQUMIIYTTJGgIIYRwmgQN4ZL69D0lnCR9Twkf\nJK2nhEvqe3OfcILc3CfcTFpPCSGEaFASNIQQQjhNgoYQQginSdAQQgjhNG8FjVuAnYAJGFhp2p+B\nfUAmcHUDp0vUkfQ95UHS95TwQd5qPZUCmIEFwCxgsz6+J7AYGATEAT8CSfp3bUnrKSGEqCN/bj2V\nCex1MP4GYAlQBmQD+4HBDZcsIYQQNfG1Oo1YINfmcy4qxyGEEMIHePIhTCuBaAfjnwKW1mE5Dsuh\n0mzulk1NTSU1NbUOixRCiMbPaDRiNBrdukxv3xG+Bvs6jSf11xf11+XAbGBjpfmkTkMIIerIn+s0\nbNluwDfAJKA50BnoDmzyRqKEc6TvKQ+SvqeED/JWTmMi8C+gLXAW2AJco097CrgHKAceAX5wML/k\nNHyE9D3lQdL3lHAzd+Q0vF08VV8SNHyEBA0PkqAh3KyxFE8JIYTwExI0hBBCOE2ChhBCCKdJ0BAu\nkb6nPEj6nhI+SCrChRCiiZCKcCGEEA1KgoYQQginSdAQQgjhNAkaQgghnCZBQ7hE+p7yIOl7Svgg\naT0lXCLdiHiQdCMi3ExaTwkhhGhQEjSEEEI4TYKGEEIIp0nQEEII4TQJGsIl0veUB0nfU8IHSesp\nIYRoIqT1lBBCiAblraBxC7ATMAEDbcYnAsWoZ4ZvAd5o8JQJIYSoVpCX1rsdmAgscDBtPzCgYZMj\nhBDCGd4KGpleWq8QQggX+GKdRmdU0ZQRGO7dpIjaSN9THiR9Twkf5MnWUyuBaAfjnwKW6u/XALOA\nzfrn5kAYkI+q6/gv0AsorLQMaT3lI6TvKQ+SvqeEm7mj9ZQni6fG1GOeUn0AFUgOAN2xBpUKaTZX\nYampqaSmptZjdUII0XgZjUaMRqNbl+nt+zTWAH8EftM/t0XlMkxAF+AnoDdQUGk+yWn4CMlpeJDk\nNISb+fN9GhOBw8AQ4Fvge338SGArqk7jc2AGVQOGEEIIL/F2TqO+JKfhIySn4UGS0xBu5s85DdFI\nSN9THiR9TwkfJDkNIYRoIiSnIYQQokFJ0BBCCOE0CRpCCCGcJkFDCCGE0yRoCJdI31MeJH1PCR8k\nraeES+Q+DQ+S+zSEm0nrKSGEEA1KgoYQQginSdAQQgjhNAkaQgghnCZBQ7hE+p7yIOl7SvggaT0l\nhBBNhLSeEkII0aAkaAghhHCaBA0hhBBOk6AhhBDCaRI0hEuk7ykPkr6nhA/yVuupl4DrgFLgAHA3\ncFaf9mfgHsAEPAyscDC/tJ7yEdL3lAdJ31PCzfy59dQKoBfQD9iLChQAPYHb9NdxwBtIbsjjjEaj\nt5PQqMj+dC/Zn77FWyfklYBZf78RiNff3wAsAcqAbGA/MLihE9fUyJ/SvWR/upfsT9/iC1fx9wDf\n6e9jgVybablAXIOnSAghhENBHlz2SiDawfingKX6+6dR9RqLa1iOFOoKIYSP8GY3ItOAe4GrgIv6\nuCf11xf11+XAbFQRlq39QFcPp08IIRqbA0A3byeiPsYBO4G2lcb3BDKA5kBn1Ab6a/9YQggh3GQf\nkANs0Yc3bKY9hcpJZAJjGz5pQgghhBBCiEZpHCqXsQ94oprv/EufvhUYUMd5mxpX9mc2sA2VG9zk\nuST6jdr2ZQrwM6qOblYd522KXNmf2cixWVlt+3MK6j++DdgA9K3DvD4rEFU0lQg0Q9Vv9Kj0nWux\nNtG9DEivw7xNjSv7EyALaOPZJPoNZ/ZlO+BS4DnsT3JybFblyv4EOTYrc2Z/DgVa6e/HUc9zpy/c\np2FrMCrx2agb/D5B3fBn63rgA/39RqA1qmmvM/M2NfXdnx1spktDBMWZfXkS+FWfXtd5mxpX9qeF\nHJtWzuzPn7F212R7U3Wdjk9fCxpxwGGbz45u7qvuO7FOzNvUuLI/Qd0j8yPqj3uvh9LoL5zZl56Y\nt7FydZ/IsWmvrvtzOtYShjrN68mb++rD2Rv55ArDOa7uz+HAEVQxwUpUmec6N6TLH7lyk6ncoFqV\nq/tkGHAUOTYt6rI/r0T1xDGsHvP6XE4jD0iw+ZyAfbcijr4Tr3/HmXmbmvruzzz9/RH99STwFU27\nHzBXji85NqtydZ8c1V/l2FSc3Z99gbdQxdL5dZzXJwWhbuhLRN3gV1vF7RCslTnOzNvUuLI/WwAR\n+vswVGuLqz2YVl9Xl+MrDfuKWzk2q3Jlf8qxWZUz+7Mjqu5iSD3m9WnXAHtQG2fpMn2GPli8pk/f\nCgysZd6mrr77swvq4MkAdiD7E2rfl9GosuGzqKu4Q0B4DfM2dfXdn3JsOlbb/nwbOI31pupNtcwr\nhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQomkwYW1rvgV105KvuwT4p5uW9SPWm9vOV5o2DZhf\nw7zXA8+4KR1CCOEXCmuYZqBx9002Cnjd5nPlfTGVmoOGAXVTXDM3p0s0Eb7W95QQ9ZGIupv1A2A7\nqu+cx1F3vG5FdUNh8bT+3XXAYqzdUxhRuQFQz67P0t8HAi/ZLOv3+vhUfZ7Pgd3ARzbrGITq2iID\n1QV1uP79pfr0MOBdfdpm1NU/QC993BZ9Xd0cbOvtwNeOdwNgHzAzsObGLgAjUJ3T/Yx0uyGEaELK\nsZ4MvwQ6oYqsLJ3WXQ0s0N8HoE7WI1BBYRsQgire2Qc8pn9vDdYuVGyDxu9RgQYgGPgFFaRSgQJU\nl/wG4H/A5ai+ew5gDUDhqMCTijVovIB6ihqo55fsQfWn9C9UUADVH1CIg23fjf3Dh2z3xRYgR1+O\nrQnAWj0dAHcDf3OwbCFq5WtdowvhjGLsH0ubiDpZWvrSuVoftuifw4DuqEDxH9TjQy8C3zixrquB\nPsDN+ueWqBxAmb4+S0/AGUBnVHHRUeA3fXzlOgfLMicAf9Q/B6PqZX5GBah4PZ37HcwbC5yx+Vx5\nX0xFPe3Oojvwd1TQMunjjqCe3CZEnUnQEI1FUaXPfwUWVhr3CPbFN7bvy7EW11a+wn8Q9cwGW6lA\nic1nE+r/5OyzCf4fKqdjKxPVy/B1qJ6HZ6ByQHVhu03hwKfA74DjNuMD6pBOIexInYZojH5APWQm\nTP8ch3pYz0/AjViLp66zmScb6xX6zTbjfwDux3qBlYQqSnJEQxU1xdgsKwJrsZDtMh+2+WzJKXRG\nFYvNR9Vb9HGwjiNAVDXrr+xd4D1U/YqtGFTOTIg6k6Ah/JGjq2TbcStRldw/o+owPkNddW9BXXlv\nRV3J/4L1yvxl4D5UxXSUzfLeBnbp47cDb2LNUThKRxlwG+rEn4EKECGVvj8X1XppG6pr7zn6+Fv1\nz1tQleKLHCx/PfbFT5XTYFlPR+AmVPC01HdY6mwGowKoEEKIOpiN/cN9/EEqKnDVVwAqmEnRtKgX\nyWmIps7fyvaNWCv16+M64AtUHY4QQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEII9/r/6JSAjRYx\n9pkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11058d9d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plots = []\n",
    "colors = ['r','g']\n",
    "\n",
    "# Plot lag-frequency spectrum\n",
    "for i in range(0,len(lags)):\n",
    "    plots += plt.plot(cross_spectrums[i].freq, lags[i], colors[i], label=str(energies[i])+'keV')\n",
    "    plt.axvline(v_cutoffs[i],color=colors[i],linestyle='--')\n",
    "    plt.axhline(h_cutoffs[i], color=colors[i], linestyle='-.')\n",
    "\n",
    "# Define axes and add labels\n",
    "plt.axis([0,0.2,-20,20])\n",
    "plt.legend(plots)\n",
    "plt.xlabel('Frequencies (Hz)')\n",
    "plt.ylabel('Lags')\n",
    "plt.title('Energy Dependent Frequency-lag Spectrum')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note:\n",
    "\n",
    "Currently, lag-energy spectrum isn't plotted and hence I am unable to verify results from Uttley et al. However, as soon as it is implemented in library project, I will test it here as well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### With same position and varying intensity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we use delta impulse responses whose position remains same but intensity varies. \n",
    "\n",
    "Again, first we define energies and then create impulse responses, and subsequently using convolution, obtain the output light curves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 615,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "energies = np.array([4.5,8.5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 616,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "h_zeros = np.zeros(int(10/lc.dt))\n",
    "responses = [np.append(h_zeros, i+1) for i in range(0,len(energies))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 617,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "delay = int(10/lc.dt)\n",
    "outputs = [signal.fftconvolve(s, h)[delay:-delay] for h in responses]\n",
    "s_mod = s[delay:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 618,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "t_mod = lc.time[delay:] \n",
    "lc_input = Lightcurve(t_mod, s_mod)\n",
    "lc_output = [Lightcurve(t_mod, output) for output in outputs]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 619,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cross_spectrums = [Crossspectrum(lc_input, lc2).rebin(0.0075) for lc2 in lc_output]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 620,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "lags = [np.angle(cross.cs)/ (2 * np.pi * cross.freq) for cross in cross_spectrums]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 621,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "v_cutoff = 1.0/(2.0*10) \n",
    "h_cutoff = lags[0][int((v_cutoff-0.0075)*1/0.0075)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 622,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEZCAYAAABrUHmEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8FHX+x/HXpoceUCCQQOhlKYKCcAgERUFB0d8polhQ\nLOid7RTPbrBwNs7z8O4Uyykq9lMBKSISQKS3kFASei9SQ0nd/f3xnU02IQmbZJPdZN/PxyOP7Mzs\nzHxmdnY+8/1+Z78DIiIiIiIiIiIiIiIiIiIiIiIiIiIiIlKNJACf+DqIAJOA9rlYgnwdQDWzHTgN\npLv9/dOXAXnIAZzExPs78DMw3KcRFc/ppeXEYba7pO9AApBNwc/zMS+tvyrx1j4vymhgA3AC2A/8\nCNSqwPXFce7PXUoQ4usAqhknMBT4pYLXEwzkenmZXYCtQH3gKuBtoD3wgpfX429sJUxzAp8Dt51j\nGUGYE1F1VdI+Ko/+wMvAIGAtEIX5/lSGkrapIr5f1YaybeUZBfwKvA4cwZygB7tNrwt8AOwFdgMv\nkv/5jAIWAX/HlASex5zcpwHHgWXAS8BC6/3/At4otP6pwMMexHkE+BS4D3jSWo+n8U0EjmGuHC8t\nxbaVtF9aAPMxV6I/AecVircX8BtwFFiDORG5JGKS3q/W/LOBBta0Bdb/Y5gSxMVF7AsbRZ9cPgL+\nA8zAlNDigSbAt8BBaxsecHt/pDXPESAFGAvscpvuAFoWWv6LbsNDrW07itnPnd2mbQcexZx0jwFf\nAOFu04dZ8x4HNmNO0DcAKwpt01+A74vY1qJ8Deyz1jcf6Og2rQHFH5eF9QAWW7GD2b5PMPsUzH54\nB/O5n8B8ns3c5m8PzAEOAxut7XKJBCZg9s8xzOcdQcHP/QTm+BlFwe9XAmdXycVRsISSiPmMFmGO\nn6mYY/Mzt21vXsx2i+TZBlxWzLRRQBamOG4DxgB73KZ/hzkRRQLnA0uBe9zmzQb+hDloIzAnhynW\n6w7ATvK/ED2sZbtOeOcBp6zlFqXwSQsg1FrnoFLE9xDmKm045ktZz8N5S9ovizEJMBToi/miT7am\nNcV8yV1JZqA17EoMiUAa0Bqzn+YBf7OmNcez6qmi6vI/sravtzUcCawEnsGU3lsAW4ArrOmvYE6u\n9YAYIBnzebkU3v//Jb+E1w04gPlMbZhSzzbM/sB6vQRojLlSXw/ca03racXpOiabAO2AMMyJtr3b\nOlcD1xW1Ezh7P4wCaloxvGnN61LScVnYJZjq3ASgDwWTHZj9fMJ6XxjwD/ITUE1M4r0d8xleAByy\n1gnmwukXINqa3staRlGf+yjO/n49z7mTRirms66DuRhIw1wsBQMfAx8Ws90iebZjrjqOuv2NtqaN\nwhxULjUwB2FDoBGQgTlYXW4iv5prFLDDbVow5kTbxm3cixS8oluPOYkC/BmYXkLcRSUNMFeTN3kY\n3x4KWgrc4uG8xe2XZpgvc6Tb9M/ITxp/dXvtMov86qR5wFNu0+4DZlqv4/AsaWSS/1kewZyEPrL+\nXC6m4OcDppTmOmm4JxCAuym5pOGeNP7D2VWEGzEJFEzSuNlt2qvWPADvYq62i/IfTCkAwI7ZttBi\n3ptA8Q3h9az4a+PZcVnYYMxV+lHMd2cC+Z/JR5gE5FITyMEk3hs5Oxm9CzxnzX+agiUylziKThqF\nP78ESk4a8zCfscsbmPYYl6EUTKbVhto0vMuJqQ4ork1jv9vr09b/WpiSQCjmJO0SRMGrUfeTzPmY\nz8593O5C65qMOWn/bP1/89zhFxBqrecI5ursXPEVTho7MFe2zTyYt7j90hBzMjlTaLmx1uvmmCqJ\nq92mh1Bw/7sv+wylb2T9krPbNJwU3N7mmG096jYumPyTWhMKflbu234uza31u1d3hVrLdCm8jdHW\n6xgKnsjcfYw5IT8D3IrZzmxgJKZKCCv+IYXmC8a0Q1yPOT4cmP1xnrXucx2Xhc2y/sBcpX8NbAIm\nWct1n/8U5nhsgtkvF1Nwn4dgjvsGmIuULedYt7td537LWQ64vc7AVE26D1dkg77PKGn4h12YK9oG\nFN+g6n4HyyHMFVcs+VfpsYXe/ymwDuiKqYbwtL7aZZi1jmWYL+C54mtaaLg58AOebVtx9mGqXGqQ\nn0yak99IuRNzNXjP2bOekyd3BDkpvsHUff6dmCv+tsW8dx8meW6whpsVmn4as40u0eSfxHZiTtLj\nPYi3sF2YqrmiLMGUCvphSn43WeM/s/6KczNwDabKawempHEEs588OS5L8ov1Z7eGbYXmr4VpY9uD\n2S/zKViCcwnCnLRbA0mFphX3uRcef5KCn0njc8RekXeY+RU1hHtfWe402Ydp7Ps7ppgfBLTCfKGL\nkgv8D1OEjsQkhVspeODuxjR2Tga+wZy4PYm7PuZq821MXfxRD+NrCDyIuQq+wYppBuYquDTb5m6H\ntQ3jrOVeQsG7az7FlDKuwFwBR2Aapd0TWHGfxyFMEmtVwvqLm7fw+GWYqpXHMZ9HMNAJuMia/hWm\nKsPVpvEABT+rNZh9HoyprnHfN+9h2nl6Wuutibn6L+kq1hXfB8AdmCv4IMx+aef2vk8wn3MW5mYC\nT9TCHEtHrFjck5knx6W7azDVTFFWzD0xNzIscXvPVZj2jjBMVddiTNL4EZOkb8EcG6GYdp/2mM/1\nQ8wxF43Zr72tZXjyuYP5TPphklZdClZFudiKeV2tKWl43zQK3tf/rTXeydlfHvfh2zAH9XrMF/Jr\n8q9uipr3z5iDeT+mquFzzJff3ceYel1Pfpi11oo3DbgTc6dVgofxgWnDaIP5Ur4I/JH8qoPSbpv7\n8M2YaogjmPrqj92m7caUiJ7CVA3sxNxJ5P4FdhZ67Ro+jbmCX2TF2ZOzFRVbUeMdmGR2AebOqUOY\n6pU61vRxmAS4DVMVM7lQjA9hkt9Ra3u/c5u2EtMG8jZmH6Rh9mdJV8yuacsxSeNNTIP4PAqWcj7B\nXNV/WsyyilrmZGtb9mAa9BcXisWT49LlqLVtqZg7jj4BXrPmca13CqZR+jDmpoBbrGnpmIuFEVYs\n+zA3OYRZ0x/DlLSXW/P+DbPP3T/3I5hjq6jP+WdMlV2StYxpRbynuGOrqOkifudVTAOqu76c3cBX\nEUZRcmOnFBRP2erQvS0Sc3fSua66y6Oo49JT/6XgrcfiB1TSqLraYX6Q5yrW30nBK9RQTGnhvcoP\nTaqI+zBVa6VpMD6Xcx2XpREwVT5ViRrCq67amGJ8E8xdHG9gbl0Ec6/6cky97D8qIZbiqnGkeL7e\nX9utGK718nJLOi5LS8eViIiIiIhIQKiSdYZdu3Z1rl279txvFBERd2sxd/mVWZVMGoDT6VRVp7fE\nj4on8aNEX4dRbSQkJJCQkODrMKoN7U/vsdlsUM7zvu6eEuZvn+/rEESkilDSEBERjylpiOm/U7wm\nPj7e1yFUK9qf/kVtGoJtnA3n89qfItWdN9o09OM+kSqufv36HD169NxvlIARFRXFkSNHKmTZShrC\n8/2f93UIUg5Hjx5FJW9xZ5UoKmbZFbbkiqXqKRGLzWZT0pACijsmdMutiIhUKiUNERHxmJKGiFSa\nhQsX0r59e1+HUcBNN93EDz/84PXlfvTRR/Tt29fryz1w4AAdO3YkK6u4Z1tVLCUNEakQcXFxzJ07\nt8C4vn37snHjxnIvOysri9GjRxMXF0edOnXo1q0bs2bNypuemJhIUFAQtWvXpnbt2sTGxnLjjTey\nYsWKAstJSkoiKSmJYcOGnbWOO++8k6CgILZu3VpsHHFxcdSoUSNvPYMHDy7Xdi1ZsoRatWpx6tSp\ns6Z169aNf//73zRq1IgBAwYwadKkcq2rrJQ0hITEBF+HINWQzWarsLt4cnJyaNasGQsWLODEiRO8\n9NJLDB8+nB078h9U2bRpU9LT00lPT2fJkiW0b9+evn378ssvv+S959133+WWW245a/m//vorW7du\nPWf8NpuN6dOn563HPXGVRa9evYiJieGbb74pMD45OZkNGzZw0003ATBy5Ejefffdcq2rrJQ0hHHz\nx/k6BAkQiYmJxMbG5g3HxcUxYcIEunbtSr169RgxYgSZmZl506dPn84FF1xAVFQUffr0Yd26dQDU\nqFGD559/nmbNzCPPhwwZQosWLVi1alWR623atCnjxo3jrrvu4q9//Wve+FmzZtG/f/8C783JyeHB\nBx9k4sSJHt2V5umda2PHjqVv376kp6dz/PhxRo8eTZMmTYiJieHZZ5/F4XAAcPvttzN58uQC806e\nPJkhQ4YQFRUFQM+ePdm6dSu7dlX+U4OVNETEZ2w2G19//TWzZ89m27ZtJCUl8dFHHwGwevVqRo8e\nzXvvvceRI0e49957ueaaa4qsyz9w4ACpqanY7fYS13fdddexatUqzpw5w6lTp9i2bRvt2rUr8J43\n33yT/v3707lzZ4+2YeTIkTRs2JBBgwaRlJR01nSn08ndd99NcnIyc+bMoXbt2owaNYqwsDC2bNnC\n6tWr+emnn3j//fcBuOWWW1iwYAG7d+8GwOFw8Pnnn3P77bfnLTMkJITWrVuzZs0aj2L0JiUNkerO\nZvPOXwV58MEHady4MVFRUVx99dV5J8JJkyZx77330qNHD2w2G7fddhvh4eEsWbKkwPzZ2dmMHDmS\nUaNG0bZt2xLX1aRJE5xOJ8eOHePYsWMA1K5dO2/6rl27mDRpEi+88IJHsU+ZMoUdO3awY8cOBgwY\nwKBBgzh+/HiB2EaMGMGxY8eYNm0aERERHDhwgJkzZ/Lmm28SGRnJ+eefz8MPP8wXX3wBQGxsLPHx\n8XzyyScAzJ07l8zMTIYMGVJg3bVr1y6wrsqipCFS3Tmd3vmrII0bN857HRkZycmTJwHYsWMHEyZM\nICoqKu9v9+7d7Nu3L+/9DoeDW2+9lYiICN5+++1zrmvPnj3YbDbq1atHvXr1AEhPT8+b/vDDD/Pc\nc89Ru3btvGqnkqqfevfuTXh4OJGRkTzxxBPUq1ePhQsX5k3fvHkz06ZN47nnniMkJCRvu7Kzs4mO\njs7brjFjxnDo0KG8+W6//fa8pPHJJ59w0003ERwcXGDd6enpedtQmXyZNGKBeUAKkAw8aI2vD8wB\nUoGfgMrfKyLiM67G52bNmvH0009z9OjRvL+TJ09y4403AuZkPnr0aA4dOsS333571km1KN999x0X\nXnghkZGR1KxZk1atWrFp06a86b/88gtjx44lOjqaJk2aACYxuEoBnsbu0qFDBz788EOuvPJKUlNT\nAVOSCA8P5/Dhw3nbdfz48bz2GjDVaLt372bevHl89913BaqmwLS7bN68ma5du3oUlzf5MmlkA48A\ndqAX8CegA/AEJmm0BeZaw1KB1PeUVJSsrCwyMjLy/nJycs45j+vK/u677+add95h2bJlOJ1OTp06\nxY8//phXErnvvvvYuHEjU6dOJTw8vMTl7dmzh3HjxvHBBx8wfvz4vGlXXXUV8+fnP4QsLS2NpKQk\n1q5dm1dNNn36dK699tqzlrtr1y4WLVqUt42vv/46hw8fpk+fPgXeN2LECMaPH8/AgQPZunUr0dHR\nXHHFFfzlL38hPT0dh8PBli1bWLBgQd48NWvW5Prrr+eOO+4gLi6O7t27F1jmsmXLiIuLK3BTQSD6\nHhgIbAQaWeMaW8OFOUXE8NfvQ1xcnNNmsxX4u+SSS5yxsbEF3jN37ty84YSEBOett96aNzxr1ixn\njx49nPXq1XNGR0c7hw8f7kxPT3du377dabPZnJGRkc5atWrl/U2ZMsXpdDqd8+bNcwYFBTlr1arl\nrFmzprNJkybOG264wbl06dICMSYnJzvtdnux2xAUFOTcsmVL3vCYMWOcY8aMcTqdTmdKSoqzS5cu\nzpo1azobNGjgHDhwoHPlypV57/3oo4+cffv2zRt+7733nM2bN3fu2LHDefz4ced9993njImJcdat\nW9fZrVs355dffllg3YmJiU6bzeZ87bXXzorr/vvvd06cOLHYuIs7JoBy1zP6S4eFccB8oBOwE4iy\nxtuAI27DLtb2i4g6LCyfkSNHMnz48CJ/4OePDh48SHx8PGvWrCEsLKzI91Rkh4X+0DV6LeBb4CEg\nvdC0YjOj+4Pm4+Pj9XQvESmTzz77zNchlErDhg1Zv369R+9NTEwkMTHRq+v3dUkjFJgOzAT+YY3b\nCMQD+4FoTGN54c5qVNIQsaikIYVV167RbcAHwHryEwbAVMB1q8DtmLYOERHxA75MGn2AW4ABwGrr\nbzDwCnA55pbbS61hqUDqe0pEPOXr6qmyUvWUF9nG2XA+r/1ZVal6SgqrrtVTIiJSxShpiIiIx5Q0\nRETEY0oaIlJp9LjX8tPjXsXn1PeUVISq+rjXSZMm0bp1a+rWrUuPHj1YtGhRiduox71KwEmIT/B1\nCFINVcXHva5Zs4ZHH32Ur7/+Ou/petddd12xd6fpca8iIhXI3x/3un79ejp27Ei3bt0AuPXWW/n9\n9985ePBgsdvk6e3OetyriEg5+dvjXvv27cu2bdtYtmwZubm5fPjhh3Tr1o1GjRoVu8xAe9yrP3RY\nKCIVyDbOO1VEFfUDUNfjXoFiH/cKcNtttzF+/HiWLFlCv3798uYv6+NeXVf27o97jY2N5aWXXsp7\nJkZUVBQzZswodnlTpkyhe/fuOBwO3nrrLQYNGsTGjRupW7duXmwjRozA4XAwbdo0QkJC8h73euzY\nMSIiIoiMjOThhx/mvffe45577inwuNcnn3zS7x73qqQhUs35+6/9Cz/ude/evYB5LOrkyZOZOHFi\n3vTs7GyvPe7VlTTS09Np0KABAFOnTmXChAls2LCB1q1bM3v2bIYOHcrq1auJjo4+a3m9e/fOe/3E\nE0/w8ccfs3DhQoYOHQqYx70mJSWxdOnSIh/36r4drqo2MFVU48eP58knn9TjXsX/qO8p8Se+fNzr\n7NmzGTJkCK1btwZg0KBBREdHs3jx4lLF7qLHvUq1NG7+OF+HINVUVXvca9euXfnxxx/Ztm0bTqeT\nOXPmkJqaSqdOnc5arh73WrUU+5hDKT0StD+rMn/9PlTFx73m5uY6x44d64yJiXHWrl3b2bFjR+en\nn36aN12Pe1Uvt4J6ua3q1Mtt+ehxr6WjpCFKGlWckoYUpq7RRUTELyhpiPqeEhGPqXpKpIpT9ZQU\npuopERHxC0oaIiLiMXUjIlLFRUVFVVgX5FI1uXrDrQhV9UhTm4aISCmpTUO8Qn1PiYinVNIQ/bhP\nJECopCEiIpVKSUNERDympCEiIh5T0hAREY8paYj6nhIRj+nuKRGRAKG7p0REpFIpaYiIiMeUNERE\nxGNKGiIi4jElDVHfUyLiMV8njQ+BA8A6t3H1gTlAKvATUM8HcQWUcfPH+ToEEakifJ00/gsMLjTu\nCUzSaAvMtYZFRMQP+DppLASOFhp3DfCx9fpj4NpKjUhERIrl66RRlEaYKius/418GIuIiLjxx6Th\nzmn9iYiIH/DHZ4QfABoD+4Fo4GBRb0pISMh7HR8fT3x8fCWEVj2p7ymR6ikxMZHExESvLtMf+p6K\nA6YBna3h14DDwKuYRvB6nN0Yrr6nKkjmqRPs3LiUpq27UaPueb4OR0S8yBt9T/k6aXwO9AfOw5Qw\nngN+AL4CmgHbgeHAsULzKWmU0ZkTR9i5aSnbt65i+9717Diyje2n97I99zDbw05zONxBRC6MjbiM\nZ5772dfhiogXVYekUVYBmTROHNpNyooZ/P77LjKzTpOZfYaMrNNkZmeQmZ1BRs4ZMnMyyczNJCMn\nk0xHFpmOLNJzz7DTeZTtYWc4FuYg9nQocbm1aR56HnG1Y4k7vw1xsZ1p3upCmrTuxpcfPML3qVP5\n6u+7fL3JIuJFAd3LrW2cDds4W5G/Zk5ITCh2fFWY77m5z+TN1/f5WK55pDEtHgsl+h+xPPjTIzy6\n/u88u3kSP6bNYNGu31h3cB07ju/gx6wUxtnm80rIEtbXPIX9/I70junNde2uZcLA17m552j+OuBZ\n0iZkMecfh3n/9U0889zPbO7QiL4r/0Szr3rx4qLx2Dv2J4VDfrdfNJ/m03zln6+8VNLwEqfDQerK\n2Rw/up/g4FCCQ0IJCQmz/ocTHGKNCw0n2DU+NILT6YdJXj2bpC2/se7IRtbl7mNTzTM0PRNKF+f5\ndK7Ths7NetCl6xW07NKf4NCwCt+WjJPHiHoliuNPpxMWWavC1ycilcMbJQ1/vHuqysjNzmLx7Pf5\nfvGHfJ+xlswgJ41zwsnFSY7NSS5Ocm1OcnCQa4NcmzVsc1rDEO6wYc+sS+fI5sQ368cDHeLpeNGV\n1IxqWGnbkZCYQEJ8Qt5wRK16NDsTSuqqOXTqc12lxSEi/k8ljVLKOHmMn7//O9+v/ZJpQWk0zgrn\n2ro9uTZ+DBf0uxFbUNWr8bONs+F8vuD+/L9HmnBjxxu48e63fBSViHibShqV5Oi+bcz44XW+T53K\nTxF7uOB0Xa5tHM9TV75Ly67xvg6vQnSq05qU3at9HYaI+BkljSI4HQ7SVs1hTuIHfL97LktrHGHA\nmcZc2/JK/j1sLOc36+DrECucPaYbX234xtdhiIifUdLAtE0kLfqWhcu+YeH+pSwM3UeYw8YAZ3Pu\nv+Auvh/2aKW2MfgDe8d4Uja+4+swRMTPBGTSyDx1guXzPmXhmh9YeHg1v0UcIjozjH4hLRnW5mre\n6HcLze19fB2mT7XtfjnbZ2WRcfIYEbX0SBMRMQImaSyZ9T4/Lv6YBcfXsbLmcdqfrkG/yPbc3fUO\nPr50VEBUORWnqL6nwiJr0fJUGJtW/UTXfsN9EJWI+KOAuXuqxWOhDAvtxFXdhtP7stup3aBJBYVW\nfQz/SyzXtr2Gm8f8y9ehiIgX6O6pUsi2OXhs1CRi2vXwdShVhr1ua1L2rvF1GCLiR6rejwrKKNcG\nwSGhvg6jSrHHdCclfauvwxARPxJAScNJcLCSRmnYOw0gJeiwr8MQET8SMEkjJwhCwiJ8HUaV0rrr\npeyOzOb08d99HYqI+ImASRq5Nqeqp4pRXE+YoRE1aH06nI0rZ1duQCLitwIoaUBwSMX3EFsVjZs/\nrthpnWyNSNm4sBKjERF/pqQhJbLXbaM7qEQkT+AkjSAq5VkU1Y292YWknNrm6zBExE8ETNLICVJJ\noyzsnS4lOfiIr8MQET8REEnD6XDgtEFQcMD8ltFrWnUdwIGIHE4e2e/rUETEDwRE0sjNySLYQZV8\nQFJlKKrvKZfg0DDanopgg+6gEhECJWlkm6QhRXN/1GtROgU1JmWT7qASkQBJGjlZGQT75umw1YI9\nqi0p+9b6OgwR8QMBkTRyc7KUNMrB3vwiUk7v8HUYIuIHAiRpZBPiqKq9wPuevfNlJIfoDioRCZCu\n0VXSKJ8WnftxOCyXE4d2U+f8GF+HI9Wc0+Fg/7Yk0lIWkrptBe1bXcwlQ+/3dVhiCZCkka2kUYKE\nxIQSG8ODgkNofzqS9Stn0WvwXZUXmFRrxw7sIHXtL6RuWUbq/hTSTu4g1XGI1BpniMyx0SazFi2C\nG/DM3k/Z1m8EkXXq+zpkIUCSRk52BsFOVU8VZ9z8cee+gyo4mpTUX5U0xGNOh4MD25PZunExW3es\nYduhVLac2EFazkHSwk+SEeykzZlI2gadT9vacQxtM4Q2LXvQpnM8UdEt8pZzzSONmfzBA9z7yGc+\n3BpxCYikYUoaShrlYY9qR8r+db4OQ/zMqaMH2bZ+EVu3rmDbvg1sPbqVrRn72Go7zvbITGrm2GiZ\nVYMWQQ1oWTOGS2L7cGeLi2jbqT+N4jp59Nupxy97jjvmPcxd2f9VV0B+IGCSRoiSRrnY43ry8/KJ\nvg5DfCTj5DE2rJhJ8vr5rNuziuTT21kXcoTfw3OJOx1GS2ddWoQ3pmW9FgxofyUtW15Ii45/oHaD\nJuVed5+rxnDe3L/y/SdP88c7X/fC1kh5BEbSyFWbRnnZuw4kZe2Lvg5DKlhudhZb1s4jed1c1u1c\nTvLxNNbZDrEjMotWp8PpREM6R7Xj3nZ30bnrFcR1uqTCu+exBQXxeNf7+dvaf/F/jlfVs4OPefJp\nDwdmASeAZ4HuwIvAqgqMy6tyc7IIRiWN8mjWoRcnQh0cO7CDeo2a+zoc8YKcrAzWL53O8tXTWb57\nGSuyd7C+5mkaZYTQKbc+nWu14rq2w3iu06W0u3AQYZG1fBbrNSNf5Il1b7Jg2kT6D3vIZ3GIZ0nj\nWeAr4BLgMuAN4D/AxRUYl1flZGepTaMEJfU95RIUHELH0zVJWTGDPkPuq4SoxJucDgdb1s5j+bLv\nWb7jN5adTmNNzXSaZoTSg6b0bNSd2zo9RKceQ/zyturg0DAeix3Ba4kvK2n4mCdJI9f6PxR4D5iO\nKWlUGbm52SpplOBcd0652EOakJK2iD4oafgzp8PBnrSVrFz6Hcs3L2DZiQ2siDxKzZwgeuQ0pGeD\nLoy7cAQXXjK8SpUab73rnzz3QgPW/fotnS/5o6/DCVieJI09wCTgcuAVIIIq9kty06ahpFFe9vrt\nSDmQ7OswxM3p47+TsuxHkjbMY+3e1SRl7CAp8gQhDuieVZ8edTrwwIX306PPDTRu2cXX4ZZLRK16\nPFhnIG98/zgfK2n4jKdtGoOB14FjQDTweEUG5W25OdmEVK0855fsLS5mxpJFvg4jIDkdDnas/42k\n1TNZu20JScc2kWQ7yM7IbNqfiqRLcDRdGnTk6gtvostFQ2nUopOvQ64QY0a/Q6t/tmLXhqXEdqgy\nNeTViidJIxyYZ72uD2QCv1RYRMZg4B9AMPA+8Gp5FqbqKe+wd72clNXnbv+Q8ju8O42lC6aweNPP\nLDm+nuWRR6iVE0SXnAZ0qdmS69oO4/nOA2l34SBCI2r4OtxKExXdgjts3fnHJ39iwvgVvg4nIHmS\nNFYBzYCj1nAUsN/6uxtY6eWYgoG3gYGYqrHlwFRgQ1kXmJOjhnBviGnXgzPBDg7vTqNBTBtfh1Nt\n5GRlkLJkKotXfM/ivUtZ4tzF/vBsepyOoledjjzc4wF69r2R85t18HWofuHhW/9F18m9eWbftgK/\nHJfK4Uk+AyHGAAAUY0lEQVTSmAN8A7ge3XYFcD3wX8xdVD29HFNPYDOw3Rr+AhhGOZJGbm4Owaqe\nKta5+p5ysQUFYT9di5RVs+inpFFmR/ZuYdHcj1icOpcl6RtYUeMYTTPC6BXUjD4xvXnsotfoePHV\n+vVzMWI7XMzV2S1454MxPPmMnihZ2TxJGr0xJQqXn4AJwD1ARRzVTYFdbsO7KeftvaqeKpknfU+5\n2EObkpL2G/14oGKDqkacDgfrl0xj+rx3mHZoEetqpHPx6fr0rteJsRc/wsX9R1K/SStfh1mlPHbt\na1wxdTiPnDxGRK16vg4noHiSNPYBf8Vc8dswDeMHMNVIFfEQVY9+u52QkJD3Oj4+nvj4+GLfm5ub\no4ZwL7Gf14Hkg7qD6lyyzpxk/vS3mbbyc6bnrCfX5uTq4I480+tx4of+WSe6cup8yR/p/m0DPv3g\nIe566GNfh+O3EhMTmT1jGieO7OX4sQNeWaYnSeNm4Hnge2t4EXATJmkM90oUBe0BYt2GYzGljQLc\nk8a5mJKGkoY32FtezPf75p37jQHo4PYUZv74D6ZtmcnP4XvokFGLq+v35vv4Z+nc5//U/YWXPR7/\nNPcsGMuduR9UeFcm/uzMiSNsWTef1NTFbN6bzM703ezMPMhOjrMzIpPMUCfN6oQRW8s7v+j3ZE8f\nAv5czLTNXomioBVAGyAO2AvciElSZZaTk02wTV9Yb7BfcAUpK5/ydRh+welwsHbh18xY8AHTDy9m\nfeRJBmY2ZWjLwfx76CM0jLP7OsRqrd/VD1Av8Wmmfvos197+N1+HU6GyM06zPWURqRsXkbZ7LalH\n0kjL3Edq6HEOROTS4lQYbZxRtImMoV2DtlzeaCjNmnWmWeuLqN+kVd4Fi+2t8lfTe5I0GmJ+l9ER\niLTGOYFLy732ouVgktRsTGnmA8rRCA6Q61CbhrdEt7qAHJu5qg7Ek+Lh3WnMmTGRWakzmB20jdq5\nwQwO7UDCH56i/5A/EV6zjq9DDBi2oCAe73wfr655m2G3vlwtSnKu3oSTUn5h3d7VbDy9i9Sgo+yo\nkU3TMyG0yalL24gmtG/QjmuajaBth0to1qEXIWERlRajJ0njM+BLTDci9wKjMKWPijTT+vMK3T1V\nMk/6nnKxBQXR6UxtUlbPDoikkZudxYpfPmHWkk+ZdWwF6yNP0v9MIwbHxPPsZZNodUFFXTuJJ669\n9WWeeOItFs14p0o9Etb9x5rrti8j6egm1nGAbTWyaH06nM62xnSOasddbS+nXfs+tOzUz28uSDxJ\nGg0wP7B7EJhv/VWpX9Xk5uYQYgv2dRh+y9M7p1zsYU1J2byYARUTjs/t35rE7JkTmbX1J+aE7qJJ\nVjiDIzvzct8E+gy622++vGJ1ZBhzA6/NfcFvk8bvuzaRsmoWyWmLWHcwmXVZu1lXI5062UF0zqlP\nl5otubrNEJ7udBntul/h98eXJ0kjy/q/H1Pa2Iv5gV+VketQScOb7Od3JPlQiq/D8Kq9aav44ttx\nTDnwM1sizjAwswmD4wby+qA/EdOuh6/DkxLcNnoiz798PuuXTKVjr2t8FseRvVtIWTmTlLTfSDmY\nTErGLlLCj5MZ5MR+pjb2sKZ0Ob8TN7W9l849hlbZ26w9SRovA/WAR4GJQB3gkYoMyttyctUQ7k32\nVr35as8cX4dRbscP7uR/Xybw2ZbvWRV5jGtzWvPqJQn0v/qBSq0jlvKJrFOfB2pdyhvfPsaHlZA0\njh/cScqKmSSnLiTlQDIpZ3aSEnac0yEOOp6phT2kCfYGHbi65U3YL7icJq27V4v2FhdPksY06/8x\nIN56XaWShto0vMvebRApy8fidDiq3Jch68xJZn49ns/WfMLsiN0MONOYMZ3vYMgNTxNZp76vw5My\num/0O7Se2IYXU1fQtO1FXlnm6eO/s375DFI2LSR57xqST28nOeQoR8Ny6Xi6Jp1CmmCv344ru16P\nvevlxLTrUeW+D2VR1pub/wK86c1AKlKuI0clDS9q2LwjwU7Yvy2J6FYX+Dqcc3Lk5rBoxjt8tvDf\nfBO0EXtGHUbGDeWd4eOqbBWBFFS/SStuoytvfXw/r728rFTz5mRlsHH5TJJT5pG8exXJ6VtIDvqd\nPZE5tDsVQaegxnSKasf97QfTqevlNLf/IaB/FxIQW57ryFVDeAk87XvKxRYUhD2jDimrf/LrpLEt\naQHvff0En2Uso3ZuCLfUj2flde/R3N7H16FJBXjklrfp/klfnj64k7oNmxX7vhOHdrNk3icsWj+L\nRcfWsazGUaIzQ+nsbEinOq25pdPNdOp8Ga0vuFTVlEUIkKShkkZJStP3lIs9PJbkLYsZWDEhlcv6\nJVN55ZtH+DFkG6OCujP16il0ueT6gKg6CGTN7X24Mrs57743hsefnpE3fuf6xSxaOIVF2+azKHMz\naTXO0P1UHfrU7sjDF/2Z3vG3qtfmUigpaZyk+H6gqlQH/qYhXCUNb7I3tLN6/xpfh1HAip8nM/7H\nJ1gUtp+H6g7kn/f8UqUeZyrlN3bYq1w5/SYiX7+eRfuWsSh4L1lBDvpkNaZPwwu5pduDdO8/grBI\n73SpEYhKShrVZq/mOnJV0vAye+vefLrrR1+HgdPhYMG0iYyf9yLrQ4/xWPQ1fHL3O9SMaujr0MQH\nuvYbzo0zX2HdoRSubDWIl/4wglZdB6iU6UUBUj2Vq5KGl9m7DyZl6cM+u4PK6XAw84sXGb9sAgeC\nM3ii5c3cMvqffv/DKKl4b/5tla9DqNYCI2k4c5Q0vOy82HZE5gaxJ21lpf74LTc7i/9NfoLxye+Q\na3PyVMe7uX7Ua2qwFKkkgZE0HLmEBClpFKc0fU+5s2fWJWXNT5WSNHKyMvh00p/527bJRDnCeKH7\nowwdOU7VDiKVLCCSRo5DJY2SlPbOKRd7RCzJW5cwyLvhFOB0OJj22bM8uXoCDRzhvNP/FeKHPaxk\nIeIjAZE0Ro+cQG5O1rnfKKVib9SJpXuXV9jyF8+cxOOzx3I0KJNXuz/OkJsTlCxEfCwgkkbjll18\nHUK1ZG/Thw93/uD15W5aPpOnPr+LZcEHeKHVbdw25h2CQyvicfQiUlq6bJMys194JetrnMLp8M6j\n4vdtWcOYxztyybdD6NmgC6nPHuSOBz5UwhDxIwFR0pCKERXdgjrZQezcsLhcXXOcOLSb19++iX9n\nLuKO0O5s+lOa+oQS8VMqaQgJiQllnteeVY/kNWXrJj3rzEn++dofafNGM3ae2suqWxfyxssrlDBE\n/JiShjBu/rgyz2uPbE7K9tL1Kgrw42cJdHg2ipn7FvDTkC/4+I0t6khQpApQ9ZSUi71xZxbu+tXj\n9586epDHXhnAzNxNvN9rPAOvf7wCoxMRb1NJQ8rF3rYPKbn7PXrvip8n0/2lGE7lZrD28a1KGCJV\nkJKGlIv9oqvYWOM0jtycYt+Tm53Fyy9ezlVzRvFC+/uY/MaWEp93ICL+S9VTUi51zo+hflYw25N/\npWXX+LOmb0tawK0fDiWcEFbdt7RS+6kSEe9TSUPK3PeUiz0niuSkgndQOR0OPv7X3fScEs91jfoz\nZ8JBJQyRakAlDSlz31Mu9hrNSdmxgmus4cO70xjz9wFs4Hd+vvYLuvYbXu4YRcQ/qKQh5WZv3IWU\no6kAzPn6Fbq+1Z6Y8IaseGG/EoZINWPzdQBl5HQ6i3sSrVS2ZT/9lztn3cfl4R35Omct/+35Mpff\n8ISvwxKRQmw2G5TzvK/qKSm3jhddyYZFmbRPP8jav2ykQUwbX4ckIhVEJQ3xim1JC4jrdIm6Lhfx\nY94oaegbLuXqe8qlRZd+ShgiAUAlDcE2zobzee1PkepOJQ0REalUShoiIuIxJQ0REfGYkoaIiHjM\nV0njBiAFyAW6F5r2JJAGbASuqOS4AlJ5+54SkcDhq7un2gMO4F3gUWCVNb4jMAXoATQFfgbaWu91\np7unRERKqSrfPbURSC1i/DDgcyAb2A5sBnpWXlgiIlISf2vTaALsdhvejSlxiIiIH6jIvqfmAI2L\nGP8UMK0Uy1E9lIiIn6jIpHF5GebZA8S6DcdY486SkJCQ9zo+Pp74+PgyrE5EpPpKTEwkMTHRq8v0\ndTci84DHgJXWsKshvCf5DeGtObu0oYZwL0pITCj3g5hExP9V5Ybw64BdQC/gR2CmNX498JX1fyZw\nP6qeqnDj5o/zdQgiUkX4uqRRVippeJE6LBQJDFW5pCEiIlWQkoaIiHhMSUNERDympCHqe0pEPKaG\ncBGRAKGGcBERqVRKGiIi4jElDRER8ZiShoiIeExJQ0hITPB1CCJSRejuKVE3IiIBQndPiYhIpVLS\nEBERjylpiIiIx5Q0RETEY0oaor6nRMRjuntKRCRA6O4pERGpVEoaIiLiMSUNERHxmJKGiIh4TElD\n1PeUiHhMd0+J+p4SCRC6e0pERCqVkoaIiHhMSUNERDympCEiIh5T0hD1PSUiHtPdUyIiAUJ3T4mI\nSKVS0hAREY8paYiIiMeUNERExGNKGqK+p0TEY7p7StT3lEiA0N1TIiJSqXyVNF4HNgBrgf8Bdd2m\nPQmkARuBKyo/NBERKY6vksZPgB3oCqRiEgVAR+BG6/9g4N+oNCQi4jd8dUKeAzis10uBGOv1MOBz\nIBvYDmwGelZ2cCIiUjR/uIq/E5hhvW4C7HabthtoWukRBRj1PSUingqpwGXPARoXMf4pYJr1+mkg\nC5hSwnKKvK0nISEh73V8fDzx8fFliVGAhPgEX4cgIhUgMTGRxMREry7Tl7fcjgLuBi4DMqxxT1j/\nX7H+zwKex1RhudMttyIipVSVb7kdDIzFtGFkuI2fCowAwoAWQBtgWaVHJyIiRarI6qmSTMQkhjnW\n8GLgfmA98JX1P8capyKFiIif0C/CRUQCRFWunhI/or6nRMRTKmmI+p4SCRAqaYiISKVS0hAREY8p\naYiIiMeUNERExGNKGqK+p0TEY7p7SkQkQOjuKRERqVRKGiIi4jElDRER8ZiShoiIeExJQ9T3lIh4\nTHdPifqeEgkQuntKvGObrwOoXrz9eM1Ap/3pX5Q0BLb7OoDqRSc579L+9C9KGiIi4jElDRER8VhV\nbQhfA3T1dRAiIlXMWuACXwchIiIiIiIiIiJSjQwGNgJpwF+Lec8/relrgW6lnDfQlGd/bgeSgNXA\nsooLsco4175sDywGMoBHSzlvICrP/tyOjs3CzrU/R2K+40nAIqBLKeb1W8HAZiAOCMU0eHco9J6r\ngBnW64uBJaWYN9CUZ3+C+dlf/YoNscrwZF+eD1wEvETBk5yOzbOVZ3+Cjs3CPNmfvYG61uvBlPHc\n6W+33PbEBL8dyAa+AIYVes81wMfW66VAPaCxh/MGmrLuz0Zu06vqHXbe5sm+PASssKaXdt5AU579\n6aJjM58n+3MxcNx6vRSIKcW8efwtaTQFdrkN77bGefKeJh7MG2jKsz8BnMDPmC/u3RUUY1Xhyb6s\niHmrq/LuEx2bBZV2f44mv4ahVPOGlDHAiuJpr3m6wvBMeffnJcBeTDXBHEyd50IvxFUVladHR/UG\nebby7pM+wD50bLqUZn8OAO7E7MPSzut3JY09QKzbcCwm65X0nhjrPZ7MG2jKuj/3WK/3Wv8PAd9h\nirGBqjzHl47Ns5V3n+yz/uvYNDzdn12A9zDV0kdLOa9fCgG2YBpkwjh3w20v8htzPJk30JRnf9YA\naluva2LutriiAmP1d6U5vhIo2HCrY/Ns5dmfOjbP5sn+bIZpu+hVhnn92pXAJszGPWmNu9f6c3nb\nmr4W6H6OeQNdWfdnS8zBswZIRvsTzr0vG2Pqho9jruJ2ArVKmDfQlXV/6tgs2rn25/vAYcxtyoVv\nVdbxKSIiIiIiIiIiIiIiIiIiIiIiIiIiIhK4csm/13w15kdL/u5C4C0vLetn8n/cdrLQtFHAxBLm\nvQZ41ktxiIhUCeklTLNRvfsmuxT4l9tw4X1xOyUnDRvmR3GhXo5LAoS/9T0lUhZxmF+zfgysw/Sd\nMxbzi9e1mG4oXJ623rsQmEJ+9xSJmNIAwHmY5zWAedbA627LuscaH2/N8zWwAfjUbR09MF1brMF0\nQV3Lev80a3pN4ENr2irM1T+A3Rq32lpX6yK29Wbgh6J3A1AwYa4hvzR2GuiL6ZxuMep2Q0QCSA75\nJ8NvgeaYKitXp3VXAO9ar4MwJ+u+mKSQBERgqnfSgL9Y75tHfhcq7knjHkyiAQgHlmOSVDxwDNMl\nvw34DfgDpu+eLeQnoFqYxBNPftIYj3mKGpjnl2zC9Kf0T0xSANMfUEQR276Bgg8fct8Xq4Ed1nLc\nXQ3Mt+IAuAN4tYhli5yTv3WNLuKJMxR8LG0c5mTp6kvnCutvtTVcE2iDSRT/wzw+NAOY6sG6rgA6\nA9dbw3UwJYBsa32unoDXAC0w1UX7gJXW+MJtDq5lXg08Zg2HY9plFmMSVIwV5+Yi5m0CHHEbLrwv\nbsc87c6lDfAaJmnlWuP2Yp7cJlJqShpSXZwqNPw3YFKhcQ9RsPrG/XUO+dW1ha/w/4x5ZoO7eCDT\nbTgX833y9NkE/4cp6bjbiOlleCim5+F7MSWg0nDfplrAl8BdwAG38UGliFOkALVpSHU0G/OQmZrW\ncFPMw3oWANeSXz011G2e7eRfoV/vNn42cD/5F1htMVVJRXFiqpqi3ZZVm/xqIfdlPug27CoptMBU\ni03EtFt0LmIde4EGxay/sA+B/2LaV9xFY0pmIqWmpCFVUVFXye7j5mAauRdj2jC+wlx1r8Zcea/F\nXMkvJ//K/A3gPkzDdAO35b0PrLfGrwP+Q36Joqg4soEbMSf+NZgEEVHo/S9i7l5KwnTtPc4aP9wa\nXo1pFJ9cxPJ/pWD1U+EYXOtpBvwRkzxd7R2uNpuemAQqIiKl8DwFH+5TFcRjEldZBWGSmaqmpUxU\n0pBAV9Xq9hPJb9Qvi6HAN5g2HBERERERERERERERERERERERERERERHxrv8Hq8ciIJ6tlA8AAAAA\nSUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1140249d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plots = []\n",
    "colors = ['r','g']\n",
    "\n",
    "# Plot lag-frequency spectrum\n",
    "for i in range(0,len(lags)):\n",
    "    plots += plt.plot(cross_spectrums[i].freq, lags[i], colors[i], label=str(energies[i])+'keV')\n",
    "\n",
    "# Draw horizontal and vertical line\n",
    "plt.axvline(v_cutoff, color='g', linestyle='--')\n",
    "plt.axhline(h_cutoff, color='g', linestyle='-.')\n",
    "\n",
    "\n",
    "# Define axis\n",
    "plt.axis([0,0.2,-25,25])\n",
    "plt.legend(plots)\n",
    "plt.xlabel('Frequencies (Hz)')\n",
    "plt.ylabel('Lags')\n",
    "plt.title('Energy Dependent Frequency-lag Spectrum')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected (and also demonstrated in Utley et al), the shape of lag-frequency spectrum for both energy channels is similar.\n"
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}