{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1A. Schechter function:\n", "$$\\phi(L) = (\\frac{\\phi_*}{L_*}) \\exp(-\\frac{L}{L_*}) (\\frac{L}{L_*})^\\alpha$$\n", "Since we will just be looking for fractions of the cumulative function, the normalization doessn't matter, and we have:\n", "$$\\phi(L) \\propto (\\frac{L}{L_*})^\\alpha \\exp(-\\frac{L}{L_*})$$\n", "If we make the approximation that a single power law represents the luminosity function (as was suggested that you can do, but see below for numerical calculation), using $\\alpha$=-1, and removing exponential term:\n", "$$\\int_{L_1/L_*}^{L_2/L_*} (\\frac{L}{L_*})^{-1} dL = \\ln \\frac{L_2}{L_*} - \\ln \\frac{L_1}{L_*}$$\n", "\n", "We have $M_{lo}=-5$, $M_{hi}=-23$, $M_*=-21$, $L/L_* = 10^{(-0.4 (M-M_*))}$\n", "\n", "$$\\int_{L_{lo}/L_*}^{L_{hi}/L_*} L^{-1} dL \\propto \\ln (\\frac{L_{hi}}{L_*}) - \\ln (\\frac{L_{lo}}{L_*})$$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "total: 16.57861266955713\n" ] } ], "source": [ "#analytic solution assuming power law, get total integral first\n", "mlo=-5\n", "mhi=-23\n", "mstar=-21\n", "lhi = 10**(-0.4*(mhi-mstar))\n", "llo = 10**(-0.4*(mlo-mstar))\n", "tot=np.log(lhi)-np.log(llo)\n", "print('total: ', tot)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now get the 25, 50, and 75th percentiles:\n", "$$\\int_{L_{lo}/L_*}^{L_p/L_*} (\\frac{L}{L_*})^{-1} dL = \\ln L_p/L_* - \\ln L_{lo}/L* = (p) (tot)$$\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.25 -9.50\n", "0.50 -14.00\n", "0.75 -18.50\n" ] } ], "source": [ "for perc in [0.25, 0.5, 0.75] :\n", " lum = np.exp(perc*tot + np.log(llo))\n", " m = -2.5*np.log10(lum) + mstar\n", " print('{:.2f} {:.2f}'.format(perc,m))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we include the exponential term, you can do the integral numerically:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.25 -8.84\n", "0.50 -12.69\n", "0.75 -16.55\n", "1.00 -23.00\n" ] } ], "source": [ "def phi(m,mstar=-21,alpha=-1) :\n", " \"\"\"Schechter function in magnitudes\"\"\"\n", " return 10.**(-0.4*(alpha+1)*(m-mstar))*np.exp(-10.**(-0.4*(m-mstar)))\n", " \n", "\n", "# simple integral over full range, using small steps\n", "tot=0.\n", "for m in np.arange(-5,-23,-0.001) :\n", " tot+=phi(m)\n", "\n", "# now do it to get the percentiles to nearest 0.001\n", "tmp=0.\n", "p=0.25\n", "for m in np.arange(-5,-23.001,-0.001) :\n", " tmp+=phi(m)\n", " if tmp/tot > p :\n", " print('{:.2f} {:.2f}'.format(p,m))\n", " p+=0.25\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, note that if you express the Schechter function in magnitudes:\n", "$$\\phi(M) \\propto 10^{-0.4 (\\alpha+1) (M-M_*)} \\exp(-10^{-0.4 (M-M_*)})$$\n", "then, for $\\alpha=-1$, you have\n", "$$\\phi(M) \\propto \\exp(-10^{-0.4 M})$$\n", "and you can integrate analytically" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1B. To get luminosity integral ($\\phi_L$), multiply number at each luminosity by the luminosity:\n", "$$\\phi_L(L) = L (\\frac{\\phi_*}{L_*}) \\exp(-\\frac{L}{L_*}) (\\frac{L}{L_*})^\\alpha$$\n", "For the case of $\\alpha = -1$, this is particularly easy, since the $L$ cancels:\n", "$$\\phi_L(L) \\propto \\exp(-\\frac{L}{L_*})$$\n", "So the integral is easy and analytic without any approximation:\n", "\n", "$$\\int_{L_{lo}/L_*}^{L_{hi}/L_*} \\exp(-\\frac{L}{L_*}) dL = -\\exp (-\\frac{L_{hi}}{L_*}) + \\exp(-\\frac{L_{lo}}{L_*})$$" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.25 -19.64\n", "0.50 -20.60\n", "0.75 -21.35\n" ] } ], "source": [ "tot=-np.exp(-lhi)+np.exp(-llo)\n", "\n", "for perc in [0.25, 0.5, 0.75] :\n", " lum = -np.log(-perc*tot + np.exp(-llo))\n", " m = -2.5*np.log10(lum) + mstar\n", " print('{:.2f} {:.2f}'.format(perc,m))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1C. These results show that the bulk of the number of galaxies are at low luminosity, but the bulk of the luminosity comes from higher luminosity galaxies\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "select galaxy properties from SDSS" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/holtz/anaconda3/lib/python3.7/site-packages/requests/__init__.py:91: RequestsDependencyWarning: urllib3 (1.26.9) or chardet (3.0.4) doesn't match a supported version!\n", " RequestsDependencyWarning)\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "62710\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/Users/holtz/anaconda3/lib/python3.7/site-packages/astroquery/sdss/core.py:877: VisibleDeprecationWarning: Reading unicode strings without specifying the encoding argument is deprecated. Set the encoding, use None for the system default.\n", " comments='#'))\n" ] } ], "source": [ "from astroquery.sdss import SDSS\n", "zmin=0.015\n", "zmax=0.05\n", "sql=' SELECT \\\n", " objid,ra,dec,petroMag_g,petroMag_r, \\\n", " extinction_g,extinction_r,z,class \\\n", " FROM SpecPhotoAll \\\n", " WHERE \\\n", " z > {:f} AND z < {:f} AND class = \"GALAXY\" \\\n", " AND (ra between 140 and 240) and (dec between 0 and 50)'.format(zmin,zmax)\n", "dr17=SDSS.query_sql(sql, data_release=17) \n", "\n", "print(len(dr17))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Determine absolute magnitudes, first correcting for foreground extinction, then using distances from Hubble law" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(-0.1, 1.2)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEHCAYAAACncpHfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXhU933v8fdPG2gBtCBjFoFAYFzsYGILL9iOF3CcNJT0psFOuoS4z32wn9u6NPFte924t3bbW7dNE4c2vbFpGj84SZNCblO7aUiCCbEBG2yBMbExm5BALMaDJATaGCR+94+ZczgazUgjaTRnRufzeh4eNMuZ+c6ZOb/v7/y2Y6y1iIhI8OT4HYCIiPhDCUBEJKCUAEREAkoJQEQkoJQAREQCKs/vAIZi8uTJtrq62u8wRESyyu7du89aaytj78+qBFBdXU1dXZ3fYYiIZBVjzLF496sJSEQkoJQAREQCSglARCSglABERAJKCUBEJKCUAEREAsrXBGCMKTXG/MAYc8AY854x5jY/4xERCRK/zwDWAj+x1l4L3AC853M8IjIKWjrCPPdKPS0dYb9DEQ/fEoAxZiLwEeBfAKy1YWvtOb/iEZHRs7Guiac3HWBjXZPfoYiHnzOB5wAh4HljzA3AbmCNtbbDx5hEZBSsrK3q879kBj+bgPKAG4FvWGs/DHQA/yv2ScaY1caYOmNMXSgUSneMIpIC5cUFPHxXDeXFBX6HIh5+JoATwAlr7a7o7R8QSQh9WGvXWWtrrbW1lZX91jISEZFh8i0BWGvfB5qMMfOjdy0F9vsVj4hI0Pi9GuijwHeNMQXAUeAhn+MREQkMXxOAtXYvUOtnDCIiQeX3PAAREfGJEoCISEApAYiIBJQSgIhIQCkBiIgElBKAiEhAKQGIiASUEoCISEApAYiIBJQSgIhIQCkBiIgElBKAiEhAKQGIiASUEoCISEApAYiIBJQSgIhIQCkBiIgElBKAiEhAKQGIiASUEoCISEApAYiIBJQSgIhIQCkBiIgElBKAiEhAKQGIiASUEoCISEApAYiIBJQSgIhIQCkBiIgElBKAiEhAKQGIiASUEoCISEApAYiIBJQSgIhIQPmeAIwxucaYt4wxP/I7FhGRIPE9AQBrgPf8DkJEJGh8TQDGmBnAJ4Bv+hmHiEgQ+X0G8DXgj4HLPschIhI4viUAY8xy4ANr7e5BnrfaGFNnjKkLhUJpik5EZOzz8wzgdmCFMaYR+D5wrzHmO7FPstaus9bWWmtrKysr0x2jiMiY5VsCsNY+bq2dYa2tBj4D/Nxa+9t+xSMiEjR+9wGIiIhP8vwOAMBa+wvgFz6HISISKDoDEBEJKCUAEZGAUgIQEQkoJQARkYBSAhARCSglABGRgFICEBEJKCUAEZGAUgIQEQkoJQARkYBSAhARCSglABGRgFICEBEJKCUAEZGAUgIQEQkoJQCRMaKlI8xzr9TT0hH2OxTJEkoAImPExromnt50gI11TX6HIlkiI64IJiIjt7K2qs//IoPRGYBIlknU1FNeXMDDd9VQXlzgU2SSbZQARLKMmnokVdQEJJJl1NQjqaIEIJJlnKYekZFSE5CISEApAYiIBJQSgIhIQCkBiIgElBKAiEhAKQGIiASUEoCISEApAYiIBJQSgIhIQCkBiIgElBKAiEhAKQGIiASUbwnAGFNljNlqjHnPGPOuMWaNX7GIiASRn6uB9gCPWWv3GGMmALuNMZuttft9jElEJDB8OwOw1p621u6J/n0BeA+Y7lc8IiJBkxF9AMaYauDDwC5/IxERCQ7fE4AxpgT4f8AfWmvPx3l8tTGmzhhTFwqF0h+gyAASXZ9XJBv4mgCMMflECv/vWmv/Pd5zrLXrrLW11traysrK9AYoMghdn1eymW+dwMYYA/wL8J619qt+xSHiaOkIs7GuiZW1VZQXFyS1ja7PK9nMzzOA24HfAe41xuyN/vtVH+ORgBtObd65Pm+yCUMkk/h2BmCt3Q4Yv94/Uwyn1imjQ7V5CRrfO4GDTm3ImUO1eQmaQc8Aom31M6y1KqFGgWqdIuKXQc8ArLUW+I80xBJIqnUOj3f4pYZiigxPsn0AO40xi621b45qNCJJcprOHM7fD99V41dIIlkn2QRwD/CwMeYY0EGk89ZaaxeOWmQiA4jXdKZmNJGhMZEWnkGeZMysePdba4+lPKIB1NbW2rq6unS+5ZgWxBFIQfzMIsaY3dba2tj7kxoFZK095vwD7vf8HThjqb05iCOQxuJnHku/SUmv4cwDeARYl+pAsoW37Tnb25uDOAJpLH7msfSblPQaTgII9OStsVSAOCOQgmQsfuax9JuU9Eq2D+CLnptFQCfQBuy21u4dpdj6UR+AiMjQjagPAKgl0vQzHSgDVgN3A/9sjPnjVAUpIiLpk2wCqAButNY+Zq19jEhCqAQ+Anx+lGKTLOFXJ6Q6P8cufbfpkWwCmAl4v4lLwCxrbRdwMeVRSVbxa2TNWBzRM9YlW7Dru02PZDuB/5XIbOAXo7d/DfieMaYY0EXcA6ylI0xnuJc1S+emvRNSnZ/ZJ9kRS6P53WouyBVJJQBr7V8aY34M3EFkFNAj1lqnN/a3Ris4yXwb65pYu+Uwj3/82rQfTGNxRM9Yl2zBPprfrYbNXpH0ctDW2t3W2rXW2q95Cn8JuJW1VTz+8WszuhaeyvbkbGubzqR4Y2vesbGla4G/bPjNpouuBxBgzkFWH2of9sGWDauZprI92Y+26aEUhrHPHWm8Qy2I4z3fuW/9aw1uLC0dYf7ge2/x9KYDrH+tsU+sj23Yy7Ov1PP0pgN8bfMhHnr+DepD7cOKf7ifI9XbZ6rcJ5980u8YkrZu3bonV69e7XcYY8YLrzfy9KYDNLV08u2dxykvLqC2ujyl79HSEeaF1xuZU1lCYUFuSl872fe9btokyosLWFlblTCGZOOcU1nivlbXpd5R+WyxsTjf0/j8HPadaIv7fs42m989wz/9oh5r4c5rKplTWcL4/BzOd/VQd6yV+VdPcLeNfZ94tx/bsJdv7zzO+PxcNv3yFE/9534Wzihlamlh3NidWA+ducBHrqmksCCXZ14+xNoth5k/ZQJ3z6/kfHcPL7zeyPYjzQBYa7lxVhlvNrTQ3BFmb1Mb1lpOt3XTcLaDg2faORpqJ3ThInMqS/rs91NtXTy2YS8zy4t46e1TlBUXsLGuqd8+cuIqLy6grrGFv//ZIcbn53BbzeSkvhPv9s4x4t1fo/VbSJWnnnrq9JNPPtlvBQffLgkp/nNOgZctmMKtc86wsrYq5R1kI21vTTae2OfFvu9g750oTud1ly2Ywsv7I/vIefy5aC11uJ8t0eeLjcX5npo7wqzdcoDmjjAV0SQUu01ZUT4Am945DQYK83MAw7ptRwHYd+IcTyxfwMv7z9DcHmbdtqN0hnu565pKPvetXbRf7OVkaxd/8evXs7Guia0HQ9wzv5LWjou8sPM4AF/csJePLpjC2yfamFVexOEPLtAZ7uVEaye5OZBjYOvBEDf+5Wbycw2XeiOTTb+z6xjj8qDrUt998EZjK/c/8wo9lyO3DbDn+DkAznf3APDOyXNsP9LMi2+doKG5k65Ll/nW9gbOtl+k18IvDoawwE/fPc2e42189WcHmVFeyPnOHnLzcrilupRbZpdzsrWTd06dB2Db4bN0hS/TfamXg2fOc8OMMh65u//ZbKKBDmNhSfKkZgJnCs0EHn1Oofb4x69NyQ95uAnF2a4z3Ot2Mg8UT2zc8QruoSSQ2Ne9Z34lWw+G+sQx0Gdr6Qiz/rUGwLBqSbX7eOw2f/3j91j36lEeqJ1B6MJFVt1WzTdeqQcsf/2phZQVFbjPf/YX9azbdpQJ43O50N3LLbPLuXVOOSsWTeff3mzipb0nef/8RfJzDJcuXzmuq8oKmTAul4bmDrouJXe85wCXk3pmZjHASEq0eL+zRMeE97sEMnpkUaKZwEoA0kemDJF7ZvNB1m45wuo7Z1NRMm7EBfhwEpq3EL/rmkr+8eeHeWL5AmoqSwYt4DvDPazdcgSAW2aX843fvgmAxzbs7ZNIfuufd7KjvpmCXAj3QlF+Dp2XIkXvxPF55OVAS2ekFjwuFy729o9zpIVeEE0qzKOtq8dNptMmjWf5wmkJzwCypaBPJFECUBOQ9DHU4XcjreEn3i6y5mBhQV7ceGK3TxR3ssMOEzXFrN1yhDVL5/KPPz/M1oMhbp1zhrLaAh791z3sqG92ty8qyGXZgin8+YvvsP1IM5+8YSq5Bnot7Gpo4fPPv8HZCxc51dZNybhcfvru++Qa2FXfjCFS+ANu4Q9Xmj8c8Qp/UOGfSH4uXIrZZ2WFeRSNy+P//PqHWP96I1sPhgBY9itTqCiJ//v1/rZS0eyXSZQAZESG28Y/2HarllRTVJCbsOBO9n2TTWjxXs/b9r71YIjqiiKWLZjCs7+odwv/G2eWsvnd0+x/v50X955k/+kLALz49uk+r7/vRJv7d/vFXvYcP+e2c8vIFeQawtG+hqICQ2fYMnFcPs2dkQ4H5yyptauH1q4e3mo6x7yrSjjyQTtNrV28ejhEY3Mn4N8ENT8oAWS5kTTZpKK5Z9mCKew82syyBVOGtN3K2io6w710hnto6Qj3e//BCu5UHYj1oXb+6kf7efTeef3Ghjujff7ge3sAaGzu5I82vs37bd2eV7Dsfz8yRPHg+xdGFIsMX2/vlfOgznDkb6cfZFJhHvdfdzUb6k4wZcI4ppcV8nr9Wd5obAWgprKY+lAH98yv9HWCmh+UALKct+a6srZqSAX6cGrvsUnj5f1n3KaRmrtKko67vLiAooJcnt50gKIEzTyDbZ+KUUV/9aP9bjPA8w/d3O/561+LDFecND6Ptu6efrX20IUwd9RUsL2+GacMKszPoetSNnahZq9eIC8Hei7Dh6sm0dp5icbmTmoqi1n3uVrKigqoqSxxBxUA3DK7jFvnTGbFomlJDRQYi5QAspy3JjzUAj3ZWrRTS3aGD8YbojjUmrhfawjFJsx5V03gUu9lnli+oM/IoZf2nqIrHBk7D3D/9Vfzen0zTa1dLJw+ibxcw6lzXTS1dvHBhe4+76HCP/3G5Rou9lpumV3OmfPdNDZ3cs/8Svc36wzfbekIs/tYC9uPNHPrnAq+cN81AEOqvIwlgU4AmTLiZSScZgqn4ILkC+Nka9FXasn7+coDi/q8h/MaLR1hntl8kNhRMYk4awitWTovLd+Bt3B34t9Y18S6bUd5/OPXUlZU4I7Q2Xm02T0rcEydVMinbpzO2i1HOHO+mzMXLnLVhEi8F3sSd8PmGLisXtoR8Y5yKhmXy3XTJrGrocV9fFZ5Ed96aDEv7T3Fi3tPujX/rzywqF+lqLy4gH/47I19RvUEWaATwFhZFGq0J1s9sXwBEDkD8Bb4z71S726z/rUGd9hjUUFuvzhih9J1hntYs3QeYAeMPTLcshGwrFoye9hJwtlHze1hDn9wgcXV5XSGe1l95xw6wz2sf63RnfT0xPIFLJxxitaOMO+caiM/N4cVi6YB8IPdJzh5LlLjP985+LIAKvxHztmFOcCvfmgqn1k8k69uPsS00vHupLWyogL2nTjXp9nHqRwB/fp2svl4T6VAJ4Cx0qOfzOcYqJAfLIHUVJb0ax/vv01k2OYdcyfHjSN21uTaLUfcwraoIG/A0T5Om+1w+gocTmd13bFm9hxvY8/xVtq6erhjbgXbjzSz+s7Z3F5TwbwpE6JbWI6ebXfb/D+xdhszKwrdwh+gO8GwTEmdcXk5XOy5TFlRPq2dl9hQd4JT57rZfuQsa5bO4+8+fQMQGZ7pJPCvPLDI/Y2rsB9YoBNApv04htok5X3+UJc68G7rjMhpbr/IM5sPsWpJtbtNolgi2/TQGe6lpSPcZ9hmvOd7RwuVFRV4mln29zlgYz+XExvYPkkimX3lfY7TWb1w+iQA2roiY+zPd11iVnkR3Zcus6O+mR31zWx57wz1oQ4At/O3u+cyh850DLiPJTVuri5j0cwy3j3Zxo76Zrei8NLek4ChK9zD9iNn8c6A8FaCsrU51w+BTgCpkqq+BKeQ3nm0uV+hONDzIXHTT7y279htV9ZWse/EuT7t3s7tznAvXeEe9p++wFOfvI6aykhnWWQUT150FE/i8foOtwCecYqigly3WWnrwRAb65r6xB/7uZyOuqF8dmchM3eiTzQBtUUXoqkqK2RWRZG7INkP3zrBwukTaO3soT7U4dY427p73A5GGZ7icbnMKit0h8s6yorymTpxHPvfb+fGmaXUziqnsCDHbeqLPa6+cN98IPLdOrPDHZlWmcsWgUwAiabxD1eq+hJW1la5NePYQjHR873/w+CLonm3dWrwz0ZPnyNry1QA1j2dBsu6bQ0A/O//eIePXFPpvrb3/Z0+gM5wD6uWzO637rvT3l7X2MKO+mZ2Hm3mieUL3EXoBvtcyXx2L6dN32mSchY2e6B2Bg1nO7htTgUP313D73xzF6faurlwsZd9JyPj+K+eOI73z1+50qkK/6G7cWYpYNlzvI3funkWj9xdw/rXGth64AP2nTzPjTMn8c1VkWZF72+lPtTOYxv2uktuxDsGVNinTiATgDPFH+J3WA5VqvoSyosL3JELybxWvAPBKfA7w73u8gSxsTlJAgxrtxzm5uoyILIsr9P8422X7wr3sv/0BRZMm9hvRMWV9zfu/7FJx2nHv2d+JTvqm6mpLHbnDozeAR4ptG+aVdonWb16KBRpS959gpqrSlh+wzTWvXrUbWsGaL/Yk/BVvSqL8wl1XBr8iVkq3hpDRXk5dEb3U+zj1RVFfGReJWXFBXSFe1m37Si3zI7U6gG+cN/8fpUD6Fsx8Y44izcvQ1IrkAnAqf2CSUkHcCprJMm8lnMG0xW+TGFBbp82e6fA7wz3JDwrcQroNUvnsmbpPHYejTSDvNHY6r7GzqPNLK4u55VDIQoLcvmHz36Y1s4wh89ciDvr19sH0NoZdtv7I7X/yIifFYumceucMyyuLucff37YfXw4zWcDnXVFLtphWLN0LquWzO7z2GMfnY/ZfJAFUye5331FcQGLq8v5258cwFpL96XL7Dt5ZemG8uJ8/uij8/nyzw7S4inwz2dpL7BTcI/LNUwYn8dZz2fyLt8cW/iXFuazYtE0Xnj9GNUVRdw8u5wDp89ztj3Msl+5ij+8b777HUaGBENejmHtliNuB7532HK879w74kxGn68JwBjzMWAtkAt801r7N+l4X297Yjbpu0TyEff+ougFKGI7eYsK8lhcXc5Dz7/Bo/fO483GFpYtmBJtjpmNUxTsamjhjrmTWTB1Ap3hXndBs7eOn+NctM28qCCPznAPWw+GOBp6k289tNjtD4jlbe93mnzWLJ1LTWUJZbVXxtvfOudMv7iTlWgJCm/bv/c6xU7CePzj1/Kd/35rvw70+lC7p28C/vzFd5hTWUJZUb7bJn3/9VNZ/1oD2w6fZc/xc0wcn0fIc4WoiePzKC/Kp7Gli1uqy3j7RBvdPZeZMC6XooI8zlyINCt5l1ouzDN0xZlHkJsDvZ75ZOPzcujuGf4Es7wcQ090TOqHpk8kLzeH/Nwc/uRj1/LKoRBd4R42vfM+Ta1dTCrM45OLpjM+L5fuSz0cPNPODVWlPBL9fqaXFrozau+ZX8m+k+eZXlYE4A4NXrVkNkUFeX2W43YMlLzjjTiT0eNbAjDG5AL/BNwHnADeNMa8ZK3d71dMmSBejTh2iWFnyOL00kJOtXW5NXXvrFrnTOKh599g68EQx1s6qQ91uH0Mzvr2a5bOY83SuYChteMi67Y18MBNM9z1UQDumFsRbedvBOBYSyd/9aPIKXpsbK8eClFTWcLtNRW0dlx0F03bdvise/q/9WCImsriaNwfsGbpvCGfiSVagsJ7EZOVtVcucLO4upx75le6CSO2EIptenCSxPrXGlj/WqPbV7Ri0XR2H2tl9Z1zWFJTwZP/+S43V5cxtbSIVUuqWf9aY2TYqjF87rZZFBbkAZa1W464HcuXiTSXNDZ3kpeXAz29lIzLpT263OeM0kKW1FSwYfcJFk6fyMTCfL5433we/d4eTp7rpnhcLh0Xe/s1wUybNJ5Tbd1MHJ9HT28vndG1/4sLcikZl0dFcT7nunrYd/K8u81XNx/kplnlFBbkcs/8q3hh5zEeXDyTP/3VX0n4u3R+C2uWzmXFouluP06ii/DEzrIdK8OvxwI/zwBuBo5Ya48CGGO+D3wSGHMJYCjNHPFqR1eabCILlnWGe9229EghfcgdFx37+p/68HR2Hm3ms4ur2FHfzKP3zuPWORUsWzCFeVOaePXQB3xw4SInz3VTVRa5zN8bjS38y+cjMyu9E7BWLamO9ge0uTVlb2y311S4QykBjDFuHHuOn3OblzbUNVEf6nCXWPbW1JOVqBCJHQ4Ye0EXZ9hp7BlEvKYHb18RRM60th0+y/YjzeTn5lBRUkBjcyefvXmm+12tWDSNH+07xa6GFnY1tLD6ztkUFuRyy+xydjW0UFVWyMevn8qDN1e5ScdpO4fIshvOxKbIZRctKxZN5+X9Z1h67RRe2HmM3/jwDOpD7eyob3aHTBbm57Bi0XT3NW+cWcqe4+eYXjoeg+HEuS4uXb7sNmEtnD6R8909zKkscedZrL5zDvfMr+TBxcnM44hcpKemssQt4GO/k0S/e3XiZg4/E8B0wHul6hPALbFPMsasBlYDzJw5Mz2Rpdhg7dXeEUnxCrbYQq2lI+wO06ypLGbO5OJ+46IdX9tymK5Ll/n61nrOdV1i4YxSd1jl4TMXeKvpSlv3ZWuZUVpIY3MnL+8/w6ol1Tz7Sj2PfLuORTPLeHBxFRUlkan0zgEdG++O+uY+hdyf/vsv2dXQwi2zy91aorPyYqJRQMlIVIjE3u+8dqSgvzLsFOhzBhGv6cHbV+TMWF5952zycw1PLF9AWVHffQCRMxPvMNL9py+w/chZbq+pAKCptYuKkgLKigpYOKOUhTMi8xKciXE76pt5eX+kc9xZLG/fiTb3bM1ZsXT9a43sqG9m0cwyKooL3KYWZ596m14az3bwRz94mz/7xALeamrFGUu/blsDH83LcWdkx+6T2H0R+1ljv7fYfT9WZtqPZX4mABPnvn4lmLV2HbAOIlcEG+2gRsNAp7zxRiQluuycd3bjVx5Y5LZ1L184tc9Sxt5tvvzpG/jihr1MLMzn3Ik2usK9bjvtE8sX0HFxHz2XLaELF2lq7eKOuZP5ndtmsWzBlD7j6N9obOXwmQvubSdG70G/akm1m5gqSiKrL37jt2/qE39sMhvJIlzJnFl544s3wmqg5BM79twZGeV9r9jLQzpnFt5lJV7ef8ZdYM6Z0OaMjHK+N6e93JsQvcnLud95b6fT3ens965fFNvJWl5cwJbH7gbg7muvAuCZzYein8D0mceRaFZ2bOGeysUGxT++XRLSGHMb8KS19v7o7ccBrLVPJ9pmLF4ScrA5CQNd0nCwyyCuvnMOhz+4QFVZES/sPMYdcydz06xS1m45wu01FdRWl7vv6V3xs6ayxH2NG2eW8n5bN1dPGs8f3DuP9a83us+JFwf0n0Ec7zOmYvJcqq9fPBLeWAYa5eJI1eRBbx+Hc8lKZ8XWgfbLUK+5LNktEy8J+SYwzxgzGzgJfAb4TR/j8cVgI5IGqkV5h1t6CxHnuT8/8AG7GlqYVBj5mm+aVcqqJbPdJoUd9c10hXupKInUzJ0mEO/Qza5Lvax79Sin2rr54sa90Tbkvss3xJt7EDsZzduW/oX7rklJ88BgNcx0rvYae2Yz2GeKfU5sAk6W8zrOWjjxziAG2i5yZjP4TG4Zm3xLANbaHmPM7wM/JTIM9FvW2nf9iidTeQt56Fu7jh25Eju0sbn9IrsaWmjr6uGe+ZWsWDSdjXVNPHrvPI6GOjjW0sn+021sP9LsFt5XmicibdJOx/CE8bm0dFyirCifrQdDrH+twW0ucAoPpznC6eSFK+3ot8wuY1dDK04rXyqaBwYraNPZBj3Sjs2RToAaagJyqEM22HydB2Ct/THwYz9j8NtgtVRvwXDrnIo+BVrsyJXYAu/Bm2fy7qnzXDd9Eo9EZ+Q6I2KOtfS9YIZ34ph3SYo75k4GYN5VE9hz/By/tnAq08uKaO4Is3ZLpNb/hfuucWuTztmFMwzTSSZ3zJ3cZ2JWOgqeTG2DjtckNtIJUCrIZTgCORN4NA212WGwEULzpkwg3HM57qiT2JErsQXey/vPsKO+mY9cU0l5dKTIzqPNrLqtGsBtbqi5q6RfJ6fTYeqMJvGOKikvLnBnenr77WOXsnA6fZ1kcue8yWldqTHZQjHdFwaK1/GvCVDiByWAFEtmlcp4F0ZJNEJo3atH3fHWiV7TkWgI5MraKupD7ax+oY76UAeXei3bj5xl4YxTfVba7Az38jeb3mPX0RaeeXBRv4k83hE7zkzPeGsMxRuxlMlXYEr3cMVUL0UiMlxKACk2WLNDvAujJJoI5X2t4dRSvQnhsQ17qQ91UFNZzIKpE/rNG/BeeAXgixv28u//4/aE7znQQnTQtyAdrCaeihr4SF4j3U1F2boUiYw9SgApNlhhF6+wWbZgSp/LK8Z7LWeYIQyvluptYy4rKui3nvrK2qroLNezANxcXd5nHkAyQxuHW5CmogY+1GspeKn9XILKt3kAwzEW5wFAcuPZ01FL9nZOAu5iX96La4/GePFUfTbvInAq0EWuyMR5AIExWAEXr+Ycu00qaqmD1bT7z3y9conH0Wwm8X62kSSDhTMmsXBGqdrVRZKU43cAQeAUvM4aNLGcAtBb4A22zXCsrK1yZ6q2dIR57pV66kPtPPdKfXQN/cQxxYtxNAz3czsja4oKcnVNWJEk6QwgDYZTex7ONoPVnuP1KThDNIfTdj4ahnumMdIzlHQPBRXJBDoDSIPh1J6Hs81gtWen1t/SEXbPBp5YvsBdKjnedt5t0mG4ZxojPUNx9t1jG/am7bOK+E0JYAzxNvHE4y3kAHcC0lceWJRwu6EklWy2srZqwESYqcbK/hd/KAGMIQPVgp0F3u6YO7lfITfQdvH6DbyFzWj0VfjBmS6P+m8AAAi3SURBVLA2UALNRGNl/4s/1Afgo8GWgk4lp5N0zdK53DlvctKF3GBzEUbS9p5p7e7ZOB8gU9c7kuygBOCjeGvCjJbY1SJH+hqOkRSaumLUyGVj0pLMoQTgo3SuCZOKsfapLmwytfaaaWcmIqNFCcBHfq0Jkyk170ytvWbK/hEZbUoAAZSpNW+H3zXwTN8/IqmiBBBAmVrzdvhdA8/0/SOSKkoAknFUAxdJD80DCKhMnkCUaF5CJscsko2UAAIqWyYQeQv9bIlZJFuoCShLjbSjNFuaWbz9AdkSs0i2UALIUiPtKB2tjs5Uj+CJncA21Jj9HlEkksmUALJUptaGUz2CZ6SJyu8RRSKZTAkgS2XqUMVMS0yZFo9IJtE1gUVExrhE1wTWKCARkYBSAhARCSglABGRgFInsKRN5AI4jYBl1ZLZGpYp4jMlAEmbyAVwDgNQVJCXkaOYRIJECUDSJnIBnF7AalimSAZQApC0iVwA5xq/wxCRKHUCi4gElBKAiEhA+ZIAjDFfNsYcMMbsM8b80BhT6kccIiJB5tcZwGbgemvtQuAQ8LhPcYiIBJYvCcBa+zNrbU/05k5ghh9xiIgEWSb0AfwusCnRg8aY1caYOmNMXSgUSmNYIiJj26gNAzXGvAxcHeehL1lrX4w+50tAD/DdRK9jrV0HrIPIaqCjEKqISCCNWgKw1i4b6HFjzCpgObDUZtOa1CIiY4QvE8GMMR8D/gS4y1rb6UcMIiJB51cfwNeBCcBmY8xeY8yzPsUhIhJYvpwBWGvn+vG+IiJyRSaMAhIRER8oAYiIBJQSgIhIQCkBiIgElBKAiEhAKQGkQEtHmOdeqaelI+x3KCIiSVMCSIGNdU08vekAG+ua/A5FRCRpuiRkCjjXt9V1bkUkmygBpEB5cQEP31XjdxgiIkOiJiARkYBSAhARCSglABGRgFICEBEJKCUAEZGAUgIQEQkoJQARkYBSAhARCSglABGRgApMAtCCbSIifQUmAWjBNhGRvgKzFpAWbBMR6SswCUALtomI9BWYJiAREelLCUBEJKCUAEREAkoJQEQkoJQAREQCSglARCSglABERALKWGv9jiFpxpgLwEG/4xiiycBZv4MYgmyLFxRzOmRbvJB9MY9mvLOstZWxd2bbRLCD1tpav4MYCmNMXTbFnG3xgmJOh2yLF7IvZj/iVROQiEhAKQGIiARUtiWAdX4HMAzZFnO2xQuKOR2yLV7IvpjTHm9WdQKLiEjqZNsZgIiIpIgSgIhIQGVFAjDGrDTGvGuMuWyMqY157HFjzBFjzEFjzP1+xeiVKF5jzH3GmN3GmF9G/7/Xzzi9BtrH0cdnGmPajTH/04/4Yg3ym1hojHk9+vgvjTHj/YrTa4DfRb4xZn001veMMY/7GaeXMebLxpgDxph9xpgfGmNKPY9l4rEXN94MP/YS7uPo46N27GVFAgDeAT4FvOq90xizAPgMcB3wMeD/GmNy0x9eP3HjJTLJ49estR8CVgHfTndgA0gUs+MZYFP6whlUot9EHvAd4BFr7XXA3cCltEcXX6J9vBIYF/1d3AQ8bIypTm9oCW0GrrfWLgQOAY9DRh97ceMls4+9RDE7Ru3Yy4qJYNba9wCMMbEPfRL4vrX2ItBgjDkC3Ay8nt4I+0oUr7X2Lc/Nd4Hxxphx0fh9NcA+xhjz68BRoCPNYSU0QLwfBfZZa9+OPq85zaElNEDMFiiOJq9CIAycT2908Vlrf+a5uRP4dPTvTD324sab4cdeon086sdetpwBJDId8F7l/UT0vmzwG8BbmfADHIgxphj4E+Apv2NJ0jWANcb81Bizxxjzx34HlIQfEDnATwPHgb+31rb4G1Jcv8uVmmg2HHveeL0y+dhzY07HsZcxZwDGmJeBq+M89CVr7YuJNotzX1rGtQ4zXmfb64C/JVJbTZthxvwU8Iy1tj3e2cFoGma8ecAdwGKgE9hijNltrd0ySmH2McyYbwZ6gWlAGbDNGPOytfboKIXZRzIxG2O+BPQA33U2i/P8jDn24sTrbJuxx16cmEf92MuYBGCtXTaMzU4AVZ7bM4BTqYloYMOMF2PMDOCHwOestfWpjWpgw4z5FuDTxpi/A0qBy8aYbmvt11MbXX8j+E28Yq09C2CM+TFwI5CWBDDMmH8T+Im19hLwgTFmB1BL5NR/1A0WszFmFbAcWGqvTBzK2GMvQbwZfewliHnUj71sbwJ6CfiMMWacMWY2MA94w+eYEor27v8X8Li1doff8STDWnuntbbaWlsNfA3463QU/iPwU2ChMaYo2qZ+F7Df55gGcxy410QUA7cCB3yOCQBjzMeINEOssNZ2eh7KyGMvUbyZfOwlijkdx15WJABjzH8zxpwAbgP+yxjzUwBr7bvABiIH+E+A37PW9voXaUSieIHfB+YCf2aM2Rv9d5VvgXoMEHNGGuA30Qp8FXgT2Avssdb+l3+RXjHAPv4noITIKKE3geettft8CjPW14EJwObo7/VZyNxjjwTxksHHHoljHnVaCkJEJKCy4gxARERSTwlARCSglABERAJKCUBEJKCUAEREAkoJQEQkoJQARAZhjLHGmG97bucZY0LGmB8NsM3no8/ZG13q9wvpiVYkeUoAIoPrAK43xhRGb98HnExiu3+z1i4Cbge+ZIypGmwDkXRSAhBJzibgE9G/Pwt8L9kNo0tSHwGmjkJcIsOmBCCSnO8TWftmPLAQ2JXshsaYmcB4IFOWdxABlABEkhJdm6eaSO3/x0lu9qAx5l0iq3qutdZ2j1J4IsOiBCCSvJeAvyf55p9/i16W8k7gK8aYeOvBi/hGCUAked8C/sJa+8uhbGStfZ3INWjXjEpUIsOkBCCSJGvtCWvt2mFu/rfAQ8aYCamMSWQktBy0iEhA6QxARCSgMuaawCLZyBjzEP3b9ndYa3/Pj3hEhkJNQCIiAaUmIBGRgFICEBEJKCUAEZGAUgIQEQmo/w8M/xqSXY+ApAAAAABJRU5ErkJggg==\n", "text/plain": [ "