diff --git a/band_gap.ipynb b/band_gap.ipynb new file mode 100644 index 0000000..efcc360 --- /dev/null +++ b/band_gap.ipynb @@ -0,0 +1,341 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Silicon Band Gap Energy" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "# We will use the scipy curve_fit function to fit a model to data.\n", + "from scipy.optimize import curve_fit" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Enter data from table as numpy arrays" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "t_cel = np.array([-15., -10., -8., -6., -4.2, -2.3, 0., 2., \n", + " 3.6, 5.8, 8.2, 10, 12.8, 16.2, 20.])\n", + "adu = np.array([13, 15, 16, 17, 19, 22, 24, 28, 32, 37, \n", + " 43, 50, 62, 89, 139])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Convert Celsius to Kelvin" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[258.15 263.15 265.15 267.15 268.95 270.85 273.15 275.15 276.75 278.95\n", + " 281.35 283.15 285.95 289.35 293.15]\n" + ] + } + ], + "source": [ + "t_kel = t_cel + 273.15\n", + "print(t_kel)\n", + "t_inverse = 1 / t_kel" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Convert ADU (counts) to electrons" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "electrons = adu * 2.3" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Convert electrons to electrons per second" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "electrons_per_sec = electrons / 120" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot dark current (e-/sec) vs. inverse temperature\n", + "Make a plot in the cell below" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnMAAAHtCAYAAAByNPigAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3XmYlXX9//Hn+7AIA2quWcjMKGguGKa4kinuC8LB1H41fnPB0FwqkxbEBcVBy1bDVLDEarLSZFNcEtdQvyq5pKUBsphLKpoLm8B8fn/M6BdxYGbgzNznzDwf13Vfw9zLOS+uq0tefe7787kjpYQkSZJKUy7rAJIkSVp3ljlJkqQSZpmTJEkqYZY5SZKkEmaZkyRJKmGWOUmSpBJmmZMkSSphljlJkqQSZpmTJEkqYR2zDtCaNt9881RZWZl1DEmSpEbNnDnzjZTSFo2d167KXGVlJY8//njWMSRJkhoVEfObcp63WSVJkkqYZU6SJKmEWeYkSZJKmGVOkiSphFnmJEmSSphlTpIkqYRZ5iRJkkqYZU6SJKmEWeYkSZJKmGVOkiSphFnmJEmSSphlTpIkqYRZ5iRJkkqYZU6SJKmEWeYkSZJKmGWuQGpqaqisrCSXy1FZWUlNTU3WkSRJUjvQMesAbUFNTQ3Dhg1j8eLFAMyfP59hw4YBUFVVlWU0SZLUxjkyVwAjR478sMh9YPHixYwcOTKjRJIkqb2wzBXAggULmrVfkiSpUCxzBVBeXt6s/ZIkSYVimSuA6upqysrKPrKvrKyM6urqjBJJkqT2wjJXAFVVVYwbN46KiooP9/30pz918oMkSWpxlrkCqaqqYt68eTzyyCMAdOvWLeNEkiSpPbDMFdgee+zBpz71KSZNmpR1FEmS1A5Y5gosl8sxePBgbr/9dpYuXZp1HEmS1MZZ5lpAPp9n0aJFTJ8+PesokiSpjbPMtYABAwaw0UYbMXHixKyjSJKkNs4y1wI6d+7MUUcdxZQpU1i5cmXWcSRJUhtmmWsh+Xye119/nYcffjjrKJIkqQ2zzLWQww8/nM6dOzurVZIktSjLXAvZaKONOOigg5g0aRIppazjSJKkNsoy14Ly+Txz5szhmWeeyTqKJElqoyxzLWjQoEFEhLdaJUlSi7HMtaCtttqKffbZxzInSZJajGWuheXzef72t7+xYMGCrKNIkqQ2yDLXwvL5PACTJ0/OOIkkSWqLLHMtbLvttmOnnXbyVqskSWoRlrlWkM/nuf/++1m4cGHWUSRJUhtjmWsFQ4YMYeXKldx2221ZR5EkSW1Mq5a5iDg2Iv4cEfMjYklEPB8Rl0XEhk24Nq1h27U1sq+P3XffnR49enirVZIkFVzHVv6+4cAC4Dzg38DngFHAgIjYN6VU28j1E4BrV9v3rwJnLLiIIJ/P8+tf/5rFixdTVlaWdSRJktRGtPZt1qNTSsenlGpSSvenlH4GfAPYCzigCde/lFJ6ZLVtcYsmLpB8Ps+SJUu4++67s44iSZLakFYtcyml1xvY/Vj9zx6tmaW17b///my88cZMnDgx6yiSJKkNKYYJEPvX//xnE879ekQsi4jFEXFPROzXksEKqVOnTgwcOJCpU6eyYsWKrONIkqQ2ItMyFxE9gEuAu1NKjzdy+u+AM4CDgWHAZsA9EXFAi4YsoCFDhrBw4UJmzJiRdRRJktRGZFbmIqI7MBlYAZzc2Pkppf9JKf0xpfRgSul3wOeBl4FLG/meYRHxeEQ8/vrrDd3lbT2HHXYYG2ywgbNaJUlSwWRS5iKiCzAF2BY4LKX07+Z+RkrpXeA2YI9GzhuXUuqXUuq3xRZbrFPeQunevTuHHHIIkyZNIqWUaRZJktQ2tHqZi4hOwJ+BPYEjU0p/X5+PA0qqFeXzeebNm8fTTz+ddRRJktQGtPaiwTmgBjgIGJxSemQ9Pmsj4CjgfwsUr1UcffTR5HI5Z7VKkqSCaO2RuauA44AfAYsiYu9Vtq0BIqIiIlZExIUfXBQRwyNifER8JSIOiIgTgRnAVsD5rfx3WC9bbrkl/fv397k5SZJUEK1d5o6o/zkSeHi17dT6YwF0WC3b88BOwJXAX4CfAHOBz6eUHmz52IWVz+d56qmnmDt3btZRJElSiWvtRYMrU0qxhm1U/TnzVv29ft/UlFL/lNLmKaVOKaXNUkqDUkqPtmb+Qhk8eDAAkydPzjiJJEkqdcWwaHC706tXL3bZZRdvtUqSpPVmmctIPp/nwQcf5I033sg6iiRJKmGWuYzk83lqa2uZOnVq1lEkSVIJs8xl5HOf+xzl5eXeapUkSevFMpeRiCCfz3PXXXexaNGirONIkqQSZZnLUD6fZ+nSpdx1111ZR5EkSSXKMpeh/fbbj0022cRbrZIkaZ1Z5jLUsWNHjj76aKZOncry5cuzjiNJkkqQZS5jQ4YM4a233uLBB0vuRRaSJKkIWOYyduihh9K1a1dvtUqSpHVimctYWVkZhx56KJMmTSKllHUcSZJUYixzRSCfz/Piiy/yxBNPZB1FkiSVGMtcERg4cCC5XM5brZIkqdksc0Vg8803Z7/99mPixIlZR5EkSSXGMlckhgwZwjPPPMPs2bOzjiJJkkqIZa5IDB48GIDJkydnnESSJJUSy1yRqKysZNddd/W5OUmS1CyWuSKSz+eZMWMGr732WtZRJElSibDMFZF8Pk9KiSlTpmQdRZIklQjLXBH57Gc/S2VlpbdaJUlSk1nmikhEMGTIEO6++27efffdrONIkqQSYJkrMvl8nmXLlnHnnXdmHUWSJJUAy1yR2Xfffdl888291SpJkprEMldkOnbsyNFHH82tt97K8uXLs44jSZKKnGWuCOXzed5++23uu+++rKNIkqQiZ5krQocccghlZWXeapUkSY2yzBWhrl27cvjhhzN58mRqa2uzjiNJkoqYZa5I5fN5XnrpJWbOnJl1FEmSVMQsc0XqqKOOokOHDt5qlSRJa2WZK1Kbbrop+++/v2VOkiStlWWuiOXzef7xj3/wr3/9K+sokiSpSFnmilg+nwdwdE6SJK2RZa6I9ezZk913390yJ0mS1sgyV+Ty+TyPPPIIr7zyStZRJElSEbLMFbl8Pk9KialTp2YdRZIkFSHLXJHbeeed6dWrFxMnTsw6iiRJKkKWuSIXEQwZMoTp06fzzjvvZB1HkiQVGctcCcjn8yxfvpzbb7896yiSJKnIWOZKwN57782WW27prFZJkvQxlrkS0KFDBwYNGsRtt93GsmXLso4jSZKKiGWuROTzed59913uu+++rKNIkqQiYpkrEQcddBDdu3d3VqskSfoIy1yJ6NKlC0cccQSTJ0+mtrY26ziSJKlIWOZKSD6f59VXX+XRRx/NOookSSoSlrkScuSRR9KxY0dntUqSpA9Z5krIJz7xCQYMGGCZkyRJH7LMlZh8Ps/zzz/PP//5z6yjSJKkImCZKzGDBw8GcHROkiQBlrmS06NHD/bcc0/LnCRJAixzJSmfz/Poo4/y0ksvZR1FkiRlzDJXgvL5PABTpkzJOIkkScqaZa4E7bDDDmy//fbeapUkSZa5UhQR5PN57rnnHv773/9mHUeSJGXIMleihgwZwooVK5g2bVrWUSRJUoYscyVqzz33ZKuttvJWqyRJ7ZxlrkTlcjkGDx7M7bffztKlS7OOI0mSMmKZK2H5fJ733nuPe+65J+sokiQpI5a5EjZgwAA23HBDb7VKktSOWeZK2AYbbMCRRx7J5MmTWblyZdZxJElSBixzJW7IkCG89tprPPLII1lHkSRJGbDMlbgjjjiCTp06eatVkqR2yjJX4jbaaCMOOuggJk6cSEop6ziSJKmVWebagHw+z5w5c/jHP/6RdRRJktTKLHNtwKBBg4gIJk6cmHUUSZLUyixzbcCnPvUp9t57b5+bkySpHbLMtRH5fJ6ZM2fy4osvZh1FkiS1IstcG5HP5wGYPHlyxkkkSVJrssy1Edtvvz077rijt1olSWpnLHNtSD6f57777uOtt97KOookSWolrVrmIuLYiPhzRMyPiCUR8XxEXBYRGzbh2i4RcUVEvFJ/7cMR8YXWyF0q8vk8K1eu5NZbb806iiRJaiWtPTI3HFgJnAccDlwNfB34S0Q0luVXwNeAC4GBwCvAnRGxa8vFLS39+vWjR48e3mqVJKkd6djK33d0Sun1VX6/PyLeBG4ADgDuaeiiiOgLfAU4JaV0ff2++4FngUuAQS0ZulTkcjkGDx7MhAkTWLJkCV27ds06kiRJamGtOjK3WpH7wGP1P3us5dJBwHLgj6t81grgD8BhEbFBwUKWuHw+z+LFi7n77ruzjiJJklpBMUyA2L/+5z/Xcs7OwNyU0uLV9j8LdAZ6t0SwUrT//vuz8cYbe6tVkqR2ItMyFxE9qLtNendK6fG1nLop0NAUzTdXOb6m7xgWEY9HxOOvv97QwGDb0rlzZ4466iimTJnCihUrso4jSZJaWGZlLiK6A5OBFcDJjZ0OpDXsX6uU0riUUr+UUr8tttii+UFL0JAhQ3jjjTd46KGHso4iSZJaWCZlLiK6AFOAbYHDUkr/buSSN2l49G2TVY6r3mGHHcYGG2zgrVZJktqBVi9zEdEJ+DOwJ3BkSunvTbjsWWCbiChbbf9OwPvA7MKmLG0bbrghBx98MJMmTSKlhgY0JUlSW9HaiwbngBrgIGBwSumRJl46BegEHLfKZ3UEvgTclVJaVuispS6fzzN37lz+/vemdGVJklSqWntk7irqCtmPgEURsfcq29YAEVERESsi4sIPLkopPUndsiQ/i4hTI+Ig6pYl2Qa4qJX/DiXh6KOPJiK81SpJUhvX2mXuiPqfI4GHV9tOrT8WQIcGsp0MXA9cCtwG9AQOTyn9rYUzl6RPfvKT7LvvvkycODHrKJIkqQU1+Q0QEbE3da/g2hv4NNAVeAN4HrgfmJRSWusb3lNKlY19T0ppHg3MUk0pLQG+Xb+pCYYMGcLw4cOZN28elZWVWceRJEktoNGRuYg4MSL+DjwEfAsoA2YB/0vd2m97AdcBL0XEhIjYpgXzqhkGDx4MwOTJkzNOIkmSWspay1xEPAVcDkwDdgc2SSl9IaX0xZTSCSmlI1NKO1K3bMjXgC2BZyPiSy0dXI3r3bs3ffr08bk5SZLasMZG5q4HtkkpfS+l9ERawzoXKaW3U0o1KaUjgX2A/xY6qNZNPp/ngQceYOHChVlHkSRJLWCtZS6l9LOU0tLmfGBK6amU0p3rF0uFks/nqa2tZerUqVlHkSRJLaDJs1kjolNEdFvDsW71iwGryOy222707NnTW62SJLVRzVma5Dpg/BqOXVu/qchEBPl8nrvuuovFixdnHUeSJBVYc8rcAGBN0yKnUPdWBxWhfD7PkiVLuOuuu7KOIkmSCqw5ZW5L4LU1HHsd+OT6x1FL2G+//dhkk0281SpJUhvUnDL3GrDLGo7tAjhdskh16tSJgQMHMnXqVFasWJF1HEmSVEDNKXO3AhdExGdX3RkRu1D3ei6nSxaxfD7Pm2++yYMPPph1FEmSVEDNKXMXUrd+3MyIeCgi/hQRM4C/AW8D57dEQBXGYYcdRpcuXbzVKklSG9PkMpdSegPYA7iMunen7lr/sxrYo/64ilS3bt049NBDmTRpEmtY+1mSJJWg5ozMkVL6b0rpwpTSPiml7VNK+6aURqWU3m6pgCqcfD7PggULePLJJ7OOIkmSCqRZZQ4gIjaPiIERcWJEbFq/r0tENPuz1LoGDhxILpfzVqskSW1Ic94AERFxBfBv6taV+zVQWX94MnWTIFTEtthiCz7/+c9b5iRJakOaM5o2AjgLuATYi7rn5T4wFRhYwFxqIfl8nqeffpoXXngh6yiSJKkAmlPmTgUuSSmNoW4G66pmA70KlkotJp/PAzg6J0lSG9GcMtcDeGQNx94Huq1/HLW0bbbZhr59+1rmJElqI5pT5l4C+qzhWF9g7vrHUWvI5/PMmDGD115b09vZJElSqWhOmbsJuDAi+q+yL0XE9sC5wB8KmkwtJp/PU1tby6233pp1FEmStJ6aU+ZGAc8BDwCz6vfdBPy9/vfLC5pMLaZv375UVFQwceLErKNIkqT11Jw3QCwBDgBOAh4C7gYeA4YBh6SU3m+BfGoBEcGQIUP4y1/+wnvvvZd1HEmStB6a+waIlSml36aUTkgpHZpS+nJK6YaU0oqWCqiWkc/nWbZsGXfeeWfWUSRJ0npozqLBm0dE+Wr7TouIX0SEa8yVmP79+7PZZps5q1WSpBLXnJG5XwPf/+CXiLgAuBr4CjA5Ir5U4GxqQR07duToo4/m1ltvZfny5VnHkSRJ66g5Za4fMH2V308HxqSUNgOuAr5dyGBqefl8nv/+97888MADWUeRJEnrqDllblPgPwAR0QfYCrih/tgk4DOFjaaWdsghh9C1a1dntUqSVMKaU+YWAlvX//lA4OWU0gdLlHRq5mepCJSVlXH44YczadIkUkpZx5EkSeugOQXsbmBURJxF3SLBqz45vwMwv5DB1Dry+TwvvfQSM2fOzDqKJElaB80pc98FXgQuA+YAF69yrAr4awFzqZUcddRRdOjQwVmtkiSVqFjb7bWI6JRSanSqY0RsBCwt9oWD+/Xrlx5//PGsYxSdAw88kNdee41nnnkm6yiSJKleRMxMKfVr7LzGRuYWRsQfI+L/1Re2BqWU3in2Iqc1y+fzPPvss8yaNavxkyVJUlFprMydBiTgGuC1iLgzIr4eEZ9u+WhqLYMHDwbwVqskSSVorWUupXRjSun/AVsAQ4C5wPnAixHxvxExIiJ2bIWcakEVFRXstttuljlJkkpQkyZApJSWp5RuTymdnlLqAXweuBf4KvBMRDwfET+IiL1bMqxaTj6f5+GHH+bVV1/NOookSWqGdVobLqX0cErp+ymlHYGdgeuB/XFGa8nK5/OklJg6dWrWUSRJUjOs90K/KaXnUkqXp5T25v8WFVaJ6dOnD9tuu623WiVJKjHrVOYiokNErIyI3Vbdn1LyHl2Jigjy+Tx3330377zzTtZxJElSE63PyFwULIWKQrdu3Xj//ff5xCc+QWVlJTU1NVlHkiRJjVifMufLPNuQmpoafvzjHwOQUmL+/PkMGzbMQidJUpFzZE4AjBw5ksWLF39k3+LFixk5cmRGiSRJUlOs62zWlcAA4PnCxlFWFixY0Kz9kiSpOKzrBIjuwDzAV3i1EeXl5c3aL0mSikOzylxEDIyIvwFvA3OAXer3XxcRX2mBfGol1dXVlJWVfWRfhw4dqK6uziiRJElqiiaXuYjIA5OBN4DvrXbtXODEwkZTa6qqqmLcuHFUVFQQEWy44YasXLmS3XffPetokiRpLZozMncRcH1K6VDgZ6sdewboU7BUykRVVRXz5s2jtraWOXPm0L17dy644IKsY0mSpLVoTpnbEfhj/Z9XX5bkLWCzgiRSUdhiiy0YPnw4N998M4899ljWcSRJ0ho0p8y9A2y+hmOVwOvrnUZF5dvf/jabb745I0aMyDqKJElag+aUub8AIyLiE6vsSxGxAXAWcHtBkylzG264Ieeffz7Tp0/n7rvvzjqOJElqQKTUtBc5REQl8Ch1t1inAV8FbgY+C2wM9EspvdwiKQukX79+6fHHH886RklZtmwZn/nMZ9h888157LHHiHCtaEmSWkNEzEwp9WvsvCaPzKWU5gG7AbcChwArgS8AjwB7FXuR07rZYIMNuOSSS5g5cyZ//vOfs44jSZJW0+SRubbAkbl1s3LlSvr27cvy5ct59tln6dixY9aRJElq8wo+Mqf2q0OHDowZM4Z//etfXH/99VnHkSRJq1hrmYuIn0fEVs35wIg4JiL+3/rFUrE5+uij2XfffRk1ahRLlizJOo4kSarX2MjctsALEfHHiBgUER9bSy4ichGxa0RcEBHPA1cDb7ZEWGUnIrj88st5+eWXGTt2bNZxJElSvUafmYuILwDDgSOoK38vU7em3DJgE6An0AV4BRgP/DSl9E4LZl5nPjO3/o466igefvhhXnjhBT7xiU80foEkSVonBXtmLqX0QEppEFAOnAJMBRYAbwMzgcuBA4HylNLFxVrkVBhjxozhrbfe4oc//GHWUSRJEs5m1Tqoqqpi4sSJzJkzh0996lNZx5EkqU1yNqtazCWXXMLy5csZPXp01lEkSWr3LHNqtl69enHaaacxfvx4Zs+enXUcSZLaNcuc1sn5559P586dufDCC7OOIklSu2aZ0zrZaqutOOecc7jxxht54oknso4jSVK7ZZnTOvvOd77DpptuynnnnZd1FEmS2q0ml7mI+GpDiwbXH9s0Ir5auFgqBRtvvDEjRozgjjvu4L777ss6jiRJ7VJzRuauB3qt4dg29cfVzpx55plsvfXWjBgxgva0zI0kScWiOWUu1nKsG7BiPbOoBHXt2pVRo0bxyCOPMGXKlKzjSJLU7qx10eCI2BXYrf7X64BqYO5qp3UF/h+wSUqpT0uELBQXDW4ZK1asoE+fPnTo0IGnn36aDh06ZB1JkqSS19RFgzs2cnwwcFH9nxMwcg3nLQSGNj2e2pKOHTtSXV3Nsccey29/+1tOOumkrCNJktRuNDYytzHwCepusb4AHAOsvg7FMuA/qQQemHJkruWklNhrr734z3/+w/PPP0+XLl2yjiRJUkkryOu8Ukpvp5Tmp5TmUTfJYVr976turzanyEXE1hHxi4h4OCIWR0SKiMomXjuv/vzVt3xTv18tIyK4/PLLWbBgAddcc03WcSRJajfWOjK3xositgQ+NvSSUlrQhGsPAP4IzAQ6AIcC29QXxsaunQc8B4xa7dDzKaW3GrvekbmWd+ihh/LEE08wZ84cNtpoo6zjSJJUsgoyMrfaB24UEddHxGLgFeomQqy+NcUDKaVPppSOBG5q6vev4o2U0iOrbY0WObWOMWPG8MYbb/CTn/wk6yiSJLULjU2AWNVVwBeBXwF/p+5ZuWZLKdWuy3UqDf369eO4447jxz/+MWeccQZbbrll1pEkSWrTmrPO3GHAd1JKZ6eUxqWUblh9a6mQqzm6/lm7ZRHxiM/LFZ/Ro0ezZMkSqqurs44iSVKb19xFg59vqSBNNBU4m7piWQUsBSZGxAlruiAihkXE4xHx+Ouvv95KMdu3z3zmM5xyyilcffXVzJs3L+s4kiS1aU2eABERvwBqU0rfLNiXR5wKjKeJEyAauL4D8AiwVUqpZ2PnOwGi9bz00kv07t2b448/nhtuaK1BW0mS2o6CT4AA7gIGRsSvI+LYiDhw9W3d466blNJK6iZRbB0Rn2rt79ea9ejRg2984xv89re/5Zlnnsk6jiRJbVZzytxk6taaOwn4E3B3/faXVX5m4YN3xhb9osXtzfe+9z022mgjzjvvvKyjSJLUZjVnNuuAFkuxjiKiI3AcsCCl9GrWefRRm266Kd/73vc477zzmDFjBv379886kiRJbc46LRq83l8acWz9Hw8CTgfOAF4HXk8p3V9/zgrghpTS0Prfv0zdu2KnAS8CnwTOBD4PfDml9IfGvtdn5lrfokWL6N27N7179+aBBx4gIhq/SJIkNfmZueaMzH3wwZsDewObAVNTSm9GRBfg/WasIbf6YsG/rP95P3BA/Z871G8fmAtsCVwBbAosBh4DDk8p3dncv4daR7du3bjooov4+te/zu23386RRx6ZdSRJktqU5sxmDeCH1C0N0pm6Z9T2SCn9LSLuBP6aUhrdYkkLwJG5bCxfvpwdd9yRbt268cQTT5DLNedRTUmS2qeWmM06AjgLuATYi/+beAB1678NbFZCtRudOnXi0ksv5emnn+bGG2/MOo4kSW1Kc8rcqcAlKaUxwN9WOzYb6FWwVGpzjj/+eHbddVcuuOAC3n///azjSJLUZjSnzPWgboHehrwPdFv/OGqrcrkcl112GXPnzmX8+PFZx5Ekqc1oTpl7CeizhmN9qZugIK3RYYcdxgEHHMDo0aN57733so4jSVKb0JwydxNwYUSsulhYiojtgXOBRpcGUfsWEVx22WX85z//4ec//3nWcSRJahOaU+ZGAc8BDwCz6vfdBPy9/vfLC5pMbdLee+9NPp/nhz/8IQsXLsw6jiRJJa/JZS6ltIS6NeBOAh6i7hVejwHDgENSSj7Vria59NJLee+997jsssuyjiJJUslr0jpzEdEJOBJ4OqVUss/Guc5c8Tj55JO58cYbmTVrFj179sw6jiRJRaeg68yllJYDfwIq1zOXBMDFF19MSomLL7446yiSJJW05jwz9wJ1r9OS1lt5eTlnnnkm119/Pc8991zWcSRJKlnNKXM/BEZGxBYtFUbty4gRI+jWrRsjR47MOookSSWrYzPOPZC6F9zPjYhHgFeoez/rB1JK6cRChlPbtsUWWzB8+HAuuugiHn30Ufbcc8+sI0mSVHKaNAECICLm8dHytrqUUtq2EKFaihMgis+7775Lr1696NOnD9OnTyciGr9IkqR2oKATIABSSpUppW3WshV1kVNx2nDDDTn//PO59957ufvuu7OOI0lSyWlSmYuIzhHx04jYo6UDqf057bTTqKysZMSIEdTW1mYdR5KkktLUpUneB04DurZsHLVHG2ywAZdccgkzZ87k5ptvzjqOJEklpTmzWZ8AdmmpIGrfvvKVr9CnTx/OP/98li9fnnUcSZJKRnPK3LnA8IgYGD6lrgLr0KEDY8aMYdasWVx//fVZx5EkqWQ0Zzbri8DGQDdgBfAaH1+apKLgCQvI2azFLaXEfvvtx9y5c5k1axZlZWVZR5IkKTNNnc3anHXmprP2pUmk9RIRXH755ey3336MHTuW7373u1lHkiSp6DV5ZK4tcGSuNAwcOJAZM2bwwgsvsMkmm2QdR5KkTBR8nTmptYwZM4a3336bH/7wh1lHkSSp6DXnmbmvNnZOSuk3652oBTkyVzpOOOEEbrnlFmbPns2nP/3prONIktTqmjoy15wyt6bVXD/8gJRSh6bFy4ZlrnS88MIL7LDDDgwdOpSrr7466ziSJLW6lrjNuk0DWz/gYmAWsNc65JQatO2223Laaacxfvx4Zs2alXUcSZKKVnPezTq/ge1vKaVLgBtUEm/UAAAgAElEQVSBb7dcTLVH559/Pl26dOGCCy7IOookSUWrUBMgHgSOKtBnSQB88pOf5JxzzuGPf/wjf/vb37KOI0lSUSpUmdsbeK9AnyV9aPjw4Wy66aacd955WUeRJKkoNXnR4Ii4sIHdnYE+1I3KjS1UKOkDG2+8Meeddx7Dhw/n3nvvZcCAAVlHkiSpqKzvbNZlwHzgD8BlKaVlBcxWcM5mLU1Lly5lu+22o0ePHjz88MP4amBJUntQ8NmsKaVcA1vXlNIOKaVRxV7kVLq6dOnCxRdfzP/+7/8yefLkrONIklRUfJ2XSsKKFSvYZZddiAiefvppOnZszmuFJUkqPQUZmYuIXhExMyIGreWcQfXnVDY/ptQ0HTt2pLq6mn/+85/89re/zTqOJElFo7HbrN8CalNKU9Z0Qv2xFcA3ChlMWt2QIUPYc889ueiii1i6dGnWcSRJKgqNlblDgF834XN+DRyx/nGkNYsILr/8cl588UVf8SVJUr3GylwF8I8mfM5zQOV6p5EaMWDAAA499FCqq6t55513so4jSVLmGitzy6lbS64xnam71Sq1uDFjxrBw4UJ+9KMfZR1FkqTMNVbmZgH9m/A5nwf+tf5xpMbtvvvuHH/88fzkJz/hP//5T9ZxJEnKVGNl7mbg7IjYdk0nREQv4CzgpkIGk9Zm9OjRLF26lOrq6qyjSJKUqcbK3M+BV4FHI+KciOgdEZ3qt94RcQ7wCPAycGVLh5U+sP322zN06FCuueYa5s6dm3UcSZIys9Yyl1JaDBwEPAn8GHgeWFq/PV+/70ngkPpzpVZz0UUXkVJil112IZfLUVlZSU1NTdaxJElqVY0uo59SehU4OCL2AA4GetYfehG4O6X0WAvmk9bo3nvvBWDRokUAzJ8/n2HDhgFQVVWVWS5JklqTr/NSyaqsrGT+/Pkf219RUcG8efNaP5AkSQVUkNd5ScVswYIFzdovSVJbZJlTySovL29wf8+ePRvcL0lSW2SZU8mqrq6mrKzsY/t33HFH2tPjA5Kk9s0yp5JVVVXFuHHjqKioICIoLy/n8MMP58477+Rb3/qWhU6S1C40Opu1qSKie0rpvUJ9ntQUVVVVH5m5mlLi3HPP5ac//Sm1tbVceeWVRESGCSVJallNLnMRcWVK6RtrONYduJOmvfpLajERwY9//GNyuRw//vGPqa2tZezYsRY6SVKb1ZyRuZMj4tWU0phVd0ZEN+AO/m/9OSlTEcEVV1xBLpfjiiuuIKXE2LFjyeV8qkCS1PY0p8wdB0yOiFdSStcDREQZdUVuG+ALLZBPWicRwQ9+8ANyuRw/+MEPqK2t5Ze//KWFTpLU5jS5zKWU7oiIrwHjI+J1YDpwO9AL2D+lNKeFMkrrJCK47LLLyOVyXHbZZdTW1nLNNddY6CRJbUqzJkCklH4TEVsBfwKeAcqBA1JKs1oinLS+IoLq6mpyuRzV1dXU1tYybtw4C50kqc1Ya5mLiIb+xfsRsDXw/4CDgH99cF5KqbbgCaX1FBGMHj2aXC7H6NGjqa2t5brrrrPQSZLahMZG5lYAa1qsK4AnV/k9NeHzpExEBJdccgm5XI6LL76YlBLXXXcdHTp0yDqaJEnrpbHydQlrLnNSyRk1ahQRwahRo6itreXXv/61hU6SVNLWWuZSSqNaKYfUai666CIigosuuoja2lomTJhgoZMklawm3RaNiM7Aq8BJKaUpLRtJankXXnghuVyOCy64gNraWm644QY6dvQpAUlS6WnSv14ppfcjYgWwtIXzSK3m/PPPJ5fLMXLkSFJK/OY3v7HQSZJKTnP+5ZoEHAvc1UJZpFZ33nnnkcvlGDFiBCklfvvb31roJEklpTn/at0OXBkRN1NX7F5htckRKaV7CphNahXf//73yeVyfO9736O2tpaamhoLnSSpZDTnX6w/1/88pn77QKJumZIE+BS5StJ3v/tdcrkc3/nOd0gpUVNTQ6dOnbKOJUlSo5pT5ga0WAqpCAwfPpxcLse5555LbW0tN954o4VOklT0mvNu1vtbMohUDL797W+Ty+U455xz+NKXvsQf/vAHOnfunHUsSZLWyPcZSav51re+xc9+9jMmTpzIl770Jd5///2sI0mStEbNeso7IvoAQ4HPAF1WO5xSSgcVKpiUpW9+85vkcjm+8Y1vcNxxx3HTTTc5QidJKkpNLnMRsRdwPzAP2A54GtgEKAf+DcxugXxSZs4++2xyuRxnnXUWxx57LDfddBMbbLBB1rEkSfqI5txmHQPcAuxM3ezVoSmlSuBg6maxXlrwdFLGzjzzTK666iqmTp3KF7/4RZYtW5Z1JEmSPqI5Ze6zwO/4v7XlOsCHa8tdClxW2GhScTjjjDO4+uqrue222zjmmGNYutQXoUiSikdzylwnYFFKqRZ4E/jUKseeB/o05UMiYuuI+EVEPBwRiyMiRURlE6/NRcSIiJgXEUsj4qmI+GIz/g7SOjn99NO59tprmTZtGkOGDLHQSZKKRnPK3BygR/2fnwZOqS9XOeBk4NUmfk5v4HjgLeDBZnw/wGhgFDAWOAJ4BLgpIo5s5udIzTZs2DDGjx/PHXfcQT6ft9BJkopCc2azTgUOAH5P3fNztwHvACuB7sA3mvg5D6SUPgkQEacChzbloojYEhgOXJ5S+lH97nsjojdwOTCtid8vrbNTTz2VXC7HqaeeyuDBg5k0aRJdu3bNOpYkqR1rzqLBo1b5890RsTdwLNAVuCOldFcTP6e2uSHrHQZ0pu65vVX9Dvh1RGyTUpq7jp8tNdkpp5xCRDB06FAGDRrE5MmTKSsryzqWJKmdWue3iaeUngCeKGCWxuwMLOPjS6A8W/9zJ8Ayp1Zx8sknk8vlOPnkkxk0aBBTpkyx0EmSMtGcdea6AP2om/iQgFeAmSml1npwaFPgvymltNr+N1c5/jERMQwYBlBeXt5y6dTunHjiiUQEJ510EgMHDmTq1Kl069Yt61iSpHam0QkQEbFBRPycutJ0P/BH4E/AA8DCiPhRRLTG0vjB/y2Lsvr+NUopjUsp9Usp9dtiiy1aJpnara9+9av85je/4f7772fgwIEsWrQo60iSpHZmrSNzERHArcCBwGTqJhksoK5A9QQGAudQd4uzpWeUvglsEhGx2ujcJqscl1rdCSecQETw1a9+lSOPPJLbbruN7t27Zx1LktRONHab9VhgAHBsSmliA8evi4hjgD9FxDEppVsKnvD/PAtsAPTio8/N7VT/8x8t+N3SWlVVVZHL5TjhhBM48sgjmTZtmoVOktQqGrvN+mXgT2socgDUF7ibgKpCBmvAHcD7DXzPCcAzzmRV1r785S/z+9//noceeojdd9+d8vJycrkclZWV1NTUZB1PktRGNTYy9zng/CZ8zq00492sEXFs/R93r/95RES8DryeUrq//pwVwA0ppaEAKaXXIuKnwIiIeBf4G/Al6m4BD27qd0st6Utf+hJ//etfGTt27If75s+fz7Bhw4C6ETxJkgqpsTK3BXXPyDVmAbBlM773ptV+/2X9z/upW5gY6t792mG180YC7wHfBLai7jVix6eUpjbju6UWNXXqx//nuHjxYkaOHGmZkyQVXGNlroy6td0a8z7QpalfmlJa6wzUNZ2TUlpJ3Qhgk0cBpda2YEHD//9nTfslSVofTVlnrkdEbNvIOVsXIozUFpSXlzN//vyP7d9qq60ySCNJausaXWcOuBmY1ci2+m1Tqd2qrq7+2NsgIoI333yTO++8M6NUkqS2qrGRuZNbJYXUhnzwXNzIkSNZsGAB5eXlDB8+nPHjxzNw4EDGjx/PSSedlG1ISVKbER9/O1bb1a9fv/T4449nHUPt1DvvvMMxxxzD9OnTGT16NCNHjqRuXW5Jkj4uImamlPo1dl5TbrNKKoCNNtqIadOmccIJJ3DBBRdw+umns2LFiqxjSZJKXFMmQEgqkM6dO/Ob3/yGrbfemssvv5yXX36ZP/zhD3Tr1i3raJKkEuXInNTKIoLLLruMq666imnTpnHggQfy+uuvZx1LklSiLHNSRs444wxuueUWnn76afbdd1/mzJmTdSRJUgmyzEkZGjx4MPfccw9vvfUW++yzD48++mjWkSRJJcYyJ2Vsn332YcaMGXTv3p0BAwZw6623Zh1JklRCLHNSEfjMZz7DQw89xI477sjgwYMZP3581pEkSSXCMicVia222or77ruPQw89lGHDhnHhhRfSntaBlCStG8ucVES6d+/OlClTOPnkkxk9ejRDhw5l+fLlWceSJBUx15mTikynTp341a9+RXl5ORdffDGvvPIKN910E927d886miSpCDkyJxWhiGDUqFGMHz+ev/zlL+y///68+uqrWceSJBUhy5xUxE499VSmTJnCc889xz777MPzzz+fdSRJUpGxzElF7sgjj+S+++5j0aJF7Lvvvjz00ENZR5IkFRHLnFQC9thjDx5++GE222wzDjroICZNmpR1JElSkbDMSSWiV69ezJgxg759+/LFL36Rk046icrKSnK5HJWVldTU1GQdUZKUgWhP61j169cvPf7441nHkNbL4sWL+fznP88TTzzxkf1lZWWMGzeOqqqqjJJJkgopImamlPo1dp4jc1KJKSsrY+HChR/bv3jxYkaOHJlBIklSlixzUgl68cUXG9y/YMGCVk4iScqaZU4qQeXl5Q3ujwjGjh3rWyMkqR2xzEklqLq6mrKyso/s69KlCzvssANnn302u+yyC1OnTvXdrpLUDljmpBJUVVXFuHHjqKioICKoqKjguuuu45lnnmHKlCkADBo0iIMPPpgnn3wy47SSpJbkbFapDVq+fDnXXnsto0aN4s033+Skk07i0ksv5dOf/nTW0SRJTeRsVqkd69SpE2eddRazZ8/m3HPPpaamhu22246LL76YRYsWZR1PklRAljmpDfvEJz7BFVdcwT//+U+OOuooRo0axfbbb8+ECROora3NOp4kqQAsc1I7sO222/KnP/2Jv/71r2y99dacfPLJ9OvXj3vvvTfraJKk9WSZk9qR/v378/DDD/P73/+ehQsXcuCBB5LP5/nXv/6VdTRJ0jqyzEntTC6X48tf/jLPPfccY8aMYfr06ey8885885vfbPDNEpKk4maZk9qprl27MmLECGbPns3QoUMZO3YsvXv35qc//Snvv/9+1vEkSU1kmZPauU9+8pNcc801PPXUU+y11158+9vfZqedduKWW25x0WFJKgGWOUkA9OnThzvuuIPbb7+dLl268MUvfpH999+fxx57LOtokqS1sMxJ+ojDDz+cJ598kmuuuYbnn3+ePffck//5n//hxRdfzDqaJKkBljlJH9OxY0dOO+00Zs2axYgRI7jpppvYfvvtOf/883n33XezjidJWoVlTtIabbTRRowZM4bnn3+eY445hurqarbbbjvGjx/PypUrs44nScIyJ6kJKioqqKmp4ZFHHqFXr14MGzaMz33uc/zlL3+hpqaGyspKcrkclZWV1NTUZB1XktqVaE+z1fr165cef/zxrGNIJS2lxJ///Ge++93vMnfuXHK53EdeDVZWVsa4ceOoqqrKMKUklb6ImJlS6tfoeZY5Seti2bJlbLXVVvz3v//92LGKigrmzZvX+qEkqQ1papnzNqukdbLBBhvw9ttvN3hswYIFrZxGktovy5ykdVZeXt7g/pQSp5xyCi+88EIrJ5Kk9scyJ2mdVVdXU1ZW9pF9Xbt25bDDDuP3v/8922+/PUOHDmXu3LkZJZSkts8yJ2mdVVVVMW7cOCoqKogIKioqGD9+PHfccQcvvPACZ555JjU1NWy//fZ87Wtf8zk6SWoBToCQ1KJeeuklfvCDHzBu3DhWrlzJySefzHnnnUdlZWXW0SSpqDkBQlJR6NGjB1deeSVz5szh9NNP54YbbmC77bbjtNNOY/78+VnHk6SSZ5mT1Cp69OjBL37xC+bMmcNpp53GhAkT2G677Tj99NOd/SpJ68EyJ6lVbb311owdO5bZs2fzta99jeuvv57evXvz9a9/3VInSevAMicpEz179uSqq65i9uzZnHrqqfzqV7+id+/enHHGGbz44otZx5OkkmGZk5Spnj178stf/pLZs2czdOhQrrvuOnr37s2ZZ57Jv//976zjSVLRs8xJKgrl5eVcffXVzJ49m5NPPpnx48fTq1cvzjrrLEudJK2FZU5SUSkvL+eaa65h1qxZnHTSSVx77bX06tWLs88+m5deeinreJJUdCxzkopSRUUF1157LbNmzeLEE0/kmmuuoVevXnzjG9/g5ZdfBqCmpobKykpyuRyVlZXU1NRknFqSWp+LBksqCfPmzaO6upoJEybQoUMHDjjgAB544AGWLFny4TllZWWMGzeOqqqqDJNKUmE0ddFgy5ykkjJ37lyqq6v51a9+1eDxiooKXxsmqU2wzDXAMie1Hblcjob++xUR1NbWZpBIkgrL13lJatPKy8sb3L/JJpuwfPnyVk4jSdmxzEkqSdXV1ZSVlX1kXy6X480336RPnz78+c9/bnDkTpLaGsucpJJUVVXFuHHjqKioICKoqKjgN7/5DVOnTqVjx44ce+yx7Lvvvjz44INZR5WkFuUzc5LanBUrVnDDDTdw4YUX8vLLLzNo0CAuv/xydtxxx6yjSVKT+cycpHarY8eODB06lFmzZjFmzBjuu+8++vTpw7Bhwz5co06S2grLnKQ2q6ysjBEjRjBnzhzOPvtsJkyYwHbbbccFF1zAO++8k3U8SSoIy5ykNm/zzTfnZz/7Gc899xyDBg3i0ksvpXfv3owdO5b3338/63iStF4sc5LajW233ZYbb7yRxx57jD59+nD22Wez0047cdNNNznzVVLJssxJanf69evH9OnTmTZtGl27duX4449n77335v777886miQ1m2VOUrsUERxxxBE8+eSTXH/99bz88ssccMABHH300Tz77LPU1NRQWVlJLpejsrKSmpqarCNLUoNcmkSSgCVLlnDllVdy2WWX8fbbb9OhQwdWrlz54fGysjLGjRtHVVVVhikltSe+m7UBljlJjVm4cCHbbLMN77777seOVVRUMG/evNYPJaldcp05SVoHm222Ge+9916DxxYsWOBECUlFp9XLXET0jIibI+LtiHgnIm6JiIbfmP3xa9Matl1bOrek9qO8vOH/JKWU2Hnnnbn66qvXWPgkqbW1apmLiDLgHmAH4ETgf4DtgHsjolsTP2YCsM9q278KHlZSu1VdXU1ZWdlH9pWVlXH66adTVlbGGWecwdZbb825557LCy+8kFFKSarT2iNzXwO2BfIppUkppcnAIKACOK2Jn/FSSumR1bbFLRVYUvtTVVXFuHHjqKioICKoqKhg3LhxXH311Tz22GM89NBDHHHEEVx55ZX07t2bwYMHM336dG/BSspEq06AiIjpQJeUUv/V9t8PkFLav5HrE1CdUjp/Xb7fCRCSCumll17immuu4dprr+X1119n55135uyzz+aEE06gW7em3myQpIYV6wSInYFnGtj/LLBTEz/j6xGxLCIWR8Q9EbFf4eJJUtP16NGD0aNHs2DBAiZMmEDnzp05/fTT6dmzJ9/97ned+SqpVbR2mdsUeKuB/W8CmzTh+t8BZwAHA8OAzYB7IuKANV0QEcMi4vGIePz1119vfmJJakSXLl048cQTmTlzJn/96185+OCD+clPfkKvXr045phjuO+++7wFK6nFZLE0SUP/RYsmXZjS/6SU/phSejCl9Dvg88DLwKVruWZcSqlfSqnfFltssW6JJakJIoL+/fvzpz/9iblz5/K9732PBx54gAEDBtC3b1+uu+46Fi/2EV9JhdXaZe4t6kbnVrcJDY/YrVVK6V3gNmCP9cwlSQXVs2dPxowZw4svvsivfvUrcrkcX/va1+jZsyff//73WbBgQdYRJbURrV3mnqXuubnV7QT8Yx0/M2h4tE+SMte1a1dOOeUUnnjiCe6//34GDBjAFVdcwbbbbstxxx3HAw884C1YSeultcvcFGDviNj2gx0RUQn0rz/WLBGxEXAU8L8FyidJLSIi+MIXvsDNN9/MCy+8wPDhw7nnnnvYf//92W233bj++uuZMGEClZWV5HI5KisrqampyTq2pBLQ2kuTdAOeApYA51M3ojYa2BD4bErpvfrzKoA5wCUppUvq9w0HPgPcS91zchXAB/sOSik92Nj3uzSJpGKyePFifv/73/Pzn/+cZ575+ET/srIyxo0bR1VVVQbpJGWtKJcmSSktAg6k7o0NvwVqgLnAgR8UuXoBdFgt3/PU3Y69EvgL8JP6az/flCInScWmrKyMU089laeffpott9zyY8cXL17M6aefztixY3nwwQd5++23M0gpqdi16shc1hyZk1Sscrlck56dq6yspG/fvuy666707duXvn37fnhrVlLb0tSRuY6tEUaStHbl5eXMnz+/wf0zZszgqaee+sg2ZcqUD8vfhhtuyGc/+9kPy13fvn3ZZZddPvZ+WUltkyNzklQEampqGDZs2EfWoVvbM3OLFy/mmWee4amnnuLJJ5/kqaee4umnn+bdd98F6iZcbLfddh8bxevRowcRTVraU1LGmjoyZ5mTpCJRU1PDyJEjWbBgAeXl5VRXVzdr8kNtbS3z5s372Cje3LlzPzxn0003/cgIXt++fdlpp53YYIMNCpZDUmFY5hpgmZPUHr399tv8/e9//3AE76mnnuKZZ55hyZIlAHTs2JEddtiBXXfdlZUrV3LLLbewbNmyD693Vq2UDctcAyxzklRn5cqVzJo16yMjeE8++SQvv/xyg+ev6Zk+SS3HMtcAy5wkrd3aZtWec845HHvssey9997OnpVaQVGuMydJKm7l5eUN7u/atStXXXUV/fv3p7y8nG9+85s8+OCDrFy5spUTSlqdZU6S9KHq6uqPLWlSVlbG+PHjee211/jd737HHnvswbXXXssXvvAFtt56a8466yzuu+8+i52UEW+zSpI+oimzWd99912mTZvGTTfdxLRp01iyZAlbbrklQ4YM4bjjjmP//fenY0eXMpXWh8/MNcAyJ0mFt2jRIm6//XZuvvlmbr31VhYtWsRmm23GkCFDOPbYYznwwAPp1KlT1jGlkmOZa4BlTpJa1uLFi7nzzju5+eabmTp1Ku+++y6bbLIJ+XyeY489loMPPpjOnTtnHVMqCU6AkCS1urKyMoYMGUJNTQ2vvfYaU6ZMYeDAgdxyyy0cddRRbLnllnz1q19l6tSpLF269MPrampqPnzHbGVlJTU1NRn+LaTS4sicJKnFLVu2jOnTp3PzzTczadIk3nrrLTbccEOOPvpottxyS6699toPFzEGFyqWwNusDbLMSVL2li9fzj333MPNN9/MxIkTWbhwYYPnVVRUMG/evNYNJxURy1wDLHOSVFxWrFhB586d17hQ8ZgxY+jfvz/9+vX72JIpUlvX1DLnvHFJUmY6duy4xleFdezYkfPOO+/DP++2227su+++7LvvvvTv359Pf/rTrR1XKkpOgJAkZWpNCxVPmDCBN954g6lTp/Kd73yHrl27cu2113L88cfTo0cPKisrqaqq4qqrruLJJ5900WK1W95mlSRlrikLFQO8//77PPXUU8yYMePD7ZVXXgGge/fu7L333h+O3O21115svPHG6/xdUtZ8Zq4BljlJaltSSixY8P/bu/vguKr7jOPfx6wd2wJhbCIje2qZlzaYxG1IQoe05SUtxTSZEBpaCDGJSwOekkwamkIgEAhM4hhKyYS0U6ceOk14MWmABGgZEgYoQ7AdF/Bgg3knYF6sSLblFzmyZcs6/ePcFev1Xe36dX12n8/MnZXunnPv0f3tan97zzn3vsnChQtZtGgRCxcuZPny5QwODiKJ6dOn79A1u3jxYmbPnk1fX9/QNjxz1g5UTuZyOJkzM2t8vb29LFmyZCi5W7x4Mb29vQCMGDGCwcHBnep45qwdiJzM5XAyZ2bWfLZv386KFStYtGgRF198ccVy06ZNY9KkSblLe3s77e3tjB49ej+23Jqdk7kcTubMzJrb1KlTc2fOHnLIIcyYMYNVq1YNLVu3bt2p3Pjx4ysmfMXliCOOqHgvWo/Xs13hS5OYmZmVmTNnTu6YuXnz5u2QVIUQ6Onp2SG5K1+ef/55Ojs7c2fRtrW17XBWb9KkSbz99tssWLCA/v5+AFauXMns2bMBnNDZHvGZOTMzayp78+zY4OAgq1ev3iHJ6+zs3Cnx6+rqyh2rBzBy5EhOP/30HRK/0seJEydSKPjcSzNyN2sOJ3NmZlYP1e50cfzxx7Nq1Sq6u7t3KiOJiRMn5iZ6pY+1JH3u5k2Lu1nNzMwOEMPd6aKjo4OlS5cCMenr6uoaOrtXepav+PPTTz9NV1dXbtJX7N7NS/iWL1/O3Llz2bx5M+Bu3kbiZM7MzGw/qDReb86cOUO/FwoFJk+ezOTJk4fd1sDAAN3d3TskeeWJ39KlS+nu7q7YvQvQ19fHRRddxKOPPsqECRM4/PDDc5dx48YxYoRvGnWgcjJnZma2HxTPfu2Nbs5CoTA0wWI4xaSvs7OTE044Ibebd/PmzTz00EOsXr16aHJGuREjRgwle8MlfaVLa2srknb5b7Nd5zFzZmZmTaDSZVmKF0wOIdDX18eaNWtyl7Vr1+au37ZtW+7+CoVC1cSv/PmDDz7YCWAJj5kzMzOzIdW6eSXR0tJCS0sLHR0dNW0zhEBvb29Nid8LL7ww9Fze5VwARo0aVVPSV7qMHTt2zw9O4pzMmZmZNYG92c1bJInW1lZaW1s56qijaqozODjIhg0bhj3bV1yWLVvGmjVr6OnpqTgTeMyYMbvU/TthwoTdvpPHgTob2N2sZmZmdkDbvn0769evHzbxKz8ruG7duorba2lpqSnpK/35rrvuyj2zOX/+/H2W0Pk6czmczPv+3TkAAAwCSURBVJmZmTWHgYEBenp6akr8isvGjRsrbk9S7tnB4pjDfcFj5szMzKxpFQoF2traaGtrq7nO1q1bc7t+165dy9VXX51b580339xbTd5tTubMzMzMiBMw2tvbaW9v3+m5W265JXc28JQpU/ZH04blKwCamZmZVTFnzpydZs6WX/S5XpzMmZmZmVUxc+ZM5s+fT0dHB5Lo6OjYp5MfdoUnQJiZmZkdgGqdAOEzc2ZmZmYJczJnZmZmljAnc2ZmZmYJczJnZmZmljAnc2ZmZmYJczJnZmZmljAnc2ZmZmYJczJnZmZmljAnc2ZmZmYJczJnZmZmljAnc2ZmZmYJczJnZmZmljAnc2ZmZmYJczJnZmZmljAnc2ZmZmYJUwih3m3YbyStBlbuYrXDgTX7oDm27zhmaXLc0uS4pccxS0dHCOG91Qo1VTK3OyQ9FUL4SL3bYbVzzNLkuKXJcUuPY9Z43M1qZmZmljAnc2ZmZmYJczJX3fx6N8B2mWOWJsctTY5behyzBuMxc2ZmZmYJ85k5MzMzs4Q5mTMzMzNLWLLJnKTfkXS3pA2SNkr6qaQpNdYdLelGSZ2SNktaLOnknHIjJH1d0huStkhaJunsnHL/KemFrB2bsnJflnTQMG34I0mDkoKkwq799elKMW6SHsviVL5csvtHIh0pxiwre5ik70l6U1K/pLcl/XC3DkKCUoubpFMrvM+Ky4l7dkTSkFrcsnJjJV0n6eVsv29JulXS1N09DraLQgjJLcBY4BXgOeAs4FPAs8BrQEsN9e8A1gMXAX8G/BTYDHywrNwcoB+4FPgY8O/AIPDxsnI/Br4IzABOB27Kyt1cYf8js/Z2AgEo1PuYOm6V4wY8BiwDTixbjqj3MXXMKsbsMGBF1u7PAycDnwH+pd7H1HHLjxvQmvMeOzGLYydwUL2Pq+NW8f22AOgDLsu2Nwt4I2v3wfU+rs2w1L0Bu9Vo+AqwHTimZN2RwADw1Sp1/4CYQF1Qsq4AvATcX7KuLXuxX1dW/xFgeQ1tvBPorfDcldmbdQ7NlcwlGTdiMvdEvY+fY7ZLMfsB8W4vrfU+ho5b7XHLKdNBTB5urPcxddzy4waMydr3nbJyZ2TtmVHv49oMS6rdrGcCvwohvFpcEUJ4HVhI/CZTre424L9K6g4Qv4HMkPSebPUMYBRwe1n924Hpko6ssp+1xBf4DiQdDVxF/Lazrco2Gk2ycWtiycVMUgvxbNwtIYSNVeo2quTiVsHnAAE/qlKuUaQYtwJwEFD+XlufPaaaZyQl1YP8fuKZrXIrgONqqPt6CKEvp+4o4JiScv3AqznlKN+PooKkcdnYg1nAd3P2Pw+4O4TweJV2NqKU43Z8NoZlm6Tlkr5Qpb2NIsWYfZh4tqArG3u0ORvvc28NH1SNIsW45fk8sDSEkPe3NKLk4hZC6AVuA/5e0sckHSzp/cCNxOEpj1Rpt+0FqQ68Hw+sy1nfQxwrs7t1i88XH9eH7HzxMOWKPgH8d/ZzAK4PIXyrtICk84GPAMdWaWOjSjJuwOPEsSgvA+PIzvpIag8hfLtKu1OXYswmZY//DDxIPGPxXmAu8JikD2QfQI0sxbjtQNJHgd8ldj02i1TjdgHwfeDRknVLgD8PIWyt0m7bC1JN5iC+qMqphnqqsW6t5Yp+CZwAHEoceHqppBBCuApA0nji4NErQwjdNbSzUSUVN4AQwjVlde6T9DPgKknfCyFsqqH9KUstZsUeh9eBzxQ/tCS9BvwKOJ94hrzRpRa3crOI3YYLqra4saQYt28T31eXAk8CU4BvAg9KOiWE8Nsa2m97INVkbh07f3uA+M0l75tJqR7iCy2vbvH54uNhyl61w5QDIISwAXgq+/URSVuBqyX9WwjhHeKLvQv4iaRxWbnR2eOhkrY0wQs+xbhVcidxttl0YHGVtqcsxZitzZ57uHR7IYQlkjYCx1dpdyNIMW5DsvFd5wAPhBDWVGlvI0kublmX6hXAhSGE/yjWk7SE2JtxIXBzlbbbHkp1zNwKYr9/ueOA52uoe6SksTl1t/LuOIIVwHuAo3PKUcN+niIe3+IYneOIH/xriW/KdcDl2XNriN14jS7FuFVS/Bbb6PfDSzFmxbE/lWIzWGV7jSDFuJU6k5hcNMvEh6IU4zY9e3yytFAI4RXiJIhpVbZne0Gqydz9wImSjiquyC5O+MfZc9XqjgT+uqRuATgXeCiE0J+t/jnxDTCzrP75wHPZDKPhnEL8MPl19vslxOvvlC7Ff1SnAd+osr1GkGLcKvks8fpNz1Ypl7rkYhZCeJv4gXO6pKGuo2wMVitlHzoNKrm4lZlF/OL7QJVtNJoU4/ab7PEPSwtJ+j3iGOPhejhsb9lf10DZmwvQQvyW8SxxuvaZxFkzv6bkAoXEaxQNANeU1f8x8czYhcQxAHcDW4APlZW7Plv/VeBU4jibQeCTJWU+kdWfRUzQzszKbQfmVfk7rqW5rjOXXNyAk4gfKF/I9vlp4L4sbpfX+5g6ZvnvtWxfA8A9wF8QJ628BbwAjKn3cXXcKv+PJF4HbRvw/XofR8etpv+RBwHPABvY8aLBxTNzU+p9XJthqXsDdrvhcWzAPcRr2/QC9wJTy8pMJX7oXlu2fgxxavVvshf0EuDUnH0cRDxjtpI4lXs58FdlZY4lXmX7raxMF/AE8VvPiCp/w7U0UTKXYtyI0/kfJH677Ac2AYuA8+p9LB2z4d9rxCTuyWy/a4FbgYn1Pp6OW9W4/UPWpg/X+xg6brXFDZhAnOD3CrHH4i3i9e7eV+/j2SyLskCYmZmZWYJSHTNnZmZmZjiZMzMzM0uakzkzMzOzhDmZMzMzM0uYkzkzMzOzhDmZMzMzM0uYkzkzS46kv5EUJB1T77bsb5J+mP3t1ZZT691WM9s/CvVugJmZ7ZJvAT8o+f1C4h1K/oR4df6iavfYNLMG4WTOzGwvkjQSGAj76IrsIYTXgNdK9ndG9uOSEMLAvtinmR3Y3M1qZg1B0mOSnpB0mqSlkvokPSfprJIy52RdkL+fU/9BSc+U/F6Q9HVJL0rql7RK0k2SRpeUmZpt74uS/knSKuKtj8ZJOkLSj7J6/ZI6Jf2PpLaS+mMl3SDpdUlbs8erJPl/s5nVzGfmzKyRHA3cDMwF1gD/CNwt6dgQwqvA/cQbgp8PfK1YSdJE4DTgipJt3Q58EriBeD/eacQuzqnA2WX7vYp4H9jZxPtebgF+Qrwh+mXEe1VOJN78fGy2zwLwC+C4bLvPAicCVwPjs7abmVXlZM7MGsnhwMkhhFcAJC0FOoFzgO+EELZIugv4rKQrQgiDWb3zAAELsnonAecCs0IIt2ZlHpbUA9wu6YMhhGfe3S1dwF+Wdq1K+ihwZQjhjpJyd5X8fB5xnNspIYTHs3WPSAL4pqQbQgjde3Y4zKwZ+FS+mTWSV4qJHECWDHUDU0rK3AZMBv60ZN3ngIdDCJ3Z72cAW4F7su7WQnYm7aHs+ZPL9ntvzhi5J4HLJH1F0nRlWVqJM4CVwKKcfYwknqUzM6vKyZyZNZKenHX9wOiS338JvEFM4JA0DfgQMckragNGAZuAbSVL8UzZhLJ9dLKzc4ndul8DlgPvSLqmZDxcG7EbdlvZ8n8V9mFmlsvdrGbWVEIIQdLtwCWSLiYmdZuAn5UUW0sc93ZShc2sKt9szn66gS8BX5L0PmAWcB2wGpiX7eN1Yhdwnjdq+XvMzJzMmVkzug34BvBpYCZwTwihr+T5nwOXA4eGEB7Z052FEF4CrpT0d8AHSvZxNrAphPDinu7DzJqXkzkzazohhJclLQGuJ46fu63s+cck3UmcCftdYtfnIHEm68eBy0MIL1favqRDgYeBO4AXid2nnwIO491xd3cAFxAnPdwELCN27R4NnAmcVZZgmpnlcjJnZs3qNuBfgXeA/815/nzgy8DfEi890k/s+vwFcfbqcLYAS4GLiOPiBoGXgJkhhPsAQgjbJM0gXg5lNnAk8FviBYEfIE7AMDOrSvvoIuVmZmZmth94NquZmZlZwpzMmZmZmSXMyZyZmZlZwpzMmZmZmSXMyZyZmZlZwpzMmZmZmSXMyZyZmZlZwpzMmZmZmSXs/wFtvMXj4/OoIAAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.figure(figsize=(10, 8))\n", + "plt.rcParams['font.size'] = 16\n", + "\n", + "plt.plot(t_inverse, electrons_per_sec,'-ko')\n", + "plt.ylabel('Dark Current (e-/sec)')\n", + "plt.xlabel('Inverse T')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Fit for the band gap energy\n", + "We will try to fit a model for dark current of the form\n", + "\n", + "$D = \\alpha e^{-e_g/2k_bT}$" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "# The Boltzmann constant\n", + "k_b = 8.6175e-5\n", + "\n", + "\n", + "def dark_current(t_k, alpha, e_g):\n", + " \"\"\"\n", + " Analytic expression for dark current as a function of temperature.\n", + " \n", + " Parameters\n", + " ----------\n", + " t_k : numpy.ndarray\n", + " Temperature in Kelvin\n", + " alpha : float\n", + " Constant coefficient in front of exponential funtion.\n", + " e_g : float\n", + " Band gap energy in eV.\n", + " \n", + " Returns\n", + " -------\n", + " dark_current : numpy.ndarray\n", + " Dark current in electrons/pixel/second.\n", + " \"\"\"\n", + " dark_current = alpha * np.exp(-e_g / (2 * k_b * t_k))\n", + " return dark_current\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Inital guesses for parameters\n", + "We need initial guesses for the values of $\\alpha$ and $e_g$." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "t_0 = t_kel[-1]\n", + "d_0 = electrons_per_sec[-1]\n", + "\n", + "alpha_0 = d_0 / np.exp(-1.1 / (2 * k_b * t_0))\n", + "e_g_0 = 1.1" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Least squares fit" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "p_opt, p_cov = curve_fit(dark_current, t_kel, electrons_per_sec, p0=[alpha_0, e_g_0])\n", + "\n", + "# Errors in the fit\n", + "sig_alpha, sig_e_g = np.sqrt(np.diag(p_cov))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Best fit values" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "77098910165.43013\n", + "1.2205226843061665\n" + ] + }, + { + "data": { + "text/plain": [ + "0.4242143350125111" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "alpha_fit, e_g_fit = p_opt\n", + "print(alpha_fit)\n", + "print(e_g_fit)\n", + "dark_current = alpha_fit * np.exp(-e_g_fit / (2 * k_b * 273.15))\n", + "dark_current" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot the data and best fit model\n", + "\n", + "Make a plot in the cell below. To plot the model, use the `dark_current` function with `alpha_fit` and `e_g_fit`." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "ename": "TypeError", + "evalue": "'numpy.float64' object is not callable", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0;31m#dark_current = alpha_fit * np.exp(-e_g_fit / (2 * k_b * t_kel))\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt_kel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdark_current\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt_kel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m77098910165.43013\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1.2205226843061665\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 7\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mylabel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Dark Current (e-/sec)'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 8\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mxlabel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Temperature (K)'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mTypeError\u001b[0m: 'numpy.float64' object is not callable" + ] + }, + { + "data": { + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.figure(figsize=(10, 8))\n", + "plt.rcParams['font.size'] = 16\n", + "\n", + "#dark_current = alpha_fit * np.exp(-e_g_fit / (2 * k_b * t_kel))\n", + "\n", + "plt.plot(t_kel, dark_current(t_kel, 77098910165.43013, 1.2205226843061665))\n", + "plt.ylabel('Dark Current (e-/sec)')\n", + "plt.xlabel('Temperature (K)')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "anaconda-cloud": {}, + "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.3" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +} diff --git a/gain.ipynb b/gain.ipynb new file mode 100644 index 0000000..f9214ed --- /dev/null +++ b/gain.ipynb @@ -0,0 +1,592 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Detector Gain\n", + "\n", + "This is a Jupyter notebook, which allows us to write and run Python code in a realtime, interactive fashion. There are two main kinds of cells\n", + "\n", + "* Markdown cells like this one. We can make notes and write equations in LaTeX (e.g, $e = mc^2$)\n", + "\n", + "* Code cells like the one below. Type Shift-Enter to run a code cell." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2\n" + ] + } + ], + "source": [ + "# This is a comment.\n", + "print(1 + 1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Importing packages\n", + "We first need to import the packages we'll use to load, analyze, and visualize the data." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# Astropy is \n", + "from astropy.io import fits\n", + "\n", + "# Numpy is a powerful package for numerical analysis.\n", + "import numpy as np\n", + "\n", + "# Matplotlib is the most popular Python plotting package.\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Working with FITS files in Python\n", + "At their most basic, FITS files are just arrays of integers corresponding to the counts in each pixel. Usually FITS files also contain metadata in a header, such as the exposure time, or transformations from pixel coordinates to sky coordinates (RA and Dec.). The astropy fits package allows us to load FITS files as Numpy arrays, which we can then perform calculations on.\n", + "\n", + "Note the call to `astype(np.int32)` at the end of `fits.getdata`. This is very important for getting the correct result for computer science reasons you can ignore if you want to. If you're interested in why this matters, see the end of the notebook.\n", + "\n", + "#### Replace 'data/Flat.15S0X1.V.14.fits' with one of your images." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(495, 657)\n", + "\n", + "[[ 669 823 1091 ... 768 636 652]\n", + " [ 914 1114 1088 ... 1388 1225 1107]\n", + " [1737 932 1323 ... 1145 1170 1377]\n", + " ...\n", + " [2911 2343 2669 ... 3179 3239 3266]\n", + " [2617 2262 2666 ... 3440 3991 3272]\n", + " [2417 2588 2637 ... 3274 3418 3347]]\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "WARNING: Header block contains null bytes instead of spaces for padding, and is not FITS-compliant. Nulls may be replaced with spaces upon writing. [astropy.io.fits.header]\n" + ] + } + ], + "source": [ + "image = fits.getdata('/Users/megankokoris1/Desktop/ccd/flat29.5.FIT').astype(np.int32)\n", + "print(image.shape)\n", + "print('')\n", + "print(image)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can see that the image is just a 1472 x 2184 array of integers." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Load Images\n", + "Using the `fits.getdata` function, Load the two flat images and two bias images that you will use to calculate the gain. Remembers to put `.astype(np.int32)` at the end" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "flat1 = fits.getdata('/Users/megankokoris1/Desktop/ccd/flat29.5.FIT').astype(np.int32)\n", + "flat2 = fits.getdata('/Users/megankokoris1/Desktop/ccd/flat29.52.FIT').astype(np.int32)\n", + "bias1 = fits.getdata('/Users/megankokoris1/Desktop/ccd/bias.FIT').astype(np.int32)\n", + "bias2 = fits.getdata('/Users/megankokoris1/Desktop/ccd/bias1.FIT').astype(np.int32)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can write Equation 3 in code much like it appears in the instructions. You can make variables that are the same as the variables in the equation.\n", + "\n", + "* For the mean: `mean_f1 = np.mean(flat1)`\n", + "* For the squared standard deviation of the image difference ($\\sigma^2$): `sigma_f1f2 = np.std(flat_1 - flat_2)`\n", + "\n", + "In the cell below, write out the equation. Use parentheses to ensure the correct order of operations." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "5958.109835032209\n", + "6739.779149178236\n", + "112.93477545623665\n", + "115.25704841412605\n", + "373.7896456416462\n", + "51.41217361743975\n", + "0.09096956267906531\n" + ] + } + ], + "source": [ + "mean_f1 = np.mean(flat1)\n", + "print(mean_f1)\n", + "\n", + "mean_f2 = np.mean(flat2)\n", + "print(mean_f2)\n", + "\n", + "mean_b1 = np.mean(bias1)\n", + "print(mean_b1)\n", + "\n", + "mean_b2 = np.mean(bias2)\n", + "print(mean_b2)\n", + "\n", + "sigma_f1f2 = np.std(flat1 - flat2)\n", + "print(sigma_f1f2)\n", + "\n", + "sigma_b1b2 = np.std(bias1 - bias2)\n", + "print(sigma_b1b2)\n", + "\n", + "gain = ((mean_f1 + mean_f2) - (mean_b1 + mean_b2))/ (sigma_f1f2**2 - sigma_b1b2**2)\n", + "print(gain)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Read noise\n", + "Calculate the read noise. The square root function is `np.sqrt()`" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "3.307098075421236\n" + ] + } + ], + "source": [ + "read_noise = (gain * sigma_b1b2) / np.sqrt(2)\n", + "\n", + "print(read_noise)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Working with Subregions\n", + "Is your gain value close to the manufacturer's specification (1.3 e-/ADU)? What happens when you calculate the gain for a subregion of the detector? In the example below, we select a 100 x 100 pixel sub region from `row0` to `row1` and `colm0` to `colm1`." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(100, 100)\n", + "\n", + "[[ 94 192 124 ... 66 195 78]\n", + " [ 46 76 93 ... 84 136 102]\n", + " [ 82 172 256 ... 130 151 121]\n", + " ...\n", + " [151 165 131 ... 91 139 125]\n", + " [ 72 61 178 ... 108 130 158]\n", + " [120 106 130 ... 60 71 107]]\n" + ] + } + ], + "source": [ + "row0 = 100\n", + "row1 = 200\n", + "colm0 = 200\n", + "colm1 = 300\n", + "\n", + "# This is called taking a 'slice' of the array, i.e. a subregion\n", + "sub_image = image[row0:row1, colm0:colm1]\n", + "\n", + "print(sub_image.shape)\n", + "print('')\n", + "print(sub_image)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now try calculating the gain in different subregions by varying `row0`, `row1`, `colm0`, and `colm1` Try different sizes, e.g. 100 x 100, 200 x 200, etc" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "6938.0038\n", + "7864.8839\n", + "115.8075\n", + "119.5273\n", + "219.73691343056132\n", + "51.41724698931284\n", + "0.31917985246234\n", + "11.604576384156436\n" + ] + } + ], + "source": [ + "row0 = 100\n", + "row1 = 200\n", + "colm0 = 200\n", + "colm1 = 300\n", + "\n", + "sub_mean_f1 = np.mean(flat1[row0:row1, colm0:colm1])\n", + "print(sub_mean_f1)\n", + "\n", + "sub_mean_f2 = np.mean(flat2[row0:row1, colm0:colm1])\n", + "print(sub_mean_f2)\n", + "\n", + "sub_mean_b1 = np.mean(bias1[row0:row1, colm0:colm1])\n", + "print(sub_mean_b1)\n", + "\n", + "sub_mean_b2 = np.mean(bias2[row0:row1, colm0:colm1])\n", + "print(sub_mean_b2)\n", + "\n", + "sub_sigma_f1f2 = np.std(flat1[row0:row1, colm0:colm1] - flat2[row0:row1, colm0:colm1])\n", + "print(sub_sigma_f1f2)\n", + "\n", + "sub_sigma_b1b2 = np.std(bias1[row0:row1, colm0:colm1] - bias2[row0:row1, colm0:colm1])\n", + "print(sub_sigma_b1b2)\n", + "\n", + "sub_gain = ((sub_mean_f1 + sub_mean_f2) - (sub_mean_b1 + sub_mean_b2))/ (sub_sigma_f1f2**2 - sub_sigma_b1b2**2)\n", + "print(sub_gain)\n", + "\n", + "sub_read_noise = (sub_gain * sub_sigma_b1b2) / np.sqrt(2)\n", + "print(sub_read_noise)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Adding Bells and Whistles\n", + "\n", + "Below is a little more advanced Python programming. You can try running it on your own data and tweaking the code to see how it works. \n", + "\n", + "You probably found varying and typing everything out by hand to be a bit tedious. Alternatively, we can write a function to calculate the gain. The text in triple quotations is called the doc(umentation)string. It tells the user what the function does, what the arguments are, and what the function returns." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "def calculate_gain(flat1, flat2, bias1, bias2):\n", + " \"\"\"\n", + " Calculate detector gain given two flat frames and two bias frames.\n", + " \n", + " Parameters\n", + " ----------\n", + " flat_1, flat_2 : numpy.array_like\n", + " The flat frames\n", + " bias_1, bias_2 : numpy.array_like\n", + " The bias frames\n", + " \n", + " Returns\n", + " -------\n", + " gain : float\n", + " The detector gain\n", + " \"\"\"\n", + " # This is Equation 3 from the assignment\n", + " numerator = (np.mean(flat1) + np.mean(flat2)) - (np.mean(bias1) + np.mean(bias2))\n", + " denominator = np.std((flat1 - flat2)) ** 2 - np.std((bias1 - bias2)) ** 2\n", + " gain = numerator / denominator\n", + " \n", + " return gain" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Replace these files with your own data to try out the code." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "ename": "FileNotFoundError", + "evalue": "[Errno 2] No such file or directory: 'data/Flat.10S0X1.V.15.fits'", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mFileNotFoundError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mflat_1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfits\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgetdata\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'data/Flat.10S0X1.V.15.fits'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mint32\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mflat_2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfits\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgetdata\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'data/Flat.15S0X1.V.14.fits'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mint32\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mbias_1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfits\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgetdata\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'data/Bias(0.0S0X1).02.fits'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mint32\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0mbias_2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfits\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgetdata\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'data/Bias(0.0S0X1).03.fits'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mint32\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/miniconda3/envs/astroconda/lib/python3.7/site-packages/astropy/io/fits/convenience.py\u001b[0m in \u001b[0;36mgetdata\u001b[0;34m(filename, header, lower, upper, view, *args, **kwargs)\u001b[0m\n\u001b[1;32m 187\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mclosed\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_get_file_mode\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 188\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 189\u001b[0;31m \u001b[0mhdulist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mextidx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_getext\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 190\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 191\u001b[0m \u001b[0mhdu\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mhdulist\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mextidx\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/miniconda3/envs/astroconda/lib/python3.7/site-packages/astropy/io/fits/convenience.py\u001b[0m in \u001b[0;36m_getext\u001b[0;34m(filename, mode, ext, extname, extver, *args, **kwargs)\u001b[0m\n\u001b[1;32m 1022\u001b[0m \u001b[0;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'extver alone cannot specify an extension.'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1023\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1024\u001b[0;31m \u001b[0mhdulist\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfitsopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 1025\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1026\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mhdulist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mext\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/miniconda3/envs/astroconda/lib/python3.7/site-packages/astropy/io/fits/hdu/hdulist.py\u001b[0m in \u001b[0;36mfitsopen\u001b[0;34m(name, mode, memmap, save_backup, cache, lazy_load_hdus, **kwargs)\u001b[0m\n\u001b[1;32m 149\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 150\u001b[0m return HDUList.fromfile(name, mode, memmap, save_backup, cache,\n\u001b[0;32m--> 151\u001b[0;31m lazy_load_hdus, **kwargs)\n\u001b[0m\u001b[1;32m 152\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 153\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/miniconda3/envs/astroconda/lib/python3.7/site-packages/astropy/io/fits/hdu/hdulist.py\u001b[0m in \u001b[0;36mfromfile\u001b[0;34m(cls, fileobj, mode, memmap, save_backup, cache, lazy_load_hdus, **kwargs)\u001b[0m\n\u001b[1;32m 388\u001b[0m return cls._readfrom(fileobj=fileobj, mode=mode, memmap=memmap,\n\u001b[1;32m 389\u001b[0m \u001b[0msave_backup\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msave_backup\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcache\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcache\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 390\u001b[0;31m lazy_load_hdus=lazy_load_hdus, **kwargs)\n\u001b[0m\u001b[1;32m 391\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 392\u001b[0m \u001b[0;34m@\u001b[0m\u001b[0mclassmethod\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/miniconda3/envs/astroconda/lib/python3.7/site-packages/astropy/io/fits/hdu/hdulist.py\u001b[0m in \u001b[0;36m_readfrom\u001b[0;34m(cls, fileobj, data, mode, memmap, save_backup, cache, lazy_load_hdus, **kwargs)\u001b[0m\n\u001b[1;32m 1037\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfileobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_File\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1038\u001b[0m \u001b[0;31m# instantiate a FITS file object (ffo)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1039\u001b[0;31m \u001b[0mfileobj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_File\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfileobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmemmap\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmemmap\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcache\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcache\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 1040\u001b[0m \u001b[0;31m# The Astropy mode is determined by the _File initializer if the\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1041\u001b[0m \u001b[0;31m# supplied mode was None\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/miniconda3/envs/astroconda/lib/python3.7/site-packages/astropy/utils/decorators.py\u001b[0m in \u001b[0;36mwrapper\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m 519\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mnew_name\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 520\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 521\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mfunction\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 522\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 523\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mwrapper\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/miniconda3/envs/astroconda/lib/python3.7/site-packages/astropy/io/fits/file.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, fileobj, mode, memmap, overwrite, cache)\u001b[0m\n\u001b[1;32m 176\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_open_fileobj\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfileobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moverwrite\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 177\u001b[0m \u001b[0;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfileobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 178\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_open_filename\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfileobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moverwrite\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 179\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 180\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_open_filelike\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfileobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moverwrite\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/miniconda3/envs/astroconda/lib/python3.7/site-packages/astropy/io/fits/file.py\u001b[0m in \u001b[0;36m_open_filename\u001b[0;34m(self, filename, mode, overwrite)\u001b[0m\n\u001b[1;32m 553\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 554\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_try_read_compressed\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmagic\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mext\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mext\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 555\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_file\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfileobj_open\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mIO_FITS_MODES\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mmode\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 556\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclose_on_error\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 557\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/miniconda3/envs/astroconda/lib/python3.7/site-packages/astropy/io/fits/util.py\u001b[0m in \u001b[0;36mfileobj_open\u001b[0;34m(filename, mode)\u001b[0m\n\u001b[1;32m 386\u001b[0m \"\"\"\n\u001b[1;32m 387\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 388\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbuffering\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 389\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 390\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'data/Flat.10S0X1.V.15.fits'" + ] + } + ], + "source": [ + "flat_1 = fits.getdata('data/Flat.10S0X1.V.15.fits').astype(np.int32)\n", + "flat_2 = fits.getdata('data/Flat.15S0X1.V.14.fits').astype(np.int32)\n", + "bias_1 = fits.getdata('data/Bias(0.0S0X1).02.fits').astype(np.int32)\n", + "bias_2 = fits.getdata('data/Bias(0.0S0X1).03.fits').astype(np.int32)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Calculating the gain in different subregions\n", + "We can divide the detector up into a grid of subregions. The `bins` are the row and column boundaries of the subregions. Effectively we are making a coarse \"map\" of the gain calculated on different parts of the detector." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "row_bins = np.linspace(0, flat_1.shape[0], 5).astype(int)\n", + "print(row_bins)\n", + "col_bins = np.linspace(0, flat_1.shape[1], 10).astype(int)\n", + "print(col_bins)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# This is an array to store the gain in each subregion.\n", + "gain_map = np.zeros((len(row_bins) - 1, len(col_bins - 1)))\n", + "\n", + "# This nested for loop goes through each subregion.\n", + "for ii in range(len(row_bins) - 1):\n", + " for jj in range(len(col_bins) - 1):\n", + " row_slice = slice(row_bins[ii], row_bins[ii + 1])\n", + " col_slice = slice(col_bins[jj], col_bins[jj + 1])\n", + " \n", + " # The `local_gain` is the gain in the subregion\n", + " local_gain = calculate_gain(flat_1[row_slice, col_slice], flat_2[row_slice, col_slice],\n", + " bias_1[row_slice, col_slice], bias_2[row_slice, col_slice])\n", + " \n", + " # Store the local gain in the `gain_map`\n", + " gain_map[ii, jj] = local_gain" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot the gain values" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=(20, 10))\n", + "\n", + "# Use the imshow function to plot one of the flat images for reference.\n", + "plt.imshow(flat_1, vmin=np.percentile(flat_1, 5), vmax=np.percentile(flat_1, 90),\n", + " origin='lower', cmap='binary_r', interpolation='nearest')\n", + "\n", + "\n", + "# Plot the boundaries of the subregions\n", + "for row in row_bins:\n", + " plt.axhline(row)\n", + "for col in col_bins:\n", + " plt.axvline(col)\n", + "\n", + "# Print the local gain value in each subregion\n", + "for ii in range(len(row_bins) - 1):\n", + " for jj in range(len(col_bins) - 1):\n", + " row_loc = (row_bins[ii] + row_bins[ii + 1]) / 2\n", + " col_loc = (col_bins[jj] + col_bins[jj + 1]) / 2\n", + " \n", + " plt.text(col_loc, row_loc, '{:.3f}'.format(gain_map[ii, jj]),\n", + " ha='center', va='center', fontsize=14, color='r')\n", + "\n", + "plt.xlim(0, col_bins[-1])\n", + "plt.ylim(0, row_bins[-1])\n", + "\n", + "# Add a colorbar\n", + "plt.colorbar()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Unsigned versus signed integers\n", + "\n", + "As their names suggest, unsigned integers can only represent non-negative numbers. That's 0 up to $2^n - 1$, where $n$ is the number of bits the variable takes up in memory. Signed integers can have values from $-(2^{n-1})$ up to $2^{n - 1} - 1$.\n", + "\n", + "Negative values have no real physical meaning on an astronomical image, so the negative values are a waste of memory. Therefore images are usually stored as unsigned integers.\n", + "\n", + "The problem is when we need to subtract two sets of unsigned integers, like we did in the gain calculation. This can give weird results if we're not careful." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "a = np.array([1], dtype=np.uint16)\n", + "b = np.array([3], dtype=np.uint16)\n", + "\n", + "print('a =', a)\n", + "print('b =', b)\n", + "print('a - b =', a - b)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can fix this by converting the unsigned 16-bit integers to signed 32-bit integers." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "ename": "NameError", + "evalue": "name 'a' is not defined", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0ma\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mint32\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mb\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mint32\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'a ='\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'b ='\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mNameError\u001b[0m: name 'a' is not defined" + ] + } + ], + "source": [ + "a = a.astype(np.int32)\n", + "b = b.astype(np.int32)\n", + "\n", + "print('a =', a)\n", + "print('b =', b)\n", + "print('a - b =', a - b)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "anaconda-cloud": {}, + "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.3" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +} diff --git a/linearity.ipynb b/linearity.ipynb new file mode 100644 index 0000000..c9aa1fb --- /dev/null +++ b/linearity.ipynb @@ -0,0 +1,258 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Detector Gain\n", + "\n", + "You could calculate the mean of each flat image manually using IRAF. Alternatively, you could calculate the means for all of the images automatically using Python." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from astropy.io import fits\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "# glob serves some of the same functions as ls in the terminal\n", + "import glob" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## FITS Headers\n", + "The headers of the FITS files contain the exposure times of the flat images. Now we use `fits.open` instead of `fits.getdata`. HDU stands for Header/Data Unit." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "29.5\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "WARNING: Header block contains null bytes instead of spaces for padding, and is not FITS-compliant. Nulls may be replaced with spaces upon writing. [astropy.io.fits.header]\n" + ] + } + ], + "source": [ + "hdu = fits.open('/Users/megankokoris1/Desktop/ccd/flat29.5.FIT')\n", + "header = hdu[0].header\n", + "print(header['exposure'])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Calculating Mean Counts\n", + "We can find all of the flat images, assuming they all have 'Flat' in the name.\n", + "#### You will need to change the path to the directory containing your data." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "['/Users/megankokoris1/Desktop/ccd/Xavier/12_sec.FIT', '/Users/megankokoris1/Desktop/ccd/Xavier/15_sec.FIT', '/Users/megankokoris1/Desktop/ccd/Xavier/18_sec.FIT', '/Users/megankokoris1/Desktop/ccd/Xavier/18_sec2.FIT', '/Users/megankokoris1/Desktop/ccd/Xavier/24_sec.FIT', '/Users/megankokoris1/Desktop/ccd/Xavier/30_sec.FIT', '/Users/megankokoris1/Desktop/ccd/Xavier/36_sec.FIT', '/Users/megankokoris1/Desktop/ccd/Xavier/3_sec.FIT', '/Users/megankokoris1/Desktop/ccd/Xavier/42_sec.FIT', '/Users/megankokoris1/Desktop/ccd/Xavier/48_sec.FIT', '/Users/megankokoris1/Desktop/ccd/Xavier/54_sec.FIT', '/Users/megankokoris1/Desktop/ccd/Xavier/60_sec.FIT', '/Users/megankokoris1/Desktop/ccd/Xavier/6_sec.FIT', '/Users/megankokoris1/Desktop/ccd/Xavier/80_sec.FIT']\n" + ] + } + ], + "source": [ + "# This is equivalent to $ ls Flat*.fits\n", + "#flat_list = glob.glob('/Users/megankokoris1/Desktop/ccd/flat/flat*.FIT')\n", + "flat_list = glob.glob('/Users/megankokoris1/Desktop/ccd/Xavier/*.FIT')\n", + "print(flat_list)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we can loop through each flat image, and keep track of the exposure time and mean counts" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Exposure time 12.0 sec\n", + "Mean counts: 18178.72\n", + "\n", + "Exposure time 15.0 sec\n", + "Mean counts: 23983.32\n", + "\n", + "Exposure time 18.0 sec\n", + "Mean counts: 28975.72\n", + "\n", + "Exposure time 18.0 sec\n", + "Mean counts: 26611.05\n", + "\n", + "Exposure time 24.0 sec\n", + "Mean counts: 36629.32\n", + "\n", + "Exposure time 30.0 sec\n", + "Mean counts: 45259.31\n", + "\n", + "Exposure time 36.0 sec\n", + "Mean counts: 52164.65\n", + "\n", + "Exposure time 3.0 sec\n", + "Mean counts: 5353.18\n", + "\n", + "Exposure time 42.0 sec\n", + "Mean counts: 56429.01\n", + "\n", + "Exposure time 48.0 sec\n", + "Mean counts: 57330.16\n", + "\n", + "Exposure time 54.0 sec\n", + "Mean counts: 57626.20\n", + "\n", + "Exposure time 60.0 sec\n", + "Mean counts: 57723.53\n", + "\n", + "Exposure time 6.0 sec\n", + "Mean counts: 9578.54\n", + "\n", + "Exposure time 80.0 sec\n", + "Mean counts: 57786.30\n", + "\n" + ] + } + ], + "source": [ + "# These are empty lists (arrays) to store the exposure times and mean counts\n", + "exp_times = []\n", + "means = []\n", + "\n", + "for filename in flat_list:\n", + " # Open the FITS file\n", + " hdu = fits.open(filename)\n", + " \n", + " exptime = hdu[0].header['exptime']\n", + " print('Exposure time {} sec'.format(exptime))\n", + " \n", + " # This will append the exposure time for each image to the array\n", + " exp_times.append(exptime)\n", + " \n", + " # Same for mean counts\n", + " mean_counts = np.mean(hdu[0].data)\n", + " print('Mean counts: {:.2f}\\n'.format(mean_counts))\n", + " means.append(mean_counts)\n", + "\n", + "# Convert to Numpy arrays so they can be sorted\n", + "exp_times = np.array(exp_times)\n", + "means = np.array(means)\n", + "\n", + "# Sort by exposure time so the plot looks correct\n", + "time_sort = np.argsort(exp_times)\n", + "exp_times = exp_times[time_sort]\n", + "means = means[time_sort]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot mean counts versus exposure time\n" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAowAAAHtCAYAAACTcy+0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xl4VOX5xvHvEwhoVBQQd5NQ0bohtU37U2u1YhW1CC7gRtW6EMsioEgVwioGRUHEBRRUqiWIYqWIdSuo1KqoQUUBl6qQuIHIIntY8vz+mAkOwyQZYDJnMrk/15VrmDPPnLmTC8nje973PebuiIiIiIhUJiPoACIiIiKS2tQwioiIiEiV1DCKiIiISJXUMIqIiIhIldQwioiIiEiV1DCKiIiISJXUMIqIiIhIlQJpGM3sHDP7j5mtMbNVZlZsZq0jXm9sZg+b2Q9mttbMZphZyxjn2c3M7jKz78xsvZm9ZWanxKjLMLO+ZrbIzDaY2Vwzu7CSbJ3N7BMzKzOzT83sL4n97kVERERql6Q3jGZ2HTANmAOcD3QEpgBZ4dcNeBY4C7geuBDIBF41s0OiTvcI0BkYCLQFvgNeMrNfRNUNBQYD9wNnA7OBKWZ2TlS2zsBDwD/Cnz8FGGNmXXb1+xYRERGprSyZd3oxs1zgY6Cvu99TSU174J9Aa3d/NXxsb2AhMNHde4SPtQI+AK529wnhY/WB+cCn7t4ufGw/4CvgDncfFPE5M4Fm7n5cxHu/BV5w9ysj6h4F2gEHuvum6r7Hfffd13Nzc+P9kYiIiIgEZs6cOT+4e7Pq6uonI0yEq4Fy4MEqatoB31Y0iwDu/qOZTQfaAz0i6jYBT0bUbTazycAtZtbQ3cuANkADYGLU50wEHjWz5u6+EDgRaBaj7u/AVcDJwKtUIzc3l+Li4urKRERERAJnZiXx1CX7kvTJwCfAJWb2hZltNrPPzaxbRM0xwLwY750PZJvZnhF1C919XYy6BkCLiLoy4PMYdQBHR9QR47Oj60RERETqlGQ3jAcBhwN3AXcAZwL/Bu43s57hmibAihjvXR5+bBxnXZOIx5W+/bX3WHXEOGd03XbMLD+8cKd46dKllZWJiIiI1ErJbhgzgL2A69x9vLu/4u5dgBeBvuEFLwbEmlhpMZ4nuo5Kaqvk7uPcPc/d85o1q3YagIiIiEitkuyGcVn48d9Rx18G9gcOJDSiF2s0r2JksWIEsLq65RGPjcPNaHV1xDhnk6jXRUREROqUZDeM8ys5XtHMlYdrjolRczRQ6u5rIs7V3MyyYtRt5Kc5i/OBhsBhMeoAFkRli/7s6DoRERGROiXZDePU8GObqONtgK/dfTGhPRgPNrNTK140s0bAueHXKjxLaH/GjhF19YGLgZfDK6QhdLl7I9Ap6jP/BMwLr5AGeAv4oZK65cAbcX6PIiIiImkl2dvqPE9oa5qHzGxf4EugA6HFL1eFa54l1LxNNLM+hC5B9yU0CnlnxYnc/QMzexK4x8wyCe3T2AVoTkTT5+7fm9koQnMkVwPvEWoqWxPapqeibpOZDSC0Ufc3wIxwzdXA9e6+MdE/DBEREZHaIKkNo7u7mZ0H3A4MITSP8BOgk7tPCteUm1lbYAQwBtiNUAN5mrt/FXXKq4BC4DZgH2AucJa7vxdVVwCsAXoCBwCfAhe5+/SofA+amQO9gT5AKdDd3cck4vsXERERqY2SeqeXuiAvL8+1cbeIiIjUBmY2x93zqqtL+r2kRURERKR2UcMoIiIiIlVSwygiIiIiVVLDKCIiIiJVUsMoIiIiIlVSwygiIiIiVVLDKCIiErCioiJyc3PJyMggNzeXoqKioCNJgFLx74MaRhERSXup+Au4QlFREfn5+ZSUlODulJSUkJ+fn1IZJX7uTnl5OeXl5WzZsoXNmzezadMmNm7cSFlZGRs2bGD9+vWsW7eOtWvXsmbNGlavXs2qVatYtWoVDz/8MJ07d065vw/auDvBtHG3iEhqqWjI1q1bt/VYVlYW48aNo1OnTtvUbtmyhU2bNrFp06atv+hr+vmYMWNYvXr1drn33HNPLr74Yty92i8grrpdeU8yPiNVc8VbX5NycnJYtGhRws8b78bdahgTTA2jiNRVRUVFFBQUUFpaSnZ2NoWFhds1ZImwceNGfvzxx60jMhV/ruxx6tSprF+/frvzZGRksNdee23TwCXzd6KZkZmZycaNGyutOeiggzCzuL4qzrkjXzv6nnT5jFTO1bt370r/vpSXl+/S37lKzhtXw5jUe0mLiEh6ih7Fq7iMBmxtGrds2cLq1au3a+jiafoiH8vKyqrN06BBA/bee28aNWoUs1kEKC8v58orryQzM5P69euTmZm5zVf0sUQ/z8gIzQrLzc2lpKRku3w1NaIkqe3ee++N+fchOzs7gDQ/0QhjgmmEUUTqmmXLlnHMMcewZMmS7V6rX78+++23H6tWrWLNmjXVnisjI2Nro1fxGPnnyh6jjzVs2HDrOVO9IduRS+aS/pL990EjjCIiklCrV69mwYIFzJs3b5uvxYsXV/qezZs3c/bZZ8fd6GVlZW29lJcohYWFMX8BFxYWJvRzdlZFE5CMy/mS+lL174NGGBNMI4wiUtuVlZXxySefbNcYRo7GZWVlccwxx3Dsscdy7LHHMnz4cL7//vvtzpVKo3ip9gtYJBVohFFERKq0efNmvvjii+0aw//9739s2bIFgMzMTI488khOPPFEOnfuvLVBrNiipsL++++f8qN4ahBFdp4aRhGRNOfulJaWbtcYfvzxx1sXkJgZhx12GMceeywdOnTY2hgeccQRZGZmVvsZqXoZTUQSQ5ekE0yXpEWkpsRzWXXJkiXbNYbz58/fZp+/Qw45ZGtDWPF11FFHkZWVlexvSUQCpn0YA6KGUURqQqyVkw0bNqRTp07sscceW5vDpUuXbn29adOmtGzZcpvG8JhjjmGfffYJ4lsQkRSkhjEgahhFpCZUtjUMhO4IEj1ieOyxx7LffvslfMWxiKQXLXoREUkjpaWlMY+bGT/++OM2C1BERBJN/8KIiKS4b775ptKFJ9nZ2WoWRaTG6V8ZEZEUNnv2bPLy8jCzbe5eAqm1bY2IpDc1jCIiKWrChAmceuqpZGVlUVxczCOPPEJOTg5mRk5Ojm4dJyJJozmMIiIpZvPmzdx0002MHj2a008/nSeffJKmTZty7LHHqkEUkUBohFFEJIUsW7aMNm3aMHr0aHr16sWLL75I06ZNg44lInWcRhhFRFLERx99RPv27fnmm2+YMGECf/7zn4OOJCICaIRRRCQlTJ06lRNPPJH169cza9YsNYsiklLUMIqIBKi8vJwhQ4ZwwQUXcMwxx1BcXMwJJ5wQdCwRkW3okrSISEDWrFnDFVdcwdSpU7niiit46KGH2G233YKOJSKyHTWMIiIB+PLLL2nfvj0LFixg1KhR9OzZU7fxE5GUpYZRRCTJZs6cyUUXXYS78+KLL3LGGWcEHUlEpEqawygikiTuzr333kubNm044IADeOedd9QsikitoIZRRCQJysrKuOaaa+jZsydt27Zl9uzZtGjRIuhYIiJxUcMoIlLDvvvuO37/+98zYcIEBgwYwDPPPMNee+0VdCwRkbhpDqOISA165513OP/881m5ciVTpkyhQ4cOQUcSEdlhGmEUEakhjz/+OKeccgoNGjTgzTffVLMoIrWWGkYRkQTbvHkzvXv35sorr+TEE0/k3XffpVWrVkHHEhHZabokLSKSQCtWrOCSSy7h5Zdfpnv37tx9991kZmYGHUtEZJeoYRQRSZAFCxbQvn17SkpKGD9+PNdee23QkUREEkINo4hIAjz77LN06tSJPfbYg1dffZXf/va3QUcSEUkYzWEUEdkF7k5hYSHnnXceP//5z3n33XfVLIpI2tEIo4jITlq7di1XXXUVU6ZMoVOnTowfP57dd9896FgiIgmnhlFEZCcsWrSI9u3bM2/ePO666y569+6NmQUdS0SkRqhhFBHZQa+99hodOnRg8+bN/Otf/+Kss84KOpKISI3SHEYRkTi5O2PGjOGMM86gWbNmvPPOO2oWRaROUMMoIhKHjRs3ct1119GtWzfatGnD7NmzOeKII4KOJSKSFGoYRUSqsWTJElq3bs348ePp27cv06ZNY++99w46lohI0mgOo4hIFebMmcN5553HsmXLmDx5MhdffHHQkUREkk4jjCIiEYqKisjNzSUjI4NmzZpxwgknkJGRwRtvvKFmUUTqLI0wioiEFRUVkZ+fz7p16wD44YcfyMjI4JZbbuH4448POJ2ISHA0wigiElZQULC1WaxQXl7O8OHDA0okIpIa1DCKiISVlpbu0HERkbpCDaOISNi+++4b83h2dnaSk4iIpBY1jCIiwBdffMGaNWvIyNj2n8WsrCwKCwsDSiUikhrUMIpInbdhwwY6duxIw4YNGTlyJDk5OZgZOTk5jBs3jk6dOgUdUUQkUFolLSJ13g033MD777/PtGnTaNeuHb169Qo6kohIStEIo4jUaZMmTeLBBx+kT58+tGvXLug4IiIpSQ2jiNRZH3/8Mfn5+Zx88smapygiUgU1jCJSJ61du5aOHTuy++67M3nyZDIzM4OOJCKSsjSHUUTqHHena9euLFiwgJdeeomDDz446EgiIilNI4wiUuc8+uijPP744wwcOJAzzjgj6DgiIilPDaOI1Clz586le/fu/OEPf2DAgAFBxxERqRWS3jCa2e/NzGN8rYyqa2xmD5vZD2a21sxmmFnLGOfbzczuMrPvzGy9mb1lZqfEqMsws75mtsjMNpjZXDO7sJKMnc3sEzMrM7NPzewvifsJiEhQVq1aRceOHWncuDFFRUXUq1cv6EgiIrVCkHMYewDvRjzfXPEHMzPgWaA5cD2wAugLvGpmv3D3ryPe9wjwR6AP8CXQDXjJzE509w8i6oYCNwEFwBzgEmCKmbV19+cjPrsz8BBwOzADOB0YY2bm7mMT8p2LSNK5O9deey1ffvklr7zyCvvtt1/QkUREao0gG8aP3X12Ja+1A04GWrv7qwBm9hawEPgroWYTM2sFXAZc7e4TwsdmAfOBW8Pnwcz2I9Qs3uHuI8Kf8aqZtQDuAJ4P19UHCoG/u3tBRN1BwFAze9jdNyXqByAiyfPAAw8wZcoU7rjjDk45ZbuLECIiUoVUncPYDvi2olkEcPcfgelA+6i6TcCTEXWbgclAGzNrGD7cBmgATIz6nIlASzNrHn5+ItAsRt3fgaaEmlgRqWXeffddbrzxRtq2bUufPn2CjiMiUusE2TAWmdkWM1tmZpPMLDvitWOAeTHeMx/INrM9I+oWuvu6GHUNgBYRdWXA5zHqAI6OqCPGZ0fXiUgtsXz5cjp27MiBBx7IY489RkZGqv5/sohI6grikvSPwEhgFrAKOB7oB7xlZse7+/dAE2BRjPcuDz82BtaE61ZUUdck4nGlu3scdcQ4Z3TdNswsH8gHyM7OjlUiIgEoLy/nyiuv5Ntvv+W///0vTZrE/E9YRESqkfSG0d3fB96PODTLzP4DvENobmJ/wIDo5o7w8ejnia6jktpKufs4YBxAXl7eDr1XRGrOyJEjee6557j33nv5zW9+E3QcEZFaKyWuzbj7e8BnwK/Dh5YTezSvcfhxRZx1yyMeG4dXX1dXR4xzNol6XURS3Ouvv07fvn3p0KED3bt3DzqOiEitlhINY1jkKOB8fppPGOlooNTd10TUNTezrBh1G/lpzuJ8oCFwWIw6gAURdcT47Og6EUlh33//PZdccgnNmzfn4YcfZvv/VxQRkR2REg2jmeUBRwBvhw89CxxsZqdG1DQCzg2/RkRdJtAxoq4+cDHwsruXhQ+/SKiB7BT10X8C5rn7wvDzt4AfKqlbDryxM9+fiCTPli1b6NSpE8uWLePpp59m7733DjqSiEitl/Q5jGZWRGg/xfeAlYQWvfQFvgHuC5c9S6h5m2hmffhp424D7qw4l7t/YGZPAveYWWb4vF0IbfjdKaLuezMbBfQ1s9Xhz74YaE3ENj3uvsnMBhDaqPsbQht3twauBq53940J/nGISILddtttzJgxg/Hjx9OqVaug44iIpIUgVknPAy4ldAeXLGAx8AwwyN1/AHD3cjNrC4wAxgC7EWogT3P3r6LOdxWhzbZvA/YB5gJnhedFRiogtLK6J3AA8ClwkbtPjyxy9wfNzIHehO4eUwp0d/cxCfjeRaQGzZgxgyFDhnD55ZdzzTXXBB1HRCRt2PY7zciuyMvL8+Li4qBjiNQ53377Lb/4xS9o1qwZ77zzDnvssUfQkUREUp6ZzXH3vOrqUmIOo4jIrti8eTOXXHIJ69at4+mnn1azKCKSYEHeS1pEJCH69+/P66+/TlFREUcddVTQcURE0o5GGEWkVnvuuecYPnw41113HZdddlnQcURE0pIaRhGptRYtWsQVV1zB8ccfzz333BN0HBGRtKWGUURqpY0bN3LRRRexZcsWpkyZwm677RZ0JBGRtKU5jCJSK9100028++67/OMf/+Cww6Jv4iQiIomkEUYRqXWmTJnCfffdR69evbjggguCjiMikvbUMIpIrfK///2Pa665hhNOOIHhw4cHHUdEpE5Qwygitcb69evp0KEDmZmZPPnkkzRo0CDoSCIidYLmMIpIrdGjRw8+/PBDnn/+ebKzs4OOIyJSZ2iEUURqhccff5yHH36Yfv36cfbZZwcdR0SkTlHDKCIpb/78+XTp0oVTTz2VIUOGBB1HRKTOUcMoIiltzZo1dOjQgb322osnnniC+vU1k0ZEJNn0L6+IpCx357rrruOzzz5jxowZHHjggUFHEhGpk9QwikjKGjduHJMmTWLo0KGcdtppQccREamzdElaRFLSe++9R48ePWjTpg39+vULOo6ISJ2mhlFEUs6PP/5Ix44dadasGRMnTiQjQ/9UiYgESZekRSSluDtXXXUVpaWlzJo1i3333TfoSCIidZ4aRhFJKaNHj2bq1KmMHDmSk046Keg4IiKCLkmLSAqZPXs2ffr0oX379txwww1BxxERkTA1jCKSEpYtW8ZFF13EoYceyt/+9jfMLOhIIiISpkvSIhK48vJyLr/8cpYsWcKbb77JPvvsE3QkERGJoIZRRAI3fPhwXnjhBcaMGcOvfvWroOOIiEgUXZIWkUC99tpr9O/fn0suuYS//OUvQccREZEY1DCKSGAWL17MpZdeyuGHH864ceM0b1FEJEWpYRSRpCoqKiI3N5eMjAxyc3NZtmwZU6ZMYa+99go6moiIVEINo4gkTVFREfn5+ZSUlODulJWVYWZ8+OGHQUcTEZEqqGEUkaQpKChg3bp12xzbuHEjBQUFASUSEZF4qGEUkaQpLS3doeMiIpIa1DCKSNJkZ2fv0HEREUkNahhFJGm6d+++3bGsrCwKCwsDSCMiIvFSwygiSeHuvPjii+y+++4ccsghmBk5OTmMGzeOTp06BR1PRESqoDu9iEhSPPXUU8ycOZMHHniArl27Bh1HRER2gLl70BnSSl5enhcXFwcdQySlrFq1iiOPPJKDDjqIt99+m3r16gUdSUREADOb4+551dVphFFEatzgwYNZvHgx//znP9UsiojUQprDKCI16sMPP+Tee+8lPz+f3/zmN0HHERGRnaCGUURqTHl5OV26dKFx48YMGzYs6DgiIrKTdElaRGrMY489xptvvsmjjz5KkyZNgo4jIiI7SSOMIlIjli9fzl//+ld++9vfcuWVVwYdR0REdoEaRhGpEf369WPFihWMGTOGjAz9UyMiUpvpX3ERSbh33nmHcePG0aNHD4477rig44iIyC5SwygiCbVlyxa6dOnCAQccwODBg4OOIyIiCaBFLyKSUA899BDvvfcekydPplGjRkHHERGRBNAIo4gkzJIlS+jXrx9/+MMfuOiii4KOIyIiCaKGUUQS5q9//Svr1q3j/vvvx8yCjiMiIgmihlFEEuI///kPjz/+OH369OHnP/950HFERCSBzN2DzpBW8vLyvLi4OOgYIkm1adMmjj/+eNasWcOCBQvIysoKOpKIiMTBzOa4e151dVr0IiK7bPTo0cyfP59p06apWRQRSUO6JC0iu+Trr79m8ODBnHvuubRr1y7oOCIiUgPUMIrILrnhhhvYsmULo0ePDjqKiIjUEDWMIrLTXnrpJZ5++mn69+9P8+bNg44jIiI1RIteEkyLXqSu2LBhAy1btiQjI4MPP/yQhg0bBh1JRER2kBa9iEiNuuuuu/j88895+eWX1SyKiKQ5XZIWkR325ZdfMmzYMC666CLOOOOMoOOIiEgNU8MoIjvE3bn++uupX78+d999d9BxREQkCXRJWkR2yLRp03j++ecZOXIkBx98cNBxREQkCbToJcG06EXS2dq1azn66KPZe++9mTNnDpmZmUFHEhGRXaBFLyKScLfddhulpaW8/vrrahZFROoQzWEUkbh8/PHHjBw5kj//+c+cfPLJQccREZEkUsMoItVyd7p168aee+7JnXfeGXQcERFJMl2SFpFqTZ48mVdffZWxY8fSrFmzoOOIiEiSaYRRRKr0448/cuONN/LrX/+azp07Bx1HREQCoBFGEanSoEGDWLJkCdOnT6devXpBxxERkQAEPsJoZi+amZvZbVHHG5vZw2b2g5mtNbMZZtYyxvt3M7O7zOw7M1tvZm+Z2Skx6jLMrK+ZLTKzDWY218wurCRTZzP7xMzKzOxTM/tL4r5jkdrjgw8+4L777qNLly7k5VW764KIiKSpQBtGM7sUaBXjuAHPAmcB1wMXApnAq2Z2SFT5I0BnYCDQFvgOeMnMfhFVNxQYDNwPnA3MBqaY2TlRn90ZeAj4R/jzpwBjzKzLTn+jIrVQeXk5Xbt2pWnTptx2223Vv0FERNJWYJekzWwfYBRwAzAp6uV2wMlAa3d/NVz/FrAQ+CvQI3ysFXAZcLW7TwgfmwXMB24Nnwcz2w+4CbjD3UeEP+NVM2sB3AE8H66rDxQCf3f3goi6g4ChZvawu29K6A9CJEVNmDCBt956i8cee4zGjRsHHUdERAIU5AjjncB8d38ixmvtgG8rmkUAd/8RmA60j6rbBDwZUbcZmAy0MbOG4cNtgAbAxKjPmQi0NLPm4ecnAs1i1P0daEqoiRVJe8uWLePmm2/md7/7HZdffnnQcUREJGCBNIxmdjJwBdC1kpJjgHkxjs8Hss1sz4i6he6+LkZdA6BFRF0Z8HmMOoCjI+qI8dnRdSJprW/fvqxcuZIxY8YQmiEiIiJ1WdIbRjPLJDRHcIS7f1pJWRNgRYzjy8OPjeOsaxLxuNK3v3F2rDpinDO6bhtmlm9mxWZWvHTp0lglIrXG7Nmzefjhh+nVqxfHHnts0HFERCQFBDHCeDOwO6G5gpUxILq5qzhe03VUUlspdx/n7nnunqdNjaU227JlC127duWggw5i0KBBQccREZEUkdRFL2aWDRQA1wINI+YYEn6+D7Ca0IherNG8ipHFihHA5UB2FXXLIx4bm5lFjTLGqiP82d9F1DWJel0kLY0dO5b333+fp556ir322ivoOCIikiKSPcL4M2A3QotKVkR8QWgV8wqgJaE5g8fEeP/RQKm7rwk/nw80N7OsGHUb+WnO4nygIXBYjDqABRF1xPjs6DqRtLN48WIKCgo488wz6dChQ9BxREQkhSS7YfwAOC3GF4SayNMINXnPAgeb2akVbzSzRsC54dcqPEtof8aOEXX1gYuBl929LHz4RUINZKeoPH8C5rn7wvDzt4AfKqlbDryxY9+uSO3Rp08fNmzYwP3336+FLiIiso2kXpJ295XAa9HHw7+cStz9tfDzZwk1bxPNrA+hkce+hOYY3hlxvg/M7EngnvBimoVAF6A5EU2fu39vZqOAvma2GniPUFPZmohtetx9k5kNILRR9zfAjHDN1cD17r4xMT8JkdTy2muvMXHiRAYMGMDhhx8edBwREUkxKXkvaXcvN7O2wAhgDKHL2G8Bp7n7V1HlVxFaQHMbsA8wFzjL3d+LqisA1gA9gQOAT4GL3H161Gc/aGYO9Ab6AKVAd3cfk8BvUSRlbNy4kW7dutG8eXP69u0bdBwREUlBtv1OM7Ir8vLyvLi4OOgYInG78847ufnmm5k+fTpt27YNOo6IiCSRmc1x97zq6gK9l7SIBOurr75iyJAhtG/fXs2iiIhUSg2jSB3Wq1cv3J3Ro0cHHUVERFJYSs5hFJGa98ILL/DMM88wbNgwcnJygo4jIiIpTHMYE0xzGKU22LBhA8ceeyyZmZnMnTuXBg0aBB1JREQCEO8cRo0witRBw4cP54svvmDmzJlqFkVEpFq7NIfRzJomKoiIJMcXX3zB7bffzqWXXkrr1q2DjiMiIrVAXA2jmXUOb6Bd8bylmX0NfG9mxWZ2QI0lFJGEcXe6d+9OgwYNGDlyZNBxRESkloh3hPF6YH3E87uBlUAvYG/g1gTnEpEaMHXqVF588UWGDh3KgQceGHQcERGpJeKdw5gNfAJgZnsDpwLnufvzZrYMuL2G8olIgqxZs4ZevXrRqlUrunXrFnQcERGpReJtGOsB5eE/nww4P90T+itgv8TGEpFEGzp0KF999RWTJ0+mfn2tdxMRkfjFe0n6f8Afw3++BHjT3deFnx8ELE90MBFJnAULFnD33Xdz9dVXc9JJJwUdR0REapl4hxlGAH83syuBxkDHiNdOAz5MdDARSQx3p2vXrjRq1Ijhw4cHHUdERGqhuBpGd59kZiXACcC77v6fiJeXANNqIpyI7LpJkyYxa9Ysxo0bx7777ht0HBERqYXiutOLmZ0CvOfua2K8tifwy6gmss7SnV4klaxcuZIjjzyS3Nxc3nzzTTIydPt4ERH5Sbx3eon3t8erwNGVvPbz8OsikiKKiorIzc2lcePGLFmyhHPPPVfNooiI7LR4f4NYFa81BLYkIIuIJEBRURH5+fmUlJRsPTZs2DCKiooCTCUiIrVZpZekzSwX+Fn46QygO+G9GCPsDlwNtHL3FjUTsXbRJWkJWm5u7jbNYoWcnBwWLVqU/EAiIpKy4r0kXdWilyuBQYT2XHTgPrYdafTw882AdgEWSRGrFjKfAAAgAElEQVSlpaU7dFxERKQ6VTWMfyO0ObcBrxBqChdE1ZQBn7m79mEUSRGHHnpozOYwOzs7gDQiIpIOKm0Y3b0EKAEws9MIrZJenaxgIrJzzjzzTB5++OFtjmVlZVFYWBhQIhERqe3iWvTi7rPULIqkvrKyMl566SUOO+wwsrOzMTNycnIYN24cnTp1CjqeiIjUUnFt3G1mDYC+wKVANqGV0ZHc3XVzWpGAPfTQQ3z11VfMnDmT1q1bBx1HRETSRLxN3l2E5jC+ADxDaO6iiKSQNWvWUFhYyOmnn65mUUREEirehrEDMMjdNQlKJEXde++9fP/995qrKCIiCRfvxt17Am/VZBAR2XkrVqzgzjvvpH379vzf//1f0HFERCTNxNswTgdOqckgIrLz7rrrLlatWsXQoUODjiIiImko3kvS9wGPm1k58Dyw3b6L7v5lIoOJSHwWL17M6NGjueyyy2jZsmXQcUREJA3F2zBWXI4eTOjuL7HU2+U0IrLDhg0bxsaNGxk8eHDQUUREJE3F2zBeTehWgCKSQkpKSnjwwQe5+uqradFCt3MXEZGaEVfD6O5/q+EcIrIThgwZQkZGBgMGDAg6ioiIpLF4F72ISIr55JNPeOyxx+jWrRuHHHJI0HFERCSNxXunl0erKXF3vyYBeUQkTgMHDiQrK4tbbrkl6CgiIpLm4p3D2Jrt5zA2AfYCVoa/RCRJ3nvvPaZMmcLAgQNp1qxZ0HFERCTNxTuHMTfWcTM7BXgQ6JTATCJSjf79+9OkSRNuvPHGoKOIiEgdsEtzGN39P8AoQvs0ikgSvP7667zwwgvccsst7L333kHHERGROiARi16+BI5PwHlEpBruTr9+/TjwwAPp1q1b0HFERKSOiHcOY0xmVh/4M/B1QtKISJVeeukl/vvf/zJmzBiysrKCjiMiInVEvKukX4lxuAFwBNAU+EsiQ4nI9srLy+nXrx/Nmzfnmmu0KYGIiCRPvCOMGWy/Sno18Aww2d1fS2QoEdneM888w/vvv8/jjz9OgwYNgo4jIiJ1iLnrjn+JlJeX58XFxUHHkDSzefNmWrZsSUZGBh9++CH16unW7SIisuvMbI6751VXt0tzGEUkOSZOnMgnn3zCM888o2ZRRESSLu5V0mbW0syeNrOlZrbZzL43s6fMrGVNBhSp68rKyhg8eDB5eXmcd955QccREZE6KN5FL78GZgHrgWeBxcABwLnAH83sFHefU2MpReqw8ePHU1JSwvjx4zGzoOOIiEgdFNccRjObATQCTnf31RHH9wJmAD+6+5k1lrIW0RxGSaS1a9dy2GGHcdRRR/HKK6+oYRQRkYRK9BzGE4DLI5tFAHdfbWbDgcd2IqOIVOP+++9nyZIlPPPMM2oWRUQkMPHOYaxuGFJLrUUSbOXKlQwfPpy2bdty0kknBR1HRETqsHgbxreBfuFL0FuZ2R7AzcDsRAcTqetGjhzJihUruO2224KOIiIidVy8l6T7Aa8BJWb2HPAdoUUvfwR2B35fE+FE6qrvv/+eUaNGcckll9CqVaug44iISB0XV8Po7u+Y2QnAQKAN0ARYDrwCDHX3j2ouokjdc/vtt7NhwwaGDBkSdBQREZH4N+529w+BDjWYRUSA0tJSxowZw5///GeOOOKIoOOIiIhUPofRzDLM7FwzO7aKmpZmdm7NRBOpm4YOHQrAwIEDA04iIiISUtWilz8BTwBrq6hZDTxhZpcmNJVIHfXZZ58xYcIEunTpQnZ2dtBxREREgOobxgnuvrCyAndfBDwCXJngXCJ10qBBg9htt93o27dv0FFERES2qqph/CXwchznmAFUu0O4iFRt7ty5TJ48mV69erH//vsHHUdERGSrqhrGvYAVcZxjRbhWRHZB//792WeffbjpppuCjiIiIrKNqhrGH4CcOM6RHa4VkZ305ptv8txzz3HzzTezzz77BB1HRERkG1U1jP8lvrmJfw7XishOcHf69evH/vvvz/XXXx90HBERke1U1TDeA5xuZqPMrEH0i2aWaWajgdbAqJoKKJLuZsyYwaxZs+jfvz977LFH0HFERES2Y+5e+YtmvYCRwDJCC2BKwi/lAGcATYHe7j66hnPWGnl5eV5cXBx0DKkl3J3f/OY3LF26lE8//ZSGDRsGHUlEROoQM5vj7tUuXq7yTi/ufo+ZvQfcApxP6L7RAOsJ3Vv6Dnd/fRezitRZ//znPykuLmbChAlqFkVEJGVVOcK4TaFZBrBv+Okyd99SY6lqMY0wSry2bNnCcccdR3l5OR999BH168d9p04REZGESMgIYyR3Lwe+36VUIrLVpEmTWLBgAVOmTFGzKCIiKa2qRS81wszamNkrZrbYzMrM7Gsze8rMjo6qO9TMnjazH81slZk9Y2bb3SvNzBqb2cNm9oOZrTWzGWbWMkbdbmZ2l5l9Z2brzewtMzslRl2GmfU1s0VmtsHM5prZhYn9KUhdt3HjRgYNGsQvf/lLLrjggqDjiIiIVCnpDSPQBJgDdAfOBPoCxwCzzSwHwMyygFeAIwlt7XM5cDjwqpltXUZqZgY8C5wFXA9cCGSG6w6J+txHgM7AQKAt8B3wkpn9IqpuKDAYuB84G5gNTDGzcxLwvYsA8Mgjj7Bw4UIKCwvJyAjiP0MREZH4xT2HsUZDmP0c+AS4yd1HmllP4G7g5+7+ebimOfA/4K/ufnf4WHvgn0Brd381fGxvYCEw0d17hI+1Aj4Arnb3CeFj9YH5wKfu3i58bD/gK0KLeQZF5JsJNHP346r7XjSHUaqzbt06WrRoQYsWLZg1axah/+8RERFJvnjnMKbK0May8OOm8GM7YHZFswjg7guBN4D2Ee9rB3xb0SyG634Epseo2wQ8GVG3GZgMtDGziuWpbYAGwMSofBOBluGmVWSXPPDAA3z33XcUFhaqWRQRkVohsIbRzOqZWQMzOxx4CFhMqIGD0CXqeTHeNh+InOtYVV22me0ZUbfQ3dfFqGsAtIioKwM+j1FH1GeL7LAff/yRO+64g7PPPpvf/e53QccRERGJS9xLM82sEXAOoXtH7xb1srv70B387LeBX4X//Dmhy8oVq7CbACtivGc50DjieRNgUSV1hGvXVHO+ivNUPK707a/TR9dtw8zygXyA7Ozt1uWIbHX33XezfPlybrvttqCjiIiIxC2uhtHMfkvoMu8+lZQ4ocUiO+JyoBHwM+Am4N9mdrK7L4o453ZRYjwPom4b7j4OGAehOYxV1UrdtXTpUu6++246duzIL3/5y6DjiIiIxC3eS9L3EBrJ+zWwm7tnRH3V29EPdveP3f1td38COB3Yk9AdZSA0GhhrNK8x244ULq+ijoja6uqWRzw2tu0nlkXXieywO+64g3Xr1nHrrbcGHUVERGSHxNswHgX0d/c57r4x0SHcfSWhy9IVcwnnE5pPGO1oYEHE86rqSt19TURd8/B2PdF1G/lpzuJ8oCFwWIw6oj5bJG5ff/01DzzwAFdccQVHHnlk0HFERER2SLwNYymhRqpGmNn+hPZc/CJ86FngBDP7WURNLvDb8GtE1B1sZqdG1DUCzo1Rlwl0jKirD1wMvOzuZeHDLxJqIDtFRfwTMC+8Ultkh912222Ul5czaNCg6otFRERSTLwN4xDglnAztkvMbKqZDTCz9mZ2mpldB8wCNgMjw2XjCV0CnxauawdMI7RH4kMRp3sWeAuYaGaXmFmb8DED7qwocvcPCG2pc4+ZXWtmpxNakd0cGBRR9z0wCuhrZjea2e/NbCzQGui3q9+7pI+ioiJyc3PJyMggNzeXoqKiSms///xzHnnkEa677jpyc3OTF1JERCRB4l0l3RbYH1hoZm+x/Vw+d/cr4zzXbOAioDehLW2+Al4Dbq9Y8OLua82sNaHm7e+EGsCZQK+Iy8y4e7mZtQVGAGMIrd5+CzjN3b+K+tyrgELgNkKLd+YCZ7n7e1F1BYRWVvcEDgA+BS5y9+lxfn+S5oqKisjPz2fdutAuTSUlJeTn5wPQqVP04DQMHjyYzMxMCgoKkppTREQkUeK604uZVXcp1t39Z9XU1Am600v6y83NpaSkZLvjOTk5LFq0aJtjH330Ea1ateLmm2/m9ttvT1JCERGR+MR7p5e4RhjdXXc4EQkrLS2N+/iAAQNo1KgRffr0qelYIiIiNSZVbg0oUmtUtjl79PG3336badOm0adPH5o0ibnnu4iISK2www2jme1nZtnRXzURTiQV9e7de7tjWVlZFBYWbnOsoKCAZs2a0bNnz2RFExERqRHx3uklg9Bikeuo/G4vO7x5t0htNG/ePDIyMsjMzKSsrIycnBwKCwu3WfAyc+ZMZs6cyT333MOee+5ZxdlERERSX7yLXm4ktP3McEKNYyFQTmi/wnLgDnd/tAZz1hpa9JLePv/8c4488ki6dOnC2rVrmTFjxnZzF92dE088kW+//ZbPPvuM3XaLvvW6iIhIaoh30Uu8l6SvAm4l1DACTHX3QYTuAPMNoEvSUicMHDiQhg0bVrlFzvTp03n77bcZNGiQmkUREUkL8TaMPwOK3X0LoQ22dwdw902E7jN9dc3EE0kdc+fO5YknnqBnz54ccMABMWvKy8spKCjg8MMP58or492aVEREJLXFu3H3j4Q2xQb4Fvg58EbEObQEVNLegAED2GeffarcImfy5MnMmzePyZMnU79+vP95iYiIpLZ4f6O9DxwNvBT+GmJm6wmNNhYC0XdLEUkrb775JtOnT2fYsGE0btw4Zs2mTZsYOHAgrVq1omPHjjFrREREaqN4G8Z7CF2WhtDil18CFTfPLQG6JziXSMpwd/r168f+++9Pjx49Kq2bMGECX3zxBc899xwZGdriVERE0ke8d3r5d8SfF5vZb4DDgCzg4/BcRpG09O9//5tZs2Zx3333sccee8SsWb9+PbfeeisnnXQS55xzTpITioiI1KydmmTlob14Pk9wFpGUUzG6mJOTQ+fOnSutGzt2LN988w1FRUWYWRITioiI1Ly4r5uZ2cFmdreZFZvZQjM7Nny8l5n9X81FFAnO1KlTmTNnDoMHD6Zhw4Yxa1atWsWwYcM488wzOfXUU5OcUEREpObF1TCa2THAR8DlhFZJZwMNwi/nALr3maSdLVu20L9/f4466iguv/zySuvuueceli1btt2tAUVERNJFvJekRwIfA22ADcDGiNfe5KcNvUXSxsSJE/n44495+umnqVcv9p0vV69ezYgRI7jgggvIy6t2o3wREZFaKd6G8WTgUndfY2bRvzmXALF3MRappcrKyhg0aBC/+tWvuOCCC7Z7vaioiKeeeoq1a9cC8H//p1kZIiKSvuKdw1hexWv7AusTkEUkZYwfP56SkhKGDRu23SKWoqIi8vPztzaLAEOGDKGoqCj6NCIiImnBQgueqykymwGscvcLwiOMm4A8d3/PzCYDWe7eroaz1gp5eXleXFwcdAzZBWvXruWwww7jqKOO4pVXXtmuYczNzaWkpGS79+Xk5LBo0aIkpRQREdl1ZjbH3audUxXvJemhwAwzexmYBDjwBzPrCZwPnLLTSUVSzL333suSJUt45plnYm6RU1paGvN9lR0XERGp7eK6JO3us4DzgObAo4ABdwC/A85z97drLKFIEq1YsYI777yTtm3bctJJJ8Wsyc7O3qHjIiIitV3c+zC6+7/c/XDgCEKLYI5y95+5+ws1lk4kye666y5WrlxZ5RY5hYWFZGVlbXMsKytL2+qIiEja2uEb3rr75+7+prt/WhOBRIKyePFiRo8ezaWXXspxxx1XaV2nTp0YN24cOTk5mBk5OTmMGzeOTp06JTGtiIhI8lS66MXMWu/Iidz9lYQkquW06KX26tGjB2PGjOGTTz6hRYsWQccRERGpcYlY9DKD0OIWCM1ZjMXDrzkQe2djkVpg0aJFPPjgg1xzzTVqFkVERKJUt0p6NfCP8NfaampFaq0hQ4aQkZHBgAEDgo4iIiKScqpqGE8DrgAuBDoCU4HHdOlZ0s2CBQt4/PHHueGGGzjkkEOCjiMiIpJyKl304u6z3P0aQrf9+wuwH/CSmZWa2e1mdlSyQorUpIEDB7LHHntwyy23BB1FREQkJVW7StrdN7j7JHc/G8gGRgPnAPPM7P6aDihSk4qLi/nHP/5B79692XfffYOOIyIikpJ2dFudZcCi8JcDjROcRySpCgoKaNq0KTfccEPQUURERFJWXLcGNLPfApcTmsvYEJgG/BH4d81FE6lZr732Gi+//DIjRoygUaNGQccRERFJWZU2jGbWglCT+CcgF/gPcBMwxd3XJCWdSA1xd/r168fBBx9M165dg44jIiKS0qoaYfwMWAU8A1wLlISP72dm+0UXu/uXiY8nUjP+9a9/8dZbb/HQQw+x++67Bx1HREQkpVV1p5fyiKexiyIL3LVxN7rTS21QXl7O8ccfz7p161iwYAGZmZlBRxIREQlEIu70clUC84ikjCeffJIPP/yQSZMmqVkUERGJQ6UjjLJzNMKY2jZt2sTRRx9NVlYW77//PhkZO7pRgIiISPqId4RRvy2lTigqKiI3N5cGDRrw+eefc8YZZ6hZFBERiZN+Y0raKyoqIj8/n5KSkq3Hxo4dS1FRUYCpREREag81jJL2CgoKWLdu3TbH1q1bR0FBQUCJREREahc1jJL2SktLd+i4iIiIbEsNo6S97OzsHTouIiIi21LDKGmvsLBwuwUuWVlZFBYWBpRIRESkdlHDKGkvJyeH8vJymjRpgpmRk5PDuHHj6NSpU9DRREREaoWqNu4WSQsjR46kadOmlJaWkpWVFXQcERGRWkcjjJLWPvvsM6ZNm0bXrl3VLIqIiOwkNYyS1kaNGkWDBg3o1q1b0FFERERqLTWMkraWLl3K3/72Ny6//HL233//oOOIiIjUWmoYJW2NHTuWDRs2cOONNwYdRUREpFZTwyhpacOGDdx///388Y9/5Kijjgo6joiISK2mhlHS0t///neWLl3KTTfdFHQUERGRWk8No6Sd8vJyRo4cyS9/+UtOPfXUoOOIiIjUetqHUdLO888/z6effsqkSZMws6DjiIiI1HoaYZS0M2LECLKzs+nQoUPQUURERNKCGkZJK8XFxcyaNYuePXuSmZkZdBwREZG0oIZR0srIkSNp1KgR1157bdBRRERE0oYaRkkbJSUlTJkyhfz8fBo1ahR0HBERkbShhlHSxujRozEzevToEXQUERGRtKKGUdLCypUrGT9+PJdccgmHHnpo0HFERETSihpGSQvjxo1jzZo19O7dO+goIiIiaUcNo9R6Gzdu5N577+X000/nF7/4RdBxRERE0o427pZa78knn+Sbb75h/PjxQUcRERFJSxphlFrN3Rk5ciRHH300Z511VtBxRERE0pJGGKVWmzlzJnPnzuXRRx/VbQBFRERqSFJHGM2sg5n9w8xKzGy9mX1qZreb2V5RdY3N7GEz+8HM1prZDDNrGeN8u5nZXWb2Xfh8b5nZKTHqMsysr5ktMrMNZjbXzC6sJGNnM/vEzMrC+f6SuJ+AJNqIESM44IADuOyyy4KOIiIikraSfUn6JmAL0A84CxgLdAH+bWYZABYaJno2/Pr1wIVAJvCqmR0Sdb5HgM7AQKAt8B3wkplFr3wYCgwG7gfOBmYDU8zsnMgiM+sMPAT8I/z5U4AxZtZlV79xSbx58+bx0ksvcf3119OwYcOg44iIiKQtc/fkfZhZM3dfGnXsCuAx4HR3f8XM2gP/BFq7+6vhmr2BhcBEd+8RPtYK+AC42t0nhI/VB+YDn7p7u/Cx/YCvgDvcfVDE584Emrn7cRHv/RZ4wd2vjKh7FGgHHOjum6r7HvPy8ry4uHgnfjqyo6666iqeeuopvvrqK5o0aRJ0HBERkVrHzOa4e151dUkdYYxuFsPeDT8eHH5sB3xb0SyG3/cjMB1oH/G+dsAm4MmIus3AZKCNmVUMObUBGgAToz53ItDSzJqHn58INItR93egKXBydd+fJM93331HUVERV199tZpFERGRGpYKq6RPDT9+HH48BpgXo24+kG1me0bULXT3dTHqGgAtIurKgM9j1AEcHVFHjM+OrpMUcN9997FlyxZ69eoVdBQREZG0F2jDaGYHA7cCM9y94jpuE2BFjPLl4cfGcdY1iXhc6dtfe49VR4xzRtdJwNasWcPYsWM5//zzOeyww4KOIyIikvYCaxjDI4XTgM3AVZEvAbEmVkbvmVITdVRSWyUzyzezYjMrXro01lV3SaQJEyawcuVKbrrppqCjiIiI1AmBNIxmthuhldA/A9q4+9cRLy8n9mhexcjiijjrlkc8NrbtN+mLVUeMczaJen077j7O3fPcPa9Zs2aVlUkCbN68mVGjRnHSSSdxwgknBB1HRESkTkh6w2hmmYS2rfkNcI67fxRVMp+f5hNGOhoodfc1EXXNzSwrRt1GfpqzOB9oCERfu6yYk7ggoo4Ynx1dJwGaOnUqCxcu1OiiiIhIEiV74+4MoAg4HWjv7rNjlD0LHGxmp0a8rxFwbvi1yLpMoGNEXX3gYuBldy8LH36RUAPZKepz/gTMc/eF4edvAT9UUrcceCPOb1NqiLszYsQIWrRoQbt27YKOIyIiUmck+9aADxBq8AqBtWYWeU3x6/Cl6WcJNW8TzawPoUvQfQnNMbyzotjdPzCzJ4F7wqOWCwltAt6ciKbP3b83s1FAXzNbDbxHqKlsTcQ2Pe6+ycwGENqo+xtgRrjmauB6d9+Y2B+F7Kg33niDd955hwceeIB69eoFHUdERKTOSPbG3YuAnEpeHuLug8N1TYARwHnAboQayBvdfW7U+XYn1HxeBuwDzAVudvfXourqEWo6OwMHAJ8Ct7r70zEyXgf0DucsBUa5+5h4v0dt3F1zzj//fF5//XVKS0vJyoqeiSAiIiI7Kt6Nu5PaMNYFahhrxmeffcaRRx5JQUEBQ4cODTqOiIhIWkjJO72I7KxRo0aRmZlJ9+7dg44iIiJS56hhlJS3dOlS/va3v3HFFVew//77Bx1HRESkzlHDKClv7NixbNiwgRtvvDHoKCIiInWSGkZJWUVFRWRnZzNo0CB233133nvvvaAjiYiI1EnJ3lZHJC5FRUXk5+ezbt06ANavX09+fj4AnTpFb5UpIiIiNUkjjJKSCgoKtjaLFdatW0dBQUFAiUREROouNYySkkpLS3fouIiIiNQcNYySkg4++OCYx7Ozs5OcRERERNQwSko6/PDDtzuWlZVFYWFhAGlERETqNjWMknJmz57Na6+9xjnnnENOTg5mRk5ODuPGjdOCFxERkQBolbSklM2bN9O1a1cOOuggJk+ezF577RV0JBERkTpPDaOklLFjx/L+++8zZcoUNYsiIiIpQpekJWV899139O/fnzZt2nDhhRcGHUdERETC1DBKyujduzdlZWXcf//9mFnQcURERCRMDaOkhJkzZ/LEE09wyy230KJFi6DjiIiISARz96AzpJW8vDwvLi4OOkatUlZWRqtWrdi8eTMfffQRu+++e9CRRERE6gQzm+PuedXVadGLBG7kyJF8+umnPP/882oWRUREUpAuSUugFi5cyNChQ7nwwgs5++yzg44jIiIiMahhlED17NmTevXqMWrUqKCjiIiISCV0SVoC8+yzzzJ9+nTuuusuDj300KDjiIiISCW06CXBtOglPmvXruWYY45hzz335P333yczMzPoSCIiInWOFr1ISissLKSkpIT//Oc/ahZFRERSnOYwStJ9/PHHjBgxgiuvvJLf/e53QccRERGRaqhhlKRyd7p168Yee+zBnXfeGXQcERERiYMuSUtS/X979x4vdVXvf/z1BgHFy0EIU1OEMjtKoRX2M/OhVKZlHtTMIinTHnFT0/BoXrapoUjeES8oaYfK7SXxrkfwEhJ2vEAoCh4pPVzEK4oixlX4/P5Y353DNHvYwN7z3Xvm/Xw8vo/Zs75rZj7fxXeYz6y1vmtuueUWJk+ezNixY9luu+3yDsfMzMyawD2MVjFLlizhlFNOYe+992bQoEF5h2NmZmZN5B5Gq5hf/vKXLFq0iAceeID27dvnHY6ZmZk1kXsYrSJmzJjBNddcw/HHH88Xv/jFvMMxMzOzDeCE0Vrc2rVrGTZsGN27d+f888/POxwzMzPbQB6SthZ3ww038PTTT3PTTTfRpUuXvMMxMzOzDeQeRmtRixYt4owzzqBfv34cffTReYdjZmZmG8EJo7WoX/ziFyxdupRrr70WSXmHY2ZmZhvBCaO1mKlTpzJ+/HhOPfVUdt9997zDMTMzs42kiMg7hqrSt2/fmD59et5h5G716tV84Qtf4P333+eFF15gyy23zDskMzMzKyLprxHRd331fNGLtYgxY8Ywa9Ys7r77bieLZmZmbZyHpK3ZLVy4kHPPPZdDDz2U/v375x2OmZmZbSInjNbshg8fzpo1axgzZowvdDEzM6sCThitWU2cOJEJEyZw9tln06tXr7zDMTMzs2bgi16aWS1f9LJixQo++9nP0r59e5577jk6deqUd0hmZmZWhi96sYq76KKLePnll3n44YedLJqZmVURD0lbs3jppZcYNWoUAwYM4MADD8w7HDMzM2tGThhtk0UEJ554Ih07duSyyy7LOxwzMzNrZh6Stk125513MmnSJEaPHs2OO+6YdzhmZmbWzHzRSzOrtYteli5dyu6770737t2ZNm0am23m7yBmZmZthS96sYoYMWIEr776KhMmTHCyaGZmVqU8h9E22qxZs7jiiisYNGgQ++yzT97hmJmZWQtxwmgbJSIYNmwYXbp0YdSoUXmHY2ZmZi3IY4i2UX7/+9/z+OOPc+ONN9KtW7e8wzEzM7MW5B5G22CLFy/mtNNOY9999+XYY4/NOxwzMzNrYU4YbYOdddZZLF68mLFjx9KunU8hMzOzaudPe9sgTz31FCyo/jAAABkPSURBVOPGjeOkk06iT58+eYdjZmZmFeCE0ZpszZo1DBs2jB122IHzzjsv73DMzMysQnzRizXZ2LFjeeaZZ7jtttvYZptt8g7HzMzMKsQ9jNYkb7zxBnV1dXzjG9/gqKOOyjscMzMzqyAnjNYkp556KitWrODqq69GUt7hmJmZWQU5YbT1mjx5MvX19Zx++unstttueYdjZmZmFaaIyDuGqtK3b9+YPn163mE0m1WrVrHnnnuycuVKZs+ezRZbbJF3SGZmZtZMJP01Ivqur54verGyLr/8cl588UUeeOABJ4tmZmY1ykPS1qj58+czYsQIjjjiCA455JC8wzEzM7OcOGG0Rp188slIYvTo0XmHYmZmZjnykLSVdN9993HPPfdw0UUX0aNHj7zDMTMzsxz5opdmVg0XvSxbtozevXvTuXNnnn32WTp06JB3SGZmZtYCfNGLbbQLL7yQefPmMWXKFCeLZmZmVvk5jJJ2knSVpCckLZMUknqWqLe5pEskvS5peVZ//xL12kk6U9I8SSskzZR0ZCOvPUjSi5JWSpojaWgj9Q6X9Ez2fPMlnS2p/aYee1swZ84cLr74Yo455hj23/9fmtvMzMxqUB4XvewKfA94F5hapt6NwCDgHOBQ4HVgkqS9iuqdD5wHXA18C3gSuF3SOpf1ShoEXA/cAXwTuB24VtKwonoHZ3WmZc93JXA2cOEGHmebExGccMIJbLnlllx88cV5h2NmZmatRMXnMEpqFxFrs79/CvwG6BUR8wrq7Ak8C/wkIv4rK9sMmA3MiYj+Wdl2wCvAryPi3ILHPwp0j4g+BY99DXgwIn5cUO+3QH9gh4hYnZU9A7wfEQcU1DuHlDT2iIg3yh1fW57DeOutt/KDH/yAa6+9lmHDhq3/AWZmZtamNXUOY8V7GBuSxfXoD6wGbit43IfArcDBkjplxQcDHYGbih5/E/A5Sb2y+18Gupeo9wegG7AfgKSdgb0aqdeB1ONYlZYsWcLw4cPp27cvgwcPzjscMzMza0Va6zqMvYG5EbGsqHw2KUHctaDeSuClEvUA9iioBzBrY+pFxFxgWUG9qnPuuefy5ptvMnbsWNq3r4npmmZmZtZErTVh7Eqa41hsccH+htv34l/H1UvVo8RzNrVeQ1nXEuVIGixpuqTpixYtKlWlVaqvr6dnz560a9eOK6+8kq9//ev07bveXmkzMzOrMa01YRRQanKlNqEejdRtar3i5/yniBgXEX0jom/37t3X8xKtQ319PYMHD2b+/Pk05Nt/+ctfqK+vzzkyMzMza21aa8K4mNK9edsW7G+43VZScTJXqh4lnrNrE+sBdCnY3+bV1dWxbNm6I/7Lly+nrq4up4jMzMystWqtCeNsoJekzkXlewCr+GjO4mygE/CpEvUAXiioBx/NUdygetk6kZ0L6rV5CxYs2KByMzMzq12tNWG8l3RV8lENBdnSON8HHoqIlVnxRFICObDo8T8EZmUXqwA8AbzdSL3FwF8AImIBMLORequBBzf+kFqXxn4f2r8bbWZmZsVy+WlASd/N/vxidvstSYuARRExJSKelXQbMFpSB2AuMAzoRUEyFxFvSboCOFPSUmAGKan8GnBYQb3Vkn5JWqj7VeCRrM5PgJ9FxKqC8M4C7pd0PXAL8HnSGoxXrm8Nxrbk6KOPZtSoUeuUde7cmZEjR+YUkZmZmbVWFV+4G0BSYy86JSL6ZXW2AEYCR5PmD84ETo+Ix4qeqz1wJulXYbYH5gAjImJCidcdAvwnsAuwALgiIq4tUe87wLnAvwNvAjcAIyNizfqOrS0s3L1y5Ur69OnDe++9R6dOnVi4cCE9evRg5MiRDBxY3LlqZmZm1aqpC3fnkjBWs7aQMI4YMYJzzz2XSZMmcdBBB+UdjpmZmeWk1f7Si+Xr73//OxdeeCEDBgxwsmhmZmZN4oSxhkQEw4YNo1OnTlx++eV5h2NmZmZtRC4XvVg+br75Zh599FGuueYadthhh7zDMTMzszbCPYw14t133+WUU07hS1/6EkOGDMk7HDMzM2tD3MNYI84880zeeecdJk2aRPv27fMOx8zMzNoQ9zDWgCeeeILrr7+ek08+mb322ivvcMzMzKyNccJY5VavXs2QIUPYaaed+NWvfpV3OGZmZtYGeUi6yo0ePZrnn3+eu+66i6222irvcMzMzKwNcg9jFZs/fz7nnXce/fv35/DDD887HDMzM2ujnDBWqYjgxBNPRBJXXXVV3uGYmZlZG+Yh6Sp19913c//993PppZfSo0ePvMMxMzOzNsy/Jd3MWsNvSS9dupTdd9+dbt26MX36dDp06JBrPGZmZtY6NfW3pN3DWIXOOeccXnvtNSZMmOBk0czMzDaZ5zBWmRkzZjBmzBiGDBnCPvvsk3c4ZmZmVgWcMFaRNWvWMHToULp3786oUaPyDsfMzMyqhIekq8h1113HtGnTuPnmm+nSpUve4ZiZmVmVcA9jlXjttdc466yzOPDAAxkwYEDe4ZiZmVkVccJYJYYPH87KlSsZO3YskvIOx8zMzKqIE8YqMHHiRP74xz9SV1fHrrvumnc4ZmZmVmW8DmMzq/Q6jMuXL6d379507NiRmTNn0qlTp4q9tpmZmbVtXoexRlxwwQXMnTuXyZMnO1k0MzOzFuEh6TbshRde4JJLLuGYY46hX79+eYdjZmZmVcoJYxu1du1ahg4dytZbb82ll16adzhmZmZWxTwk3UaNHz+eqVOncsMNN9C9e/e8wzEzM7Mq5h7GNqS+vp6ePXvSrl07Bg0axG677cZxxx2Xd1hmZmZW5ZwwthH19fUMHjyY+fPnExGsXbuWBQsWcMstt+QdmpmZmVU5J4xtRF1dHcuWLVunbMWKFdTV1eUUkZmZmdUKJ4xtxIIFCzao3MzMzKy5OGFsI3r06LFB5WZmZmbNxQljGzFy5Eg6d+68Tlnnzp0ZOXJkThGZmZlZrXDC2EYMHDiQcePGscsuuyCJXXbZhXHjxjFw4MC8QzMzM7Mq59+SbmaV/i1pMzMzs43V1N+Sdg+jmZmZmZXlhNHMzMzMynLCaGZmZmZlOWE0MzMzs7KcMJqZmZlZWU4YzczMzKwsJ4xmZmZmVpYTRjMzMzMrywmjmZmZmZXlhNHMzMzMynLCaGZmZmZlOWE0MzMzs7KcMJqZmZlZWU4YzczMzKwsJ4xmZmZmVpYiIu8YqoqkRcD8Ers+Brxd4XBaG7dB4nZI3A6J2yFxOyRuh8TtkFSiHXaJiO7rq+SEsUIkTY+IvnnHkSe3QeJ2SNwOidshcTskbofE7ZC0pnbwkLSZmZmZleWE0czMzMzKcsJYOePyDqAVcBskbofE7ZC4HRK3Q+J2SNwOSatpB89hNDMzM7Oy3MNoZmZmZmU5YTQzMzOzspwwthBJO0uaIGmJpPcl3SmpR95xtSRJO0m6StITkpZJCkk9S9TbXNIlkl6XtDyrv3/lI25+kr4r6Q5J87NjmyNplKSti+ptK+kGSW9L+oekRyR9Lq+4m5ukgyX9SdIbklZKWijpj5L2KKpXi++Tidl744Ki8qo9JyT1y465eHuvqF7VtkEhSYdI+rOkD7LzfrqkrxXsr+p2kPRYI+dDSJpYUK+q2wFA0lckPSTprexcmCHpJ0V1WsVnphPGFiCpM/An4N+BHwM/Aj4NTJa0ZZ6xtbBdge8B7wJTy9S7ERgEnAMcCrwOTJK0V4tH2PJOBdYAZwHfBMYCw4CHJbUDkCTg3mz/z4AjgQ6k82OnPIJuAV2BvwInAgcBZwK9gScl7QK1+T6R9ANgzxLltXBOAJwEfLlgO7BhR620gaQhwD2k98cRwFHA7UDnbH8ttMPxrHsefBk4Jdt3L9RGO0jqAzxCOq5BpGOcBtwoaVhB1dbxmRkR3pp5A04mJQ27FpT1Aj4ETsk7vhY87nYFf/8UCKBnUZ09s/LjCso2A+YA9+Z9DM3QBt1LlB2THfPXsvuHZfe/WlDn34DFwJi8j6EF2+Yz2XH/Z3a/pt4nQBfgDeAHWTtcULCvqs8JoF92fAeWqVPVbZAdT09gOfDzWm6HRo77RmAl0LVW2gG4EFgFbFVU/iTwRPZ3q/nMdA9jy+gPPBkRLzUURMRc4C+kN0FVioi1TajWH1gN3FbwuA+BW4GDJXVqofAqIiIWlSielt1+IrvtD7wWEZMLHrcEuI8qPj+Ad7Lb1dltrb1PLgZmR8QtJfbV6jlRqBba4CfAWuC6MnVqoR3WIWkLUk/rfRGxOCuuhXboSPr/cHlR+Xt8NALcaj4znTC2jN7ArBLls4E9SpTXkt7A3IhYVlQ+m/Tm2bXyIbW4A7Lb/81uy50fPSRtVZGoKkBSe0kdJX0auJ7Uw3Zrtrtm3ieS9iP1NB/fSJVaOSfqJa2R9I6km4vmq9ZCG+wHvAgMkPSypA8lvSTphII6tdAOxb4DbA38rqCsFtphfHY7RtKOkrpIGgR8Hbgi29dqPjOdMLaMrqR5fMUWA9tWOJbWplzbNOyvGpI+AYwAHomI6Vnx+tqgms6Rp0jDTH8D+pCG5d/K9tXE+0RSB1KyfGlEzGmkWrWfE0uAy0hTVb4GnE+av/iEpO2yOtXeBgA7kubpXgL8mjS/92HgakknZ3VqoR2KHQO8BTxYUFb17RARs0jTNQ4DXiUd7zXA0Iho+GLdaj4zN6vUC9WgUiuiq+JRtD6iRtom+wZ8D2lO3nGFu6iRNiBdyLIN8EnSBUEPS9ovIuZl+2uhHU4HtgBGlqlT1edERDwDPFNQNEXSn4GnSRfCnE2Vt0GmHakn7diIuDMr+5PSahJnShpDbbTDP0nakfTl4cpsqPWfu6jydshGXu4g9RYOJQ1NHwZcJ2lFRNTTitrBCWPLeJfSWf+2lP6mUEsWA6WWTdm2YH+bJ2lz0hV+nwQOiIiFBbsX0/j5AVV0jkREwzD8U5IeBOYBZ5D+c6z690k25FpH6lnrVDTfqJOkLsBSauicaBARMyT9Ddg7K6qFNniH1MP4cFH5Q6SrgXegNtqh0A9JifTvisproR0uJM1PPDQiGuZ2PyqpG3ClpFtoRZ+ZHpJuGbNJ8w6K7QG8UOFYWpvZQK9sSZVCe5CuFnvpXx/StmRDkHcAXwIOiYjni6qUOz8WRMQHLRxiLiLiPdK/b8Ocm1p4n3wS2By4ifQB17BB6nF9F/gcNXpOsG7vSS20wexGyht6i9ZSG+1Q6BhgZkTMLCqvhXb4HOnYVxeVPw10A7ajFX1mOmFsGfcC+0j6ZENBNuTwlWxfLbuXtObUUQ0FkjYDvg88FBEr8wqsOWRrLdaTJi0fFhFPlqh2L/AJSQcUPG4b4D+o4vND0sdJay6+nBXVwvvkWeCrJTZISeRXSf/h19w5IakvsBtpnivURhvcld0eXFR+MLAwIt6gNtoB+Oc50Jt/7V2E2miHN4C9JHUsKv9/wApS72Gr+cxUtqaPNaNs0eGZpPkIZ5O+QZ9PmrvSp0q+GZUk6bvZn18nDTseDywCFkXElKzOraT/IE8D5pIWtj4U2DciZlQ86GYkaSzpuEcC9xftXhgRC7Ok8nFgZ1IbvEta2LoPsGdEvFLBkFuEpLuAGcBzwPukxGA4sD3wpYj4W42/TwIYGRFnZ/er+pyQVE96r88gLRnyedLxLQO+EBFvV3sbwD8Xo36UtLZeHfB/wHdJizIfFxHja6EdGmRzNocBO0XEm0X7qr4dss/L20lTEq4l/V/YHzgBuCIiTsnqtY7PzEou+lhLG2nOwR2kD8ulwN0ULWJdjRvpQ7/U9lhBnS2Ay0nfrlaQehj65R17Mx3/vDJtcF5Bva7Ab0nfIJeRfYjkHX8ztsPppF+yeC87vjmkK4V7FtWr5ffJBUVlVXtOkD7onyNdLb0aeAUYB+xQK21QcIzbkK6EfZM0pPgccHQNtkMHUmfCfWXq1EI7fAt4LGuLpaRRieOB9gV1WsVnpnsYzczMzKwsz2E0MzMzs7KcMJqZmZlZWU4YzczMzKwsJ4xmZmZmVpYTRjMzMzMrywmjmZmZmZXlhNHMWpykYyVFI9t7ecfXWknqWabdCrfHsvqPNfzdWki6T9JVOb32cEnPZYtAm9km2CzvAMysphwFLCwq+zCPQNqI14EvF5U9AYwnLYTe4P3s9vgKxNRkkvYHvgF8KqcQriMtIv9j4L9yisGsKjhhNLNKejYiXso7iEqS1Ck28vdes8et83vk6dfleDVK/E55RLywUUG2nNNIv+Txah4vHhHLJf0eOBUnjGabxN30ZtYqSGqXDanOk/RvBeWfk7Rc0iUFZfMk3SRpkKSXJK2QNEPSV0s87w8lzczqvC3pD5J2KKpztKRnJH0gaYmk5yUNKdhfcqg3i2N8wf2Goff9Jd2eDbc/VbD/AEmPSloq6R+SJkn67CY0W3E868QpqV8Wz+GSrpe0WNK7kq6Q1F7S3pIez2KZLengEs+5UTFL2pH0s2c3F5VvL+l3kl6TtFLS65Lul7RdQZ3Oki6SNFfSquy2rnhoWVJ3SddKeiV7rleyf99OBdVuBfaQtG+TG9LM/oUTRjOrpPaSNiva2gFExFrgh8DWZMOtkrYgfeDPBuqKnusA4JSsfACwEnhQ0mcaKkgaDPwB+F/gO8AZwMHAFElbZXX2A24CpgCHk4bNfwN02YTjrAfmAt/NXhNJ3yb9Fu4H2XEenR3rVEk7b8JrNcVo4B/A94GrgZ9nZb8n/Vbvd0i/13unpI81PGgTY/4G0B54vKj8D6Rh9tOyOieRpil0zl5zM2AS8FPgSlLSeQPwS6DwS8O2wP9kx3Q5cAjwC9JvFHcseL1nSUP231xPvGZWTt4/vO3Nm7fq34BjgWhku7+o7hFZ+XHAOFKysltRnXnAKqBHQdnWpKTnD9n99sCbwOSix+6XPf9J2f1TgcXrif8x4LES5fOA8SWO84oSdV8CHi0q2wZ4Gxi9AW0ZwAVNiRPol9X/bVG9GVn5fgVlfbKyHzdHzMBY0tB5cfkHDW3fyON+lMWxf1F5XfZvvl12fwSwBvh8E9psKvBQ3u8Db97a8uYeRjOrpCOAvYu2nxdWiIi7SD2MY4FBwM8i4m8lnuvJiFhQ8LilwAN8dJHIZ4DtSL19hc//ODCf1EMJMA3YNhviPlTSpvQsNrir8I6kT5Mu/Kgv7F0FlpEuYtm/GV6znAeL7r8I/CNri8IygJ2bKeYdgUUlyqcBp0k6OZtuoKL93yT9+/xP0es+ROo93CerdxAwLSKeWU8cZHHs2IR6ZtYIJ4xmVkmzImJ60VbqIpjfAZ2AtyiaA1fgzUbKPpH93TW7fb1EvTca9kfEFNIw9M6kRG+RpEck9WnSEZVW/JoN8/NuBFYXbYcC3TbhtZri3aL7q4B1ljOKiFXZn5tnt5sa8+akaQLFvg/cSxo+fg54VdI5BfMTtwN2KfGaT2f7uxXcFl9x35jlwBZNrGtmJfgqaTNrVSR1Js2rmwV8Gvg1MLxE1Y83UtZwRe7i7Hb7EvW2B6Y33ImICcCEbF5jP+AiYKKknSLNrVxBGoot1rVEGaQh1ULvZLdnAo+UqL+qRFneNjXmd4BexYUR8RZwAnBCNt/0x8CvSL2AY7PHzQW+18jzzstu3+ajLwfr0zWrb2YbyQmjmbU2V5ISgb1IPVmjJU2KiIlF9faRtHNEvAIgaWvg26RhaYA5pB7HAaReMrJ6+5J6sC4rfuGI+AC4X9Inszi6kRKZ+cCRkjo29MQprTG4dROPaQ4p0ekdEb9u4mPytqkxvwgcIWmziCi51mZEzAHOkjQUaLjyeiJwJPBBRLxY6nGZh4CzJe0ZETPXE0svPuqhNLON4ITRzCppr8KrcAtMj4gPJR1Jujr2RxHxf8AYSQcB4yX1yXqnGrwJPCTpPNLQ5+nAlsD5ABGxRtI5wPWSbiJdCf0JYCTwd7J1+SSNIPVMTgZeA3YiXbn7bEQ0zMG7FRgM/DZbRqcX6QrtJU056IgISScA90jqCPyR1OP1cWBfYEFEXN6U56qUZoj5z6Sewz6ki2xQWi7pEdK80hdJQ82HAduSEkCyfccBj0q6DJhJuur5U0B/4PCIWAZcQbpq+xFJFwDPAx/Lnm9oNqeVbE7qbsClm9omZrXMCaOZVdLtjZR3z5bQ+Q1QHxE3Few7jjTXbbykb0dEw3DvFNJVwReSkrwXgG8VXiATEeMkLSMt4XIP6Qrd/wZ+kfUmQlon8SRSAtKVNG/yIdIyLg3PMznrBTuV1Pv1DGmZmTuaeuAR8d9Zr2QdaZmYLUhzKZ8Ebmvq81TSJsY8lZSA/wdZwkga2p9BuphpF2AtqSdzYETck73m6mw9yDNISXov0pJAL5N6j1dl9d6T9BXggqxuN9KXiD+x7nD5t7P761yIZGYbRh/932tm1jZImgc8HhE/zDsWa1zW+zuQtCxSLh82kh4E3o6IH+Xx+mbVwldJm5lZS7mCtAD6kXm8uKS9gK+ShsbNbBM4YTQzsxYREUtIC3F3XF/dFrI9cFwjSzeZ2QbwkLSZmZmZleUeRjMzMzMrywmjmZmZmZXlhNHMzMzMynLCaGZmZmZlOWE0MzMzs7L+P5+hsKuZ45q4AAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.figure(figsize=(10, 8))\n", + "plt.rcParams['font.size'] = 16\n", + "\n", + "plt.plot(exp_times, means, '-ko')\n", + "plt.xlabel('Exposure Time (sec)')\n", + "plt.ylabel('Mean Counts')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "anaconda-cloud": {}, + "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.3" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +}