{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "

Data simulation and cross correlation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this notebook, we are going to start with the stellar simulation from last, where we took a template spectrum and convolved it with a smoothing kernel. Now we are going to finish it by adding noise.\n", "

\n", "After that, we'll simulate a spectrum with a radial velocity, and see if we can recover it using cross-correlation, experimenting with different S/N." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Spectrum simulation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's start with our simulation of a stellar spectrum from last time." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "minimum dispersion: 0.300000 maximum dispersion: 0.300000\n", "using dispersion: 0.300000\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAD4CAYAAADCb7BPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2deXgV1fnHv++9WYAQEiABAgECEsCwIwKKWlRWsYJWq9aKtVbaSluX+rOorbst1VotVm1REXdFqxUFFwSsKLLKviZAgLAGwhIg+z2/P+7MvXPnztxZ7tz9/TxPnsycOXPOeye57zvnfc95DwkhwDAMwzBWccVaAIZhGCYxYQPCMAzD2IINCMMwDGMLNiAMwzCMLdiAMAzDMLZIi7UA0SIvL08UFRXFWgyGYZiEYvXq1UeEEPla11LGgBQVFWHVqlWxFoNhGCahIKLdetfYhcUwDMPYgg0IwzAMYws2IAzDMIwt2IAwDMMwtmADwjAMw9iCDQjDMAxjCzYgDMMwjC3YgDCMDp+s34/jZ+pjLQbDxC1sQBhGg/3Ha/Cbt9bgtje/j7UoDBO3sAFhGA0amjwAgIpjNTGWhGHiFzYgDKOBiwgA4OEdOxlGFzYgDKOBZD/g8bABYRg92IAwjAbyCKTJgRHIbW+uxifr94fdjlVeWrITbyzbjUc/2QwhfY73V1dg/3F2yzHOkDLZeBnGLB6PwDFp9lW49kMIgfkbDmL+hoO4vH9HB6Qzz2PztviO7xhVjMw0N+5+bx06tGqGZfddGlVZmOSERyBMyvBt2RFMffN739u4Hk8t2IYJM74BAByursP7qytw3l8W4kRNg+U+7dzjBLUNTapzDxo93okBB0/WxkIkJgnhEQiTMtz48nJ4BPD3pgHITHMHXf9k/X6cld8SX2w6FFB+93vrAAAbKk7gguK8oPvmrNoLj0fguqFdgq6dqmt0SHpzPPPldjzzZWlQeU19EzLT+X2RcRY2IEzKIMfD6xs9mgbkN2+tAQD0ap+teb8cWFdzz/vrAQDnn5WHdq0y0Szd33a0J3E9/9UOzfKahiZkN+l/3T/dcAAVx2pw60XdIyUak4SwAWFSDqOJVXqGotEjUFldh/zsTM3rFz25GACw4aExyG6WHo6ItklzEbTWzu+oPIWKY2d07/u1tGCSDQhjBR7TMimHUQyEdCzI/723Duc+/iW2HjyJic99q+ue2n7olKIv+3I6yW1vfo9bXnV2S+eHP96Ef/9Pe8TDpAZsQJiUQ63Uj52ux39WV/jOdQYgOFxdBwB47JMtWLf3OL4pPaJZT7n4UMB/7PEINDR5UH7ktD3BDej5x09xpr7JuCKAt5bv0Sy/5l9LA9a+7K06g7V7j2Px1sNBn/eVb8vxl0+34g/vr8feKv3RDZO8sAFhUo4mIQKU+F1z1uL3UqAcAKpOh06g+E2ZV5HWNWor6xkLS1E0bR6EEAHG6tlFZSi+/1OM/NtXOHjC+ZlQ9Y0e03Xv+3CD5iLJleXHcO8HG1A0bR5WlVfhwicWY9Jz3+Lm2Svx05eXY2V5VdA9767a63t+35QewcZ9J7B462F4PAI97/8Ury4tx+GTtZi/4UDAfUIINDaZl5mJPzgGwqQcry4tx7OLyvCP6wZi1je7sK7iRMB1s9NcK47V4OZXVgSVL5He1L8pO4LC1i185U9/ud13vLPyFDrkNLMjviZGbjktTtc3asZq3l21FwDw8brgxY+b9p3AuUVtgsplQ/DTl5f7yt6ZMhz1TR48OHcTHpy7CQDwwg2DMbJXO9z34QYcrq7Ft2VH8eYvhmFIUWvNiQ1MfMMjECblWFV+DABw+ztrg4yHFZ78fBsWb6vUvX6qtlFXsf/kpeWa5Xaxk3HlZG3oKcaNGo0e11nXotX/gRPBK95//eb3mLfhAD5csw/flh0FANzw0nI88N9NJiRm4g02IEzK8d3Oo1Hph4gQrRi6naSP1bUNqG1owrsrteMhWjRJlkJtGLUM5SkdA3W3wl0oM2/DASzaeghX/PMbdmslEOzCYpgIQRS9WVh2+tl2sBrz1x/AjEVlmtdX7z4WVPbsojIQES7t3S6gfF3FCfx+TqBhqGkwF9AHgAt65OGOd9biZG0jqs7Uo112MyzeehgnahowaVAn0+0w0YUNCMNEiLZZGVhSqu/icor9x2vwxGdbLd/3xeZDaJuVoXt968FqzfIZC0vx/OJgo/Of7ysCzr/ffdy0LB4hfC61mvomzFt/AFPf8q5N+WDNPny93fscV/1xFPJaaq/DYaIPu7AYJkL8++udePjjzRHv5/zpi/Dftdaz/c5bfwBpLnsqQCs+ouazTQdNt/fFZn/6mJM1jT7jAcBnPABgymuBa1lO1jZgeZRckkwwbEAYJkLsitB6DyfRm4ocS46F2If++z3H0f+hz+HxCJyqa0T/h77AtTOX4cWvd+KKf34TlESSiSzswmKYCGEU2D5R04Cc5rFJeSJz5FRdTPvX4kx96NlhJ2sb8ef5W/DK0nJf2ePzvanr9x2vwVn5LSMpHqOARyBM0vPxuv3YevBk1Ps1CmxP+8/6oLIZC0vx27fXREiiYKwsPowWR06FXsgJAC99swsdc4PX0VQbTE1mnIUNCJP0/PbtNRj3zJKo92vkwvp0Y2CMQAiBvy/YrrmAL1KEWscSK/74342m6pFG0plFWw+j6nQ9HvhoY5A7a87KvRHJAJDKsAuLYWLIxn0n0LdTDgDgpldWxliaxGKPRv6tGQtLMWOhdz+UkoJWvj1aTtY24B5pxPfUNQPwo3MKoydoEsMjEIaJIZc/+43vWDnbKFXo3UF77xUnePLzbfhsozf/Vo0iyeTv31uH7Ye0pygz1mADwjAxZExJ+1iLEFMiudDy6Ol6/OqN7zH1re9R1xAY69lZeQpbD54MmIX2/uoKVFbH36SCeIYNCMPEEDPrKRKZkoJWuGt0T93rHiHwwOUlEZVh3voDmPXtroCyo6frMe6ZJej1x8+weNthzFt/AHe/tw6TnvsWn200v34l1WEDwjAxpEEn79OFTyyKsiSRoWNuM90dHAFvav0BnXMiLsdsxZRfIDB+8svXVvsWLu47XoNfvbEaOyr9m4Jt3n8Sc9ft10x/n+qYNiBE5CaiNUT0iXTejYiWE1EpEb1LRBlSeaZ0XiZdL1K0ca9Uvo2IxirKx0llZUQ0TVFuuQ+GSST0ptHurQrOZJuIGLmo6hs9aJERei7PH8b1xtZHxzkoFfDv/+30y6BhxC996n/ocd98FE2bh8tmLMHv3l6D/g9/gX3Hk+Pv4hRWRiC3A9iiOP8rgKeFEMUAjgG4RSq/BcAxIUQPAE9L9UBEJQCuA9AHwDgAz0tGyQ3gOQDjAZQAuF6qa7kPhkk0lu+qwmmdrXHjjXF9OjjeZn2jB93yskLWKW7XEs3StfcK+c3FPRyXSUbtXjxV14gR0xdhzZ5jeG5xGV5ftht7jqb2ToymDAgRFQKYAOAl6ZwAXALgfanKqwAmSccTpXNI1y+V6k8E8I4Qok4IsQtAGYCh0k+ZEGKnEKIewDsAJtrsg0kRznl0AX4+OzmmvfZ58POEcI90advCuJIGWivyb7+0GADQvzAX6e7QaijUk2nbMgN9O7WyJZddrnx+KZ78fBv+9N+NuPW1VahtaMI/F5Wm5Cw6syOQZwDcA0Ae67UFcFwIIb86VQCQcy53ArAXAKTrJ6T6vnLVPXrldvoIgIimENEqIlpVWZl6f9xk5ujpeizaejjWYjiG1lax8caAwlzda1cN1k+53kxjp8EBnXMw9zcjMOP6gXC7rL379WjnT1VCAJ65dpCl+51k26Fq/OTFZfjbF9sxedYKjP/HErz+XbmtHSITEUMDQkSXAzgshFitLNaoKgyuOVVu1L+/QIiZQoghQogh+fn5GrcwTHxgZe+MWDGhf4HuNZeOA0AAmvt5uF0u9C/MNYx/AMGbVeUq8ocJwLIBcprv9/jT1m85cBJ/+mgTJj2/NIYSRQ8zI5ARAK4gonJ43UuXwDsiySUi+a9fCEDOv1ABoDMASNdzAFQpy1X36JUfsdEHwwSgXEAWz2h5YD0egTvfXYv1Ffr7avx5/hbda/GAECJAwbdq5v06p1lQ+uo3Q+WjOl3XCLkpt4swoZ++kYsm6/Yex0Nz/dv0Lt52GH9fsB1jn/46qVxdhgZECHGvEKJQCFEEbxB8kRDiBgCLAVwtVbsJwEfS8VzpHNL1RcL7CjEXwHXSDKpuAIoBrACwEkCxNOMqQ+pjrnSP1T6YJEcIgQWKvSPmrNqLy58NznN10ROL8dDcTXjtu3JL7Y+2ubDvwuI8/PYS+wFdLXV6qLoWH67Zhymvrda46mXm1zt1r0WTIovxESujBvU3W5kD65P1B3yjHwIwrq/zgX67yFOHi6bNw82vrMSMhaXYdqgak2etgBACNfVNaGzyxGVKfbOEsw7kDwDuIqIyeOMPL0vlLwNoK5XfBWAaAAghNgGYA2AzgM8ATBVCNEkxjN8A+BzeWV5zpLqW+2CSn6+2V+JWxaZC97y/Hhv3BWfa3VN1BrOXlqPJ4nuFHNy1yogeefj9mF7oajPQrOUCkhWlnf3Oo82153YxVU/edVDPgJzTtbVGaZAF8dEptzlcUlvxOI1G77324Y834+wHPkOP+z9Frz9+lrCjEkvJFIUQXwH4SjreCe8MKnWdWgDX6Nz/OIDHNcrnA5ivUW65Dya5OXZaO9V3k8frKmls8mDpDv8Odev3nrDUvpZiu3F4V7y+bLep+/ViAUbEo/KzgtU4hJ5NlF1coVD21Dorw+fC0srO6xR3j+mJv32x3fJ9z3xZqlmuXtg4edYK/Oz8Itw/4Wxs3HcC/TrlIM1gdlo8EP8SMowCPUXbJE2DnbGoDJNnrfCV66301kPLN3/d0M4aNbWx60nVMjyNHq/s4QaJtd/qnUUp4saHx+pXlNB7TmaUpvJRjevTwf/sKPSU33DI1JhJZgalu9WI2UvLMeiRBbjy+aUY/48luPu9dSiXtgT4fNNBvPZdedylo2cDwiQU+rN9vKpj99HAPTisun9cGsrazJvtMmlf7natgjc5sotsFK2Man4+oltQWYsMe8rPCspn1DJTfxQxoNCbtkTPtZju9rbTtW0LjO3jjUfJVd+ZMhyf33GRr+47U4ZjVEl7KOxHWEwa2FH3WkaaPVV5oqbBUv1T0qLS0sOn8P7qCoz821d4+Ztd+OXrq/HAR5tw3vSFuG7md3hr+R7d0Xg0YQPCJBR660VlJaNWtk0qPXVuUei3cbdG+y4T35LjZ7yKol2IvE+h0Hojb5CEP3jS/Fun1gu8XbeaFUjnGcmfami3NhjZK9+/olzHrrulh00ABnXx/q06tW4OABjevS16dcjGXaN7Ia9lBvp09C4g9AXRw/iYNwzrgnvG9da9btWAFOR4XySsGhAtHv1kMwDvRI3fXVKMQyfrcN+HG3Du41/i57NX4sM1FTHLZsAbSjEJhZ6O0DMg6hXeXdpkYWX5Md32tdxFZkYg7Vt5DUeTzRXlWm/ksgurySPwbdkRjOiRZ9iOW8PaKT9TmosikgFYd2QodTXnl+cBAK7993ch20n3BcQJUy7sjkt7t0Nx+8A9Q4Z2a4NVfxxt2LcVhnZrg8wQRiLDZDyie34WHpvUFy9+vRMHTtT6RhROMKhLa9w5uifuGFWMjftOYu66ffh43QEs2noYzdI34NLe7fHDAR1xce982y43q/AIhEko9HSF7KpS63+1CysjLbSycbsIX9x5UUCZGf0kZ5y9c3RPnF1gPbWGlk5X7mFxw0vLgxIvao1atPTc2r3+dSRWprka5aiSue+y3kHP/ZcXddes++TVA3DVoE4YUtRG8zoppuS6XBRkPLRQjhrtxqDSXK6QsSazI5BOuc1x/ll5EYnFyKMMIkK/whzcP6EES6ddgnenDMePBhdi6Y4j+NUbqzHk0S9x93vr8PX2SjRajAFahQ0Ik1DojQbkL2yQC0ulmdMM/FFuF6FtVoaqz2AuLA4cDch6q2f7bHx6+4Uh+9BCKxeWMqU4ADzyyaaAcy2jo+WCq1L4yq28rb84eYjpuuq/y/CzgjILAfDm0/r7tQONFbKFQYXsPgtnFpbbRZrxL5lQoxN1O4D23yZcDmm4Ml0uwrDubfH4lf2w8v5RmH3zuRjdpz0+23gQk2etwLA/L8Sf/rsRpRHagZFdWExCYTQCURNkQNyhlYyLKEjJasVdHpvUFz948quQbVmh9HDwF/yuOesCzrcdDKyj9Zm1XFhKrEzoslI36BHZVKB2AuJy3XAmq6W5gv/uSsyOQOQXlEisazZao5TmdmFkr3YY2asdahuasHjrYXy8fj/mrNqL/oU5pkZzVmEDwiQUekpCCGDrwZNBezuolaxR5le3i4LuUeqVWy/shsnnFaFzm8AFg+Gqi883GU/3XFl+DM98uR03Du+KNlkZOgYkdBtWklZb+Uxq5Su/ieco8laZQW7FjpzhxEKItEdvMkYjV3+98OMxz98wGLe9+X1A2UM/LLFkAJqluzG+XwHG9yvA8TP1pnKO2YFdWEyCof0FPSZtUfrhmn0B5WpXgtEX3O2ioIC28o4+HXOCjAcQ/t7eq3frB/aVPPNlKc557EvMXlqu2afRCMSKetN6i87WmaKrfqwX9MjD70f3xCMT+5juT55VBViTU1b8Z3e0n9bdI7Rn22VLCxuN7MLk87p6ZZFGuGb+H/SeZbbGYspu+S01apojt0WG7WnIRrABYRIKvZfEMzpJE9UuLDMjkPyWmbhzlH8fb+XbsPqtX/6yX6mRcTaSLNp62NYIxIpm1lKCn94RHN8hUNCIweUi/PbSYuS2yAiqr8WWR8bhw9tG+F1YFuTMykzD27cOtxSzUSOECBrBlE+fgMLW3pcFI3sgj7TSfTEQYwuiV0dd/NYvhuEHPeMzmzi7sJiEQk+v6L14B7uwDEYg5FWGt48qxpvLd+NwdV3A26fSHn1190jkZWeGXDgXKZo8wtYIJJSbJqgPjQ4KW7fAiB5t8W3ZUXTPz8LOSu/CzXA9N82lxY5yINxqQPw8naC9WQS8o9OfDu+CN5bt8ZdLz8BoMaZsQN2+GIh2vXenDMe1M5f5+pS5sDgPS0qP+MpfmjwE6WkutMhw41ydGWvxAI9AmITCaCGhmuAguvEIJKhPhTJTunWK8rJiYjwAYOmOo9ojEAO9ayVO4DGYAapsyakNQe2MQKyitaZDCAEiwmOT+qnKvb+bGxkQ6XeawQhkWHe/oRt1tj/z81TF1rxCCIwqaY8f9MyPa+MBsAFhEoxQQXQtjp0JTPdgJgaiRqnM4ikvruY0XgMDaTIWDMC7piEUWkbjrtE9NWpax65BMhN7+O/UEUFZk/XuG9nb6zpqm+XPMKA1GpHFlacCh5JDnhL85DX98ZNh3izGHXP8zzqe/seMYAPCJBR6ekXofO0OnawLODeKgSjth1ZfeoHPWKAV5DaeBWReMee0SEf59Am619VdlU+fgN/ZTIcv41Req1DIcS6tftXcM7Y3lt93qW+hKACsf3BMUD15ZFcn7SwZKgay9oEx2PLIOGSmufHYxL5Ydu+l6NK2BS6QMg1EI/WMU7ABYRIKPd+4kbtFxmgdiNabr3KB2dg+8bNh0bwNB4LKjGIcTu7+Gpn06eHntTLD8z8djEcn9vElnywpyNGs53YR2qsSZGopeHmx5gfSLEBllUFdAveSb57h9rnEXC5CBylv1iMT+2BMSXsM6xbfbisl8fM6xTBm0FEsZjeOMhqBGHUZarVytLn/w43BhQ7GQMwSCWVvt02zaWTaZTfDjecVAQAe+GFJwLWV94/SHdHqyabeNuDvPx6Iq15YisrqOtNmtnt+S8wMYyZZLOARCBM17v9wA866L2jfMEvoJbUzm8RQbxZW83Q3Zt54jua1BPIoGBqIOLJ/mvhdWPYE7dUhG5sfCb0fidHfMz87E+2y9dPya41S1S8mndu0wKPSGhinJhjEI2xAmKjx5vI9trPVGmF23w+9FcVd2rTAmDhyTxmR11I7bbyRgXBSmYV6S7eLfyW6/TZaZKRh1s+i+yYvy1uiGAEls+GQYRcWk1Do2Ql1plo99EYgifZd1zMURrsXmkm10TIzDT8cUBBUFk3C/XNc0ru9cSUHkf8vrz6n0Ffm26ckqpJEFzYgTNyyYlcVTtY0YFSJXxnovfXe8NJyU22azWkE+N0oQgAv3DDYvxmSDZzeh0PP4Bm99U69uAde+mZXyDpf3vUDX2AXAN7/1Xm+FdkBfUVANZITQxADcltYy89lBaXY8e4udAI2IEzc8mNp86FQU0mtkq7ICTSmpD3OO6stHv54s6bilY2Viwjj+xUEXbfC9UO74PVlu8NqQ4me8g6ltAZ1ydXMsxTUtqoNvb07IuPCitxb+/bHxuNwdW3I+EYonry6P95asUfzmuxCDZhw4cBOifEOGxAmoQg3aWG6QsPOuH4Qyg5799zQ+o7LfTmhAJx+G9WTKVQQnWDOL29VVCd9/ZFciZ6R5tIcSZnlmiGdcc2QzprX/P8ryoVE/sMFd16EFnG0hsgpku8TMXGPnDbC1r1h9q2MEQR81zXEkftyQpc5HVDVzQkWopt0t8vcZzGoFIGtLqyKEHfIo7FAF5Y/r1ck9uKIB3gWFhN1wgkFhLtRj3Idh5EPXzhoQZx+o9YzSFojkEcn9QXgzVob7+4UfwgkzgVV4RuBKMp8/2qJ9VEswQaEiTpmp9xqEe7Lr1LBEoV2U/38giIAQKtm4QddnV7AZyWI3koR9zDnwjJnWCMTRPe2mWgBaN//peL5JlJKEruwC4uJOk0eAdsTmsK0IOrvtM/1oKEMbxvZA7eN7BFUbqtfR1pRtKcbA/EfvzNlOJaWHUFWhHaj61+Yg22HqtErAu6ZyKRJiRxaI5Bo5PWKNWxAmKgTjhfKyZk/FKYslvpy2oWlNwtLYUF6d8jG8O5t8eXmQ9I9znLloE6YenEPFOVlOdbmyZoGAMCK8irH2owOwTGQRDOCdmAXFhN1zOatknl83mbsO14DwFmFT0ToXZCNET3a4s9X9jO+Icy+nG1Pu1yZTDEaCsxJ4wEAWw5WO9petNBy60VimnO8wQaEiTpWYyAvLtmF3729BoDDBgRAZpobb/5iOPoVamdjdbIvJ9HzrwesRLfZaSxd9zZyXcacdDf5/i8DYjcOTgOPVxLwz8UkCnurzuAfX5YGzZxat/e45bbkVCXhBODVxNsX22i7XSV6NY2mJjtBJN+sE8Xtc1a+d+S1/qExWPvAGM1pvP5JfInxmezAMRAmYtz62ipsPViNSYM6omtbv6tj5tc7cWFxvqW25C9oOKorIy3wfSnepop6RxXmPqGe6F3a+BfKRfzTRaCDOPuT6PLh1BGoOlXvm6EXyZlp8QyPQJiIUevbnS2wXJ2Rd9nOo9gvxTjW7T2OZxeWBrUlf0HDWQcSzaC5ZucA+nXKQU5z7WnBVqZ96hm/bMWU40gZyJg9wziiVbP0gPiPxz/c8JEKpoQNCBMxZAWmdjupkwpeN3MZLnnqKwDAxOe+xVMLtuNUXWNAHb8BsS9PPMzL71eYg+752oHnGsngAsC943vbal8rBGL1kaWbTDgZkXUgjrcYHfzTwVMLw/8UImpGRCuIaB0RbSKih6XybkS0nIhKiehdIsqQyjOl8zLpepGirXul8m1ENFZRPk4qKyOiaYpyy30w8YOsr9VKX2sUUdsQmI5dnZ5dviOcVezehYPReX1ecd+l2jKYvL9/oX8b1C2PjAu6LufwUqNeKGmVHu1aIscgW62TOcKCiAMjbwutXFgpgJlXjToAlwghBgAYCGAcEQ0H8FcATwshigEcA3CLVP8WAMeEED0APC3VAxGVALgOQB8A4wA8T0RuInIDeA7AeAAlAK6X6sJqH0x8ISsztdKWT78pPYIvNh3UvLdRtcm53EY4AdynrhkQlgGyQrtW2hlfQ+mXZ68f5DvOz/ZvGCXvn20GrXUIZlXav346GJ/dfqFhvT9dXoIBnXMxoDDXsK5VIqF+L+9fgIIcexl4zRIq600yT+c1DKIL7zdXft1Jl34EgEsA/EQqfxXAQwBeADBROgaA9wH8k7xmeSKAd4QQdQB2EVEZgKFSvTIhxE4AIKJ3AEwkoi1W+xDRer1kTCG7U4JiINKf6acvB+7hsXHfCd9xY5P2nzIcAzC+X4Ev1hIrQrl9+nT07maXmebyGYJOuc0ttR/OCKR/YS7STMyj7VeYg4+mjrDWeAz5508GR7wPWfVoPfNkDqybcnZKI4W1AA4DWABgB4DjQgjZUV0BoJN03AnAXgCQrp8A0FZZrrpHr7ytjT7Uck8holVEtKqystLMR2UcRP7iqGMgQmi7ki5/9hvf8bz1B4Lu8f4O7x2ho0WF7BS+0UAIXaJ0fzgZrzH7xMzsFRJpEtUDJD9j5d8tFd5mTRkQIUSTEGIggEJ4Rw1na1WTfmuP4pwrD9VHYIEQM4UQQ4QQQ/LzrU0bZcJHGQPZp3jzX7v3OLrdOz/kve+vrgg4903jTYJvpd5nIJ1jK4Qagegp52nje2PVH0cFzOBKBi7okRe1viIaF4pjLL1yCCGOE9FXAIYDyCWiNGkEUAhgv1StAkBnABVElAYgB0CVolxGeY9W+REbfTBxhEsxC8uq60jtN95+6BSqTtfbWkg4vHsbXDWo0LhiFAilX7QWAVpVSIGzsAJv1nt0aS5CXstM7YsJzKyfnYu6xibjig7g5ALXRMLMLKx8IsqVjpsDGAVgC4DFAK6Wqt0E4CPpeK50Dun6Iik2MRfAddIMqm4AigGsALASQLE04yoD3kD7XOkeq30wcYRyBOLEi9kf/rPe1ghkZK92+PG52jvJRZtQs3SUCt+My8uofaNssOP7drDWeBRw8gU+I80VtVGVbEDSFFOg5cD98O5B3vWkwcwIpADAq9JsKReAOUKIT4hoM4B3iOgxAGsAvCzVfxnA61KQvApegwAhxCYimgNgM4BGAFOFEE0AQES/AfA5ADeAWUKITVJbf7DSBxNfVJ2ut32vlqGorm2w9aYXD16FUIbgk99egLyWmWho8vjq2p25Y2UdSHudmWKxJIhGwz4AACAASURBVFGnwWZneg1VC8WMue75LfG//xuJzmFsoxvvmJmFtR7AII3ynfDPolKW1wK4RqetxwE8rlE+H0CQU9xOH0x8cN+HG3DgRC0ArzK0qhcuKM5DqWqtgxD2ApPxsIAwFG2yMtAhpxkqjp0BIE8y8F6zOoMncARCqmva98T780kE/vTDEvTskI2RvQJjrcoUPskIr0RnHKW+0YNbZq/EW8v3+Mq8U2+tKSmttzYBe7Ow1Prx+qHRd2cpg6zqT+CPdzic8t1QJv2pp7EijkSxRMvMNNxyQbeEHUHZhQ0I4yibD5zEwq2HA8rsKH3lmhB/Q/ZmYam/1H+5qj/Kp0+w3lAYhNr5UMZp1WOkyxzc8t0xUkz/JjxsQBhH0TIW936wwbJi+GDNvuC2IWwtJIwLnRRimqdWwNyn3MMQXjacegbcPyqKiycEILkX3SUjbECYiLP1YDWqaxuNKxrgjYGE78KKNXriBOxmJ2l3Z+MTqmm9GntYxJx4koUxhA0I4yh66t2JefIC9lKZxFOQmKD/jJRiehyMT+j15w/Uxw/xJAtjDBsQxlH07IQTy3Q8QjgSRI8FoaSWt6ENWIkuCZ2dGX56Ea1V0mP7tMfQbm0AAGcXtAq7D6eIh78VY57YJ79hUgInlnkKh4LoscDnklIs1Mhulobq2kZ/tl3F6vPueVm4Z1wvXDmok7opO70rmw+YQDC8e9u4XA/CJAZsQBiH0dbw6l0I7bacSAsJv7v3Etz7wQZ8ta0SmWnuIFmev2Ewmqe70VIaZQSsRCfCbSN7OCqPlh1l48GEA7uwGEfRsxNO5JkRQtgagTQ2eYwrRYCCnOb410/PwWd3XIiWiky38mgku1k6hhS18ZWHGiiNOrud7la4RiRSkh+ehZVYsAFhHEVPWXmcGIEIeyOQd1buNa4UIZqlu9G7Qyv/cyH/Z3CpdKV/A67gdl666Vyse3CMLRn8XbNyZpyFDQjjKHoKvt6BUYDd/FAnahrC7jtclLLrpSkxo95fvyUos49x3wmUajwRZGT8sAFhHEVvgNCgs8OgFTweeyMQJ+Iv4aI0GnoK3UwK9wuL83HVYHuB9URQzokgI+OHDQjjKA06Iw0n4hDeXFjW78tMj59/c2UurCAD4tvDPLQWteqKSuY9uZnYEj/fLCbh2bz/JCbPWqF5rcGBUUBdY5OthYT9C3PD7tspCN4dANtmZaB7XsvgiyZQx04Ab0bfZ68PSpoNwH5mX4YxgqfxMo7x/Z5jutf2HD0ddvs7K08n1DReJcoFkD/omY/VfxodVMfsLoRa1+8YVYwfDuioWb9vpxwAwIT+BeaEZRiTsAFhosKLS3Y50k6jjVhKPCwklAklCql+66GVmiWUXe2Wl4Vdf7ksrp6DHok05ZhhFxaTYDz95XbL94wuaR8BSaxhxo3kj40YxEA0rhuNzBLBeDCJBxsQxjHi9eXxCh3XTjQxk57dbHJDrTbiYKIZk4KwAWGcg/0Pulh6NAYWRCuI7kSyynggST5GysAGhHEMK2/Bg7o4NzPqoR+WaJa3zcrAoxP7ONZPOAhVQkMtZMPQJisjZFtabjAn0uXHAzzlOLHgIDrjGFaU2LVDOmPNnuOO9Nu1bZZm+QXFebjxvCJH+ggX2bi6tIYPEtnN0vH4lX3xg575IdvSaiKe9jxhUgc2IIxjWHkJdodQpFZJBN0p5wIzUvQ3DOtq2FZ1nX93Rxd5jVM/aapuopMkA6mUgQ0I4xhWRiBpbr8iTXdTWKlO9JRyrw7Zttt0GvnZOGE4P/jev1/8+L4FePCKErTL5rTsTPThGAjjGFYMiNKPf+25ncPqV8t+fHDb+fjVRWeF1a6TjO3TAQBwca92jrfNxoOJFTwCYRzDShBdqfTdNn1QGWku1Dd6NIPKg7u0ttVmpBjQOTdgJ0CGSQbYgDABXPjEIrTJysRHU0dYvtfuTCC7i9w++PX5OHiiVjerbbJClLyxgiT9WEkLu7CYAPZW1WDdXnuzo6woNaXRsLuGoUNOM4wqaR80/lhw5w9stZcopMKMq4w0Vk2JAP+VGMewkqdKqQLtrqKmoAPpNMn1q4MT2OKWl28aEmsRGBOwAWEco9Fjb88Pu4vH5FGMOgaS7G/oSZ3XSvpXSHezakoEOAbCOIaVbWuVJiNcf37wxkzJTcAIROPDvnXrMOysDD99PsMYwQaEcQwzLqy2WRk4errekdxNeoYimV/QAdWIS+Mxnn9WHs4/Ky96AjEpC48TU5ANFSdQNG0elpYdCSgPV6m/u3KvYZ2Zk89B+fQJAaMOu736NmBSlyf5GCTZDSSTOLABSUGW7TwKAFi49XBAeWOYOcFPKVJs6KFlo+zaLd8e4iqNmuwKlnRPkgDp8yTrNOVkw9CAEFFnIlpMRFuIaBMR3S6VtyGiBURUKv1uLZUTEc0gojIiWk9EgxVt3STVLyWimxTl5xDRBumeGSRpBDt9pBLLdx5FXWOTqbq1DU245KmvUHqoWlfBRjOjqzJwHu7IJ/XWgST5B2QSBjMjkEYAvxdCnA1gOICpRFQCYBqAhUKIYgALpXMAGA+gWPqZAuAFwGsMADwIYBiAoQAelA2CVGeK4r5xUrmlPlKJ7Yeqce3MZXjsky2m6v989krsrDyNm2atwI7KUwCC3/JsTqKyhNylsu+mMOfxBrmwklzBJvenYxIJQwMihDgghPheOq4GsAVAJwATAbwqVXsVwCTpeCKA14SXZQByiagAwFgAC4QQVUKIYwAWABgnXWslhPhOeF9FX1O1ZaWPlOHY6XoAwNaDJ03VX7rD67baf6IWb6/wxirU6wmaojACkbtQdmXbflDgb1+5veYShiS3j0wCYSkGQkRFAAYBWA6gvRDiAOA1MgDkLHGdACijqRVSWajyCo1y2OgjZbD6lq2VBVa9N0V0XVh+uudr7+dhnhSLgST7B2QSBtMGhIhaAvgPgDuEEKFee7X+u4WN8pDimLmHiKYQ0SoiWlVZWWnQZGJiVufLGW8nn+ffb0KthzxR2Fhb7lNprPp0bIUv7rzIeluqNv3lya1gU8F+pMJnTAZMGRAiSofXeLwphPhAKj4ku42k3/KUngoAyvzchQD2G5QXapTb6SMAIcRMIcQQIcSQ/PzQu7zFgle+3YWDJ2pt3St/wcyq/DYtvNuk5rfM9JWpV2xHwX5oqnYiQnsbKcn9K9EDSfZUH0n+8ZgEwswsLALwMoAtQoi/Ky7NBSDPpLoJwEeK8snSTKnhAE5I7qfPAYwhotZS8HwMgM+la9VENFzqa7KqLSt9JAwVx87g4Y834xevrbR1v1UlopUuJCgGYsGCrNhVheNn6i1KoXizdMBY+UcgqRUEYRcWEy+YWYk+AsCNADYQ0Vqp7D4A0wHMIaJbAOwBcI10bT6AywCUATgD4GYAEEJUEdGjAGSN+YgQoko6/jWA2QCaA/hU+oHVPhIJecbT8TMNEe3nsU82g8ifW+iVpeW+a++urMBtI3sgK9P7b2A2BtLQ5MGP//0d+nRshXm/u9CiRF7lpzRo4arDVFtImMx0zPGORFtkuGMsCWMGQwMihPgG+t/xSzXqCwBTddqaBWCWRvkqAH01yo9a7SNRkF8i5dlUdjFaQ/HSN7sAALeN9O7OV6Xo78ipOjw0dxOevGYAAPMGZPa35QCATfvNzQBT4nO9qbqyk75bdsGl3DqQWAsQQR6/sh8u7t0O/QtzYy0KYwJeiR5jTtebWwioxqqS1DMNlafqfMdmXVjKe6zSuXULb18KC0IENM9wY/pV/Sy15U9lopqFZVu6xCCZDWRWZhomDkypCZUJDSdTjDLrK46jvtGDgtzmjrQXbihBaTSisZAwP9sbxNeyVbkt0m21qVaoyZ7OPflNJJMosAGJMlf881sAwNJpl0S0n8v+scSUQla6kuau2xdBiQJRThmWRxBWl6Ho2Ylktx/Kz5fkH5WJc9iAxIhwlZyRst18IDA+oVdfGff42xfbTfWtNkyb9p9AhsUNgLTcZVZnF/mTKWqXJyvKT8c5B5lYwgYkQdHKKRW6vnZFq6vPT9Q0+BR0lzbeeMaEGd8A8AbCrzu3M177brdhOx5VDMQOejGQJLcfKeCiYxIFNiAxItxFe5YTEOqOQKw1c+2/v8PWg9UAgDR3oCJzE6FZurnpl1ryW83Kq7sSPcn1a7J/PiZx4FlYMSLcFOZWRw5O5bmSjQegkc1XCNPKTWk/SKPMDKQzjded5BrWYEdbhokabEBiRLj63JfV1mR93e3Kw5BDPYoQwrx75UZFTi5fe9KHym5mbmDsG4Go1Giyu3h4JToTL7ABiRHhGhCPVl50M/XVcoRhQZo8ImAk5RHCdB6qlplpGN69TUBZYWvv1OaBnc0tIkvVWVgMEy+wAYkR4biUPlxTgc82HgRgfgARrstMC48QKD96JuDc1tu/dMvgLq0x73cX4NYLuwPwGpS1D4zWv02nLzYgDBMd2IDEiHDU+Z3vrsOby/d425Ea+ueiUmxTxCfUWN0s6tmFpahtCL1KPngEYs29oiVSn445AUYoV8oibIVkn8bLMPECG5AYYWdEcPxMPU7UBCZfFBCob/Tgb19sx1XPf6t7r9UA9VMLtuO5xWUh62iNOOykUg9KRRKm/tfaPIthGOfhabwxws403oGPLAgqEwI4cKIGQOi8WnoGK5QdO1kTOlOw1mew4sLS61ov4aIR3fKyMOeX5yW9AYmEO5Jh7MAjkBjx5ZZDvuNwFIIQwB3vrjWsZ2bdSOc21vJzeURwCN5FwCMT++CNW4YZNyDd7NRK8nQ3+XJtMQwTediAxIjpn271HVccq7HdjoA3LbsR2w+dMqxT0CrQgBjFM9QxEPmeyecV4YLiPMP+ZPMTtJ+HCfux8eGxxpWSlOd/ek6sRWAYAGxA4oKlO47YvlcIgb1VxgZo7d7j2vcrjtUzw2YrNp/SwuMRQW6sNAvuI+EbgVhPx94yM3W9rz3bt/Qd85oQJpawAYkDnHRph+MOszq1uEkIHK4O3NPdSvxB7i04FYk9pZgqoQGeZcbEC2xA4oBw9N5W1dRdq0pUPQ3XCh4P8H/vrQ8oS7eQlVfuW8+FpbYjetucptpLOKdzZ+IFNiBJRjgLFK2OXpqEwOn6xoAyKyOQzDSvQXC5zLmwLCeQTFKcyGTMME6Quo7kOMKpRIeAV6nb/aNa1c/eIHpgWbrbvEZ7+tqBeH1ZOQaq9r+WXVjqtlPFRWVEsuf6YhIHHoHEAc7GQOzf64Qhc7vM/0t1yGmG/xvbO3gEoqMfnTS0iUyzdDcevqIPAKBPx1YxloZJZXgEEgf89bOtWL6rCtOv6odTdY1o36qZ7bb2VJ1Bz/bZpusHzsKy3p/a7WVxY0JN9N6v2YD4uen8IpzTtTVKCtiAMLGDRyBxQHVtIz5etx9XPv8thv15YVhtjXn6a0v1lTrZzgwu9R2hRiD3XdbbVJt6s7A4BBJI3045QaM3hokmbEDiCL3Ffoera33pSiKJnTf8X17UPeBcuZnT9sfGB1y7+pzOptpklcgwiQG7sOKYQydr0S47E0MfD29UYoYVu6pQWW28ol2NekGf0oWVkRb4fmLWMOjFQPKzM23JyDBMZOARSIT4ensljp2ut33/jspTGPbnhXhpyS7L99pxRf3439/h2JnQyRO1ULuVrATR9dCbZfTBr88PeR97uBgmurABiQA19U2YPGsFfvbKCttt7K3ybtT0dWml5Xut2I9wla7a7RUqiB7u7NPObVpotxteswzD2IQNSASQN28qPWycwFAP+S3cTlxiR6WFfkVwQkQrqG91Yo2Cuo3cFumhZQi7R4Zh7MAxEAfYUXkK7bIzkd3Mq+hk9Re8EM68qvOl87Dxfv3zV1dqlhe1bRGwBa1MOLOb1AndndiLQ22DFv1+JI6f0XcHyo+VJyQxTHRhA+IAlz71P5QUtML82y8EoD9qqDIZE9lbdcZnOHYdOW1ZnsYm7f475DRDVmYaNu0/GVAeToqQJk/guTvECMSsMXT5VqJ75WqTlYE2Wd6tbWf9bAjyWwauk5GD9e2y7a+fYRjGOmxAHGLzAb9S1tPHf1+w3VRbR07VQY5F7ztuffqu3kDHRaTpYgpngV5wDETbSBABLTK1kyFq1dXjkt7tg8q65WXhiR/1x+iS4GsMw0QOjoFEAPnNuaahCYu3HvaVv7l8j6n7m6W7w4ol6BmEhiYPNuw7EVAmYG0E0kG1Sn730cARkp4B2fWXCaYz9cpNWDFrPz63M1pLoxSGYaIDG5Aw0YprKHcYfHz+FsttprtdYea00i6va/RoljcZdPbbt9fgjWW7AcDnSpKZs6oi4NyZldEczGCYRIBdWGGi9fb+yMebQ143QggRkbTsdQ3aBsRjIOPH6/bj43X78eKSnditEYRXEioGYhYOhjNMYmA4AiGiWUR0mIg2KsraENECIiqVfreWyomIZhBRGRGtJ6LBintukuqXEtFNivJziGiDdM8MkhIh2ekjFjRqKV+FAtx15HRQHKNTbnOEwiOsG56hRW0U95u/V1joy8h4AE7NwmILwjCJgBkX1mwA41Rl0wAsFEIUA1gonQPAeADF0s8UAC8AXmMA4EEAwwAMBfCgbBCkOlMU942z00esUAbPZdQjgAv+uijgXG9nPRmPEIZuJTVpin04LO/r4WCWWycMiC8Gwgs8GCauMTQgQoivAVSpiicCeFU6fhXAJEX5a8LLMgC5RFQAYCyABUKIKiHEMQALAIyTrrUSQnwnvFr3NVVbVvqICSt2+R9NTX0TiqbNw8ryYwF1rC628whh6FZSo1TclkYgEPBoe7Zs4cgIhGMgDJMQ2A2itxdCHAAA6Xc7qbwTgL2KehVSWajyCo1yO30EQURTiGgVEa2qrLSeEsSI2oamAFV3uLrW1H1GHhorbiWZDMUMJ6vGJ9QIJNQCPi2cWIlOvllYPARhmHjG6SC6lvYQNsrt9BFcKMRMADMBYMiQIY5ro95/+izgPM2J3ZQgjUAs+m9qG5t8x3oGQU8hhzI46ypO6F7TIi0CK9EZholP7BqQQ0RUIIQ4ILmP5MUOFQCUmz4UAtgvlY9UlX8llRdq1LfTR8wZMX2RcSUTNDR5dKfc6nGqzm9ArLikNu47qbtC3kVAu+xMS3KoXViPTOwTtHbECJfOnugMw8QXdl+Z5wKQZ1LdBOAjRflkaabUcAAnJPfT5wDGEFFrKXg+BsDn0rVqIhouzb6arGrLSh8Jw/VDu4S8/qMXvsPt76y11GamYvTTaDGoMfG5bzXL7cwGU68DmXxeEcb06WCpDR6BMExiYGYa79sAvgPQi4gqiOgWANMBjCaiUgCjpXMAmA9gJ4AyAC8CuA0AhBBVAB4FsFL6eUQqA4BfA3hJumcHgE+lckt9JBJXDdYM2YSHQuk6ufWr5jTlEDixDoSD6AyTGBi6sIQQ1+tculSjrgAwVaedWQBmaZSvAtBXo/yo1T6iiZ0khzJOzFRSEymV26DOlihx1aBO+GDNvqDySGTjZRgmPuFUJjY4VdeIi//2laV7ztQ3+o5dRHj8yr4YdXa7EHdYI1JKt0EnFqOXGNFJ48ghEIaJb9iA2KDeYoAbAPZW+Veju12EG4Z1xb2Xne2kWBGhXmcEoufZcsaFxTBMIsAGxAZ2XrL3K9KZuFX7XUQLO9016OwtomdEHdgS3Tfy4FlYDBPfsAGxgZ0g9c2z/bsEyjOVnAx2m1G2RMATV/cPWWdg59yAc3W6dpma+ibNcidcWLzDIMMkBmxAbBDODn5KYvGGXWiQyLFX++yA88fmaaej//2Ynpr3OGFA5EWUkZhswDCMc7ABsYAQIuxU6wHtORgmvmJgR3MVDXSy2c/WPb8lSgpaBd3jRAxE3qK2WGXMGIaJL9iAWKDbvfNx62urHTMgTo1kAKCkoBX+dHlJyDpm1ldYkUiOdwzt5k8l78SoIa9lJt64ZRj++ZNBYbfFMEzkYANikS+3HIpLF5bbRY7EDKzIJBukvp1y/GUOzSe+oDgPrZqlO9IWwzCRgQ2IDZxKf96nYyvf8Y8GF4aoaYyLzK3fNjJ+AgJXDTK3Ul42WE64rRiGSTx4S1sbXPTkYkfaISLs/PNlqGloQlZmGv7zfYXxTTq4XWT49k8EnNO1NQZ3ycX3e45r1lmw6RA2PDwWFcdqsKJcvQ2MRoNwah90hmESDR6BRJmv/+/igHOXi5CVGb4dTzPpwmqRkYYPbhuhe726zrti3swuhb4RCP8XMUxKwl/9KJPTPDy/fu8O2jOTXC4KSOlul7F92gMAitu1NKwr2yu3E6sHGYZJOPibH2VC6drpV/UzvP/i3tr5s9xEeHP5brti+XjhhnMAAPeOD06z0iYrI+Bc3reDYyAMk5pwDMQBiMzPXkoLYUHOVUyH1UOprLMz03wuJ+8srPAVuRzPyEwPltMjBO67rDfOPyvPW1c2IPwawjApCX/1TRIqb1Wn3OYYdXZ7U+2EGoGYMQDKgHWn1s2RLcVPnA5ka63naPIITLnoLN+0XVlcJwwXwzCJBxsQkxw4Uat7ze0iXNbP3K57oUYgZlxByj3HR53d3heIcBMZpnS3skZDSxb1FGDyBdHZgDBMKsIGxCRPfLZV95qLCFeZXMcRStcqdfbPzi/Crr9cFlxH+n3NOYW4a7Q/H5XbZW0fv04aObFm3niOX04NQdW7E/pdWGxAGCYVYQNiklBuGll/tsjQ3mRJSahRgFpph6rbrlUmXC7yxV4y0lyWRhhL7rk4qKynQe4p9QhEfibswmKY1IQNiAPICjTcHFlmXuRlIyN3JfeZ5iLcb7BBldL9pTXC0BtJyCMT9eeT7QZv28EwqQkbEAeQFW+4KU6sTIeVlfYdo4oBAJlpLgxQ7eUhI6dJ+cmwLiHb1OteniDQs13gCIUcMpwMwyQmPI3XJKGUpNYI5OyCVthy4KSlPqy4oOSuplx0FqZcdJYkR2CdaeN7Y/qnW3H90M546scDDNtUj0Bm33wuNu0/CZeL8MYtw9C7QGVAfLKwAWGYVIRHICYJpSLliVWjS/xTeR+b1Md3PLJXPoaZWONhxoXl8cibLWlcUwk55cLu+OS3F2BIkXHf3v4DBRjZqx2mXtwDgDc7bl7LzIDrcndOZeBlGCax4BGIAQ1NHhTf/2noOo1eVfrMdQPxwOl6FOQEznCaffNQNHmEYSZcM8Hoi3u3w1MLtmNMSfC04byWGXCR35C4XBSQat0Iu3aAzQfDpCY8AjHg2Ol6wzrbDlUDADLT3EHGQ8btIt9Oe3poGZDJ53XFvxXTa/t2ykH59Ama8Q4iwn0GgfRQZLqNZ5EpUbqunrl2IB6Z2CdEbYZhkg0egcQRyjWGzdK9yvyRiX0ttTGsW1sAwCQTW9xO6F+AeesPAAAendQXOS2sJXrsnpeFJaVH0LpFBkb20s7RxTBM8sIGxICGMHYf/NUPzkJ+dqZxRQnlCOQGgxlTevQr9I5QzPDEj/rjD2N7AwC6tG1hua/7JpyNkb3b6c7+YhgmuWEDYkBDo/25udPG97ZUXzYg6W5C5zaBCj3dTZjQr8C2LFpkZaaFtRdJZpobF/PIg2FSFjYgBjQ0ObR/rQW0YiGljwenNWEYhoklHEQ3oC7ECKSwtXbA3C7pbq/huOn8IkfbZRiGiQQ8AjFAPQKZ0K8A8zYcQIsMN/p0bIWKYzWO9ZXmdqHs8fGcnJBhmISARyAGNDQFBtH/enV/AN6V4E/8yHh1t1XS3NaSIjIMw8QKHoEYcLq+MeA8K8ONm87riqsGFyKnRTpG9spHt7ysGEnHMAwTO9iAGHD0VOBCQiLCw4q1GbNvHhptkRiGYeKChHVhEdE4ItpGRGVENC1S/TQ0edC5jbPBcoZhmGQgIUcgROQG8ByA0QAqAKwkorlCiM1O93X90C64fmgXvLl8N3p3CL3hEsMwTCqRkAYEwFAAZUKInQBARO8AmAjAcQMic8OwrpFqmmEYJiFJVBdWJwB7FecVUlkARDSFiFYR0arKysqoCccwDJMKJKoB0ZrnGpS0SggxUwgxRAgxJD8/PwpiMQzDpA6JakAqAHRWnBcC2B8jWRiGYVKSRDUgKwEUE1E3IsoAcB2AuTGWiWEYJqVIyCC6EKKRiH4D4HMAbgCzhBCbYiwWwzBMSpGQBgQAhBDzAcyPtRwMwzCpSqK6sBiGYZgYwwaEYRiGsQUJYX/L1kSCiCoB7LZ5ex6AIw6K4xQsl3XiVTaWyxoslzXCkaurEEJzHUTKGJBwIKJVQoghsZZDDctlnXiVjeWyBstljUjJxS4shmEYxhZsQBiGYRhbsAExx8xYC6ADy2WdeJWN5bIGy2WNiMjFMRCGYRjGFjwCYRiGYWzBBoRhGIaxRcoaECJqRkQriGgdEW0iooel8m5EtJyISonoXSlZI4goUzovk64XKdq6VyrfRkRjIyTXbCLaRURrpZ+BUjkR0Qyp//VENFjR1k3S5yglopvCkUtqz01Ea4joE+k8ps8qhFwxf1ZSm+VEtEGSYZVU1oaIFkj9LCCi1tGWTUeuh4hon+KZXaaor/k3I4e3lSaiXCJ6n4i2EtEWIjovTp6XllwxfV5E1EvR91oiOklEd0T9eQkhUvIH3j1FWkrH6QCWAxgOYA6A66TyfwH4tXR8G4B/ScfXAXhXOi4BsA5AJoBuAHYAcEdArtkArtaofxmAT6X7hgNYLpW3AbBT+t1aOm4d5jO7C8BbAD6RzmP6rELIFfNnJbVbDiBPVfYEgGnS8TQAf43B31FLrocA3K1RV/NvJv3sANAdQIZUpyRMuV4F8AvpOANAbpw8Ly25Yv68FH26ARwE0DXazytlRyDCyynpNF36EQAuAfC+VP4qgEnS8UTpHNL1S4mIpPJ3hBB1QohdAMrg3XLXabn0mAjgNem+ZQByiagAwFgAC4QQVUKIYwAWABhnVy4iKgQwAcBL0jkhxs9KOxaG0gAAA29JREFUSy4DovKsTMggPxv1M4u1bHryav3NfNtKCyHqAcjbStuCiFoBuAjAywAghKgXQhxHjJ9XCLn0iMrzUnEpgB1CiN2I8vNKWQMC+FwfawEchvfB7QBwXAjRKFVRbpXr20ZXun4CQFuY3F43HLmEEMulS49Lw8+niShTLZeqf6flegbAPQA80nlbxMGz0pBLJpbPSkYA+IKIVhPRFKmsvRDiAABIv9vFQDYtuQDgN9IzmyW7PqIoV3cAlQBeIa878iUiykLsn5eeXEBsn5eS6wC8LR1H9XmltAERQjQJIQbCu6PhUABna1WTfutto2tqe91w5CKivgDuBdAbwLnwDjf/EC25iOhyAIeFEKuVxSHaj8qz0pELiOGzUjFCCDEYwHgAU4noohB1oymbllwvADgLwEAABwA8FWW50gAMBvCCEGIQgNPwumD0iLVcsX5e3s68cccrALxnVDUScqW0AZGRhqRfwesbzCUieZ8U5Va5vm10pes5AKoQwe11FXKNE0IckIafdQBegd/1o9e/k3KNAHAFEZXDO/S+BN43/1g/qyC5iOiNGD8rH0KI/dLvwwA+lOQ4JLkOIP0+HG3ZtOQSQhySXlw8AF5E9J9ZBYAKxWj7fXgVd6yfl6ZccfC8ZMYD+F4IcUg6j+7zMhssSbYfAPkAcqXj5gCWALgcXkuuDAzfJh1PRWBgeI503AeBQbOdCC+IridXgVRG8Crv6dL5BAQGx1YIf3BsF7yBsdbScRsHnttI+IPVMX1WIeSK+bMCkAUgW3G8FF7f8pMIDHI+EU3ZQshVoKhzJ7x+fN2/Gbxv5julMjko3CfMZ7YEQC/p+CHpWcX0eYWQK+bPS+rvHQA3K86j+rzC/uIm6g+A/gDWAFgPYCOAB6Ty7gBWwBv8eg9AplTeTDovk653V7R1P7zxk20AxkdIrkUANkhlb8A/U4sAPCf1vwHAEEVbP5fkLVP+k4Up30j4FXVMn1UIuWL+rKRns0762QTgfqm8LYCFAEql322iKVsIuV6X+l0PYC4CFaTm3wzemT3bpWv3O/DMBgJYJcnwX3gVWkyfVwi54uF5tQBwFECOoiyqz4tTmTAMwzC24BgIwzAMYws2IAzDMIwt2IAwDMMwtmADwjAMw9iCDQjDMAxjCzYgDMMwjC3YgDAMwzC2+H9iKiaR51CGSQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "#read in stellar model, plot and determine dispersion\n", "from astropy.io import ascii\n", "a=ascii.read('star.txt')\n", "\n", "# copy the wavelength array into a variable for convenience\n", "wave=a['WAVE']\n", "\n", "#plot the model spectrum\n", "plt.plot(wave,a['FLUX'])\n", "\n", "# see what the dispersion is, checking to see if it is constant\n", "dispersion=wave[1:]-wave[0:-1]\n", "print('minimum dispersion: {:f} maximum dispersion: {:f}'.format(dispersion.min(),dispersion.max()))\n", "disp=dispersion.mean()\n", "print('using dispersion: {:f}'.format(disp))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As before, we'll simulate a low resolution spectrum with $\\delta\\lambda=2.5$ Angstrom. Create the smoothing kernel for this resolution, in the pixel units of our model spectrum." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dlambda = # delta(lambda) in Angtroms\n", "sig = # corresponding sigma in Angstroms\n", "sig_pixels = sig / disp # corresponding sigma in pixels, given disp\n", "\n", "def gauss(x,mean=0,sig=1) :\n", " \"\"\" Return normalized Gaussian function here\"\"\"\n", " return 1/sig/np.sqrt(2*np.pi)*np.exp(-(x-mean)**2/2/sig**2) # add code here\n", "\n", "# create a gaussian kernel with this sigma. Let this extend from -5*sig to 5*sig. Note that we want this to be an odd \n", "# number to prevent the convolution from shifting the spectrum.\n", "x=np.arange(-5*sig_pixels,5*sig_pixels+1)\n", "kernel= gauss(x,0,sig_pixels) #use your gauss routine to fill the kernel array\n", "\n", "# plot it\n", "plt.plot(kernel)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "OK, now convolve our spectrum with the smoothing kernel using numpy.convolve() " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sim=np.convolve(a['FLUX'],kernel,mode='same')\n", "plt.clf()\n", "plt.plot(wave,a['FLUX'])\n", "plt.plot(wave,sim)\n", "plt.xlim(5000,5050)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "OK, let's add noise to the image. Suppose we want to generate a spectrum with some specified median S/N. For a target S/N value, what do we need to scale the median flux level to, given Poisson statistics?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sntarg= # Choose a S/N here\n", "fluxtarg= # what does the flux need to be to reach this target S/N?\n", "\n", "# scale the model spectrum to the desired target flux\n", "sim=sim/np.median(sim)*fluxtarg\n", "# plot should have a median of the desired level\n", "plt.plot(sim)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now add Poisson noise using np.random.poisson() " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "noisesim = # get a Poisson sampling of the spectrum\n", "\n", "#plot both the spectrum with noise and the noiseless one\n", "plt.plot(wave,noisesim)\n", "plt.plot(wave,sim)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we will add in some readout noise, which is represented by a Gaussian distribution with zero mean and a $\\sigma$ corresponding to the readout noise level of the detector, let's say 10 units here." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "rn=10 # desired readout noise\n", "\n", "noisesim2= # add in Gaussian noise\n", "plt.plot(wave,noisesim2,label='spectrum+Poisson+readout')\n", "plt.plot(wave,noisesim,label='spectrum+Poisson')\n", "plt.plot(wave,sim,label='spectrum')\n", "plt.legend()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For later use, let's write a routine that scales an input spectrum and adds the noise" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def addnoise(spec,sn=100,rn=10) :\n", " \"\"\" Add noise to an input spectrum to get desired median S/N and input readout noise\n", " \"\"\"\n", " # calculate target flux given input S/N\n", " fluxtarg=sn**2 \n", " # scale the model spectrum to the desired target flux\n", " sim=spec/np.median(spec)*fluxtarg\n", " \n", " return np.random.poisson(sim)+np.random.normal(0,rn,len(sim))\n", "\n", "# plot results with this routine\n", "plt.plot(wave,addnoise(sim,sn=100,rn=10),label='spectrum+Poisson+readout')\n", "plt.plot(wave,addnoise(sim,sn=100,rn=0),label='spectrum+Poisson')\n", "plt.plot(wave,sim,label='spectrum')\n", "plt.legend()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's simulate a spectrum with a radial velocity of 200 km/s. To do this, you can just shift the wavelengths by the Doppler shift:\n", "$$\\lambda_{shifted} - \\lambda_{rest} = \\lambda_{rest} \\frac{v}{c}$$\n", "To do this, we will initialize an interpolator with these shifted wavelengths using scipy.interpolate() , which we can then use to get a simulated spectrum at any desired wavelengths. We are also going to want to interpolate the rest spectrum to resampled wavelengths as well, so we can cross correlate them, so initialize an interpolator for that as well." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "rv= # set desired radial velocity here (in km/s)\n", "\n", "#instantiate an interpolation for the rest spectrum\n", "interp=scipy.interpolate.CubicSpline(wave,sim)\n", "\n", "#instatiate a interpolator given the input (shifted) wavelengths and the spectrum\n", "rvinterp=scipy.interpolate.CubicSpline(wave+wave*rv/3.e5,sim)\n", "\n", "# plot the rest spectrum (here the interpolator is doing nothing!)\n", "plt.plot(wave,interp(wave),color='g',label='rest')\n", "\n", "# plot the shifted spectrum\n", "plt.plot(wave,rvinterp(wave),color='r',label='shifted')\n", "plt.legend()\n", "\n", "#plot a small segment of both\n", "plt.figure()\n", "plt.plot(wave,interp(wave),color='g',label='rest')\n", "plt.plot(wave,rvinterp(wave),color='r',label='shifted')\n", "plt.legend()\n", "plt.xlim(5000,5050)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "OK, to get the radial velocity by cross correlation, we need to resample the spectra to a $\\log\\lambda$ scale, so that the shift is independent of the wavelengths. Let's choose a $\\log\\lambda$ array with the same number of pixels and the same wavelength range as the original wavelengths." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import scipy\n", "\n", "# set up a log(wavelength) array with the same number of pixels as the original\n", "#numpy.linspace returns a uniformly spaced array between a starting and ending values with the \n", "# specified number of pixels\n", "logwave=np.linspace(np.log(wave[0]),np.log(wave[-1]),len(wave))\n", "\n", "#dlog_lambda is the spacing in log(lambda)\n", "dlog_lambda=logwave[1]-logwave[0]\n", "\n", "# interpolate to the new wavelength scale. Remember that we need to give the interpolator the wavelengths\n", "# not the log(wavelengths)\n", "sim_logwave = # use the interp interpolator with the new wavelengths\n", "\n", "# plot the old and new, demonstrate they are similar, i.e. interpolator works well here\n", "plt.plot(wave,sim)\n", "plt.plot(np.exp(logwave),sim_logwave)\n", "plt.xlim(5000,6000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we'll be a little unrealistic (cheat!) and interpolate the rv shifted spectrum first, then add noise" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# interpolate the RV-shifted spectrum to the log(lambda) grid\n", "rvsim_logwave = # use the rvinterp interpolator with the new wavelengths\n", "\n", "# need to be careful around edges that cubic spline interpolation doesn't give negative values,\n", "# because they will cause the Poisson routine to fail (why?). Set negatives to 0.\n", "\n", "rvsim_logwave[rvsim_logwave<0]= 0.\n", "\n", "# now add noise using the addnoise routine from above\n", "rvsim_noise= # add an appropriate call to addnoise()\n", "plt.clf()\n", "plt.plot(np.exp(logwave),sim_logwave,label='noiseless rest')\n", "plt.plot(np.exp(logwave),rvsim_noise,label='RV shifted with noise')\n", "plt.plot(np.exp(logwave),rvsim_logwave,label='RV shifted, no noise')\n", "plt.legend()\n", "plt.xlim(5000,5050)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "OK, now we will see if we can recover the radial velocity using cross correlation.\n", "

\n", "Apart from sampling on a $\\log\\lambda$ grid, another thing we need to be careful of in cross correlation is the continuum level. If the signal has a mean value greater than 0, the cross correlation will be \"diluted\" by the continuum contribution, and edge effects will be very significant; imagine that the signal drops to zero off the edges, which will lead to a big change in cross-correlation as the edges slide relative to each other. As a result, the mean level is typically removed before cross-correlation. Here, the continuum level is clearly a function of wavelength. We'll remove an approximate continuum by subtracting a \"median-filtered\" version of the spectrum: imagine running a window across the spectrum and taking the median value in the window around each pixel. We will use scipy.signal.medfilt() to do this, although the principle is simple.\n", "\n", "We will build this continuum subtraction into a short cross-correlation function, using numpy.correlate to do the actual cross-correlation of the continuum-removed spectra. numpy.correlate() will cross-correlate two spectra across a large range of shifts, depending on the mode keyword, but even the smallest of these is likely overkill for the RV problem, if we expect that the shift is small (pixels or tens of pixels). Here, we use it with the mode='same' option, which will produce a cross-correlation function of the same size as the first array, which means it is computing shifts from the second array being shifted one half of its length to the left up through one half of its length to the right: zero lag will be in the middle of cross-correlation array.\n", "\n", "It would be very straightforward, and possibly simpler to understand, to write your own (perhaps one-line!) cross correlation routine that returned the cross-correlation function as a function of shift (lag)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import scipy.signal\n", "def xcorr(a,b,filter=401) :\n", " \"\"\" Cross correlate two spectra, but remove the median levels first\n", " \n", " filter gives the width for the median filtering\n", " \n", " returns the lag array and the cross correlation function\n", " \"\"\"\n", " xc=np.correlate(a-scipy.signal.medfilt(a,filter),b-scipy.signal.medfilt(b,filter),mode='same')\n", " # with mode='same', the zero shift cross-correlation is in the middle of the output array\n", " if len(xc)%2 == 0 :\n", " lag=np.arange(len(xc))-len(xc)/2\n", " else :\n", " lag=np.arange(len(xc))-(len(xc)-1)/2\n", "\n", " return lag, xc \n", "\n", "# cross correlate template with itself as a check\n", "lag,xc=xcorr(sim_logwave,sim_logwave)\n", "\n", "# cross correlate template with object\n", "lag,xc2= # add the call to xcorr with the arrays you want to cross correlate\n", "\n", "# plot cross correlation as a function of lag\n", "plt.plot(lag,xc/xc.max())\n", "plt.plot(lag,xc/xc.max(),'bo', label='autocorrelation')\n", "plt.plot(lag,xc2/xc2.max(),'ro',label='cross correlation of shifted and rest')\n", "plt.legend()\n", "plt.xlim(-20,20)\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can look for the peak of the cross correlation, which will give us the pixel shift in the $\\log\\lambda$ data arrays. To turn that into a velocity, we need to determine the velocity shift per pixel, which is just $d\\log\\lambda \\times c$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print('dloglambda: {:f} velocity per pixel: {:f}'.format(dlog_lambda,dlog_lambda*3.e5))\n", "# get the peak pixel, the lag corresponding to this pixel, and the corresping RV\n", "imax=xc2.argmax()\n", "print('peak pixel: {:d} corresponding lag: {:f}'.format(imax,lag[imax]))\n", "print('corresponding RV: {:f}'.format(lag[imax]*dlog_lambda*3.e5))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How well does this match the input radial velocity? Note that we can do (much) better by fitting for the peak: if you look carefully you can see that the cross-correlation is not symmetric around the peak pixel. Write this as a function so we can use it later (always good to write things as function!)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot around peak\n", "plt.plot(lag,xc2/xc2.max(),'ro',label='cross correlation of shifted and rest')\n", "plt.xlim(lag[imax-5],lag[imax+5])\n", "plt.ylim(0.8,1.05)\n", "print(lag[imax])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can see by eye that the maximum may not be exactly on an integer pixel. Getting a little ahead of ourselves, let's do a quick quadratic fit to the peak." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def fitpeak(lag,xcorr,plot=True) :\n", " \"\"\" Do a quadratic fit to the peak of a cross correlation function\n", " \"\"\"\n", " # get the peak pixel\n", " imax=xcorr.argmax()\n", " \n", " # set up x and y arrays +/- 5 pixels around peak\n", " x=lag[imax-5:imax+6]\n", " y=xcorr[imax-5:imax+6]/xcorr.max()\n", " \n", " # do a 2nd order (quadratic)polynomial fit. polyfit() returns highest order coefficient first\n", " fit=np.polyfit(x,y,2)\n", " \n", " if plot :\n", " # plot the cross correlation peak and fit\n", " plt.plot(x,y,'ro',label='cross correlation of shifted and rest')\n", " plt.xlim(lag[imax-5],lag[imax+5])\n", " xx=np.linspace(x[0],x[-1],100)\n", " plt.plot(xx,fit[0]*xx**2+fit[1]*xx+fit[2])\n", "\n", " # determine the peak of the quadratic\n", " #for a quadratic fit y = ax**2 + bx + c, the peak is where the derivative, 2ax +b, equals 0, so xpeak=-b/2a\n", " return -fit[1]/2/fit[0]\n", "\n", "peak=fitpeak(lag,xc2)\n", "\n", "#convert this to a velocity as above\n", "print('peak lag: {:f} corresponding RV: {:f}'.format(peak,peak*dlog_lambda*3.e5))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Better? Try doing this at different S/N for the input spectrum, and report your results. See how compact it is if we use functions for the different steps!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for sn in [5,10,50,100,1000] :\n", " # get the simulated spectrum for this S/N\n", " out=addnoise(rvsim_logwave,sn=sn,rn=10)\n", " \n", " # plot the spectra\n", " plt.figure()\n", " plt.plot(out)\n", " \n", " # do the cross correlation including continuum removal\n", " lag,xc2=xcorr(out,sim_logwave)\n", " \n", " # fit the cross correlation peak and report results\n", " peak=fitpeak(lag,xc2,plot=False)\n", "\n", " #convert this to a velocity as above\n", " print('S/N: {:f} peak lag: {:f} corresponding RV: {:f}'.format(sn,peak,peak*dlog_lambda*3.e5))\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Comments or discussion?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.11" } }, "nbformat": 4, "nbformat_minor": 2 }